A new ellipsoidal gravimetric, satellite altimetry and astronomic boundary value problem, a case study: The geoid of Iran

Journal of Geodynamics 39 (2005) 545–568

A new ellipsoidal gravimetric, satellite altimetry and astronomic boundary value problem, a case study: The geoid of Iran Abdolreza Safari a , Alireza A. Ardalan a,∗ , Erik W. Grafarend a,b a

Department of Surveying and Geomatics Engineering, University of Tehran, P.O. Box 11365-4563, Tehran, Iran b Department of Geodesy and GeoInformatics, Stuttgart University Geschwister-Scholl-Str. 24 D, D-70174 Stuttgart, Germany Accepted 18 April 2005

Abstract A new gravimetric, satellite altimetry, astronomical ellipsoidal boundary value problem for geoid computations has been developed and successfully tested. This boundary value problem has been constructed for gravity observables of the type (i) gravity potential, (ii) gravity intensity (i.e. modulus of gravity acceleration), (iii) astronomical longitude, (iv) astronomical latitude and (v) satellite altimetry observations. The ellipsoidal coordinates of the observation points have been considered as known quantities in the set-up of the problem in the light of availability of GPS coordinates. The developed boundary value problem is ellipsoidal by nature and as such takes advantage of high precision GPS observations in the set-up. The algorithmic steps of the solution of the boundary value problem are as follows: - Application of the ellipsoidal harmonic expansion complete up to degree and order 360 and of the ellipsoidal centrifugal field for the removal of the effect of global gravity and the isostasy field from the gravity intensity and the astronomical observations at the surface of the Earth. - Application of the ellipsoidal Newton integral on the multi-cylindrical equal-area map projection surface for the removal from the gravity intensity and the astronomical observations at the surface of the Earth the effect of the residual masses at the radius of up to 55 km from the computational point. - Application of the ellipsoidal harmonic expansion complete up to degree and order 360 and ellipsoidal centrifugal field for the removal from the geoidal undulations derived from satellite altimetry the effect of the global gravity and isostasy on the geoidal undulations.

∗

Corresponding author. Tel.: +49 711 121 3389; fax: +49 711 121 3285. E-mail address: [email protected] (A.A. Ardalan).

0264-3707/$ – see front matter © 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.jog.2005.04.009

546

A. Safari et al. / Journal of Geodynamics 39 (2005) 545–568

- Application of the ellipsoidal Newton integral on the multi-cylindrical equal-area map projection surface for the removal from the geoidal undulations derived from satellite altimetry the effect of the water masses outside the reference ellipsoid within a radius of 55 km around the computational point. - Least squares solution of the observation equations of the incremental quantities derived from aforementioned steps in order to obtain the incremental gravity potential at the surface of the reference ellipsoid. - The removed effects at the application points are restored on the surface of reference ellipsoid. - Application of the ellipsoidal Bruns’ formula for converting the potential values on the surface of the reference ellipsoid into the geoidal heights with respect to the reference ellipsoid. - Computation of the geoid of Iran has successfully tested this new methodology. © 2005 Elsevier Ltd. All rights reserved. Keywords: Geoid computations; Ellipsoidal approximation; Ellipsoidal boundary value problem; Ellipsoidal Bruns’ formula; Satellite altimetry; Astronomical observations

1. Introduction During the past decade the geodetic community has seen tremendous progress in the accuracy of the gravity observations on the earth’s surface and slightly above (aerogravimetry). Nowadays, sub-microgal accuracy for the land gravity observations is obtainable. Recent satellite gravimetry missions like CHAMP, GOCE, GRACE, have supplied the geodesists with huge amount of high quality global gravity data. The following list summarizes the contributions of the geodetic community to the studies related CHAMP, GOCE, GRACE satellite missions within the year 2002 (Albertella et al., 2002; Moreaux and Balmino, 2002; Iorio, 2002; Iorio et al., 2002; Jakowski et al., 2002; Knudsen and Andersen, 2002; Konig et al., 2002; Leuliette et al., 2002; Wickert et al., 2002; Oberndorfer et al., 2002; Petrovskaya et al., 2002; Reigber et al., 2002; Rodell and Famiglietti, 2002; Schroter et al., 2002; Luhr et al., 2002). The availability of the GPS coordinates for surface points in relation to the boundary value problems for geoid computations resulted in that the outer boundary (i.e. the surface of the Earth) can be regarded as a known and fixed. In addition to the above-mentioned progresses in the quality of input data, geodesists and geophysicists put big effort into refining the theory of geoid computations. The goal is to be able to cope with the high accuracy of modern observations, and to make the theory compatible in its accuracy with the input data. Some of the recent contributions are Ardalan and Grafarend (2001, 2004), Blitzkow (1999), Denker and Torge (1993), Featherstone and Kirby (2000), Featherstone et al. (1998), Forsberg (1993), Hwang (1997), Kearsley et al. (1998), Novak et al. (2001a, 2001b), Vajda and Vanicek (1999), Vanicek et al. (2001), Tziavos and Andritsanos (1999) and Tziavos (1996). Ardalan and Grafarend (Ardalan and Grafarend, 2004), and (Ardalan, 1999) presented a new methodology for geoid computations based on the gravity observables of the type gravity potential and the gravity intensity, via the solution of the so called ellipsoidal fixed-free two-boundary value problem. In this paper we have further developed the fixed-free two-boundary value problem to include boundary data of the type geoidal heights from the satellite altimetry, astronomical longitude and astronomical latitude. This new boundary value problem is now operational according to the following flowchart:

A. Safari et al. / Journal of Geodynamics 39 (2005) 545–568

547

Chart 1. Geoid computations procedure based on fixed-free two-boundary value problem.

Chart 1 Implementation of the new types of observables introduced to the boundary value problem has following advantages: (i) Adds to the degree of freedom of the solution. (ii) The use of the satellite altimetry data in the set-up of the problem helps the problem of geoid computations in the coastal areas. As an example, we refer to Japan where dense and accurate gravity data are available on islands and low accurate gravity data are available at sea where strong geophysical signatures are present. Satellite altimetry data combined with gravity data over land areas can properly solve this problem, since the satellite altimetry data can replace the gravity data at sea. One should keep in mind that the satellite altimetry observations having centimeter accuracy at the open ocean, which is equivalent to microgal level of accuracy in the gravity space, and at the coastal areas, where we have degraded decimeter accuracy, we still have the information equivalent to tens of microgal. This is of course more accurate than sea gravity observations. See, for example, Anzenhofer et al. (1999) for more details on the accuracy of satellite altimetry observations. (iii) The use of the vast amount of astronomical data in the set-up of the problem. Combination of the GPS and the total station observations (GPS/LPS observations) is a new data source for measuring the astronomical longitude and the astronomical latitude (see Grafarend and Awange, 2000; Awange, 1999 for more details).

