Ren~ FORSBERG

Geodetic Institute
Gamlehave Alle 22
DK-2920 Charlottenlund, Denmark

G R A V I T Y FIELD T E R R A I N EFFECT COMPUTATIONS BY FFT

Abstract

The widespread availability of detailed gridded topographic and bathymetric

data for many areas of the earth has resulted in a need for efficient terrain effect
computation techniques, especially for applications in gravity fieM modelling. Compared
to conventional integration techniques, Fourier transform methods provide extremely
efficient computations due to the speed o f the Fast Fourier Transform (FFT). The
Fourier techniques rely on linearization and series expansions of the basically unlinear
terrain effect integrals, typically involving transformation of the heights/dep ths and their
squares. The FFT methods will especially be suited for terrain reduction o f land gravity
data and satellite altimetry geoid datat
In the paper the basic formulas will be outlined, and special emphasis will be
put on the practial implementation, where a special coarse~detailed grid pair
formulation must be used in order to minimize the unavoidable edge effects o f FFT, and
the special properties of FFT are utilized to limit the actual number o f data
transformations needed. Actual results are presented for gravity and geoid terrain effects
in test areas of the USA, Greenland and the North Atlantic. The results are evaluated
against a conventional integration program : thus, e.g., in an area of East Greenland
(with terrain corrections up to lOmgal), the accuracy of FFT-computed terrain
corrections in actual gravity stations showed an r.m.s, error o f 0.25 mgal, using height
data from a detailed photogrammetric digital terrain model. Similarly, isostatic ocean
geoid effects in the Faeroe Islands region were found to be computed with r.m.s, errors
around 0.03 m .

Introduction
Gravitational effects of the topography and bathymetry are the dominant
sources of local gravity field variations in many areas of the earth. By computationally
removing the effects of topography, bathymetry and possibly the associated isostatic
compensation, a more smooth residual field will usually be obtained, more suitable for
gravity field modelling and enhancing effects of other density anomalies in the earth's
crust and upper mantle.

Bull. G6od. 39 (1985) pp. 3 4 2 - 3 6 0 .

342
 R. F O R S B E R G

In the present context, the term "terrain effect" is used to express the
gravitational effects on gravity anomalies, deflections of the vertical, etc., from a mass
body representing topographic/isostatic masses in part or in whole. Examples of possible
terrain effects are "topographic effects," the direct influence of the visible topography
in mountainous areas, "isostatic effects," where the effect of a hypothesized isostatic
compensation is also taken into account, and finally, the "residual terrain model (RTM)
effects," where only local high-frequent topographic irregularities are taken into account
by referring all elevations to a smooth mean elevation surface, defined, e.g., through a
high-degree and order spherical harmonic expansion of topographic heights, cf. Fig. 1D.
For more details and application examples, see Forsberg (1984) or Forsberg and
Tscherning (1981 ).

A} T O # ~ q A , ~ I Y ,) r e m , u . co 'r,o.

c) ISOSTA TIC D) RESIDUAL TERRAIN MODEL

t
D "3;r lem
mean e l e ~ t ~
~ce

1 ,, l,
,,,, ,., , ~ ~ ~ ~ l i ~' ~,-',', ~ . - - 0 . d l

",l l ~ ~ l ~ l' ,,l 'l t I l ~ "

%1 l l ~ la *',b . o

Fig. 1 - Examples of density distributions associated with various types of terrain

effects.
A : topographic effect, Le., "refined" Bouguer reduction (consisting o f
a Bouguer plate effect minus the terrain correction "B"),
C : conventional Airy-isostasy,
D : residual terrain model (R TM), mean elevation surface, given, e.g., through
a 180 x 180 spherical harmonic expansion of global topography.

Of special importance is the conventional gravity terrain correction. It is not a

"terrain effect" in the sense of the above definition, since different mass models are

343
 GRAVITY FIELD TERRAIN EFFECT.....

taken into account, depending on the elevation of the computation point. Rather, the
terrain correction should be looked upon as an auxiliary quantity, taking into account
the "unlinear" part of the total effect of the topography, i.e.,

AgBC = 2 ~ r G p h - t c (1)

