Application of variance components estimation to calibrate geoid error models

The method of using Global Positioning System-leveling data to obtain orthometric heights has been well studied. A simple formulation for the weighted least squares problem has been presented in an earlier work. This formulation allows one directly employing the errors-in-variables models which completely descript the covariance matrices of the observables. However, an important question that what accuracy level can be achieved has not yet to be satisfactorily solved by this traditional formulation. One of the main reasons for this is the incorrectness of the stochastic models in the adjustment, which in turn allows improving the stochastic models of measurement noises. Therefore the issue of determining the stochastic modeling of observables in the combined adjustment with heterogeneous height types will be a main focus point in this paper. Firstly, the well-known method of variance component estimation is employed to calibrate the errors of heterogeneous height data in a combined least square adjustment of ellipsoidal, orthometric and gravimetric geoid. Specifically, the iterative algorithms of minimum norm quadratic unbiased estimation are used to estimate the variance components for each of heterogeneous observations. Secondly, two different statistical models are presented to illustrate the theory. The first method directly uses the errors-in-variables as a priori covariance matrices and the second method analyzes the biases of variance components and then proposes bias-corrected variance component estimators. Several numerical test results show the capability and effectiveness of the variance components estimation procedure in combined adjustment for calibrating geoid error model.

where h is ellipsoidal height obtained from GPS measurements, H is orthometric height obtained from leveling survey, N is gravimetric geoid data calculated from global geopotential models or regional gravimetric geoid. Based on this inherent appeal of the relationship between ellipsoidal, orthometric heights and gravimetric geoid data, if any two of the heights are given, the third can be derived with Eq. (1).
In practical terms, the performance of Eq.
(1) is more complex due to the adaptation of the parametric models for the systematic errors (e.g. long-wavelength systematic errors in N, distortions in the orthometric height due to an over constrained adjustment of leveling network, etc.) and the correctness for the stochastic model of the observable noises due to the covariance (CV) matrices for the height types obtained from separate adjustments of the individual height types.
Up until today, the former problem has been a topic of interest and well studied by different kinds of parametric functions to absorb the systematic errors of ellipsoidal, orthometric heights and gravimetric geoidal undulations. Many models have been performed ranging from a simple linear regression to more complicated models, such as polynomial fit, multiquadric function fit, Spline fit method and so on (Sansò and Sideris 2013). Since it is evident that the suitability of the parameter model depends on the density, distribution and quality of height network data, there is no universal model applicable to any situation. Here the method of multi-surface function is used after evaluation and assessment of parametric models in this paper.
As for the latter issue, the problem of stochastic modeling for observables in the combined adjustment of heterogeneous height types has not yet to be satisfactorily solved by this traditional formulation. This problem will be the main focal point of this paper. In order to reach the best unbiased estimators of the unknown parameters, a proper CV matrix of the observables is required. Here the well-known method of variance component estimation (VCE) is employed to the combined least-squares (LS) adjustment in Eq. (1) (Teunissen and Amiri-Simkooei 2008; Amiri-Simkooei 2013; Wensch et al. 2013). There are many reasons for performing VCE. For example, VCE is an effective statistical tool to test the noise level. This method is very flexible and easily understood. It can be conducted to estimate variance and CV components for linear and nonlinear stochastic model. The main idea of this paper is to give a detailed analysis of the combination adjustment and to estimate all errors of the heterogeneous heights data and to develop a method for calibrating the geoid error models using the method of VCE. In practice, the iterative minimum norm quadratic unbiased estimation scheme is implemented via a combined adjustment using existing ellipsoidal, orthometric heights and gravimetric geoid data.