548

A. Safari et al. / Journal of Geodynamics 39 (2005) 545–568

The idea of using satellite altimetry data in the problem of geoid computations has already been proposed by Svensson (1983), Sacerdote and Sanso (1983, 1987), Grafarend and Keller (1995), Jinghai and Xiaoping (1997) and Lehmann (1999). Here, we differ from aforementioned contributions, by the type of the observation equation that we introduce via application of the inverse Bruns’ formula. This observation equation is then added to the observation equations for other type of surface gravity data and solved via a combined least squares adjustment. The application of the astronomical observations in the problem of geoid computations goes back to Helmert (1880) and his “astro-geodetic geoid computation algorithm” for the computation of the geoidal height differences. Here, in contrast to the classical implementation of the astronomical observations we are using these data as boundary values along with other type of boundary observables, e.g. the gravity observations. In this way, various types of boundary values will be used in a unified manner for the geoid computations. Section 1 covers the details of the set-up of our gravimetry, the satellite altimetry and the astronomic boundary value problem. The fixed-free two-boundary value problem for the gravity observables of the type modulus of the gravity intensity is a nonlinear problem. The linearization of the problem will be discussed in Section 2. For the linearization, we have taken advantage of prior information of the global gravity field by means of the ellipsoidal harmonic expansion complete up to degree and order 360 as well as the ellipsoidal centrifugal field. Also, the effect of the local terrain masses up to distance of 55 km around the computational points. In Section 3, we will introduce the Abel–Poisson integral, which plays the central role in the solution of our boundary value problem. Section 4 has been devoted to the introduction of the ellipsoidal Bruns’ formula of the Somigliana–Pizzetti type, which we have used in its direct and inverse forms. Finally, in Section 5, the proposed boundary value problem has been numerically tested via computations of the geoid of Iran.

2. Formulation of the over-determined, nonlinear, fixed-free two-boundary value problem for the geoid computations We define the fixed-free two-boundary value problem based on the field equation of the type Laplace–Poisson partial differential equation with the boundary data of the type modulus of the gravity intensity (i.e. the gravity acceleration γ(x)) the gravity potential w(x), and the astronomical observations of the longitude and the latitude (Λ(x), Φ(x)) on the surface of the Earth. Furthermore, the satellite altimetry data over the sea areas were included. Thanks to Global Positioning System (GPS), the surface of the earth Mh2 can be assumed a fixed boundary. The potential value of the geoid w0 , comprises the boundary data on the free boundary, i.e. the geoid Mg2 . The surface of the sea Ms2 also serves as a fixed boundary, since its geometry is point-wise known from the satellite altimetry measurements. A brief definition of the fixed-free two-boundary value problem is given in Table 2.1. The fixed-free two-boundary value problem as defined above is an over-determined, nonlinear, oblique boundary value problem. It is an over-determined boundary value problem because there is more than minimum one boundary data (see Ardalan, 1999, p. 13). The fixed-free two-boundary value problem is nonlinear since boundary operators for the modulus of the gravity intensity, the astronomical longitude, and the astronomical latitude are nonlinear operators. The problem is an oblique boundary value problem because the direction of the gravity vector is in general not perpendicular to the boundary, i.e. to the earth surface.

A. Safari et al. / Journal of Geodynamics 39 (2005) 545–568

549

Table 2.1 Nonlinear fixed-free two boundary value problem 1. 2. 3. 4. 5. 6. 7. 8. 9.

div grad w(x) = 2ω2 (outside the Earth’s masses) div grad w(x) = −4πGσ + 2ω2 (inside the surface of the Earth) E{||grad w||} = µλ (boundary value data of the type modulus of gravity) E{w} = µw (boundary value data of they type gravity potential) E{h} = µh (boundary value data of the type satellite altimetry) E{Φ} = µ (boundary value data of the type astronomical latitude) E{Λ} = µ (boundary value data of the type astronomical longitude) w(x) = w0 (equipotential value at the level
of thegeoid close to mean sea level) w(x) = 21 ω2 ||x − x|eω eω ||22 +

gm ||x||2

+ Ow

1 ||x||32

∀x ∈ R3 /D ∪ ∂G+ e ∀x ∈ D ∪ ∂G− e ∀x ∈ ∂Ge = Mh2 ∀x ∈ ∂Ge = Mh2 ∀x ∈ ∂Ge = Ms2 ∀x ∈ ∂Ge = Mh2 ∀x ∈ ∂Ge = Mg2 ∀x ∈ ∂Ge = Mg2

||x||2 → ∞ (regularity condition at infinity)

3. Linear fixed-free two-boundary value problem The first step towards the solution of the proposed fixed-free two-boundary value problem is the linearization of the problem. Here, we apply the Euler δ-perturbation theory (Ardalan, 1999). Indeed, we take advantage of prior information of the terrestrial gravity field by means of the normal/reference gravity field of the ellipsoidal harmonic expansion complete up to degree and order 360, the ellipsoidal centrifugal potential, and the effect of the local terrain in the vicinity of the computational point. The concept of Euler δ-increment for the observable of the type gravity intensity ||γ(x)||2 = γ, the mass density σ and the angular velocity ω is summarized in Table 3.1. In Table 3.2, we have also expanded the modulus of the gravity intensity into Taylor series up to the order of magnitude of 
|δ
2 δΓ 2 . Table 3.3 is devoted to the linear part of Taylor series expansion of the boundary operator of the modulus of the gravity intensity. Note that e |δ
 Table 3.3 is the directional derivative of δ
in the direction of the reference gravity vector
. In order to linearize the astronomical observables (Λ,Φ), we use the relation between the reference longitude and the latitude and the astronomical longitude and latitude as done by Grafarend et al. (2003). Table 3.4 shows the relation between the astronomical coordinates (Λ, Φ) and the reference coordinates (λ , φ ). λ is the reference longitude and φ is the reference latitude. Table 3.5 shows how incremental astronomical observables δΛ, δΦ, can be constructed. The reference latitude and the reference longitude can be computed as explained in Table 3.6. To be able to implement the satellite altimetry observations in a linearized boundary value problem, we will go through the various stages outlined in Table 3.7 as follows. In step (i), we subtract the sea surface Table 3.1 Euler δ-increments of the gravity intensity γ, the mass density σ, and the angular velocity ω γ = Γ + δΓ

