ERRORS MEASUREMENT OF INTERPOLATION METHODS FOR GEOID MODELS: STUDY CASE IN THE BRAZILIAN REGION Mensuração do erro de métodos de interpolação para modelos geoidais: Estudo de caso na região brasileira

The geoid is an equipotential surface regarded as the altimetric reference for geodetic surveys and it therefore, has several practical applications for engineers. In recent decades the geodetic community has concentrated efforts on the development of highly accurate geoid models through modern techniques. These models are supplied through regular grids which users need to make interpolations. Yet, little information can be obtained regarding the most appropriate interpolation method to extract information from the regular grid of geoidal models. The use of an interpolator that does not represent the geoid surface appropriately can impair the quality of geoid undulations and consequently the height transformation. This work aims to quantify the magnitude of error that comes from a regular mesh of geoid models. The analysis consisted of performing a comparison between the interpolation of the MAPGEO2015 program and three interpolation methods: bilinear, cubic spline and neural networks Radial Basis Function. As a result of the experiments, it was concluded that 2.5 cm of the 18 cm error of the MAPGEO2015 validation is caused by the use of interpolations in the 5'x5' grid.


Introduction
Obtaining geopotential models that best represent the equipotential surface of the Earth's gravitational field has aroused interest among Geodetic Scientists, especially in applications that need to set the fluid flow direction.Some users of the geoid models, such as oceanographers and geophysicists, needs quality of decimeter level for geoid heights (Kuroishi 2009); however, certain geodetic works and high-accuracy mapping applications, such as GNSS/leveling, global unification of altitude systems and deglaciation studies, require an accuracy of about 1 cm (Plag and Pearlman 2009).Thus, it is essential to seek methods that provide accurate results for geoid models.
Modified Stokes integral and spherical approximations are usually used to obtain geoid models.As a result, a regular grid of geoid heights is available to users in order to perform interpolations to obtain the geoid height in certain positions.Although the latest geoid models use a modern method for obtaining a regular grid of high quality, the extraction of information should be done by users with interpolations which faithfully represent the surface of the mesh.This task may not be simple, because each interpolator represents the surface according to its own characteristics.
Despite progress in developing geoid models in recent decades, little attention has been given to the choice of the interpolating method used for the extraction of the mesh information.Against this background, two questions are addressed in this work: (1) what is the magnitude of the interpolation error on a geoid model?(2) is there a significant difference in information extracted by users using different interpolations?In order to discuss these two issues in the Brazilian reality, this work presents an analysis for comparing some interpolation methods that can be used to obtain heights from the geoid model of Brazil in its latest version, MAPGEO2015.There are three specific objectives: quality analysis of interpolation using meshes with different spatial resolutions, assessing the interpolation quality across the standard interpolation of MAPGEO2015 and the validation of interpolating using real data.For this purpose, the next chapter presents a brief history of MAPGEO and Chapter 3 presents the methodology used to conduct the analysis of the different interpolators analyzed: the bilinear, the cubic spline and neural networks Radial based Function (RBF).Chapter 4 shows the results and Chapter 5 presents the conclusions.

