GEOPHYSICS, VOL. 34, NO. 1 (FEBRUARY 1969), P. 39-53, 9 FIGS.

, 1 TABLE
Downloaded 01/14/15 to 157.211.3.15. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

THE EQUIVALENT SOURCE TECHNlQlJEt

C. N. G. DAMPNEY*

The inherent ambiguity of potential field inter- and above the terrain.
pretation can be put to advantage. Bouguer Once the equivalent source is obtained, the
anomaly measurements on an irregular grid and projection of the Bouguer anomaly onto a reg-
at a variety of elevations can be synthesized by ularly gridded horizontal plane is easily done.
an equivalent source of discrete point masses on a In addition, the equivalent source can then be
plane of arbitrary depth below the surface. By economically used to carry out vertical continua-
keeping the depth of the plane within certain tion. The technique is illustrated by a hypothet-
limits relative to the station spacing, we can ical example and a case history of a local gravity
ensure that the synthesized field closely approxi- survey in precipitous topography.
mates the true gravity field in the region close to

INTRODUCTION of equation (1.1) any value of g, (x, y, 2) for z <h

fundamental property of a potential field is A
is similarly defined uniquely.
the inherent ambiguity in determining its source Thus we fiild a surface density contrast distri-
from values of the field outside the source region. bution on an arbitrary plane which synthesizes a
For instance, given the vertical anomalous gravi- known gravity field. It is then possible to cal-
tational field intensity, g, (x, y, z) over a hori- culate from gravity measurements at one height
zontal plane, it is impossible to find the anomalous the gravity field at any other point in space.
mass distribution uniquely. The projection and interpolation of the field from
Roy (1962) has discussed in detail the ambi- known points of measurement may thus be
guity of the relationship between g, (x, y, z) and carried out.
c (c~,/3, h) in the equation

gz(x,y, 2) = K
J --oo J --oo I<%- o)* + (y - P)” + (2 - h)Zj3’?