(2)

σ = p + δp

(3)

ω2 = Ω2 + 2Ω|δΩ + δΩ2

(4)

Where γ is the norm of the gravity vector, Γ the reference gravity intensity, δΓ the incremental gravity intensity, p the reference mass density, δp the incremental mass density, Ω the reference angular velocity and δΩ is the incremental angular velocity.

A. Safari et al. / Journal of Geodynamics 39 (2005) 545–568

550

Table 3.2 Taylor series expansion of the nonlinear boundary operator γ(x) up to the order 5 in δΓ (i) Nonlinear boundary operator γ(x) = ||grad w(x)|| = ||grad W(x) + grad δW(x)|| = grad W(x) + grad δW(x)|grad W(x) + grad δW(x) 1/2 = ( grad δW(x) 2 + grad W(x) 2 + 2grad W(x)|grad δW(x))  1/2 2 1 = Γ (x) 1 + 2 
|δ
 + 2 δΓ 2 Γ Γ

(5)

Subject to Taylor expansion of (1 + x)1/2 where |x| < 1, i.e. 1 1 1 (1 + x)1/2 = 1 + x − x2 + x3 + O(x4 ) 2 8 16

(6)

We have (ii) Expansion of the nonlinear boundary operator γ(x) γ(x) = Γ (x) +

1 1 1 1 1 
|δ
 + δΓ 2 − 
|δ
2 − 
|δ
δΓ 2 − δΓ 4 + Oγ {δΓ 5 } Γ 2Γ 2Γ 3 2Γ 3 8Γ 3

(7)

(iii) Expansion of the nonlinear boundary operator δΓ (x)



⇒

δΓ (x) = γ(x) − Γ (x) =

1 1 1 1 1 
|δ
 + δΓ 2 − 
|δ
2 − 
|δ
δΓ 2 − δΓ 4 +Oγ {δΓ 5 } 3 3 Γ 2Γ 2Γ 2Γ 8Γ 3

(8)

Table 3.3 Linearized boundary operator δΓ According to Table 2.1, the linearized boundary operator δΓ can be constructed as follows: 1 δΓ (x) = γ(x) − Γ (x) = 
|δ
 + O(δΓ 2 ) Γ  
|δ
+ O(δΓ 2 ) = Γ = e |δ
 + O(δΓ 2 ) = e |grad δW  + O(δΓ 2 ) = ∇e δW + O(δΓ 2 )

(9)

Where ∇e is the directional derivative in the direction of e , the unit vector along the reference gravity vector. Table 3.4 Euler δ-increment for the astronomical coordinates (Λ, Φ) and the Jacobi ellipsoidal coordinates (λ, φ) (see Appendix A) Λ = λ + δΛ,

(10)

Φ = φ + δΦ,

(11)

λ = λ + δλ,

(12)

φ = φ + δφ.

(13)

A. Safari et al. / Journal of Geodynamics 39 (2005) 545–568

551

Table 3.5 The linearized boundary operators δΛ and δΦ δΛ = −

1 Γ cos φ

1 =− Γ cos φ 1 δΦ = − Γ =−

1 Γ

 

 

1 1 1 √ Dλ δW + sin φδλ √ Dφ δW − cos φδλ √ Dη δW gλλ gφφ gηη 1 1 1 √ Dλ + sin φδλ √ Dφ − cos φδλ √ Dη gλλ gφφ gηη

1 1 − sin φδλ √ Dλ δW + √ Dφ δW − gλλ gφφ 1 1 − sin φδλ √ Dλ + √ Dφ − gλλ gφφ





δφ +





δW

1 sin 2φ δφ + 4 sinh2 η

1 sin 2φ 4 sinh2 η



(14)



1 √ Dη δW gηη

1 √ Dη gηη





δW

(15)

{λ, φ, η} are the Jacobi ellipsoidal coordinates as defined in Appendix A, and {Dλ , Dφ , Dη } are the partial derivatives with respect to the Jacobi ellipsoidal coordinates {λ, φ, η}.

topography hSST from the mean sea level heights hMSL derived from satellite altimetry. The result is the geoid undulation N. The sea surface topography hSST is computed by the POCM 4B model (Rapp, 1998). In step (ii), we convert the geoid undulation N to its equivalent potential value dW(X) = w0 − W(X). In 2 step (iii), we compute the actual gravity potential w(X) at the surface of the ellipsoid Ea,b . In step (iv), we remove the reference gravitational field of the ellipsoidal harmonic expansion complete up to degree and order 360 and the centrifugal potential. In step (v), we reduce the residual effect of the local sea masses above the reference ellipsoid up to distance 55 km around the computational point. Finally, the linearized fixed-free two-boundary value problem is given in Table 3.8.

4. Abel–Poisson integral The incremental observables of the types (i) the gravity potential, (ii) the gravity intensity, (iii) the astronomical longitude and (iv) the astronomical latitude (v) the satellite altimetry can now be downward continued to the surface of the reference ellipsoid via the inverse solution to the Abel–Poisson integral, and its derivatives (see Table 4.1). However, we must first discretize the integral equations and build a system of linear equations of the type y + i = Ax, where y stands for the five incremental observations mentioned above, x is the vector of unknowns (i.e. the incremental gravity potential at the surface of reference ellipsoid), and i is the vector of the inconsistencies caused by the random error of the observations. The following systems of equations present our discretized observation equations:  y1 + i1 = A1 x        y2 + i2 = A2 x

y +i =A x

3 3 3     y + i = A 4 4x  4  

y5 + i5 = A5 x

(1)

552

Table 3.6 The reference gravity vector
and the reference longitude λΓ and latitude φΓ (i) The reference gravity vector

= grad

360 n 

unm

n=0 m=−n

Q∗n|m| (sinh η) Q∗n|m| (sinh η0 )

 enm (λ, φ)

 + grad

1 2 2 ω ε cosh2 η cos2 φ 2

 + grad(Uterrain )

1 ∂(V (φ, η)) 1 ∂(Uterrain (λ, φ, η)) 1 ∂(Uterrain (λ, φ, η)) 1 ∂(Uterrain (λ, φ, η)) eη + √ eλ + √ eφ + √ eη +√ gηη ∂η gλλ ∂λ gφφ ∂φ gηη ∂η

 =