Combined LS adjustment of GPS-leveling heights and gravimetric geoid data
To establish the LS variance component estimation, a combined adjustment is carried though using the following functional model with orthometric, ellipsoidal heights and gravimetric geoid data. Considering the Gauss-Markov model, the system of observation equation and the solution can be expressed as, where A is the m × t matrix of known coefficients depends on the parametric model (m is the number of observational equations, t is the number of unknown parameters), x is t × 1 vector of the unknown parameters depending on the parametric model, E{•} is the mathematical expectation operator, Q v is the m × m CV symmetric cofactor matrix of the observables, σ is the variance or CV components, and the observations vector l consists of the height 'misclosure' at the GPS-leveling benchmark as follows, Here, for the design matrix A corresponding to the parametric model the model of Multi-surface function fit is being used, δ, ϕ and λ are latitude, longitude respectively, and a i are the coefficients.
v is a (m × 1) vector of unobservable random error with zero mean, for each of heterogeneous height data types are given by, where B is the block-structured matrix B = [I−I−I], such that I is an m × m unit matrix.
The corresponding CV matrix is described in Eq. (2-c), which can be written as, where C is a positive-definite symmetric matrix. Here we assume that there are no correlations between the heterogeneous height types. Q(⋅) are known positive cofactor matrices for ellipsoidal, orthometric heights and gravimetric geoid data, and σ 2 h , σ 2 H , σ 2 N are the corresponding variance components.
The vector of unknown parameters is computed by using LS minimization principle, where the corresponding weight matrix P associated with the observations take the forms, where P(⋅) is assumed diagonal with elements which is the reciprocal of the error variances. According to LS adjustment, the unknown parameters and adjusted residuals of observations can be easily solved. The unknown parameters is given as follows, The adjusted residuals of the GPS, leveling and gravimetric geoid height observations can be obtained by and with And the mean square error and the accuracy of parameters estimations are given as follow, According to this combined LS adjustment approach, the solution can be achieved depending on two issues, one is the appropriateness of the parametric models, which refers to the correction for the data inconsistencies and the systematic errors, and another is the residuals of the height data types which allow for the calibration of data covariance matrices. Since the former problem has been well studied, the main work of next section focus on a description of the implement of VCE schemes to the combined LS height adjustment.

Application of MINQUE to the Combined LS Adjustment
There are many methods available to implement VCE within the LS adjustment (Helmert 1924;Rao 1970). The first solution was proposed by Helmert. And an independent solution was proposed by Rao (1970) who put forward a method known as minimum norm quadratic unbiased estimation (MINQUE). Assuming the observations are normally distributed the approach of Helmert's and Rao's are equivalent. In this paper the MINQUE procedure is employed since this procedure does not require distributional assumptions.
Here the general MINQUE algorithm is modified to estimate the variance components in the combined LS adjustment with heterogeneous observation types, where θ is a vector composed of the unknown variance components σ 2 h , σ 2 H and σ 2 N . The matrix S is denoted by where S ij in the matrix is derived by where tr(·) is the trace operator, i, j = h, H , N, and Q(·) is the cofactor matrix for observables. Considering the matrix S may not be of full rank, the algorithm of pseudo-inverse should be used to solve Eq. (15). R is a matrix denoted by where A is an appropriate matrix corresponding to the parametric model as in Eq. (4) and C 1 is the CV matrix of the observables. A CV matrix model for the heterogeneous observation types can be written as the following linear model, The vector q can be expressed as where v i = Q i Rw, which is a vector composed of the residuals for heterogeneous observations. Substituting the formulations above into Eq. (15), we can obtain the explicit expression by It is evident from the expression in Eq. (18) that an iterative process should be used because the unknown variance components σ 2 h , σ 2 H and σ 2 N are embedded in C 1 . So, the estimation for the variance components must be conducted with a convergence criterion. Here the iterative performance of MINQUE is employed to determine the variance components. In practice, a convergence criterion should be specified to terminate the computation. In this paper, the computation should not be stopped until all the values are equal.

