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

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

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

Fig. 2 Spline kernel used in the combined GRACE-GOCE simulation scenario: developed until

degree n = 300, formal errors of simulated spline solution up to degree n = 150, above n = 150

Kaula’s rule

The basis functions were located at the nodes of a triangular grid with an average

nodal point distance of about 74 km. This results in a resolution slightly higher

than a spherical harmonic degree of n = 300. From the setting described above

regional solutions were calculated. The results for an individual regional reﬁnement

are displayed in Fig. 3 for the area of the Andes. The ﬁgure shows the differences

to the pseudo-real ﬁeld EGM96 in terms of geoid heights.

The comparison was only performed up to degree n = 240, because t he higher

degrees are too strongly corrupted by noise. The left part of the ﬁgure shows the

result with one uniform regularization parameter determined for the region, the

Fig. 3 Differences between spline solution from combined GRACE/GOCE analysis (simula-

tion study) and the pseudo-real ﬁeld EGM96: one uniform regularization parameter (left,RMS:

9.24 cm), two adapted regularization parameters for land and ocean (right,RMS:8.08cm)

262 A. Eicker et al.

right part demonstrates the improvement that was achieved using the regionally

reﬁned regularization procedure described in Sect. 2.2 with adapted regularization

for ocean and continental area. The RMS value is r educed from 9.24 to 8.08 cm.

This corresponds to an improvement of 12.5%.

Subsequently, regional reﬁnements were calculated covering the complete sur-

face of the Earth, with an overlapping border of 10

◦

having been taken into account

at the boundaries to avoid truncation effects. The individual patches have a size of

λ = 40

◦

× ϕ = 50

◦

, complemented by two spherical caps with an aperture

angle of 30

◦

covering the poles. The pattern of the reﬁnement regions is illus-

trated in Fig. 4. It leads to a number of about 5,000–9,000 spline parameters to

be determined for each region, the size of the patches being limited by storage

restrictions.

From the individual spline solutions the gravity ﬁeld functionals were predicted

to the nodes of a Gaussian grid with a spacing of λ = 0.5

◦

. From the global ﬁeld,

a spherical harmonics expansion was calculated using the Gauss-Legendre quadra-

ture, for the speciﬁc grid and the respective quadrature method see, for example,

Stroud and Secrest (1966) or Sneeuw (1994). For the global solution, again com-

pared up to degree n = 240, the global RMS amounts to 6.51 cm, including the

poles. Only the solution with the adapted land/ocean regularization is presented in

Fig. 4, as on a global scale the differences between the use of one or two regulariza-

tion parameters are not easily detectable. The problem arises that for many of the

patches the choice of continental and oceanic regions does not pose an ideal option.

Here a more tailored adjustment of the regularization areas would be advisable.

Fig. 4 Individual reﬁnement patches, differences in geoid heights between combined

GRACE/GOCE solution and EGM96, RMS for global coverage: 6.51 cm

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

4 Conclusions

The concept of taking a global solution and calculating regional reﬁnements to

the global solution presents an elegant method of combining different data sets.

The lower resolution data set can serve as global reference, and observations with

assumed higher resolution can be exploited for the regional reﬁnements. The reg-

ularization parameters are determined individually for different geographical areas

by variance component estimation. With this procedure the regularization parame-

ter is optimally chosen according to the given signal-to-noise ratio in the respective

regularization area. In this way the dampening can be adjusted in accordance with

the residual signal content (i.e. the information additional to the global reference

solution).

It shall be pointed out that the direct computation of the spherical harmonic

coefﬁcients by solving the improperly posed downward continuation cannot provide

a stable solution up to an (arbitrarily) high degree as it can be achieved by means

of the quadrature method. Finally, the regional recovery procedure offers a chance

to deal with the polar gap problem in a tailored way, as in regions without any data

the regional reﬁnement can either be skipped or the regularization parameter can be

adjusted accordingly.

Acknowledgements The support by BMBF (Bundesministerium fuer Bildung und Forschung)

and DFG (Deutsche Forschungsgemeinschaft) within the frame of the Geotechnologien-Programm

is gratefully acknowledged.

References

Eicker A, Mayer-Gürr T, Ilk KH (2004) Global gravity ﬁeld solutions based on a simulation sce-

