引力位函数定义及球谐计算

Scientific Technical Report 09/02 Deutsches GeoForschungsZentrum GFZ Contents 1 Introduction 2 2 Definitions 2 2.1 The Potential and the Geoid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.2 The Height Anomaly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.3 The Gravity Disturbance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.4 The Gravity Anomaly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.4.1 The Classical Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.4.2 The Modern Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.4.3 The Topography-Reduced Gravity Anomaly . . . . . . . . . . . . . . . . . . . . . . 7 3 Approximation and Calculation 8 3.1 The Geoid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 3.2 The Height Anomaly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.3 The Difference: Geoid - Height Anomaly . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.4 The Gravity Disturbance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.5 The Gravity Anomaly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.5.1 The Classical Gravity Anomaly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.5.2 The Modern Gravity Anomaly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3.5.3 The Topography-Reduced Gravity Anomaly . . . . . . . . . . . . . . . . . . . . . . 17 4 Calculation from Spherical Harmonics 17 4.1 Spherical Harmonics and the Gravity Field . . . . . . . . . . . . . . . . . . . . . . . . . . 17 4.2 The Geoid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 4.3 The Height Anomaly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 4.4 The Gravity Disturbance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 4.5 The Gravity Anomaly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 5 Practical calculations using the model EIGEN-GL04C 24 5.1 Geoid and Height Anomaly . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 5.2 Gravity Disturbance and Gravity Anomaly . . . . . . . . . . . . . . . . . . . . . . . . . . 24 1 Scientific Technical Report 09/02 Deutsches GeoForschungsZentrum GFZ 1 Introduction The intention of this article is to present the definitions of different functionals of the Earth’s gravity field and possibilities for their approximative calculation from a mathematical representation of the outer potential. In history this topic has usually been treated in connection with the boundary value problems of geodesy, i.e. starting from measurements at the Earth’s surface and their use to derive a mathematical representation of the geopotential. Nowadays global gravity field models, mainly derived from satellite measurements, become more and more detailed and accurate and, additionally, the global topography can be determined by modern satellite methods independently from the gravity field. On the one hand the accuracy of these gravity field models has to be evaluated and on the other hand they should be combined with classical (e.g. gravity anomalies) or recent (e.g. GPS-levelling-derived or altimetry-derived geoid heights) data. Furthermore, an important task of geodesy is to make the gravity field functionals available to other geosciences. For all these purposes it is necessary to calculate the corresponding functionals as accurately as possible or, at least, with a well-defined accuracy from a given global gravity field model and, if required, with simultaneous consideration of the topography model. We will start from the potential, formulate the definition of some functionals and derive the formulas for the calculation. In doing so we assume that the Earth’s gravity potential is known outside the masses, the normal potential outside the ellipsoid and that mathematical representations are available for both. Here we neglect time variations and deal with the stationary part of the potential only. Approximate calculation formulas with different accuracies are formulated and specified for the case that the mathematical representation of the potential is in terms of spherical harmonics. The accuracies of the formulas are demonstrated by practical calculations using the gravity field model EIGEN-GL04C (Fo¨rste et al., 2006). More or less, what is compiled here is well-known in physical geodesy but distributed over a lot of articles and books which are not cited here. In the first instance this text is targeted at non-geodesists and it should be “stand-alone readable”. Textbooks for further study of physical geodesy are (Heiskanen & Moritz, 1967; Pick et al., 1973; Van´ıcˇek & Krakiwsky, 1982; Torge, 1991; Moritz, 1989; Hofmann-Wellenhof & Moritz, 2005). 2 Definitions 2.1 The Potential and the Geoid As it is well-known, according to Newton’s law of gravitation, the potential Wa of an attractive body with mass density ρ is the integral (written in cartesian coordinates x, y, z) Wa(x, y, z) = G ∫∫ v ∫ ρ(x′, y′, z′)√ (x− x′)2 + (y − y′)2 + (z − z′)2 dx ′dy′dz′ (1) over the volume v of the body, where G is the Newtonian gravitational constant, and dv = dx′dy′dz′ is the element of volume. For √ (x− x′)2 + (y − y′)2 + (z − z′)2 → ∞ the potential Wa behaves like the potential of a point mass located at the bodies centre of mass with the total mass of the body. It can be shown that Wa satisfies Poisson’s equation ∇2Wa = −4piGρ (2) where ∇ is the Nabla operator and ∇2 is called the Laplace operator (e.g. Bronshtein et al., 2004). Outside the masses the density ρ is zero and Wa satisfies Laplace’s equation ∇2Wa = 0 (3) 2 Scientific Technical Report 09/02 Deutsches GeoForschungsZentrum GFZ thus Wa is a harmonic function in empty space (e.g. Blakely, 1995). On the rotating Earth, additionally to the attracting force, also the centrifugal force is acting which can be described by its (non-harmonic) centrifugal potential Φ(x, y, z) = 1 2 ω2d 2z (4) where ω is the angular velocity of the Earth and dz = √ x2 + y2 is the distance to the rotational (z-) axis. Hence, the potential W associated with the rotating Earth (e.g. in an Earth-fixed rotating coordinate system) is the sum of the attraction potential Wa and the centrifugal potential Φ W =Wa +Φ (5) The associated force vector ~g acting on a unit mass, the gravity vector, is the gradient of the potential ~g = ∇W (6) and the magnitude g = |∇W | (7) is called gravity. Potentials can be described (and intuitively visualised) by its equipotential surfaces. From the theory of harmonic functions it is known, that the knowledge of one equipotential surface is sufficient to define the whole harmonic function outside this surface. For the Earth one equipotential surface is of particular importance: the geoid. Among all equipotential surfaces, the geoid is the one which coincides with the undisturbed sea surface (i.e. sea in static equilib- rium) and its fictitious continuation below the continents as sketched in Fig. 1 (e.g. Van´ıcˇek & Christou, 1994, Van´ıcˇek & Krakiwsky, 1982 or Hofmann-Wellenhof & Moritz, 2005). Being an equipotential sur- ellipsoid gravity v ector topography geoid N U = Uo W = Uoht H Figure 1: The ellipsoid, the geoid and the topography face, the geoid is a surface to which the force of gravity is everywhere perpendicular (but not equal in magnitude!). To define the geoid surface in space, simply the correct value W0 of the potential has to be chosen: W (x, y, z) =W0 = constant (8) As usual we split the potential W into the normal potential U and the disturbing potential T W (x, y, z) = U(x, y, z) + T (x, y, z) (9) 3 Scientific Technical Report 09/02 Deutsches GeoForschungsZentrum GFZ and define “shape” and “strengths” of the normal potential as follows: (a) The equipotential surfaces (U(x, y, z) = constant) of the normal potential should have the shapes of ellipsoids of revolution and (b) the equipotential surface for which holds U(x, y, z) = W0 (see eq. 8) should approximate the geoid, i.e. the undisturbed sea surface, as good as possible (i.e. in a least squares fit sense). It is advantageous to define ellipsoidal coordinates (h, λ, φ) with respect to this level ellipsoid U(h = 0) = U0 = W0, where h is the height above ellipsoid (measured along the ellipsoidal normal), λ is the ellipsoidal longitude and φ the ellipsoidal latitude. Thus eq. (9) writes (note that the normal potential U does not depend on λ): W (h, λ, φ) = U(h, φ) + T (h, λ, φ) (10) and the geoid, in ellipsoidal coordinates, is the equipotential surface for which holds W ( h = N(λ, φ), λ, φ ) = U ( (h = 0), φ ) = U0 (11) where N(λ, φ) is the usual representation of the geoid as heights N with respect to the ellipsoid (U = U0) as a function of the coordinates λ and φ. Thus N are the undulations of the geoidal surface with respect to the ellipsoid. This geometrical ellipsoid together with the normal ellipsoidal potential is called Geodetic Reference System (e.g. NIMA, 2000 or Moritz, 1980). Now, with the ellipsoid and the geoid, we have two reference surfaces with respect to which the height of a point can be given. We will denote the height of the Earth’s surface, i.e. the height of the topography, with respect to the ellipsoid by ht, and with respect to the geoid by H, hence it is (see fig. 1): ht(λ, φ) = N(λ, φ) +H(λ, φ) (12) Here H is assumed to be measured along the ellipsoidal normal and not along the real plumb line, hence it is not exactly the orthometric height. A discussion of this problem can be found in (Jekeli, 2000). Like the potential W (eq. 5) the normal potential also consists of an attractive part Ua and the centrifugal potential Φ U = Ua +Φ (13) and obviously, the disturbing potential T (h, λ, φ) =Wa(h, λ, φ)− Ua(h, φ) (14) does not contain the centrifugal potential and is harmonic outside the masses. The gradient of the normal potential ~γ = ∇U (15) is called normal gravity vector and the magnitude γ = |∇U | (16) is the normal gravity. For functions generated by mass distributions like Wa from eq. (1), which are harmonic outside the masses, there are harmonic or analytic continuations W ca which are equal to Wa outside the masses and are (unlike Wa) also harmonic inside the generating masses. But the domain where Wa is harmonic, i.e. satisfies Laplace’s equation (eq. 3), can not be extended completely into its generating masses because there must be singularities somewhere to generate the potential, at least, as is well known, at one point at the centre if the mass distribution is spherically symmetric. How these singularities look like (point-, line-, or surface-singularities) and how they are distributed depends on the structure of the function Wa outside the masses, i.e. (due to eq. 1) on the density distribution ρ of the masses. Generally one can say that these singularities are deeper, e.g. closer to the centre of mass, if the potential Wa is “smoother”. For further study of this topic see (Zidarov, 1990) or (Moritz, 1989). 4 Scientific Technical Report 09/02 Deutsches GeoForschungsZentrum GFZ Due to the fact that the height H = ht −N of the topography with respect to the geoid is small compared to the mean radius of the Earth and that in practise the spatial resolution (i.e. the roughness) of the approximative model for the potential Wa will be limited (e.g. finite number of coefficients or finite number of sampling points), we expect that the singularities of the downward continuation of Wa lie deeper than the geoid and assume that W ca exists without singularities down to the geoid so that we can define (ht is the ellipsoidal height of the Earth’s surface, see eq. 12): W ca(h, λ, φ) = Wa(h, λ, φ) for h ≥ ht ∇2W ca = 0 for h ≥ min(N,ht) W c(h, λ, φ) = W ca(h, λ, φ) + Φ(h, φ) (17) However, this can not be guaranteed and has to be verified, at least numerically, in practical applications. From its definition the normal potential Ua is harmonic outside the normal ellipsoid and it is known that a harmonic downward continuation U ca exists down to a singular disk in the centre of the flattened rotational ellipsoid (e.g. Zidarov, 1990). Thus, downward continuation of the normal potential is no problem and we can define U ca(h, φ) = Ua(h, φ) for h ≥ 0 ∇2U ca = 0 for h ≥ min(N, 0, ht) U c(h, φ) = U ca(h, φ) + Φ(h, φ) (18) and hence T c(h, λ, φ) = W ca(h, λ, φ)− U ca(h, φ) ∇2T c = 0 for h ≥ min(N,ht) (19) 2.2 The Height Anomaly The height anomaly ζ(λ, φ), the well known approximation of the geoid undulation according to Molo- densky’s theory, can be defined by the distance from the Earth’s surface to the point where the normal potential U has the same value as the geopotential W at the Earth’s surface (Molodensky et al., 1962; Hofmann-Wellenhof & Moritz, 2005; Moritz, 1989): W (ht, λ, φ) = U(ht − ζ, λ, φ) (20) where ht is the ellipsoidal height of the Earth’s surface (eq. 12). An illustration of the geometrical situation is given in fig. (2). The surface with the height ζ = ζ(λ, φ) with respect to the ellipsoid is often called quasigeoid (not shown in fig. 2) and the surface ht − ζ is called telluroid. It should be emphasised, that the quasigeoid has no physical meaning but is an approximation of the geoid as we will see. In areas where ht = N (or H = 0) i.e. over sea, the quasigeoid coincides with the geoid as can be seen easily from the definition in eq. (20) if we use eq. (12): W (N +H,λ, φ) = U(N +H − ζ, λ, φ) (21) set H = 0 and get W (N,λ, φ) = U(N − ζ, φ) (22) and use eq. (11), the definition of the geoid to write U(0, φ) = U(N − ζ, φ) (23) from which follows N = ζ for H = 0 (24) 5 Scientific Technical Report 09/02 Deutsches GeoForschungsZentrum GFZ N ζ ζ tW(h ) = W = Uo U = Uo N ζU(h − )t geoid ellipsoid topography telluroid th = N N =ζ if equipotential surfaces surfac esequipo tential th t h Figure 2: The ellipsoid, the geoid, and the height anomaly ζ In the history of geodesy the great importance of the height anomaly was that it can be calculated from gravity measurements carried out at the Earth’s surface without knowledge of the potential inside the masses, i.e. without any hypothesis about the mass densities. The definition of eq. (20) is not restricted to heights h = ht on the Earth’s surface, thus a generalised height anomaly ζg = ζg(h, λ, φ) for arbitrary heights h can be defined by: W (h, λ, φ) = U(h− ζg, φ) (25) 2.3 The Gravity Disturbance The gradient of the disturbing potential T is called the gravity disturbance vector and is usually denoted by ~δg: ~δg(h, λ, φ) = ∇T (h, λ, φ) = ∇W (h, λ, φ)−∇U(h, φ) (26) The gravity disturbance δg is not the magnitude of the gravity disturbance vector (as one could guess) but defined as the difference of the magnitudes (Hofmann-Wellenhof & Moritz, 2005): δg(h, λ, φ) = ∣∣∇W (h, λ, φ)∣∣− ∣∣∇U(h, φ)∣∣ (27) In principle, herewith δg is defined for any height h if the potentials W and U are defined there. Ad- ditionally, with the downward continuations W ca and U ca (eqs. 17 and 18), we can define a “harmonic downward continued” gravity disturbance δgc(h, λ, φ) = ∣∣∇W c(h, λ, φ)∣∣− ∣∣∇U c(h, φ)∣∣ (28) With the notations from eqs. (7) and (16) we can write the gravity disturbance in its common form: δg(h, λ, φ) = g(h, λ, φ)− γ(h, φ) (29) The reason for this definition is the practical measurement process, where the gravimeter measures only |∇W |, the magnitude of the gravity, and not the direction of the plumb line. 6 Scientific Technical Report 09/02 Deutsches GeoForschungsZentrum GFZ 2.4 The Gravity Anomaly The term gravity anomaly is used with numerous different meanings in geodesy and geophysics and, moreover, there are different practical realisations (cf. Hackney & Featherstone, 2003). Here we will confine ourselves to the classical free air gravity anomaly, to the gravity anomaly according to Molodensky’s theory and to the topography-reduced gravity anomaly. 2.4.1 The Classical Definition The classical (historical) definition in geodesy is the following (cf. Hofmann-Wellenhof & Moritz, 2005): The gravity anomaly ∆gcl (subscript “cl” stands for “classical”) is the magnitude of the downward continued gravity |∇W c| (eq. 17) onto the geoid minus the normal gravity |∇U | on the ellipsoid at the same ellipsoidal longitude λ and latitude φ: ∆gcl(λ, φ) = ∣∣∇W c(N,λ, φ)∣∣− ∣∣∇U(0, φ)∣∣ (30) The origin of this definition is the (historical) geodetic practise where the altitude of the gravity mea- surement was known only with respect to the geoid from levelling but not with respect to the ellipsoid. The geoid height N was unknown and should be determined just by these measurements. The classical formulation of this problem is the Stokes’ integral (e.g. Hofmann-Wellenhof & Moritz, 2005; Martinec, 1998). For this purpose the measured gravity |∇W (ht, λ, φ)| has to be reduced somehow down onto the geoid and the exact way to do so is the harmonic downward continuation of the attraction potential Wa (eq. 17). This is the reason for the definition of the classical gravity anomaly in eq. (30). In practise the so-called “free air reduction” has been or is used to get |∇W c(N,λ, φ)| approximately. Thus the classical gravity anomaly depends on longitude and latitude only and is not a function in space. 2.4.2 The Modern Definition The generalised gravity anomaly ∆g according to Molodensky’s theory (Molodensky et al., 1962; Hofmann-Wellenhof & Moritz, 2005; Moritz, 1989) is the magnitude of the gravity at a given point (h, λ, φ) minus the normal gravity at the same ellipsoidal longitude λ and latitude φ but at the ellipsoidal height h− ζg, where ζg is the generalised height anomaly from definition (25): ∆g(h, λ, φ) = ∣∣∇W (h, λ, φ)∣∣− ∣∣∇U(h− ζg, φ)∣∣, for h ≥ ht (31) or in its common form: ∆g(h, λ, φ) = g(h, λ, φ)− γ(h− ζg, φ) (32) Here the height h is assumed on or outside the Earth’s surface, i.e. h ≥ ht, hence with this definition the gravity anomaly is a function in the space outside the masses. The advantage of this definition is that the measured gravity |∇W | at the Earth’s surface can be used without downward continuation or any reduction. If geodesists nowadays speak about gravity anomalies, they usually have in mind this definition with h = ht, i.e. on the Earth’s surface. 2.4.3 The Topography-Reduced Gravity Anomaly For many purposes a functional of the gravitational potential is needed which is the difference between the real gravity and the gravity of the reference potential and which, additionally, does not contain the effect of the topographical masses above the geoid. The well-known “Bouguer anomaly” or “Refined Bouguer anomaly” (e.g. Hofmann-Wellenhof & Moritz, 2005) are commonly used in this connection. However, they are defined by reduction formulas and not as functionals of the potential. The problems arising when using the concepts of the Bouguer plate or the Bouguer shell are discussed in (Van´ıcˇek et al., 2001) and (Van´ıcˇek et al., 2004). 7 Scientific Technical Report 09/02 Deutsches GeoForschungsZentrum GFZ Thus, let us define the gravity potential of the topography Vt, i.e. the potential induced by all masses lying above the geoid. Analogously to eq. (27), we can now define a gravity disturbance δgtr which does not contain the gravity effect of the topography: δgtr(h, λ, φ) = ∣∣∇[W (h, λ, φ)− Vt(h, λ, φ)]∣∣− ∣∣∇U(h, φ)∣∣ (33) and, analogously to eq. (31), a topography-reduced gravity anomaly ∆gtr: ∆gtr(h, λ, φ) = ∣∣∇[W (h, λ, φ)− Vt(h, λ, φ)]∣∣− ∣∣∇U(h− ζ, φ)∣∣ (34) where, consequently, W − Vt is the gravity potential of the Earth without the masses