508
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING,
VOL 32.
NO
3,
MAY
1994
Fast Reduction
of
Potential Fields Measured Over an
Uneven Surface to a
Plane
Surface
P.
S.
Naidu and M.
P.
Mathew
Abstract-The present work is aimed at rapid reduction of
the gravity and magnetic fields observed over an uneven sur-
face to a horizontal plane. The approach suggested here is to
estimate the Fourier transform
of
the potential field over an
imaginary horizontal plane lying entirely above the ground sur-
face and impose boundary conditions; namely, the solution must
satisfy the observed field over the ground surface and vanish
over an infinite hemisphere. The desired Fourier transform is
obtained in an iterating manner.
A
2-D
FFT algorithm can con-
siderably reduce the computational burden. The FFT approach
cannot be used unless the discrete data is available on a rectan-
gular grid. If the observations are scattered, interpolation to
the nearest grid point will have to be carried out. Interpolation
introduces marginal increase in the rms error. The iterating
approach is about
10
times faster than the least squares ap-
proach.
I.
INTRODUCTION
N
geophysical surveys, and in particular, in gravity and
I
magnetic surveys, the observation stations are often
controlled by topography, accessibility, and the distance
to the base station to which one has to return for frequent
calibration. These factors naturally control the distribu-
tion of the observation stations; mostly along the existing
motorable tracks. In aerial low-altitude mineral surveys
the aircraft is often flown at a constant height (e.g.,
50
m) above the topographic surface. Consequently, the
flight lines, specially in hilly terrain, are undulating in the
vertical plane, though practically parallel in the horizontal
plane. Thus, the observation stations may be considered
as unevenly distributed on an irregular topographic sur-
face or on flight lines lying on an irregular surface parallel
to the topographic surface below. Such observations of
the gravity and magnetic field are inadequate for modern
quantitative methods that require discrete observations
over a square grid on a horizontal plane. The earliest work
on the subject of the continuation of the potential field
measured on a curved line to a horizontal line above the
curve line was by Strakhov and Devitsyn
[l],
who used
the method of successive approximations for the solution
of an integral equation of first order. Their method works
well when the surface undulations are weak. Tsirul’skiy
[2] reduced the problem to solving of an integral equation
Manuscript received November
3
1,
1993.
P.
S.
Naidu is with the Department
of
Electrical Communication Engi-
M.
P.
Mathew is with Airborne Mineral Survey and Exploration, Geo-
IEEE
Log Number
9400439.
neering, Indian Institute
of
Science, Bangalore
560012,
India.
logical Survey
of
India, Bangalore
560001,
India.
of the second kind in the frequency domain. An algorithm
to solve the integral equation by a method of successive
approximation is given in
[3].
Bhattacharya and Chen
[4]
have formulated the problem in terms of the Fredholm
integral equation and gave a series solution in the space
domain
[5].
Hansen and Miyazaki
[6]
modified the method
given in
[4]
for better accuracy. Henderson and Cordell
[7]
have used finite harmonic series expansion of the ob-
served data and then, using the estimated coefficients, they
were able to reduce the field to a plane surface. The
French workers used the Backus and Gilbert inverse prob-
lem formalism for the continuation of the potential field
from an uneven surface to a plane surface
[8],
[9].
Parker
and Klitford
[
101 have used the Schwartz-Christoffel
transformation to map an uneven track of a bottom-towed
magnetometer to a straight line.
The present paper addresses itself to the problem of re-
ducing the field observed over an irregular surface onto a
horizontal plane above the surface of observation. The ap-
proach adopted here is similar to that in
[l]
and [2], but
it
is
valid both for two- and three-dimensional problems
and it uses a fast computational algorithm. The compu-
tations require the repeated application of a
2-D
fast Fou-
rier transform, thus making the task of reduction very fast
compared
to
the space-domain methods described in
[4]-
[lo]
and currently used in industry. However, the major
limitation of the present approach lies in the requirement
that the
x
and
y
coordinates of a measurement station must
lie on a square grid, though the
z
coordinate may take any
arbitrary value (see Fig.
1).
Only then can one exploit the
high speed of the FFT algorithm. From a practical point
of view this does not pose a serious problem as it is a
common practice to generate gridded data from the scat-
tered observations. Using such an approach we have ob-
tained accurate reduction even when the observation sta-
tions were scattered, but a sufficiently large number of
them were available. The present paper is divided into
five sections, In the next section we briefly review the
least-squares inversion method and study its performance
when the available data is randomly scattered. In Section
I11 we describe the theory of an iterating FFT-based
method and in Section IV we provide a few illustrations
as well as a synthetic example. In Section V we provide
a practical example. We shall also compare the accuracy
and speed of the present method with that of the least-
squares method.
0196/2892/94$04.00
0
1994
IEEE
NAIDU
AND
MATHEW:
REDUCTION
OF
FIELDS
OVER
UNEVEN
TO
Plane of
reduction
A
Number
of
2500 2000
observations.
rms
error
(%)
1.29 1.49
CPU
Time (min)
125
105
on
Vax
750.
PLANE
1500
1000
750
600 450
1.78 1.81 3.20 6.10 127
89 62
50
42 33
Fig.
1.
Profile
of
observation stations over
a
hilly terrain.
The
observed
field
is
now reduced to a horizontal plane just above the highest point in
the area
of
survey.
11.
LEAST-SQUARES
INVERSION
We shall assume that all observation points are scat-
tered in a three-dimensional space but that they all lie on
a smooth plane representing a topographic surface. The
data points when projected onto a horizontal plane are
scattered over the plane. Let
f(x,
y,
H)
be the potential
field on a horizontal plane
H
units above the mean obser-
vation surface and
F(u,
v,
H)
be its
2-D
Fourier trans-
form (see Fig.
1).
We shall continue the unknown field
on the horizontal plane downward to every point of ob-
servation,
f(xi,
yi,
H
+
Azi)
.
lrca
nm
d~
dv,
i
=
0,
1,
-
*
*
N
-
1.
(1)
.
e
jWi
+
Uyi)
In above equation the unknown quantity is
F(u,
v,
H),
which may
be
discretized over a finite domain with the
assumption that it vanishes outside this domain. For the
sake of simplicity we shall assume the domain is rectan-
gular,
2a
X
2b,
and it is divided into
(2p
+
1)
x
(2q
+
1)
cells. We shall discretize the integral in
(1)
and write
f@i,
yi,
H
-
Azi)
lP
=-
C
C
F
PQ
k=
-p
I=
-4
where
P
=
2p
+
1
and
Q
=
2q
+
1.
Define following
matrices:
f
=
CO~
{f(~;,
yi,
H
-
Azi),
i
=
1,
2,
3
*
N}
Size:
N
x
1
1
k= -p, -p
+
1
**'p
-
1,p
l=-q,-q+l*..
9-
17q
Size:
N
x
(2p
+
1)(2q
+
1)
SURFACE
509
Table
I
and
k=-p,-p+l.**p-l,p
I=
-9,
-q+
1
**.
4-
1,q
Size:
(2p
+
1)(2q
+
1)
Equation
(2)
in matrix version becomes
f
=
AF.
We
should be able to solve the above system
of
linear equa-
tions for the unknowns contained in vector
F.
Naturally,
we must have
N
1
(2p
+
1)
(2q
+
1);
in the presence of
noise a reliable solution is obtained only when
N
>>
(2p
+
1)
(2q
+
1).
A
least-squares solution of
(2)
is given by
F
=
(AHA)-'AHf
(3)
where the superscript
H
stands for the hermitian trans-
pose. We have evaluated
(3)
on a Vax
750.
The compu-
tation time required for different data sizes is shown in
Table
I.
A.
Numerical
Example
A ground magnetic survey over an area of
140
X
140
units (arbitrary units) over an undulating terrain was sim-
ulated. The field at randomly scattered stations uniformly
covering the entire area was calculated. The magnetic field
is caused by four semiinfinite prisms. A map of the un-
dulating surface including the stations and an outline of
magnetic targets are shown in Fig.
2.
The undulating sur-
face was generated mathematically using the equation
1
X2
[Jx2
+
y2
+
16
A.z(x,
y)
=
4.0
+
2.0
COS
Jx'
+
y2
+
16
y2
+
4
+
sin
-
The height of the observation stations varied from
1.1
to
7.0
units. First, the scattered field measurements were
used for reduction. The least-squares inversion approach
was used to reconstruct the field at a height of
7.5
units,
just above the highest point. The mean-square error be-
tween the reconstructed field and the actual (computed)
field is tabulated in Table I.
The error is expressed in percentage relative to the peak-
to-trough anomaly amplitude of
65.5
gamma. The
rms
er-
ror is found to increase rapidly when the number of sta-
tions is reduced below
600.
I
510
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING,
VOL. 32,
NO.
3,
MAY
1994
Fig.
2.
A
map
of
the undulating surface including the observation stations
(+)
and an outline
of
magnetic targets.
111.
FFT-BASED METHOD
Letfo(x,
y,
Az)
be the potential field (gravity or mag-
netic) observed over a known but arbitrary surface
Az(x,
y).
We shall assume that the
x
and
y
coordinates of each
observation station lie on a regular grid; however, the
z
coordinate is arbitrary. Where this is not true we shall
assume that through interpolation, gridded digital mag-
netic data has been prepared. The data points, when
pro-
jected onto a horizontal plane, fall on a regular grid of
points. Only the height of a data point is irregular. Let
f(x,
y,
h)
be the desired field on a horizontal plane Hunits
above the mean observation surface
[H
>
Az
(x,
y)].
Now
we shall continue the field downward to the observation
surface and require that the downward continued field
must satisfy the observed field
(4)
.
ej(u+vY)
du dv
where
F(u,
u,
h)
is the Fourier transform of the desired
fieldf(x,
y,
h)
and
s
=
v.
Note that the repre-
sentation of the potential field given in
(4)
satisfies the
Laplace equation and the vanishing boundary condition at
+
00.
Next, we express
Az
(x,
y)lh
=
€17
(x,
y)
where
E
<
1
and
r]
(x,
y)
is a normalized surface, that is,
1
xy
-
c
c
$(x,
y)
=
1,
x,
Y
-+
00
4XY
-x
-Y
Use the following expansions:
k=O
k!
IXI
F(u,
u,
h)
=
c
eiFi(u,
u,
h)
i=O
in
(4).
Express the right-hand side
of
(4)
as a polynomial
in
E.
Then,
(4)
becomes
1
h(x,
Y,
Az)
=
47~2
1:
A@,
u,
r](x,
Y))
.
esh
e
l(u
+
VY)
du dv
(5)
r]
(x,
Y)
*s2FI(~,u,h)
+-
1
*
sF2(u,
u,
h)
-
F3(u,
u,
h)
I
+
-
*
-
.
We set the coefficients
of
E
to zero giving us the following
set
of
interconnected equations from which we can obtain
Fk(u,
u,
h),
k
=
0,
1,
2,
-
.
For example,
Fo(u,
v,
h)
=
e-sh
c
XY
cfo(x,
y,
Az)e-i(ux+vy)
(.+
PT
Fo(u,
u,
h)sesheei("x+vy)du du
s-
s,
(6)
and
so
on. The above system
of
iterating equations is eas-
ily evaluated using the FFT algorithm. The convergence
is rapid, usually requiring about
10
iterations, depending
upon the amplitude of the surface undulations.
IV.
COMPUTER SIMULATION
For the purposes of testing the new method of reduction
to
a horizontal surface we have a synthetically generated
field along a set
of
16
nearly parallel lines placed on the
NAIDU AND
MATHEW:
REDUCTION
OF
FIELDS
OVER
UNEVEN
TO
PLANE
SURFACE
511
Fig.
3.
The contoured maps
of
measured field are shown in (a) (parallel
flight lines) and
(b)
(scattered stations). The results
of
field reduction to a
horizontal plane at
7.5
units above using the new method are shown in (c)
and (d). Compare the reduced field with the theoretical field shown in Fig.
4
Fig.
4.
The theoretically computed field on the plane
of
reduction.
undulating surface, as well as at randomly scattered
(2500)
points (see Fig.
2).
The parallel lines would correspond
to undulating flight paths of an aircraft over a hilly terrain.
On each line
160
samples were collected, giving a total
of
2560
samples. The data thus created was gridded using
a commercially available software package (e.g., Data
Plotting Services Inc., Toronto, Ont., Canada). The con-
toured maps are shown in Fig. 3(a) and (b). The result of
field reduction to a horizontal plane at
7.5
units above
using the new method is shown in Fig. 3(c) and (d). The
theoretically computed field on the plane of reduction is
shown in Fig.
4.
Table I1 gives the
rms
difference be-
?.I'
.d
??*
*.'
Fig.
5.
Aeromagnetic field over a hilly terrain flown at
low
altitude (nom-
inal height
50
m above ground). Contour interval is
207.
Table
I1
Parallel
Profiles
Grid data 1.44
tween the reduced and the theoretical field. We note that
in comparison with the least-squares methods there
is
a
marginal increase in the
rms
error; however, there is a
considerable savings in computer time, by a factor of
10
to
12.
The increase in the
rms
error
is
largely due to the
interpolation error. In support of this statement we have
carried out the reduction to a horizontal plane using data
on a regular grid requiring no gridding. From the last row
in Table
I1
we note that the
rms
error has come down and
is
of
the same magnitude as in the least-squares approach.
V.
APPLICATION
TO
REAL
DATA
We have applied our new method of reduction to a real
data set obtained in a low altitude (nominal flight height
50
m) aeromagnetic survey carried over a part of the dia-
mond-bearing rocks, consisting of granites and granite
gneisses belonging to the south Indian shield. The survey
was carried out by the
AMSE
division
of
the Geological
Survey
of
India using the Scintrex system mounted on a
Twin Otter aircraft. In addition to the total magnetic field,
the system provided continuous measurements of baro-
metric height and air column thickness, which were used
I
I
512
IEEE TRANSACTIONS
ON
GEOSCIENCE AND REMOTE SENSING,
VOL.
32,
NO.
3,
MAY
1994
*
Id
50‘
Fig. 6. Reconstructed flight surface from the barometric data along with
the flight paths. Contour interval is 20 m.
Fig. 7. Total field reduced to a horizontal plane at a height of 678 m above
MSL, The plane of reduction
is
just above the highest point in the survey
area.
to reconstruct the actual flight surface. The navigation
system consisted of Doppler as well as a visual method
based on aerial photostrips. It is believed that the posi-
tional accuracy is better than
5
m. The total measured
field after the application of the diurnal and instrumental
drift corrections and subtracting a base level of
41
loOy
is
shown in Fig.
5
and the reconstructed flight surface is
shown in Fig.
6.
On the flight surface we have also su-
perposed actual flight lines. The contour values are in me-
ters, representing the height above
MSL
plus
297
m. In
Fig.
7
we show the reduced field, which now corresponds
to what would have been observed on a horizontal plane
at a height of
678
m above
MSL
or
381
m above the local
base level.
ACKNOWLEDGMENT
We wish to thank the Director-General
of
the Geolog-
ical Survey
of
India, Calcutta, for permission to use their
aeromagnetic data. We also wish to thank
M.
R. Nair,
Deputy DG,
AMSE,
for his interest in this work and
M.
V.
Ramanmurthy, R. Sharma, Dr.
S.
V.
Anand,
P.
V.
Ani1 Kumar,
M.
Singh, and
B.
Singh, all from
AMSE,
for their help in preparing the maps.
REFERENCES
[l]
V. N. Strakhov and V. M. Devitsyn, “The reduction of observed
values of a potential field to values at a constant level,”
Bull.
(Izv.)
Acud.
Sci.
(I.S.S.R.,
Earth Physics, no. 4, pp. 256-261, 1965.
[2] A. V. Tsirul’skiy, “On the reduction of observed values of potential
fields to a single level,”
Bull.
(Izv.)
Acud.
Sci.
U.S.S.
R.,
Earth Phys-
ics, no.
3,
1968.
[3] A. V. Tsirul’skiy and L. Ya. Ospishcheva, “An algorithm for the
reduction of observed potential field values to a common level,”
Bull.
(Izv.)
Acud.
Sci.
U.S.S.R.,
Earth Physics, no. 4, pp. 266-270, 1968.
[4] B.
K.
Bhattacharya and
K.
C. Chan, “Reduction of magnetic and
gravity data on an arbitrary surface acquired in a region of high to-
pographic relief,”
Geophys.,
vol. 42, pp. 141 1-1430, 1977.
[5] J. C. Wynn and B. K. Bhattacharya, “Reduction of terrain induced
aeromagnetic anomalies by parallel surface continuation: A case his-
tory from southern San Juan mountains of Colorado,”
Geophys.,
vol.
42, pp. 1431-1449, 1977.
[6]
R.
0.
Hansen and Y. Miyazaki, “Continuation of potential fields be-
tween arbitrary surfaces,”
Geophys.,
vol. 49, pp. 787-795, 1984.
[7] R.
G.
Henderson and L. Cordell, “Reduction of unevenly spaced
potential field data to a horizontal plane by means of finite harmonic
series,”
Geophys.,
vol. 36, pp. 856-866, 1971.
[8]
V. Courtillot,
J.
Ducruix, and
J.
L.
Le Mouel, “Le prolongement
d’un champ de potentiel d’un contour quelconque sur un contour hor-
izontal: Une application del la method de Backus et Gilbert,”
Ann.
Geophys.,
vol. 29, pp. 361-366, 1973.
[9]
J.
Ducruix,
J.
L. Le Mouel, and V. Courtillot, “Continuation of three
dimensional potential fields measured on an uneven surface,”
Geo-
phys.
J.
R.
Asrr.
Soc.,
vol. 38, pp. 299-314, 1974.
[lo]
R.
L.
Parker and
K.
D. Klitgord, “Magnetic upward continuation
from an uneven track,”
Geophys.,
vol. 37, pp. 662-668, 1972.
P.
S.
Naidu
was born in 1937. He received the
M.S.
degree from the Indian Institute of Technol-
ogy, Kharagpur, and the Ph.D. degree from the
University of British Columbia, Vancouver, Can-
ada.
Currently, he is a Professor of Electrical Com-
munication Engineering at the Indian Institute
of
Science, Bangalore. He has published over sev-
enty research papers in geophysical and under-
water signal processing.
Dr.
Naidu was a reciuient
of
a Humboldt Fel-
lowship (1979) and an NRC Senior Associate (1988).
M.
P.
Mathew
was born in 1947. He received the
M.S. degree in mathematics and statistics from the
Kerala University.
He joined the Geological Survey of India in
1970 and performed gravity and magnetic
sur-
veys. Currently, he is with the Airborne Mineral
Survey and Exploration wing of the
GSI
where he
is engaged in aeromagnetic data processing.