nario of GRACE SST data and regional reﬁnements by GOCE SGG observations. Proceedings

of the International Conference Gravity, Geoid and Space Missions – GGSM2004, 2004 August

30–September 3, Porto, Portugal.

Eicker A (2008) Gravity ﬁeld reﬁnement by radial basis functions from in-situ satellite data.

Dissertation University of Bonn. URN: urn:nbn:de:hbz:5N–13754, URL: http://hss.ulb.uni-

bonn.de/diss_online/landw_fak/2008/eicker_annette.

ESA (1999) Gravity ﬁeld and steady-state ocean circulation mission. Reports for Mission

Selection, ESA, SP-1233(1). European Space Agency Publications Division, Nordwijk.

Freeden W, Gervens T, Schreiner M (1998) Constructive Approximation on the Sphere. Oxford

University Press, Oxford.

Koch KR, Kusche J (2001) Regularization of geopotential determination from satellite data by

variance components. J. Geod. 76(5), 259–268.

Lemoine FG, Kenyon SC, Factor JK, Trimmer RG, Pavlis DS, Chinn NK, Cox CM, Klosko SM,

Luthcke SB, Torrence MH, Wang YM, Williamson RG, Pavlis EC, Rapp RH, Olson TR (1998)

The Development of the Joint NASA GSFC and NIMA Geopotential Model EGM96. NASA

Goddard Space Flight Cent., Greenbelt, MD.

Mayer-Gürr T (2006) Gravitationsfeldbestimmung aus der Analyse kurzer Bahnbögen am

Beispiel der Satellitenmissionen CHAMP und GRACE. Dissertation at the University

of Bonn. URN: urn:nbn:de:hbz:5N–09047, URL: http://hss.ulb.uni-bonn.de/diss_online/

landw_fak/2006/mayerguerr_torsten

264 A. Eicker et al.

Rapp RH, Wang YM, Pavlis NK (1991) The Ohio State 1991 geopotential and sea surface topog-

raphy harmonic coefﬁcient models. Number 410 in Reports of the Department of Geodetic

Science. Ohio State University (OSU), Columbus, OH.

Reigber C, Schwintzer P, Lühr H (1999) The CHAMP geopotential mission. Boll. Geof. Teor.

Appl. 40, 285–289.

Schmidt M, Fengler M, Mayer-Gürr T, Eicker A, Kusche J, Sanchez L, Han SC (2007) Regional

gravity modeling in terms of spherical base functions. J. Geod. 81(1), 17–38.

Sneeuw N (1994) Global spherical harmonic analysis by least squares and numerical quadrature

methods in historical perspective. Geophys. J. Int. 118, 707–716.

Stroud AH, Secrest D (1966) Gaussian Quadrature Formulas. Prentice Hall, Englewood Cliffs, NJ.

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

experiment: mission overview and early results. Geophys. Res. Lett. 31, L09607, doi:

10.1029/2004GL019920.

Quality Evaluation of GOCE Gradients

Jürgen Müller, Focke Jarecki, Insa Wolf, and Phillip Brieden

1 Cross-Over Analysis

The ﬁrst ESA Earth Explorer Core Mission GOCE (Gravity Field and Steady-

State Ocean Circulation Explorer) entered the operational measurement phase in

September 2009. Before gravity ﬁeld processing, the quality of the GOCE gradi-

ents has to be assessed. Here, two procedures have been developed in Hanover,

the mutual comparison and analysis of observed gradients in s atellite track cross-

overs and the application of terrestrial gravity data which are upward continued

and transformed into reference gradients for the GOCE observations. As shown in

Jarecki et al. (2006), comparison of the gravity tensor elements in satellite ground

track cross-overs (“XO analysis”) offers a good opportunity to validate the gra-

dients measured by a satellite gravity gradiometer. There, closed loop results are

presented, with accuracies ﬁtting well into the margins of the overall GOCE gra-

diometer performance. The radial tensor component can be reproduced with an

accuracy of better than 3 mE, just using the un-ﬁltered measurements in correspond-

ing track parts and a recent high-resolution geopotential model like the combined

GRACE/terrestrial model (e.g. EIGEN-GL04C, Förste et al., 2006). The XOs are

identiﬁed by a measurement oriented search algorithm (independent of the predicted

orbit characteristics) resulting in the actual position of the ground track crossing.

