Unbiased least-squares modification of Stokes’ formula

As the KTH method for geoid determination by combining Stokes integration of gravity data in a spherical cap around the computation point and a series of spherical harmonics suffers from a bias due to truncation of the data sets, this method is based on minimizing the global mean square error (MSE) of the estimator. However, if the harmonic series is increased to a sufficiently high degree, the truncation error can be considered as negligible, and the optimization based on the local variance of the geoid estimator makes fair sense. Such unbiased types of estimators, derived in this article, have the advantage to the MSE solutions not to rely on the imperfectly known gravity signal degree variances, but only the local error covariance matrices of the observables come to play. Obviously, the geoid solution defined by the local least variance is generally superior to the solution based on the global MSE. It is also shown, at least theoretically, that the unbiased geoid solutions based on the KTH method and remove–compute–restore technique with modification of Stokes formula are the same.


Introduction
Today ultra-high Earth Gravitational Models (EGMs) allow detailed geoid determination all over the Earth. However, the higher-order harmonics of the EGMs are typically much less accurate than the low-to-medium wavelengths, which calls for improving the geoid estimator by using additional terrestrial gravity data in a combined local/regional solution. This is the case for some versions of the remove-compute-restore (RCR) method (e.g., Forsberg 1993;Sansò and Sideris 2013) and the least-squares modification of Stokes' formula (LSMSF; Sjöberg 1980Sjöberg , 1984aSjöberg , b, 1991Sjöberg , 2003Sjöberg , 2005a. Due to truncations of the EGM series and the area of Stokes integration, such solutions have an inherent bias, and, in case of the LSMSF technique, the optimum solution is provided for the minimum of the expected global mean square error (MSE). However, Sjöberg (2005a) set forth the idea of optimizing the LSMSF method by minimizing the local MSE. While the LSMSF method implies that the EGM and gravity anomaly data are combined by spectral weighting at spherical harmonic degrees, Sjöberg (2011, Sect. 3) presented a method for local geoid determination by spectral combination and weighting by degree and order, being unbiased in case the truncation error is negligible. Similarly, Klees et al. (2018) tested the use of single-and multi-scale spherical radial base functions for local spectral combination of different sets of gravity-related data in geoid determination, but the result shows that the method needs further test and development before being suitable for practical application.
Here we will consider the LSMSF method in the case that Stokes' integral covers a sufficiently large region, such that the remote zone effect becomes negligible, in which case the solution can be regarded as unbiased (provided that all data are unbiased) from a statistical point of view. This implies that each solution (not necessarily the least-squares combination) will be (at least practically) unbiased, so that the MSE solution can be replaced by a local minimum variance solution. Sjöberg and Bagherbandi (2017, Sect. 4.4.6) shortly discussed the possibility of modifying Stokes' formula with a very high-degree EGM, but the strategy was still to include the estimation bias in the model and minimizing the MSE. In contrast, here we will assume that the bias is negligible, implying that the minimum variance solution is optimal. Hence, both the unbiased LSMSF and Remove-Compute-Restore (RCR) solutions will be derived in Sects. 2 and 3, and the results are shortly discussed in Sect. 5. Section 6 concludes the study.
Although all of the above methods need direct and indirect corrections for topography, atmosphere and ellipsoidal shape of the Earth as well as for downward continuation of gravity observations to the sphere of integration, only the topographic effects are shortly discussed in Sect. 4 (but are more extensively treated in Sjöberg and Bagherbandi 2017, Sect. 6.2).

Spectral combination
Consider the geoid height model where is the unit sphere, c = R∕(2 ) , being normal gravity on the reference ellipsoid, R is the mean Earth radius, Δg T is the terrestrial gravity anomaly observation, Δg EGM n is the EGM derived gravity anomaly Laplace harmonic of degree n, and L is the upper degree of the EGM as well as upper degree of the modification of Stokes' kernel function S( ) by parameters s k , i.e.: where P k are Legendre's polynomials and is the geocentric angle of integration w.r.t. the computation point.
The spectral form of Eq. (1a) reads: where n = 2∕(n − 1) , Δg n is the true gravity anomaly, T n and EGM n are the random errors of Δg T n and Δg EGM n , respectively, and Assuming that all error components are uncorrelated with vanishing expectations, and by subtracting the true geoid height from Eq. (2a), the following error and variance of Ñ are obtained: and where E{} is the statistical expectation operator, and 2 n and dc 2 n are the error degree variances of the terrestrial and EGM gravity anomalies, respectively.
From Eq. (4), one notices that the geoid error is unbiased, and from Eq. (5) follows that the minimum variance related with the least-squares solution N becomes for the modification parameters The spectral weights given by Eq. (7) were first presented by Sjöberg (1980) and and by Wenzel (1981), who named this technique spectral combination.
We now generalize the initial geoid model (1a) by assuming that the terrestrial gravity and also the EGM errors are internally correlated with covariance matrices and , respectively, (but gravity and EGM data are assumed as mutually uncorrelated). Then the geoid model (1a) is still unbiased, but now its variance becomes: or where we have decomposed the infinite matrix into four parts with LL and ∞∞ on the diagonal and L∞ = T ∞L off-diagonal, of dimensions (L,L), (∞, ∞) and ( L, ∞ ), respectively. (In practice, the infinite dimensions must be approximated by some finite degree.) Differentiating 2 N w.r.t. and equating to zero one arrives at the least-squares choice of : with the variance of the geoid height estimator