History of MAPGEO
MAPGEO92 was one of the first versions of geoid model in Brazil.This model showed an accuracy of ±3.0 m (Blitzkow et al. 1993).Subsequently, MAPGEO98 emerged with an accuracy of ± 1.5 m (Diniz et al. 2004).The dataset used to generate the model in the earliest versions was not disclosed, however, since the third version (MAPGEO2004) came out, dataset information has become available.According to Lobianco, Blitzkow and Matos (2005) this model was calculated using the Fast Fourier Transformation (FFT), modified Stokes kernel and remove-restore technique.Also, the dataset used was the Earth Gravitational Model 1996 (EGM96) which considered up to degree and order 180, a Digital Terrain Model (DTM) of 1' developed by Escola Politécnica da Universidade de São Paulo (EPUSP) and Helmert's gravity anomalies with 10' resolution.The geoid model accuracy was evaluated by Matos et al. (2012), showing an accuracy of ± 0.68 m.
The successor geoid model in the MAPGEO series was released in 2010, named MAPGEO2010, whose regular grid was made available with a resolution of 5'.This model was generated using the same configurations as MAPGEO2004, such as FFT, modified Stokes kernel and remove-restore technique, but the dataset used was the EGM08 which considered up to degree and order 150, 925,878 gravimetric points and also the free air gravity anomalies model named DNSC08GRA (Andersen et al. 2009).This national model series was developed by a partnership of the Instituto Brasileiro de Geografia e Estatística (IBGE) and EPUSP.The most recent geoid model was released in November 2015, named MAPGEO2015.Although MAPGEO2015 has an equal spatial resolution of 5', it presents a 20% improvement in the root mean square difference (RMSD) compared to MAPGEO2010.Among the novelties of MAPGEO2015, we highlight the use of Artificial Neural Networks (ANN) in South and Southeast of Brazil for obtaining Helmert anomalies in regions with gravimetric gaps (Machado, Blitzkow and Matos, 2013).The same methodology as before was kept in the generation of MAPGEO2015, i. e., the use of FFT and modified Stokes kernel linked to a removerestore technique.More details about the MAPGEO can be obtained from Blitzkow (2003), Guimarães, Matos and Blitzkow (2013), IBGE (2015) and Blitzkow et al. (2016).Figure 1 shows the grid of the geoid height model.
Bulletin of Geodetic Sciences, 24(1): 44-57, Jan-Mar, 2018 Together with MAPGEO2015 regular mesh, IBGE also provides users with an executable, the last version is mapgeo2015_v1, capable of calculating values of geoid height between longitudes 35° 34' 0.12 ' W and 75° 25' 1.19' W and latitudes 5° 34' 59.89'' N and 35° 25' 1.19 '' S. According to IBGE (Electronic correspondence replied by Technical Team of the Gravimetric Network and Geoid Development of IBGE on May 13, 2016 -Customer Service Number: #90996/2016 -3#) the mapgeo2015_v1 uses a bilinear interpolation in a 5'x5' grid, allowing users to determine geoid height values for distinct points in the Brazilian region.It is worth mentioning that the latitudes and longitudes of the mesh are provided to the seventh decimal place in degrees, while the geoid undulations are provided to the second decimal place.

Methods of Investigation
The methods of interpolation evaluated in this work were: bilinear, cubic spline and RBF network interpolations.Those methods were selected in order to verify whether there is any difference in using robust or simple methods for geoid interpolation.For more details about those interpolators, the authors suggest Wright (2003) or Haykin (2009) for RBF networks; Pedrini and Schwartz (2008) for bilinear; and Ruggiero and Lopes (1998) or Subbotin (2011) for splines.
The analysis of the methods was divided into three experiments.The first consisted initially of a 10'x10' grid generation from the MAPGEO2015 5'x5' grid.The interpolators were applied to the 10'x10' grid (control points) and the height values of the original grid (validation points) were used as a reference for the calculation of the interpolator errors.The interpolator errors were calculated comparing the validation points with the interpolated ones, ignoring the common points in both grids.In addition, the validation points located within 1° close of the boundary were discarded in the experiments.Figure 2 presents an example of the resample achieved.The second experiment aimed to analyze the mapgeo2015_v1 interpolation method and the bilinear, spline and RBF interpolations.The MAPGEO2015 5'x5' grid was used to interpolate geoid heights each 2.5' in latitude and longitude, generating a 2,5'x2,5' grid of interpolated values.These values were compared with the bilinear, spline and RBF interpolation values.In this case, all the geoid heights were interpolated from the original grid and the validation was performed by calculating the differences between the interpolators.RMSD (Eq.02) and absolute maximum difference were used as quality parameters for the evaluations.For the third experiment, the original grid of MAPGEO2015 was also used and all the interpolators, including the mapgeo2015_v1 interpolation.However, the geoid heights were obtained for 539 levelling points of the Altimetric Fundamental Brazilian Network tracked by GNSS. Figure 3 shows the distribution of the levelling points used in the third experiment.
Bulletin of Geodetic Sciences, 24(1): 44-57, Jan-Mar, 2018 The reports of the geodetic stations from the IBGE database were used to obtain the geometric (h) and orthometric (H) heights of the validation stations.Using the height information, it was possible to obtain the geoid undulation (N) in each station by the equation. ( where N is the geoid height; h is the geometric height, and H is the orthometric height of the stations.Once calculated, the geoid undulations of the stations were compared with the values resulting from the interpolation methods.The validation was performed by calculating the RMSD (Eq.4) and the absolute maximum differences obtained with each interpolation. ( where    are the calculated geoid heights (Eq.03) and    are interpolated geoid heights obtained from MAPGEO2015.

Results
The experiment results are analyzed and discussed with the purpose of estimating the magnitude of the interpolation error on a geoid model to check whether there is a significant difference in the use of various interpolators, considering the geoid model of Brazil in its latest version, the MAPGEO2015.

Interpolators Analysis with Resample Grids
In this first experiment, the errors measurements of selected interpolators were performed for 10'x10 ', 15'x15', 20'x20' and 25'x25' grids and also an estimation of expected error for the MAPGEO2015 original 5'x5' grid.Figure 4 shows a profile from a side view in longitude of geoid height differences, contemplating all latitudes of the study area, obtained for each interpolator (row) and mesh grid (column).And Figure 5 shows the same differences between interpolated and expected height values from the original grid but in map form.It can be seen from Figure 4 that the interpolators show larger errors with a larger resolution.In addition, most part of the differences in the center of the study area range from 15 to 50 cm while larger differences can be observed at the edges of the study area (~ 1 m).
In the 10'x10' resample (first column of Figure 5) the interpolators show similar behavior since the differences are almost the same.However, while the mesh resample is described with a larger resolution, the bilinear interpolation shows expressive dissimilarity compared to RBF and spline interpolators.In general, the larger differences are located in the Andes Mountains and coastal surroundings of Brazil.The largest interpolated errors in the Andes occurred at the border of the mountains, which is coincident with the location where the geoid shows most variation.This geoid variability can be justified by the presence of large mountain ranges.In the coastal region of Brazil, the biggest differences are also related to geoid variability, being occasioned by ocean relief behavior.
Figure 5: Map differences of geoid height on validation points between the interpolated values (bilinear, spline and RBF) and those from the MAPGEO2015 grid, for each resample.
In the 10'x10' resample (first column of Figure 5) the interpolators show similar behavior since the differences are almost the same.However, while the mesh resample is described with a larger resolution, the bilinear interpolation shows expressive dissimilarity compared to RBF and spline interpolators.In general, the larger differences are located in the Andes Mountains and coastal surroundings of Brazil.The largest interpolated errors in the Andes occurred at the border of the mountains, which is coincident with the location where the geoid shows most variation.This geoid variability can be justified by the presence of large mountain ranges.In the coastal region of Brazil, the biggest differences are also related to geoid variability, being occasioned by ocean relief behavior.
Figure 6 shows the RMSE calculated for each interpolator and mesh grid, taking into account the number of grid points used for interpolations.The difference between interpolations gets smaller with the reduction of the mesh resolution and the increase of points used during the interpolation.,An RMSE of 2.5 cm is estimated for the three interpolators and approximately 230.000 points in the study area from the data extrapolation for a 5'x5' mesh.Thus, when using the MAPGEO2015 5'x5' grid, the bilinear interpolation will not show expressive difference with robust interpolators such as spline or RBF.Table 1 shows the RMSD between the interpolators for each resample in order to analyze their consistency.The RMSD of bilinear-spline and bilinear-RBF is ~ 13 cm in the 25'x25' grid and it reduces to ~ 4 cm in the 10'x10' grid.Also, a small deviation can be observed between the spline and the RBF interpolations, which remains around 3 cm and 2 cm.These values were expected as the neural network used a cubic basis function, with features similar to the cubic spline interpolation.

Comparing interpolators with mapgeo2015_v1 interpolator
In this experiment the comparison was made in order to evaluate the behavior of various interpolators compared to bilinear interpolation from mapgeo2015_v1, which is used by the geodetic community.Figure 7 shows the deviation obtained in the comparison.
Figure 7: Difference between MAPGEO2015 interpolation and bilinear, spline and RBF in a 2.5'x2.5'mesh.The larger differences can be observed in the same regions where the geoid model has more variability, i. e., the Andes Mountains and the coast of Brazil.In this experiment, RMSD was 1.11 cm, 1.44 cm and 1.54 cm respectively for bilinear, spline and RBF interpolation; and the maximum differences obtained were 20.47 cm, 36.75 cm and 32.51 cm.The bilinear interpolator thus showed more similarity with the MAPGEO2015 executable.In addition, the RMSD obtained is lower than the estimated error in the interpolations in section 4.1.
It is worth mentioning that IBGE informs the MAPGEO2015 executable carries out a bilinear interpolation which corroborates with the results of the bilinear interpolator, showing a lower RMSD.However, the 1.11 cm RMSD and the 20.47 maximum deviation were seen to be significant.We may therefore inquire about how the interpolator is implemented in the mapgeo2015_v1.Three bilinear interpolators were confronted with the MAPGEO2015 program.The first one was implemented by the authors according to Pedrini and Schwartz (2008), the second one comes from the functions already provided by Matlab and the third one comes from Python's libraries.In addition, the interpolator implemented by the authors obtained identical results to that provided by Matlab and Python.

Interpolator Analyses on Leveling Benchmarks
This experiment aimed to evaluate the interpolators in real measurable points.The geoid heights were calculated considering the geometric and orthometric heights available in 539 leveling benchmarks tracked by GNSS and provided by IBGE.The comparison of those geoid heights was made with the three interpolators in the original MAPGEO grid (5'x5').Table 2 shows the results of RMSD and maximum differences for each interpolator.
The lowest RMSD and maximum difference values were obtained by mapgeo2015_v1.to, The RMSD of bilinear, spline and RBF showed a difference of 0.03 cm, 0.01 cm and 0.01 cm respectively; and the maximum errors differ up to 1.61 cm, 2.46 cm and 2.41 cm.Therefore, the interpolators are similar for a centimetric level.(2015) show an RMSD of 17 cm and the maximum error of 49 cm, which is quite similar to the results of Table 2 for RMSD, but the maximum difference may be justified by outliers that affected our results because we used distinct stations in comparison to the MAPGEO2015 validation.
In this experiment, the lowest RMSD and maximum difference values were obtained by mapgeo2015_v1.The RMSD of bilinear, spline and RBF showed a difference of 0.26 cm, 0.11 cm and 0.14 cm respectively; and the maximum errors differ up to 0.39 cm, 0.46 cm and 0.41 cm.Therefore, the interpolators are similar at a centimeter level.The Figure 8 shows respectively the histogram of bilinear, spline, RBF and mapgeo2015_v1 interpolations.We performed statistical analyses in a search of an answer to the second question of this work "is there a significant difference in information extracted by users using different interpolations?".For this analyses it is assumed all samples are independent and normally distributed.We first applied the F-test to check whether the interpolators have equal (null hypothesis) or different (alternative hypothesis) variance.The largest difference, in terms of standard deviation values, was obtained between mapgeo2015_v1 and bilinear.Consequently, we presume if they show statistically equal variance other interpolator pairs will have the same behavior.Otherwise, all pair combination must be checked.The calculated F-test value for the interpolators was 0.974.Since the theoretical F-test values with 5% significance level (SL) are 0.844 and 1.184, for samples with same degree of freedom (DF) 538, we can affirm the interpolation variances are not significantly different.
In addition, we applied the Fisher method, also known as least significant difference (LSD) test, to check whether all pairs of interpolators arithmetic means are significantly different.This test consists of two stages.In the first stage, a multiple-comparison ANOVA F-test is performed to check whether there are two or more groups showing different means.The second stage can be applied only if this F-test does not reject the null hypothesis of equality (Hayter 1986).Following the Fisher method stage 1, the ANOVA test was carried out.The theoretical F-test (for 5% SL, 3 DF between-group and 2152 DF within-group) is 2.609 and the calculated F-test value is 0.021.Based on this statistical analyze, it is possible to affirm the interpolators arithmetic means are not significantly different.
Since the interpolators did not show a significant difference, only one interpolator was chosen to be shown in Figure 9, which represents the errors of MAPGEO2015 obtained using the bilinear interpolator in leveling benchmarks tracked by GNSS.(2015).Lower accordance between the geoid model and GNSS positioning occurred in the North region due to the water load of the Amazon Basin as well as the greater distances between that the local vertical data, located in the South of Brazil (Imbituba, SC), and the North stations.More study is therefore required to improve this accordance.

Conclusion
Two main questions were addressed in this work.The first refers to the error magnitude on the use of interpolation methods in regular grids that are available as a product of a geoid model in Bulletin of Geodetic Sciences, 24(1): 44-57, Jan-Mar, 2018 Brazil.The second inquires if there is a significant difference in using the three interpolators to extract information from a regular grid that resulted from a geoid model.
It was observed in the first experiment that typically the interpolators showed larger errors with lower spatial resolution.However, bigger differences are located where the geoid shows more variability, such as the Andes Mountains and the coast of Brazil.In these areas non-linear interpolators showed better results.Studies for quantification of errors in areas with high and low variability are strongly suggested.We estimated the error is 2.5 cm for the three selected interpolators applied in geoid height when using a grid 5'x5'.This error must be taken into consideration in the 18 cm accuracy of MAPGEO2015.Another interpolator may not show big difference when using a grid of 5'x5' for the study area.This methodology of estimation geoid interpolator error can be applied to other interpolation methods in order to compare our results.
In the second experiment the results showed unexpected differences between mapgeo2015_v1 and bilinear interpolation.These differences may indicate the MAPGEO2015 executable uses bilinear interpolation with additional modifications.The proposed methodology cannot be applied directly using the executable, however since in the first experiment the interpolators showed similar results we can conclude the mapgeo2015_v1 may shows same error magnitude.
The calculated statistical values showed high similarity in the third experiment.In order to verify whether these statistics are significantly different two statistical tests were performed, allowing the affirmation the selected interpolation methods are not significantly different to obtained results from the MAPGEO2015 interpolator.
In a general conclusion, the interpolator methods are consistent with geoid applications that require a centimetric level.However, the interpolators error needs to be taken into consideration in high accuracy application, such as in the definition of reference systems and geodynamics, where millimeter-lever is required.In these cases, further analysis should be carried out to define a more efficient interpolator and the proposed methodology can be used as an underlying base.

Figure 2 :
Figure 2: Resample example (10'x10') and choice of validation points.Another three grids were generated from the MAPGEO2015 5'x5' grid: 15'x15', 20'x20' and 25'x25'.As the resolution gets larger, more validation points from the original grid in 5'x5' are used as reference values in the calculation of errors.To assess the quality of the interpolation for new grids, the validation was performed through the analysis of the root mean square error (RMSE, Eq. 01) and the absolute maximum error obtained for each interpolation and new grid.
obtained by other interpolators and  2015  are mapgeo2015_v1 interpolated geoid heights.

Figure 3 :
Figure 3: Location of the 539 stations used as validation in the third experiment.

Figure 4 :
Figure 4: Plot differences of geoid height on validation points between the interpolated values (bilinear, spline and RBF) and those from the MAPGEO2015 grid, for each resample.

Figure 6 :
Figure 6: RMSE calculated for 25'x25', 20'x20', 15'x15' and 10'x10' resample and estimated to the 5'x5' MAPGEO2015 mesh.It is possible to verify from Figures 4 and 5 that maximum deviations reach a magnitude of just a few meters, despite an RMSE below 18 cm for interpolators.The maximum errors on the 10'x10', 15'x15', 20'x20' and 25'x25' grids were, respectively, 0.87 m, 1.63 m, 2.36 m and 2.87 m for bilinear interpolation; 1.32 m, 1.48 m, 1.50 m and 1.88 m for spline interpolation; and 0.73 m, 1.40 m, 1.45 m and 1,85 m for interpolation by RBF networks.It can be seen that the reduction of maximum errors is the result of the decrease of mesh distance.The RBF network was shown to be the most efficient interpolator in terms of RMSE and maximum difference.

Figure 9 :
Figure 9: Geoid height difference of MAPGEO2015 (5'x5') on leveling benchmarks obtained using the bilinear interpolation.The larger errors of MAPGEO2015 are located in the North of Brazil, which agrees with Blitzkow et al. (2016) and the Technical Report of IBGE (2015).Lower accordance between the geoid model and GNSS positioning occurred in the North region due to the water load of the Amazon Basin as well as the greater distances between that the local vertical data, located in the South of Brazil (Imbituba, SC), and the North stations.More study is therefore required to improve this accordance.

Table 2 :
RMSD and absolute maximum errors of the comparison of known geoid heights in benchmarks and interpolated heights with bilinear, spline, RBF and mapgeo2015_v1