The Best Robust Estimation Method to Determine Local Surface

Geoid/quasi-geoid modelling in a local (small) area is often made with polynomials. The optimal version of a polynomial is found by testing its successive versions using different statistical parameters, i.e. with different degree and number of coefficients. In this publication, the authors presented 3 approaches to search for an optimal version of the polynomial. Two of them were based on one of the robust estimation methods, i.e. the Danish method. Therefore, the question arises whether this method is suitable for this type of work. The publication analyses 6 different methods of robust estimation and the method of least squares. The control parameters of the methods were selected in such a way as to obtain a comparable accuracy of fitting the polynomial into the points each time. Then, the actual accuracy obtained for each of them was compared. The result of the study is a ranking of methods in terms of the most accurate fitting results.


Introduction
In most cases, it is not possible to measure processes over a large area with high accuracy. For this reason, at selected locations, a process is precisely determined and then approximated to the desired surface. There are two related issues: selection of the appropriate model, i.e. the approximation function and its optimisation. Optimisation should be understood as the selection and development of function parameters in such a way that they approximate the studied process as accurately as possible. Therefore, it is important to remove or minimise gross errors that disrupt the model. In geodesy, prominent fields of study are connected with the above question, i.e. modelling the shape of the earth (a geoid). Despite the development of surface measurements, including satellite measurements (Godah and Kryński, 2011;Pail et al., 2010), ground measurements (gravimetry, GNSS) are still the basis for its precise determination.
The above principle is even more visible when we model a quasi-geoid or a figure of Earth according to the theory of Molodensky (Molodensky, 1958;Molodensky et al., 1962). In this case, surface measurements are not necessary, but only precise measurements in a finite number of points. This is related to the renouncement of Bouger's gravimetric correction (large-surface measurements) in favour of freeair (measurements in points) (Czarnecki, 2010). The discussion on which of the reference surfaces is "better" is still ongoing (Vaníček et al., 2012;Sjoberg, 2013) and is not the subject of the article.
Another type of surface requiring optimisation is the model of transformation of height systems. Currently (2019)(2020) in Poland there is a change from Kronstadt'86 (zero point -Kronstadt, Russian Federation) to PL-EVRF2007-NH (Amsterdam, the Netherlands) (Ordinance, 2012).