1 ∂(U(λ, φ, η) + Uterrain (λ, φ, η)) √ ∂λ gλλ

 ×

eλ eφ eη

1 ∂(U(λ, φ, η) + V (φ, η) + Uterrain (λ, φ, η)) √ ∂φ gφφ

1 ∂(U(λ, φ, η) + V (φ, η) + Uterrain (λ, φ, η)) √ ∂η gηη



 (16)

where U(λ, φ, η) =

360  n 

unm

n=0 m=−n

Q∗n|m| (sinh η) Q∗n|m| (sinh η0 )

enm (λ, φ)

(17)

is the ellipsoidal harmonic expansion of the gravitational potential field of the Earth to degree and order 360, and V (φ, η) =

1 2 2 ω ε cosh2 η cos2 φ 2

(18)

is the centrifugal potential in terms of Jacobi ellipsoidal coordinates, and {eλ , eφ , eη } are base vectors of the local moving frame. (ii) Reference gravity vector
in terms of the reference longitude λ and the reference latitude φ .


= [ Γ cos φ cos λ

Γ cos φ sin λ

Γ sin φ ]

e1 e2 e3



where {e1 , e2 , e3 } are the base vectors of the global Cartesian frame.

(19)

A. Safari et al. / Journal of Geodynamics 39 (2005) 545–568

1 ∂(U(λ, φ, η)) 1 ∂(U(λ, φ, η)) 1 ∂(V (φ, η)) 1 ∂(U(λ, φ, η)) eλ + √ eφ + √ eη + √ eφ = √ gλλ ∂λ gφφ ∂φ gηη ∂η gφφ ∂φ

Table 3.6 (Continued ) (iii) The relation between the local moving frame {eλ , eφ , eη } with the global Cartesian frame {e1 , e2 , e3 }: eλ eφ eη





e1 = T e2 e3

 (20)

where the rotation matrix T is defined as follows:



− sin λ − cosh η 

cos λ sin φ  2 2 T =  sinh η + sin φ sinh η   cos λ cos φ sinh2 η + sin2 φ

cos λ

 

− cosh η

sinh η + sin φ sinh η 2



0 2

sinh2 η + sin2 φ

sin λ sin φ



sin λ cos φ



sinh η

sinh η + sin φ sinh η 2

2

sinh2 η + sin2 φ

   

cos φ  sin φ

(21)

(iv) Inserting (20) in (16) we get


=

1 ∂(U(λ, φ, η) + Uterrain (λ, φ, η)) √ gλλ ∂λ

 ×T

e1 e2 e3

1 ∂(U(λ, φ, η) + V (φ, η) + Uterrain (λ, φ, η)) √ gφφ ∂φ

1 ∂(U(λ, φ, η) + V (φ, η) + Uterrain (λ, φ, η)) √ gηη ∂η



 (22)

A. Safari et al. / Journal of Geodynamics 39 (2005) 545–568



(v) From equality of Eqs. (22) and (19) the reference longitude λΓ and the reference latitude φΓ can be obtained.

553

A. Safari et al. / Journal of Geodynamics 39 (2005) 545–568

554

Table 3.7 2 Conversion of the satellite altimetry to the incremental gravity potential δW(X)at the surface of the reference ellipsoid Ea,b (i) Combine the mean sea level (MSL) hMSL data from satellite altimetry and the sea surface topography (SST) hSST data in order to get the geoidal undulation N: N = hMSL − hSST

(23)

2 via (ii) Conversion of the geoidal undulation N to the equivalent potential dW(X) at the surface of the reference ellipsoid Ea,b the nonlinear ellipsoid Bruns’ formula: 1/2

N = dW(X)

ε2 cosh η(cos2 hη − cos2 φ) gm

(24)

2 (iii) Computation of the actual gravity potential at the surface of the ellipsoid Ea,b :

w(X) = w0 + dW(X)

(25)

(iv) Removal of the reference potential field W(X) (the expansion of the gravitational potential to the degree and order 360 and the centrifugal potential) δW(X) = w(X) − W(X)

(26)

2 (v) Reduction of the potential WS (X) generated by the seawater masses above the reference ellipsoid Ea,b

δW(X) = w(X) − W(X) − WS (X)

(27)

1/2

N = (W(X) + WS (X) + δW(X) − w0 )

ε2 cosh η(cosh2 η − cos2 φ) gm

1/2

N = (W(X) + WS (X) − w0 )

ε2 cosh η(cosh2 η − cos2 φ) gm

(28)

1/2

+ δW(X)

N = hR (X) + δh (X)

ε2 cosh η(cosh2 η − cos2 φ) gm

(29)

(30)

Where 1/2

hR (X) = (W(X) + WS (X) − w0 )

ε2 cosh η(cosh2 η − cos2 φ) gm

(31)

and δh is the incremental geoid height 1/2

δh(X) = δW(X)

ε2 cosh η(cosh2 η − cos2 φ) gm

(32)

div grad δW(x) = 0, ∀x ∈ R3 /D ∪ ∂Gg div grad δW(x) = −4πδp(x), ∀x ∈ D ∪ ∂Gg

 field diff. Eq.

δΓ (x) = ∇e δW(x) 1/2 ε2 cosh η(cosh2 η − cos2 φ) δW(X) δh(X) = gm δΛ(x) = −

1 Γ cos φ





1 1 1 √ Dλ + sin φδλ √ Dφ − cos φδλ √ Dη gλλ gφφ gηη

1 1 1 δΦ(x) = − − sin φδλ √ Dλ + √ Dφ − Γ gλλ gφφ δW(x) = δW(x) w0 = W(x) + δW(x) δW(x) =

gm ||x||2

+ Oδω




1 ||x||32



for||x||2 → ∞

    2   ∀X ∈ Ea,b        ∀x ∈ ∂Ge = M 2  ∀x ∈ ∂Ge = Mh2  





1 δφ + sin 2φ 4 sinh2 η



regularity at infinity

δW(x)



1 √ Dη gηη

h



    δW(x) ∀x ∈ ∂Ge = Mh2      2 ∀x ∈ ∂Ge = Mh    2   ∀x ∈ ∂Gi = Mg

Boundary Value

A. Safari et al. / Journal of Geodynamics 39 (2005) 545–568

Table 3.8 Linearized fixed-free two-boundary value problem

555

556

Table 4.1 The ellipsoidal Abel–Poisson integral for the incremental observable δw, δΓ , (δΛ, δΦ)

