Flechtner F.M., Gruber Th., G?ntner A., Mandea M., Rothacher M., Sch?ne T., Wickert J. (Eds.) System Earth via Geodetic-Geophysical Space Techniques

Подождите немного. Документ загружается.

250 O. Baur et al.

The eigenvalues  (0 ≤  ≤ 1), denoted as concentration factors, are a measure

for the spectral concentration of the eigenvectors g in B. They are dependent on the

kernel D, representing the area of base function concentration. In the special case of

an axisymmetric double polar cap its elements become

lml





π−θ



sin θ dθ ⇔ D

= 

.(5)

Due to its block-diagonal structure, the eigenvalue problem in Eq. (5) can be

solved for each order m separately. The eigenvectors related to large concen-

tration factors, i.e.,  ≈ 1, allow to set up a system of (normalized) band-

limited base functions

( cos θ). They result from the product between the

eigenvectors g

and the (normalized) associated Legendre functions

( cos θ)

subject to

( cos θ) =



l=m

( cos θ), m ≤ j ≤ L .(6)

The Slepian functions are orthogonal within the spherical belt B, orthonormal on

the sphere respectively. For constant order m they are constructed by a linear

combination of Legendre functions of the same order. The order m is identi-

cal for both Legendre and Slepian functions. This does not hold for the ﬁrst

index. In the following we refer j to as the Slepian index to distinguish from the

degree l.

4.1 Solving the Eigenvalue Problem

The problem Albertella et al. (1999) were faced with is that the EVP according to

Eq. (5) can not be solved uniquely. The parity of most of the concentration factors,

cf. the left panel in Fig. 6, only allows the non-structured computation of the eigen-

vectors, hence Slepian functions. Grünbaum et al. (1982) gave the key impulse to

overcome this computational problem. For the axisymmetric double polar cap con-

ﬁguration they found a commuting EVP to the original one. As outlined in Fig. 5,

the commuting matrices D

and T

, subject to D

= T

, generate the same

Fig. 5 Slepian functions

computation via commuting

operators

Spectral Approaches to Solving the Polar Gap Problem 251

10 20 30 40 50

−3000

−2000

−1000

Grünbaum eigenvalues

10 20 30 40 50

0.5

Concentration factors

Fig. 6 Spectrum of concentration factors (left panel) and Grünbaum eigenvalues ( right panel);

exemplary for order m = 0 and spectral resolution L = 50

Co−latitude [deg]

Function index

45 90 135 180

−5

Co−latitude [deg]

Function index

0 45 90 135 180

−5

Fig. 7 Spatial oscillation of Legendre functions (left panel) and Slepian functions (right panel);

exemplary for order m = 0 and spectral resolution L = 50

eigenvectors but different eigenvalues. Hence, in terms of Slepian functions compu-

tation both the original and the commuting EVP are equivalent: D

= 

↔

= X

The spectrum of the so-called Grünbaum eigenvalues X

, cf. the right panel

in Fig. 6, allows the unique estimation of the eigenvectors and hence Slepian func-

tions. Moreover, as T

turns out to be tridiagonal with its elements following simple

analytical expressions, the computation of the eigenvectors is highly efﬁcient.

Figure 7 displays Slepian functions opposed to Legendre functions, exemplary

for the order m =0 and a spectral resolution of L =50. As already mentioned before,

Legendre functions exhibit no spatial concentration. Independent of the function

index l, they oscillate over the whole sphere. Slepian functions, on the other hand,

clearly demonstrate their band-limited nature. Low function indices j represent base

functions that are exclusively concentrated in the polar caps, followed by a small

transmission zone, with signal content in both the polar caps and the spherical belt.

From the seventh base function on there is negligible signal in the polar caps. Hence

(for m = 0 and L = 50), these Slepian functions constitute a system of orthogonal

band-limited base functions that are optimally concentrated on the spherical belt of

GOCE data coverage.

The Slepian parameterization deﬁnitively solves the PGP in the “Slepian world”.

In this context, it may be interpreted as a method for regional gravity ﬁeld determi-