Solutions for a small integration cap size
We now allow for a high-degree maximum (M) of the EGM, such that M ≥ L , and we limit the Stokes integration in Eq. (1a) to a cap 0 of spherical radius 0 . Then the general estimator becomes: where and are the so-called Molodensky truncation coefficients with In this case, the geoid estimator can be written in the spectral form as and the geoid height error becomes Here it will be assumed that M, the maximum degree of the EGM, is so large, that the bias, the last term of Eq. (12), is negligible at the cm-level already for integral cap sizes of a few degrees (see Note 1 at the end of this section). Sjöberg (2005b, Sect. 5) reported that by selecting the modification parameters according to Molodensky's method (Molodensky et al. 1962), which is specifically designed to minimize the truncation error, with L = M = 360, the RMS truncation error is within 1 mm already for a cap size of 2°. However, in the LSMSF technique the bias will be larger as it is balanced with other error components in a minimum MSE solution.
Disregarding the bias and assuming that the errors of Δg T and Δg EGM are uncorrelated, the following variance of the general geoid estimator of Eq. (10a) Note 2 If the aim were to find the minimum MSE solution, the square of the bias must be added to the variance in Eq. (13) for differentiation w.r.t. s. However, this is not the goal of this study.

The RCR solution
Using the remove-compute-restore technique, the geoid estimator corresponding to Eq. (1a), where L = M, can be written where and the spectral form of this estimator becomes which, after some rearrangements of terms, equals Eq. (2a). Similarly, using the spectral representation, the more general Eq. (10a) for the KTH approach with L ≤ M can be rewritten in its RCR-form as That is, from a theoretical point of view the RCR technique yields the same set of estimators as the straightforward modification of Stokes' formula. In particular, the least-squares modification parameters and geoid variance become the same for the two methods.

Topographic corrections
As stated at the end of Sect. 1, the corrections for topography are treated only sparsely in this study. The KTH approach uses analytical continuation of the external disturbing potential to point level and geoid level in applications for quasiand geoid determination, respectively, and there is no further topographic correction in quasigeoid estimation, while geoid estimates need also corrections for the topographic bias. In the RCR method, firstly one removes the effect of the topography on the observables as the direct topographic effect on the Molodensky type of surface gravity anomalies ( Δg T and Δg EGM . Then the gravity anomalies are downward continued in one way or another to the sphere of computation, a secondary indirect topographic effect on Δg T is applied, and after Stokes integration, the first indirect topographic effect on the geoid is added (e.g., Sjöberg 2018). The interested reader is also referred to Sjöberg and Bagherbandi (2017, Sect. 6.2) to find that also these corrections theoretically agree in the KTH and RCR approaches when properly applied. See also Sjöberg (2005b).

Discussion
It must be emphasized that the above modification parameters and least-squares solutions are local, implying that the modification parameters change with position on the mean Earth sphere. However, in practice, one can expect small changes in a small region, implying that only one set of modification parameters is needed there. As all solutions are unbiased, any choice of modification parameters yields an unbiased solution (if the truncation bias vanishes). Hence, the weight relation between terrestrial and EGM gravity data will not be critical for the solution, but, as is known from any least-squares adjustment, the error estimation will be more dependent on the choice of error covariance models. An important advantage of the unbiased solutions is that the gravity anomaly signal degree variance model needed in the MSE solutions is not required. A similar approach for geoid determination was outlined in Sjöberg (2005a), but the principle difference is that there the geoid bias term was included in the optimization, which leads to a least-squares solution based on the (local) MSE. Hence, if the bias term is small, the two solutions should be practically the same.

Concluding remarks
Unbiased minimum variance estimators are usually preferred to biased/MSE solutions, demanding less a priori information for practical applications. The spectral combination without assuming correlations is the most simple solution, but not very realistic as there are mostly correlations among the data that should be taken into account. Doing so in the strict sense lead to huge, even infinite dimensional, covariance matrices, which, of course, in practice must be restricted to a suitable resolution/dimension. Here one should also remember that approximate/simplified covariance models also lead to unbiased, even if not optimal, solutions.
We have shown that the unbiased LSMSF and RCR solutions are theoretically the same, although there may be some differences in their practical implementations.