where the terrain effect AgBC is the conventional refined Bouguer correction and
27rGph the simple Bouguer correction (p : density, G : gravitational constant,
h : point elevation). The terrain correction tc is always positive, and typically one
order of magnitude smaller than the simple Bouguer term.
The computation of terrain effects traditionally involves time-consuming
numerical integrations, where the effects are summed up using rectangular prisms of
uniform density as basic "building blocks." However, since digital elevation data are
nearly always made available in regular grids, application of methods able to utilize
FFT (Fast Fourier Transform) would certainly be attractive in terms of speeding
computations.
In the present paper such methods will be outlined, especially applicable for
computing land gravity terrain corrections and ocean isostatic geoid effects, two types
of reductions where typically lots of data need to be processed. The FFT methods are
based on various approximations, allowing the basically unlinear terrain effect integrals
to be replaced by a sequence of convolutions involving the elevations and their powers.
The accuracy of these approximations will later in the paper be evaluated by comparing
results obtained by the FFT techniques with similar results obtained using a prism-type
space domain integration program, using actual data from a variety of areas.
Recently 3ideris (1985) has published results for gravity terrain corrections
computed by FFT. The terrain correction results of the present paper, based on a more
elaborate two-stage approach and more detailed height data, support his general
conclusion that FFT-terrain corrections may obtain accuracies at a fraction of a regal.
Results for topographic/isostatic reductions are based on an integral expansion technique,
essentially analogous to the methods of Dorman and Lewis (1974) and Parker (1972).
Generally good results are obtained for gravity anomalies and geoid undulations, while
deflections of the vertical are less accurate when using only the first terms of the integral
expansions.

2. FFT and Terrain Effect Integrals

In order to be able to utilize the computational speed of FFT, the terrain effect
integrals must be approximated by two-dimensional convolutions, i.e., expressions of the
form

t(x,y) = ff(x-x, ,y-y')h(x' ,y')dH' = f(x,y)*h(x,y) (2)

rl
where 11 is a reference x - y plane, dH = dxdy the plane integration element, t the
result, h the data and f a kernel function determined by the type of gravity field
quantity and terrain reduction wanted. The convolution (2) becomes a simple product in
the frequency domain

344
 R. FORSBERG

T(u,v) = 7 ( u , v ) ~(u ,v) (3)

where the Fourier transform pair is given by

t(u,v) =;t(x,y) e -i(ux+vy) dxdy (4)

H

t (x , y) = -~-~21r2
f "~(u , v)ei(Ux +VY) dudv (5,

For actual computations, the Fourier transforms are approximated by the

discrete Fourier transform, FFT, resulting in t being expressed by a discrete periodic
convolution
n m
t (xi ,yj) = Z Z f (Xi_k, yj_s h (Xk ,ys Ax A y (6)
k=l s

where x i , i = ] . . . . . n and yj , j = ] . . . . . m are the coordinates of a given data grid

with spacing ( A x , A y ) . The functions in (6) are to be viewed as periodic, i.e.,
f(x i ,yj) = f(Xi+n,Yj+m) for all i and j . The above should be familiar to most
geodesists and geophysicists ; more details might be found in standard textbooks such as
Kanasewich (1975).
The simplest example of a convolution terrain effect integral is the so-called
condensation approximation. Consider Fig. I A . The potential at P for uniform
topographic density p will be

T(P) = G rh(XQ,yQ) dflQ , r = X/(Xp _yQ)2+(yp_yQ)2 +hp (71

II
if the masses are "condensed" as a mass layer at sea level. The function l / r has a simple
analytical Fourier transform when hp ~ 0 , (see Eorsberg,1984, sec. 6.2). However, to
avoid "edge effects" coming from the basic periodicity of EFT, it turns out to be often
more advantageous to transform both the elevation data (h) and the kernel function
( l / r ) numerically, truncating the kernel at a certain distance dma x . Setting the kernel
to zero at distances greeter than dma x is equivalent to computing the terrain effect out
to a distance of only dma x from the computation point. Since the FFT approach
implies a computation at all points of the original data grid, the use of truncation thus
means that all points in the grid with a distance greater than dma x from the edges of
the grid will be computed without errors due to the periodicity (Fig. 2). When using
results from this "reliable" central area only, there is no need for the arbitrary windowing
of data otherwise so common for EFT applications.
The influence of the distant zones, i.e., regions outside draax, may if
necessary be computed analogously by FFT, using less detailed (mean) elevation data,
now truncating the kernel both below dma x and a new, sufficiently large outer radius.
This "hierarchical" approach is especially important for gravity terrain corrections,

345
 GRAVITY FIELD TERRAIN E F F E C T .....

"7":=
CO

HEIGHT DATA KEItNEL RESULT

Fig. 2 - Schematic convolution giving terrain effects, using a periodical//extended

kernel truncated at distance d max "

which demand very detailed elevation data to compute accurate inner zone effects. (In
actual implementations, a rectangular rather than circular "filter domain" should be used
in order to avoid overlaps and gaps between detailed and coarse grid elements along the
circle of radius dma x .)
By approximating (7) with a discrete convolution (6), a discretization error is
introduced. Basically the continuous terrain effect integral is replaced by a sum of effects
of vertical "mass lines" (in (7) degenerated to mass points) at the nodes of the grid.
However, since the result usually needs to be interpolated from the computation grid to
actual "station" locations, the relative influence of the discretization error is diminished,
and improved prism-type terrain representations would probably not give significantly
better results.
In the sequel, results of interpolation to actual station locations have been
tested by both simple bilinear interpolation and windowed bicubic spline interpolation,
using a 5 x 5 or 9 x 9 subgrid surrounding the point to be interpolated. As expected,
the 9 x 9 spline algorithm yields the best results.

3. Gravity terrain corrections by FFT

The gravity terrain correction due to a terrain of uniform density p may be
written as a triple integral (cf. Fig. 1B)

tc(P) = Gofj~ah dzdIIQ, r = X / ( X p - X Q ) 2 + ( y p - y Q ) 2 + ( z - b p ) 2 (8)

1-1 hp
If the maximal terrain slopes immediately surrounding the computation point
are assumed to be small, then

! 1_._
r"T- %3 , r 0 = X / ( X p ' - x Q ) 2 +(yp _ y Q ) 2 (9)

and

tc (P) =~ ~-1 Gp f (hQ - h3p' ) 2 dl-IQ

ro
(10)
II

346
 R. FORSBERG

This is the so-called "linear approximation" to the terrain correction ; see Moritz(1968a).
The approximation is closely related to the approximations involved in the familiar
Molodensky boundary value problem of geodesy. The approximation errors are typically
a fraction of a mgal, but may become large; see Forsberg (1984, sec. ?.4). The formula
(10) may be written as a combination of convolutions

tc (P) = ~ Gp d]IQ + drlQ

. j r~ hp ro
lI II rl
where it is noted that the second integral in (11 ) is a constant. For practical computation
by FFT, the kernel function ] / r o3 is evaluated in a grid corresponding to the data grid
and truncated at the maximal computation distance dma x . For x = y = 0 a zero value
may be substituted since it is easily shown that in the discrete convolution expression
]
tc(P) = ~ Gp [h2 *f+h~ fo-2hl, ( h , f ) ] , f= (x2+ y2) -3/2 (12)