nation, as the polar caps are totally ignored within the analysis process. We have to

252 O. Baur et al.

point out that the Slepian approach has no beneﬁt compared to the Legendre param-

eterization of the geopotential for GOCE data analysis on a global scale (Baur and

Sneeuw, 2007).

5 Conclusions

The GOCE satellite’s ground tracks leave out a double polar cap with a radius of

more than six degrees. As a consequence, GOCE data analysis based on spherical

harmonics is confronted with numerical instabilities, as the spatial data gaps map

into non-resolvable low-order constituents. Spectral approaches are efﬁcient and

powerful to solve the PGP. We stabilize the normal equations system by tailored

regularization. In terms of a combination approach, we incorporate EGM08 to pro-

vide a global GOCE-based gravity ﬁeld model. In this context we recommend the

mixed estimator with additional information as stochastic linear restrictions.

Besides stabilization with external a priori information, the parameterization of

the Earth gravity ﬁeld by means of band-limited base functions, or Slepian func-

tions, ultimately solves to PGP. The transformation of base functions is done by

solving an algebraic eigenvalue problem. Commuting matrices allow the unique

estimation of the new base. The philosophy of the Slepian approach is not to adapt

the observation scenario to the modeling but instead, to adapt the parameterization

of the geopotential to GOCE observations. Hence it can be interpreted as method for

regional gravity ﬁeld determination, meeting the ESA requirement for the recovery

of the Earth’s gravity ﬁeld from GOCE data only.

Acknowledgments This is publication no. GEOTECH-1243 of the programme

GEOTECHNOLOGIEN, Grant 03F0421C. We gratefully acknowledge the support of

the BMBF (Bundesministerium für Bildung und Forschung) and t he DFG (Deutsche

ForschungsGemeinschaft). Furthermore, we kindly acknowledge helpful comments by W.-D.

Schuh and an anonymous reviewer of the manuscript.

References

Albertella A, Sansò F, Sneeuw N (1999) Band-limited functions on a bounded spherical domain:

The Slepian problem on the sphere. J. Geod. 73, 436–447.

Baur O, Sneeuw N (2007) Slepian approach revisited: New studies to overcome the polar gap.

Proceedings of GOCE User Workshop, ESA SP-627.

Baur O, Sneeuw N, Grafarend EW (2008) Methodology and use of tensor invariants for satellite

gravity gradiometry. J. Geod., doi: 10.1007/s00190-007-0178-5.

Cai J (2004) Statistical inference of the eigenspace components of a symmetric random deforma-

tion tensor. DGK Series C 577, Munich.

Ditmar P, Klees R, Kostenko F (2003) Fast and accurate computation of spherical harmonic

coefﬁcients from satellite gravity gradiometry data. J. Geod., doi: 10.1007/s00190-002-0298-x.

ESA (1999) The four candidate Earth explorer core missions – gravity ﬁeld and steady-state ocean

circulation mission. ESA Report SP-1233(1), Granada.

Forsberg R, Kenyon S (2004) Gravity and geoid in the arctic region – The Northern polar gap now

ﬁlled. Proceedings of GOCE User Workshop, ESA-ESRIN.

Spectral Approaches to Solving the Polar Gap Problem 253

van Gelderen M, Koop R (1997) The use of degree variances in satellite gradiometry. J. Geod. 71,

337–343.

Grünbaum FA, Longhi L, Perlstadt M (1982) Differential operators commuting with ﬁnite

convolution integral operators: some non-abelian examples. SIAM J. Appl. Math. 42, 941–955.

Kusche J (2003) A Monte-Carlo technique for weight estimation in satellite geodesy. J. Geod., doi:

10.1007/s00190-002-0302-5.

Kusche J, Klees R (2002) Regularization of gravity ﬁeld estimation from satellite gravity gradients.

J. Geod. 76, 359–368.

Metzler B, Pail R (2005) GOCE data processing: The spherical cap regularization approach. Stud.

Geophys. Geod. 49, 441–462.

Miranian L (2004) Slepian functions on the sphere, generalized Gaussian quadrature rule. Inv.