Description of data
Here a rugged area bounded by 39°N and 41°N, and 82°W and 85°W, is chosen for testing the theory. Several numerical studies are implemented with the two datasets composed of GPS-leveling heights, gravimetric geoid data obtained from gravimetrical data and the initial cofactor matrices for ellipsoidal, orthometric heights and gravimetric geoid. A data set of Free Air (FA) gravity anomalies used for geoid computation are consisted of 4,574 ground data irregularly distributed in the test area and are referred to GRS80, and the data set of FA gravity anomalies is referred to the NAD83. All data have been removed the duplicate points and validated Least squares collocation (LSC). The error ratio used in this method is about 8 % and the remaining points are used for test (Tscherning 1991;Gil et al. 1993). Figure 1 gives the isoline map of FA gravity anomalous gridded by Kriging interpolation.
This gravity anomaly is used to estimate the geoid by application of Stokes's formulae. Before implementing Stokes's formula, the gravimetric geoid computations require that the topography effect should be reduced by removing step of the remove-compute-restore (RCR) technique for the purpose of making the gravity anomalous small and easy to be gridded (Guo and Xu 2011b). Here a method of Residual Terrain Model (RTM) is employed. With this technology only the topographic irregularities relative to a smooth reference surface (Omang and Forsberg 2000) are taken into account. Table 1 shows the statistical results of the FA gravity data and the residual gravity anomaly, where Δg fa is the FA anomaly, Δg ref is the effect of the reference field of EGM2008 on the gravity anomalies (Pavlis et al. 2008), Δg rtm is the terrain effect of RTM reduction with the smooth reference surface, e.g. of resolution 120 km, and Δg res is the residual gravity anomalous obtained by subtracting the effect of a global Earth gravity field model Δg ref and the RTM effect of the topography Δg rtm from Δg fa , g res = g fa − g ref − g rtm . Fig. 1 The FA gravity anomalous of tested area (unit: mGal).
From Table 1, it is obvious that the reduced gravity anomaly is significantly smoother than the FA gravity anomalous.
Then, the residual gravity anomalous are gridded on a 2.5′ × 2.5′ grid using Kriging interpolation in order to meet the requirements of geoid determination employing the technique of Fast Fourier transform (FFT) (Haagmans et al. 1993). By modifying Stokes formula as convolution integrals, the FFT techniques are used to estimate the residual geoidal undulation with the reduced gravity anomalous here. Finally, the gravimetric geoid can be obtained after applying the restoring step. The gravimetric geoid with a resolution of 2.5′ × 2.5′ is shown in Fig. 2.
The GPS-leveling heights are employed to estimate the absolute and relative accuracies of gravimetric geoid. A total number of 170 GPS-leveling benchmarks distributed throughout the test area were used. Figure 3 gives the distribution of GPS-leveling stations. The GPS observations were processed with the Bernese GPS software version 4.2 with observation periods between 4 and 16 h. All the GPS heights used in a LS adjustment are given with respect to the GRS80 reference ellipsoid, the reference frame is ITRF2005. Geodetic leveling observations are given with respect to the North American Vertical Datum of 1988 (NAVD88).

Table 1 Statistics of the gravity anomalous (unit: mGal)
Δg fa FA anomalies, Δg ref effect of the reference field of EGM2008 on the gravity anomalies, Δg rtm terrain effect of RTM reduction, Δg res residual gravity anomalies.

Min
Max Mean STD Fully populated initial CV matrices are computed for ellipsoidal, orthometric heights and gravimetric geoid data. The CV matrix for the ellipsoidal heights is obtained from the results of the GPS post-processing software package being used. These are the direct results without considering the spatial, temporal and physical correlations between GPS phase observables. Similarly, the CV matrix for the orthometric height are obtained from the rigorous adjustment of leveling measurements and the CV matrix for the gravimetric geoid heights at the GPS benchmarks are derived in a straightforward manner by applying the law of error propagation to the geoid solution. Table 2 summarizes the characteristics of the 'a-priori' CV matrices of all the previous data types for comparison purposes.

Iterative MINQUE scheme
When estimating variance components, convergence depends on the given initial values of weight of the measurements. Here we will study an algorithm to improve the convergence speed and behavior of estimating variance components.
Considering the availability of fully populated CV matrices is a luxury, the iterative application of MINQUE is used to estimate the variance components in practice. First of all, suppose that by using a conventional technique of variance component estimation, we have the three estimates at the jth step σ j h , σ j H , σ j N . To continue the next iteration, we Table 2 Initial CV matrix characteristics (unit: cm)

