Comparison of remove-compute-restore and least squares modification of Stokes' formula techniques to quasi-geoid determination over the Auvergne test area

Comparison of remove-compute-restore and least squares modification of Stokes' formula techniques to quasi-geoid determination over the Auvergne test area The remove-compute-restore (RCR) technique for regional geoid determination implies that both topography and low-degree global geopotential model signals are removed before computation and restored after Stokes' integration or Least Squares Collocation (LSC) solution. The Least Squares Modification of Stokes' Formula (LSMS) technique not requiring gravity reductions is implemented here with a Residual Terrain Modelling based interpolation of gravity data. The 2-D Spherical Fast Fourier Transform (FFT) and the LSC methods applying the RCR technique and the LSMS method are tested over the Auvergne test area. All methods showed a reasonable agreement with GPS-levelling data, in the order of a 3-3.5 cm in the central region having relatively smooth topography, which is consistent with the accuracies of GPS and levelling. When a 1-parameter fit is used, the FFT method using kernel modification performs best with 3.0 cm r.m.s difference with GPS-levelling while the LSMS method gives the best agreement with GPS-levelling with 2.4 cm r.m.s after a 4-parameter fit is used. However, the quasi-geoid models derived using two techniques differed from each other up to 33 cm in the high mountains near the Alps. Comparison of quasi-geoid models with EGM2008 showed that the LSMS method agreed best in term of r.m.s.


Introduction
Nowadays, a major goal for physical geodesy is the determination of the geoid with an accuracy on the cm level, matching the accuracy of GPS height determination. Although the cm-geoid * E-mail: hasan.yildiz@hgk.msb.gov.tr, Phone: +90 312 5952216 The manuscript solely reflects the personal views of the author and does not necessarily represent the views, positions, strategies or opinions of Turkish Armed Forces has been demonstrated in numerous cases in lowland areas with dense gravity data coverage, so far no convincing case of attaining a cm-geoid in mountainous regions has been reported. This is likely a consequence of insufficient gravity data coverage, theoretical shortcomings and insufficient quality of the levelling data, used to compute "ground truth" geoid (or quasi-geoid) values.
The most commonly adopted and applied approach to regional gravimetric geoid determination today is probably the removecompute-restore (RCR) technique (e.g., Schwarz et al., 1990). In the remove step, a long-wavelength part (predicted by a global Brought to you by | DTU -Technical Information Center of Denmark (DTIC) Authenticated | 192.38.67.112 Download Date | 3/4/13 3:23 PM gravity field model) and a short-wavelength part (predicted by topography) are removed from the original gravity data. In the compute step, the obtained band-pass filtered gravity anomalies are transformed into quasi-geoid heights either using Stokes' integration methods, e.g. 2-D spherical Fast Fourier Transform (FFT), or using Least Squares Collocation (LSC). After having carried out the compute step, the long-wavelength and the short-wavelength parts are restored to the quasi-geoid. One advantage of the LSC method is that it does not require gridded gravity anomalies as the FFT method and provides error estimates for the resulting quasi-geoid models.
Alternatively, least squares modification of Stokes' formula (LSMS), developed at the Royal Institute of Technology, does not require gravity reductions and includes a least squares kernel modification with additive corrections for the topography, downward continuation, the atmosphere and the ellipsoidal shape of the Earth (Sjöberg 2003). A detailed theoretical discussion about RCR and LSMS techniques can be found in Sjöberg (2005).
This study aims to numerically compare the RCR and LSMS techniques in the Auvergne test area. A number of studies have already reported results using different quasi-geoid determination methods in the same test area (Barzaghi and Sanso, 2009;Ågren et al., 2009b;Forsberg, 2010). However, this study differs from the previous studies because identical input data have been used for all methods, and the quasi-geoid models are not only compared with the GPS-levelling data, but also among each other and with the EGM2008 (Pavlis et al., 2008). Geodetic Gravity Field Modelling Programs (GRAVSOFT) (Forsberg and Tscherning, 2008) and Scientific Software for Precise Geoid Determination Based on the Least-Squares Modification of Stokes' Formula (LSMS-GEOLAB) (Kiamehr and Sjöberg, 2010) are used for the practical implementations of the RCR and LSMS techniques, respectively.
Section 2 briefly describes the test area, data and the quasigeoid determination methods. Section 3 gives the details of the computation of regional quasi-geoid models and comparison with GPS-levelling data. Quasi-geoid difference analyses are given in Section 4. Finally, the results are discussed and some conclusions are presented in Section 5.

