Print this page

Functional and stochastic modelling of satellite gravity data

PoG 67, Jasper van Loon, Functional and stochastic modelling of satellite gravity data

Jasper van Loon

Publications on Geodesy 67
Delft, 2008. 248 pagina's. ISBN: 978 90 6132 307 5.


Summary

Functional and stochastic modelling of satellite gravity data

The gravitational field of the Earth is caused by the heterogeneous mass distribution of the Earth. Due to this field, the mean sea level (at rest) shows fluctuations (geoid heights) up to 100 meters with respect to a reference ellipsoid. The time-varying gravity signal is a source of information on the melting of the polar ice, the hydrological cycle and other mass variations within the Earth. Moreover, the gravity field of the Earth is needed to derive the physical (orthometric) heights and to estimate the ocean at rest, which is necessary for oceanographic research. The gravitational field is an harmonic field and can be expressed in a model of spherical harmonic coefficients. A higher resolution can be obtained by an increase in this number of coefficients; e.g., more than 130,000 coefficients are needed for a resolution of 111 kilometers. Traditional techniques for the determination of this gravity field are terrestrial measurements, marine measurements and airborne gravimetry. Since the launch of the first satellite in 1957, satellite tracking measurements are used to derive the gravitational field, together with altimetry measurements to estimate the sea surface at rest. With the launch of the dedicated satellite missions CHAMP (2000), GRACE (2002) and the upcoming GOCE mission (2008), millions of observations have become available, which will (and already have) improved the gravitational field by several orders. Moreover, the measurements of the GRACE mission enable us to derive monthly estimates of the mass changes in a thin shell near the Earth's surface, mainly due to hydrology, post-glacial rebound and the melting of the polar ice. The combination of different observation groups (satellite / surface gravity data) raises a number of research questions.

How to deal with systematic effects in the data?

The (linear) relationship between the observations and the unknown parameters is defined in the functional model. These parameters consist of the 'global' coefficients of the gravitational field and some 'local' parameters, such as systematic effects. The residuals, i.e., observations minus model, should only consist of observation noise. Local parameters need to be tested for significance. A correct stochastic model is needed for such an hypothesis testing. As unmodelled systematic effects enter into the stochastic model, which will change this stochastic model with each improvement of the vector of unknowns, this testing should be done in an iterative process.

How to weight the different observation groups?

The quality of the satellite measurements vary considerably in time due to changing GPS constellations and instrument conditions. The weighted least-squares approach will give more weight to high-quality observations. The precision of the observations is expressed in the variance-covariance matrix, as part of the stochastic model. A-priori estimates of these stochastic properties are often not reliable and in general too optimistic. To achieve precise and reliable estimates of the unknown parameters, it is therefore necessary to calibrate the stochastic model. In the 1970s and 1980s such a calibration was mainly done externally, i.e., by comparing the result to independent data. A breakthrough was the implementation of the subset solution method, which is an internal calibration method suggested by F.J. Lerch [1989]. However, with the launch of the CHAMP and GRACE missions, which have brought millions of gravity-related observations to the geodetic community, this community has made a step backwards with respect to the validation of the stochastic model. This is probably due to the huge dimensions of the least-squares problem. In this thesis, we have shown that even for large amounts of data (and unknown parameters), it is still possible to compute unbiased variance components and consequently observation weights. With the use of Monte Carlo simulations, we have re-written the equations in such a way that every least-squares software package can be used to compute the new weights of the different data sets. In this way, the stochastic model can absorb the varying quality of the satellite data and it enables an improved combination of different observation types.

How to decrease the influence of the blunders in the data?

Outliers are observations that do not follow the assumed functional and stochastic model. The residual of an outlier is much higher than one could expect from the stochastic model. The conventional approach is to test a bias in an observation for significance (w-test statistic) and consequently remove the observation with the highest test statistic (data snooping). An alternative method is to assume a different distribution than the normal distribution as the probability density function of the observations (M-estimation). In this way, more probability is assumed for observations with a high test statistic. Therefore, the influence of these observations is highly reduced. With a sufficiently high redundancy in the observations, it is even possible to estimate this probability density function from the observations and consequently re-weight all observations accordingly. This cost function estimation(CFE) produced the best results in a test setup with real CHAMP satellite gravity data. It is recommended to further refine this method and to use it on different sources of data, including the very precise measurements of the GOCE satellite.

