Computers& GeosciencesVol. 23, No. 1. pp. 91-101, 1997 0 1996Elsevier Science Ltd Printed in Great Britain. All rights reserved PII: SCtW8-3004(%)00064-7 0098-3004/97317.00+ 0.00
Pergamon
GravMap
AND PFproc: SOFTWARE FOR FILTERING GEOPHYSICAL MAP DATA G. R. J. COOPER
Departments of Geophysics and Geology, University of the Witwatersrand, Johannesburg 2050, South Africa (e-mail:
[email protected]) (Received 15 December 1995; revised 22 July 1996) Abstract-A pair of programs have been written for the filtering of map data, in both the space and frequency domains. Although the programs are designed primarily for geophysical data, several of the filters are applicable to any type of two-dimensional dataset. The programs are suitable for application to small datasets, or for teaching purposes. 0 1997 Elsevier Science Ltd. All rights reserved
Key Words: Filtering, Magnetics, Gravity, Frequency domain, Space domain.
INTRODUCTION
cessed. Before the dataset can have its FFT calculated it must consist of a grid of 2N x 2M points (M not necessarily equal to N). A regional field comprising the average value of the map first is subtracted from the data, then all the edges are padded out to 2N points using a cosine bell. After the data have been filtered and the inverse FFT applied, the border is removed to restore the original grid size, and (if necessary) the regional field is restored. This process helps minimise edge effects in the map being introduced by the filtering process.
GravMap and PFproc have been written to perform basic filtering operations on small map datasets of up to 5000 observation points. This limit is set by the 640K memory limit of DOS. GravMap takes data values from scattered locations (such as gravity or magnetic stations) and grids and contours them to produce a map. This map may then be filtered using vertical continuation and polynomial surface fitting. In addition, linear profiles may be generated at any angle across the map. PFpr oc operates in a similar manner to GravMap, and the filters that may be applied to the data are strike filtering, vertical derivatives, edge enhancement, sun-shading, lowpass/highpass filtering, and pole reduction. Maps also may be displayed in a histogram-equalised form to enhance detail. Both programs allow combinations of the filters to be applied to the data. For example, the polereduced data may be lowpass filtered. The maximum grid size is 60 x 60 points, again due to the memory restrictions of DOS. A sample aeromagnetic dataset, consisting of just over 4000 points, has been gridded and contoured in Figure 1. The data consisted of flight lines 1 km apart, flown in a North-South direction. Hardcopy output can be produced on dot matrix and laser printers, and the filtered data can be saved to disk in ASCII format.
Vertical continuation
Upward continuation is applied frequently to magnetic data to enable surveys collected at different heights (such as ground and aeromagnetic data) to be compared. Upward continuation has the effect of smoothing the data, whereas downward continuation increases detail but also makes the map more noisy. Figure 2 shows the map in Figure 1 upward and downward continued. Vertical continuation is applied in the frequency domain using the following function (Gunn, 1975): A(& v) = A(u,
v,e-z-,
where u and v are orthogonal frequencies present in the map, A(u,v) is the amplitude present at those frequencies, and z is the vertical distance that the data are to be continued. Strike filtering
FREQUENCY DOMAIN FILTERS Several of the filters are performed in quency domain and therefore involve the fast Fourier transform (FFT), which places restrictions upon the gridded dataset to
Features on a map that have a particular orientation may be enhanced by using a strike filter. This is often useful to remove “geological noise”, such as dyke anomalies cutting across a feature of interest. Two parameters describe the filter: the direction
the freuse of a various be pro91
92
G. R. 3. Cooper
2024.8
2060.9
% 0‘
Figure 1. Contour map of aeromagnetic data from South Africa. Contour interval 100 nT. Distances .__.. which is to be enhanced, and the severity of the filter. The filter works by upward continuing the map, except for features that lie in the direction to be
enhanced. Hence the severity is a distance which is measured in the same units as the map itself. The algorithm is given by Syberg (1972): A(& r) = A(& “)e-h@+JsB+dl, where h is the continuation height, 4 is tan-‘(&), and 0 is the direction that is to be enhanced. Figure 3 shows the map in Figure 1 strike filtered. Vertical derivatives
The vertical gradient operator is applied often to magnetic data to sharpen the edges of anomalies and other features. The order of the derivative may be chosen. The higher the order entered, the noisier the map will become. The algorithm is given by Gunn (1975): A&, v) = A(u, v>(dm/,‘, where n is the order of the derivative. See Figure 4 for the first vertical derivative map of Figure 1. Pole reduction
Reducing a map of magnetic data to that which would have been measured had the observed geology been at the magnetic pole is a useful aid to interpretation since asymetric anomalies become symmetric. The pole reduction assumes that there is
no remanent magnetisation present. The inclination of the geomagnetic field over the map must be entered, as well as the bearing of the map. The algorithm is given by Spector and Parker (1979): A(u, v) =
A(u, 4 Sin I + iCos I Sin@ + v)2 ’
where Z and D are the inclination and declination of the geomagnetic field, and cp is tan-‘(u/v). Figure 5 shows the data in Figure 1 after pole reduction has been applied.
SPACE DOMAIN FILTERING
OPTIONS
Several of the filters that can be applied do not require the use of FFT, since they operate on the data directly. To distinguish these filters from those requiring FFT, they are termed space domain filters. Polynomial surface jitting
Orthogonal polynomial surfaces of up to fourth order may be calculated from the data using Gr avMap. The program reports the fit of the polynomial surface to the data, and stores both the surface itself and the residuals. To calculate the surface, the normal equations for the surface are calculated, then solved using matrix algebra. For example, a linear trend surface is given by:
Grdvbfap
1808.0
1845.3
and PFproc:
1882.7
software
1920.0
for filtering
1957.3
geophysical
1994.7
map data
2032.0
93
2069.3
0
od
2
t %
1808.0
1845.3
1882.7
1920.0
1957.3
1994.7
2032.0
2069.3
w
*10*
B
*lo2 o 1808.0
1845.3
1882.7
1920.0
1957.3
1994.7
2032.0
2069.3
1845.3
1882.7
1920.0
1957.3
1994.7
2032.0
2069.3
od 2
9 s 1808.0
*10* Figure
2. Vertically
continued
data
applied to map in Figure 1. A, Upward Downward continued 300 m.
continued
500 m. B.
G. R. J. Cooper
94
*1d 1808.0
Figure
1844.1
3. Strike filtering.
1880.3
Features
1916.4
1952.5
1988.6 I
lying at angle of 60” clockwise
from North
in Figure
1 are enhanced,
*lo2
1808.0
9 %
1845.3
1882.7
1920.0
1957.3
1994.7
2032.0
2069.3
1882.7
1920.0
1957.3
1994.7
2032.0
2069.3
9 CT 1808.0
1845.3
*IO2
Figure
4. First vertical
derivative
of data in Figure
1.
GravMap and PFproc: software for filtering geophysical map data
95
*lo2 1812.0
1848.8
1885.6
1922.4
188k.6
19i2.4
1959.2
19i9.2
1996.0
2032.8
2069.6
2Oi9.6
2Oi2.8
1956.0
N
*lo2
Figure 5. Pole reduction of data in Figure 1.
where Y is the data in the map and Xi and Xz are the two map coordinates. The normal equations which must be solved are then (Davis, 1986, p. 402): c
Y = bon + bl c
XI +
b2 1
x2
E-W enhancement [ -101
-101
1 1 1
0 0 0
-1 -1 -1
1 1 0
1 0
0 -! -1 1
N-S enhancement
[ NW-SE where there are n points and all summations extend from 1 to n. The equations are solved for the unknowns bo, h,, and b2. The value of the trend surface can then be calculated at any point on the map and subtracted from the original data to give a residual map. Figure 6 shows both a first order polynomial surface fit to the contour map in Figure 1, and the residual surface remaining after the polynomial surface has been subtracted from the original map. Edge enhancement This produces a map where the edges of features have been “sharpened up” using convolution kernels, a process routinely applied in the image processing of geophysical and satellite data (Richards, 1993, p. 123). Four kernels of size 3 x 3 grid points are applied to each grid point, and the greatest value resulting is retained. The kernels are:
-101
enhancement
[ NE-SW
1
-1
enhancement 0
-1 [ -1
1 0 -1
1 1 0
1
Figure 7 shows the effect of edge enhancement on the map in Figure 1. Note that since the kernels cannot be applied to the outer pixels of the map, edge effects are present there. Sun-shading Sun-shading is another technique for enhancing directional features of maps and images. It considers the map as if it were a 3D topographic surface, then calculates the length of the shadows that
G. R. J. Cooper
96
A
*Id
O 1808.0 d 0‘
1845.3
\o g-
1882.7
!k * 2
1957.3 I
1994.7 2032.0 I ^ I 033.G
803.3 796.7
803.3 796.7
790.2 783.7
790.2 783.7
777.1
777.1
I
2 c!
1920.0 1
,
,_.”
2069.3 1
o i E
I ‘0 ri ‘0
._._
764.1 757.5
764.1 757.5--
751.0 744.5
751 .o 744.5
738.0 731.4
731.4 738.0
“! zi B l
-g
*10*
a 1808.0 0s 0‘
1845.3
1882.7
1920.0
1957.3
1845.3
18i2.7
19io.o
19i7.3
1994.7
2032.0
2069.3
9
3
1808.0
*1d
Figure 6. Polynomial surface fitting. A, First-order trend surface fitted to data in Figure map generated by subtracting the trend surface from the original data.
would he cast by a “Sun” which occupied a particular elevation and azimuth. These values are then contoured. The lower the Sun elevation, the greater the “shadows” cast, and hence the greater the enhancing effect. Features lying normal to the sun azimuth are enhanced, whereas features parallel are reduced. The reflectance is calculated using an equation from Pelton (1987):
1.
B, Residual
1 +P.Po+Q.Qo
R=
J(1+P2+Q2)+,/m' where
PO = - cos q5tan 0, Q. = -sin
is the sun elevation and 4 is the azimuth East). The gradients
4 tan 0; tl
(measured
from
(measured
anticlockwise
are calculated
the vertical), from
from the follow-
GravMap
and PFproc:
software
for filtering
geophysical
map data
97
*IO2 1808.0
18445
1880.3
1916.4
1952.5
1988.6
2024.8
2060.9
1844.5
1SSCi.3
191b.4
195;s
1988.4
2024.8
2060.9
Figure
7. Edge enhancement
‘-
1.
of data in Figure
*1d 1808.0 ro d 2
9 ;s 1808.0
1844.1
1880.3
1844.1
1880.3
1916.4
1916.4
1952.5
1952.5
1988.6
1988.6
2024.8
2024.8
2060.9
2060.9
*1d Figure
8. Figure
1 sunshaded
from
sun direction of 60” clockwise from horizontal.
from
North
and elevation n of 30”
Cl. R. J. Cooper
98
A
1808.0
*lo2
1844.1
1880.3
1916.4
1952.5
1998.6
2024.8
206b.9
s 9
1808.0
1844.1
18803
1914.4
1952.5
1988.6
2060.9
2024.8
*1d
B 1808.0
1808.0
*1@ 1844.1
1844.1
1880.3
1880.3
1916.4
19 16.4
1952.5
1952.5
1988.6
1988.6
2024.8
2060.9
2024.8
2060.9
*1d Figure
9. A, Moving
average
low-pass
filtered map of data in Figure map of same data.
1. B, Laplacian
high-pass
filtered
GravMap
and PFproc:
software
for filtering
geophysical
map data
99
A Modified histogram
Original histogram
150
150 -
75
I
I
32
16
0
Class
Class
Old contour values
New contour values
14-
14
4-
4
1
47
“0 l
s s s 0
:--/ -6
I 17
1
1 33
-
I 49
I
;3
I;
1
44
Class
Class B Equalised map 1808.0
1808.0
1844.1
1844.1
1880.3
1916.4
1916.4
1880.3
1952.5
1952.5
1988.6
2060.9
2024.8
1988.6
Z&O.9
2624.8
*1@
Figure
10. Histogram
equalisation.
A,
Original and equalised Equalised map.
histogram
and
contour
values.
B,
G. R. J. Cooper
loo
Gr avMap uses the data from two grid cells around each interpolation point in a second-order polynomial interpolation to calculate the field value at the interpolation point (Press and others, 1992, p. 118).
ing kernels: -1 P=-2 -1
Q=
0 0 0
1 2 1
2 0
:,
:, -1
Histogram equalisation
-2
-1
See Figure 8 for a map of the data in Figure sunshading.
1 after
Low-pass and high-passjltering Low-pass filters are used to smooth noisy datasets, whereas high-pass filters are used to bring out detail in smooth maps (at the risk of enhancing noise). There are two techniques offered for lowpass and high-pass filtering: moving average and Laplacian. The moving average option replaces each grid point of the map with the average of itself and its eight neighbouring grid cells; this is a lowpass filter. The resultant high-pass filtered map is merely the difference between the original map and the low-pass map. The Laplacian operator produces a high-pass map and is given by Hord (1982, p. 78): L = V=Z = d’Z/dx=
+ d=Z/dy’
(8)
The low-pass map option is the difference between the original map and this high-pass map. The filters are implemented using the kernels shown: Lowpass filter
This is a modification in the way that the contour maps are displayed. Instead of the contour interval being constant across the map, it is changed so that there are an equal number of data points in each interval. This has the benefit of representing detail in the map equally at all data values (Richards, 1993, p. 97). To perform the equalisation, the data values are first sorted using the quicksort method (Press and others, 1992, p. 323). The number of data points in each contour interval is then: N = Npts/NI, where Npts is the number of data points in the map, and NZ is the number of contour intervals desired. Since the points have been sorted, the value of the first contour interval is simply the value of the data point after N points have been counted. The second interval is the value of the data point after another N points have been counted, and so on. Figure 10 shows the map histogram before and after equalisation, along with a graph of the contour values themselves. The equalised map is plotted also, using 25 contour levels.
CONCLUSIONS
[i Laplacian
high-pass
i
i]
filter
Figure 9 shows the results of moving average lowpass and Laplacian high-pass filtering the map in Figure 1. As with the edge detection and sun shading filters, the fact that the kernels cannot be applied to the outer pixels of the map causes edge effects.
The filtering of small datasets can be performed easily and quickly using Gr avMap and PFp ro C. The programs are undemanding in the hardware that they require to operate, even running on 8086 based computers, and so are well suited to field operations and teaching laboratories where state of the art computer hardware may not be available. They can be obtained via anonymous FTP from ftp.cs.wits.ac.za in directory /pub/general/geophys, or from the server at IAMG.ORG.
REFERENCES
Davis, J. C., 1986, Statistics
OTHER
FEATURES
Projile extraction across the map Profiles may be extracted at any angle across the contour map using GravMap, for example for modelling purposes. The profiles can contain up to 128 points, and may be taken across any of the different maps produced by the program, such as polynomial surface maps, vertically continued maps, and so on. The profile must be a straight line that starts and ends on an edge of the map.
and data analysis in geology (2nd ed.): John Wiley & Sons, New York, 646 p. Gunn, P. J., 1975, Linear transformations of gravity and magnetic fields: Geophysical Prospecting, v. 23, no. 2, p. 300-312. Hord, R. M. 1982, Digital image processing of remotely sensed data: Academic Press, New York, 256 p. Pelton, C., 1987, A computer program for hill-shading digital topographic data sets: Computers & Geosciences, v. 13, no. 5, p. 545-548.
Press, W. H., Teukolsky, S.‘ A., Vetterling, W. T., and Flannerv. B. P., 1992, Numerical recipes, the art of scientific computing (2nd Press, New York, 963 p.
ed.): Cambridge
University
GravMap and PFproc: software for filtering geophysical map data Richards, J. A., 1993, Remote sensing digital image analysis, an introduction (2nd ed.): Springer-Verlag, New York 340 p. Spector, A., and Parker, W., 1979, Computer compilation and interpretation of geophysical data, in Hood, P. J., ed., Geophysics and Geochemistry in the Search for
101
Metallic Ores: Geol. SUIT. Canada, Economic Geology Report 3 1, p. 527-544. Syberg, F. J. R., 1972, A Fourier method for the regionalresidual problem of potential fields: Geophysical Prospecting, v. 20, p. 47.
