Three-dimensional gravity modeling in all space

30 Pages·0.84 MB·English

by  Xiong L.I., Michel Chouteau

THREE-DIMENSIONALGRAVITYMODELINGINALLSPACE
XIONGLI
BHPResearch–MelbourneLaboratories,P.O.Box264,ClaytonSouth,Victoria3169,Australia
E-mail:[email protected]
MICHELCHOUTEAU
DépartementdeGénieMinéral,ÉcolePolytechniquedeMontréalC.P.6079,succ.Centre-ville,
Montréal,QuébecH3C3A7,Canada
E-mail:[email protected]
Abstract. We review available analytical algorithms for the gravity effect and gravity gradients
especiallytheverticalgravitygradientduetoarightrectangularprism,arightpolygonalprism,and
a polyhedron. The emphasis is placed on an investigation of validity, consistency, and especially
singularitiesof different algorithms, which havebeen traditionallyproposed for calculationof the
gravityeffectonground(oroutsideanomalousbodies),whentheyareappliedtoallpointsinspace.
Theroundingerrorduetothecomputerfloatingpointprecisionisestimated.Thegravityeffectand
verticalgradientofgravityinthreedimensionscausedbyacubicmodelarecalculatedbydifferent
typesofalgorithms.Thereliabilityofalgorithmsforthecalculationofgravityofarightpolygonal
prismandapolyhedronisfurtherverifiedbyusingaregularpolygonalprismapproximatingavertical
cylinderandaregularpolyhedronapproximatingasphere,respectively.ByhighlightingHaáz-Jung-
PlouffandOkabe-Steiner-Zilahi-Sebess’formulaeforarightrectangularprism,Plouff’salgorithm
for aright polygonal prism, and Götzeand Lahmeyer’salgorithm for apolyhedron andremoving
their singularities, we demonstrate that these formulae and algorithms can be used to model the
gravityanomalyanditsverticalgradientatallpossiblecomputationpositions.
Keywords: three-dimensional body, gravity,gravitygradient,forwardmodeling,allspace,singu-
larity,roundingerror
1. Introduction
Gravitymetersmeasuretheverticalcomponentofthegravitational attraction,sim-
ply and normally called the gravity and denoted by g. For convenience, we adopt
two Cartesian systems of coordinates: .(cid:24);(cid:17);(cid:16)/for body coordinates and .x;y;z/
for fieldcoordinates. Thevertical (cid:16) axis and z axis are taken to be positive down-
ward and the (cid:24) and (cid:17) axes and the x and y axes are arranged into a right-handed
system. The gravity due to a three-dimensional body with density (cid:26).(cid:24);(cid:17);(cid:16)/ and
arbitrary shapeobserved atpoint.x;y;z/is
Z Z Z
z(cid:0)(cid:16)
g.x;y;z/ D(cid:0)G (cid:26).(cid:24);(cid:17);(cid:16)/ d(cid:24)d(cid:17)d(cid:16) (1)
r3
SurveysinGeophysics 19: 339–368,1998.
©1998KluwerAcademicPublishers. PrintedintheNetherlands.
340 XIONGLIANDMICHELCHOUTEAU
whereGistheuniversalgravitational constant, and
p
r D .x (cid:0)(cid:24)/2C.y (cid:0)(cid:17)/2C.z(cid:0)(cid:16)/2 (2)
The simple formula for sphere, cylinder, and other geometrical models can
serve very well to give the approximate magnitudes of calculated gravity effects
forcomparison withtheobserved values. Thereare, however, manyinterpretation
problems in which more detailed calculations are useful. Such problems arise,
for instance, when it is desired to determine the gravity effect from a geological
situation which is partially defined by known stratigraphy and drilling and where
thegravityexpressioniscomplexbutwelldefined.Three-dimensionalmethodsare
preferred, because in general, the anomalies are irregular or do not have a strong
elongation inonedirection.
There are a number of ways to solve three-dimensional density problems. The
popular three categories according to the basic geometrical model used are the
stack of right rectangular prisms, the stack of polygonal prisms, and the polyhe-
dron. Traditionally, amethod has been developed for surface gravity observations
where the computation point is generally outside the body. In the calculation of
thegravityeffectinaboreholeoratunnel,theformulaisespecially requiredtobe
applicable to all possible computation positions, outside, on the boundary of, and
inside the three-dimensional body. An analytical solution is preferred in compar-
ison with the surface gravity, because a numerical quadrature solution, generally,
cannotproduceareasonably accurate gravityfieldinsideornearthebody.
In Sections 2, 3 and 4 we review the developments of analytical algorithms
for the gravity effect due to a right rectangular prism, a right polygonal prism,
and a polyhedron, and explain the techniques to remove singularities when the
computation point is inside or on the boundary of a prism or a polyhedron. Some
cautionsinprogrammingaresuggested. Wehavedesignedastandardcubicmodel
combining a three-dimensional observation array so as to check the correctness
and consistency of different algorithms. Plouff’s (1975, 1976) algorithm for the
calculation of gravity of a right polygonal prism is further tested by using an
approximation of a 16-side right regular polygonal prism to a vertical cylinder.
SoisGötzeandLahmeyer’s(1988)algorithmforapolyhedronbyusinga24-facet
regularpolyhedron approximating asphere.
Section5isdevotedtothedescriptionoftheroundingerrorduetothecomputer
floating point precision. We use a cubic model to examine the operational ranges
ofsingleanddoubleprecision codes.
Section 6 is a general discussion beyond the calculation of the vertical com-
ponent of the gravitational attraction, including the calculation of the gravity gra-
dient tensor. In Section 7, we present distributions of the vertical gravity gradient
in three-dimensional space caused by the standard cubic model and the regular
polyhedron.
THREE-DIMENSIONALGRAVITYMODELINGINALLSPACE 341
X
(x, y, z)
O Y
(x , h , z )
1 1 1
(x , h , z )
2 2 2
Z
Figure1.Therightrectangularprismmodel.
2. TheRightRectangularPrism
Acollectionofrectangularprismsprovidesasimplewaytoapproximateavolume
of mass. If small enough, each prism can be assumed to have constant density.
Then by the principle of superposition, the gravitational anomaly of the body at
anypointcouldbeapproximated bysummingtheeffectsofalltheprisms.
According to Equation (1), a rectangular parallelopiped, which sides of the
prism are parallel to x, y, z axes, with uniform density (cid:26) and with dimensions
described by the limits (cid:24) (cid:20) (cid:24) (cid:20) (cid:24) ;(cid:17) (cid:20) (cid:17) (cid:20) (cid:17) , and (cid:16) (cid:20) (cid:16) (cid:20) (cid:16) has avertical
1 2 1 2 1 2
attraction atpoint.x;y;z/(Figure1)givenby
Z Z Z
(cid:16)2 (cid:17)2 (cid:24)2 z(cid:0)(cid:16)
g D(cid:0)G(cid:26) d(cid:24)d(cid:17)d(cid:16) (3)
r3
(cid:16)1 (cid:17)1 (cid:24)1
342 XIONGLIANDMICHELCHOUTEAU
2.1. DIFFERENT DERIVATIONS
Theexpression forEquation(3)hasbeenderived bynumerous researchers. These
analyticalformulaehaveseveraldifferentforms.Theearliestonemightbederived
by Everest (1830, p. 94–97). He calculated closed-form expressions for the hor-
izontal and vertical gravitational effects of a rectangular parallelopiped and used
theseequations toestimatethetopographic deflection oftheplumbbobduetothe
Satpura Range in India. MacMillan (1930, p. 72–79) gave an analytical formula
for computing the potential of this body and discussed the derivation of the three
components of gravitational attraction. Another solution was presented by Ansel
(1936) both in a closed-form and in the form of an infinite series. He pointed out
that for numerical evaluation, a combined use of these two forms is preferred. In
fact,hisanalyticalformularepresentsthegravityeffectofarectangularparallelop-
ipedinthecaseofy(cid:0)(cid:17) D0.Forthenormalcases,toobtainthegravityeffectof
1
anyrectangular parallelopiped weneedtoobtain thedifference oftworectangular
parallelopipeds withacommonsideinthexz-plane. Kolbenheyer (1963, equation
(11))expressedtheverticalcomponentoftheattractionofarightprismaticalbody
(a right polygonal prism) as the potential difference between the top and bottom
boundary surfaces. Applying his solution to the right rectangular prism, a closed
expression for gravity could be derived. Motivated to access errors in computa-
tionofgravity terraincorrections, Müller(1963, p.20)gaveananalytical solution
which is valid for all possible positions of the right rectangular prism with square
topandbottomsurfaces.
ItisfoundthatSorokin(1951,p.370,equation(426))publishedanunequivocal
solution tothisproblem. Hisformulacanbewritteninto
X2 X2 X2
g D (cid:0)G(cid:26) (cid:22)
ijk
(cid:20) iD1 jD1 kD1 (cid:21)
(cid:0) (cid:1) (cid:0) (cid:1)
z r
(cid:2) x ln y Cr Cy ln x Cr Cz arctan k ijk (4)
i j ijk j i ijk k
x y
i j
where
x Dx (cid:0)(cid:24) ; y Dy (cid:0)(cid:17) ; z Dz(cid:0)(cid:16) (5)
i i j j k k
q
r D x2 Cy2 Cz2 (6)
ijk i j k
and
(cid:22) D.(cid:0)1/i.(cid:0)1/j.(cid:0)1/k (7)
ijk
Severalformulae similarto(4) werederived byother researchers. Their differ-
ence is only in the last term of (4). Haáz (1953, p. 62) applied Euler’s theorem of
THREE-DIMENSIONALGRAVITYMODELINGINALLSPACE 343
homogeneous functions to the second derivative of the potential of a prism, thus
obtaining the first derivative and the potential itself without integration. Hisresult
fortheverticalcomponent oftheattraction ofarightrectangular prismis
X2 X2 X2
g D (cid:0)G(cid:26) (cid:22)
ijk
(cid:20) iD1 jD1 kD1 (cid:21)
(cid:0) (cid:1) (cid:0) (cid:1)
x y
(cid:2) x ln y Cr Cy ln x Cr (cid:0)z arctan i j (8)
i j ijk j i ijk k
z r
k ijk
Jung(1961) andPlouff(1966, 1975, 1976) derivedindependently Equation (8)
whichPlouffusedforcomparisonwithterraincorrection approximation formulae.
The equivalence of Equation (4) and Equation (8) can be easily demonstrated.
It is well known arctan.a/ D (cid:25)=2 (cid:0) arctan.1=a/, the (cid:25)=2 terms cancel by the
summations.
Nagy (1966) carried out the integration (3) in a different way from Sorokin
(1951), inwhichheusedarcsineratherthanarctangent expressions, andobtained
"
X2 X2 X2 (cid:0) (cid:1) (cid:0) (cid:1)
g D(cid:0)G(cid:26) (cid:22) x ln y Cr Cy ln x Cr
ijk i j ijk j i ijk
iD1 jD1 kD1 #
y2 Cz2Cy r
(cid:0)z arcsin (cid:0) j k(cid:1)q j ijk (9)
k
y Cr y2 Cz2
j ijk j k
Okabe (1979, equation (46)) and Steiner and Zilahi-Sebess (1988) presented
another closed-form formula which, like Sorokin’s and Haáz, Jung and Plouff’s,
hasanarctangent term.Itis
"
X2 X2 X2 (cid:0) (cid:1) (cid:0) (cid:1)
g D(cid:0)G(cid:26) (cid:22) x ln y Cr Cy ln x Cr
ijk i j ijk j i ijk
iD1 jD1 kD1 #
x Cy Cr
C2z arctan i j ijk (10)
k
z
k
2.2. COMPUTATIONAL ASPECTS
Nagy (1966) wrote, when either the x or y or both axes is crossed (that the x axis
is crossed means (cid:24) (cid:20) x coordinate of a computation point (cid:20) (cid:24) ), the integration
1 2
must be carried out from the lower limit to the axis, then from the axis to the
upper limit, the sum of these integration giving the required effect. This means
that for an interior point, the prism needs to be divided into four subprisms. This
subdivision processing and computational difficulty are due to the use of arcsines
344 XIONGLIANDMICHELCHOUTEAU
in the expression. It can be avoided by the use of arctangents, e.g., by Sorokin,
Haáz-Jung-Plouff, andOkabe-Steiner- Zilahi-Sebess’formulae.
Either Sorokin’s formula (Equation (4)) orHaáz-Jung-Plouff’s formula (Equa-
tion (8)) or Okabe-Steiner-Zilahi-Sebess’ formula (Equation (10)) contains 24
terms: 16 logarithms, and 8arctangents. Their difference is only in the arctangent
terms. As a cautionary note, the function ATAN2.a;b/ instead of ATAN(a/b) is
recommended touseinFORTRANprogramming, whenonecomputestheinverse
tangent functions in the formulae. Note that the function ATAN2 accepts two ar-
gumentswhicharetakentobe,respectively, theordinateandabscissaofapointin
anyofthefourquadrants.Withthisadditionalinformation,ATAN2candistinguish
the four quadrants and returns an angle in the range from (cid:0)(cid:25) to (cid:25), while ATAN
function providesaresultintherangeof(cid:0)(cid:25)/2to(cid:25)/2only.
We have designed a cubic model with a unit density (1 g/cm3) to check the
validity and consistency of Haáz-Jung-Plouff and Okabe-Steiner-Zilahi-Sebess’
formulae.Thelowerlimitsofthecubeare(cid:24) D(cid:17) D (cid:16) D(cid:0)10m,theupperlimits
1 1 1
(cid:24) D (cid:17) D (cid:16) D 10 m.The computation points are placed on athree-dimensional
2 2 2
grid with a lower bound corner ((cid:0)20, (cid:0)20, (cid:0)20), an upper bound corner (20, 20,
20), and a grid interval (the same in all three directions) 2 m. This arrangement
produces an observation array composed of9261 points. Thismodel isnamed the
standard cubicmodel.
The calculation results show (Figure 2) that the gravity is zero at the z D 0
plane, is antisymmetrical about the z D 0 plane, is symmetrical about the x D 0
and the y D 0 planes, and reaches an absolute maximum (346.426 (cid:22)Gal) at the
center of the upper or lower boundary surface of the prism. The magnitude of
gravity decreases smoothly, when the computation point is away from the body.
The FORTRAN single precision codes result in a maximum absolute difference
between the two formulae 0.00056 (cid:22)Gal at point ((cid:0)18, (cid:0)8, (cid:0)18) where the cal-
culated gravity is 51.203 (cid:22)Gal. Translating into the relative difference it is about
11permillion.Themaximumrelativedifferenceis87.5%whichoccursat(18,10,
0)andatwhich thecalculated gravity fields bytwoformulae are0.0001018 (cid:22)Gal
and0.0000127(cid:22)Gal,respectively.Thesedifferencesrepresentthelimitationofthe
floatingpointprecision. Therefore,thevalueof87.5%hereismeaningless. Bythe
way,allthecalculations reported inthispaperweredoneonaSUNSparc5.
Inconclusion, althoughHaáz-Jung-Plouff’s formula(Equation(8))andOkabe-
Steiner-Zilahi-Sebess’ formula (Equation (10)) do not seem to agree, itis verified
thattheyareidenticalandvalidoutside,inside,ontheboundaryof,andevenatthe
verticesoftherightrectangular prism.
Itisworthtopointoutthatbychangingthesignofx ;y ;z inEquation(5),in
i j k
otherwords,bydenoting
x D(cid:24) (cid:0)x; y D(cid:17) (cid:0)y; z D(cid:16) (cid:0)z (11)
i i j j k k
without any other modifications, the gravity calculated by Equations (4), (8), and
(10) will not change, will remain the same magnitude and the same sign. This
THREE-DIMENSIONALGRAVITYMODELINGINALLSPACE 345
300
−20
200
−10
100
m) 0
Z( 0
10
−100
20
20 −200
10 −20
−10
0
0 −300
−10
10
Y(m) −20 20
X(m)
Figure 2. The gravity effect in three dimensions caused by the standard cubic model with a unit
density (1 g/cm3), and with the lower limits (cid:24) D (cid:17) D (cid:16) D (cid:0)10 m and the upper limits
1 1 1
(cid:24) D (cid:17) D (cid:16) D 10 m. The gravity is zero at the z D 0 plane, is antisymmetrical about the
2 2 2
z D 0plane,issymmetricalaboutthex D 0andy D 0planes,andreachesanabsolutemaximum
(346.426(cid:22)Gal)atthecenterofthetoporbottomsurfaceofthecube.
fact can be explained from the structure of formulae (Goodacre, 1973). Since the
integrandinEquation(3)isanevenfunctionofx(cid:0)(cid:24) andy(cid:0)(cid:17)andanoddfunction
ofz(cid:0)(cid:16),theintegratedresultshouldbeanoddfunctionofx(cid:0)(cid:24) andy(cid:0)(cid:17)(butthe
product oftwooddfunctions isanevenfunction incaseofasimultaneous change
of their signs) and an even function of z (cid:0) (cid:16). At a first glance, neither Equation
(4)nor(8)nor(10)seemstosatisfytheseconditions, but,forexample,(8)maybe
rewrittenas
" 0 s 1
X2 X2 X2 y y2
g D(cid:0)G(cid:26) (cid:22) x ln@q j C 1C j A
iD10jD1 kD1 ijk si xi2 Cz1k2 xi2 Czk2 #
x x2 x y
Cy ln@q i C 1C i A(cid:0)z arctan i j (12)
j y2 Cz2 yj2 Czk2 k zkrijk
j k
346 XIONGLIANDMICHELCHOUTEAU
since
(cid:18) q (cid:19)
X2 X2 X2
(cid:22) x ln x2 Cz2 D0 (13)
ijk i i k
iD1 jD1 kD1
and
X2 X2 X2 (cid:16) q (cid:17)
(cid:22) y ln y2Cz2 D0 (14)
ijk j j k
iD1 jD1 kD1
Because
0 s 1
a a2 a
ln@p C 1C A Darcsin p (15)
b2 Cc2 b2Cc2 b2Cc2
itisclear that (12) has the desired symmetry property. Equation (9)does not have
this property because of the form of the argument in the arcsine term. In practice,
westilluse(8)ratherthan(12)because(8)isequivalent toandsimplerthan(12).
3. TheRightPolygonalPrism
Althoughconceptuallystraightforward,therightrectangularprismapproachwould
be cumbersome in practice. Geological bodies are often difficult to be modeled
withrectangularblocks.Moreover,thecomputationdoesnottakeadvantageofthe
fact that if the densities of neighboring prisms are identical, there is no need to
includetheirmutualinterface inthecalculation.
Talwani and Ewing (1960) and Talwani (1973) described a more practical
method. Their technique approximates a body by a stack of infinitely thin, hor-
izontal laminas. The shape of each lamina is approximated by a polygon. The
key formula in this method is to calculate the gravity effect of a horizontal poly-
gon. Grant and West(1965, equations (10–16)) and Golizdra (1981, equation (4))
improvedTalwaniandEwing’sformulabyinvolvinganevaluationoflesstranscen-
dental functions. Integration over the vertical coordinate is done with numerical
quadrature technique. The shorter the height of observation point compared to
the thickness represented by each lamina, the more inaccurate the result of the
numerical integration will be. The gravity field within or near a lamina cannot be
accurately determined bythelaminamethod(Plouff,1976).
Plouff (1975, 1976) took this representation one step further. He derived an
analytical formula forthegravity effect ofalayer offinitethickness, withvertical
sides and with top and bottom surfaces approximated by polygons. Analogous to
themethodofTalwaniandEwing(1960),thesepolygonallayerscanbestackedon
THREE-DIMENSIONALGRAVITYMODELINGINALLSPACE 347
X
X
O D y
Y (x, y)
1 1
R
k1 b
r r 1 d1 D s D x
1
dv
z
1 d
P 2
z (x, y)
2 r 2 2
A 2
b
2
Y
O
Z
(a) (b)
Figure3.Geometricalelementsinvolvedincalculationofgravityanomaliescausedbyarightpolyg-
onalprism.TheobservationpointislocatedattheoriginoftheCartesiancoordinatesystem.(a)a
polygonalprism;(b)planviewofonesideofprism.(AfterPlouff,1976).
top of one another in order to approximate three-dimensional bodies of arbitrary
shape.
The geometrical symbols used in Plouff’s formula are shown in Figure 3. For
the ith edge, the subscript k in R refers to the endpoint (or corner) of that edge,
kj
andthesubscriptj thefaceoftheprism.k D1indicatesthefirstendpointfoundin
theclockwiseprogression,k D2thelatterofthetwoendpoints.j D1referstothe
faceclosertothecomputationpoint,j D 2themoredistantface.Thegravityeffect
ofapolygonal prism areexpressed intermsofthesummation ofthecontributions
from the individual edges of a n-sided prism. The final Plouff’s formula for the
verticalcomponent ofgravitycaused byan-sidedrightpolygonal prismis
"
(cid:18) (cid:19)
Xn
z d z d
g D(cid:0)G(cid:26)s s A.z (cid:0)z /Cz arctan 2 1 (cid:0)arctan 2 2
m p 2 1 2
PR PR
iD1 12 22 #
(cid:18) (cid:19) (cid:18) (cid:19)
z d z d R Cd R Cd
(cid:0)z arctan 1 1 (cid:0)arctan 1 2 (cid:0)P ln 22 2 11 1 (16)
1 PR PR R Cd R Cd
11 21 12 1 21 2
whereP istheperpendicular distancefromthecomputationpointtoaverticalside
ofthepolygonal prism
x y (cid:0)x y
P D 1 2 2 1 (17)
1s
p
1s D .1x/2 C.1y/2 Djd (cid:0)d j (18)
1 2
348 XIONGLIANDMICHELCHOUTEAU
and
x x Cy y
A Darccos 1 2 1 2 (19)
r r
1 2
The symbol s D 1 or (cid:0)1 depending on if the center of mass of the prism is
m
below orabove the computation point. Thesymbol s D 1or (cid:0)1depending on if
p
P ispositiveornegative.ForthespecialcaseofP D0,thegravityvalueandarea
subtendedbythecorrespondingedgeoftheprism,asviewedfromthecomputation
point,arezero.ComputertimecanbereducedbydeterminingthesumofanglesA
outsidetheindicatedsummation.Thesumis2(cid:25) forthecomputationpointslocated
overtheinteriorofthepolygon, zeroforexteriorpoints,(cid:25) overanedge,andequal
totheinteriorangleifthecomputationpointislocatedovertheintersection oftwo
edgesofthepolygon.
3.1. TEST USING THE STANDARD CUBIC MODEL
Plouff(1976)alreadyprovedthatgravityvaluesoutsidethebodyobtainedbyusing
Equation(16),foranyx(cid:0)y orientation ofarightrectangular prism,agreeexactly
with those obtained by using Equation (8). By using the standard cubic model,
our calculation results show that the twoFORTRANsingle precision codes based
on Equation (16) and Equation (8) result in a rms deviation 0.000088 (cid:22)Gal. The
maximum absolute difference is 0.000410 (cid:22)Gal at (20, 16, (cid:0)18) at which the
calculated gravity is31.580(cid:22)Gal.Orsay, therelative difference is13permillion.
The maximum relative difference 50 per million is located at (18, 20, (cid:0)2) where
the gravity has a value of 5.283 (cid:22)Gal. Again, these differences are caused by the
floatingpointprecision limitation.
3.2. APPROXIMATION TO A VERTICAL CYLINDER
Furthermore, a vertical cylinder is used to check Plouff’s formula. The cylinder
has a height of 20 m, a radius of 10 m, and a density of 1 g/cm3. The gravity
response of the cylinder can be analytically estimated. A right regular polygonal
prismwith16sidesisinscribedinthecylinder. Theprismhasthesamedensityas,
hence, a smaller mass than the cylinder. Its gravity anomaly calculated according
to Plouff’s formula is shown in Figure 4. The maximum absolute errors due to
the model approximation occur along the vertical axis of the cylinder and are less
than3.1(cid:22)Gal,about1%ofthemaximumgravityanomaly(317.18(cid:22)Gal)occuring
at the center of the top or bottom surface of the prism (Figure 5). Larger relative
errors occur, for example, nearly 2% at the vertical location of 20 m and (cid:0)20
m along the axis. A better gravity approximation can be expected by improving
the model approximation, especially by increasing the number of sides of the
polygonal prism.
Plouff’s formula is valid for any observation positions, outside, inside, on the
surfaceof,andatthecorners ofthepolygonal prism.