Gradiometer measurements and auxiliary data are interpolated into the position of

the XO along each of the crossing orbit tracks. Gradient differences from differing

orbit altitudes and orientations of the measurement reference frame are reduced by

applying the geopotential model (and auxiliary data). These closed loop tests show,

that residual gradient differences in the XO can be assigned to measurement errors,

as they should vanish in the absence of noise and systematic errors. Such noise

and systematic errors have been modelled for the integrated gradiometer (single

components of the system, like speciﬁc accelerometers, are not treated separately)

J. Müller (B)

Institut für Erdmessung, Leibniz Universität Hannover, 30167 Hannover, Germany

e-mail: mueller@ife.uni-hannover.de

265

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_21,



Springer-Verlag Berlin Heidelberg 2010

266 J. Müller et al.

by Koop et al. (2002), with model reﬁnements by Denker et al. (2003), as

meas

(t) = V

real

(t)+V

∗

+o(t) V

∗





k=1

cos (kωt)+b

sin (kωt)+ε(t). (1)

In the XOs, differences of the measured gradients are further investigated. But

not all modelled systematic errors project into those XO differences according to

Eq. (2). Temporally constant or, due to the relation to a discrete point on the earth

sphere, locally constant features V

∗

vanish. The ﬁrst vanishes because of the rela-

tive nature of this internal validation approach, the latter could also be characterised

as some kind of geographical signal (and therefore not as measurement error).

Constant biases V

∗

, which are not present in the whole time series investigated

o(t) = const, can project into the XO differences. A linear trend V



projects lin-

early (with the measurement time difference t

= t

–t

in the XO) into the

investigated value V

meas

. Fourier harmonics a

cos(kωt)+b

sin(kωt) sum up

to combined terms and are still present in the gradient differences. Contrary to the

trend (and polynomial parameters of higher degree) where only time differences

t

show up, they carry absolute time information into the differences, modelled

as t

and t

. The noise parameter (t) represents the stochastic noise. Thus, the

following systematic parts of the gradients’ differences remain in the XOs:

V

meas



o(t

) −o(t

)



∗

+ V



t



k=1



cos (kωt

) −cos (kωt

)



+ b



sin (kωt

) −sin (kωt

)



(2)

Consequently, the XO approach is (only) capable to detect three kinds of sys-

tematic errors, namely linear trends (and polynomial drifts), Fourier-type periodic

disturbing effects and biases which are only present in short parts of the data (called

Short Term Biases). The latter may reach from single outliers (affecting just one

special measurement) to larger portions of the data set under investigation.

In the following sections, it will be discussed if and under which conditions those

systematic errors can be detected, and which approaches might be used to get robust,

accurate and reliable results. The studies are shown for the tensor component V

for the ease of formulation and readability, but similar results are achieved for the

other diagonal and off-diagonal tensor components.

1.1 Short Term Biases

Biases which are not present in the whole data set (called sequential outliers) project

into XO gradient differences in contrast to biases affecting the whole data set. In the

ﬁrst case there is some data left which can be considered as “reference data”. The

amount and the coverage (and, of course, the quality) of reference data is the main

criterion to decide, whether the detection of outliers and short term biases works: On

the one hand, there has to be enough reference data to intersect with the erroneous

part of the track at all. On the other hand, shorter erroneous parts of the track (leaving

more correct reference data) might be hidden between the single XOs, especially

Quality Evaluation of GOCE Gradients 267

Table 1 Single and sequential outliers as detected with the XO method when a threshold is a pplied

to the XO differences of the tensor component V

. The 30 d 1 Hz input data have been superposed

by 4,400 single outliers of 10 mE affecting 1,633 XOs and by three series of sequential outliers (10

and 15 mE) affecting 6,045 XOs. Type 1 error indicates the rate of theoretically detectable outliers

not detected, while type 2 error means false alarms

No. of alarms Error (%)

Threshold Total Correct Type 1 Type2

Single outliers

1.5 mE 1,773 1,399 14.3 22.9

2 mE 1,508 1,152 29.5 21.8

Sequ. outliers

1.5 mE 6,018 6,001 0.72 0.28

2 mE 6,007 6,001 0.72 0.09

in the mid latitudes where the XO coverage is sparse. The latter problem becomes