where u ((Y, p, h) is the surface contrast density

distribution over the plane z= h, K is the gravita-
tional constant and the positive direction is By using the equivalent source u (a, p, h) as an
down. We note that the depth h of the apparent integral part of gravity field computations, two
source may take any value. advantages result: (1) All available information
To perform the inversion of equation (1.1) it is contained in the measurements of g, (x, y, z) is
sufficient to have a complete knowledge of g, on a used, and (2) we ensure the correct analytical
plane z=zr zi <h. Therefore a unique corre- gravity potential field behavior of the projected
spondence exists between the function g, (x, y, field values intrinsically as they are calculated
z=zi) and (T (x, y, h) ; and furthermore as a result from a gravitating equivalent source.

t Presented at the 37th Annual International SEG Meeting, Oklahoma City, Oklahoma, November 1, 1967.
Manuscript receivedby the Editor December4, 1967;revised manuscriptreceivedOctober 30, 1968.
* Dept of Geology,University of Tasmania, Australia; Contribution No. 187. Now with GeophysicsLaboratories,
University of Toronto, Canada.
39
 40 Dampney

REVIEW uous surface density contrast distributiona(c+,h)

There are, in the literature, many data-pro- on any plane zl=h. The equivalent source
cessing techniques (e.g. Henderson, 1960; Bhat- technique is based on approximating this con-
Downloaded 01/14/15 to 157.211.3.15. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

tacharyya, 1965, 1966; Keidell, 1966; Sax 1966) tinuous distribution by a series of discrete masses.
which require data to be put on a gridded hori- If we have N data points, we can calculate the
zontal plane. In practice, however, points of values of N point masses at a suitable depth which
measurement may be scattered horizontally and will then constitute the equivalent source from the
vertically over normal rough topography. equations (using the principle of superposition).
i\s Naudy and Neumann (1965) emphasize, the
gl = aim + 61121132 + . .
Bouguer anomaly must be regarded as being the
vertical gravity field at points on the terrain due
+ alkmk + . + (IlNmN
to the anomalous massesin the ground. The points
are not situated on the geoid or some other g2 = azlml + a2?m2 + . . .
reference surface. The observed values must not
only be interpolated horizontally from their orig- + a?kmk + ’ ’ ’ + a?NmN

inal coordinates of measurement, but be ver- . . . . . . . . . . . . . .

. (3.1)
tically projected onto the plane of the grid as
well. gi = Uii?YZi+ Uigm2 + . . .
Projection of Bouguer anomalies onto a hori-
+ aikmk + . . . + Ui.vWZ.v
zontal plane has received very little attention in
the literature. Strakhov and Devitsyn (1965) and . . . . . . . . . . . . . . .
Grant and West (1965, p. 269) have proposed the
only two solutions known to the author. In this
paper we set up an equivalent source to represent
and thereby interpolate and project the observed
Bouguer anomaly values.

uik
= {(xi- cYk)*+ (ri - @k)*+ (zi - /4*)3/z

The idea of an equivalent layer has been ex- and z=h is the horizontal plane containing the
ploited by Danes (1961) in gravity interpretation point masses mk at ((Yk,@k, h). The position of g,
incorporating borehole data and recently by Bott is (Xi, Yz, &I.
(1967) to interpret magnetic anomaly profiles. This can be written in matrix form
This demonstrates that representing potential
g = Am (3.3)
fields by an equivalent source may have many
applications. In fact, as Zidarov (1965) showed, which represents N simultaneous equations in N
the idea can be applied to electrical fields and to unknowns and is thus solvable.
the general Dirichlet-Neumann inverse problem Suppose the N data points gi have an average
as well as to gravity and magnetics. station spacing of Ax. The approximation of the
Zidarov’s (1960, 1965) papers give the general discrete masses mk to a continuous distribution
outline of representing potential fields by equiva- will be a valid representation of g if the mk are
lent sources. sufficiently far below the surface relative to Ax.
Make the average separation of the N masses
THE EQUIVALENT SOURCE TECHNIQUE
AZ so as to correspond to the average station
From equation (1.1) we see that the anomalous spacing of gi. Consider the anomaly gz (x, y, s)
gravitational field intensity g, (x, y, z) can be at (x, y, z) resulting from the discrete masses
represented by or be synthesized from a contin- mk at (ffk, Pk, h)

N
mk(h- 4
k=l ((x - ak)2 + (y - ok)’ + (2 - h)2f3’2 ’
 Equivalent Source Technique 41

Taking the Fourier transform of g, we obtain:

G,(u, z’, z) = 27rK esp - e 4-x)’ + (zlAx); \I X 2 mkek, (3.5)

Downloaded 01/14/15 to 157.211.3.15. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

I k=l

where tion of a gravity field, due to a line source, below

the source depth induced violent oscillations in
ek = exp (-&U - &v)
the field. As the gravitational intensity and sur-
and face contrast density are linked by the equation

u = 21rjz; v = 21rjU,

where the fz and .f, are spatial frequencies mea-

sured in cpm in the mks system. the equivalent source would also experience this
effect if placed below the level of a point or line
1 CwkckI 2 c mk as / ck 1 2 1. (3.6)
source. In the synthetic example following, the
Hence G, (p, Y, z) is asymptotically dominated by depth of the equivalent source was taken at two-
the term thirds the depth of the anomalous gravitating
sphere.
In the case of a local survey, its area1 extent
may also limit the depth of the plane Zi. If (h-z<)
is large relative to the dimensions of the survey,
.d[(uAx):! + (vAx)~] . (3.7) the coefficients sib tend to approach a common
1 value a, where

(h - zi)
a=
I: ((Xi - (.yk)’ + (y; - /3k)* + (zi - A)“] 312

= im l/(z - h)*. (3.9)

h+co

hccording to sampling theory (Blackman and Thus the matrix .4 becomes ill-conditioned and
Tukey, 1959) the maximum frequency at which its solution unreliable if the equivalent source is
Gz(u, 7’, z) can be specified from the set of values too far below the surface; that is if
g, (xi, yi, z;) is the aliasing or folding frequency
(h - Zi)! 1(Xj - (ok)’ + (yi - fik)‘) “2
given by urnax= r/Ax. In fact, in a well-designed
survey, the amplitude of G,(u, D, z) computed is too large in equation (3.9).
from the g, will approach zero at the aliasing fre- From equation (3.5), a lower limit of (h-z,) =
quency. 2Ax gives
Hence, in spectral terms, the discrete equiv-
alent source representation of gz must also (unl,,, X’lKl,, z) = 2 X lo-* G(0, 0, z) (3.10)
satisfy this requirement. From equation (3.7) sufficient to make G negligible beyond the aliasing
Gz(u, V, z) is seen to be negligible at frequency.
rt,,,, = x/Ax; Tln.1,= T/AK In a test on part of the data in the case history
(3.8)
following, it was found that a value of (h-zi)
providing (h-z)/Ax is sufficiently large. This ?SAx (Table 1) did not make the matrix A ill-
condition places an upper limit on the plane zi of conditioned, demonstrating that the limits on
the equivalent source. (h-ZJ are sufficiently broad to cover the case of
;\n effect which places a lower limit on the rough topography. Over the entire survey (h-
plan: zi can be deduced from Bullard and Cooper zi)=2.5Ax.
(1948). They noted that the downward continua- In summary the values (/z-z,) should satisfy
 42 Dampney

2.5Ax < (h - zi) < 6Ax, (3.11) forces the solution through regions in hyperspace
where the hyperellipsoid of R becomes very
where the upper bound is based empirically on the elongated. This makes the method over-relax
case history.
Downloaded 01/14/15 to 157.211.3.15. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

raising the possibility of a divergent iteration, in

Equation (3.1) was solved by using Zidarov’s which case the iteration is repeated with succes-
application of the “method of steepest descent” sively halved values of a until it does converge.
which is based on the geometric notion that the While theoretically it is possible to reduce R to
equation can be solved by minimizing R with zero, it is more expedient to reduce R to a value at
respect to mk. which it is mainly composed of the random errors
R = (g - Am)r(g - Am), (3.12) in g,. As R follows the path of steepest descent in
hyperspace, the random errors e, of gi do not
where the superscript T denotes the transpose influence the gradient vk =dR/&nk or A. In ac-
operator. cordance with normal practice (e.g. Kempthorne,
This is done by choosing a unit vector v, along 1962, p. 129) assume e, normal with mean zero
the line of maximum change R and along which and variance (T.
the initial vector m(l) is moved a maximum dis-
tance x(l) torn@). “Distance” is meant in the sense
of Ralston (1965, page 44) for a hyperspace of
dimension N, where N is the number of variables
m,. =ZgeiE+Z$,f:E (3.17)
In general k k

m(j+l) = m(j) + X(j)v(j)_ where

(3.13)
4‘
Naturally, the reduction of R to zero by m(j) f; = g: - c
would give the solution of equation (3.1). k=l

The geometry underlying this method is dis- mk(h - Zi)

cussed in the appendix together with the deriva-
tion of Zidarov’s equation (2). Dropping the
’ {(xi - ak)* + (yi - @k)’ $ (zi - h)‘) 3’2
superscript (j) we have and g”iis the value of gi stripped of its random
error ei. That is gi= g:+ ei, and therefore /i=
,fi+Ei. As (afi/&&) <0 for all i and k, then

(3.14)

where as N--+a from the property that ei has a zero

mean.
Thus, in the solution of g=Am, the error cri-
terion R must be reduced until it is less than E,
and the total estimated summed, squared noise in the

-
N
mk(h- Zi)
fi = gi - t (3.15)
k-1 ((Xi - Lyk)* f (J’i - Pk)2 + (Zi - h)2j3’*

where (x,, yi, zi) are the points of measurement measurements of g. This advantage of the tech-
and (&, Pk, h) iS the pOSitiOn Of mk. nique in reducing the presence of random errors in
For practical applications the introduction of a the solution of the mk is passed onto the Bouguer
factor a> 1 into equation (3.13), anomaly values later computed from the equiv-
alent source.
m(sti) = m(j) + &,(j) v(j) (3.16) E was worked out from
 Equivalent Source Technique 43

A’ Table 1. The numerical convergence of R

R=CfI using the steepest descent procedure
i=l
Valuezf
Downloaded 01/14/15 to 157.211.3.15. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

Itera-
= g (Ei+fY R=z
tion
cycle
Conver-
gence?
i-1
(8
(gi- gP)
asft+O while R-+0 despite the influence of the
a) Derby-Winnaleah Area (h-z;) =2.5&
random numbers pi, therefore R+xLN, (cJ2 as 3887.459 yes
f&O. 1562.7364 ; yes
Thus there is no point in reducing R below the yes
1063.3290
739.07411 : yes
value Nu2 where c is the variance of errors. The Note slow 556.31239 5 yes
process of iteration should continue until convergence 556.31232 6 yes

Region of 1 2896.7126 7 no
R < E = _17ag. (3.19) nonconvergence\ 556.31234 8 no

Nate fast yes

I3 deveioping an algorithm to solve the matrix convergence yes
[equation (3.1)] by the method of steepest de- yes
scent, the large number of coefficientsaikgenerated yes
44.441891 13 yes
cause computer storage problems. This may be ~ ~~
overcome by either storing the values tempo- b) Test Area (h-z;)=S.OAx
1197.9390 1 yes
rarily outside core-storage if sufficient locations are
yes
not available or by recomputing the coefficients 214.14105
384.77229 : yes
when required. The latter method was chosen in 118.76362 4 yes
Note slow ( 84.163410 - yes
this case. The algorithm requires N2/5000 set per convergence \ 84.163409 :, yes
iteration on an I.B.M. 360/65 where N is the
number of data points. Thus 1000 points would Region of 1 467.92798 7 no
nonconvergence 84.163406 8 yes
require three minutes per iteration. 84.163423 9 no
Table 1 gives the values of R for successive
iterations in analyzing the data discussed in the yes
convergence
Note fast 19.190139
30.586830 :: yes
Case History section below.
The algorithm may be speeded up by improving
the initial approximations to mk. The survey is
broken up into a number of interconnecting Figure 1C compares the theoretically correct
blocks of data sufficiently small so that the ma- field and the field projected from the terrain sur-
trix which solves for their equivalent source fits face values to a common datum plane by use of
within computer core storage. This vastly in- the equivalent source technique.
creases algorithm efficiency as the operation of The difference between the two fields, which
recalculating coefficients is not required. For ex- is the error of the projected field, is plotted in
ample, the equivalent source of 100 points can be Figure 1C. The maximum absolute error, which
calculated in less than a second per iteration on the is seen to occur over the sphere, represents a rela-
360/65. The values of the mk found for the indi- tive error of only three percent.
vidual blocks can then be put together as the
initial values of the equivalent source for the
entire survey.
SYNTHETIC TEST

The technique was tested on synthetic data.

The gravity field of a sphere at a depth of 100
units was computed at points 25 units apart over
an idealized terrain representing a valley and
plateau separated by a 25 unit high cliff (Fig-
FIG. 1A. The synthetic test of an anomalous
ure 1A). This field is contoured in Figure 1B. gravitating body and terrain.
 44 Dam pney

50 100 150 200 250 300 phism. It was extensively intruded and folded by
I I I I I I granite to form hills which rose out of the pre-
5OC
Tertiary landscape. A long period of erosion which
Downloaded 01/14/15 to 157.211.3.15. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

followed the granite intrusions, cut deeply into the

topography. Veins in the granite were leached of
tin which was subsequently deposited as low
density (2.0 gm/cm2) alluvium along depressions
in the plane levelled out of the Mathinna Sand-
stone.
\ The topography at the beginning of the Ter-
tiary had thus been formed by the eroding in-
fluence of the Ringarooma River. The river’s
path was controlled by the resistant hills of
granite and sandstone, and by a relative depres-
250 sion in the land surface to the nearby sea which
lowered the river’s slope, changing it to a series of
small lakes and estuaries.
2oc
Then, during the Tertiary, this quiescent scene
was disrupted by the outpourings of basaltic
150 lava flows. The Ringarooma River channel was

IO0

5C

FIG. 1B. The Bouguer anomaly measured over the

terrain of Figure 1A. Contour valuesin mgal.

CASE HISTORY-THE DERBY W’INNALEAH

GRAVITY SURVEY

The equivalent source technique was tested on

data collected by Howland-Rose (1966) on behalf
of the Australian Commonwealth Bureau of
Mineral Resources at Derby, Tasmania, Aus-
tralia. The location of the survey is shown in
Figure 2.
The survey of about 4 sq mi was chosen for this
study as its area is extremely precipitous due to
the influence of the Ringarooma River on the
topography.
The geology of the Ringarooma Valley, which
includes the survey area, has been described by
Nye (1925) (Figure 3).
The oldest rock in the area is the Silurian
Mathinna Group composed of slates and sand-
stones. During the course of geological time this
group has been extensively folded and faulted. FIG. 1C. The absoluteerrorsof the Bougueranomaly
During the Devonian, the Mathinna Group values projected from Figure 1B using the equivalent
sourcetechnique.Dashed contoursrepresentnegative
suffered its most extensive period of diastro- errors.Contour interval= 0.01 mgal.
 Equivalent Source Technique 45

BASS STRAIT
Downloaded 01/14/15 to 157.211.3.15. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

SCALE IN M\LES

AREA COVERED 6Y
1964 SURVEY

REFERENCE TO AUSTRALIA t250@0

STANBARD MAP SERIES

FIG. 2. Derby-Winnaleah area, Tasmania, locality map.

buried by lava and its course changed at Derby manner following Hammer (1963) so that the
as it was forced along the granite basaltic lava final topographically corrected Bouguer anoma-
contact in the survey area. The ancient river that lies took into account the simple model of the
had previously flowed west of the Mt. Cameron geology shown (see cross section) down to 150
Range (Figure 2) was diverted to the east. m above mean sea level (“base” level in Figure
The modern Ringarooma River has eroded its 4). Bott’s (1961) method for calculating the topo-
way down the gorge now containing the Briseis graphic effect was employed. The topographic
Mine (Figure 3) in the survey area. A rough correction was defined as in Grant and West
north-south cross section is shown in Figure 4 of a (1965, p. 239) and hence, following Bott, the
simplified interpretation of the area’s geology. various densities of the rocks were used to com-
The gravity observations were reduced in a pute the topographic correction. The Bouguer
 46
Downloaded 01/14/15 to 157.211.3.15. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

- Elevationcontous AMSL Eevonian Granite 2.6g/cm3 0

80 II II
__-- l&W&tnaC&oup 26g/cm3 m
(inferredfrom air photos)
- Ringorooma River Tertiary Gravek ZOg/cm’O

FIG. 3. Geology and topography of Derby. Elevation values are in units of 100 ft.
Northings and Eastings are in units of 1000 yards.

HILL

PLATEAU
------- 250 METERS
“Y”“YYY”YVY”V
“,“““Y”V” “Y” YYYY ““YV”
Ovvvv BASALTvvv vvv vvv vv
vvYvvYvvv”vvvv V”VVV”V
- ---- yv v v v V_V v-v v-v _v_ S&EgVEL v v 150 METERS
ALLUVIUM
- -~- 3 MILES .-. L

FIG. 4. Diagrammatic cross section through plateau and valley, north to south. All heights
are with reference to mean sea level.
 Equivalent SOI~rce Technique 47

anomalies are shown in Figure 5. As absolute the precipitous topography and the associated
Bouguer anomaly values are not required in a difficulty of making exact topographic correc-
local survey, the zero contour has been set so that tions. The plane of the equivalent source was
Downloaded 01/14/15 to 157.211.3.15. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

it passes through the center of the map. taken at a height equal to mean sea level which
The equivalent source of this Bouguer anomaly satisfies the previously discussed limits.
map was found. The noise level parameter E was The Bouguer anomalies at regular grid points
worked out from equation (3.19) to be 50 (mgal)* on a horizontal plane about the same height as the
for the 860 stations involved. The random error of basalt plateau (250 m = 800 ft) were then com-
each Bouguer anomaly value was taken to be the puted from the equivalent source and are plotted
order of 0.25 mgal. This is reasonable in view of in Figure 6. The grid itself is not shown, but it has

l:IG. 5. TopographicallycorrectedBougueranomaliesfor the Derby-Winnaleah area.

Contour interval = 1.0 mgal.
 48 Dampney
Downloaded 01/14/15 to 157.211.3.15. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

FIG. 6. The topographicallycorrected Bouguer anomalies at 250 m above mean sea level.
Contour interval = 0.5 mgal.

a mesh interval of 100 m and extends over the as the plateau values, allows a direct comparison
area contoured. Comparison of Figures 5 and 6 of the Bouguer anomalies across the 300 ft eleva-
shows that the technique’s projection of the orig- tion difference between the two regions. Evenly
inal Bouguer anomaly values has eliminated spaced grid points also allow more objective con-
minor random fluctuations present in the con- touring than do long survey lines.
tours of Figure 5. Projection of the values in the The Bouguer anomalies are now in a satisfac-
valley region (delimited to the North by the -3.0 tory form for processing by techniques requiring
contour in Figure 5) to the same reference plane gridded data.
 Equivalent Source Technique 49

However, the simplicity of directly calculating major negative trend across the map apparently
the field from the equivalent source suggests us- reveals the low density alluvium marking the
ing the technique for vertical continuation. Figure buried channel of the ancient Ringarooma River.
7 shows the dominant regional trend of the
Downloaded 01/14/15 to 157.211.3.15. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

The positive anomaly in the southeast corner is

Bouguer anomalies, computed in this way, on a attributed to the water filled Briseis mine work-
plane 500 m above sea level. ings taken as being empty in the topographic
The final residual map, Figure 8, was found by corrections. The two positive anomalies in the
subtracting the “regional” values (defined by northeast and northwest corners, corresponding
Figure 7) from the values at 250 m (Figure 6). The to Mathinna sandstone outcrops, mark the limits

\
\
\ \/

FIG. 7. The topographicallycorrectedBougueranomaliesat 500 m above mean sea level.

Contour interval=O.S mgal.
 50 Dampney
Downloaded 01/14/15 to 157.211.3.15. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

FIG. 8. The topographicallycorrectedresidualBougueranomaliesat 250 m above mean sea level.

Contour interval=0.2 mgal.

of the alluvial plane imposed by these resistant zontal reference plane. The technique has the im-
pre-Tertiary hills of sandstone and granite. portant application of objective preparation of
Bouguer anomalies for processing by methods
CONCLUSION
mentioned in the review. It should also be possible
The equivalent source technique as demon- easily to extend this technique to other potential
strated in Synthetic Test section and the Case fields, particularly magnetic fields.
History offers a convenient and accurate way to However, limitations have to be realized. Pro-
interpolate gravity data onto a grid. It can be jecting gravity data onto a horizontal plane in-
used to make the final correction to the Bouguer volves vertical continuation and so a large hori-
.anomaly by projecting measurements onto a hori- zontal coverage of field values may be required.
 Equivalent Source Technique 51

This can be seen from the lateral extent of the - 1966b, Geophysical studies in Tasmania; Part
A, Interpretation of the gravitational potential field
theoretically perfect coefficient set for vertically using the frequency domain: M.Sc. Thesis (unpubl.),
continuing a potential field at a height equal to Dept. of Geology, Univ. of Tasmania.
Downloaded 01/14/15 to 157.211.3.15. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

the station spacing. (Dampney, 1966a). Danes, 2. F., 1961, Structure calculations from gravity
data and density logs: Mining Trans, v. 223, p. 23-29.
This technique offers a low-error-level, objec- Grant. F. S.. and West. G. F.. 1965. interpretation
tive approach to three-dimensional interpolation theory in applied geophysics: ‘New York, McGraw-
and the related problem of contouring data. As Hill Book Co., Inc.
Hammer, S., 1963, Deep gravity interpretation by
an added bonus it is very economical to use the strippmg: Geophysics,v. 28, p. 369-378.
equivalent source (once calculated) for reason- Henderson,R. G., 1960,Comprehensivesystemof auto-
ably accurate vertical continuation. matic computationin magneticand gravity measure-
ments: Geophysics, v. 25, p. X4-585.
ACKNOWLEDGEMENTS Howland-Rose, A., 1966, Derby-Winnaleah Survey,
Tasmania 1964: Dept. of National Development,
This research formed part of an M.Sc thesis Bureau of Mineral Resources, Geology and Geo-
program carried out in the Department of Geol- physics, Record No. 1966/10.
ogy, University of Tasmania (Dampney, 1966b). Kempthorne, O., 1962, Design and analysis of experi-
ment: New York, John Wiley and Son.
I gratefully acknowledge the help and inspiring Naudy, H., and Neumann, R., 1965, Sur la dL:finition de
supervision of Dr. R. Green. I’anomalie de Bouguer et ses consequencespratiques:
Sincere thanks are extended to Dr. G. F. West, Geophys. Prosp., v. 13, p. l-11.
Neidell! N. S., 1966, Spectral studies of marine geo-
Geophysics Laboratories, University of Toronto, physical profiles: Geophysics, v. 31, p. 122-134.
for his useful comments made during the prepara- Nye, P., 1925, The sub-basaltic tin deposits of the
tion of this manuscript. Ringarooma Valley: Tas. Dept. of -Mines, Geological
Survey Bull. No. 35.
The author is also indebted to Mr. J. Boothroyd
Ralston, A., 1965, A first course on numerical analysis:
for programming advice and to the Director, New York, McGraw-Hill Book Co. Inc.
Australian Commonwealth Bureau of Mineral Roy, A., 1962, Ambiguity in geophysical interpreta-
tion: Geophysics, v. 27, p. 90-99.
Resources for permission to publish the Derby
Sax, R. L., 1966, Applications of filter theory and infor-
Winnaleah Case History. mation theory to the interpretation of gravity mea-
surements: Geophysics, v. 31, p. 570-575.
REFERENCES
Strakhov, V. N., and Devitsyn, V. M., 1965, The reduc-
Bhattacharyya, B. K., 1965, Two-dimensional har- tion of observed values of a potentral field to values at
monic analysis as a tool to magnetic interpretation: a constant level: Bull. Acad. Sci. USSR., Geophys.
Geophysics, v. 30, p. 829-857. Ser. (English Transl.), p. 250-255.
- 1966, Continuous spectrum of the total magnet- Zidarov, D., 1960, Determination des champs de gravi-
ic field anomaly due to a rectangular prismatic body: tation (ou magnetiques) locaux et regionaux: Comp-
Geophysics, v. 31, p. 97-121. tes rendus de 1’AcadCmie bulgare des Sciences, v. 13,
Blackman, R. B., and Tukey, J. W., 19.59,The mea- p. 531-534.
surementof power spectra: New York, Dover Pub- - 1965, Solutions of some inverse problems of ap-
lications,Inc. plied geophysics: Geophys. Prosp., v. 13, p. 240-246.
Bott, M. H. P., 1961,The useof electronicdigital com-
puters for the evaluation of gravimetric terrain cor- APPENDIX
rections:Geophys.Prosp., v. 7, p. 45-54.
___ 1967!Solutionof the linear inverse problem in Derivation of equation (3.12)
magneticinterpretation with applicationsto oceanic
magneticanomalies:Geophys J. R. Astr. Sot., v. 12,
p. 151-162. Minimize R = (g - AmjT(g - Am) (1)
Bullard, E. C., and Cooper, R. I. B., 1948, The deter-
mination of the masses necessarv to oroduce a eiven to find the solution of m in equation (3.12). The
gravitat :ional field: Proc. Roy. Sot. (London) Se;. A.,
194, p. 332-347. direction of maximum change of R with respect to
D& pney, C. N. G., 1966a, Three criteria for the judge- rni is given by dR/dm, the ith component of VR
ment of vertical continuation and derivative methods
with respect to m.
of geophysical interpretation: Geoexploration, v. 4. p.
3-24. Thus following Zidarov (1965)

R= 5
i=l
(gi(Mk)- gY(mJ f, (2)
(h - zi)
(3)
 52 Dampney

gi(Mk) is the measured field due to the true masses Developing 4(h) as a power series in X from
Mk, and gij (mk) is the jth approximation of Taylor’s theorem, and taking into consideration
g, from the approximate masses mhj’at (CQ,&, h). only the first two terms (as d3R/ami=0 for all k
Downloaded 01/14/15 to 157.211.3.15. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

The positions of the rnk are not restricted, but from equation 4), we obtain
for convenience place all the mk at the same
dR dR
height. In the case history, each discrete mass was +(A) = R - X dm ;;
positioned vertically below each of the N data
points making up the survey. However, any other
arrangement is valid within restrictions imposed
(8)
by equation (3.12) and the nature of the matrix
A. In the solution we assume nothing about the d?R
elements of A except that they are real and that __ zz
A is square. dm2
Therefore,
Kow d2j/dm2=0 from equation (a), therefore

f#a(X)= f? - 4Xf2
where
df 4
S
+ 4X2f? _ (10)
k = gk(Mi) - c i=l
( dm ) ’
+(A) will have a maximum when d+/dX=O.
mi(h - zk) Thus
.___-
I (xk - (Yj)’ + (m - pi)’ + (Zk - h)?J3’2 ’
that is
4fz(E)2- Shfl($)” = 0. (11)
S
R = zj;. Hence
(3
z=l

We find following equation (3.11) that

(i) (j-1,
-_x
(j-1)
-__
dR
mk = ‘?& (6)
arn(i-l)
k

where m&j) is a closer approximation to the true

value of mk than rnt-‘).
Therefore,
(12)
aR
R = R m;j’ _ A(‘) -, .’‘
dm(i) ’
which is Zidarov’s (1965) equation (2) with ??I!+
substituted for 8jk and R for LT.
This method assumes R has only one minimum.
The condition for one minimum that d3R/dm3=0
is less than R=R(nz(,“, . , mf$‘).
is seen to be true.
To reduce R as quickly as possible, we maxi-
mize + (A) with respect to X: Geometric considerations

(j) aR Ralston (1965, p. 439-442) derives an equation

+(A) = R&j) - X for X for the case of a symmetric matrix. He also
yjq .‘.)
discusses the geometry of the method of steepest
(I 1 8R descent.
m.v - p_ .
(7) The geometry for the unsymmetric matrix is
am,:: )>
similar. For ideal convergence the magnitude of
 Equivalent Source Technique
Downloaded 01/14/15 to 157.211.3.15. Redistribution subject to SEG license or copyright; see Terms of Use at http://library.seg.org/

FIG. 9. The path of the steepestdescentsolution.

the successive iterations X(l), X(*) . . from the and shallow and the effect described by Bullard
first approximation is shown in Figure 9. and Cooper (1943) discussed in section 3 occurs.
As the solution will work its way towards
c=g-Am
elongated regions of the hyperellipsoid, the in-
in the figure and t represents the major and fluence of small (relative to the station spacing)
minor axes of the hyperellipsoid R. shallow sources will be the last to be extracted
Physically one sees in the solution of g=Am from R. Hence as R is not reduced to zero, the
that if the true value of WZ~is greatly different shallow sources will be treated as noise in the
from the other m’s, the hyperellipsoid is very data. As small shallow sources arc by definition
elongated along the mk axis. This may result if inadequately sampled, their elimination as noise
mk is in a position where the actual source is small is in accordance with normal practice.