Prob. 20, 877–892.

Pail R, Plank G, Schuh W-D (2001) Spatially restricted data distributions on the sphere: The

method of orthonormalized functions and applications. J. Geod. 75, 44–56.

Phillips DL (1962) A technique for the numerical solution of certain integral equations of the ﬁrst

kind. J. ACM 9, 54–97.

Rao CR, Toutenburg H (1999) Linear Models, Springer, Berlin, Heidelberg.

Reigber C (1989) Gravity ﬁeld recovery from satellite tracking data, theory of satellite geodesy

and gravity ﬁeld determination. In: Sansò F, Rummel R (eds.), Lecture Notes in Earth Science,

Vol. 25, Springer, Berlin, Heidelberg, pp. 197–234.

Rummel R, Sansò F, van Gelderen M, Brovelli M, Koop R, Miggliaccio F, Schrama E,

Scerdote F (1993) Spherical harmonic analysis of satellite gradiometry. Netherlands Geodetic

Commission, New Series 39.

Scheinert M (2005) The Antarctic geoid project: Status report and next activities. In: Jekeli C,

Bastos L, Fernandes J (eds.), Gravity, Geoid and Space Missions, Vol. 129, Springer, Berlin,

Heidelberg, pp. 137–142.

Simons FJ, Dahlen FA, Wieczorek MA (2006) Spatiospectral l ocalization on a sphere. SIAM Rev.

48, 504–536.

Slepian D (1978) Prolate spheroidal wave functions, Fourier analysis and uncertainty. V. The

discrete case. Bell. Syst. Tech. J. 57, 1371–1430.

Sneeuw N, van Gelderen M (1997) The polar gap. In: Sansò F, Rummel R (eds.), Geodetic

Boundary Value Problems in View of the One Centimeter Geoid. Lecture Notes in Earth

Science, Vol. 65, Springer, Berlin, pp. 559–568.

Tapley BD, Bettadpur S, Watkins M, Reigber C (2004) The gravity recovery and climate experi-

ment: Mission overview and early results. Geophys. Res. Lett., doi: 10.1029/2004GL019920.

Tikhonov AN (1963) Regularization of incorrectly posed problems. Sov. Mat. Dokl. 4, 1035–1038.

Regionally Reﬁned Gravity Field Models

from In-Situ Satellite Data

Annette Eicker, Torsten Mayer-Gürr, Karl-Heinz Ilk,

and Enrico Kurtenbach

1 Introduction

The third satellite in the sequence of dedicated gravity satellite missions after

CHAMP (CHAllenging Minisatellite Payload – Reigber et al. 1999) and GRACE

(Gravity Recovery and Climate Experiment – Tapley et al., 2004) will be GOCE

(Gravity ﬁeld and steady-state Ocean Circulation Explorer – ESA, 1999). The

gravity ﬁeld will thereby be measured with unprecedented accuracy by a gravity

gradiometer consisting of six three-axes accelerometers which determine in-orbit

gravity gradients in three spatial directions.

In the past decades several gravity ﬁeld analysis techniques tailored to this new

kind of observables have been proposed and tested based on various realistic sim-

ulation scenarios. One problem is the processing of the extremely large amount of

observations and in solving the huge unstable adjustment problem. Another problem

seems to be related to the representation of the gravity ﬁeld by an appropriate set

of basis functions and the associated gravity ﬁeld parameters. The usual way is to

model the gravity ﬁeld by spherical harmonics up to a certain degree. This degree is

limited by the signiﬁcance of the gravity ﬁeld signal in the observations and its ratio

to the observation noise. Because of the inhomogeneous gravity ﬁeld of the Earth

the signal content varies in the s pace domain. The gravity ﬁeld in regions with rough

gravity ﬁeld features should be recovered up to a higher degree than within regions

with smooth gravity ﬁeld features. This means that the degree of instability of the

gravity ﬁeld recovery process, which increases with increasing high frequencies in

the gravity ﬁeld signal, varies with the speciﬁc gravity ﬁeld characteristics in dif-

ferent regions. A gravity ﬁeld representation with basis functions of global support

