Royal Institute of Technology
Complementary studies in Satellite Gravity Gradiometry
Mehdi Eshagh
Post-doctoral report in Geodesy
Royal Institute of Technology (KTH) Division of Geodesy
10044 Stockholm Sweden December 2009
Acknowledgment
I would like to appreciate Professor Lars E. Sjöberg for providing the opportunity of working with him and the post-doc/research associate position in Division of Geodesy at Royal Institute of Technology. The Swedish National Space Board (SNSB) is cordially acknowledged for the financial support of project no. 63:07:1.
Contents
Summary………...1
1.Background………1
1.1 Topographic and atmospheric effects……….………1
1.2 Local gravity field determination………2
1.3 Variance component estimation………..5
1.4 Validation of satellite gradiometry data and least-squares modification…………6
1.5 The polar gaps……….7
1.6 Special studies of geoid and Moho discontinuity………...8
2. The Author’s contributions……….11
3. Obtained results and conclusions………15
References………...18
Paper A:
Eshagh M. (2010) On the convergence of spherical harmonic expansion of topographic and atmospheric biases in gradiometry, Contributions in Geophysics and Geodesy, (Accepted).
Paper B:
Eshagh M. (2009) Contribution of 1st-3rd order terms of a binomial expansion of topographic heights in topographic and atmospheric effects on satellite gravity gradiometric data, Artificial Satellites, Vol. 44, No. 1, P. 21-31.
Paper C:
Eshagh M. (2009) The effect of lateral density variation of crustal and topographic masses on GOCE gradiometric data: A study in Iran and Fennoscandia, Acta Geodaetica et Geophysica Hungarica, Vol. 44, No. 4, P. 399-418.
density in topographic effect on satellite gravity gradiometric data, Acta Geophysica, (Accepted). (0.3)
Paper E:
Eshagh M. and Sjöberg L. E. (2010) Determination of gravity anomaly at sea level from inversion of satellite gravity gradiometric data, Journal of Geodesy (Submitted).
Paper F:
Eshagh M. (2010) spatially restricted integrals in gradiometric boundary value problem, Artificial Satellites. (submitted)
Paper G:
Eshagh M. (2010) On integral approach to regional gravity field modelling from satellite gradiometric data, Acta Geophysica (Submitted).
Paper H:
Eshagh M. (2010) Variance component estimation in linear ill-posed problems:
TSVD issue, Acta Geodaetica et Geophysica Hungarica, Vol. 45, (Accepted).
Paper I:
Eshagh M. (2010) Error calibration of quasi-geoidal, normal and ellipsoidal heights of Sweden using variance component estimation, Contributions to Geophysics and Geodesy (Submitted).
Paper J:
Eshagh M. (2010) Least-squares modification of extended Stokes’ formula and its second-order radial derivative for validation of satellite gravity gradiometry data, Journal of Geodynamics (Accepted).
Paper K:
Eshagh M. (2010) Least-squares modification of second-order radial derivative of Abel-Poisson’s formula for generating satellite gravity gradiometry data, Computational Geosciences (Submitted).
Eshagh M. (2009) Least-squares modification of Stokes' formula with EGM08, Geodesy and Cartography, Vol. 35, No. 4 (in press).
Paper M:
Eshagh M. (2009) The effect of polar gaps on the solutions of gradiometric boundary value problems, Artificial Satellites, Vol. 43, No. 3, P. 97-108.
Paper N:
Eshagh M. (2010) Gravity anomaly recovery in polar gaps from full tensors of gravitation at satellite level, Acta Geodaetica et Geophysica Hungarica (Submitted).
Paper O:
Sjöberg L.E. and Eshagh M. (2010) Considering data gaps in geoid modelling by modifying Stokes's formula, Acta Geodaetica et Geophysica Hungarica, Vol. 45, (Accepted).
Paper P:
Sjöberg L.E. and Eshagh M. (2009) A geoid solution for airborne gravity data, Studia Geophysica et Geodaetica. Vol. 53, P. 359-374.
Paper Q:
Eshagh M., Bagherbandi M. and Sjöberg L.E. (2010) A combined global Moho model based on seismic and gravimetric data, Journal of Geodynamics (Submitted).
Summary
The gravity field and steady-state ocean circulation explorer (GOCE) (ESA 1999, Albertella et al. 2002, Balmino et al. 1998 and 2001) was finally launched on 17th March in 2009. In this satellite mission the second-order derivatives of the Earth’s gravitational potential are measured based on differential accelerometer at a satellite- borne gradiometer. It is expected to recover the geopotential coefficients to higher degrees and orders than those were obtained from the former satellite missions; say up to degree and order 300. Such an Earth’s gravity model will have an accuracy of 1 cm in global geoid height and of 1 mGal for the gravity anomalies, which are extremely good accuracies of the long-wavelength structure of the gravity field.
The present report is a summary of the studies of Mehdi Eshagh when a post- doct/research associate position in division of Geodesy was available for him. The research work consists of the studies continuing his thesis work and completing project no. 63/07:1, funded by the Swedish National Space Board (SNSB). The report is of collective papers type with a series of Papers A-Q. In the following we present summaries of his complementary studies. We just present the background of each study and the author’s contributions in comparing his research to others’ works.
The details about the methodology, theory, numerical investigations and conclusions are given in the corresponding papers of each subject.
1. Background
Generally, the research work is divided into six parts: topographic and atmospheric effects, local gravity field determination, validation, variance component estimation, the polar gaps and some special studies on geoid and Moho discontinuity. The background of each part is given in the following.
1.1 Studies in Topographic and atmospheric effects
The satellite gravity gradiometry (SGG) data are affected by the topographic masses.
The topographic effect has been considered by several geodesists (e.g. Martinec et al.
1993, Martinec and Vaníček 1994, Sjöberg 1998a, Sjöberg and Nahavandchi 1999, Tsoulis 2001, Heck 2003, Seitz and Heck 2003, Sjöberg 2000 and 2007) to determine a precise geoid. However, Wild and Heck (2004a and 2004b), Makhloof
and Ilk (2005 and 2006), Makhloof (2007), and Novák and Grafarend (2006), Eshagh and Sjöberg (2008a and 2009a) investigated the topographic effect on the SGG data.
Atmospheric effect was investigated by Ecker and Mittermayer (1969), Wallace and Hobbs (1977), Sjöberg (1993, 1998b, 1999, 2001 and 2006), Sjöberg and Nahavandchi (2000), Novak (2000), Tenzer et al. (2006). Eshagh and Sjöberg (2008, 2009a, 2009b and 2009c) and Eshagh (2009a and 2009b) have investigated the atmospheric effect on the SGG data. Sjöberg (2007) presented the topographic bias in analytical continuation and concluded that for considering the topographic effect on geoid, it is enough to consider the topographic Bouguer shell effect. After that, Vermeer (2008) commented on his paper and he used the spherical harmonics to prove that Sjöberg’s (2007) topographic bias is an approximation of the total topographic effect. Sjöberg (2008), in response of Vermeer (2008), mentioned that the spherical harmonic expansion using a binomial expansion of topography to fourth term is not convergent and the spherical harmonic expansion for the topographic bias will not be a correct way to disprove his results. Later on Sjöberg (2009a and 2009b) theoretically proved that the topographic bias is exact and even he proved that there is no need to consider the terrain correction in geoid computation.
Considering a constant density (2.67 g/cm3) for the topographic masses introduces errors in the computational harmonic space before downward continuation of the SGG data to geoid. Most of the efforts in this subject was done for geoid determination aspects, see e.g., Martinec (1993), Martinec et al. (1995), Kuehtreiber (1998), Pagiatakis and Armenkis (1999), Tziavos and Featherstone (2000), Huang et al. (2001), Hunegnaw (2001). Sjöberg (2004) showed that the total effect on the geoid due to the lateral density variations inside the masses at deepest and highest parts of the Earth is about 1.5 cm and 1.78 m, respectively. Kiamehr (2006) showed that the lateral density variation effect can reach to 30 cm in Iran. All the efforts of considering the effect of variations were in geoid determination process to recover high frequencies of the geoid. However, the high frequencies are considerably attenuated at satellite level, therefore considering general patterns of the lateral variations of the masses could be useful to obtain some ideas about the magnitude of their effect. Makhloof (2007) showed this attenuation numerically and concluded that the spherical harmonic expansion of the topographic potential is suitable for considering the topographic effect at satellite level.
1.2 Local gravity field determination
The second-order partial derivatives of the extended Stokes formula can be used to recover the gravity anomaly. In fact, the integral formulas are solved in reverse case by considering them as integral equations, namely the gravitational gradients are observed at satellite level and the unknowns are the gravity anomalies at sea level. In practice, these integral equations have to be discretized prior to computations. (The solutions of the equations yield the downward continued gravity anomaly.) The solutions are ill-posed, which implies (among other properties) that they are sensitive to the errors of the data. Any small error in the data can vary the solution considerably, and regularization is one method to stabilize such a system of equations.
Reed (1973) was one of the earliest geodesists who suggested direct downward continuation of the SGG data in a local gravitational field determination. He derived the corresponding kernels of the second-order partial derivatives of the extended Stokes formula and inverted the integrals to estimate the gravity anomaly at sea level. The results of the numerical studies were very promising. The ill-posed problems of downward continuation in his simulations were solved by defining some constraints for the gravity anomaly being estimated to the system of equations. After that, Rummel (1976) theoretically investigated the downward continuation of satellite data in stochastic and deterministic views. Zielinski (1975) considered least- squares collocation (Krarup 1969) for downward continuation and concluded that the covariance model is not local but global in statistical characters, and it is impossible to get rigorous results. Rummel et al. (1979) compared and presented the relation between regularization and least-squares collocation. Krarup and Tscherning (1984) evaluated the isotropic covariance function of torsion balance observations.
Tscherning (1988) investigated the recovery of gravity anomaly and geoid from SGG data at different altitudes, resolutions and noise levels and concluded that the geoid can be computed with an error of 10 cm. In a similar study, Tscherning (1989) studied the effect of sampling rate of the SGG data and concluded that it was possible to eliminate the effect of shortening the mission duration by improving the data sampling rate, the number of gravity gradients components and decreasing the noise standard deviation. Arabelos and Tscherning (1990) used some covariance functions to determine the local gravity field from SGG data. They concluded that the gravity anomaly can be recovered within 8.3 mGal and the geoid within 0.33 m.
Tscherning et al. (1990) investigated different methods of regional gravity field modelling from SGG data. They studied the Fourier, integral and least-squares collocation methods and concluded that removing the long-wavelength structure of the gravity field will smooth and decorrelate the data and thereby the error of the estimated qualities. Visser (1992) investigated downward continuation of satellite accelerations. He used Tikhonov regularization and least-squares collocation as deterministic and stochastic methods and concluded that collocation is superior to Tikhonov regularization. Tscherning (1993) derived the covariance of derivatives of the disturbing potential in a local frame and Arabelos and Tscherning (1993 and 1995) used these covariance functions to recover the gravity field locally from the observed SGG data in the orbital frame. Arabelos and Tscherning (1999) found 25%
improvement in recovering the gravity anomaly when the correlated errors are used rather than the uncorrelated errors in local gravity field determination from airborne gradiometric data. Least-squares collocation relies on a priori information of the unknowns, and Xu (1991) investigated the robustness of this method against incorrect prior information of data. Xu (1992) discussed the determination of the gravity anomaly from the SGG data without a priori information by using the Tikhonov regularization. He presented a regularization method minimizing the mean square error of the estimated anomalies. Xu and Rummel (1994) investigated the smoothness of the regional gravity field and concluded that the biased estimators improve the least-squares solution. They proposed an iterative ridge estimator. The truncated singular value decomposition was used by Xu (1998). van Gelderen and Rummel (2001) and (2002) presented the closed form expressions of Green’s functions of the gradiometric boundary value problem and removed the zero- and first-degree harmonics from their derivations. Martinec (2003) stated that one can
consider these degrees without problem for the vertical-vertical component and first- degree harmonics for the vertical-horizontal components. Both Green’s functions presented by van Gelderen and Rummel (2001) and (2002) and Martinec (2003) can be used to upward or downward continue the gravitational gradient to the mean orbital sphere. However, they cannot be used for downward continuation of the gradients to some related quantities at sea level. Toth et al. (2004) and (2006) investigated the downward and upward continuation of the gravitational gradient.
The difference between their approach and others is that they downward/upward continue the gradients to the same type of gradients at the mean orbital sphere. They applied a second-order Taylor expansion for this continuation and derived the corresponding kernel formulas in the spectral forms. Toth et al. (2007) investigated the practical aspects of this continuation using band-limited kernels and concluded that a cap size of at least 8D is needed to sufficiently suppress the omission errors for both upward and downward continuation with band pass filtered kernels. Kotsakis (2007) presented a covariance-adaptive method for regularization and used this method for determining the gravity anomaly at sea level from the second-order radial derivative of the gravitational gradients. This method is sensitive to the choice of covariance-adaptive matrix, which should be known a priori. Xu (2009) used an iterative generalized cross-validation method for simultaneous estimation of the regularization and variance components with an example of recovering the gravity anomaly from the SGG data. Janak et al. (2009) carried out the inversion of the SGG data in a spherical geocentric frame (see e.g. Koop 2003, Novak and Grafarend 2006 or Eshagh 2008) to the gravity anomaly at sea level using the truncated singular value decomposition for regularizing the system of equations being inverted.
In the gradiometric boundary value problem three combinations for the tensor elements (partial second-order derivatives) are constructed so that tensor spherical harmonics can be used to make integral formula for the function (the Earth’s gravity field); see for instance Rummel (1997), van Gelderen and Rummel (2001) and (2002) and Martinec (2003). If we assume that the boundary is a sphere, then the integration should be performed on it. However, if the integration is performed in a small part of the sphere, the function will not correctly be obtained and the solution is biased. Such a problem was considered for satellite altimetry data, as we do not have a full global coverage of satellite data. The integration will not produce the true function because the spherical harmonics are not locally orthogonal. Mainville (1986) was one of the earliest persons who presented recursive formulas for the integral product of the associated Legendre functions and used them in satellite altimetry problem. Later on Hwang (1991) presented another recursive formula for generating the gaunt coefficients (Varshalovich et. al. 1989) and used them in solving the same problem. Pail et al. (2001) also used some other numerical techniques to orthonormalized the based function which have global orthogonality support for local gravity field determination. Another idea which is of interest today is the Slepian method (Slepian 1983). In Slepian’s method the ratio of the signal energy in the local and global sense are compared using eignevalue analysis. The eighenvectors maximizing this ratio are selected as the orthogonal base functions with local concentration. Slepian (1983) presented the method for one-dimensional signals but later it was generalized to three-dimensions (Weizerok and Simons 2005) and further developed by Simons et al. (2005).
1.3 Variance components estimation
The estimation of unknown parameters of any least-squares adjustment problem is independent of the choice of a priori variance factor. However if the estimated accuracy of the parameters is desired this variance factor, which plays role of a scaling factor for the variance-covariance matrix of the unknowns, is important. In an ordinary adjustment involved with heterogeneous data estimating one factor for this scale is not relevant and one variance factor should be estimated for each set of observations, which is called variance component. The variance component estimation is a well-known topic in geodetic and related sciences. Different methods exist for computing the variance components. One of the most famous methods is the MInimum Norm Quadratic Unbiased Estimator (MINQUE), presented by Rao (1971). Helmert’s method is another approach presented by Kelm (1978) and by Grafarend and Schaffrin (1979). The method of LaMotte (1973) and Pukelsheim (1981), which is very popular in statistics, was generalized by Schaffrin (1981) for applications in geodesy. Sjöberg (1985) presented a variance component estimator for the adjustment model with a singular covariance matrix. Variance components estimation and applications can also be found in Rao and Kleffe (1988). Maximum likelihood estimation of variance components was first presented for geodetic applications by Kubik (1970), Patterson and Thompson (1971) and (1975), Persson (1980) and Koch (1986). Searle et al. (1992) and Koch (1999) provided useful discussions on the concepts of variance components. As an alternative approach a Monte-Carlo algorithm can be used for variance component estimation (e.g., Kusche 2003) and Fotopoulous (2003) and (2005). In Grafarend (2006) many details of the Gauss-Markov and Gauss-Helmert models are presented, and variance component estimation is treated both theoretically and numerically. Xu et al. (2006) presented a method for computing variance components in linear ill-posed models. They have considered a zero-order Tikhonov regularization method (Tikhonov 1963) for estimating the components. They also found out that the simultaneous estimation of the regularization parameter and variance components is advantageous. Xu et al.
(2007) discussed the estimability of the variance components and proved (as could be expected) that they are not estimable for a fully-unknown variance-covariance matrix. Amiri-Simkooei (2007) and Teunissen and Amiri-Simkooei (2008) explored a new type of least-squares estimator to the components. Eshagh and Sjöberg (2008) presented the modified best quadratic unbiased non-negative estimator of variance components.
Recently a quasi-geoid model was presented over Sweden by Ågren et al. (2009) and named the KTH08 model. This model was computed based on Sjöberg’s theory which was developed during last 24 years. This method is known as the least-squares modification of Stokes’ formula with additive corrections. More details about the method and its implementations can be found in Sjöberg (1984a, 1991 and 2003) and Ellmann (2004 and 2005), Ågren (2004), Kiamehr (2006). The internal error of the KTH08 quasi-geoid model is 19 mm (Ågren et al. 2009). This error was not calibrated and it was estimated based on propagation of the random errors from the gravitational signal and the error degree variance models used in computing the KTH08. The error of the KTH08 model can be calibrated with the global navigation
satellite system (GNSS) data and a variance component estimation process through a combined adjustment of the GNSS heights, normal heights and the KTH08.
Fotopoulos (2003 and 2005) was one of the persons who started using variance components estimation to calibrate the error of geoid models and leveling data. She considered different corrective surfaces to model the discrepancies between geoid models and the geoid computed from the leveling data. Full variance-covariance matrices of geoidal, orthometric and ellipsoidal heights were considered in her analyses. Such matrices can be obtained in different ways; for more details see Fotopoulos (2003 and 2005). Another similar work was carried out by Kiamehr and Eshagh (2008) for calibrating the error of the gravimetric geoid of Iran. They did not use the full variance-covariance matrices of the heights, but they empirically estimated the error of the geoid, orthometric and ellipsoidal heights.
1.4 Validation of satellite gravity gradiometry data and least-squares modification
The geodetic data should be validated prior to use and there is no exception for SGG data either. Therefore we should establish a network over the desired area prior to regionally recovering the gravity field. Different methods of validating SGG data have been proposed. A simple way could be the direct comparison of the real SGG data with the synthesized gravitational gradients using an existing Earth’s gravitational model. Another idea is to use regional gravity data to generate the gradients at satellite level. Haagmans et al. (2002) and Kern and Haagmans (2004) used the extended Stokes and extended Hotine formulas to generate the gravitational gradients using terrestrial gravity data. Denker (2002) used the least-squares spectral combination technique to generate and validate the gravitational gradients. Bouman et al. (2003) has set up a calibration model based on instrument (gradiometer) characteristics to validate the measurements. Mueller et al. (2004) used the terrestrial gravity anomalies to generate the gravitational gradients, and after that Wolf (2007) investigated the deterministic approaches to modify the integrals and validate the SGG data. In fact, the spectral weighting scheme (Sjöberg 1980 and 1981 and Wenzel 1981) was used by Wolf (2007). Stochastic methods of modifying Stokes’
formula, or in other words least-squares modification can be used for the extended Stokes formula as well; see Sjöberg (1984a), (1991) and (2003). Least-squares collocation can be used for validation purposes. Tscherning et al. (2006) considered this method and concluded that the gradients can be predicted with an error of 2-3 mE in the case of an optimal size of the collection area and optimal resolution of data. Zielinski and Petrovskaya (2003) proposed a balloon-borne gradiometer to fly at 20-40 km altitude simultaneously with satellite mission and proposed downward continuation of satellite data and comparing them with balloon-borne data. Bouman and Koop (2003) presented an along-track interpolation method to detect the outliers.
Their idea is to compare the along-tack interpolated gradients with measured gradients. If the interpolation error is small enough the differences should be predicted reasonably by an error model. Pail (2003) proposed a combined adjustment method supporting high quality gravity field information within the well-surveyed test area for continuation of local gravity field upward and validating the SGG data.
Bouman et al. (2004) stated that there are some limitations in generating the
gravitational gradients using terrestrial gravimetry data and geopotential model.
When a geopotential model is used, high degrees and orders should be taken into account and the recent models seem to be able to remove the greater part of the systematic errors. In their regional approach they concluded that the bias of the gradients can accurately be recovered using least-squares collocation. Also, they concluded that the method of validation using high-low satellite-to-satellite tracking data fails unless a higher resolution geopotential model is available. Kern and Haagmans (2004) and Kern et al. (2005) presented an algorithm for detecting the outliers in the SGG data in the time domain.
The second-order radial derivative of extended Stokes’ kernel is isotropic and azimuth-independent. The isotropy is an important property in modifying extended Stokes’ formula otherwise it will not be an easy task. Two methods of generating the SGG data are investigated in Paper J, in the first method takes derivative prior to modification and the second method modifies the formulas prior to taking derivative.
Modification of extended Stokes formula and its second-order order radial derivative based on the biased, unbiased and optimum methods of least-squares modification is the main subject of this study which is a new issue in the scope of SGG. Obviously, these two methods will not deliver the same results, but we are going to test in which cases these methods are comparable. We select the second-order order radial derivative of the extended Stokes formula as its kernel function is isotropic, in such a case, we can use both methods to generate the second-order radial gradient and compare the results. The importance of this study is mostly related to the second method although the first one is new as well. If we can find the cases, where the second method performs well, the method can be used to some how modify the horizontal derivatives of extended Stokes formula having non-isotropic kernels, to generate the other gradients. A similar study was done by Wolf (2007) but just based on deterministic approaches.
1.5 The polar gaps
Polar orbits are suitable for applications of dynamic satellite geodesy because of observing geodetic data with a global coverage. However, the sun-synchronized orbits are preferred for SGG missions. Since the inclination of low orbiters is not equal to 90D, then the polar areas are not covered with satellite observations. The reason is that the gradiometer is a sensitive instrument that one cannot allow movable solar panels. Instead, the solar panels are mounted onto the satellite body and in order to guarantee a constant power supply, the orbit has to be sun-synchronous such that the normal vector of the orbital plane points to the sun (Sneeuw 2009, personal communication).
The satellite orbits used in practice for the gravity field recovery missions are inclined causing the polar gaps in the data coverage. If the inclination of low orbiters is not equal to 90D, then the polar areas are not covered with satellite observations.
The polar gaps problem was described by Rummel et al. (1993) and Koop (1993) with concentration on determining the spherical harmonic coefficients using the least-squares method. Sneeuw and van Gelderen (1997) also studied this problem and explained why the near-zonal harmonics are weakly determined in geopotential
modeling. They have considered both space-wise and time-wise approaches in their discussions. Tscherning et al. (2000) investigated the polar gaps problem using least- squares collocation. They concluded that the inclusion of gravity data in the gaps will improve the estimates of the geopotential coefficients. Tscherning (2001) mentioned that the near-zonal harmonics were typically two times larger than the other harmonics if the poles were not covered. He found that if ground gravity data is used to fill-in the gaps, the data should have a resolution twice the second-order radial derivative T . He concluded that the available airborne and sub-marine gravimetry zz data would not improve the gravity field and steady-state ocean circulation explorer (GOCE) solution if data from CHAMP (Challenging Minisatellite Payload) (Reigber et al. 1999 and 2004), and GRACE (Gravity Recovery and Climate Experiments) (Tapley et al. 2005) are available. Pail et al. (2001) used an orthonormalization scheme to orthonomalize the base functions with reference to Hwang (1991) and (1993). They recovered the geopotential coefficients without filling-in the gaps.
Rudolph et al. (2002) considered the polar gaps filling-in problem from the satellite gradiometric data, they conclude that 1 1D× D gravity anomaly blocks can be recovered. They also mentioned that the error of the recovered gravity anomalies reduced from 3 to 10-12 mGal from the rim to the centre of gap. They found also the cross-track direction contain the most information. Metzler and Pail (2005) presented a spherical cap regularization method considering an analytical function, which is defined in the polar region using a geopotential model. Simons and Dahlen (2006) used the spherical Slepian function for the polar gaps, and they proposed a new method that expands the source field in terms of a truncated basis set of spherical Slepian function and compared it a with least-squares method. Siemes et al. (2007) incorporated three strategies: augmenting data in the gaps, introducing a priori information for regularization and the Slepian approach again.
1.6 Special studies on geoid and Moho
Most methods of geoid determination combine different types of gravimetric data.
As it is well–known in physical geodesy a spherical harmonic expansion can be considered for expressing the geoid height, but the problem is that this series of expansion will be restricted to a certain frequency. This is why that Stokes integral is combined with a spherical harmonic expansion of the geoid. However, these approaches need global coverage of gravimetric data with very high resolution to consider the higher frequency portion of the geoid which is not practical.
Molodensky et al. (1962) presented a method to modify Stokes formula in such way that minimizes the effect of far zones data, and the integration is limited to a small cap around the computation point. Following his idea different geodesists worked on modification of Stokes formula; see e.g. Wong and Gore (1969), Meissl (1971), Jekeli (1980 and 1981), Sjöberg (1984), Heck and Grueninger (1987), Vanicek and Kleusberg (1987) and Featherstone et al. (1998).
Sjöberg (1984) proposed a method based on spectral combination (Sjöberg 1980 and 1981, Wenzel 1981 and 1983) to minimize error of terrestrial and global data as well as the truncation error in a least square sense. This method is well-known as least-squares modification. Sjöberg (1991) presented refined least-squares modification and discussed about biasedness and unbiasedness of his geoid
estimators. Sjöberg (2003) also proposed a general model of modifying Stokes formula in a least-squares sense and presented an optimum geoid estimator that the bias, unbiased geoid estimators are especial cases of that estimator. Ellman (2004) and Ågren (2004) both worked on Sjöberg’s geoid estimators practically to determine the geoid of Baltic countries and Sweden, respectively. Ellman (2005) compared the biased, unbiased and optimum way of modifying Stokes formulas and concluded that the results of these three types of modification are very similar.
Kiamehr (2006) successfully used the unbiased type modification to determine a geoid for Iran.
Nowadays airborne gravity mapping is performed more or less as a standard survey. The primary data are semi-continuous gravity data along flight paths, and various methods exist in converting this information to geoid heights (e.g., Forsberg and Kenyon, 1994, Tscherning et al. 1997, Forsberg et al., 2002 and Jekeli, 2001, Ch.10). A primary idea in aerial gravimetry is to use a two step procedure to compute the geoid, where in the first step the airborne gravity data is analytically continued downward to sea level, and in the second step Stokes’ or Hotine’s integral formula is used to convert these data to a geoid model. Numerous scientists have worked on related problems, especially on the problem of downward continuation of the data.
For example, Jekeli (1987) presented a Molodensky type of technique, which first solves for a surface density layer on the Earth’s surface, and, once the surface layer has been solved for, it can be used to determine any gravity related function on or above the Earth’s surface. Tscherning et al. (1997) used least squares collocation.
Also, as a more efficient approach than the two-step procedures, Novak (2003) proposed to directly solve for the geoid by a one-step Fredholm integral equation of the first kind. Tenzer and Novak (2008) compared these types of solutions further, and once again they concluded that the one-step solution is more stable and efficient.
Alberts and Klees (2004) numerically compared several different approaches for the inversion of airborne gravity data, and they concluded that least squares collocation performs slightly better than the integral based methods. Forsberg et al. (2000) and Forsberg (2002) reported on successful aerial gravity surveys over Skagerrak, Greenland and Svalbard. They estimated the accuracy of the resulting geoid solution for the former survey to the order of 10 cm. For the latter survey they used collocation and fast Fourier transform to downward continue airborne data, and they concluded that the merging of high altitude airborne gravity with surface data with the joint utilization of digital terrain models seems essential for obtaining reasonable geoid models in very rugged topography regions. Hwang et al. (2006) also used the fast Fourier transform in downward continuing scalar airborne gravimetry data in Taiwan. Serpas and Jekeli (2005) proposed a method for local geoid determination from airborne vector gravimetry. They concluded that the horizontal components of the gravity vector improved the quality of the geoid, which can be determined to accuracy better than 10 cm.
As the above literature review shows, inversion of integral formulas or least squares collocation play a key role for downward continuation of airborne gravimetry. In contrast, in this study we propose solving for the geoid by using a two-step approach: first computing the disturbing potential at flight level by a combination of an Earth gravity model and Hotine’s integral formula and then the disturbing potential is continued analytically downward to geoid level by a Taylor
expansion. Topographic corrections are applied separately. This scheme implies the following computational formula for the determination of the geoid height. The details of this study will be presented in Paper P.
The discontinuity between the Earth’s mantle and crust is called the Mohorovičič discontinuity (Moho). This surface is of great importance in Geophysics to characterize the overall structure of the crust and it can be related to geologic and tectonic evaluations. Having a good knowledge about the crustal thickness is essential in heat flow determination and seismology. There are two well-known methods for determining this surface: a) the seismic method and b) the gravimetric method. In Seismology, Moho is defined as a surface where the velocity of the seismic wave changes, and we call this surface the seismic Moho. In the gravimetric method and gravimetric data is used based on an isostatic hypothesis, which is based on isostatic equilibrium of the crust on the mantle. We name this model of Moho the gravimetric Moho. Three main hypotheses in this respect are due to Pratt-Hayford (Pratt 1855), Airy-Heiskanen (Airy 1855) and Vening-Meinesz (Heiskanen and Vening Meinesz 1958 p. 131). The last hypothesis agrees better with reality than the others, as it states that the isostatic compensation is regional rather than local. The seismic and gravimetric methods of Moho determination will not yield the same surface, as they are theoretically and practically different. Consequently, the combination of the results of these methods should be a better solution in statistical point of view and hopefully more realistic than the individual ones. The seismic method is based on the travel-time observation of the reflected seismic signal. The advantages of using artificial explosions rather than earthquakes are that the time and position of the shot are accurately known. There are two limitations in the seismic surveys to estimate of the Moho depth. First, the seismic method suffers by lack of global coverage of seismic data, while according to the latest developments in satellite technology and gravimetry, the problem of data gap was solved in the gravimetric method. However, the latter method suffers by the isostatic hypothesis, which is used in formulating the Moho model, because it is not clear how realistic the applied hypothesis is, and we do not know if the real Moho follows this hypothesis. Second, the seismic data acquisition is expensive and the amount of data collected in a survey can rapidly become overwhelming. Therefore, data reduction and processing can be time consuming, requiring sophisticated computer hardware, and demand the considerable expertise.
Vening-Meinesz (1931) modified the Airy-Heiskanen theory to estimate the crustal thickness by considering a regional instead of local compensation. The Vening-Meinesz hypothesis assumes a constant density contrast and a variable Moho depth in defining the Moho. Parker (1972) presented a practical gravimetric method based on a constant density contrast and a varying Moho depth, similar to Vening- Meinesz’s hypothesis, in the Fourier domain. Oldenburg (1974) used the method of Parker (1972) and presented a method to compute the topographic density contrast from the inversion of gravity anomaly in a two-dimensional frame. Parker- Oldenburg’s method was generalized by Gomes-Oritz and Agarward (2005) and Shin et al. (2006). Kiamehr and Gomes-Ortiz (2009) estimated the Moho depth in Iran based on terrestrial gravimetric data and the Earth gravity model EGM08 (Pavlis et al. 2006) by Parker-Oldenburg’s method. Čadak and Martinec (1991) presented a
model for Moho in terms of spherical harmonics to degree and order 30. Martinec (1992) and (1994) investigated the gravitational effect of the Moho discontinuity to determine the density contrast between the mantle and crust by minimizing the sum of the gravitational potentials induced by the Earth’s topographic masses and the discontinuity of the Moho.
Moritz (1990) improved the Vening-Meinesz hypothesis by developing it to a spherical Earth model. Sjöberg (2009b) presented a new solution for the theory of Moritz, named the Vening Meinesz-Moritz (VMM) inverse problem in isostasy. He presented some iterative as well as approximate spherical harmonic solutions. The gravimetric Moho model, which is discussed here, is based on the spherical harmonic solution. Bagherbandi and Sjöberg (2009) investigated the regional and local compensations in the inverse isostatic problem. They compared the traditional isostatic hypothesis i.e. the Airy-Heiskanen Moho model and the gravimetric Moho model (based on the Vening-Meinesz hypothesis) in some geophysically critical areas. They found that the Vening-Meinesz solution has a good agreement with the simple mean, but there are notable differences in some special areas, e.g.
Fennoscandia and Canada due to glacial isostasy.
Sjöberg (1980) and (1981) and Wenzel (1981) proposed a method of combining integral formulas in the spectral domain. Sjöberg (1984) develop this technique to the spectral form of Stokes integral and presented the least-squares modification of Stokes’ formula. Bjerhammar (1983) suggested using variance component estimation to mix discrete boundary value problems, and he used this technique to combine solutions of Fredholm integral equations of the first kind.
2.The author’s contributions
1. The topographic and atmospheric biases can be considered for SGG data when the gravity gradients are directly continued downward to sea level for local gravity field determination. It is promising to use the gravity field and steady-state ocean circulation explorer (GOCE) data for this goal, namely direct analytical downward continuation of the GOCE data and removing the topographic and atmospheric biases from the downward continued data. In Paper A we review the topographic bias and present the atmospheric biases in spherical harmonics which is a new subject. After that the biases are considered in gravity gradiometry and the mathematical models of the biases are presented for each gravitational gradient which were not investigated until now. The convergence of topographic and atmospheric biases on the gradients is numerically studied. It was shown the spherical harmonics of expansions of the topographic and atmospheric biases are convergent to second-order topography terms while they are asymptotically convergent to third term. This means that as long as the gradient are measured at satellite level we have not problem in convergence of the harmonic series, but when the altitude reduced higher degrees of the effect are being important and the contribution of the third term increases as well.
2. The above studies miss an investigation of the contribution of the topographic terms obtained based of a binomial expansion of the topographic heights in formulating either the topographic or atmospheric effect on the SGG data. Paper B will investigate this matter. Some studies has been done in convergence of the binomial expansion of the topographic heights by Sun and Sjöberg (2001) and they concluded that truncation of the binomial expansion to the third order is a good approximation as along as the degree of the spherical harmonic expansion is not so high. Here we take advantage of the real topographic data and we study the contribution of each term of the binomial expansion in the topographic and atmospheric effects on the SGG data at 250 km level, which is a new study in the scope of SGG.
This investigation is presented in Paper B. Also the effect of each topographic term of the spherical harmonic expansions of the topographic and atmospheric effect was also investigated. It was concluded that the atmospheric effect to linear approximation of the binomial expansion of the height is enough for SGG data. The topographic effect to second-order topographic term is meaningful for the SGG data.
3. Paper C used Sjöberg’s method of considering the lateral density variations of the topographic masses for the topography and the crustal layers. Paper D differs from the reviewed articles above, as we are going to compare two different methods of considering the laterally varying density of the topographic masses in formulation of their effect on the SGG data. Since the topographic features are the upper part of the upper crust layer we can assume that the topographic masses have the same lateral density variations as the upper crust. For this goal we use CRUST2.0 (Bassin et al. 2000) which is the updated version of CRUST5.0 (Mooney et al. 1998). This paper will investigate the Sjöberg method and Novák and Grafarend’s method theoretically and numerically which is also a new study in the scope of SGG.
It should be mentioned, that the method of Sjöberg was originally presented for geoid determination goals and the method presented by Novak and Grafarend (2006) was only an idea which was not developed. In this part of the studies, we developed the method of Sjöberg to satellite gradiometric data and develop the idea of Novak and Grafarend (2006) and used it in practice; for more details please see Paper D.
4. Reed (1973), Xu (1992, 1998 and 2009), Kotsakis (2007) and Janak (2009) used the inversion of the extended Stokes integral formulas for recovering the gravity anomaly at sea level. Xu (1992, 1998 and 2009) and Kotsakis (2007) considered the second-order radial derivatives of the gravitational potential. The works of Reed (1973) and Janak et al. (2009) are similar to those presented in Paper E, but they used a geocentric spherical frame in their simulation, while we use the local north-oriented frame. Janak et al.
(2009) took advantage of some existing methods of regularization without using a priori information. No effort was done for estimating the error of the recovered gravity anomalies and removing the bias of regularization as well as Joint inversion of the gradients. The use of variance component
estimation in linear ill-posed problems is also a new subject which is presented in Paper E.
5. Until now most of the studies are related to a spatially restricted integral product of scalar spherical harmonics. However, having such a product is not sufficient for localizing the gradiometric boundary value problem. The spherical Slepian functions in gradiometric boundary value problems involve the integral products of the vector and tensor spherical harmonics. These integral products are spatially restricted and complicated. Paper F will simplify these integrals in terms of integral of associated Legendre functions and combinations of the gaunt coefficients.
6. As the above literature review shows, most of the scientists used the inverse solution of integral formulas, least-squares solution or least-squares collocation for regional gravity field determination. The method which was not investigated is related to direct integral formulas. Tscherning et al.
(1990) proposed such a method in regional gravity field modelling. They just considered the second-order radial derivative of disturbing potential to express the idea. Their idea was not developed further until now. Paper G tries to use and develop this idea to the solutions of gradiometric boundary value problems. The beauty of the integral approach is its direct solution, which means that we do not invert integral formulas and the satellite data are continued downward and integrated simultaneously to estimate the disturbing potential at sea level.
7. As the above review shows most of the variance component estimation problems are outlined for the well-posed adjustment problems. Only Xu et al. (2006) has investigated it in the ill-posed problems which are solved using the Tikhonov stabilizer. They showed how to estimate the bias of the variance components due to the regularization and proposed a biased- corrected estimator of variance components. The main idea of Paper H is to theoretically study the variance component estimation in the ill-posed problems which are solved by the truncated singular value decomposition which is one of the famous regularization methods. The bias and the biased- corrected estimator are presented. A non-negative biased-corrected estimator is developed, which is totally new. This estimator was studied by Eshagh (2009a) for the Tikhonov regulator (Tikhonov 1963) but not for the truncated singular value decomposition in Paper H. Some studies on this subject were done in Paper E, in local gravity field determination based on direct inversion of the SGG data. The ill-conditioned system of equations related to the inversions was regularized using Tikhonov regularization.
Later on, we developed this technique further and presented a method of estimating the variance components in a linear ill-posed problem which are stabilized by truncated singular value decomposition for more details see Paper H. In this paper we are going to perform similar error calibration for the KTH08 model. However, this study and its involved problems are somehow different with those performed by others. First of all, we have not full covariance matrices of the heights and the quasi-geoid model. Second,
the errors of the quasi-geoid model as well as the normal heights are constant for all the points. This yields singularity in the system of equations from which the variance components are estimated. Also the error of the zero-, first- and the second-order GNSS height networks are available.
Although the theory of variance component estimation and error calibration of geoid is not new but more practical considerations are needed to do similar job for the Swedish data. In this respect different stochastic models will be defined and used in the variance component estimation process; see Paper I for more details.
8. However we concentrate on generating the second-order radial gradient based on the least-squares modification approaches; for more details see Paper J. The second-order radial gradient can also be generated using the Abel-Poisson integral. This matter is investigated in Paper K of the present report. Both studies have shows than the quality of the gravity field model which is used plays an important role in validation. Some studies have been done based on the EGM08 (Pavlis et al. 2006) geopotential model and the traditional Tscherning/Rapp (1974) model. The conclusion of this study is that the spectra of the gravity do not improve the least-squares modification while the error spectra of EGM08 are very important to obtain a qualified solution. Such an investigation is given in Paper L.
9. In Paper M, it is investigated which parts of the gravitational signal are not well estimated in the spectral solutions of gradiometric boundary value problems. Similar investigation was done by Sneeuw and van Gelderen (1997) only for the vertical-vertical solution, which is simplest one among the other solutions of gradiometric boundary value problems. The other solutions have not been treated so far. Here, we want to present the polar gaps effects on these solutions. In Paper N we consider full tensors of gravitation at satellite level, in other words, all types of SGG data are tested and investigated for filling-in the polar gaps by gravity anomaly. Such an investigation was not already done. Another new issue is the joint inversion of SGG data for gravity anomaly recovery in the gaps. Propagation of random errors and removing the biases of regularization method are two other new issues which will be investigated in Paper N.
10. All the efforts were done by the geodesists was to modify Stokes formula based on data inside the cap around the computation point. However, in Paper O an attempt is made to modify Stokes formulas in a reverse situation based on the data outside the cap or inside a ring (between two concentric caps) around the computation point. This type of modification is useful to approximate the geoid for the regions with data gaps which is a considerable problem of most developing countries.
11. We also present an airborne solution for geoid. In this method the gravity disturbance are integrated based on the modified Hotine’s formula at airborne level to obtain the disturbing potential. After that the potential in continued downward to see level; for more details see Paper P
12. In Paper Q we first estimate an optimal value for each spectrum of the Moho model. We construct condition adjustment models and take advantage of the variance component estimation process to get a statistically optimal spectral combination of the Moho models. Eshagh (2009a) and (2009c) successfully developed and used this method to combine the integral solutions of the gradiometric boundary value problem in the spectral domain. He defined the degree-order variance components and used the variance components to update the weight of observations. Here this type of estimation is used in combining the seismic and gravimetric Moho models to estimate a combined degree-order variance component estimation of Moho model after transferring them into the spectral domain. The simple mean and the weighted mean in the spectral domain are also considered and compared with the combined model.
3. Obtained results and conclusions
The harmonics of the atmospheric bias of V is convergent just to the first zz topographic term. It means that for computing the atmospheric bias the harmonics should be generated to this term independent on the altitude of gradiometer and application. There is no convergence problem on Vxy and it can be continued downward to sea level analytically and the biases can be removed from the results;
and it is not important at which altitude the gravity gradients are measured. The harmonics of the atmospheric bias of the other gradients are convergent to the second order topographic term and they are asymptotically convergent to the third term. It means that the harmonics should be generated to the second term when the gravity gradients are measured near the Earth’s surface but for the SGG missions, e.g., GOCE the third term can also be used as the GOCE data cannot sense higher frequencies of the Earth gravity field and the maximum degree of the gravity field is smaller than the convergence degree. The numerical results show that the contribution of the second and third terms of the topographic effects reaches to 80 mE and 1 mE, respectively. The effect of these terms for the atmospheric effect is below 1 mE for the exponential, Sjöberg models and Eshagh/Sjöberg’s models and may be negligible. Among the gradients, V shows more sensitivity with respect to zz the topographic terms than the others. The topographic effects derived based on the methods of Sjöberg and Novák and Grafarend differ less than 1 E. The differences are small in comparing them with the topographic effects (less than 8%). The reason of such a difference is due to the interpolation error of the geometric and density functions in the Novák and Grafarend method. However these methods are equivalent theoretically. Therefore it is suggested to see which method is fitted to the original data to avoid further approximations in practice.
For local gravity field determination Tzz is the best gradient and also Txx, Tyy and T are more suited than the rest. The systems of equations related with xz T and zz
Tyz filter out the random noise more efficiently than the others. The propagated errors of T , xx T , yy T and zz T are decreasing from north to south and increasing for xz T yz and Txy. The numerical study showed that the removal of the regularization bias from the estimated variance of unit weight will not improve the propagated errors.
By removing the bias from the recovered gravity anomaly (the bias-correction process) one can reduce the error by 3 mGal. The recovery of the gravity anomaly with an error of 1 mGal by inverting a single gradient is not easy, unless the resolution of the grid of gradients increases considerably. The study shows that the gravity anomaly can be recovered with 2 and 5 mGal errors at 1 mE and 10 mE noise levels, respectively. The bias-correction procedure reduces these errors to 1.6 and 3 mGal, respectively. It is expected to reach to 1 mGal accuracy, if the resolution of the SGG data increases.
In order to locally concentrate the solution of the gradiometric boundary value problem, the spherical Slepian function can be used to maximize the signal energy in the desired area. The eigenvalues yielding the maximum ratio of the local and global signal is selected and their corresponding eigen functions are selected as the based functions. In the numerator of these ratios there exist spatial restricted integral products of scalar, vector or tensor spherical harmonics which have complicated forms. However this paper presented and simplified these integral products in terms of combinations of the gaunt coefficients and integral of associated Legendre functions. The presented formulas for the spatially restricted integrals of these spherical harmonics can be used in recovering the Earth’s gravity field locally from SGG data.
In the numerical studies on the integral approach to regional gravity field modeling we found the degrees between 150 and 250 based on 1 and 5 mE Gaussian noise. The kernels of the integral are not well-behave and even after removing the long wavelength structure of the field, the contribution of fat zone data remains relatively large. Here we suggest band-limiting the kernels and removing the truncation bias of the integral formulas. The truncation bias of the integral related to T decreases well, but this is not the case for the other integrals. In global gravity zz
field modelling the truncating the spherical harmonics expansion to the degree of resolution is very similar to that case where the Weiner filter is used (according to corresponding noise levels). However in regional gravity field modelling, it is not true. The numerical study on the filtering showed that this process can improve the results about 50%. In regional gravity field modelling the integrals involving T is zz suited and after that Txz yz, and Txx yy− ,2xy. Numerical studies showed that by integrating T with 1 mE noise and in a cap size of 7zz D, the geoid can be recovered with an error of 12 cm after the filtering process. Similarly, the errors of the recovered geoids from Txz yz, and Txx yy− ,2xy are 13 and 21 cm, respectively.
The variance components can be estimated in the ill-posed problems which are solved based on the truncated singular values decomposition. In this study we presented an unbiased estimator of the a posteriori variance factor and analyzed the