any actual value of f at zero drops out. In (12) the constant fo may be obtained from
the DC-component of f , while the convolutions h 2 * f and h * f may be obtained in
one operation, using as data the complex quantity h + i h 2 , which after F F T -
convolution with f yields the required convolutions in the real and imaginary parts of
the complex result respectively.
Terrain corrections for outer zones (distances greater than dma x) may be
computed analogously from a more coarse mean elevation grid, now letting the filter f
only be nonzero between the inner zone and outer zone radii. However, terrain
corrections for remote zones are known to have a systematic dependency of the actual
elevation of the computation point. Since the elevations of the computation points are
in principle given by the detailed grid, it is necessary to interpolate the convolution
results h 2 * f and h * f to the detailed grid knots before computing the final outer
zone corrections by (12), using hp values from the detailed grid. The reason for this
approach is that the coarse grid does not resolve the topography sufficiently. Using (12)
directly for the outer zones without the convolution interpolation has been found to
yield unnecessary errors of several regals in the examples of this section.
In the following examples, results of FFT-based terrain correction effect
computations have been compared to results from a conventional integration program
(Forsberg and Tscherning, 1981). The used program builds up a model of the terrain
masses using rectangular prisms, using a spline-densification of the given elevation data in
the immediate neighborhood of the computation point to reproduce sloping surfaces.
Results of this program are used asa standard against which the FFT results are evaluated.
However, since the "prism" method by its very nature obviously also is an
approximation, comparisons between the methods are not necessarily equivalent to the
"real" terrain correction computation accuracy.

347
 GRAVITY FIELD TERRAIN EFFECT.....

Example 9 Terrain corrections in Jameson Land, Central East Greenland

Terrain corrections in 42 actual gravity stations in a small 9 x 9 km area around