must be regularized globally. However, a global regularization causes an overall ﬁl-

tering of the observations leading to a mean dampening of the global gravity ﬁeld

features. As a consequence, the high frequent gravity ﬁeld signal in the observations

is lost again in some geographical regions.

A. Eicker (B)

Institute of Geodesy and Geoinformation, University of Bonn, 53115 Bonn, Germany

e-mail: annette@geod.uni-bonn.de

255

F. Flechtner et al. (eds.), System Earth via Geodetic-Geophysical Space Techniques,

Advanced Technologies in Earth Sciences, DOI 10.1007/978-3-642-10228-8_20,



Springer-Verlag Berlin Heidelberg 2010

256 A. Eicker et al.

An alternative approach, presented here, is to determine a global gravity ﬁeld

solution with high long and medium wavelength accuracy and improve this global

solution in regions with characteristic gravity ﬁeld features by an adapted regional

recovery procedure. The global solution is parameterized by spherical harmonic

coefﬁcients up to a moderate degree and the regional solutions are represented by

space localizing basis functions, e.g., by spherical splines. This procedure provides

several advantages. The regional approach allows the exploitation of the individ-

ual signal content in the observations and a tailored regularization for regions

with different gravity ﬁeld characteristics. The advantage compared to a uniform

global regularization is that the regularization factor is selected for each region indi-

vidually. By an individually adapted regularization it is possible to extract more

information out of the given data than with a global gravity ﬁeld determination

process. Regions with a smooth gravity ﬁeld signal, for example, can be regularized

stronger without dampening the signal. In addition, the resolution of the gravity ﬁeld

determination can be chosen for each region individually according to the spectral

behavior of the signal in the speciﬁc region. Furthermore, the regional approach

has the advantage of dealing with regions with different data coverages more eas-

ily. If no data at all is available (e.g. the polar gap) the regional reﬁnement can

be skipped. For regions with sparse data coverage a coarser parameterization can be

selected. Another aspect, especially relevant for the GOCE mission with its potential

to recover the high resolution gravity ﬁeld, is the fact that regional solutions contain

less unknown parameters and therefore the computation procedure is simpliﬁed.

This enables a signiﬁcant reduction of computation costs.

If desired, several regional solutions with global coverage can be merged by

means of quadrature methods to obtain a global solution, in principle, up to an arbi-

trary degree. Simulation results are presented to demonstrate this approach. Due to

the regionally adapted strategy this method provides better results than calculating

a spherical harmonic solution by recovering the potential coefﬁcients directly (see

also Eicker et al., 2004 and Eicker, 2008).

2 Mathematical Model

The gradiometer as primary observation device in case of the GOCE mission for

the determination of the high frequency part of the gravity ﬁeld, measures the grav-

ity gradient consisting of the second derivatives of the gravitational potential. To

use these measurements in the gravity ﬁeld determination process, a mathematical

model has to be established that links the particular observations to the gravity ﬁeld

parameters. The unknown gravity ﬁeld parameters are then determined by a least-

squares adjustment procedure. Therefore, the gravity ﬁeld has to be parameterized in

terms of basis functions. In the approach presented here, a global reference solution,

parameterized by spherical harmonics is reﬁned by regional solutions modeled by

spherical splines as space localizing basis functions. The design of these basis func-

tions will be described in Sect. 2.1. When solving t he system of normal equations,

special emphasis is put on the ill-posedness of the downward continuation process

Regionally Reﬁned Gravity Field Models from In-Situ Satellite Data 257

that is inevitable when dealing with satellite data and aiming at the determination of

the gravity ﬁeld on the Earth’s surface. The use of space localizing basis functions

offers the opportunity to regionally adapt the necessary regularization process as

will be explained in Sect. 2.2.

The observations are given in the reference frame of the gradiometer (GRF),

oriented along the three axes of the gradiometer instrument with origin in the

nominal intersection of the three one-axis gradiometers. In contrast to that, the

gravity ﬁeld representation is formulated in an Earth-ﬁxed reference frame (ITRF)

