- Email: [email protected]
Terrain correction program for regional gravity surveys
Pergamon
Computers & Geosciences Vol. 20, No. 6, pp. 961-972, 1994 Copyright © 1994 Elsevier ScienceLtd 0098-3004(94)E0015-L Printed in Great Britain. All rights reserved 0098-3004/94 $7.00 + 0.00
TERRAIN CORRECTION PROGRAM FOR REGIONAL GRAVITY SURVEYS X. Q. MA* and D. R. WATTS Department of Geology and Applied Geology, University of Glasgow, Glasgow GI2 8QQ, U.K. (Received and accepted 3 November 1993)
Abstract--An algorithm is presented for the computation of gravity terrain corrections from a digitized grid representing topography. Different approximations for topography at various distances from the gravity station are used to make the computation more efficient. For the distant zone (30 < r < 50 km) topography is approximated using a vertical line mass distribution. In the intermediate zone (2 < r < 30 km) a segment of a hollow vertical cylinder is used. The near zone (0.5 < r < 2 km) is approximated using vertical rectangular prisms with sloping upper faces. A new method is used to compute the correction for the inner zone which is the digitized square where the gravity station is located. The gravity effect of topography in the inner zone is calculated from four traingular prisms with vertices made of the station location and the corners of the square. The computer program has been tested comprehensively with a regional gravity database from the Southern Uplands of Scotland in which all terrain corrections were carried out originally using the Hammer chart method. Key Words: Terrain correction, Gravity survey, FORTRAN 77.
INTRODUCTION Determining the terrain correction for a gravity anomaly is a tedious but potentially important task, especially in areas of rugged terrain. The calculation must be done efficiently as well as accurately and this may be expedited by using different approximations of topography depending on the distance from the gravity station. In this paper, we present a F O R T R A N - 7 7 computer program which utilizes various features of previous methods given by Bott (1959), Karlemo (1963), Nagy (1966), Ketelaar (1976), Blais and Ferland (1984), and Lopez (1990). A new method for the inner zone correction has been developed. Topography is represented by a digitized grid of equal squares of convenient size as shown in Figure 1. In the example considered here, topography is divided into a grid of 1 km squares and the digitized point is the center of each of the squares. If a square contains a gravity station, the four corners of that square are added to the database as well as the position and elevation of the gravity station. Note that the gravity station does not have to be centered within a square of the digitized grid. The zones are defined by the distance (r) of the digitized points from the gravity station being considered. The far distant zone (r > 50 km) is neglected. The distant zone is defined as (30 < r < 50 kin). The intermediate zone is between (2 < r < 30 km). The near zone is between
(0.5 < r < 2 km) and the inner zone is the digitized square with the gravity station inside. The inner, near, and a part of the intermediate zones are shown in Figure 1. DISTANT ZONE CONTRIBUTION To reduce the calculation time without losing accuracy, the vertical line mass Equation (1) is used to determine the contribution of the distant zone (30 < r < 50 km) topography (Bott, 1959): g = GpA
x / ( r 2 + z2) 3
(1) If h is less than about one-quarter of r a good approximation to (I) is given by the following simplified formula: g =
GpA hz 2 x r3.
(2)
where G is the gravitational constant, r is distance from station to center of square, h is difference in elevation between the square center and station, A is the area of square, and p is density of the topography. INTERMEDIATE ZONE CONTRIBUTION In this zone (2 < r < 30km) it is acceptable to approximate a prism as a segment of a hollow vertical cylinder (Bott, 1959). The vertical component of gravitational attraction from a segment of a hollow
*Present address: Department of Earth Sciences, University of Liverpool, P.O. Box 147, Liverpool L69 3BX U.K. 961
962
X.Q. MA and D. R. WATTS
cylinder with angle ~, inner radius r~, and outer radius r 2 is given by
g=GpA(r2-r~+x/r~+h2-
r~2+h2).
expression for becomes:
~ / 7 + h2 = r 1 + 2 r2].
(4)
If we now replace r 2 and r~ by r + p and r - p, ~ by A/2rp (A is area of the square, r is distance from station to center of square, and p is half length of square side), the terrain correction for a prism can be approximated as g
GpA h2 2 r(r2-p 2)
(5)
where G, A, h, are as given in (1). A considerable amount of computation time is saved if squares beyond a certain distance from a station are grouped and treated as a single square (Bott, 1959). In the area where 2 0 < r < 30km, sixteen squares are treated as one square (4 × 4 km). Where 15 < r < 20 km, four squares are treated as a single square (2 x 2 km), and where 2 < r < 15 km, a single square (1 x 1 km) is used to calculate the terrain correction using formula (5).
NEAR ZONE CONTRIBUTION
The near zone is shown as the stippled region in Figure 1. This area is represented by the rectangular prisms with the horizontal lower faces and sloping upper faces determined by the heights at the centers of the prisms and heights of neighboring locations (Ketelaar, 1976). The gravitational effect of a rectangular prism with the horizontal upper and lower faces and vertical sides is given by Nagy (1966). A modification by Ketelaar (1976) includes a top face with a constant slope towards the origin. The sides of a prism have a unit length D, the length of the unit digitized square (in this paper D = 1 km). The gravitational attraction of the prism, which has its vertical axis at x = iD and y =jK, is given as in (Ketelaar, 1976):
g=Gpfff~dv=Gp(1-cos~)DK(i,j),
terrain correction therefore
(3)
Expression (3) may be approximated further with the square root expanded in a power series and the terms with higher powers neglected,
( l_h2
the
INNER ZONE CONTRIBUTION
This zone is the square containing the gravity station. We develop a new method, which not only achieves the correction automatically, but provides acceptable accuracy. The basic idea is to construct four triangular prisms within the square containing the station with four elevation values at four corners of the square, and that of the station. Figure 2 shows the four triangular prisms with horizontal lower faces and sloping upper faces so constructed. The gravitational attraction of individual prisms can be obtained by integrating gravity over the volume of each triangular prism, and the effects of the four prisms are summed together, with the result added to the total terrain correction for that station. This method is a generalization of that suggested by Lopez (1990). Suppose that four elevations at the four corners of the square are all higher than the station elevation. Based on this assumption, the terrain effects of the four triangular prisms are considered individually as follows. Let us take the first triangular prism OP~P2Z2Z~ as an example to show how the terrain effect is obtained. The basic calculation is a simple integration over volume. The integration limits are shown in Figures 2A and 2B. The equations of lines OPt, OP2, and plane O Z l Z 2 are y = k~x, y =k2x, and z =a~x +b~y respectively. Therefore, the integration can be written as:
gl = Gp
dx
~x
dy
x/(x: +y2 + z2)3 (8)
(6)
where ~ is an angle showing the slope of the upper face [ ~ = a r c t g (h/r)]. K(i, j) is an eightterm algorithm which is a function of indices i and j only, that is K(i, j) depends uniquely on the position of the prism, and may be termed the position matrix (Ketelaar, 1976). A simplification has been made to the position matrix K(i, j). For prisms in this near zone, it is determined that K(i, j) is almost equal to (i 2 +j2)-05 (Ketelaar, 1987). The
Figure 1. Digitized topography and several zones with respect to gravity station at O. Inner zone is square containing station. Near zone is stippled and portion of intermediate zone is shown as white squares outside of stippled area.
963
Terrain correction program for regional gravity surveys ,
zs
'f
Y
P3 (x3, Y 3 ) ~
Y
_
P2
(x2, Y2)
_
X
74
P4
I
'%
PI
(x4, Y4
I I
(Xl, Yl ) P4
A
B
Figure 2. A--Perspective view of inner zone with four corners above stauon level. B--Plan view of inner zone. The result from the integration (8) and those for other three triangular prisms are given by: gl = G p x l f ( m ,
n, u, v)],, = k.... k~......... b,
(9)
g2 = G p Y 2 f ( m , n, u, U)]m=k'3,n=k~,u=b 2......
(10)
gs = G p x 3 f ( m , n, u, V)[m=k.... k3. . . . . . . . . b3
(1 1)
g4 = G p Y 4 f ( m , n, u, v)[,, = ka,,=ki,, =b,,~,=,4
(12)
where f ( m , n, u, v) = In
u/3
n + x/1 + n 2 m +x/l+m
~/(
1
2
X/1 +/32
u/3 )2
l + u2
n+l~v2+
n+l--~vZ
+l+v
m + l~-~v2 +
m + i--~v2
+l+vZ
u2/32
2
(l+v2) z
.In//2/32
(1 + V2)2 (13)
and zty2--z2yl a I =
a2 =
xIZ2--X2z 1 b I
x~y2--x2yl z2y 3 -- z3y 2
,
- -
xjy2--x2yl X223 -- X3Z 2
b2 -
x2y 3 -- x3y 2
x2y3
z3Y4-- z4y3 a 3 -=
x 3 z 4 - - X4 z3 ,
b3 -
,
Xay 4 -- x4y 3 z4Yl
--
x3Y4 -- x4y 3
z~y4
a 4 -=
,
b4
x4z~ -- x~ z4 --
x4y I -- xlY 4
kl = Y_L ,
k2 =
Y_! ,
x4Yl
k3 =
Y.-! ,
X1
XI
X3
I
1
1
ki=~,
k~=~,
,
-- x3y 2
k;=~,
-- Xly 4
k4 =
Y4 , X3
1
The terrain correction of the inner zone is obtained by summing four incremental contributions for the four triangular prisms. Thus we get for the total correction for the inner zone, g = gt +g2 +g3 + g 4 .
(14)
As stated, the value g is determined on the assumption that the four elevation values: Z1 P~, Z2P2, Z3P3, and Z4P 4, are all positive. However, in reality, these four positions in the square can be either higher or lower than the gravity station. One example is shown in Figure 3 in which Z1P~ is negative and the others remain positive. Within a triangular area OP~ P2, the triangular prism OQI P2 Z2 represents a hill whereas the triangular prism OQ1PLZ~ represents a valley. Although the vertical component of gravitational attraction at the gravity station from these two prisms are opposite in sign, in determining Bougour anomaly, the terrain corrections for both a hill and a valley must be added to the observed gravity. Therefore we need to only consider the absolute value of gravity. Because of symmetry, the volume of the prism OQ~ Pt Z~ is equal to that of the prism OQ~ Pl Z~, that is, integrating over the volume OQ~ P~ Z~ is equivalent to integrating over the volume OQ~ P~ Z l . Therefore, the terrain effect (w~) of two triangular prisms OQI PI Z'~ and OQ~ P2 22 can be obtained by subtracting the effect produced by the central prism OQ~ZIZ2 from the overall effect produced by the triangular prism OP~ PzZ2ZI. The same treatment should be applied to the triangular area OPI P4. To consider sixteen cases which may occur in a real situation (Fig. 4), the gravitational attractions (w2, w3, w4) from other three central prisms (OQ2Z2Z3, OQ3Z3Z4, and OQ4Z4ZL) also must be determined. The results from integration are given by: wl = gl - G p x l [ f ( m ,
k'4=k~. +f(m,
n, u, v)lm=, ..... ;. . . . . . . . . . h,
n, u, v)l,,_ ;.... k. . . . . '2.,'= a2]
(15)
964
X . Q . MA and D. R. WARTS
Z2 P2
(x2, Y2)
~'Y
Q2
(x3, Y3)
P3 ~ y - - k 3 x
~
/
Z4 P4
y=k2x
'"1~ X
6
Q]
y=14x
P1
P4 (X4' Y4)
Q4
(Xl, Yl)
A
B
Figure 3. A--Perspective view of inner zone with one corner at Pt below station level. B--Plan view of A.
where
w2 = gz -- G p y z [ f ( m , n, u, v)lm=rz,,,=k~,.=a3 ...... + f ( m , n, u, V)lm=k~,,,=t~,,,=a4. . . . . ]
(16)
Cl =
zlY2-J-z2y
1
--XIZ2--X2Z
,
d~
I
-
xly2--x2yl
x[y2-x2yl
w3 = g3 - G p x 3 [ f ( m , n, u, v )l,,,=t3,,,=l,3 . . . . 5.~=as x~z2-Fx2z~
--zzy2--z2y 1
+ f ( m , n, u, v)],,, = k. . . . t3......... a6]
(17)
C2 - -
W4 = g4 -- G p Y 4 [ f ( m, n, u. v)l,.=k'4,,,=r4,.=a7 ...... c3 + f ( m , n, u, v)lm=r.,. = k't,,,= as,,,= c,] Z3
Z4
g-w3-w 4
g - w 1 -W 2-W3 "W 4
(18)
~
x~yE-xzyl z2y3--~z3y2 x2Y3--x3Y2
,
d2
d3 =
x~y2-x2Y~
,
--X2z3--X3z2 x2y3--x3Y2
Z2
g'w2"w 3
g - w ] "W2
g-W 1 -W 3
g-W 1 -W2-W 3-W 4
g'W2*W 3
g-w2-w 4
~-W2-W 4
g
g - w 1 -W 3
g-W 1 -W4
g - W 1 -W 2
g-W 1 - w 4
g'w3"w4
g
Figure 4. Sixteen examples of possible terrain configuration in inner zone and corresponding equations for terrain corrections.
965
Terrain correction program for regional gravity surveys Table 1. Example of gravity station data file. N o is gravity station number; (x, y) and h are coordinates and elevations of gravity stations respectivley; z~, z~, z 3, and z4 are elevations at four corners of square containing station No
x (m)
y (m)
h (m)
z I (m)
z2 (m)
z3 (m)
z4 (m)
1
273350
625460
516.9
426.7
434.3
414.5
472.4
2
274640
627040
403.9
335.3
329.2
278.9
350.5
3
276230
628160
249.3
243.8
379.5
365.8
292.6
4
277500
628020
263.5
213.4
259.1
378.0
243.8
5
278650
628440
235.3
228.6
320.0
259.1
228.6
6
279920
627770
217.3
236.3
213.4
243.8
269.7
7
279540
629700
344.4
297.2
317.0
365.8
320.0
8
276990
625680
398.4
341.4
362.7
405.4
408.4
9
276180
626520
450.2
362.7
396.2
378.0
406.9
10
275370
625240
467.6
406.9
406.9
396.2
442.0
--z2y 3 -- z3y 2 C4 ~
d4 =
, x2y 3 -- x3y 2
x2Y3 -- x3Y 2
z3y 4 "b z4y 3
c5=
,
ds=
x3y 4 -- x4y 3
--z3y 4 -- z4y 3
x3Y4
d7
,
x4y I -- xly 4
x4Yl -- Xly 4
--z4y I -- 21y4
,
ds :
x4y I -- Xly 4
Cl
x4zl .-]- X12 4 X4Yl
C3
4' 1
t2' =
A F O R T R A N - 7 7 p r o g r a m (see the Appendix) was developed on a Sun S P A R C 4 W o r k s t a t i o n for the evaluation o f terrain corrections using the expressions as described. The p r o g r a m requires two i n p u t files; a gravity station file ( s t a t n . d a t ) and a terrain file ( b l o c k . d a t ) . The station file is illustrated in Table 1, consisting o f station numbers, grid coordinates (x, y ) (corresponding with the U K N a t i o n a l Grid), elevation o f the station, elevations at the four corners of a square (z~, z2, z3, and z4). The terrain file is illustrated in Table 2 a n d consists of block numbers, coordinates of block centers (x, y), average density of block, a n d sixteen elevation values of the
-- x4y 3
- - X 4 Z l - - X 1Z4 ,
= E'
FORTRAN-77 PROGRAM
d6 =
z 4 y I -.]- Z l y 4
c8=
,
-- x4y 3
X 3 Z 4 -~- X 4 Z 3
, x3y 4 -- x4y 3
C7 =
--X3Z 4 -- X4Z 3 x3Y4
C6 =
E q u a t i o n s for the sixteen cases are shown in Figure 4 a n d are used to calculate terrain corrections for the inner zone.
X 2 Z 3 At- X 3 Z 2
1
-- xly 4
C5
t3=
C7
ds' 1
t; =- t3'
1 t; = ?, .
Table 2. Example of terrain data file. centers; p is average density of block;
N o is ht-h~6
block number; (x, y) are coordinates of block are elevations of sixteen squares within block
NO
x (m)
y (m)
I00
226000
560000
p(g/crn3) 2.73
88.4
91.4
........
85.3
101
230000
560000
2.73
83.8
76.2
........
70.1
102
234000
560000
2.73
61.0
45.7
........
53.3
103
238000
560000
2.73
41.1
54.9
........
45.7
104
222000
564000
2.73
152.4
149.4 . . . . . . . .
152.4
105
226000
564000
2.73
121.9
96.0
........
106.7
106
230000
564000
2.73
99.1
121.9 . . . . . . . .
106.7
107
234000
564000
2.73
61.0
53.3
........
67.1
108
238000
564000
2.73
65.5
70.1
........
76.2
109
222000
568000
2.73
144.8
167.6
h I (m)
h 2 (m)
........
hi6
(In)
152.4
X. Q. MA and D. R. WATTS
966
Table 3. Comparison of terrain corrections carried out using our program with those done by Hammer zone chart method. No is gravity station number; (x, y) are coordinates of gravity stations; gtneWis terrain correction carried out by our program; gtom is terrain correction done by Hammer's method x (m)
y (m)
g~m (mg~)
gt~, (mg~)
e=or (%)
107
295280
634380
13.25
13.44
1.4
168
306530
631480
10.07
10.00
0.7
487
264240
601040
5.61
5.72
1.9
1061
314160
608940
5.18
5.17
0.2
1122
314690
611100
7.41
7.36
0.7
2159
233520
591080
5.24
5.58
6.0
No
centers of the individual squares. The order of the digitized squares within a block is shown in Figure 5. The program reads the station and terrain data and calculates the distance from the station for every square. The squares and blocks are classified according to their distance from the station so the appropriate algorithm can be applied. For the inner zone, many calculations and checks of elevation are required. Because the Equations (9)-(12) and (15)-(18) have similar patterns, a subroutine termed innerzone is designed so that the main program can call it many times, each with different coefficient values for the evaluation of the inner zone correction.
TESTING THE PROGRAM
The program was tested using the British Geological Survey gravity data set for the Southern Uplands of Scotland. This data set was selected becasue all of the stations have a terrain correction that was computed originally using Hammer's method (Hammer, 1939). For this purpose, 4526 gravity stations in the western Southern Uplands, which occupies 11,000 km 2, were used. The area was digitized with a total of 1 !,328
elevation points read from the Ordnance Survey (OS) map for 708 sixteen square blocks. In addition, 18,104 elevation points also were read from the OS map in order to calculate the terrain correction for the inner zone. The densities used for the terrain file, block.dat, were taken from Mansfield and Kennett (1963), Bott and Masson-Smith (1960), and Parslow and Randall (1973). The program was tested by comparing the new terrain corrections with those given in the British Geological Survey data computed using the Hammer chart method. It was determined that most of the new terrain corrections are within a few percent of the original values. We arbitrarily select six stations. The new and old terrain corrections are listed in Table 3. However, there were a few stations where new terrain corrections were substantially different from the original ones by up to a few mgals. To examine the problem, we computed the terrain correction for these stations using the Hammer chart method. In each situation, we determined that our Hammer method calculations produced values that agreed closely with those given by our program, indicating that some terrain corrections in the original data set are in error.
SUMMARY AND CONCLUSIONS 11
12
15
16
10
13
14
4
7
8
2
5
6
Figure 5. Terrain block illustrating order of digitization for file block.dat.
A program for gravity terrain correction has been developed for regional gravity surveys. The program is made more efficient by dividing the digitized terrain into different zones. Within each zone, different equations with approximations appropriate to the distance of the topography from the station are applied. A new method is used for the computation of the inner-zone correction, which involves the use of triangular prisms. The station is not restricted to the center of the inner square. The program was tested using gravity data from the Southern Uplands of Scotland. The results have shown that the computerized terrain corrections for most stations agree closely with the original corrections made by the Hammer chart method. Discrepancies were shown to be errors in the application of the Hammer method in the original data.
967
Terrain correction program for regional gravity surveys REFERENCES
Blais, J. A. R., and Ferland, R., 1984, Optimization in gravimetric terrain corrections: Can. Jour. Earth Sci., v. 21, no. 5, p. 505 515. Bott, M. H. P., 1959, The use of electronic digital computers for the evaluation of gravimetric terrain corrections: Geophysical Prospecting, v. 8, no. 1, p. 45 54. Bott, M. H. P., and Masson-Smith, D., 1960, A gravity survey of the Criffell granodiorite and the New Red Sandstone deposits near Dumfries: Proc. Yorkshire Geol. Society, v. 32, Pt. 3, no. 13, p. 317-332. Hammer, S., 1939, Terrain corrections for gravimeter stations: Geophysics, v. 4, no. 3, p. 184 194. Karlemo, B., 1963, Calculation of terrain corrections in gravity studies using the electric computer: Geoexploration, v. 1, no. 1, p. 56 66.
Ketelaar, A. C. R., 1976, A system for computer-calculation of the terrain correction: Geoexploration, v. 14, no. 1, p. 57~5. Ketelaar, A. C. R., 1987, Terrain correction for gravity measurements, using a digital terrain model (DTM): Geoexploration, v. 24, no. 2, p. 109-124. Lopez, H. R. B., 1990, F O R T R A N program for automatic terrain correction of gravity measurements: Computers & Geosciences, v. 16, no. 2, p. 237-244. Mansfield, J., and Kennett. P., 1963, A gravity survey of the Stranraer sedimentary basin: Proc. Yorkshire Geol. Society v. 34, Pt. 2, no. 9, p. 139 151. Nagy, D., 1966, The gravitational attraction of a right rectangular prism: Geophysics, v. 31, no. 2, p. 62-371. Parslow, G. R., and Randall, B. A. O., 1973, A gravity survey of the Cairnsmore of Fleet granite and its environs: Scott. Jour. Geology, v. 9, no. 3, p. 219-231.
APPENDIX Program TERRAIN ++++++++++++++++++++++++ ++++++++++++++++++++++++++++++++++÷++++++++ A F o r t r a n - 7 7 c o m p u t e r p r o g r a m for m a k i n g t e r r a i n c o r r e c t i o n s + for r e g i o n a l g r a v i t y surveys. D e s i g n e d and w r i t t e n by + X i n - Q u a n Ma & Doyle R. Watts, + D e p a r t m e n t of G e o l o g y & A p p l i e d Geology, + U n i v e r s i t y of G l a s g o w Scotland. + D e c e m b e r 1993 + +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
i n t e g e r mm, nn parameter(mm=708,nn=4526) real e x ( n n ) , e y ( n n ) , h ( n n ) , z l ( n n ) , z i ( n n ) , z 3 ( n n ) , z 4 ( n n ) real xx (mm) ,yy (mm) ,d (ram) ,hl (mm) ,h2 (mm) ,h3 (mm) ,h4 (mm) ,h5 (irma) real h 6 ( m m ) , h 7 ( m m ) , h S ( m m ) , h 9 ( m m ) , h l 0 ( m ~ ) , h l l ( m m ) , h l i ( m m ) real h l 3 ( m m ) , h l 4 ( m m ) , h l 5 ( m m ) , h l 6 ( m m ) , x ( 1 6 ) , y ( 1 6 ) , r ( 1 6 ) , b ( 1 6 ) real a l , a i , a 3 , a 4 , b l , b i , b 3 , b 4 , c l , c i , c 3 , c 4 , c 5 , c 6 , r t , c 7 , c 8 , d l , d 2 real k l , k 2 , k 3 , k 4 , k l l , k 2 2 , k 3 3 , k 4 4 , 1 1 , 1 2 , 1 3 , 1 4 , 1 1 1 , 1 2 2 , 1 3 3 , 1 4 4 real h t 4 , z l l ( n n ) , z 2 2 ( n n ) , z 3 3 ( n n ) , z 4 4 ( n n ) , g , g l , g 2 , g 3 , g 4 , g c , p real x l , x 2 , x 3 , x 4 , y l , y 2 , y 3 , y 4 , w l l , w 1 2 , w 2 I , w 2 2 , w 3 1 , w 3 2 , w 4 1 , w 4 2 real d 3 , d 4 , d 5 , d 6 , d 7 , d 8 , h t l , h t i , h t 3 , g g ( 4 ) , r 0 , g t ( n n ) , w l , w i , w 3 , w 4 open(unit=l,file='statn.dat') open(unit=2,file='block.dat') open(unit=3,file='outpt.dat') +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
ram n u m b e r of total blocks nn n u m b e r of total g r a v i t y stations ex x c o o r d i n a t e of g r a v i t y station (in metres) ey y c o o r d i n a t e of g r a v i t y station (in metres) h e l e v a t i o n of g r a v i t y station (in metres) z l l , z 2 2 , z 3 3 , z 4 4 e l e v a t i o n s of four corners of the square w h i c h c o n t a i n s the g r a v i t y station (in metres) xx x c o o r d i n a t e of a b l o c k centre (in metres) yy y c o o r d i n a t e of a b l o c k centre (in metres) hl-hl6 e l e v a t i o n s of 16 prisms w i t h i n the b l o c k (in metres) d a v e r a g e d e n s i t y of the b l o c k in g/cm**3 gt total t e r r a i n c o r r e c t i o n for one s t a t i o n in mgal gc is the g r a v i t a t i o n a l constant (gc=6.67x10**(-ll) N m**2/kg**2) +++++++++++++++++++++++++++++++++++÷+++++++++++++++++++++++++++++++
i00 150
print*,' print*,' ' + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ' print*,' + +' print*,' + Automatic Terrain Correction Program +' print*,' + +' print*,' ' + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ' print*,' ' gc = 6 . 6 7 3 2 E - I I do 150 k = i, nn read(l,100) e x ( k ) , e y ( k ) , h ( k ) , z l l ( k ) , z 2 2 ( k ) , z 3 3 ( k ) , z 4 4 ( k ) format(14x, if8,5f7.l,f5.2) continue do 300 i = i, m m read(2,200) x x ( i ) , y y ( i ) , d ( i ) , h l ( i ) , h i ( i ) , h 3 ( i ) , + h4(i),h5(i),h6(i) read(2,250) h 7 ( i ) , h 8 ( i ) , h 9 ( i ) , h l 0 ( i ) , h l l ( i ) , h l i ( i ) , + hl3(i),hl4(i),hl5(i),hl6(i)
+ + + + + + + + + + + + + +
968
X . Q . MA and D. R. WATTS 200 250 300
f o r m a t ( 5 x , 2 x , f7,2x, f7,2x, f4.2,6(2x, f6.1)) f o r m a t ( 1 0 ( 2 x , f6.1)) continue
C
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c c c c c c c c
k is i n d e x for the g r a v i t y s t a t i o n w h i l e i is for b l o c k number. T h e p r o g r a m p r o c e s s e s d a t a b l o c k by block; a b l o c k is c l a s s i f i e d b y its d i s t a n c e to the g r a v i t y s t a t i o n a c c o r d i n g to w h i c h d i f f e r e n t f o r m u l a s are applied. In this p r o g r a m the t e r r a i n is d i v i d e d into 7 zones a c c o r d i n g to d i s t a n c e r0: I. r 0 < 0 . 5 k m ; 2. 0.550km. N o t e the g r a v i t y e f f e c t of a r e g i o n r0>50 k m is n e g l e c t e d .
C
+ + + + + + + +
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
d o 5000 k = i, nn gt(k) = 0.0 d o 4000 i = i, m m b ( 1 ) = h l (i) b(2)=h2(i) b ( 3 ) = h 3 (i) b(4)=h4(i) b(5)=h5(i) b(6)=h6(i) b(7)=h7(i) b(8)=h8(i) b(9)=h9(i) b(10)=hl0(i) b(ll)=hll(i) b(12)=hl2(i) b(13)=hl3(i) b(14)=hl4(i) b(15)=hl5(i) b(16)=hl6(i) C
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c c
c a l c u l a t e the d i s t a n c e b e t w e e n g r a v i t y s t a t i o n a n d the b l o c k c e n t r e (xx(i),yy(i)):
C
(ex(k),ey(k))
+ +
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
r0=sqrt((xx(i)-ex(k))**2+(yy(i)-ey(k))**2) i f ( r 0 . 1 e . 1 5 0 0 0 ) g o t o 600 i f ( r 0 . 1 e . 2 0 0 0 0 ) g o t o 500 i f ( r 0 . 1 t . 3 0 0 0 0 ) g o t o 450 i f ( r 0 . 1 t . 5 0 0 0 0 ) g o t o 400 g o t o 4000 C
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c c c c c c c c
T o c a l c u l a t e the t e r r a i n c o r r e c t i o n for e a c h 4x4 k m s q u a r e T h e t e r r a i n in this r e g i o n (30
c
+ + + + + + + +
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
400
ha=(hl(i)+h2(i)+h3(i)+h4(i)+h5(i)+h6(i)+h7(i)+h8(i)+h9(i)+ + hl0(i)+hll(i)+hl2(i)+hl3(i)+hl4(i)+hl5(i)+hl6(i))/16.0 tt=abs(O.5*gc*l.
OE5*(d(i)*lOOO.O)*(l.6EO7)*(ha-h(k))**2/rO**3)
gt(k)=gt(k)+tt g o t o 4000 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ T o c a l c u l a t e the t e r r a i n c o r r e c t i o n for e a c h 4x4 s q u a r e + T h e t e r r a i n in t h i s r e g i o n (20
c c c c c c c C
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
450
ha=(hl(i)+h2(i)+h3(i)+h4(i)+h5(i)+h6(i)+h7(i)+h8(i)+h9(i)+ + hl0(i)+hll(i)+hl2(i)+hl3(i)+hl4(i)+hl5(i)+hl6(i))/16.0 p = 2000.0 tt=abs(0.5*gc*l.0E5*(d(i)*l.0E3)*l.6E7*(ha-h(k))**2 + /(r0**3-r0*p**2)) gt(k)=gt(k)+tt goto 4000
C
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c c c c c c
T o c a l c u l a t e the t e r r a i n c o r r e c t i o n for the r e g i o n (15
c
+++++++++++++++++++++++++++++++++++++++++++++++++++÷+++++++++++++++
500
r(1)=sqrt((xx(i)-1000-ex(k))**2+(yy(i)-1000-ey(k))**2)
+ + + + + +
Terrain correction program for regional gravity surveys r(2)=sqrt((xx(i)+1000-ex(k))**2+(yy(i)-1000-ey(k))**2) r(3)=sqrt((xx(i)-1000-ex(k))**2+(yy(i)+1000-ey(k))**2) r(4)=sqrt((xx(i)+1000-ex(k))**2+(yy(i)+1000-ey(k))**2) gg(1)=(hl(i)+h2(i)+h3(i)+h4(i))/4.0 gg(2)=(h5(i)+h6(i)+h7(i)+h8(i))/4.0 gg(3)=(h9(i)+hl0(i)+hll(i)+hl2(i))/4.0 gg(4)=(hl3(i)+hl4(i)+hl5(i)+hl6(i))/4.0 p = I000.0 do 550 1 = i, 4 tt=abs(0.5*gc*l.0E5*(d(i)*l.0E3)*4.0E6*(gg(1)-h(k))**2 + /(r(1)**3-r(1)*p**2)) gt(k)=gt(k)+tt continue 550 g o t o 4000 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c To c a l c u l a t e terrain corrections for the region (02 km. +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ do 3000 j = i, 16 if(ex(k).gt.(x(j)-500.0) .and. ex(k).it.(x(j)+500.0) .and. + ey(k).gt.(y(j)-500.0) .and .ey(k).lt.(y(j)+500.0)) goto 850 if(r(j).ge.2000) goto 800
969
+ + + + + + + +
+ + +
970
X.Q. MA and D. R. WATTS i f ( r ( j ) . g t . 5 0 0 . 0 .and. r(j).it.2000) goto 750 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ To c a l c u l a t e t e r r a i n c o r r e c t i o n in the near zone (0.5
c c c c c c c 75O
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
c c c c c c c 8O0
c c c c c c c c c c c c c 85O
+ + + + +
tt=gc*l.0E5*(d(i)*l.0E3)*(l.0E6)*(i/r(j)i/sqrt((b(j)-h(k))**2+r(j)**2)) gt(k)=gt(k)+tt g o t o 3000 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ To c a l c u l a t e t e r r a i n c o r r e c t i o n for the region (2
p = 500.0 tt=0.5*gc*l.0E5*(d(i)*l.0E3)*(l.0E6)*(b(j)-h(k))**2 + /(r(j)**3-r(j)*p**2) gt(k)=gt(k)+tt g o t o 3000 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ To c a l c u l a t e t e r r a i n c o r r e c t i o n s for the inner zone (r0<0.5km) by d i v i d i n g square into four t r i a n g u l a r prisms. zl, z2, z3 and z4 are h e i g h t s above the g r a v i t y s t a t i o n whate v e r the p o l a r i t i e s of zll,z22,z33 and z44 are. Next 8 lines (if) are u s e d to avoid u n d e r f l o w and overflow, w h i c h is c a u s e d by d i v i s i o n by 0.0, b e c a u s e for some stations z l = z 2 = z 3 = z 4 - i.e. the t e r r a i n is totally flat, if the triang u l a t i o n is a p p l i e d to this square, it may cause a r i t h e m a t i c problems, to a v o i d this, an a r t i f i c i a l small h e i g h t is a s s u m e d w h i c h s t a b i l i z e s the s o l u t i o n w i t h o u t i n t r o d u c i n g s i g n i f i c a n t errors.
+ + + + +
+ + + + + + + + + + +
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
i f ( z l l ( k ) . e q . h ( k ) ) zll(k) = z l l ( k ) + 0 . 0 5 i f ( z 2 2 ( k ) . e q . h ( k ) ) z22(k) = z 2 2 ( k ) + 0 . 0 6 if z33(k).eq.h(k)) z33(k) = z33(k)+0.07 if z44(k).eq.h(k)) z44(k) = z 4 4 ( k ) + 0 . 0 8 if z l l ( k ) . e q . z 2 2 ( k ) ) zll(k) = zll(k)+0.09 if z 2 2 ( k ) . e q . z 3 3 ( k ) ) z22(k) = z22(k)+0.10 if z 3 3 ( k ) . e q . z 4 4 ( k ) ) z33(k) = z33(k)+0.11 if z 4 4 ( k ) . e q . z l l ( k ) ) z44(k) = z44(k)+0.12 zl(k)=abs(zll(k)-h(k)) z2(k)=abs(z22(k)-h(k)) z3(k)=abs(z33(k)-h(k)) z4(k)=abs(z44(k)-h(k)) +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ To e s t a b l i s h a local c o o r d i n a t e system w i t h o r i g i n at the g r a v i t y + station, the c o o r d i n a t e s of 4 corners of the square are: + (xl, yl),(x2, y2),(x3, y3), (x4, y4). + +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ xl=x(j)-ex(k)+500.0 yl=y(j)-ey(k)-500.0 x2=xl y2=500.0+y(j)-ey(k) x3=x(j)-ex(k)-500.0 y3=y2 x4=x3 y4=yl next lines are to c o m p u t e the c o e f f i c i e n t s of 12 plane e q u a t i o n s d e f i n e d as z = ax + by. al=(zl(k)*y2-z2(k)*yl /(xl*y2-x2*yl bl=(xl*z2(k)-x2*zl(k) /(xl*y2-x2*yl a2=(z2(k)*y3-z3(k)*y2 /(x2*y3-x3*y2 b2=(x2*z3(k)-x3*z2(k) /(x2*y3-x3*y2 a 3 = ( z 3 ( k ) * y 4 - z 4 ( k ) * y 3 /(x3*y4-x4*y3 b3=(x3*z4(k)-x4*z3(k) /(x3*y4-x4*y3 a4=(z4(k)*yl-zl(k)*y4)/(x4*yl-xl*y4) b4=(x4*zl(k)-xl*z4(k))/(x4*yl-xl*y4) cl=(zl(k)*y2+z2(k)*yl)/(xl*y2-x2*yl) dl=(-xl*z2(k)-x2*zl(k))/(xl*y2-x2*yl) c2=(-zl(k)*y2-z2(k)*yl)/(xl*y2-x2*yl) d2=(xl*z2(k)+x2*zl(k))/(xl*y2-x2*yl) c3=(z2(k)*y3+z3(k)*y2)/(x2*y3-x3*y2) d3=(-x2*z3(k)-x3*z2(k))/(x2*y3-x3*y2) c4=(-z2(k)*y3-z3(k)*y2)/(x2*y3-x3*y2) d4=(x2*z3(k)+x3*z2(k))/(x2*y3-x3*y2) c5=(z3(k)*y4+z4(k)*y3)/(x3*y4-x4*y3)
Terrain correction program for regional gravity surveys d5=(-x3*z4(k)-x4*z3(k))/(x3*y4-x4*y3) c6=(-z3(k)*y4-z4(k)*y3)/(x3*y4-x4*y3) d6=(x3*z4(k)+x4*z3(k))/(x3*y4-x4*y3) cT=(z4(k)*yl+zl(k)*y4)/(x4*yl-xl*y4) d7=(-x4*zl(k)-xl*z4(k))/(x4*yl-xl*y4) c8=(-z4(k)*yl-zl(k)*y4)/(x4*yl-xl*y4) d8=(x4*zl(k)+xl*z4(k))/(x4*yl-xl*y4) next lines are to compute the coefficients of 8 line e q u a t i o n s d e f i n e d as y = kx kl=yl/xl k2=y2/xl k3=y2/x3 k4=y4/x3 kll=i/kl k22=i/k2 k33=i/k3 k44=I/k4 ll=-cl/dl 12=-c3/d3 13=-c5/d5 14=-c7/d7 iii=i/Ii 122=1/12 133=1/13 144=1/14 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ F o r m u l a listed in the paper are u s e d to calculate g r a v i t a t i o n a l a t t r a c t i o n s from 4 triangular prisms. The results are gl, g2 g3 and g4. Note that a factor G*rho has not been multiplied. Users are also w a r n e d that a lower limit and upper limit in i n t e g r a t i o n must be chosen properly to ensure the results are all positive. +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ call i n n e r z o n e ( k l , k 2 , a l , b l , x l , g l ) call innerzone(k33,k22,b2,a2,y2,g2) call i n n e r z o n e ( k 4 , k 3 , a 3 , b 3 , x 3 , g 3 ) call innerzone(k44,kll,b4,a4,y4,g4) g = gl+g2+g3+g4 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ The next routines are to calculate the gravity effects from four central p r i s m s (see text for locations) in order to consider a real s i t u a t i o n where the polarities of zzl, zz2, zz3 and zz4 could be arbitrary. +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ call innerzone (kl,ll,cl,dl,xl,wll) call innerzone (ll,k2,c2,d2,xl,wl2) wl = gl - (wll+wl2) call innerzone (122,k22,d3,c3,y2,w21) call innerzone (k33,122,d4,c4,y2,w22) w2 = g2 - (w21+w22) call innerzone (13,k3,c5,d5,x3,w31) call innerzone (k4,13,c6,d6,x3,w32) w3 = g3 - (w31+w32) call innerzone (k44,144,d7,c7,y4,w41) call innerzone (144,kll,d8,c8,y4,w42) w4 = g4 - (w41+w42) +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ To compare the four heights of p r i s m corners w i t h the station h e i g h t a n d to decide as w h i c h formula should be used to give a correct terrain correction. +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ htl=zll(k)-h(k) ht2=z22(k)-h(k) ht3=z33(k)-h(k) ht4=z44(k)-h(k) if(htl.gt 0 . a n d . h t 2 . g t . 0 . a n d . h t 3 . g t 0.and.ht4.1t.0) goto II00 if(htl.gt 0 . a n d . h t 2 . g t . 0 . a n d . h t 3 . 1 t 0.and.ht4.gt.0) goto 1200 if(htl.gt 0 . a n d . h t 2 . g t . 0 . a n d . h t 3 . 1 t 0.and.ht4.1t.0) goto 1300 if(htl.gt 0 . a n d . h t 2 . g t . 0 . a n d . h t 3 . g t 0.and.ht4.gt.0) goto 1400 if(htl.gt 0 . a n d . h t 2 . 1 t . 0 . a n d . h t 3 . g t 0.and.ht4.1t.0) goto 1500 if(htl.gt 0 . a n d . h t 2 . 1 t . 0 . a n d . h t 3 . g t 0.and.ht4.gt.0) goto 1600 if(htl.gt 0 . a n d . h t 2 . 1 t . 0 . a n d . h t 3 . 1 t 0.and.ht4.gt.0) goto 1700 if(htl.gt 0 . a n d . h t 2 . 1 t . 0 . a n d . h t 3 . 1 t 0.and.ht4.1t.0) goto 1800 if(htl.lt 0 . a n d . h t 2 . g t . 0 . a n d . h t 3 . g t 0.and.ht4.1t.0) goto 1900 if(htl.lt 0 . a n d . h t 2 . g t . 0 . a n d . h t 3 . 1 t 0.and.ht4.gt.0) goto 2000 if(htl.lt 0 . a n d . h t 2 . g t . 0 . a n d . h t 3 . 1 t 0.and.ht4.1t.0) goto 2100 if(htl.lt 0 . a n d . h t 2 . g t . 0 . a n d . h t 3 . g t 0.and.ht4.gt.0) goto 2200 if(htl.lt 0 . a n d . h t 2 . 1 t . 0 . a n d . h t 3 . g t 0.and.ht4.1t.0) goto 2300 if(htl.lt .0.and.ht2.1t.0.and.ht3.gt 0.and.ht4.gt.0) goto 2400 if(htl.lt .0.and.ht2.1t.0.and.ht3.1t 0.and.ht4.gt.0) goto 2500 if(htl.lt .0.and.ht2.1t.0.and.ht3.1t 0.and.ht4.1t.0) goto 2600
971
+ + + + + +
+ + + +
+ + +
972
X.Q. MA and D. R. WAtts ii00 1200 1300 1400 1500 1600 1700 1800 1900 2000 2100 2200 2300 2400 2500 2600 2800 3000 4000
5000
c c c
rt = g - w3 - w4 goto 2800 rt = g - w2 - w3 goto 2800 rt = g - w2 - w4 goto 2800 rt = g goto 2800 rt = g - wl - w2 - w3 - w4 goto 2800 rt = g - wl - w2 goto 2800 rt = g - wl - w3 goto 2800 rt = g - wl - w4 goto 2800 rt = g - wl - w3 goto 2800 rt = g - w3 - w4 goto 2800 rt = g - wl - w2 - w3 - w4 goto 2800 rt = g - wl - w4 goto 2800 rt = g - w2 - w3 goto 2800 rt = g - w2 - w4 goto 2800 rt = g - w3 - w4 goto 2800 rt = g tt=abs(gc*l.0E5*(d(i)*l.0E3)*rt) gt(k)=gt(k)+tt continue goto 4000 continue write(3,*), k, e x ( k ) , e y ( k ) , g t ( k ) print*,'Station = ',k,' continue close(l) close(2) stop end
gt
=
',gt(k)
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine innerzone (m,n,u,v,z,f) This subroutine evaluate the function f(m,n,u,v) given the values m,n,u,v,z. real u,v,m,n,p,gl,g2,g3,hl,h2,h3,t,f,z gl = n+sqrt(l+n**2) hl = m+sqrt(l+m**2) p = I/(l+v**2) g2 = n+u*v*p h2 = m+u*v*p t = I/sqrt(l+v**2) g3 = sqrt((n+u*v*p)**2+(l+u**2)*p-(u*v*p)**2) h3 = sqrt((m+u*v*p)**2+(l+u**2)*p-(u*v*p)**2) f = z*(log(gl/hl)-t*log((g2+g3)/(h2+h3))) return end ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Computers & Geosciences Vol. 20, No. 6, pp. 961-972, 1994 Copyright © 1994 Elsevier ScienceLtd 0098-3004(94)E0015-L Printed in Great Britain. All rights reserved 0098-3004/94 $7.00 + 0.00
TERRAIN CORRECTION PROGRAM FOR REGIONAL GRAVITY SURVEYS X. Q. MA* and D. R. WATTS Department of Geology and Applied Geology, University of Glasgow, Glasgow GI2 8QQ, U.K. (Received and accepted 3 November 1993)
Abstract--An algorithm is presented for the computation of gravity terrain corrections from a digitized grid representing topography. Different approximations for topography at various distances from the gravity station are used to make the computation more efficient. For the distant zone (30 < r < 50 km) topography is approximated using a vertical line mass distribution. In the intermediate zone (2 < r < 30 km) a segment of a hollow vertical cylinder is used. The near zone (0.5 < r < 2 km) is approximated using vertical rectangular prisms with sloping upper faces. A new method is used to compute the correction for the inner zone which is the digitized square where the gravity station is located. The gravity effect of topography in the inner zone is calculated from four traingular prisms with vertices made of the station location and the corners of the square. The computer program has been tested comprehensively with a regional gravity database from the Southern Uplands of Scotland in which all terrain corrections were carried out originally using the Hammer chart method. Key Words: Terrain correction, Gravity survey, FORTRAN 77.
INTRODUCTION Determining the terrain correction for a gravity anomaly is a tedious but potentially important task, especially in areas of rugged terrain. The calculation must be done efficiently as well as accurately and this may be expedited by using different approximations of topography depending on the distance from the gravity station. In this paper, we present a F O R T R A N - 7 7 computer program which utilizes various features of previous methods given by Bott (1959), Karlemo (1963), Nagy (1966), Ketelaar (1976), Blais and Ferland (1984), and Lopez (1990). A new method for the inner zone correction has been developed. Topography is represented by a digitized grid of equal squares of convenient size as shown in Figure 1. In the example considered here, topography is divided into a grid of 1 km squares and the digitized point is the center of each of the squares. If a square contains a gravity station, the four corners of that square are added to the database as well as the position and elevation of the gravity station. Note that the gravity station does not have to be centered within a square of the digitized grid. The zones are defined by the distance (r) of the digitized points from the gravity station being considered. The far distant zone (r > 50 km) is neglected. The distant zone is defined as (30 < r < 50 kin). The intermediate zone is between (2 < r < 30 km). The near zone is between
(0.5 < r < 2 km) and the inner zone is the digitized square with the gravity station inside. The inner, near, and a part of the intermediate zones are shown in Figure 1. DISTANT ZONE CONTRIBUTION To reduce the calculation time without losing accuracy, the vertical line mass Equation (1) is used to determine the contribution of the distant zone (30 < r < 50 km) topography (Bott, 1959): g = GpA
x / ( r 2 + z2) 3
(1) If h is less than about one-quarter of r a good approximation to (I) is given by the following simplified formula: g =
GpA hz 2 x r3.
(2)
where G is the gravitational constant, r is distance from station to center of square, h is difference in elevation between the square center and station, A is the area of square, and p is density of the topography. INTERMEDIATE ZONE CONTRIBUTION In this zone (2 < r < 30km) it is acceptable to approximate a prism as a segment of a hollow vertical cylinder (Bott, 1959). The vertical component of gravitational attraction from a segment of a hollow
*Present address: Department of Earth Sciences, University of Liverpool, P.O. Box 147, Liverpool L69 3BX U.K. 961
962
X.Q. MA and D. R. WATTS
cylinder with angle ~, inner radius r~, and outer radius r 2 is given by
g=GpA(r2-r~+x/r~+h2-
r~2+h2).
expression for becomes:
~ / 7 + h2 = r 1 + 2 r2].
(4)
If we now replace r 2 and r~ by r + p and r - p, ~ by A/2rp (A is area of the square, r is distance from station to center of square, and p is half length of square side), the terrain correction for a prism can be approximated as g
GpA h2 2 r(r2-p 2)
(5)
where G, A, h, are as given in (1). A considerable amount of computation time is saved if squares beyond a certain distance from a station are grouped and treated as a single square (Bott, 1959). In the area where 2 0 < r < 30km, sixteen squares are treated as one square (4 × 4 km). Where 15 < r < 20 km, four squares are treated as a single square (2 x 2 km), and where 2 < r < 15 km, a single square (1 x 1 km) is used to calculate the terrain correction using formula (5).
NEAR ZONE CONTRIBUTION
The near zone is shown as the stippled region in Figure 1. This area is represented by the rectangular prisms with the horizontal lower faces and sloping upper faces determined by the heights at the centers of the prisms and heights of neighboring locations (Ketelaar, 1976). The gravitational effect of a rectangular prism with the horizontal upper and lower faces and vertical sides is given by Nagy (1966). A modification by Ketelaar (1976) includes a top face with a constant slope towards the origin. The sides of a prism have a unit length D, the length of the unit digitized square (in this paper D = 1 km). The gravitational attraction of the prism, which has its vertical axis at x = iD and y =jK, is given as in (Ketelaar, 1976):
g=Gpfff~dv=Gp(1-cos~)DK(i,j),
terrain correction therefore
(3)
Expression (3) may be approximated further with the square root expanded in a power series and the terms with higher powers neglected,
( l_h2
the
INNER ZONE CONTRIBUTION
This zone is the square containing the gravity station. We develop a new method, which not only achieves the correction automatically, but provides acceptable accuracy. The basic idea is to construct four triangular prisms within the square containing the station with four elevation values at four corners of the square, and that of the station. Figure 2 shows the four triangular prisms with horizontal lower faces and sloping upper faces so constructed. The gravitational attraction of individual prisms can be obtained by integrating gravity over the volume of each triangular prism, and the effects of the four prisms are summed together, with the result added to the total terrain correction for that station. This method is a generalization of that suggested by Lopez (1990). Suppose that four elevations at the four corners of the square are all higher than the station elevation. Based on this assumption, the terrain effects of the four triangular prisms are considered individually as follows. Let us take the first triangular prism OP~P2Z2Z~ as an example to show how the terrain effect is obtained. The basic calculation is a simple integration over volume. The integration limits are shown in Figures 2A and 2B. The equations of lines OPt, OP2, and plane O Z l Z 2 are y = k~x, y =k2x, and z =a~x +b~y respectively. Therefore, the integration can be written as:
gl = Gp
dx
~x
dy
x/(x: +y2 + z2)3 (8)
(6)
where ~ is an angle showing the slope of the upper face [ ~ = a r c t g (h/r)]. K(i, j) is an eightterm algorithm which is a function of indices i and j only, that is K(i, j) depends uniquely on the position of the prism, and may be termed the position matrix (Ketelaar, 1976). A simplification has been made to the position matrix K(i, j). For prisms in this near zone, it is determined that K(i, j) is almost equal to (i 2 +j2)-05 (Ketelaar, 1987). The
Figure 1. Digitized topography and several zones with respect to gravity station at O. Inner zone is square containing station. Near zone is stippled and portion of intermediate zone is shown as white squares outside of stippled area.
963
Terrain correction program for regional gravity surveys ,
zs
'f
Y
P3 (x3, Y 3 ) ~
Y
_
P2
(x2, Y2)
_
X
74
P4
I
'%
PI
(x4, Y4
I I
(Xl, Yl ) P4
A
B
Figure 2. A--Perspective view of inner zone with four corners above stauon level. B--Plan view of inner zone. The result from the integration (8) and those for other three triangular prisms are given by: gl = G p x l f ( m ,
n, u, v)],, = k.... k~......... b,
(9)
g2 = G p Y 2 f ( m , n, u, U)]m=k'3,n=k~,u=b 2......
(10)
gs = G p x 3 f ( m , n, u, V)[m=k.... k3. . . . . . . . . b3
(1 1)
g4 = G p Y 4 f ( m , n, u, v)[,, = ka,,=ki,, =b,,~,=,4
(12)
where f ( m , n, u, v) = In
u/3
n + x/1 + n 2 m +x/l+m
~/(
1
2
X/1 +/32
u/3 )2
l + u2
n+l~v2+
n+l--~vZ
+l+v
m + l~-~v2 +
m + i--~v2
+l+vZ
u2/32
2
(l+v2) z
.In//2/32
(1 + V2)2 (13)
and zty2--z2yl a I =
a2 =
xIZ2--X2z 1 b I
x~y2--x2yl z2y 3 -- z3y 2
,
- -
xjy2--x2yl X223 -- X3Z 2
b2 -
x2y 3 -- x3y 2
x2y3
z3Y4-- z4y3 a 3 -=
x 3 z 4 - - X4 z3 ,
b3 -
,
Xay 4 -- x4y 3 z4Yl
--
x3Y4 -- x4y 3
z~y4
a 4 -=
,
b4
x4z~ -- x~ z4 --
x4y I -- xlY 4
kl = Y_L ,
k2 =
Y_! ,
x4Yl
k3 =
Y.-! ,
X1
XI
X3
I
1
1
ki=~,
k~=~,
,
-- x3y 2
k;=~,
-- Xly 4
k4 =
Y4 , X3
1
The terrain correction of the inner zone is obtained by summing four incremental contributions for the four triangular prisms. Thus we get for the total correction for the inner zone, g = gt +g2 +g3 + g 4 .
(14)
As stated, the value g is determined on the assumption that the four elevation values: Z1 P~, Z2P2, Z3P3, and Z4P 4, are all positive. However, in reality, these four positions in the square can be either higher or lower than the gravity station. One example is shown in Figure 3 in which Z1P~ is negative and the others remain positive. Within a triangular area OP~ P2, the triangular prism OQI P2 Z2 represents a hill whereas the triangular prism OQ1PLZ~ represents a valley. Although the vertical component of gravitational attraction at the gravity station from these two prisms are opposite in sign, in determining Bougour anomaly, the terrain corrections for both a hill and a valley must be added to the observed gravity. Therefore we need to only consider the absolute value of gravity. Because of symmetry, the volume of the prism OQ~ Pt Z~ is equal to that of the prism OQ~ Pl Z~, that is, integrating over the volume OQ~ P~ Z~ is equivalent to integrating over the volume OQ~ P~ Z l . Therefore, the terrain effect (w~) of two triangular prisms OQI PI Z'~ and OQ~ P2 22 can be obtained by subtracting the effect produced by the central prism OQ~ZIZ2 from the overall effect produced by the triangular prism OP~ PzZ2ZI. The same treatment should be applied to the triangular area OPI P4. To consider sixteen cases which may occur in a real situation (Fig. 4), the gravitational attractions (w2, w3, w4) from other three central prisms (OQ2Z2Z3, OQ3Z3Z4, and OQ4Z4ZL) also must be determined. The results from integration are given by: wl = gl - G p x l [ f ( m ,
k'4=k~. +f(m,
n, u, v)lm=, ..... ;. . . . . . . . . . h,
n, u, v)l,,_ ;.... k. . . . . '2.,'= a2]
(15)
964
X . Q . MA and D. R. WARTS
Z2 P2
(x2, Y2)
~'Y
Q2
(x3, Y3)
P3 ~ y - - k 3 x
~
/
Z4 P4
y=k2x
'"1~ X
6
Q]
y=14x
P1
P4 (X4' Y4)
Q4
(Xl, Yl)
A
B
Figure 3. A--Perspective view of inner zone with one corner at Pt below station level. B--Plan view of A.
where
w2 = gz -- G p y z [ f ( m , n, u, v)lm=rz,,,=k~,.=a3 ...... + f ( m , n, u, V)lm=k~,,,=t~,,,=a4. . . . . ]
(16)
Cl =
zlY2-J-z2y
1
--XIZ2--X2Z
,
d~
I
-
xly2--x2yl
x[y2-x2yl
w3 = g3 - G p x 3 [ f ( m , n, u, v )l,,,=t3,,,=l,3 . . . . 5.~=as x~z2-Fx2z~
--zzy2--z2y 1
+ f ( m , n, u, v)],,, = k. . . . t3......... a6]
(17)
C2 - -
W4 = g4 -- G p Y 4 [ f ( m, n, u. v)l,.=k'4,,,=r4,.=a7 ...... c3 + f ( m , n, u, v)lm=r.,. = k't,,,= as,,,= c,] Z3
Z4
g-w3-w 4
g - w 1 -W 2-W3 "W 4
(18)
~
x~yE-xzyl z2y3--~z3y2 x2Y3--x3Y2
,
d2
d3 =
x~y2-x2Y~
,
--X2z3--X3z2 x2y3--x3Y2
Z2
g'w2"w 3
g - w ] "W2
g-W 1 -W 3
g-W 1 -W2-W 3-W 4
g'W2*W 3
g-w2-w 4
~-W2-W 4
g
g - w 1 -W 3
g-W 1 -W4
g - W 1 -W 2
g-W 1 - w 4
g'w3"w4
g
Figure 4. Sixteen examples of possible terrain configuration in inner zone and corresponding equations for terrain corrections.
965
Terrain correction program for regional gravity surveys Table 1. Example of gravity station data file. N o is gravity station number; (x, y) and h are coordinates and elevations of gravity stations respectivley; z~, z~, z 3, and z4 are elevations at four corners of square containing station No
x (m)
y (m)
h (m)
z I (m)
z2 (m)
z3 (m)
z4 (m)
1
273350
625460
516.9
426.7
434.3
414.5
472.4
2
274640
627040
403.9
335.3
329.2
278.9
350.5
3
276230
628160
249.3
243.8
379.5
365.8
292.6
4
277500
628020
263.5
213.4
259.1
378.0
243.8
5
278650
628440
235.3
228.6
320.0
259.1
228.6
6
279920
627770
217.3
236.3
213.4
243.8
269.7
7
279540
629700
344.4
297.2
317.0
365.8
320.0
8
276990
625680
398.4
341.4
362.7
405.4
408.4
9
276180
626520
450.2
362.7
396.2
378.0
406.9
10
275370
625240
467.6
406.9
406.9
396.2
442.0
--z2y 3 -- z3y 2 C4 ~
d4 =
, x2y 3 -- x3y 2
x2Y3 -- x3Y 2
z3y 4 "b z4y 3
c5=
,
ds=
x3y 4 -- x4y 3
--z3y 4 -- z4y 3
x3Y4
d7
,
x4y I -- xly 4
x4Yl -- Xly 4
--z4y I -- 21y4
,
ds :
x4y I -- Xly 4
Cl
x4zl .-]- X12 4 X4Yl
C3
4' 1
t2' =
A F O R T R A N - 7 7 p r o g r a m (see the Appendix) was developed on a Sun S P A R C 4 W o r k s t a t i o n for the evaluation o f terrain corrections using the expressions as described. The p r o g r a m requires two i n p u t files; a gravity station file ( s t a t n . d a t ) and a terrain file ( b l o c k . d a t ) . The station file is illustrated in Table 1, consisting o f station numbers, grid coordinates (x, y ) (corresponding with the U K N a t i o n a l Grid), elevation o f the station, elevations at the four corners of a square (z~, z2, z3, and z4). The terrain file is illustrated in Table 2 a n d consists of block numbers, coordinates of block centers (x, y), average density of block, a n d sixteen elevation values of the
-- x4y 3
- - X 4 Z l - - X 1Z4 ,
= E'
FORTRAN-77 PROGRAM
d6 =
z 4 y I -.]- Z l y 4
c8=
,
-- x4y 3
X 3 Z 4 -~- X 4 Z 3
, x3y 4 -- x4y 3
C7 =
--X3Z 4 -- X4Z 3 x3Y4
C6 =
E q u a t i o n s for the sixteen cases are shown in Figure 4 a n d are used to calculate terrain corrections for the inner zone.
X 2 Z 3 At- X 3 Z 2
1
-- xly 4
C5
t3=
C7
ds' 1
t; =- t3'
1 t; = ?, .
Table 2. Example of terrain data file. centers; p is average density of block;
N o is ht-h~6
block number; (x, y) are coordinates of block are elevations of sixteen squares within block
NO
x (m)
y (m)
I00
226000
560000
p(g/crn3) 2.73
88.4
91.4
........
85.3
101
230000
560000
2.73
83.8
76.2
........
70.1
102
234000
560000
2.73
61.0
45.7
........
53.3
103
238000
560000
2.73
41.1
54.9
........
45.7
104
222000
564000
2.73
152.4
149.4 . . . . . . . .
152.4
105
226000
564000
2.73
121.9
96.0
........
106.7
106
230000
564000
2.73
99.1
121.9 . . . . . . . .
106.7
107
234000
564000
2.73
61.0
53.3
........
67.1
108
238000
564000
2.73
65.5
70.1
........
76.2
109
222000
568000
2.73
144.8
167.6
h I (m)
h 2 (m)
........
hi6
(In)
152.4
X. Q. MA and D. R. WATTS
966
Table 3. Comparison of terrain corrections carried out using our program with those done by Hammer zone chart method. No is gravity station number; (x, y) are coordinates of gravity stations; gtneWis terrain correction carried out by our program; gtom is terrain correction done by Hammer's method x (m)
y (m)
g~m (mg~)
gt~, (mg~)
e=or (%)
107
295280
634380
13.25
13.44
1.4
168
306530
631480
10.07
10.00
0.7
487
264240
601040
5.61
5.72
1.9
1061
314160
608940
5.18
5.17
0.2
1122
314690
611100
7.41
7.36
0.7
2159
233520
591080
5.24
5.58
6.0
No
centers of the individual squares. The order of the digitized squares within a block is shown in Figure 5. The program reads the station and terrain data and calculates the distance from the station for every square. The squares and blocks are classified according to their distance from the station so the appropriate algorithm can be applied. For the inner zone, many calculations and checks of elevation are required. Because the Equations (9)-(12) and (15)-(18) have similar patterns, a subroutine termed innerzone is designed so that the main program can call it many times, each with different coefficient values for the evaluation of the inner zone correction.
TESTING THE PROGRAM
The program was tested using the British Geological Survey gravity data set for the Southern Uplands of Scotland. This data set was selected becasue all of the stations have a terrain correction that was computed originally using Hammer's method (Hammer, 1939). For this purpose, 4526 gravity stations in the western Southern Uplands, which occupies 11,000 km 2, were used. The area was digitized with a total of 1 !,328
elevation points read from the Ordnance Survey (OS) map for 708 sixteen square blocks. In addition, 18,104 elevation points also were read from the OS map in order to calculate the terrain correction for the inner zone. The densities used for the terrain file, block.dat, were taken from Mansfield and Kennett (1963), Bott and Masson-Smith (1960), and Parslow and Randall (1973). The program was tested by comparing the new terrain corrections with those given in the British Geological Survey data computed using the Hammer chart method. It was determined that most of the new terrain corrections are within a few percent of the original values. We arbitrarily select six stations. The new and old terrain corrections are listed in Table 3. However, there were a few stations where new terrain corrections were substantially different from the original ones by up to a few mgals. To examine the problem, we computed the terrain correction for these stations using the Hammer chart method. In each situation, we determined that our Hammer method calculations produced values that agreed closely with those given by our program, indicating that some terrain corrections in the original data set are in error.
SUMMARY AND CONCLUSIONS 11
12
15
16
10
13
14
4
7
8
2
5
6
Figure 5. Terrain block illustrating order of digitization for file block.dat.
A program for gravity terrain correction has been developed for regional gravity surveys. The program is made more efficient by dividing the digitized terrain into different zones. Within each zone, different equations with approximations appropriate to the distance of the topography from the station are applied. A new method is used for the computation of the inner-zone correction, which involves the use of triangular prisms. The station is not restricted to the center of the inner square. The program was tested using gravity data from the Southern Uplands of Scotland. The results have shown that the computerized terrain corrections for most stations agree closely with the original corrections made by the Hammer chart method. Discrepancies were shown to be errors in the application of the Hammer method in the original data.
967
Terrain correction program for regional gravity surveys REFERENCES
Blais, J. A. R., and Ferland, R., 1984, Optimization in gravimetric terrain corrections: Can. Jour. Earth Sci., v. 21, no. 5, p. 505 515. Bott, M. H. P., 1959, The use of electronic digital computers for the evaluation of gravimetric terrain corrections: Geophysical Prospecting, v. 8, no. 1, p. 45 54. Bott, M. H. P., and Masson-Smith, D., 1960, A gravity survey of the Criffell granodiorite and the New Red Sandstone deposits near Dumfries: Proc. Yorkshire Geol. Society, v. 32, Pt. 3, no. 13, p. 317-332. Hammer, S., 1939, Terrain corrections for gravimeter stations: Geophysics, v. 4, no. 3, p. 184 194. Karlemo, B., 1963, Calculation of terrain corrections in gravity studies using the electric computer: Geoexploration, v. 1, no. 1, p. 56 66.
Ketelaar, A. C. R., 1976, A system for computer-calculation of the terrain correction: Geoexploration, v. 14, no. 1, p. 57~5. Ketelaar, A. C. R., 1987, Terrain correction for gravity measurements, using a digital terrain model (DTM): Geoexploration, v. 24, no. 2, p. 109-124. Lopez, H. R. B., 1990, F O R T R A N program for automatic terrain correction of gravity measurements: Computers & Geosciences, v. 16, no. 2, p. 237-244. Mansfield, J., and Kennett. P., 1963, A gravity survey of the Stranraer sedimentary basin: Proc. Yorkshire Geol. Society v. 34, Pt. 2, no. 9, p. 139 151. Nagy, D., 1966, The gravitational attraction of a right rectangular prism: Geophysics, v. 31, no. 2, p. 62-371. Parslow, G. R., and Randall, B. A. O., 1973, A gravity survey of the Cairnsmore of Fleet granite and its environs: Scott. Jour. Geology, v. 9, no. 3, p. 219-231.
APPENDIX Program TERRAIN ++++++++++++++++++++++++ ++++++++++++++++++++++++++++++++++÷++++++++ A F o r t r a n - 7 7 c o m p u t e r p r o g r a m for m a k i n g t e r r a i n c o r r e c t i o n s + for r e g i o n a l g r a v i t y surveys. D e s i g n e d and w r i t t e n by + X i n - Q u a n Ma & Doyle R. Watts, + D e p a r t m e n t of G e o l o g y & A p p l i e d Geology, + U n i v e r s i t y of G l a s g o w Scotland. + D e c e m b e r 1993 + +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
i n t e g e r mm, nn parameter(mm=708,nn=4526) real e x ( n n ) , e y ( n n ) , h ( n n ) , z l ( n n ) , z i ( n n ) , z 3 ( n n ) , z 4 ( n n ) real xx (mm) ,yy (mm) ,d (ram) ,hl (mm) ,h2 (mm) ,h3 (mm) ,h4 (mm) ,h5 (irma) real h 6 ( m m ) , h 7 ( m m ) , h S ( m m ) , h 9 ( m m ) , h l 0 ( m ~ ) , h l l ( m m ) , h l i ( m m ) real h l 3 ( m m ) , h l 4 ( m m ) , h l 5 ( m m ) , h l 6 ( m m ) , x ( 1 6 ) , y ( 1 6 ) , r ( 1 6 ) , b ( 1 6 ) real a l , a i , a 3 , a 4 , b l , b i , b 3 , b 4 , c l , c i , c 3 , c 4 , c 5 , c 6 , r t , c 7 , c 8 , d l , d 2 real k l , k 2 , k 3 , k 4 , k l l , k 2 2 , k 3 3 , k 4 4 , 1 1 , 1 2 , 1 3 , 1 4 , 1 1 1 , 1 2 2 , 1 3 3 , 1 4 4 real h t 4 , z l l ( n n ) , z 2 2 ( n n ) , z 3 3 ( n n ) , z 4 4 ( n n ) , g , g l , g 2 , g 3 , g 4 , g c , p real x l , x 2 , x 3 , x 4 , y l , y 2 , y 3 , y 4 , w l l , w 1 2 , w 2 I , w 2 2 , w 3 1 , w 3 2 , w 4 1 , w 4 2 real d 3 , d 4 , d 5 , d 6 , d 7 , d 8 , h t l , h t i , h t 3 , g g ( 4 ) , r 0 , g t ( n n ) , w l , w i , w 3 , w 4 open(unit=l,file='statn.dat') open(unit=2,file='block.dat') open(unit=3,file='outpt.dat') +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
ram n u m b e r of total blocks nn n u m b e r of total g r a v i t y stations ex x c o o r d i n a t e of g r a v i t y station (in metres) ey y c o o r d i n a t e of g r a v i t y station (in metres) h e l e v a t i o n of g r a v i t y station (in metres) z l l , z 2 2 , z 3 3 , z 4 4 e l e v a t i o n s of four corners of the square w h i c h c o n t a i n s the g r a v i t y station (in metres) xx x c o o r d i n a t e of a b l o c k centre (in metres) yy y c o o r d i n a t e of a b l o c k centre (in metres) hl-hl6 e l e v a t i o n s of 16 prisms w i t h i n the b l o c k (in metres) d a v e r a g e d e n s i t y of the b l o c k in g/cm**3 gt total t e r r a i n c o r r e c t i o n for one s t a t i o n in mgal gc is the g r a v i t a t i o n a l constant (gc=6.67x10**(-ll) N m**2/kg**2) +++++++++++++++++++++++++++++++++++÷+++++++++++++++++++++++++++++++
i00 150
print*,' print*,' ' + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ' print*,' + +' print*,' + Automatic Terrain Correction Program +' print*,' + +' print*,' ' + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ' print*,' ' gc = 6 . 6 7 3 2 E - I I do 150 k = i, nn read(l,100) e x ( k ) , e y ( k ) , h ( k ) , z l l ( k ) , z 2 2 ( k ) , z 3 3 ( k ) , z 4 4 ( k ) format(14x, if8,5f7.l,f5.2) continue do 300 i = i, m m read(2,200) x x ( i ) , y y ( i ) , d ( i ) , h l ( i ) , h i ( i ) , h 3 ( i ) , + h4(i),h5(i),h6(i) read(2,250) h 7 ( i ) , h 8 ( i ) , h 9 ( i ) , h l 0 ( i ) , h l l ( i ) , h l i ( i ) , + hl3(i),hl4(i),hl5(i),hl6(i)
+ + + + + + + + + + + + + +
968
X . Q . MA and D. R. WATTS 200 250 300
f o r m a t ( 5 x , 2 x , f7,2x, f7,2x, f4.2,6(2x, f6.1)) f o r m a t ( 1 0 ( 2 x , f6.1)) continue
C
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c c c c c c c c
k is i n d e x for the g r a v i t y s t a t i o n w h i l e i is for b l o c k number. T h e p r o g r a m p r o c e s s e s d a t a b l o c k by block; a b l o c k is c l a s s i f i e d b y its d i s t a n c e to the g r a v i t y s t a t i o n a c c o r d i n g to w h i c h d i f f e r e n t f o r m u l a s are applied. In this p r o g r a m the t e r r a i n is d i v i d e d into 7 zones a c c o r d i n g to d i s t a n c e r0: I. r 0 < 0 . 5 k m ; 2. 0.5
C
+ + + + + + + +
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
d o 5000 k = i, nn gt(k) = 0.0 d o 4000 i = i, m m b ( 1 ) = h l (i) b(2)=h2(i) b ( 3 ) = h 3 (i) b(4)=h4(i) b(5)=h5(i) b(6)=h6(i) b(7)=h7(i) b(8)=h8(i) b(9)=h9(i) b(10)=hl0(i) b(ll)=hll(i) b(12)=hl2(i) b(13)=hl3(i) b(14)=hl4(i) b(15)=hl5(i) b(16)=hl6(i) C
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c c
c a l c u l a t e the d i s t a n c e b e t w e e n g r a v i t y s t a t i o n a n d the b l o c k c e n t r e (xx(i),yy(i)):
C
(ex(k),ey(k))
+ +
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
r0=sqrt((xx(i)-ex(k))**2+(yy(i)-ey(k))**2) i f ( r 0 . 1 e . 1 5 0 0 0 ) g o t o 600 i f ( r 0 . 1 e . 2 0 0 0 0 ) g o t o 500 i f ( r 0 . 1 t . 3 0 0 0 0 ) g o t o 450 i f ( r 0 . 1 t . 5 0 0 0 0 ) g o t o 400 g o t o 4000 C
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c c c c c c c c
T o c a l c u l a t e the t e r r a i n c o r r e c t i o n for e a c h 4x4 k m s q u a r e T h e t e r r a i n in this r e g i o n (30
c
+ + + + + + + +
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
400
ha=(hl(i)+h2(i)+h3(i)+h4(i)+h5(i)+h6(i)+h7(i)+h8(i)+h9(i)+ + hl0(i)+hll(i)+hl2(i)+hl3(i)+hl4(i)+hl5(i)+hl6(i))/16.0 tt=abs(O.5*gc*l.
OE5*(d(i)*lOOO.O)*(l.6EO7)*(ha-h(k))**2/rO**3)
gt(k)=gt(k)+tt g o t o 4000 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ T o c a l c u l a t e the t e r r a i n c o r r e c t i o n for e a c h 4x4 s q u a r e + T h e t e r r a i n in t h i s r e g i o n (20
c c c c c c c C
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
450
ha=(hl(i)+h2(i)+h3(i)+h4(i)+h5(i)+h6(i)+h7(i)+h8(i)+h9(i)+ + hl0(i)+hll(i)+hl2(i)+hl3(i)+hl4(i)+hl5(i)+hl6(i))/16.0 p = 2000.0 tt=abs(0.5*gc*l.0E5*(d(i)*l.0E3)*l.6E7*(ha-h(k))**2 + /(r0**3-r0*p**2)) gt(k)=gt(k)+tt goto 4000
C
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c c c c c c
T o c a l c u l a t e the t e r r a i n c o r r e c t i o n for the r e g i o n (15
c
+++++++++++++++++++++++++++++++++++++++++++++++++++÷+++++++++++++++
500
r(1)=sqrt((xx(i)-1000-ex(k))**2+(yy(i)-1000-ey(k))**2)
+ + + + + +
Terrain correction program for regional gravity surveys r(2)=sqrt((xx(i)+1000-ex(k))**2+(yy(i)-1000-ey(k))**2) r(3)=sqrt((xx(i)-1000-ex(k))**2+(yy(i)+1000-ey(k))**2) r(4)=sqrt((xx(i)+1000-ex(k))**2+(yy(i)+1000-ey(k))**2) gg(1)=(hl(i)+h2(i)+h3(i)+h4(i))/4.0 gg(2)=(h5(i)+h6(i)+h7(i)+h8(i))/4.0 gg(3)=(h9(i)+hl0(i)+hll(i)+hl2(i))/4.0 gg(4)=(hl3(i)+hl4(i)+hl5(i)+hl6(i))/4.0 p = I000.0 do 550 1 = i, 4 tt=abs(0.5*gc*l.0E5*(d(i)*l.0E3)*4.0E6*(gg(1)-h(k))**2 + /(r(1)**3-r(1)*p**2)) gt(k)=gt(k)+tt continue 550 g o t o 4000 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c To c a l c u l a t e terrain corrections for the region (0
969
+ + + + + + + +
+ + +
970
X.Q. MA and D. R. WATTS i f ( r ( j ) . g t . 5 0 0 . 0 .and. r(j).it.2000) goto 750 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ To c a l c u l a t e t e r r a i n c o r r e c t i o n in the near zone (0.5
c c c c c c c 75O
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
c c c c c c c 8O0
c c c c c c c c c c c c c 85O
+ + + + +
tt=gc*l.0E5*(d(i)*l.0E3)*(l.0E6)*(i/r(j)i/sqrt((b(j)-h(k))**2+r(j)**2)) gt(k)=gt(k)+tt g o t o 3000 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ To c a l c u l a t e t e r r a i n c o r r e c t i o n for the region (2
p = 500.0 tt=0.5*gc*l.0E5*(d(i)*l.0E3)*(l.0E6)*(b(j)-h(k))**2 + /(r(j)**3-r(j)*p**2) gt(k)=gt(k)+tt g o t o 3000 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ To c a l c u l a t e t e r r a i n c o r r e c t i o n s for the inner zone (r0<0.5km) by d i v i d i n g square into four t r i a n g u l a r prisms. zl, z2, z3 and z4 are h e i g h t s above the g r a v i t y s t a t i o n whate v e r the p o l a r i t i e s of zll,z22,z33 and z44 are. Next 8 lines (if) are u s e d to avoid u n d e r f l o w and overflow, w h i c h is c a u s e d by d i v i s i o n by 0.0, b e c a u s e for some stations z l = z 2 = z 3 = z 4 - i.e. the t e r r a i n is totally flat, if the triang u l a t i o n is a p p l i e d to this square, it may cause a r i t h e m a t i c problems, to a v o i d this, an a r t i f i c i a l small h e i g h t is a s s u m e d w h i c h s t a b i l i z e s the s o l u t i o n w i t h o u t i n t r o d u c i n g s i g n i f i c a n t errors.
+ + + + +
+ + + + + + + + + + +
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
i f ( z l l ( k ) . e q . h ( k ) ) zll(k) = z l l ( k ) + 0 . 0 5 i f ( z 2 2 ( k ) . e q . h ( k ) ) z22(k) = z 2 2 ( k ) + 0 . 0 6 if z33(k).eq.h(k)) z33(k) = z33(k)+0.07 if z44(k).eq.h(k)) z44(k) = z 4 4 ( k ) + 0 . 0 8 if z l l ( k ) . e q . z 2 2 ( k ) ) zll(k) = zll(k)+0.09 if z 2 2 ( k ) . e q . z 3 3 ( k ) ) z22(k) = z22(k)+0.10 if z 3 3 ( k ) . e q . z 4 4 ( k ) ) z33(k) = z33(k)+0.11 if z 4 4 ( k ) . e q . z l l ( k ) ) z44(k) = z44(k)+0.12 zl(k)=abs(zll(k)-h(k)) z2(k)=abs(z22(k)-h(k)) z3(k)=abs(z33(k)-h(k)) z4(k)=abs(z44(k)-h(k)) +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ To e s t a b l i s h a local c o o r d i n a t e system w i t h o r i g i n at the g r a v i t y + station, the c o o r d i n a t e s of 4 corners of the square are: + (xl, yl),(x2, y2),(x3, y3), (x4, y4). + +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ xl=x(j)-ex(k)+500.0 yl=y(j)-ey(k)-500.0 x2=xl y2=500.0+y(j)-ey(k) x3=x(j)-ex(k)-500.0 y3=y2 x4=x3 y4=yl next lines are to c o m p u t e the c o e f f i c i e n t s of 12 plane e q u a t i o n s d e f i n e d as z = ax + by. al=(zl(k)*y2-z2(k)*yl /(xl*y2-x2*yl bl=(xl*z2(k)-x2*zl(k) /(xl*y2-x2*yl a2=(z2(k)*y3-z3(k)*y2 /(x2*y3-x3*y2 b2=(x2*z3(k)-x3*z2(k) /(x2*y3-x3*y2 a 3 = ( z 3 ( k ) * y 4 - z 4 ( k ) * y 3 /(x3*y4-x4*y3 b3=(x3*z4(k)-x4*z3(k) /(x3*y4-x4*y3 a4=(z4(k)*yl-zl(k)*y4)/(x4*yl-xl*y4) b4=(x4*zl(k)-xl*z4(k))/(x4*yl-xl*y4) cl=(zl(k)*y2+z2(k)*yl)/(xl*y2-x2*yl) dl=(-xl*z2(k)-x2*zl(k))/(xl*y2-x2*yl) c2=(-zl(k)*y2-z2(k)*yl)/(xl*y2-x2*yl) d2=(xl*z2(k)+x2*zl(k))/(xl*y2-x2*yl) c3=(z2(k)*y3+z3(k)*y2)/(x2*y3-x3*y2) d3=(-x2*z3(k)-x3*z2(k))/(x2*y3-x3*y2) c4=(-z2(k)*y3-z3(k)*y2)/(x2*y3-x3*y2) d4=(x2*z3(k)+x3*z2(k))/(x2*y3-x3*y2) c5=(z3(k)*y4+z4(k)*y3)/(x3*y4-x4*y3)
Terrain correction program for regional gravity surveys d5=(-x3*z4(k)-x4*z3(k))/(x3*y4-x4*y3) c6=(-z3(k)*y4-z4(k)*y3)/(x3*y4-x4*y3) d6=(x3*z4(k)+x4*z3(k))/(x3*y4-x4*y3) cT=(z4(k)*yl+zl(k)*y4)/(x4*yl-xl*y4) d7=(-x4*zl(k)-xl*z4(k))/(x4*yl-xl*y4) c8=(-z4(k)*yl-zl(k)*y4)/(x4*yl-xl*y4) d8=(x4*zl(k)+xl*z4(k))/(x4*yl-xl*y4) next lines are to compute the coefficients of 8 line e q u a t i o n s d e f i n e d as y = kx kl=yl/xl k2=y2/xl k3=y2/x3 k4=y4/x3 kll=i/kl k22=i/k2 k33=i/k3 k44=I/k4 ll=-cl/dl 12=-c3/d3 13=-c5/d5 14=-c7/d7 iii=i/Ii 122=1/12 133=1/13 144=1/14 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ F o r m u l a listed in the paper are u s e d to calculate g r a v i t a t i o n a l a t t r a c t i o n s from 4 triangular prisms. The results are gl, g2 g3 and g4. Note that a factor G*rho has not been multiplied. Users are also w a r n e d that a lower limit and upper limit in i n t e g r a t i o n must be chosen properly to ensure the results are all positive. +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ call i n n e r z o n e ( k l , k 2 , a l , b l , x l , g l ) call innerzone(k33,k22,b2,a2,y2,g2) call i n n e r z o n e ( k 4 , k 3 , a 3 , b 3 , x 3 , g 3 ) call innerzone(k44,kll,b4,a4,y4,g4) g = gl+g2+g3+g4 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ The next routines are to calculate the gravity effects from four central p r i s m s (see text for locations) in order to consider a real s i t u a t i o n where the polarities of zzl, zz2, zz3 and zz4 could be arbitrary. +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ call innerzone (kl,ll,cl,dl,xl,wll) call innerzone (ll,k2,c2,d2,xl,wl2) wl = gl - (wll+wl2) call innerzone (122,k22,d3,c3,y2,w21) call innerzone (k33,122,d4,c4,y2,w22) w2 = g2 - (w21+w22) call innerzone (13,k3,c5,d5,x3,w31) call innerzone (k4,13,c6,d6,x3,w32) w3 = g3 - (w31+w32) call innerzone (k44,144,d7,c7,y4,w41) call innerzone (144,kll,d8,c8,y4,w42) w4 = g4 - (w41+w42) +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ To compare the four heights of p r i s m corners w i t h the station h e i g h t a n d to decide as w h i c h formula should be used to give a correct terrain correction. +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ htl=zll(k)-h(k) ht2=z22(k)-h(k) ht3=z33(k)-h(k) ht4=z44(k)-h(k) if(htl.gt 0 . a n d . h t 2 . g t . 0 . a n d . h t 3 . g t 0.and.ht4.1t.0) goto II00 if(htl.gt 0 . a n d . h t 2 . g t . 0 . a n d . h t 3 . 1 t 0.and.ht4.gt.0) goto 1200 if(htl.gt 0 . a n d . h t 2 . g t . 0 . a n d . h t 3 . 1 t 0.and.ht4.1t.0) goto 1300 if(htl.gt 0 . a n d . h t 2 . g t . 0 . a n d . h t 3 . g t 0.and.ht4.gt.0) goto 1400 if(htl.gt 0 . a n d . h t 2 . 1 t . 0 . a n d . h t 3 . g t 0.and.ht4.1t.0) goto 1500 if(htl.gt 0 . a n d . h t 2 . 1 t . 0 . a n d . h t 3 . g t 0.and.ht4.gt.0) goto 1600 if(htl.gt 0 . a n d . h t 2 . 1 t . 0 . a n d . h t 3 . 1 t 0.and.ht4.gt.0) goto 1700 if(htl.gt 0 . a n d . h t 2 . 1 t . 0 . a n d . h t 3 . 1 t 0.and.ht4.1t.0) goto 1800 if(htl.lt 0 . a n d . h t 2 . g t . 0 . a n d . h t 3 . g t 0.and.ht4.1t.0) goto 1900 if(htl.lt 0 . a n d . h t 2 . g t . 0 . a n d . h t 3 . 1 t 0.and.ht4.gt.0) goto 2000 if(htl.lt 0 . a n d . h t 2 . g t . 0 . a n d . h t 3 . 1 t 0.and.ht4.1t.0) goto 2100 if(htl.lt 0 . a n d . h t 2 . g t . 0 . a n d . h t 3 . g t 0.and.ht4.gt.0) goto 2200 if(htl.lt 0 . a n d . h t 2 . 1 t . 0 . a n d . h t 3 . g t 0.and.ht4.1t.0) goto 2300 if(htl.lt .0.and.ht2.1t.0.and.ht3.gt 0.and.ht4.gt.0) goto 2400 if(htl.lt .0.and.ht2.1t.0.and.ht3.1t 0.and.ht4.gt.0) goto 2500 if(htl.lt .0.and.ht2.1t.0.and.ht3.1t 0.and.ht4.1t.0) goto 2600
971
+ + + + + +
+ + + +
+ + +
972
X.Q. MA and D. R. WAtts ii00 1200 1300 1400 1500 1600 1700 1800 1900 2000 2100 2200 2300 2400 2500 2600 2800 3000 4000
5000
c c c
rt = g - w3 - w4 goto 2800 rt = g - w2 - w3 goto 2800 rt = g - w2 - w4 goto 2800 rt = g goto 2800 rt = g - wl - w2 - w3 - w4 goto 2800 rt = g - wl - w2 goto 2800 rt = g - wl - w3 goto 2800 rt = g - wl - w4 goto 2800 rt = g - wl - w3 goto 2800 rt = g - w3 - w4 goto 2800 rt = g - wl - w2 - w3 - w4 goto 2800 rt = g - wl - w4 goto 2800 rt = g - w2 - w3 goto 2800 rt = g - w2 - w4 goto 2800 rt = g - w3 - w4 goto 2800 rt = g tt=abs(gc*l.0E5*(d(i)*l.0E3)*rt) gt(k)=gt(k)+tt continue goto 4000 continue write(3,*), k, e x ( k ) , e y ( k ) , g t ( k ) print*,'Station = ',k,' continue close(l) close(2) stop end
gt
=
',gt(k)
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine innerzone (m,n,u,v,z,f) This subroutine evaluate the function f(m,n,u,v) given the values m,n,u,v,z. real u,v,m,n,p,gl,g2,g3,hl,h2,h3,t,f,z gl = n+sqrt(l+n**2) hl = m+sqrt(l+m**2) p = I/(l+v**2) g2 = n+u*v*p h2 = m+u*v*p t = I/sqrt(l+v**2) g3 = sqrt((n+u*v*p)**2+(l+u**2)*p-(u*v*p)**2) h3 = sqrt((m+u*v*p)**2+(l+u**2)*p-(u*v*p)**2) f = z*(log(gl/hl)-t*log((g2+g3)/(h2+h3))) return end ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++