J.P. Koch Field, Jameson Land, central east Greenland, are used here to illustrate the
accuracy of the FFT method versus height spacing and interpolation method. The
topography of the test area consists of Mesozoic sediments, with a total elevation range
from ca. 300 m to 908 m at the summit of J.P. Koch Fjeld, as outlined in the map,
Fig. 3, which additionally shows the location of the gravity stations (surveyed in 1982 in
order to investigate local geologic structures associated with J.P. Koch Field, the
mountain in the center of the plot).

mA~ Lm

u]
ZOd
2~

7841

Fig. 3 - Map of J.P. Koch Fjeld area, Jameson Land, East Greenland. Contour
interval 50 m . Gravity stations shown with dots.

Available elevation data are a photogrammetric ]00 m x ]00 m height model,

covering the area with a ca. 3 km wide margin, and an approximately | kin x ! k m
height model of the remote zones, derived from existing maps. Terrain corrections,
computed with density 2.67 g/cm 3 to a distance of 50 k m , are shown in Fig. 4 and

348
 R. FORSBERG

64 L km L55 km

Fig. 4 - Terrain corrections computed by FFT, J.P. Koch Field. Contour interval
I mgal.

Table 1. The Fourier transformations involved 150 x 150 and 128 x 120 points in the
inner and outer zones respectively. The "radius" of the used square inner zone was
3 k m . Interpolation to station locations was done by a 9 x 9 point spline procedure.

Table 1

Terrain corrections (mgal) computed to 50 km by FFT

(J.P. Koch Field, Greenland)

mean standard
min max value deviation

Computed terrain corrections 1.26 9.68 3.66 2.00

"Difference "prism" - " F F T " - 0.29 0.92 0.14 0.21

349
 GRAVITY FIELD TERRAIN EFFECT.....

Compared to the prism-results, the F F T approach is seen to give an r.m.s, error of

around 0.25 toga|, with a slight bias toward too low values. This error level is
satisfactory for most purposes.
To give insight into the influence of the elevation and spacing, the inner zone
effects have been computed using various elevation grid spacings, the coarser grids being
produced by averaging the 100 x ]00 m model. The results, outlined in Table 2, clearly
illustrate the need for very detailed height data. In the present area even a 400 m grid
spacing is too coarse, especially for the summit station. For the detailed grid, the FFT
method gives very good inner zone results, with a maximal error of only 0.]3 mga]. As
additionally shown in the table, this error estimate increases slightly when using less
smooth grid interpolation procedures.

Table 2
Inner zone terrain corrections (to 3 kin) computed by FFT,
as a function of height grid spacing and interpolation type
(J.P. Koch Fjeld, Greenland)

grid std.
(unit : mgal) spacing min max mean dev.

Reference ("prism m e t h o d " ) 100 m 0.63 5.79 2.17 1.07

Difference to F F T m e t h o d 100 m - 0.12 0.13 0.00 0.06

200 m - 0.02 1.25 0.22 0.24
(ref. - F F T , 9 x 9 spline) 400 m - 0.03 2.95 0.60 0.54
800 m 0.41 3.76 1.14 0.69

Interpolation type : bilinear 100 m - 0.24 0.20 - 0.01 0.07

5 x 5 spline 100 m - 0.14 0.16 0.00 0.06
9 x 9 spline 100 m - 0.12 0.13 0.00 0.06

Example : Terrain corrections, white Sands test area, New Mexico

102 gravity stations from New Mexico were used for an additional test using
more coarse elevation data. The topography of the area, distribution of gravity stations
and computed terrain corrections are shown in Fig. 5. The used set of gravity stations is
the same as the stations used in Forsberg and Tscherning (1981), northern block. Terrain
corrections were computed to 50 km in a two-step procedure, using 0.5' x 0.5' heights
in an inner zone out to ]5 k m , and 2' x 2' mean heights for the outer zone (EFT grids
were 160 x ]60 and 80 x 80 points respectively). Results compared to the prism
method (Table 3) are surprisingly good, much better than the expected overall error
originating from lack of resolution of the elevation data.

350
 R. FORSBERG

-107 ~ -tO6 ~
34~ 1 34 ~

v eo.~
o0.4 g o ,2
5 OOI~O.L

O 0 . 2 {~ O 0 . 2 91
~600
O0.L eo.2
eo.l eo.
eO.~o. ~ eo.2 #
eo.i eo.3
"&o., eo.5
o 0 =oo L ..eo 9 o 4

og.l
O0.2
go )2.5