Examination of the geoid
Geoid/quasi-geoid models can be developed for local, regional, or global areas. For local one's undulation is determined by set of points using satellite levelling (1a, 1b), but other methods are also available. Orthometric (H o ) or normal (H n ) heights are taken from geodetic vertical network. The ellipsoidal heights are taken from satellite measurements. The difference (undulation) can be determined at points on the ground, acting as ground survey markers. It can also be determined in place of GNSS permanent stations, but this requires the determination of their orthometric or normal heights. The use of direct levelling (spirit levelling) is often difficult or even impossible in this case. Hence, indirect methods are applied (Borowski and Banasik, 2011). Although the ellipsoidal heights of stations change over time (Haritonova et al., 2015;Maciuk and Szombara, 2018), the undulation values are considered relatively constant (Banasik, 1999;Czarnecki, 2010).
,undulation or distance between quasi-geoid and ellipsoid, , orthometric, normal height, ellipsoid heights. Surface model is done by fitting in those points approximate function. Sometimes an interpolation as a second way of fitting is mentioned. Because it is a special case of approximation (Fortuna et al., 1993), this issue is being omitted. Different types of functions can be selected for approximation. Commonly used are orthogonal and nonorthogonal polynomials. A non-orthogonal polynomial for ∀n ∈ Ν according to the variables x and y can be written as follows: (2) Where: searched value e.g. undulation, polynomial coefficient, coordinates of the points in planar system. For relatively local (small) areas plan coordinates are used. Inaccurate information on the location of a point (even up to 100 m) from a certain XY value does not significantly influence such a model (Banasik et al., 2012). The impact of mapping distortion is therefore negligible. The polynomial itself can be of any degree, but for local areas the polynomial of degree 3 and lower is commonly used. For larger areas or more differentiation surface the area is divided into samples and polynomial spline is used (Walo, 2000). In table 1 the full form of non-orthogonal polynomial of 1st, 2nd and 3rd degree is presented. The determination of the coefficients is performed on the basis of the correction equations in the form of (3) and (4), through the solution commonly used in geodesy with the least squares method (5). The weights for individual points can be taken in many ways. The simplest way is to take a fixed value for all points ( = ). Another way is to estimate a priori the weight of each point, e.g. depending on the length of the GNSS session or the class of a point (Pażus et al., 2002). Another possibility is to use the robust estimation method, which allows testing each point for a gross error. Those points, for which the model residuals do not meet the assumed conditions, have their weights reduced in the iterative process, and the calculation process starts from the beginning. This is the main part of the publication, which is presented in more detail in Chapter 4.
Where: corrections vector, undulation vector from satellite levelling (1), coefficients matrix with unknown -part x and y of equation (2), vector of polynomial coefficients .
̂= ( ) − • (5) Where: ̂calculated (estimated) polynomial coefficients , matrix of weights of individual points. Since the baseline model does not have to explain best the course of a geoid/quasi-geoid in a given area, further examination is possible. The model is then recalculated in a different form of polynomial and compared to its predecessor. Different types of statistical tests are used. The selected options are presented below. The authors omitted the hybrid adjustment approach from geometric and physical observations, developed e.g. in DFHRS project (Ameti and Jager, 2016;Jäger, 2006Jäger, , 2013, as a too far away from the article scope.

The Walo method
In Walo's paper, the total application of satellite levelling and vertical line deviation was assumed (Walo, 2000). Thus, the system of corrections (3) was extended to two systems of linear equations (6), (7) which after some transformations can be described as equations of errors of observed values (undulations and plumb line deflections) and fictional errors of observation. The main part of the method is the way of selection of polynomial. Walo adopted the study from full forms of the 1st-degree polynomial to "3+" (splines based on 3rd-degree polynomial). The Danish (Krarup et al., 1980) method was used to minimise gross errors. The parameter comparing subsequent forms of the polynomial was F-test (8), consisting in the verification of the hypothesis on the equality of variance of two populations of H 0 (̂1 2 = ̂2 2 ), of different degrees of freedom (v1, v2) with the assumed level of confidence p=0.95. It is assumed that the estimators are subject to chi-square distribution.
Where: ( , )the height of the geoid at P i point of the coordinates (x k ,y l ), ( , )the height of the geoid from the reference model (usually gravimetric), Θ ( , )the vertical deviation in azimuth A at the point P j of the coordinates (x j ,y j ), the azimuth of vertical deviation at the P j point, the parameter of the approximation polynomial, the degree of polynomial, ̂1 2 ,̂2 2unit variance of subsequent polynomial models.
The tests were carried out for 4 samples (areas) of the surface area ranging from 44.2x47.8km (A1), 96x128km, 63x83km to 152x67km. The test results showed that in every test the appropriate approximation function is a spline based on a polynomial of 3rd degree. The only exception was in the smallest area, for which additional tests related to the appropriate density of points were performed. In the case of the lowest studied density of points, F-test indicated the polynomial of 3rd degree. This was due to the relatively low number of freedom degrees (only 2 for "3+"), and thus relatively high critical value of F min .

The Zhong method
A different method is presented in Zhong's paper (Zhong, 1997). In this case, an output full version of the polynomial, e.g. of 2nd degree, is assumed. The system of equations (3) (4) is solved with the use of the Danish method. Then the significance of each polynomial coefficient (̂) is determined using a test based on its cofactor (̂̂,) and the estimate of unit weight variance in the complete model (̂0) (9).
Another (new) model is developed without the minimum significant ratio F min (9). It is then compared with the output by using 5 parameters (10-14). The procedure is repeated until the parameters reach the minimum value.
where: RMSresidual mean squares variance of interpolation biases, AICthe so-called "Akaike Information Amount" a specially defined statistic variable, mean mean-square interpolation error, nnumber of the data points, tnumber of the unknown model parameters, ̂̂the sum of the weighted residual square. ̂0 2the estimate of unit weight variance in the complete model.
The tests were carried out for a relatively small area of 3.0x7.0 km containing 13 points and 12 points for method validation. The initial model was the full polynomial of the 2nd degree (6 coefficients). In other studies the method was used for 20x30 km (Banasik, 1999;Borowski, 2015).

Tusat and Mikailsoy investigation
The search for suitable parameters to develop a polynomial can be found in the paper of Tusat and Mikailsoy (Tusat and Mikailsoy, 2018). The authors tested 5 different models of the polynomial: three full models (of 1st, 2nd and 3rd degree) and two intermediate models (of 2nd and 3rd degree reduced). To compare the subsequent models 10 statistical parameters were used: sum of squared errors (ESS), root mean-squared error (RMSE), coefficient of determination (R 2 ), adjusted R-squared (R adj 2 ), mean absolute percentage error (MAPE), Willmott's index agreement (D), confidence index (C), Theil's forecast accuracy coefficient (UI), Snedecor F-test (F) and Akaike information criterion (AIC). Due to the significant volume of equations, only selected parameters were presented, i.e. R adj 2 (15), D (16), C (17) and UI (18). Optimisation of each model was made by iterative Levenberg-Marquardt algorithm (LMA) with step length Δ=0.1. The models were developed on 33 points, spread over an area of about 25x30 km. Finally, 5 parameters indicated a model of 3rd degree (ESS, R 2 , D, C, UI) and 4 parameters for no. 2 (RMSE, R adj 2 , AIC, F). The authors found the results to be ambiguous.
There is C parameter that indicates a 3rd degree, which is dependent on D and R^2. Decreasing the number of points to 25 and reworking the models indicated a polynomial of 2nd degree. Where: , sum of squared errors, total sum of squares; ,number of points, number of coefficients.

Robust estimation methods
The methods used in this paper to identify the coefficients of the method model are derived from a broad class of M-estimations (Jäger and González, 2006;Jäger et al., 2005). Huber, Hampel, Danish, Gaździcki, IGG3 and Andrews methods were used for calculations. The estimation of polynomial coefficients was carried out with the use of the least squares method with iterative change of weights. In order to make it possible to compare the results, the control parameters were adopted in such a way as to obtain 95% efficiency in relation to the standard normal distribution (Banaś and Ligas, 2014).

The Huber method
A very popular and one of the first methods used for the robust adjustment of observations is the Huber method presented for the first time by P.J. Huber in his paper (Huber, 1964). It was developed as a result of the combination of the least squares method and the minimum average deviation method. The iterative modification of weights follows the formula below: Where: p iweight of the i-th observation from the previous iteration (in the first calculation step taken from the least squares method), p imodified weight used in the next iteration, v icalculated adjustment to the observation, fcontrol parameter, determining the range of acceptable values of adjustments v. The value of this parameter for the method assuming the above-mentioned efficiency is 1.345.

Fig. 1. Weight function diagram for the Huber method
The weight of the observations in the iterative process is unchanged if the standardised adjustment estimator v ̅ i is within the limits set by parameter f. Otherwise the value of the weight is reduced according to the damping function (Fig. 1).

The Hampel method
The Hampel method was developed as an extension of the Huber method and proposed in Hampel's paper (Hampel, 1974). It also leads to a reduction in weight for observations with a gross error. In comparison to the Huber method, 2 intermediate intervals (left and right of the permissible range) were additionally introduced. In this method the weight function has a form: Where: p iweight of the i-th observation from the previous iteration (in the first calculation step taken from the least squares method), p imodified weight used in the next iteration, v i -calculated adjustment to the observation, f, g, hcontrol parameters determining the limits of additional compartments characterised by a different way of weight modification. In order to achieve the assumed efficiency, the following control parameters were assumed in sequence, i.e. 1.7, 3.4, 8.5.

The IGG3 method
This method was proposed by Yang (Yang, 1991(Yang, , 1994. It is a modification of the Hampel method. The weight function of the IGG3 method is based on three ranges. A standardised adjustment qualifying for the first range does not change the weight from the previous iteration. In the second range there is a smooth weight damping. If the adjustment exceeds the value of the second control parameter, the observation to which it refers will get a weight equal to 0 in the next iteration. Values of control parameters are constant and depend on the current distribution. In this paper, in order to obtain the assumed efficiency, the parameter g=1.75 and h=8 have been adopted.
Where: p iweight of the i-th observation from the previous iteration, p imodified weight, v istandardised observation deviation, g, hcontrol parameters.

The Gaździcki method
It was created as a development of the Danish method (Gaździcki, 1985). Calculations are performed with the modified least squares method using an equivalent weight matrix (eq. 29). The damping function in this method is in the following form: )density function of the standardised normal distribution P Gprobability that the value of adjustment v i is not caused by a gross error in another observation ∫ f(v) g f dvprobability that v i will assume a value between <f;g>, and the scales are modified according to the formula: Where: t(v i )damping function, p imodified weight, v istandardised residual of the model, f, g -control parameters determining the boundaries of additional intervals characterised by a different way of modifying the weights. Values of these parameters were adopted as f=1.5, g=3.0, P G =0.5.

The Danish method
The idea of the Danish method belongs to a Danish mathematician, physicist and geodesist, i.e. Torben Krarup. The weighting function in this method is as follows (Krarup et al., 1980): Control parameters l and g are usually within a range from 0.01 to 0.1 and g = 2.
In this work, the value was set as l=0.1, g=2, k=2.

The Andrews method
The weight function of the Andrews method is widely used. It was introduced in 1972 ( Andrews et al., 1972) in a well-known study called the Pronceton Robustness Study. The form of the weight function is based on the trigonometric sine function and can be presented as follows: In order to obtain the efficiency at the level of 95% in relation to the standard normal distribution, the c parameter of 1.34 was adopted.

The test assumptions
To perform the test of robust estimation methods, 5 areas were selected, being the second level of administrative division of Poland (powiats/counties). The points were taken from the national basic levelling network, for which the distance between the ellipsoid and the geoid (ζ) was taken from the quasi-geoid model PL-geoid-2011 (Kadaj and Świętoń, 2014). Then polynomial of the 2nd degree was fitted in the points. On its basis the values of the distance ( − ) without unknown random and gross errors were determined. They were replaced by random and gross errors (26). Each of the counties is characterised by a different number of points (from 52 to 225) and their spatial distribution (Fig. 7). The noise (m ne ) was introduced for each of the points, using a single sampling based on a generator built in Python 3.6 (BSD license -permissive free software license, which is compatible with the GNU General Public License (GPL)). Gross errors (m ge ) were introduced gradually from 5% to 40% of all points. In this way 8 degrees of error load were obtained for each of 5 counties. The values of the gross errors were increased with the successive levels of sample load (Fig. 8).
( , )value of the distance between the ellipsoid and the quasi-geoid loaded with errors, ,coordinates of points in a two-dimensional coordinate system 1992, − distance values "cleaned" of thick and systematic errors, random error (noise) from 〈− , 〉 mm, gross error with values from 〈− * , * 〉 mm, where 〈 , 〉.

Results
The tests covered 5 areas with 8 levels of frequencies of gross errors. For the purpose of comparison, LS and 6 methods of robust estimation were used. In total, thus, 280 models were developed. The main aim of the study is to indicate which methods give the results as close as possible to the actual ones. Two sets of data were analysed. Variant I contains a comparison of the model with error-loaded data ( Where: the distance from the ellipsoid to the quasi-geoid between the model (model) and the data: error-loaded (ne+ge) or primary (q-adj); n, j, lnumber (in order): points in the area, areas, degrees of input errors; The SAAE distribution for variant I is shown in Figure 9. In the methods, the control parameters were selected in such a way as to ensure that the accuracy of the obtained results was similar (Chapter 5.1). The SAAE merging for all methods is clearly visible (Fig. 9).
The differences between them are at the level of 2 orders of magnitude below the relevant results. They can therefore be considered negligible. The same applies to methods that use root mean square error (RMSE) as a measure of error.
It was demonstrated with variant I that all methods were brought together to ensure the same level of accuracy of fitting the functions. In case of lack of information on actual errors in the points, it could be concluded that the tested methods fit the functions equally well to the data. Information on these errors is included in variant II. The distribution of SAAE is shown in Figure 10. In this variant, the results obtained do not overlap. Clear differences are visible, especially between the extreme values -the largest in least squares and the smallest errors in the Huber method. The largest errors in robust estimation methods are obtained when using the Danish method. This is not a surprise -similar results have already been obtained (Borowski and Banaś, 2018). Between them there are 4 remaining methods. The ranking of these methods is less clear; mainly due to the Andrews method, which performs worse with fewer gross errors. However, with the increase in numbers, its effectiveness increases in comparison to other methods (Hampel,Gaździcki and IGG3). It is possible that the reason for this is the feature, i.e. the weighting of all points. In this method, each observation characterised by an adjustment from the range modifying the weight obtains it. This method is the only one of the tested groups based on trigonometric function, i.e. sine function. Thus, the ranking of these 4 methods is the result of the selection of the ranking criterion. If TSAAE (29) is used, the Hampel method achieves a better result, then at a comparable level the IGG3 and the Andrews methods, and then the Gaździcki method (Table 2). With the increase in the number of gross errors, the Andrews method came out ahead of the Hampel method. In the tested sample, the limit value is 30% of the gross errors in the sample.
It should be noted that the methods gave different results with a 5% error level, and they gave similar results with a 10% error level. The first level, i.e. 5%, is loaded with a negligible number of errors of significant value. In most cases, the error is at the maximum noise level. For the level of 10%, there are also higher values of gross errors (3-4 times the maximum noise value). It can, therefore be pointed out that for a small number of outliers, similar values are obtained regardless of the used method. The exception is LS, which is the result of the lack of damping of these observations. As the number (and thus the magnitude) of gross errors increases, the methods give different results. The convergence of the Hampel, IGG3 and Andrews methods was obtained for 30% of errors.
A summary and ranking of the tested methods is given in Table 2. There was a slight difference in TSAAE between IGG3 and Andrews, so it was qualified in one position.

Conclusions
The calculations and thus the analysis of the methods were carried out under the assumption that the parameters were selected in such a way as to maintain 95% efficiency in relation to the standard normal distribution. For different "steering parameters" final conclusions might be different. It is important to mention that usually various values are used depending on authors experiences.
As a general conclusion, the authors recommend using robust estimation methods to express surface rather than "standard" least-square method. In performed tests almost every time the results were better when robust estimation was used. Especially when number of gross errors was relatively high. Standard least-squares was better than some robust methods only when number of gross was relatively small and more important its value comparable to introduced "noise" errors. This is not a surprise. Other analysis confirm the same conclusion (Borowski and Banaś, 2018).
The Danish method was applied in two cases in the quoted studies on local geoid/quasigeoid modelling. In our analyses this method turned out to be the worst (except the least squares method), i.e. it was used to obtain the highest sum of real errors. Compared to the best method that was analysed (Huber), it gives an error greater by almost 50%. The ranking of the effectiveness of the methods is presented in Table 3.
The last conclusion is the suggestion to analyse the methods based on the control parameters recommended by their authors. This will create a situation where the assumed effectiveness of estimators at the level of 95% will not be maintained, but it will enable comparison of the effects of the methods functioning in real conditions applied on a daily basis.