Fig. 3 Distribution of GPS benchmarks.
compute the new set of two variance components, from which, the weights for the next iteration can be further computed, The iterative procedures are given in Fig. 4. From Fig. 4 we can see that the iterative procedures need a priori values for the estimation of variance components which have been shown in Table 2, and need to specify a convergence criterion to determine when to terminate the process. In this scheme, the process should be stopped until all the values satisfy the converge of In order to study the correlations among observations, numerical tests are performed with diagonal CV matrices. Figure 5 shows the convergence behaviors of the estimated variance components by MINQUE in each of successive iterations for observable. The results demonstrate the computational efficiency.

Numerical examples for applying MINQUE to combined LS adjustment
A popular approach is to compute the difference between the gravimetrical geoidal undulation N, with ellipsoidal height h from GPS measurement and orthometric height H from leveling measurement N = h−H − N. The discrepancies ΔN denote an evaluation for the accuracy of gravimetric geoid. In order to estimate the accuracy of gravimetric geoid, the different between gravimetric geoid solutions are compared with 165 GPS-leveling heights.
Each of the height data types refers to different reference surfaces which results in that biases are introduced into the geoid model. In order to absorb the errors due to the datum inconsistencies between the gravimetric geoid and GPS-leveling heights, a parametric model has been used. In general, the process of choosing the parametric model suffers from a high degree of arbitrariness. In order to evaluate the availability of parametric models, the semi-automated assessment procedure combined with F-significance test is used. Figure 6 gives the computational process of semi-automated assessment procedure. The detailed procedures are not elaborated in this paper. The interested reader should refer to the work by Dermanis and Rossikopoulos (1991). According to the results, multi-surface function is selected for computing the adjusted residuals here.
Since the main idea of this paper is to discuss the influence of a priori CV matrices on the final computed variance components, two numerical examples are conducted to assess the gravimetric geoid accuracy in this section. In the scheme of case study II, combined height adjustment is conducted by using the a priori CV matrices as given in Table 2. In the scheme of case study I, combined height adjustment is conducted by using variance-component estimated by MINQUE method. The widely used approach in case study I is the comparison of different geoid model computed by the scheme of case study II over the same region, which is an external method to evaluate the accuracy of gravimetric geoid.  Such a scheme selects the best GPS benchmarks to evaluate the accuracy of geoid with minimum range and standard deviation criteria. The selected GPS benchmarks are used to fit geoid using the multi-surface function fit technique. 150 GPS benchmarks are selected to fit geoid models in this paper. The residual values at the all the rest GPS benchmarks denote the relative accuracy of the geoid. The statistical results of the residuals at the 15 GPS benchmarks are shown in Table 3. We can see that the accuracy of geoid fitting with GPS-leveling heights based on VCE is with minimum range and standard deviation.

Conclusions
The main focus point of this study is to give a procedure of geoid determined with available high precise gravimetric geoid and GPS and leveling heights and to improve the accuracy of geoid. In the context, a precise formula for the geoid computation is derived, and the relationship between ellipsoidal, orthometric heights and gravimetric geoid data is conducted. And a detailed algorithm uses the complete description of the CV matrices of the observation vector and of the coefficient matrix, possibly with unknown components of each. In the practice of this context, the iterative application of MINQUE is employed to test the variance components matrices for heterogeneous height data. Numerical case studies shows that the best results were obtained combined with the variance components estimation, with an outside precision of 1.91 cm and an inside precision of 0.29 cm, while the classical method with a-prior CV matrix is a bit worse, with an outside precision of 4.70 cm and an inside precision of 1.25 cm. The results of numerical examples show the capability and effectiveness of variance components estimation procedure in combined adjustment for calibrating geoid error model.

Authors' contributions
DG has made substantial contributions to conception and design, or acquisition of data, or analysis and interpretation of data and involved in drafting the manuscript. HX has given final approval of the version to be published; and agree to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. Both authors read and approved the final manuscript.