,o;~ 01 ,o

(DO .3 O0.3 L
0 o .4 9
eO .2
O| .2 9
DO.I e0.7
9 0., O0.2 110.4

DO. Le O .2 ~.2_ 9 .e0.6

33 ~ ~ J 33 ~
-107 ~ o~ ~ -106 ~
o

Fig. 5 - Topography (200 m contour interval) and gravity stations with computed
terrain corrections (regal). Northern part o f White Sands test area, New
Mexico.

Table 3

Terrain corrections (mgal) computed by FFT (New Mexico)

min max mean std. dev.

Computed corrections 0.06 32.48 0.96 3.34

Difference "prism" - " F F T " - 0.16 1.33 0.07 0.15

351
 G R A V I T Y F I E L D T E R R A I N E F F E C T .....

4. Isostatic gravity effects

When terrain corrections have been computed, they may be used together with
the simple Bouguer correction to provide a complete topographic reduction. Also, RTM-
reductions may usually be obtained sufficiently accurate by applying a simple Bouguer
reduction to the reference elevation level (cf. Forsberg, 1984, sec. 7.]). For isostatic
reductions, however, the gravity effect of the isostatic compensation (the "isostatic
correction") must be computed separately. This may be done by FFT through an integral
expansion technique.

z i
Y}

Fig. 6

Consider a point P at sea level" let the compensation depth be D and the
depth to the base of the isostatic compensation be d ( x , y ) , where

d(x,y) = P--P--h ( x , y ) (13)

Ap

(p : density of topography, Ap : crust-mantle density contrast); cf. Fig. 6. Then the

isostatic correction will be

D+d

gi(P) = G A p ffz
H D
~ - dzdll (14)

By expanding z/r a to first order around a suitable reference level d o (e.g.,

the mean depth to the compensating masses)

z ~do ro - 3 d o (z-d o) ro = ~/(XQ -Xp) 2 + (yQ_yr)2 + d 2o (15)

r ~ ' ~ r oa + ro
s

and integrating with respect to z , one obtains after a few derivations

gi(P)= GAp f d + r-3do

o
r os
(D-d o )] +! f r -3do
dH
2
GAp o
ro
dl-[ (16)
rl ii

352
 R. FORSBERG

which are in the form of two convolutions involving d and d 2 . These convolutions
may again be performed most efficiently by Fourier transforming complex data
(d + id 2 ) and a complex kernel as well. From the computed spectra, the spectra of each
originally real component may be isolated through the familiar conjugate symmetries of
FFT; for details, see Kanasewich (1975), sec. 3.3. Thus a total of only three FFTsub-
routine calls are needed for isostatic corrections.
Since the FFT method gives results on a constant level, a small error is made if
this method is used for land gravity stations of varying elevations. This error is, however,
expected to be very small, as the compensation depth is typically at least an order of
magnitude larger than the variation in the station elevations. The insignificance of the
error is indicated by the following example.

Example : Isostatic corrections, central east Greenland

Isostatic corrections were computed at 25 stations with a ]00 km spacing in

central east Greenland (Fig. 7), covering the area from the coastline to the edge of the
inland ice, where elevations are above 2 000 m . Gravity stations were located on the
terrain surface and the comparison between rigorously computed "prism" isostatic
corrections and the FFT-method are shown in Table 4. The obtained accuracies at the
regal-level are fully satisfactory, considering the large magnitude of the effects. The
maximum difference was only - 1.8 mgal.

2000 ~"~

.lr
sr

Fig. 7 - Perspective view of central east Greenland region, from the Norvegian-
Greenland sea to the ice cap margin, 5' x 10' mean heights.

The isostatic gravity effect formula (16) may be viewed as the space domain
equivalent of Parker's frequency domain expansion formula (Parker, 1972), expanding
only to the second order. Letting ~o = (u 2 + v 2 ) t / ~ denote the radial wave number,
the full frequency domain expansion becomes

,~i(U,V ) - 2~'Gp ~ -OO

- n e-~ao 'tin ( u , v ) (17)
n=l n !

353
 GRAVITY FIELD TERRAIN EFFECT.....

Table 4

Isostatic corrections (mgal) computed by FFT

(Central East Greenland)

mean std. dev.

Computed effects - 103.74 47.93
Difference "prism"- "FFT" - 0.88 0.66

As shown by Parker, this series converges at least as rapidly as ~(dmax/do)n , where