especially clear, when dealing with single outliers. Those single erroneous gradients

cannot be detected by the XO method, if they appear in parts of the track, where no

crossing arc is present. Moreover, even if they affect an XO, it cannot be determined

which of the four (in case of linear interpolation) measurements contributing to the

XO is the erroneous one.

The capability of the XO approach can simply be tested by applying a threshold

to gradient differences, where artiﬁcial outliers in single values of V

are added or

is disturbed constantly during single revolutions of the satellite track. In a test

case ( see Table 1) roughly 4,400 single outliers were artiﬁcially introduced in a sim-

ulated30d1HzV

GOCE test data set. Of the 1,633 affected XOs, nearly 30%

remained undiscovered or led to false alarms. While the errors of type 1, undiscov-

ered outliers, can be reduced by carefully choosing the threshold ( and most likely by

using a more sophisticated test method), the type 2 errors, or false alarms, seem to

remain in the same range when changing the threshold. This is due to the threshold

being chosen close to the accuracy of the measurements themselves. The detection

of sequential outliers works much more effectively and is not as sensitive to the

choice of the threshold. This is mainly due to the structure of the outliers, which

affect all the measurements used in the XO interpolation and therefore severely

project into the XO difference.

For single outlier detection, the detection mode should be more sophisticated to

signiﬁcantly identify most outliers. Another important improvement would be the

identiﬁcation of a single erroneous measurement rather than a suspicious XO, e.g.,

by changing the interpolation algorithm systematically. Those reﬁnements were not

planned for the present study, as other methods, like along-track interpolation (see

Bouman and Koop, 2003) show complementary behaviour and are recommended

(and adopted in the ofﬁcial GOCE processing) for use in outlier detection.

1.2 Trend

According to Eq. (2) a linear trend projects linearly into the XO gradient differ-

ences with respect to the measurement time differences and independently of the

268 J. Müller et al.

Table 2 Statistics of estimated trend factors for all 475 ascending arcs from least-squares and

from robust estimation on the basis of XO differences from different time spans

XOs with

tracks from

Mean trend

(mE/d)

Abs. error

(mE) Std (mE/d)

Min trend

(mE/d)

Max trend

(mE/d)

30 days Least squares 6.7979 0.0005 0.0014 6.7033 6.9407

Robust 6.7974 0.0000 0.0001 6.7967 6.7980

10 days Least squares 6.7994 0.0020 0.0330 6.6339 7.1368

Robust 6.7974 0.0000 0.0001 6.7971 6.7978

1 day Least squares 6.7992 0.0018 0.0168 6.7742 6.9150

Robust 6.7975 0.0001 0.0001 6.7972 6.7977

actual measurement epochs. Therefore, all XOs contribute to the determination of a

single (constant) trend factor. Consequently, its determination is possible with high

accuracy and reliability, when large data sets are used. Dealing with shorter data

sets, single outliers might severely affect the standard least-squares estimation of a

trend parameter. This is a typical example for the application of robust estimation

methods (Huber, 1981). The formal application of the Huber approach for trend

estimation from XO differences has been described in Jarecki and Müller (2009).

There, a robust estimator is successfully applied to estimate (artiﬁcial) trend fac-

tors for single and half revolution GOCE data sets. Table 2 summarises the beneﬁts

of robust trend estimation with respect to short reference data sets. Least-squares

estimation gives reasonable (mean) results and a good standard deviation for the

trend factor (introduced as 6.7974 mE/d) when using all 30 d of test data. Taking

just the XOs from the surrounding day(s) produces larger deviations and absolute

estimation errors. However, trend factors from robust estimation ﬁt very well and

deviate less. In Fig. 1, the estimated trend factors for each of the 475 ascending

arcs of a 30 d test data set is shown, as estimated with the least-squares and robust

Fig. 1 Estimated t rend factor for each of the 475 ascending arcs of a 30 d test data set with an

artiﬁcial trend of 6.7974 mE/d superposed. The estimates are based on XO differences from 1 day

of surrounding data. Robust results are shown in black, least-squares results in grey

Quality Evaluation of GOCE Gradients 269

approach, using 16 XOs. While several least-squares estimates obviously suffer

from outliers in the reference data, the robust estimates look smooth and correct.

A change in trend factors would easily be observable in the timeline of arc-wise

estimates.