with its origin located in the geocenter, z-axis directed to the pole, x-axis ﬁxed

in the equatorial plane in the direction of the Greenwich meridian, and y-axis as

completion to a right-handed system. Therefore, the observation equation has to be

transformed from the ITRF into the GRF.

In a real data analysis, the inﬂuence of non-gravitational forces have to be

corrected for. These disturbing forces, such as tides (direct tidal forces of sun

and moon, ocean tides, Earth tides polar tides) and short periodic variations of

atmospheric masses and the reaction of the ocean to theses masses, can be approx-

imated by appropriate models. In anticipation of being able to process real GOCE

observations in the near future, these background models have been implemented

in the software and have been used for real-data GRACE analysis as well. In

the simulation scenario described in Chap. 3 they have not been considered for

simpliﬁcation.

2.1 Basis Functions

As mentioned above, the parameterization of the gravitational potential is performed

in terms of spherical harmonics for the global solution, in case of a regional gravity

ﬁeld recovery the potential has to be parameterized by space localizing basis func-

tions. The gravitational potential V(r) at a given ﬁeld point r can be modeled as a

sum of basis functions as follows,

V(r) =



i=1



(r,r

), (1)

with the ﬁeld parameters a

arranged in a column matrix and the basis functions

located at the nodal point r

being denoted by 

(r,r

). These are chosen as isotropic

and homogeneous harmonic spline functions (Freeden et al., 1998) and can therefore

be expanded into a series of Legendre polynomials,



(

r,r

)



n=2





n+1

(

r,r

)

.(2)

R is the mean equator radius of the Earth, r the distance of a ﬁeld point from the

geo-centre, and P

(r,r

) are the normalized Legendre polynomials of degree n.

The coefﬁcients k

deﬁne the shape of the basis functions. In our approach these

shape coefﬁcients are chosen to ensure that the basis functions become decorrelated

258 A. Eicker et al.

with respect to the inner product in a reproducing kernel Hilbert space deﬁned by

the smoothness characteristics of the gravity ﬁeld as reproducing kernel, for more

details see Eicker (2008). This leads to the following choice for the coefﬁcients k

with σ

describing the degree variances of the gravity ﬁeld spectrum,

√

2n +1

with σ



m=−n

.(3)

The nodal points as locations for the basis functions are generated on a grid by

a uniform subdivision of an icosahedron of 20 equal-area spherical triangles. In

this way the global pattern of spline nodal points r

shows approximately uni-

form nodal point distribution as illustrated in Eicker (2008). For an overview of

regional gravity ﬁeld parameterizations used in satellite geodesy see for example

Schmidt et al. (2007).

The observation equation for gradiometer measurements is obtained by dif-

ferentiating the potential twice. It reads in case of a regional gravity ﬁeld

reﬁnement,

∇∇V(r) =



i=1

∇∇

(r,r

). (4)

2.2 Regionally Adapted Regularization

Due to the ill-posedness of the downward continuation process the solution has to

be regularized. Here the Tikhonov regularization has been applied and the regular-

ization parameter has been determined by means of variance component estimation

as proposed by Koch and Kusche (2001).

The process of variance component estimation yields the optimal regulariza-

tion parameter under consideration of the given signal-to-noise ratio. In case of

a regional gravity ﬁeld determination, this results in one regularization parame-

ter tailored optimally to the respective recovery region. This is an improvement in

comparison to a global gravity ﬁeld parameterization which allows only one regu-

larization factor for the complete Earth, resulting in an overall mean damping of the

gravity ﬁeld features. But even within smaller geographical areas, the gravity ﬁeld

features may vary signiﬁcantly. Therefore, it seems reasonable to further adapt the

regularization procedure. The proposed approach does not take into account only

one regularization matrix with one associated regularization parameter per region,

but allows several matrices with respective parameters

N =

PA +

+···+

.(5)

Here, N denotes the combined normal equation matrix. The ﬁrst term on the right

hand side contains the accumulated observation equations with the corresponding

Regionally Reﬁned Gravity Field Models from In-Situ Satellite Data 259

variance factor σ