Data
The Auvergne data set was distributed in 2008 by Institut Géographique National (IGN), France, on behalf of the International Geoid Service (IGeS), as a "ground truth" example for precise geoid determination methods (Duquenne, 2007).  Figure 1. a) The data coverage for quasi-geoid determination (Duquenne, 2007); b) Distribution of GPS-levelling points (black circles) in the central region with the heights (in metre) of quasi-geoid computational points selected at 0.02 • ×0.025 • resolution (approximately at 2 km spacing) from 3" Shuttle Radar Topography Mission (SRTM) height data over the quasi-geoid comparison area levelling connections (see Fig 1(b)), and a quoted GPS ellipsoidal height accuracy of 2−3 cm (RBF points) or "slightly better" (NIVAG points) (Duquenne, 2007). Duquenne (2007) suggested that these GPS-levelling points were linked to the French national levelling network (NGF-IGN69) by precise levelling with redundant observations, with the total standard deviation of the difference in heights between points in the GPS-levelling area better than 2 cm, including the uncertainties of the basic network and of the local ties. The NGF-IGN69 uses normal height system and is tied to Marseille tide gauge. Considering the suggested accuracy of GPS ellipsoidal heights as 2.5 cm and the accuracy of precise levelling measurements as 2 cm, an accuracy of 3.2 cm can be attained for the quasi-geoid heights of the GPS-levelling points through the error propogation law assuming uncorrelated observations.
The elevations of the GPS-levelling points range from 206 to 1235 m, and the highest mountain in the central area is 1886 m; it is therefore a relatively moderate mountainous area (Fig 1(a), 1(b)).

Journal of Geodetic Science 55
Shuttle Radar Topography Mission (SRTM) height data over the quasi-geoid comparison area are shown in Fig 1(b)). The elevations of the quasi-geoid computational points range from 21 m to 1823 m whereas the elevations of the free-air gravity anomalies range from 19 m to 1715 m over the quasi-geoid comparison area (Fig 1(b)).