dma x is the maximal thickness I d - D I of the compensating mass body. The effect of
neglecting higher-order terms in (16) may therefore be assessed fairly easily. Thus, in the
present Greenland example, maximal root thickness would be roughly 9 km which
combined with a d o value around 34 km yields a relative error bound for the third-
order term at 1.8~ . This number compares reasonably with the statistics of Table 4,
suggesting even better results may be obtained by higher-order expansions, although the
nonconstancy of the actual station elevations ultimately limits the obtainable accuracy.

5. Terrain reduction of ocean and land geoid data

With the availability of a large number of oceanic geoid undulations from

SEASAT and GEOS-3 altimetry, efficient FFT-based isostatic reduction methods will
prove to be very valuable, both in order to enhance "geologic" geoid anomalies and for
the study of actual ocean isostasy. A recently released global bathymetric data base
"SYNBAPS", covering the oceans of the world t o latitude 75~ with 5' x 5' depth
values, is especially suited for such studies of satellite altimetry results. In the present
section SYNBAPS bathymetric data will be used for evaluating accuracies of the isostatic
integral expansion method. In addition also FFT methods for land geoicl terrain reduction
(RTM) will be briefly discussed in the sequel.
The isostatic geoid effect consists of the influence of the water body (density
anomaly p "" - ] . 6 4 ) and the associated isostatic compensation. Referring again to
Fig. 6, the geoid effect associated with the bathymetry may be expressed as

D+d

Nb (P) - Tb(P)3' = - - G O l f 3 ' -rl dzdII (18)

FI D
where D in this case is the computation height above sea level, and 3' normal gravity.
A first-order expansion of r -1 around a mean depth d o

1 d
l(z) ~ r ~ (Z-do), r o = X/(XQ-Xp)' + ( y Q - y p ) = + d2o (19)
r ro o

results in

354
 R. FORSBERG

17
ro
s
ro
d]-[ - - ,i,,,o
2
Gp

17
_-3- d l ]
ro
(20)

This expression is again a sum of two convolutions, which as in the gravity case may be
implemented in practice by just one "set" of Fourier transforms, packing the depths
( d , d 2 ) and the two convolution kernels as complex arrays.
The effect of the isostatic compensation may be expressed analogously, and the
results added to give the complete bathymetric/isostatic effect. However, for a "pure"
oceanic area, the thickness of the isostatic compensation mass body ( d ' ) will be
proportional to the depth d through

d' = p d = a d (21)
Ap

SO the total isostatic effect may be obtained as a single sum of convolution as

#

FI
ro re' 'Oro3 lro3
~ dll

(22)

'o I'21 '~ '1

H
ro ro
dl-I

where the primed quantities refer to the isostatic compensation layer. Hence, total
isostatic effect requires also only a single "set" of convolution Fourier transforms.
For practical evaluation of the integral kernels of (20) or (22) at sea level
(b=0), the value of the kernel at x = y = 0 needs modification, since it might be
singular (if d o = 0 ) . What in effect needs to be done is to replace the inherent mass
line representation of FFT with a spatial density distribution (e.g., a small cylinder)
immediately below the computation point. This is obtained when the kernel value at
• = y = 0 is expressed as the integral (mean value) over a small circle C, with area equal
to the grid element area A x A y . Use of the simple integrals

f l d [ I = 27r V~s2 + d 2o - d o (23)

C
and

[+ d ~ d l ] = 21r [ 1 ,o] (24)

where S _-__ 0.56 ~/ A x A y is the radius of C , gives these "zero-distance" geoid

contributions directly for the bathymetry (20). Neglect of these corrections might result

355
 GRAVITY FIELD TERRAIN EFFECT.....

in large errors (several decimeters), especially when the grid spacing is large comparedto
the expansion level d o .

Example : isostatic #eoid effects, Faeroe Islands area, North Atlantic

The accuracy of the FFT method for oceanic isostatic geoid effects has been
tested in a 6 ~ x ]2 ~ area of the North Atlantic around the Faeroe Islands. The area
bathymetry is quite varied, containing a continental shelf area, part of the Faeroe-lsland
ridge and some major banks, with a depth range from 0 to below 3 000 m (Fig. 8 ) .
The bathymetry data used were 5' x 10' SYNBAPS derived depths, and the comparison
to effects computed by the "prism"-integration was done in 49 points on a ] o x 2 ~
grid. The isostatic computation was carried out to 200 km distance, Fourier
transforming a larger area o f totally ]28 x ]28 depth points. The computed isostatic
geoid is shown in Fig. 9. The comparison results (Table 5) show more than sufficient
accuracy is obtained by the FFT method, compared to the accuracy of satellite
altimetry.