The long-term trend estimation assesses the overall coherence of the data set, but

neglects time-variable physical effects on the measurements. Repeated determina-

tion of trend factors for short parts of the data and their comparison is a fast and easy

method to monitor t he behaviour of the measurement system. Robust trend estima-

tion provides those features and will be applied in the GOCE project as part of the

validation activities of the CalVal-Team (ESA, 2006).

1.3 Fourier Coefﬁcients

Periodic systematic errors do not project directly into the gradient differences, i.e.

their effect is not independent from the measurement epoch and therefore they give

a different feature in XO gradient differences. Nevertheless, they can be represented

in a linear model (assuming ﬁxed frequencies) and estimated as well by standard

least-squares or with robust methods. Thus, the results from the trend part regarding

overall data quality assessment and instrument monitoring also hold for low-degree

Fourier coefﬁcients, e.g., up to degree four (Bouman et al., 2004). Medium and

high degree coefﬁcients on the other hand, which would be valuable to gain infor-

mation about the gradiometer’s performance in its highly accurate measurement

bandwidth MBW (5–100 mHz), are not accessible via this way. As all coefﬁcients

up to the desired degree (100 mHz means approx. n = 500 based on complete

cycles) have to be estimated in the complete Fourier decomposition, large num-

bers of XO differences had to enter the computation. Short data sets as adopted

for the revolution-wise trend monitoring do not deliver such a big amount of data.

Consequently, a repeated estimation of Fourier coefﬁcients in the MBW is not pos-

sible. But in principle, global assessment of the gradiometer performance in the

MBW is possible from rather long data sets. This is illustrated in Fig. 2, where the

amplitudes for the (revolution based) Fourier frequencies 1–110 cpr are plotted, as

estimated from 500 or 20,000 XO differences, respectively. In a pre-processing step,

both, the simulated measurements of V

, which have been superposed by a pertur-

bation of 100 cpr which is in the MBW, and the applied reductions have been ﬁltered

with a bandpass butterworth ﬁlter to cut out the signal content in the MBW. To ﬁl-

ter the time-independent reductions, they have been computed from differences of

another time series of simulated V

. This is a straightforward method to deal with

the problem of ﬁltering the reductions. Figure 2 shows the disturbing signal clearly

and signiﬁcantly detected even with less “reference data”. But to rely on the abso-

lute amplitude estimated and to suppress wrongly estimated signals in the longer

wavelengths, much more data is needed. Finally, even the amplitude calculated from

20,000 XO differences is off by about 10% but at least the (non-existing) signals in

the other frequencies vanish.

270 J. Müller et al.

Fig. 2 Fourier amplitudes (1–110 cpr) estimated from 500 (grey) and 20,000 (black) V

XO dif-

ferences with standard least squares (dotted) and robust methods (solid lines). Measurement and

reduction data have been ﬁltered t o cover the MBW. A superposed disturbing signal of 100 cpr is

detected signiﬁcantly, even in case of sparse XO data while the lower frequencies are misﬁts (note

the logarithmic scale)

2 Accuracy Analysis of External Reference Gradients

in the Frequency Domain

To prove the planned accuracy level, the gradiometer measurements will be cali-

brated and validated internally and externally. Internal calibration will be performed

by industry and ESA. External evaluation is done by comparing measured gradients

with reference data, where we worked in the frequency domain, as it easily allows a

comparison with the given baseline accuracy. Here, the highest accuracy level of the

gradiometer within the MBW is 11 mE/

√

Hz, valid for the sum of the main diagonal

components of the gravitational tensor (1 mE = 10

–12

1/s

). We took this value as

worst case assumption for single tensor components, too.

One strategy of an external evaluation includes the use of gravity data on ground

upward continued to satellite altitude (e.g., Bouman et al., 2004; Denker, 2003;

Pail, 2003; Wolf, 2007). The evaluation can only be applied regionally, because

sufﬁciently accurate ground data are only available for selected areas. Therefore,

a global potential model, representing the long wavelength part, is introduced in

a remove-restore procedure. The error estimation for the external reference data is

carried out empirically in a synthetic environment permitting closed-loop validation

in all points.

2.1 Spectral Combination Method

In the present study, the radial tensor component V

is computed using integral

formulas, see Wolf (2007) for the other ﬁve components. In the following the