Remove Compute Restore Technique
In the RCR technique, the anomalous potential T is split into three parts: where T EGM is the contribution of an Earth Geopotential Model (EGM). T RT M are the terrain effects from Residual Terrain Modelling (RTM), and T RES the residual gravity field. T is treated as a spatial function. One reason for subtracting an EGM is to represent the gravity field outside the area covered with data.
For the terrain-reduction, the terrain effects are reduced relative to a mean elevation surface. The terrain potential is subtracted from the observations using a prism integration, i.e. representing the mass between the actual topography and the mean elevation surface as mass prisms of either positive or negative density, nominally 2670 kg/m 3 . The prism implementation of the RTM method has an inherent problem: the method leaves a point above the mean elevation surface in the mass-free domain, whereas a point below the mean elevation surface after the reduction would correspond to the value inside the reference topography mass. As all geodetic gravity field modeling methods require observations derived from a harmonic function, i.e. in a mass-free environment above the geoid, the harmonic correction is applied to the gravity anomaly points below the mean elevation surface; for details see Forsberg (1984). The geoid "restore" signal is computed by Fourier methods (for details see Forsberg, 1984Forsberg, , 1985. In RTM method the resolution of the mean elevation surface is controlled by the user, through a suitable low-pass filter. Ideally such a filter should correspond to the equivalent resolution of the highest spherical harmonic expansion degree in the reference potential.
The Stokes/Molodensky' formula (Moritz, 1980, p. 387 where ζ is quasi-geoid, R is the mean Earth radius, γ is the normal gravity at the normal height, σ is the unit sphere, S(ψ) is the Stokes' function with argument ψ as the geocentric angle, ∆ is the surface gravity anomalies, 1 is the first term in the Molodensky expansion. Equation (2) is modified using reduced surface gravity anomalies (∆ ), obtained removing the contribution of an EGM and terrain effects from the surface gravity anomalies (∆ ), instead of ∆ to obtain the residual quasi-geoid contribution (ζ ) by the formula, When the RTM reduction is used, the Molodensky 1 term will generally be insignificant (Forsberg and Sideris, 1989;Schwarz et al., 1990) and the formula converts into the conventional Stokes' formula for the geoid, in principle applied to gravity field quantities at the geoid. The evaluation of Stokes' formula is done using 2-D spherical FFT for details of the methods see Schwarz et al. (1990) and Forsberg and Sideris (1993). In the 2-D spherical FFT method, Stokes' integral is evaluated over the whole gravity anomaly area by a series of convolutions, each accurate around a certain reference latitude. To keep the inherently highly accurate Gravity Recovery And Climate Experiment (GRACE) gravity field information in EGM2008 (Pavlis et al., 2008) from being overruled by the influence from terrestrial gravity data, the integral should use modified Stokes' kernels, e.g., as the modified Wong-Gore kernel (Wong and Gore, 1969) described by Forsberg and Olesen (2010).
The LSC method also uses the RCR technique (Forsberg and Tscherning, 1981), as the data to be used for covariance function estimation and the subsequent LSC step are required to be smooth with small variance in order to properly interpret the error estimates. The LSC method can handle heterogeneous observations to estimate gravity field components and their standard errors, such as geoid heights (Tscherning, 1982). The LSC method takes into account data located at different altitudes through the use of a spatial covariance function. In the LSC method, a limitation, however, has been that as many equations as the number of data have to be solved. Therefore, we used a limited number of gravity data selected at approximately 2 km resolution for the impementation of the LSC method concerning computational capabilities, which are also used by the 2-D FFT and LSMS methods as identical data. First, emprical covariances are computed and subsequently these values are fitted to the model covariance function of Tscherning-Rapp (Tscherning and Rapp, 1974). The residual quasi-geoid undulations are determined by the LSC method, where the required auto-and cross-covariance functions are computed by covariance propagation from the analytically modelled local covariance function represented as follows: from the geocenter, R B is the radius of Bjerhammar sphere and σ 2 is the error degree variance. The covariance parameters (scale parameter), A (a constant parameter in units of (m/sn) 4 ) and the Bjerhammar radius R B are determined using an iterative non-linear adjustment (Knudsen, 1987).

Least Squares Modification of Stokes' Formula (LSMS)
The least squares modification of Stokes' formula (LSMS) (Sjöberg, 2003) to compute a gravimetric model, includes a least squares (stochastic) kernel modification with additive corrections for the topography, downward continuation, the atmosphere and the ellipsoidal shape of the Earth. The LSMS method has been presented under several different versions during the years (see e.g. Sjöberg, 2003). Ågren et al. (2009a) clearly describes the LSMS method for the estimation of geoid heights, as well as the transformation of the method to estimate height anomalies. As we aim to determine height anomalies, we used the version applied by Ågren et al. (2009a) by using the combined estimator for the height anomaly: where σ 0 is the spherical cap, R is the mean Earth radius, γ is mean normal gravity, S M (ψ) is the modified Stokes' function, are the modification parameters, M is the maximum degree of the EGM, Q M are the Molodensky truncation coefficients and ∆ EGM is the Laplace surface harmonic of the gravity anomaly determined by the EGM of degree . The four additive corrections are shown to the right in equation (5), where the combined topographic effect δζ C OMB = 0 (Sjöberg, 2000) and δζ DW C , δζ C OMB AT M and δζ ELL represent the downward continuation effect, atmospheric and ellipsoidal corrections, respectively. See Ågren et al. (2009a) for the formulas of the atmospheric and ellipsoidal corrections.
The downward continuation effect δζ DW C is where P is the computation point, P = R + H P ζ 0 P is an approximate value of the quasi-geoid height and Q is the running point in Stokes' integral, σ 0 is the spherical cap, R is the mean Earth radius, γ is mean normal gravity, S M (ψ) is the modified Stokes' function, are the modification parameters, M is the maximum degree of the EGM, Q M are the Molodensky truncation coefficients and ∆ EGM is the Laplace surface harmonic of the gravity anomaly determined by the EGM of degree (Ågren et al., 2009a). Ågren et al. (2009a) points out that Equations (5) and (6) are equivalent to analytical continuation to point level using the 1 term in Moritz (1980, p. 387) except that they differ from Moritz' in that least squares modification of Stokes' formula is utilised with improved atmospheric and ellipsoidal corrections.
One problem with using the combined quasi-geoid estimator in Eq. (5) is that Stokes' quadrature is made on the rough surface gravity anomaly, which results in large discretisation errors (Ågren et al., 2009a). However, by taking advantage of the removecompute-restore philosophy for the gridding of a comparatively dense gravity anomaly grid using a smoothing topographic correction, such errors can be counteracted (Sjöberg, 2003). This makes it possible to take advantage of the high-frequency information available in the Digital Elevation Model (DEM). However, a practical drawback here is that dense gravity anomaly grids are required in rough mountain areas (Ågren et al., 2009a). One advantage of the LSMS method is that the "real" magnitude of the corrections becomes apparent, i.e. the sum of the additive corrections is equal to the quasi-geoid errors obtained in case no corrections are applied at all (Ågren et al., 2009a).
Choosing a suitable EGM and modification parameters are essential steps in the determination of a quasi-geoid model using the LSMS formula. We need to estimate the signal and error degree variances of the EGM and terrestrial gravity to be able to estimate the least-squares modification coefficients. The signal degree variances for the degrees above the EGM are generated using the Tscherning and Rapp (1974)

Computation of regional quasi-geoid models
As previously mentioned, concerning the computational capabilities for the implementation of the LSC method, A limited number of 59097 gravity points are randomly selected at 0.02 • ×0.025 • resolution (pixel binning to approximately 2 km resolution) (Fig. 2), to be used as identical input data for all three quasi-geoid determination methods. This computational strategy can be critized, as in reality only LSC needs to be limited by the number of observations.

Regional Quasi-Geoid Model Determination Using the RCR Technique
The EGM2008 is computed to degree 360. The terrain effects are computed using 3" SRTM data with respect to a mean eleva-  The 2-D spherical FFT conversion was done using zero-padding in a zero-padded grid of dimension 600×640 data points, using 3 reference latitude parallels. The RTM geoid terrain effects are subsequently restored using height anomaly terrain effects computed by FFT in a similar grid, and finally the gravimetric quasi-geoid was obtained by adding the geoid EGM2008 effects, computed at the level of the topographic surface. We first computed the quasi-geoid without using any modification of Stokes' formula.
In this case, we found a standard deviation of 3.9 cm between the 2-D FFT derived quasi-geoid model and GPS-levelling data using a 1-parameter fit (Table 3). Subsequently, the modification parameters are investigated in comparison with the GPS-levelling data which results in a selection of an optimal degree band of 80-90, in good agreement with the estimated accuracy of GRACE (Tapley et al., 2005). In this case, a standard deviation of 3.0 cm is found between estimated gravimetric quasi-geoid model and the GPS-levelling data using a 1-parameter fit (Table 3).
For the computation of quasi-geoid model by the LSC method, residual gravity anomalies are used for the empirical covariance function estimation. Subsequently, the model covariance function is estimated using the covariance function model of Tscherning and Rapp (1974) through the program the COVFIT (Knudsen, 1987). It is found to be a quite good compromise for the covariance function Figure 3. Signal empirical (blue) and model (red) regional covariance function of reduced gravity anomalies after the removal of the contribution of EGM2008 to degree 360 and terrain effects with respect to a mean elevation surface at 30' resolution of the entire area as it is evident from Fig. 3. The estimated covariance function parameters are shown in Table 1.
These parameters (Table 1) are used as input for GEOCOL17 program of GRAVSOFT package (Forsberg and Tscherning, 2008).
In addition, the observation standard error of the reduced gravity anomalies is set to 1.0 mgal as Duquenne (2007)  A standard deviation of 3.4 cm is found between the quasi-geoid model determined by the LSC method and the GPS-levelling data using a 1-parameter fit (Table 3).

Regional Quasi-Geoid Model Determination Using the LSMS Method
A remove-grid-restore strategy is utilised to reduce the discretisation errors during the gridding operation while preparing the input data for the LSMS method (Ågren et al., 2009a (Forsberg and Olesen, 2010). We found that K =85 yielded a standard devitation of 3.5 cm between LSMS derived quasi-geoid model and GPS-levelling using a 1-parameter fit, indicating that K =85 is realistic taking into account the error characteristics of GRACE (Tapley et al., 2005). The corresponding degree variance models are specified in detail in   Tscherning and Rapp (1974). Degree variances rescaled using the factor 0.5 2

EGM error
Published ones for EGM2008 (Pavlis et al. 2008). Degree variances rescaled using the factor 0.4 2 Terrestrial gravity error Combination of the reciprocal distance and white noise models (Ågren et al., 2009a).
The reciprocal distance part is specified by the standard deviation 1 mGal and the cor- The white noise part is specified by the standard deviation 1 mGal and the Nyquist degree 10800   Table 3. Comparison of quasi-geoid models with GPS-levelling. Statistics for the 75 GPS-levelling residuals after a 1-and a 4-parameter fit. Note that max, min vales are given after the mean value is subtracted from differences for the 1-parameter fit.
(1) denotes that no modication of Stokes' integral is applied for the FFT method, while (2) indicates that the FFT method applies the modification. Unit is in metre  Table 5. Comparison of quasi-geoid models with EGM2008. Note that max, min vales are given after the mean value is subtracted from differences.
(1) denotes that no modication of Stokes' integral is applied for the FFT method, while (2) indicates that the FFT method applied the modification. Unit is in metre Quasi-geoid models Mean Max. Min. Std. dev.

Quasi-geoid defference analysis
In order to test the performance of the three different quasigeoid determination methods, their corresponding quasi-geoid undulations are compared with the quasi-geoid heights of the 75 GPS-levelling points situated in the central area ( Fig. 1(b)). Except a mean difference between the gravimetric quasi-geoid models and the quasi-geoid heights of the GPS-levelling points, all three methods agree well with the GPS-levelling data in the order of 3−3.5 cm using a 1-parameter fit provided that the modification of the Stokes' integral is implemented by the FFT method (Table 3).
There is a large difference between LSMS and RCR methods in terms of mean difference between the gravimetric quasi-geoid models and GPS-levelling quasi-geoid undulations. This is due to the differences in the coordinate reference system used. The agreement between FFT derived quasi-geoid model and the GPS-levelling quasi-geoid heights drastically improves when the modification is used (FFT (2) ). In case of using a 1-parameter fit to the differences between the gravimetric and GPS-levelling quasi-geoid heights, the FFT method using kernel modification (FFT (2) ) performs best with 3.0 cm r.m.s whereas the LSC method, the LSMS method and EGM2008 quasi-geoid heights shows 3.4 cm, 3.5 cm and 3.6 cm (Table 3).
In order to investigate the tendencies of the residuals between the gravimetric quasi-geoid models and quasi-geoid heights of GPSlevelling points after a 1-parameter fit, the differences between the quasi-geoid models and the quasi-geoid heights of GPS-levelling in the observation points are shown in Figs. 6(a)− 6(d), where it can be seen that the agreement is good. The pattern of the residuals for the FFT (Fig. 6(a)) and the LSC (Fig. 6(b)) methods are similar, not showing any systematic tendency. On the other hand, the residuals for the LSMS method ( Fig. 6(c)) and EGM2008 model (Fig. 6(d)) have a similar pattern, showing a remaining significant systematic slope in the direction from north-east to south-west, which is not the case in the FFT residuals ( Fig. 6(a)). Therefore, we use 4-parameter  Table 3. The 4-parameter fit corresponds to estimating the zero and first spherical harmonic degree terms using the observation equation (Section 2−18 in Heiskanen and Moritz, 1967), where 1 , 2 , 3 and 4 are the four parameters, φ is the latitude and λ is the longitude. In the 1-parameter fit only the x1 term is estimated. Using a 4-parameter fit, the LSMS method performs best with a 2.4 cm, whereas the FFT method, LSC method and EGM2008 show 2.9, 3.1 and 3.1 cm, respectively, in terms of the standard deviation between gravimetric and GPS-levelling quasigeoid heights at observation points. The significant systematic slope shown on Fig. 6(c) is absorbed in the 4-parameter fit, leading to better fit for LSMS than for FFT.
The mean differences and the standard deviations of the quasigeoid differences between the methods are shown in Table 4 Unit is in metre is 7 cm.
Furthermore, we used EGM2008 as a 'standard' over the study area, against which the performance of quasi-geoid models based on the three different methods are compared. Table 5 shows that the best agreement between the quasi-geoid models and the EGM2008 is observed for the LSMS method, showing a standard deviation of 3 cm. The worst agreement is for the FFT (1) method with a standard deviation of 8 cm while the LSC method gives a standard deviation of 6 cm. Obviously, this says nothing about the real accuracy of the models, as EGM2008 also has errors (but actually performs remarkably well, cf. Table 3).
In order to investigate the spatial variation of the differences between the quasi-geoid models obtained the by three methods with and without applying modification of Stokes' function for the FTT method are shown in Figs. 7(a)− 7(d). The differences between quasi-geoid models obtained from the FFT method without the use of modification (FFT (1) ) and from the LSC method has a pattern indicating a slope in the north-south direction ( Fig. 7(a)). This may potentially be due to the inherent periodicity assumption in the FFT method, while collocation implicitly assumes zero values outside the data area. Fig. 7(b)) shows that the differences between the FFT and the LSC methods increases when a modification is implemented by the FFT method (FFT (2) ). There is a large differerence between the FFT (2) and the LSMS methods reaching to 33 cm in a mountainous area ( Fig. 1(b) and Fig. 7(c)). Fig. 7 Figure 8. Comparison of the quasi-geoid models with EGM2008 quasi-geoid heights. a) FFT (1) minus EGM2008; b) FFT (2) minus EGM2008; c) LSC minus EGM2008; d) LSMS minus EGM2008. Unit is in metre reduced surface gravity anomalies ignoring the 1 term in Eq. (2).
However, this term is known to be small when terrain-reduced data are used (Forsberg and Sideris, 1989;Schwarz et al., 1990). Another issue is the fact that the LSMS method requires a comparatively dense surface gravity anomaly grid, which is computed from the observed gravity anomalies using a remove-compute-restore method for the gridding; see further Section 2.2.2. Here the grid resolution in question was chosen to approximately 2 km, which is likely too sparse in the high mountains. Notice that this does not mean that the LSMS method requires denser gravity observations than the other two methods.
Furthermore, we compare the quasi-geoid model EGM2008 to degree 2190 with the quasi-geoid models by the three methods in Figs. 8(a)− 8(d). The differences between the quasi-geoid models obtained by the FFT method, when no modification is used, and EGM2008 are relatively large both in flat and mountainous areas ( Fig. 8(a)). When the FFT method is used with the use of modification (FFT (2) ), the differences decrease in the flat areas but still exist over the mountains (Fig. 8(b)). The differences between the LSC derived quasi-geoid heights and the EGM2008 have a similar pattern with the differences between FFT (1) and EGM2008 ( Fig. 8(c)). The differences between the LSMS derived quasi-geoid heights and the EGM2008 quasi-geoid heights are relatively small, except in the mountainous areas where the differences reach to 20 cm ( Fig. 8(d)). However, these differences are mostly of a high frequency nature, corresponding to the different resolutions of the two models (0.02 • for LSMS and ≈0.08 • for EGM2008).

Conclusions
In this study, two methods using the RCR technique (2-D Spherical Journal of Geodetic Science 63 using a 4-parameter fit, the LSMS method gives the best agreement with GPS-levelling with 2.4 cm standard deviation. The reason for the remarkable improvement of LSMS after using a 4-parameter fit instead of 1-parameter fit is that the significant north-east to south-west slope in the residuals of the 1-parameter fit is absorbed by the 4-parameter fit.
These results provide only marginally better agreement with GPSlevelling than EGM2008 shows. This illustrates the remarkably high quality of EGM2008 in well-covered gravity data areas like France.
When the quasi-geoid models are compared with each other, although they agree well over the more central areas, they differ up to 33 cm over the mountainous areas near the Alps. Comparison of quasi-geoid models with EGM2008 quasi-geoid heights suggest that the LSMS method gives the best agreement with EGM2008 in term of the standard deviation over the mountainous region near the Alps. It would therefore be useful to extend the present investigations to the Alps, if gravity and precise GPS-levelling data could be released from this region as well in order to properly understand the large differences between the methods which is a very important future work.