. For each regional regularization group i, the regularization

matrix R

is a diagonal matrix featuring a “one” for each regional spline param-

eter located inside the corresponding region and a “zero” for parameters belonging

to basis functions outside the regularization group,

(j, j) =



1forj inside i

0forj outside i,

(6)

+···+R

= I.(7)

The possibility of adapting the regularization procedure in this particular way is

a unique feature of a ﬁeld parameterization by space localizing basis functions, as

each unknown parameter is related to a particular geographical location. This is an

inevitable premise when the elements of the regularization matrix are supposed to be

assigned to a certain region. The approach assumes that the original regularization

matrix used for a global calculation can be approximated by the unit matrix. This is

achieved by the speciﬁc choice of the space localizing basis functions as described

in Sect. 2.1 and the orthogonality relations emerging from it, as explained in Eicker

(2008). The separation of a geographical region into different regularization areas is

exemplarily illustrated in Fig. 1. Here the two regularization areas are t he continental

and the oceanic regions. The separation into land and ocean areas can propose a

reasonable choice in certain regions, where the gravity ﬁeld information over the

oceans is signiﬁcantly less rough compared to the continental areas. In the presence

Fig. 1 Two different regularization areas and the location of the unknown parameters belonging

to each area

260 A. Eicker et al.

of deep sea trenches or rough ocean bottom topography, however, this might not

always be a valid assumption. In this case, different choices for the regularization

areas are inevitable.

3 Simulation Scenario

Due to the lack of real GOCE observations, the following calculations were per-

formed in the frame of a simulation scenario. The simulation example is supposed

to demonstrate that a regional reﬁnement strategy cannot only be applied to reﬁne

global solutions, but it can be utilized to combine different data sets as well. In

the scenario at hand, GRACE-observations supposed to be combined with those of

GOCE. The GRACE solution serves as the global reference ﬁeld, as GRACE pro-

vides excellent results especially in the long and medium wavelength part of the

gravity ﬁeld spectrum. The regional reﬁnements of this solution are then calculated

on the basis of GOCE-like observations, as they cover the short periodic part of the

gravity ﬁeld spectrum with superior accuracy. As additional reference ﬁeld for the

higher frequencies, the OSU91 (Rapp et al., 1991) was introduced from degree n =

150 to n = 300. An extension of this combination approach is suitable to include

terrestrial or airborne gravity data as well.

In order to achieve a consistent data set, both the GRACE and the GOCE solution

were calculated from simulated observations on the basis of the EGM96 (Lemoine

et al., 1998) up to degree n = 300. The observations for GRACE as well as for

GOCE were simulated for a period of 30 days with a sampling rate of 5 s. In case

of the GRACE satellites, the simulated orbit positions were corrupted by a white

noise with a standard deviation of 3 cm and the intersatellite ranges between the

GRACE twin satellites with a white noise of 10 μm. Regarding the GOCE satellite,

the gradiometer observations were corrupted by a colored noise model with a stan-

dard deviation of 1.2 mE. In this model, the power spectrum is assumed as constant

in the measurement band of the gradiometer and features increasing energy in the

long wavelength part of the spectrum, as the gradiometer is not able to recover these

low frequencies. For the satellite positions again an accuracy of 3 cm was assumed.

From this simulated data set, ﬁrst the GRACE solution up to a harmonic degree of

n = 150 was derived. The analysis of the GRACE data was conducted using the

integral equation approach as described in Mayer-Gürr (2006). The regional reﬁne-

ments from GOCE observations were then calculated as residual ﬁelds to this global

solution. Using a parameterization by splines as space localizing basis functions, the

spline kernels are to represent the gravity ﬁeld features to be determined. Therefore,

up to degree n = 150, the error degree variances of the (simulated) GRACE solu-

tion were applied in Eq. (3) as coefﬁcients in the series expansion of the spline

kernel. These error degree variances represent the signal which is still in the data in

addition to the GRACE solution. Above degree n = 150, the degree variances were

approximated by Kaula’s rule of thumb. The resulting basis function is displayed

in Fig. 2.