_-14 ~ ..o .,+ -2 ~

6C 50 ~

58 ~ 58 ~
--I"9 --2 ~
L, I 0 ~ -6 ~

Fig. 8 - Depth map (km) of the North Atlantic test area, contour interval 200 m
(an apparent error in the SYNBAPS bathymetric data base at 6 ~
plays no role in the present study}.

356
 R. FORSBERG

_14 ~ _ _^ _^ -2 ~

6 60 ~

58' 58 ~
_2 ~
-"~ ~ ~l~ 7" -6 ~

Fig. 9 - Isostatic geoid of the North Atlantic test area, computed by FFT, contour
interval 0.2 m . (Airy isostatic model, normal crustal thickness 32 k m ,
density contrast 0.4 g/cm 3 ).

Table 5

C o m p u t a t i o n o f ocean geoid effects (m) by FFT,

Faeroe Islands area

min max mean std. dev.

Computed isostatic effects - 4.61 - 0.02 - 1.66 1.21

Difference "prism" - "FFT" - 0.03 0.07 0.02 0.02

357
 GRAVITY FIELD TERRAIN EFFECT.....

Example : RTM geoMeffects, Colorado area

Geoid data on land determined through a combination of GPS and levelling

refers to the actual surface of the earth. (Such data are, in fact, height anomalies rather
than geoid undulations in the classical 9sense.) When geoid terrain effects are wanted at
the various elevations of actual topography, the previous expansion method can in
principle not be used, since it provides values on a plane of constant elevation. However,
geoid effects associated with R T M - t y p e reductions (Fig. ID) are small, typically below
the meter level, so for such effects it might make sense using the integral method, now
keeping only the first term, i.e., using the condensation approximation for the residual
topographic irregularities.
The accuracy of the condensation approximation, implemented using FFT, has
been tested in an extremely rugged topography, a 2 ~ x 2 ~ area of the Rocky Mountains
of Colorado, located west of Denver. The area, located between latitude 38 ~ - 40 ~ N ,
longitude 108 ~ - 106 ~ W, has maximal topographic elevations above 4000 m . As the
reference surface, elevations from the spherical harmonic topography expansion of Rapp
(1982) were used, and comparison with prism integration results was done in 25 points
on a 3 0 ' x 3 0 ' grid, using 2 . 5 ' x 2 . 5 ' mean heights. The results, shown in Table 6,
indicate that even the simple FFT-condensation approach may produce useful results for
many applications, although the accuracy might not be sufficient for analysis of the most
precise GPS geoid data.

Table 6

RTM geoid effects (m), computed by FFT condensation approximation

(Colorado)

min max mean std. dev.

~Computed geoid effect - 1.93 2.69 0.47 0.98

9Difference "prism" - " F F T " - 0.10 0.08 0.00 0.03

6. A note on deflections of the vertical

Astrogeodetic deflection stations are, like gravity stations, situated at the actual
surface of the earth, thus making application of the integral expansion methods for
deflections less straightforward (except such cases as predicting deflections aloft or at sea
for use by inertial navigation systems).
On land, topographic effects on deflections correspond to terrain corrections,
since a Bouguer plate has no effect on deflection values. Therefore, a method similar to
the gravity terrain correction "linear approximation" might be attempted. In the plane
approximation, deflection terrain corrections are given by

hQ }

r/t (P) 1 o ff{ yQ-y

_

rl :hp
XQ-Xp
ra
dz dIlQ , r = ~/(XQ - X p ) 2 + (yQ _yp)2 + ( z _ h p ) 2

(25)

358
 R, FORSBERG

In the "linear approximation" of sec. 3, r -3 ~ r o 3 was put outside the

z-integral. For deflections, this approximation becomes identical to the condensation
approximation, evident since the integral (25) will be independent on the actual value
of hp. The condensation approximation is usually n o t sufficiently accurate for
deflection computations. Experiments in a rugged area of central Norway, using 1 km
elevation data, have yielded r.m.s, errors in EFT-computed deflection terrain corrections
at the ] " level. To get improved results, it will thus be necessary to expand
the integral (25),
This might be done as outlined in the sequel. Expanding r -3 around the level
of P