The testing of the local parameters, the validation of the stochastic model and the detection and treatment of the outliers are all connected to each other. It is therefore recommended to perform these algorithms in an iterative process. The algorithms have been tested on several different applications in this thesis. In the first application a global gravity field model is derived from CHAMP pseudo-observations, using the energy balance approach. In the second test setup, we have computed a joint inversion of weekly GPS site displacements with monthly GRACE models to obtain estimates of the surface mass redistributions. Simulated GRACE observations are used in the third application to quantify and reduce the effect of temporal aliasing. In the last application a corrector surface to the gravimetric geoid is estimated for Switzerland using a combination of GPS/levelling data and the gravimetric geoid.


Contents

  1. Introduction  1
  2. Estimation of the Earth's gravity field  9
  3. Augmentation of the functional model  35
  4. Stochastic model validation  49
  5. Monte Carlo implementation  83
  6. Outlier detection and robust estimation  97
  7. Application 1: CHAMP satellite gravity data  115
  8. Application 2: Joint inversion of global GPS time-series with GRACE gravity models  141
  9. Application 3: Temporal aliasing of hydrological signals in a simulated GRACE recovery  165
  10. Application 4: The computation of a height reference surface in Switzerland  177
  11. Conclusions and recommendations  191

References  197

  1. Series expansion into spherical harmonics  217
  2. Matrix algebra and matrix analysis  219
  3. Some standard distributions  221

Summary  223
Samenvatting  227
Curriculum Vitae  231


Samenvatting

Functional and stochastic modelling of satellite gravity data

Het zwaartekrachtsveld van de Aarde wordt veroorzaakt door de massaverdeling binnen de Aarde en heeft tot gevolg dat het zeeoppervlak in rust afwijkingen (geoïdehoogten) vertoont tot wel 100 meter ten opzichte van een referentie-ellipsoïde. Veranderingen in dit zwaartekrachtsveld geven ons zo informatie over de smelting van de poolkappen, de hydrologische cyclus en de massaveranderingen in het binnenste van de Aarde. Bovendien is het zwaartekrachtsveld essentieel in de bepaling van de fysische (orthometrische) hoogte en de definiëring van het zeeoppervlak in rust, hetgeen nodig is voor oceanografisch onderzoek. Dit zwaartekrachtsveld is een harmonisch veld en kan uitgedrukt worden in een model van spherisch harmonische coëfficiënten. Een verhoging van dit aantal coëfficiënten leidt tot een verbetering van de resolutie. Zo heeft men voor een resolutie van 111 kilometer meer dan 130.000 coëfficiënten nodig. Traditionele technieken voor de bepaling van dit zwaartekrachtsveld zijn zwaartekrachtmetingen op land, mariene metingen en vliegtuiggravimetrie. Sinds het begin van het ruimtetijdperk (1957) worden ook afstandsmetingen naar satellieten uitgevoerd, alsmede altimetriemetingen van het zeeoppervlak. Met de lancering van de satellietmissies CHAMP (2000), GRACE (2002) en de toekomstige GOCE missie (2008) komen vele miljoenen metingen beschikbaar, die ervoor zorgen (en al hebben gezorgd) dat het zwaartekrachtsveld met enkele ordes beter bepaald kan worden. De waarnemingen van de GRACE-missie maken het bovendien mogelijk om voor het eerst op een maandelijkse basis de massaveranderingen in de bovenste laag van de Aarde te berekenen. Deze massaverplaatsingen zijn voornamelijk toe te schrijven aan hydrologie, post-glaciale opheffing en het smelten van de poolkappen. De combinatie van verschillende waarnemingsgroepen (satellietmetingen / terrestrische metingen) levert een aantal onderzoeksvragen op.

Hoe om te gaan met de systematische effecten per data set?

Het functioneel model geeft de relatie weer tussen de waarnemingen en de onbekende parameters. Deze onbekende parameters bestaan uit de 'globale' coëfficiënten van het zwaartekrachtsveld en enkele 'locale' parameters, zoals systematische fouten. De residuen, oftewel waarnemingen minus model, dienen idealiter slechts de ruis te bevatten. De locale coëfficiënten dienen getest te worden op significantie. Hiervoor is echter ook het stochastische model nodig. Aangezien ongemodelleerde systematische fouten geabsorbeerd worden door het stochastic model, zal dit model bij iedere schatting van de onbekende parameters aangepast moeten worden. Hierdoor is het testen van de locale parameters een iteratief proces.