δw(λ, φ, u) = w(φ ) = √

 2 Ea,b

w(φ ) S

dS




a b2 +ε2 sin2 φ¯

2 S = area(Ea,b ) = 4πa

!

×

∞  n 



Qn|m| (i sinh η) Qn|m| (i sinh η0 )

× enm (λ , φ )enm (λ, φ)]δw(λ , φ ) =

1 S

 2 Ea,b

n=0 m=−n 1 2

"

+ 1 2

1 b2 4 aε

+

ln

1 b2 4 aε

a+ε a−ε

ln

a+ε a−ε

dS w(φ )K(λ, φ, η, λ , φ , η0 )δw(λ , φ )



(ii) The ellipsoidal Abel–Poisson integral for the of the incremental gravitational intensityδΓ  modulus



,φ ,η ) Γφ 1 Γλ 1

∂K(λ,φ,η,λ ,φ ,η0 ) 0 dS ω(φ ) δW(λ , φ ) + √g1φφ ||
|| dS ω(φ ) ∂K(λ,φ,η,λ δW(λ , φ ) δΓ (x) = γ(x) − Γ (x) = e |δΓ  = √g1λλ ||
|| S E2 ∂λ S E2 ∂φ +

Γη 1 1 √ gηη ||
|| S



2 Ea,b

a,b

a,b

,φ ,η0 ) dS ω(φ ) ∂K(λ,φ,η,λ δW(λ , φ ) ∂η

(iii) The ellipsoidal # Abel–Poisson integral for the astronomical observables(δΛ, δΦ) δΛ = − Γ

1 cos φ

1 √ gλλ

− cos φδλ √g1ηη





δΦ =

− Γ1






1 4 sinh2 η

) ∂K(λ,φ,η,λ ,φ ,η0 ) dS ω(φ δW(λ , φ ) + sin φδλ √g1φφ S ∂λ



2 Ea,b

− sin φδλ √g1λλ

− δφ +



2 Ea,b

$

) ∂K(λ,φ,η,λ ,φ ,η0 ) dS ω(φ δW(λ , φ ) S ∂η

%

2 Ea,b



sin 2φ

) ∂K(λ,φ,η,λ ,φ ,η ) 0 dS ω(φ δW(λ , φ ) S ∂λ

% 1 √ gηη

2 Ea,b



% +

1 √ gφφ

2 Ea,b



) ∂K(λ,φ,η,λ ,φ ,η ) 0 dS ω(φ δW(λ , φ ) ∂η S

2 Ea,b





) ∂K(λ,φ,η,λ ,φ ,η0 ) dS ω(φ δW(λ , φ ) S ∂φ





) ∂K(λ,φ,η,λ ,φ ,η0 ) dS ω(φ δW(λ , φ ) S ∂φ

A. Safari et al. / Journal of Geodynamics 39 (2005) 545–568

(i) The ellipsoidal Abel–Poisson integral for the gravity potential δW

A. Safari et al. / Journal of Geodynamics 39 (2005) 545–568

557

y1 stands for the incremental gravity potential (from the precise leveling), y2 stands for the incremental gravity intensity (from gravimetry), y3 stands for the incremental astronomical longitude (from geodetic astronomy or GPS/LPS observations), y4 stands for the incremental astronomical latitude (from the geodetic astronomy or GPS/LPS observations), and y5 stands for the incremental geoid height δh(x) (from the satellite altimetry). Determination of the incremental gravitational potential values at the surface of reference ellipsoid via the inverse Abel–Poisson integral is an improperly posed problem and as such must be regularized. Here, we have applied the Phillips–Tikhonov regularization method (Phillips, 1962; Tikhonov, 1963). For the details of the regularization process for the Abel–Poisson integral, we refer to Ardalan (1999). Having derived the incremental gravity potential values at the surface of the reference ellipsoid we can restore the effect of the removed (i) centrifugal potential, (ii) the reference gravitational potential of the ellipsoidal harmonic expansion up to degree and order 360 and (iii) the gravitational potential of topographical masses up to a distance of 55 km from the computational point at the surface of the reference ellipsoid. 5. The ellipsoidal Bruns’ formula Table 5.1 summarizes the common reference gravitational fields, associated with a surface of a sphere or an ellipsoid of revolution. The geoid can be defined as deviation with respect to any of the reference equipotential surfaces presented in Table 5.1. However, only the fourth model in Table 5.1, i.e. the reference field of Somigliana–Pizzetti (see Pizzetti, 1894; Somigliana, 1930), produces a reference equipotential surface which is compatible with the resolution No. 7 of IUGG (see the Geodesist’s Handbook, Moritz, 2000). World Geodetic Datum 2000 (Grafarend and Ardalan, 1999) is a reference equipotential surface of the Somigliana–Pizzetti type based on the current best estimates of the fundamental geodetic parameters {GM, W0 , J2 , Ω}. The computation of the directional derivatives of the actual gravity potential along the normal direction to the surface of the reference equipotential surface should be done based on the reference field. The ellipsoidal Bruns’ formulas associated with these reference fields are listed in Table 5.2. Ardalan and Grafarend (2001) have done a numerical verification for the accuracy of different variants of the ellipsoidal Bruns’ formula. The conclusion is that the Bruns’ formula computed based on the reference fields (i) and (iii) does not provide centimeter accuracy unless the nonlinear terms up to degree 3 in terms of the Table 5.1 Common types of the reference equipotential gravitational fields and their corresponding equipotential surface Reference potential fields

Ref. equipotential surfaces

Spherical gm r GM r

+

1 2 2 r ω 3

+

Ellipsoidal & ' gm arccot uε e gm arccot ε

&u' ε




1 R5 ω 2 √ 3 5 π2

−

r 2√ ω2 3 5

√

5 (3 2

Reference sphere sin φ − 1) 2

& u2 ' 3 +1 arccot( uε )−3 uε 1 2 2 & ε2 ' & b ' b (3 sin2 φ − 1) + 21 ω2 (u2 + ε2 ) cos2 φ + 6w a u2 3 2 +1 arccot ε −3 ε ε

Bjerhammar sphere Ref. ellipsoid Somigliana Pizzetti ref. Ellipsoid, e.g. WGD2000

A. Safari et al. / Journal of Geodynamics 39 (2005) 545–568

558

Table 5.2 Different variants of the ellipsoidal Bruns’ formulas, up to the order of accuracy of O(δW(x0 )2 ) Reference Field

Bruns’ formula

Spherical

& gm ' R2  &

gm r

h = dW(R)/

Bjerhammar reference field

h = −dW(X0 )

Ellipsoidal & ' gm arccot uε ε

a h gm

Somigliana–Pizzetti reference field



h=−

− gm r2

+

2 rω2 3

'

 +

 −

R5 ω 2 √

5r 4

−

2 2rω √ 3 5

√ 5 (3 sin2 2

−1 φ − 1)

b2 + ε2 sin2 φδW(X0 ) + O(δW(X0 )2 )




gm

− 2 ω2 a

β2 +ε2 sin2 φ a 6(b2 +ε2 )b arccot bε −3bε+ε2 =3ε (3b2 +ε2 )arccot bε −3bε

( ) ( )



δW(X0 ) (3 sin2 φ+1)+ω2 b cos φ

disturbing gravity potential are incorporated. The Bruns’ formula associated with the Somigliana–Pizzetti reference field offers millimeter accuracy even for the linear part. In addition, the geoidal undulations computed using the Somigliana–Pizzetti field have the advantage of fulfilling the Gauss criterion of zero value for the global mean of the geoidal undulations over the reference equipotential surface (Heiskanen and Moritz, 1967). Therefore, Ardalan and Grafarend (2001) have recommended the ellipsoidal Bruns’ formula based on the Somigliana–Pizzetti reference field for the modern sub-millimeter geoid computations.

6. Case study: the geoid of Iran In this section, the results of our case study, i.e. the construction of the geoid of Iran based on the methodology explained in the previous section will be presented. As the observables, we will use the gravity intensity values for the test area. Figs. 6.1 and 6.2 show the location of the gravity stations and the astronomical stations. Fig. 6.3, which has been produced using the 1 km × 1 km DTM file that we used for the terrain reduction, gives a picture to the topography of the test region. The maximum and

Fig. 6.1. Coverage map of 8483 gravity intensity stations in Iran (from BGI data bank).

A. Safari et al. / Journal of Geodynamics 39 (2005) 545–568

559

Fig. 6.2. Coverage map of 10 first order astronomical stations in the geographical region of Iran.

minimum elevations in this region, according to the DTM file, are 4010 m and the −64 m, respectively. In Figs. 6.4–6.6, we have plotted the surface gravity observations of (i) the type modulus of the gravity intensity, (ii) the satellite altimetry data and (iii) the sea surface topography in the Persian Gulf and in the Oman Sea. Next, we have removed the effect of the reference gravitational field of the ellipsoidal harmonic expansion complete up to degree and order 360 and the ellipsoidal centrifugal potential as well as the residual effect of the local terrain gravitational field up to a distance 55 km from the computational

Fig. 6.3. Color plot of the topography of the test region based on 1 km × 1 km DTM file.

560

A. Safari et al. / Journal of Geodynamics 39 (2005) 545–568

Fig. 6.4. Color plot presentation of the variation of the modulus of the gravity intensity in Iran based on the 8483 gravity data.

points. Fig. 6.7 shows the remaining signals for the gravity intensity after the remove step. The remaining signal for the astronomical observations after the remove step is given in Table 6.1. For the computation of terrain correction, we have first transformed the ellipsoidal mass elements using a multi-cylindrical equal-area projection and then applied the closed form solution of the Newton integral in Cartesian

Fig. 6.5. Color plot presentation of mean sea level (MSL) variation at Persian Gulf and Oman Sea based on satellite altimetry data of Topex-Posseidon.

A. Safari et al. / Journal of Geodynamics 39 (2005) 545–568

561

Fig. 6.6. Color plot presentation of sea surface topography (SST) at Persian Gulf and Oman Sea based on the SST data file.

Fig. 6.7. Incremental gravitational potential after removal of global reference field of the ellipsoidal harmonic expansion to degree/order 360/360 plus centrifugal acceleration and the residual effect of local terrain masses in a radius of 55 km around the computational point (at 8483 gravity stations).

562

A. Safari et al. / Journal of Geodynamics 39 (2005) 545–568

Table 6.1 The astronomical longitude and the astronomical altitude of the 10 first order Laplace stations in Iran after the reduction for the effects of the reference values and the terrain correction Incremental longitude, δ∆

Incremental longitude, δΦ

0◦ 0 03.98

0◦ 0 20.26

0◦ 0 20.81

0◦ 0 02.11

0◦ 0 04.42

0◦ 0 12.32

0◦ 0 15.70

0◦ 0 14.64

0◦ 0 11.12

0◦ 10 28.14

0◦ 10 07.25

0◦ 10 16.24

0◦ 0 23.31

0◦ 10 05.74

0◦ 09 35.83

0◦ 09 01.92

0◦ 08 18.19

0◦ 07 30.75

coordinates. For details on this, we refer to Ardalan and Safari (2004). Subsequently, the incremental gravitational intensity is downward continued to the gravitational potential values at the surface of the reference ellipsoid of WGD2000. Fig. 6.8 shows the variation of the trace of the mean square error (MSE) matrix of the estimated parameters versus the regularization parameter α. From this figure, the optimum value of the regularization factor has been selected as 7.22 × 10−6 . Fig. 6.9 shows the result

Fig. 6.8. Variations of trace of mean square error (MSL) matrix of estimated parameters versus the regularization parameter α.

A. Safari et al. / Journal of Geodynamics 39 (2005) 545–568

563

Fig. 6.9. Downward continued residual gravity plus the residual potential values derived from satellite altimetry data.

Fig. 6.10. Gravity potential at the surface of the reference ellipsoid of WGD2000 after restoring the reference gravity field of ellipsoidal harmonic expansion to degree/order 360/360 plus the centrifugal acceleration and the residual effect of local terrain masses in a radius of 55 km around the computational point at the surface of reference ellipsoid.

564

A. Safari et al. / Journal of Geodynamics 39 (2005) 545–568

Fig. 6.11. Computed geoid of Iran, with respect to World Geodetic Datum 2000 based upon collocation of linearized observational functionals of the type GPS, gravity intensity and satellite altimetry data.

of the downward continuation after Phillips–Tikhonov regularization. Now, if we restore the reference ellipsoidal gravitational potential field up to the degree and order 360 and the centrifugal field as well as the residual gravitational potential associated with the local terrain masses we will arrive to the gravity potential values at the surface of the reference ellipsoid. Fig. 6.10 is the result of the restoring process for the gravity potential values obtained from the downward continuation of the gravity intensity. Next, the gravity potential on the surface of the reference ellipsoid of WGD2000 is converted to the geoidal undulation with respect to the reference ellipsoid WGD2000 via Bruns’ formula (see Table 5.2). Fig. 6.11 represents the computed geoid of Iran. Finally, the computed geoid has been tested on the GPS/leveling stations. Table 6.2 summarizes the statistics of the difference between the computed geoid and the geoid derived from GPS/leveling at 51 stations along the first-order leveling network of Iran. It is important to

Table 6.2 The statistics of the difference between the computed geoid heights and 51 GPS/leveling stations along the first-order leveling network of Iran Maximum (m) Minimum (m) Mean (m) STD (m)

2.066 −2.148 0.161 1.068

A. Safari et al. / Journal of Geodynamics 39 (2005) 545–568

565

add that this difference includes the error in the zero point of the datum, and the effect of the orthometric correction, which is not applied to the precise leveling network of Iran. The order of magnitude for the difference is consistent with the sum of the aforementioned errors.

7. Conclusions We can summarize the highlights of our approach as follows: (1) Formulation of the fixed-free two-boundary value problem with respect to the reference ellipsoid. (2) Application of a high degree and order expansion of the global geopotential model, which provides us with the global gravity and isostasy at the remove step. (3) Implementation of the observables of the type the gravity intensity (from gravimetry), the astronomical longitude, and the astronomical latitude (from the classical geodesy or from combination of GPS and total station observations), and the mean sea level data (from satellite altimetry) in a combined model for geoid computation. (4) Ellipsoidal approximation has been consistently preserved throughout the computations. (5) Application of the ellipsoidal Bruns’ formula of the Somigliana–Pizzetti type to convert the incremental gravity potential at the surface of the reference ellipsoid to the geoidal undulations. (6) Computation of the terrain correction using a close form solution of the Newton’s integral in Cartesian coordinates at a multi-cylindrical equal-area map projection surface.

Acknowledgements The authors would like to thank University of Tehran for the financial support provided to this work by Grant number 8151007/1/01. The constructive review of Professor G. Strykowski and Professor M. Vermeer is gratefully acknowledged. Their comments and corrections helped us to improve the initial version of the paper.

Appendix A. Conversion of Cartesian coordinates {x, y, z} into ellipsoidal coordinates {λ, φ, η} Forward transformation of the ellipsoidal coordinates {λ, φ, η} into Cartesian coordinates {x, y, z}: x = ε cosh η cos φ cos λ y = ε cosh η cos φ sin λ z = ε sinh η sin φ √ where ε = a2 − b2 is the linear eccentricity of the reference ellipsoid with the major semi axis a and the minor semi axis b.

A. Safari et al. / Journal of Geodynamics 39 (2005) 545–568

566

The inverse transformation of the Cartesian coordinates {x, y, z} into the ellipsoidal coordinates {λ, φ, η}  x  arctan , for x > 0 and y ≥ 0    y    y    arctan + π, for x < 0 and y = 0   x   y λ = arctan + 2π, for x > 0 and y < 0  x      π, for x = 0 and y > 0    2    π  3 , for x = 0 and y < 0 2 

φ = {sgn z} arcsin

1 2 [ε − (x2 − +y2 + z2 ) + 2ε2



1 2 η = arccosh [x + y2 + z2 + ε2 + 2ε2

(

(x2

+

1/2

(

y2

(x2 + y2 + z2 − ε2 )2 + 4ε2 z2 ]

+

z2

−

ε2 )2

1/2

−

4ε2 (x2

+

y2 ]

Valid for λ ∈ {λ ∈ R|0 < λ < 2π}   π π φ ∈ φ ∈ R| − < φ < + 2 2 h ∈ {η ∈ R|η > 0}

References Albertella, A., Migliaccio, F., Sanso, F., 2002. GOCE: the Earth gravity field by space gradiometry. Celestial Mech. Dyn. Astronomy 83, 1–15. Anzenhofer, M., Shum, C.K., Rentsh, M., 1999. Coastal Altimetry and Applications. Department of Civil and Environmental Engineering and Geodetic Science of the Ohio State University, Columbus Ohio, Report No. 464. Ardalan, A.A., 1999. High resolution regional geoid computation in the world geodetic datum 2000. Based upon collocation of linearized observational functionals of the type GPS, gravity potential and gravity intensity. Ph.D. thesis. Stuttgart University. Ardalan, A.A., Grafarend, E.W., 2001. Ellipsoidal geoidal undulations (ellipsoidal Bruns’ formula): case studies. J. Geod. 75, 544–552. Ardalan, A.A., Grafarend, E.W., 2004. High-resolution geoid computation without applying Stokes’s formula; case study: highresolution geoid of Iran. J. Geod. 78, 138–156. Ardalan, A.A., Safari, A., 2004. Terrain correction on the multi-cylindrical equal-area map projection of the surface of the reference ellipsoid. J. Geod. 78, 114–123. Awange, J., 1999. www.uni-stuttgart.de/gi/research/schriftenreihe/quo vadis/pdf/awange.pdf. Blitzkow, D., 1999. Toward a 10’ resolution geoid for South America: a comparison study. Phys. Chem. Earth Part A Solid Earth Geod. 24, 33–39. Denker, H., Torge, W., 1993. Present state and future developments of the European geoid. Surv. Geophys. 14, 433–447. Featherstone, W.E., Evans, J.D., Olliver, J.G., 1998. A Meissl-modified Vanicek and Kleusberg kernel to reduce the truncation error in gravimetric geoid computations. J. Geod. 72, 154–160.

A. Safari et al. / Journal of Geodynamics 39 (2005) 545–568

567

Featherstone, W.E., Kirby, J.F., 2000. The reduction of aliasing in gravity anomalies and geoid heights using digital terrain data. Geophys. J. Int. 141, 204–212. Forsberg, R., 1993. Modeling the fine structure of the geoid—methods, data requirements and some results. Surv. Rev. 14, 403–418. Grafarend, E.W., Ardalan, A.A., 2000. World Geodetic Datum. J. Geod. 73, 611–623. Grafarend, E.W., Awange, J., 2000. Determination of vertical deflections by GPS/LPS measurements. Z. Vermessungswesen 125, 279–288. Grafarend, E.W., Finn, G., Ardalan, A.A., 2003. Ellipsoidal vertical deflections and ellipsoidal gravity disturbances, case studies. J. Geod. (Under review). Grafarend, E.W., Keller, W., 1995. Setup of observational functionals in gravity space as well as in geometry space. Manuscr. Geod. 20, 301–305. Helmert, F.R., 1880. Die mathematischen und physikalischen Theorien der h¨oheren Geode¨asie, vol. I. Minerva GMBH reprint 1962. Hwang, C., 1997. Analysis of some systematic errors affecting altimeter-derived sea surface gradient with application to geoid determination over Taiwan. J. Geod. 71, 113–130. Heiskanen, W., Moritz, H., 1967. Physical Geodesy. W.H. Freeman and Co., San Francisco. Iorio, L., 2002. A critical approach to the concept of a polar, low-altitude LARES satellite. Classical Quantum Gravity 19, 175–183. Iorio, L., Ciufolini, I., Pavlis, E.C., 2002. Measuring the relativistic perigee advance with satellite laser ranging. Classical Quantum Gravity 19, L4301–L4309. Jakowski, N., Wehrenpfennig, A., Heise, S., Reigber, C., Luhr, H., Grunwaldt, L., Meehan, T.K., 2002. GPS radio occultation measurements of the ionosphere from CHAMP: early results. Geophys. Res. Lett. 29 (Art. no. 1457). Jinghai, Y., Xiaoping, W., 1997. The solution of mixed boundary value problems with the reference ellipsoid as boundary. J. Geod. 71, 454–460. Kearsley, A.H.W., Forsberg, R., Olesen, A., Bastos, L., Hehl, K., Meyer, U., Gidskehaug, A., 1998. Airborne gravimetry used in precise geoid computations by ring integration. J. Geod. 72, 600–605. Knudsen, P., Andersen, O., 2002. Correcting GRACE gravity fields for ocean tide effects. Geophys. Res. Lett. 29 (Art. no. 1178). Konig, R., Zhu, S., Reigber, C., Neumayer, K.H., Meixner, H., Galas, R., Baustert, G., Schwintzer, P., 2002. Champ rapid orbit determination for GPS atmospheric limb sounding. New Trends Space Geod. 30, 289–293. Lehmann, R., 1999. Boundary-value problems in the complex world of geodetic measurments. J. Geod. 73, 491–500. Leuliette, E.W., Nerem, R.S., Russell, G.L., 2002. Detecting time variations in gravity associated with climate change. J. Geophys. Res. Solid Earth 107 (Art. no. 2112). Luhr, H., Maus, S., Rother, M., Cooke, D., 2002. First in-situ observation of night-time F region currents with the CHAMP satellite. Geophys. Res. Lett. 29 (Art. no. 1489). Moreaux, G., Balmino, G., 2002. Impact of some land hydrological phenomena on GOCE mission. Geophys. Res. Lett. 29 (Art. no. 1205). Moritz, H., 1980. Geodetic reference system. J. Geod. 74, 128–133. Novak, P., Vanicek, P., Martinec, Z., Veronneau, M., 2001a. Effects of the spherical terrain on gravity and the geoid. J. Geod. 75, 491–504. Novak, P., Vanicek, P., Veronneau, M., Holmes, S., Featherstone, W., 2001b. On the accuracy of modified Stokes’s integration in high-frequency gravimetric geoid determination. J. Geod. 74, 644–654. Oberndorfer, H., Muller, J., Rummel, R., Sneeuw, N., 2002. A simulation tool for the new gravity field satellite missions. New Trends Space Geod. 30, 227–232. Petrovskaya, M.S., Vershkov, A.N., Zielinski, J.B., 2002. Recovering the Earth’s potential spectral characteristics from GOCE mission. New Trends Space Geod. 30, 221–226. Phillips, D.L., 1962. A technique for the numerical solution of certain integral equations of the first kind. J. Ass. Comput. Mach. 9, 84–96. Pizzetti, P., 1894. Geodesia—Sulla es pressione della gravita alla superficie del geoide, supposto ellissoidico. Atti Reale Accademia dei Lincei 3, 166–172. Rapp, R.H., 1998. The Development of a degree 360 expansion of the dynamic ocean topography of the POCM 4B global circulation model, NASA/CR-1998-206877. Greenbelt Maryland 20771.

568

A. Safari et al. / Journal of Geodynamics 39 (2005) 545–568

Reigber, C., Luhr, H., Schwintzer, P., 2002. CHAMP mission status. New Trends Space Geod. 30, 129–134. Rodell, M., Famiglietti, J.S., 2002. The potential for satellite-based monitoring of groundwater storage changes using GRACE: the High Plains aquifer, Central US. J. Hydrol. 263, 245–256. Sacerdote, F., Sanso, F., 1983. A contribution to the analysis of the altimetry–gravimetriy problem. Bull. Geod. 57, 257–272. Sacerdote, F., Sanso, F., 1987. Further remarks on the altimetry–gravimetriy problems. Bull. Geod. 61, 65–82. Schroter, J., Losch, M., Sloyan, B., 2002. Impact of the gravity field and steady-state ocean circulation explorer (GOCE) mission on ocean circulation estimates. 2. Volume and heat fluxes across hydrographic sections of unequally spaced stations. J. Geophys. Res. Oceans 107 (C2) (Art. no. 3012). Somigliana, C., 1930. Geofisica—Sul campo gravitazionale esterno del geoide ellissoidico. Atti della Reale Academia Nazionale dei Lincei Rendiconti 6, 237–243. Svensson, S.L., 1983. Solution of the altimetry–gravimetry problem. Bull. Geod. 57, 332–353. Tikhonov, A.N., 1963. The regularization of incorrectly posed problems. Soviet Math. Doklady 4, 1624–1627. Tziavos, I.N., Andritsanos, V.D., 1999. Recent geoid computations for the Hellenic area. Phys. Chem. Earth Part A Solid Earth Geod. 24, 91–96. Tziavos, I.N., 1996. Comparisons of spectral techniques for geoid computations over large regions. J. Geod. 70, 357–373. Vajda, P., Vanicek, P., 1999. Truncated geoid and gravity inversion for one point-mass anomaly. J. Geod. 73, 58–66. Vanicek, P., Novak, P., Martinec, Z., 2001. Geoid, topography, and the Bouguer plate or shell. J. Geod. 75, 210–215. Wickert, J., Beyerle, G., Hajj, G.A., Schwieger, V., Reigber, C., 2002. GPS radio occultation with CHAMP: atmospheric profiling utilizing the space-based single difference technique. Geophys. Res. Lett. 29 (8) (Art. no. 1187).