3
r -3 ~ ro 3 --7ro-5 ( z - h p ) 2 , r o = X/ (XQ- Xp)2+ (yq'yp)2" (26)

the integral (25) then becomes

~t(P) l= _, YQ-Yp
r/t (P)) rl XQ -Xp
I[7o1 ( h Q - h p ) - 7l, r7o( h q - h p ) a dIIq (27)

By introducing the functions

f = I Y l r-~o ' g= l:il l

.T I 'o = l'%//'x"2"-+--~y21 (28)
x ro
the deflection terrain corrections (27) may finally be written as a linear combination
of four convolutions (of which the first represents the simple condensation
approximation)

~t(P) I = _GP_-
f , h - 3h~ (g,h)+ 3hp (g*h 2 ) - g*h 3 (29)
~t (P) 7

This expression involves transforming two filters (f and g) as well as the elevation
data h and their second and third powers. Since deflections are both sensitive to local
and more distant topography, a hierarchical two-step procedure, analogous to the EFT-
computation of gravity terrain corrections, will usually be needed. The EFT-
computation of deflections of the vertical therefore becomes a rather complicated
process, involving many transforms. Considering that usually only quite a few deflection
data are available for gravity field modelling (e.g., geoid prediction by collocation), it
might often be easiest simply to compute needed deflection terrain effects by a
conventional "prism" integration program.

7. Conclusions

In the previous sections it has been outlined how the traditionally time-
consuming processes of computing terrain corrections, isostatic effects, etc. may be
formulated in terms of convolution integrals, which may be computed very fast by the
Fast Fourier Transform. Especially for such often massive problems as computing terrain

359
 GRAVITY FIELD TERRAIN EFFECT.....

corrections for continental gravity data, or isostatic analysis of satellite altimetry geoid
data, the FFT methods might provide a reduction in computing time by several orders of
magnitude, compared to space domain integration methods.
In the numerous numerical examples of the paper, the FFT methods were
demonstrated to yield satisfactory accuracies in general. Typical gravity terrain
correction errors associated with the FFT approximations were at a fraction of a mgal,
isostatic geoid errors at a few cm level. The used expansions of terrain effect integrals
were all fairly simple, involving only the elevation/depth data and their squares. Only
for deflections of the vertical do higher powers of the elevations seem to be needed.
Other gravity field quantities (second-order gradients) and corresponding
integrals for magnetic terrain effects may be treated along similar lines. It would also be
possible to use most of the presented methods in a straightforward fashion, given the
realistic case of regional variations in the average topographic density. .
The effectiveness of the FFT terrain effect computations have already proven
their use for major gravity field projects undertaken at the Danish Geodetic Institute.
They are presently used as auxiliary tools for analysis of SFASAT/GEOS3 data in the
North Atlantic, as well as for geoid prediction in Scandinavia and Greenland (Tscherning,
1985).

O O

REFERENCES

L.M. DORMAN and T.R. LEWIS : The use of nonlinear functional expansions in calculation of the
terrain effect in airborne and marine gravity and gradiometry. Geophysics,39, No. 1, pp. 33-
38, 1974.
R. FORSBERG and C.C. TSCHERNING : The Use of Height Data in Gravity Field Approximation
by Collocation. Journ. Geoph. Res.,86, No. B9, pp. 7843-7854, 1981.
R. FORSBERG : A study of terrain reductions, density anomaliesand geophysical inversion methods
in gravity field modelling. Rep. 355, Dept. Geod. Sci. Surv., Ohio State Univ., Columbus, 1984.
E.R. KANASEWICH : Time Sequence Analysis in Geophysics. The University of Alberta Press,
second ed., Edmonton 1975.
H. MORITZ : On the use of the terrain correction in solving Molodensky's problem. Rep. 108, Dept.
Geod. Sci., Ohio State Univ., Columbus, 1968.
R.L. PARKER : The rapid calculation of potential anomalies. Geophys.J. R. Astr. Soc., 3], pp. 447-
455, 1972.
R.H. RAPP : Degree variances of the earth's potential, topography and its isostatic compensation.
Bull. Geod., 56, pp. 84-94, 1982
M.G. SIDERIS : A Fast Fourier Transform method for computing terrain corrections. Manuscripta
Geodaetica, ]0, pp. 66-73, 1985.
C.C. TSCHERNING : Geoid Modeling Using Collocation in Scandinavia and Greenland. Marine
Geodesy, 9, No. 1, pp. 1-16, 1985.

Received : 21.12.1984
Accepted : 25.06.1985

360