Hoe moeten de verschillende waarnemingsgroepen gewogen worden?

De kwaliteit van satellietmetingen verandert sterk in de tijd door veranderende GPS-constellaties en het functioneren van de meetapparatuur. In een gewogen kleinstekwadraten oplossing krijgen metingen met een hogere kwaliteit een hoger gewicht in de schatting van de onbekende parameters. De kwaliteit wordt uitgedrukt in de variantiecovariantie matrix, als onderdeel van het stochastisch model van de waarnemingen. De a priori schattingen van deze stochastische eigenschappen zijn meestal niet betrouwbaar en over het algemeen te positief. Een calibratie van het stochastisch model is daardoor noodzakelijk voor een precieze en betrouwbare schatting van de onbekende parameters. In de jaren '70 en '80 werd deze calibratie voornamelijk uitgevoerd door het eindproduct te vergelijken met onafhankelijke gegevens. Een doorbraak was de 'subset solution'-methode, een interne calibratiemethode, voorgesteld door F.J. Lerch [1989]. Met de komst van miljoenen waarnemingen van de nieuwe satellietmissies CHAMP en GRACE heeft de geodetische wereld echter een stap teruggedaan voor wat betreft de validatie van het stochastische model. Dit is waarschijnlijk toe te schrijven aan de enorme dimensiteit van het kleinstekwadraten probleem. In dit proefschrift laten we echter zien dat ook voor grote hoeveelheden data het mogelijk is om de zuivere variantiecomponentenschatting toe te passen. Met behulp van Monte Carlo simulaties zijn de vergelijkingen zo opgeschreven dat ieder kleinstekwadraten softwarepakket gebruikt kan worden voor de berekening van de nieuwe gewichten van de verschillende data sets. Op deze manier kan de variërende kwaliteit van de satellietgegevens worden opgevangen door het stochastisch model en is een betere combinatie mogelijk tussen de verschillende soorten waarnemingen.

Hoe kan de invloed van blunders zo veel mogelijk beperkt worden?

Blunders zijn waarnemingen, die niet overeenkomen met de aangenomen functionele en stochastische modellen. De residuen zijn veel groter dan verwacht zou worden, uitgaande van het stochastisch model. De conventionele aanpak is om een afwijking in de waarneming te testen voor significantie (w-toets) en vervolgens de waarneming met de hoogste testgrootheid te verwijderen uit de data set (data snooping). Een alternatieve methode is echter om een andere verdeling voor de waarnemingen aan te nemen dan de normaalverdeling. Hierdoor vermoed men meer kans op een blunder in de waarnemingen en verkleint men zo de invloed van zo'n blunder. Nog een stap verder is om deze kansverdeling rechtstreeks te schatten uit de waarnemingen. Deze 'cost function estimation (CFE)' levert in een testopzet met echte CHAMP-waarnemingen de beste oplossing. Het verdient aanbeveling om deze methode te verfijnen en toe te passen op verschillende soorten data, waaronder de nieuwe, zeer nauwkeurige metingen van de GOCE-satelliet.

Het testen van locale parameters, het verbeteren van het stochastisch model en de behandeling van blunders in de data zijn gerelateerd aan elkaar. Het is daarom aan te raden dit in een iteratief proces op te lossen. Dit is in dit proefschrift gedaan voor meerdere applicaties betreffende zwaartekrachtgegevens. Zo is het globale zwaartekrachtsveld gemodelleerd, gebruikmakend van satellietmetingen van de CHAMP-missie. De tweede applicatie vormt de combinatie van wekelijkse GPS-stationsverplaatsingen met maandelijkse GRACE-modellen voor de bepaling van massaverplaatsingen in de opperste laag van de aarde. De derde toepassing richt zich op het reduceren van het 'temporal aliasing'-effect. Hierbij is gebruik gemaakt van gesimuleerde GRACE-waarnemingen. De laatste toepassing is de bepaling van het correctieoppervlak in Zwitserland door gebruik te maken van GPS-waarnemingen, waterpasmetingen en de gravimetrische geoïde.