Abstract:
This paper aimed to evaluate the use of Global Navigation Satellite System (GNSS) ellipsoidal heights in conjunction with height anomalies provided by Global Geopotential Model (GGM) XGM2019e, refined by Residual Terrain Modelling (RTM) technique, to obtain normal heights in Brazil, referred to the Imbituba Brazilian Vertical Datum (IBVD) and the International Height Reference System (IHRS). For this purpose, a local modelling approach has been analyzed in contrast to the national modeling one on the reference geopotential value. For this, a methodology based on geopotential space was adapted. In the local modeling, two study subregions were defined using the spatial clustering analysis of IBVD and GGM/RTM height anomalies differences outliers. The parameters have been estimated using three different configurations. In the parameters validation step, Root Mean Square Errors (RMSE) of the discrepancies between transformed and Brazilian official normal heights were calculated. In both subregions more accurate results have been obtained with the local modeling. In the SP1 subregion the accuracy increased tenfold (0.97m to 0.10m) and SP2 improved from 0.39m to 0.17m. For the linkage to the future realization of IHRS, the accuracy analysis was not possible. However, discrepancies between calculated normal heights and Brazilian official normal heights have been analyzed.
Keywords:
Normal heights; GNSS; Global Geopotential Model; Vertical datum linkage; IHRS
1. Introduction
In the wake of a more global view of the Earth’s intrinsic processes related to Geodesy, Drewes (2006Drewes, H. 2006. The changing objectives in geodetic research. Zeitschrift für Geodäsie, GeoInformation und Landmanagement, 131 (5), pp.17.) stated that the classic goals of Geodetic Sciences were added to others that place Geodesy at the center of the discussion on the challenge of living on a dynamic planet. This more recent interpretation of Geodesy’s role as an avantgarde science brought about the need not only for measurement, but also for an understanding of physical entities directly linked to the way the Earth’s surface is organized and modified over time. This more global view considers the whole as the Earth System, needing continuous and accurate monitoring over time.
In this context, initiatives promoted by the International Association of Geodesy (IAG) such as the creation of the Global Geodetic Observing System (GGOS) and the search for a global connection of geodetic observations, confirm this purpose. In terms of global unification, the state of the art also points to the definition of a Height Reference System (HRS), which was indicated, in terms of definition and implementation, by IAG Resolution No. 01 of 2015 (Ihde et al. 2017Ihde, J., Sánchez, L., Barzaghi, R., Drewes, H., Foerste, C., Gruber, T., Liebsch, G., Marti, U., Pail, R., and Sideris, M. 2017. Definition and proposed realization of the International Height Reference System (IHRS). Surveys in Geophysics, 38 (3), pp.549570. https://doi.org/10.1007/s1071201794093.
https://doi.org/10.1007/s107120179409...
). This is the IHRS and it is future realization International Height Reference Frame (IHRF).
Regionally, as a first step towards the realization of a unified vertical frame, according to de Freitas (2015de Freitas, S.R.C. 2015. SIRGASWGIII activities for unifying height systems in Latin America. Revista Cartográfica, (91). DOI:10.35424/rcarto.i91.452
https://doi.org/10.35424/rcarto.i91.452...
), the Working Group III of Geodetic Reference System for the Americas (SIRGAS), has operated in the area for more than two decades, seeking to standardize rules, procedures and activities among Latin American countries, aiming to meet international requirements.
In the case of Brazil, after the Readjustment of Altimetric Network with Geopotential Numbers, carried out in 2018 (REALT2018  Reajustamento da Rede Altimétrica com Números Geopotenciais in Portuguese), normal heights have been adopted for the High Precision Altimetric Network (HPAN) of the Brazilian Geodetic System (BGS). According to IBGE (2019), IBGE 2019. Reajustamento da Rede Altimétrica com Números Geopotenciais REALT2018 ENT#091;Readjustment of the Altimetric Network with Geopotential Numbers REALT2018ENT#093; ENT#091;onlineENT#093;. Available from: https://biblioteca.ibge.gov.br/visualizacao/livros/liv101666.pdf.
https://biblioteca.ibge.gov.br/visualiza...
this type of height was adopted for being the most adequate to modern concepts and scientific methods and it is in line with the resolutions established by the IAG with aim to define the IHRS. Additionally, the adjustment was carried out using geopotential numbers, in order to prepare the HPAN for future actions under the IHRS; and new spirit leveling and gravimetry observations, as well as correction of eventually detected inconsistencies have been included. However, the new normal heights remain referenced to the Brazilian Vertical Datums (BVD) of Imbituba and Santana, current in Brazil. These datums, in each case, were defined from a single value of the Mean Sea Level, which was classified with data collected in only one tide gauge station, therefore, of a local approach. Thus, given the impossibility of connecting the two leveling networks, it appears that there are two vertical datums in the country. Also, according to IBGE (2019), IBGE 2019. Reajustamento da Rede Altimétrica com Números Geopotenciais REALT2018 ENT#091;Readjustment of the Altimetric Network with Geopotential Numbers REALT2018ENT#093; ENT#091;onlineENT#093;. Available from: https://biblioteca.ibge.gov.br/visualizacao/livros/liv101666.pdf.
https://biblioteca.ibge.gov.br/visualiza...
it is targeted that the future realization of the vertical component of the BGS, in a few years, will already be related to the IHRS/IHRF.
To obtain the normal height value referred to any of the BVD, at a given point in the Brazilian territory, it is possible to carry out spirit levelling from a benchmark (RN  “Referência de Nível” in Portuguese) station and interpolate the real gravity value from Brazilian gravimetric network data. However, with the expressive technological development experienced in recent decades, new methods for obtaining altimetric data have gained projection, supported by the improvement of related observations and techniques. In this context innumerable researches have proposed the use of GNSS ellipsoidal heights in conjunction with data provided by Global Geopotential Model (GGM), which can be refined by Residual Terrain Modelling (RTM) technique using Digital Terrain Models (DTM) data, to obtain physical heights. For example, Hirt, Featherstone, and Marti (2010Hirt, C., Featherstone, W.E., and Marti, U. 2010b. Combining EGM2008 and SRTM/DTM2006.0 residual terrain model data to improve quasigeoid computations in mountainous areas devoid of gravity data. Journal of Geodesy , 84 (9), pp.557567. https://doi.org/10.1007/s0019001003951.
https://doi.org/10.1007/s001900100395...
)like EGM2008, is not capable of representing the highfrequency components of Earth’s gravity field. This is known as the omission error. In mountainous terrain, omission errors in EGM2008, even when expanded to degree 2,190, may reach amplitudes of 10 cm and more for height anomalies. The present paper proposes the utilisation of highresolution residual terrain model (RTM combined residual terrain model data to improve quasigeoid computations in mountainous areas devoid of gravity data, and Castro, Jamur, and de Freitas (2013)Castro, H.D.M., Jamur, K.P., and Freitas, S.R.C. de 2012. Exploration of the residual topography effect on some functionals of the gravity field. Revista Brasileira de Geofísica, 30 (4), pp.431. https://doi.org/10.22564/rbgf.v30i4.230.
https://doi.org/10.22564/rbgf.v30i4.230...
considering their distribution and variations in density, have a measurable influence on Functionals of the Gravity Field (FGF explored of the residual topography effect on some functionals of the gravity field.
In the context of normal heights obtained from GNSS/GGM/RTM integration, there is the advantage of avoiding the costly and timeconsuming spirit levelling procedure from some benchmark, especially in regions with lower station density. However, there are two issues to consider. The first is the normal height accuracy. While GNSS ellipsoidal height accuracy can be obtained at millimeter level or better, the higher order modern GGM height anomalies, even if refined using DTM, are accurate in the order of magnitude of centimeter or decimeter. The second is the question of the datum, associated with the GGM/RTM zeroheight surface, which has a global definition approach.
For the first question, the solution is to limit the use of such a methodology to certain types of applications in which the orders of magnitude of the accuracies may be permissible. In the second question, the solution is, somehow, to estimate a transformation parameter or set of transformation parameters between Height Reference Frames (HRFs), and, therefore, establish the linkage of the normal height obtained from GNSS/GGM/RTM to one of the BVD. For this purpose, some methods have been proposed over the last years. For example, the methodology based on geopotential space as indicated in Ihde et al. (2017Ihde, J., Sánchez, L., Barzaghi, R., Drewes, H., Foerste, C., Gruber, T., Liebsch, G., Marti, U., Pail, R., and Sideris, M. 2017. Definition and proposed realization of the International Height Reference System (IHRS). Surveys in Geophysics, 38 (3), pp.549570. https://doi.org/10.1007/s1071201794093.
https://doi.org/10.1007/s107120179409...
) and the geodetic boundary value problem (GBVP) method proposed initially by Rummel and Teunissen (1988Rummel, R. and Teunissen, P. 1988. Height datum definition, height datum connection and the role of the geodetic boundary value problem. Bulletin Géodésique , 62 (4), pp.477498. https://doi.org/10.1007/BF02520239.
https://doi.org/10.1007/BF02520239...
).
Given the above, this research aimed to evaluate the use of GNSS ellipsoidal height in conjunction with height anomalies provided by a GGM, refined by RTM technique, for the determination of normal heights in Brazil, referred to the IBVD and to the future IHRS. For the linkage, an adaptation of the method based on geopotential space indicated in Ihde et al. (2017Ihde, J., Sánchez, L., Barzaghi, R., Drewes, H., Foerste, C., Gruber, T., Liebsch, G., Marti, U., Pail, R., and Sideris, M. 2017. Definition and proposed realization of the International Height Reference System (IHRS). Surveys in Geophysics, 38 (3), pp.549570. https://doi.org/10.1007/s1071201794093.
https://doi.org/10.1007/s107120179409...
) has been used. In the case of linkage at the IBVD, local estimated parameters and the one obtained by Sánchez and Sideris (2017Sánchez, L. and Sideris, M.G. 2017. Vertical datum unification for the International Height Reference System (IHRS). Geophysical Journal International, pp.ggx025. https://doi.org/10.1093/gji/ggx025.
https://doi.org/10.1093/gji/ggx025...
) at national level have been analyzed. In the second case, the IHRS geopotential value at the geoid W _{0} has been used for the linkage.
It is noteworthy that, according to IBGE (2019)IBGE 2019. Reajustamento da Rede Altimétrica com Números Geopotenciais REALT2018 ENT#091;Readjustment of the Altimetric Network with Geopotential Numbers REALT2018ENT#093; ENT#091;onlineENT#093;. Available from: https://biblioteca.ibge.gov.br/visualizacao/livros/liv101666.pdf.
https://biblioteca.ibge.gov.br/visualiza...
, in Brazil, six stations belonging to the Brazilian Network for Continuous Monitoring of GNSS Systems (PPTE, IMBT, CUIB, MABA, CEFT and BRAZ) have been preselected to integrate the IHRF, being the first in the study region of this research and within the set of stations used in this research. In this sense, studies are necessary about the densification of gravimetric and the linkage of vertical datums. This paper aims to contribute to the second topic.
2. Methods
2.1 Data and Study Region
Normal heights obtained from GNSS/GGM/RTM are referred to a zeroheight surface associated with the zerolevel geopotential W _{0} value of the GGM/RTM solution used. In order to link such heights at the BVD, as mentioned, a transformation between HRFs must be performed (Rummel and Teunissen 1988Rummel, R. and Teunissen, P. 1988. Height datum definition, height datum connection and the role of the geodetic boundary value problem. Bulletin Géodésique , 62 (4), pp.477498. https://doi.org/10.1007/BF02520239.
https://doi.org/10.1007/BF02520239...
). As will be presented in item 2.2, in the methodology indicated in Ihde et al. (2017Ihde, J., Sánchez, L., Barzaghi, R., Drewes, H., Foerste, C., Gruber, T., Liebsch, G., Marti, U., Pail, R., and Sideris, M. 2017. Definition and proposed realization of the International Height Reference System (IHRS). Surveys in Geophysics, 38 (3), pp.549570. https://doi.org/10.1007/s1071201794093.
https://doi.org/10.1007/s107120179409...
), the zerolevel geopotential value in a local datum is estimated as parameter. For this purpose, the differences between BGS and GGM/RTM height anomalies are used, considering a set of points. A given BGS height anomaly, at a given RN/Global Positioning System (GPS) connection station, can be obtained subtracting the BGS ellipsoidal height and the BGS normal height, referred to one of the BVDs. On the other hand, a GGM/RTM height anomaly can be obtained from the International Centre for Global Earth Models (ICGEM, http://icgem.gfzpotsdam.de/, last access: 6 August 2021) (Ince et al. 2019Ince, E.S., Barthelmes, F., Reißland, S., Elger, K., Förste, C., Flechtner, F., and Schuh, H. 2019. ICGEM  15 years of successful collection and distribution of global gravitational models, associated services, and future plans. Earth System Science Data, 11 (2), pp.647674. https://doi.org/10.5194/essd116472019.
https://doi.org/10.5194/essd116472019...
) and calculation of the RTM effect using one or more DTMs. It should be noted that currently there are several options for obtaining global DTMs for free on the web.
Considering all available RN/GPS connection stations in the Brazilian territory at the time of the study, Sánchez and Sideris (2017Sánchez, L. and Sideris, M.G. 2017. Vertical datum unification for the International Height Reference System (IHRS). Geophysical Journal International, pp.ggx025. https://doi.org/10.1093/gji/ggx025.
https://doi.org/10.1093/gji/ggx025...
) obtained the mean (W _{0} ) values for IBVD and Santana Brazilian Vertical Datum (SBVD), whose values were 62,636,849.61±0.18 m^{2}/s^{2} and 62,636,852.93 ± 1.51 m^{2}/s^{2} respectively. From these values, it is possible to carry out the transformation associated with each datum for the entire Brazilian territory. However, in this research, with the purpose of investigating the problem of transformation at the local level, the value of (W _{0} ) was calculated considering a smaller study region, using the subset of available RN/GPS data. The defined study region comprises the state of São Paulo, an administrative region with 102 RN/GPS stations distributed throughout the entire territory. In the Figure 1 it is shown the spatial distribution of the RN/GPS stations.
Spatial distribution of RN/GPS connection stations in São Paulo state and the 4X RN station in Imbituba.
All benchmarks in this region have normal height referred to IBVD, materialized through the station code “4X”, which is the fundamental benchmark of the IBVD. For the experiments, all 102 stations have been used, being a part separated to estimate the transformation parameter (zerolevel geopotential value in local datum), and the complementary part for it is evaluation, as will be seen below. All data were obtained free of charge from the BGS Geodetic Database (BGD).
For linking the normal heights determined by GNSS/GGM/RTM to the future IHRS, a transformation between HRFs is also required. However, by using the same transformation strategy for the previous case, an estimate of W _{0} is not necessary. This is because the IHRS zerolevel geopotential value can be directly used. In this case, only the GGM/RTM height anomalies have been used to calculate the geopotential value at the point. It is noteworthy that, for both linkage cases, to compare the values of normal heights obtained from GNSS/GGM/RTM, the same RN/GPS stations set in the study region has been used.
During the process of choosing the GGM, the study developed by Nicacio, Dalazoana, and Freitas (2018Nicacio, E., Dalazoana, R., and Freitas, S.R.C. de 2018. Evaluation of recent combined global geopotential models in Brazil. Journal of Geodetic Science, 8 (1), pp.7282. https://doi.org/10.1515/jogs20180008.
https://doi.org/10.1515/jogs20180008...
)was taken as reference, which verified the employability of some of the most recent GGMs for the extraction of normalorthometric heights in Brazilian territory. It is noteworthy that, this research was carried out in 2017, when normalorthometric heights were in use. According to IBGE (2019)IBGE 2019. Reajustamento da Rede Altimétrica com Números Geopotenciais REALT2018 ENT#091;Readjustment of the Altimetric Network with Geopotential Numbers REALT2018ENT#093; ENT#091;onlineENT#093;. Available from: https://biblioteca.ibge.gov.br/visualizacao/livros/liv101666.pdf.
https://biblioteca.ibge.gov.br/visualiza...
, the differences between previous and current HPAN heights in the study region vary almost entirely between 0.01 m and 0.01 m. Therefore, the results obtained by Nicacio, Dalazoana, and Freitas (2018)Nicacio, E., Dalazoana, R., and Freitas, S.R.C. de 2018. Evaluation of recent combined global geopotential models in Brazil. Journal of Geodetic Science, 8 (1), pp.7282. https://doi.org/10.1515/jogs20180008.
https://doi.org/10.1515/jogs20180008...
remain valid after REALT2018. As GGM XGM2016 (Pail et al. 2018Pail, R., Fecher, T., Barnes, D., Factor, J.F., Holmes, S.A., Gruber, T., and Zingerle, P. 2018. Short note: the experimental geopotential model XGM2016. Journal of Geodesy , 92 (4), pp.443451. https://doi.org/10.1007/s0019001710706.
https://doi.org/10.1007/s001900171070...
) had the best performance and has been updated to XGM2019 (Zingerle et al. 2020Zingerle, P., Pail, R., Gruber, T., and Oikonomidou, X. 2020. The combined global gravity field model XGM2019e. Journal of Geodesy , 94 (7), pp.66. https://doi.org/10.1007/s00190020013980.
https://doi.org/10.1007/s0019002001398...
), this has been selected for the development of this research. It is noteworthy that, this GGM is more accurate than its predecessor, by adding topographically derived gravimetric information over the terrain. The ICGEM calculation tool was used with data from the GGM XGM2019e up to degree and order 2190. The zerodegree term was selected. The reference ellipsoid used was the GRS80 and the mean tide system was selected.
To minimize the omission error in the GGM data and in order to increase the spectral resolution, the topography data can be used. This information is extracted from DTMs. As mentioned, this procedure proposes a refinement in the GGM height anomalies. For this, the RTM technique has been applied, by using the prisms method, with a radius of approximately 200 km, in agreement with Hirt, Featherstone, and Marti (2010Hirt, C., Featherstone, W.E., and Marti, U. 2010b. Combining EGM2008 and SRTM/DTM2006.0 residual terrain model data to improve quasigeoid computations in mountainous areas devoid of gravity data. Journal of Geodesy , 84 (9), pp.557567. https://doi.org/10.1007/s0019001003951.
https://doi.org/10.1007/s001900100395...
). It should be pointed out that harmonic corrections have been not taken in account. For a given point P, after calculating the residual height anomaly (${\zeta}_{P}^{RTM}$) by the RTM technique using a grid of rectangular prisms and planar approximation (MacMillan 1930MacMillan WD (1930): Theoretical Mechanics, vol 2: the Theory of the potential. McGrawHill, New York (reprinted by Dover Publications, New York 1958).) (Heck and Seitz 2007Heck, B. and Seitz, K. 2007. A comparison of the tesseroid, prism and pointmass approaches for mass reductions in gravity field modelling. Journal of Geodesy, 81 (2), pp.121136. https://doi.org/10.1007/s0019000600940.
https://doi.org/10.1007/s001900060094...
)the topographic surface of the Earth is often divided with respect to geographical or mapgrid lines, and the topographic heights are averaged over the respective grid elements. The bodies bounded by surfaces of constant (ellipsoidal, the refinement of the GGM height anomaly (${\zeta}_{P}^{GGM}$) can be done by:
where ${\zeta}_{P}^{GGM/RTM}$ is the refined GGM height anomaly.
The selected DTMs were those that had the best cost/benefit ratio based on factors such as accuracy, computational cost and resolution in previously tests performed by the authors. It was chosen the ETOPO1 (Amante and Eakins 2009Amante, C and Eakins, B W 2009. ETOPO1 Global Relief Model converted to PanMap layer format. application/zip. https://doi.org/10.1594/PANGAEA.769615.
https://doi.org/10.1594/PANGAEA.769615...
) model with 1’ spatial resolution and the SRTM15_PLUS (Tozer et al. 2019Tozer, B., Sandwell, D.T., Smith, W.H.F., Olson, C., Beale, J.R., and Wessel, P. 2019. Global Bathymetry and Topography at 15 Arc Sec: SRTM15+. Earth and Space Science, 6 (10), pp.18471864. https://doi.org/10.1029/2019EA000658.
https://doi.org/10.1029/2019EA000658...
) 15” spatial resolution, obtained, respectively, from the ICGEM and the ERDDAP data server of the National Oceanic and Atmospheric Administration. For the application of the RTM technique, the ETOPO1 model has been used as reference surface, while SRTM15_PLUS model has been used as detailed surface. Both models cover oceanic areas derived from bathymetric data.
In the data integration, the differences between permanent tide effect have been taken in account. While the BGS ellipsoidal heights are referred to tide free system, the BGS normal heights and DTM heights are referred to mean tide system and zero tide system, respectively. For the compatibilization, the mean tide system has been adopted and the Rapp (1989Rapp, R.H. 1989. The treatment of permanent tidal effects in the analysis of satellite altimeter data for sea surface topography. Manuscripta Geodaetica.) formulation has been used, as recommended by the IAG in the conventions for the adoption of the IHRS (Ihde et al. 2017Ihde, J., Sánchez, L., Barzaghi, R., Drewes, H., Foerste, C., Gruber, T., Liebsch, G., Marti, U., Pail, R., and Sideris, M. 2017. Definition and proposed realization of the International Height Reference System (IHRS). Surveys in Geophysics, 38 (3), pp.549570. https://doi.org/10.1007/s1071201794093.
https://doi.org/10.1007/s107120179409...
).
2.2 Strategies for normal heights linkage
After defining the study region, obtaining and performing the compatibilization of all data, the linkage of the GNSS/GGM/RTM normal heights to IBVD and to the future IHRS were carried out. For this, the methodology presented by Ihde et al. (2017Ihde, J., Sánchez, L., Barzaghi, R., Drewes, H., Foerste, C., Gruber, T., Liebsch, G., Marti, U., Pail, R., and Sideris, M. 2017. Definition and proposed realization of the International Height Reference System (IHRS). Surveys in Geophysics, 38 (3), pp.549570. https://doi.org/10.1007/s1071201794093.
https://doi.org/10.1007/s107120179409...
) to define gravitational potential values for regional datums and their respective differences from a global reference has been adapted. This methodology is based on obtaining zerolevel geopotential value (W _{0k} ) of a Local Vertical Datum (LVD) k, and it is difference (∆W _{0k} ) in relation to the zerolevel geopotential value (W _{0} ) of a global HRS. Using the HRS unification principle using geopotential source values (Ihde et al. 2017Ihde, J., Sánchez, L., Barzaghi, R., Drewes, H., Foerste, C., Gruber, T., Liebsch, G., Marti, U., Pail, R., and Sideris, M. 2017. Definition and proposed realization of the International Height Reference System (IHRS). Surveys in Geophysics, 38 (3), pp.549570. https://doi.org/10.1007/s1071201794093.
https://doi.org/10.1007/s107120179409...
), quantities related to a LVD and a LVD can be obtained. The proposed adaptation in this research is to use the value of directly W _{0k} to obtain the normal height already referenced to LVD.
As previously mentioned, for linkage to IBVD, besides analyzing the value of ${W}_{0}^{IBVD\_SS}$ obtained by Sánchez and Sideris (2017Sánchez, L. and Sideris, M.G. 2017. Vertical datum unification for the International Height Reference System (IHRS). Geophysical Journal International, pp.ggx025. https://doi.org/10.1093/gji/ggx025.
https://doi.org/10.1093/gji/ggx025...
) for the entire Brazilian region, a more adjusted local model for the study region was also searched. Thus, it was necessary to obtain the mean value of ${\stackrel{}{{W}_{0}}}^{IBVD}$ from the height anomaly data subset.
In this context, in a given P point, BGS RN/GPS station, the anomalous potential can be obtained from the GGM (${T}_{P}^{GGM}$), and the h _{P} ellipsoidal height can be obtained from BGS/BDG, referenced to SIRGAS2000 and transformed from tidefree system to the mean tide system. With these quantities available and knowing the zerolevel spheropotential (U _{0} ) of the adopted Reference Ellipsoid (GRS80) (Moritz 1980Moritz, H. 1980. Geodetic reference system 1980. Bulletin Géodésique, 54 (3), pp.395405. https://doi.org/10.1007/BF02521480.
https://doi.org/10.1007/BF02521480...
), the geopotential value for point P (W _{P} ) can be calculated by:
or yet:
where U _{P} is the spheropotential value for point P; $\frac{\partial {U}_{0}}{\partial h}$ is the normal gradient of U _{0} ; γ _{P} is the magnitude of the acceleration of normal gravity at point P; and ${\zeta}_{P}^{GGM}$ is the height anomaly obtained via GGM or GGM/RTM (${\zeta}_{P}^{GGM/RTM}$) for point P. In this research, the GGM refinement has been performed, according to Equation (1).
On the other hand, the same value of W _{P} is also related to the ${W}_{0}^{IBVD}$ value and to the geopotential number value of point P referred to IBVD (${C}_{P}^{IBVD}$) by:
From the combination of Equations (2) and (4), a value of ${W}_{0}^{IBVD}$ can be obtained for point P by:
or yet:
where ${H}_{P}^{NIBVD}$ is the normal height referred to the IBVD at point P, obtained from BGS/BGD; and ${\zeta}_{P}^{BGS}$ is BGS the height anomaly.
Knowing the values of ${W}_{0}^{IBVD,P}$ for a set of n points, the mean value of the geopotentials can be assumed as the local benchmark:
Finally, the normal height obtained from GNSS/GGM/RTM, at a given point Q, referred to IBVD, is:
where W _{Q} can be calculated with the form of Equation (2) with the data of point Q. It is noteworthy that, in this strategy, the quantity ${\stackrel{}{{W}_{0}}}^{IBVD}$ is considered as a transformation parameter between the HRFs.
In the second approach of this research, for the case of linkage to the future IHRS, the same functional mathematical modeling has been used. The difference, however, is that there is no $\stackrel{}{{W}_{0}}$ parameter to be estimated, since IHRS zerolevel geopotential value (${W}_{0}^{IHRS}$) has already been defined by IAG Resolution 01/2015. In this case, the normal height can be obtained by:
2.3 Validation of linkage and detection of outliers
In the context of linkage to IBVD, in order to estimate the local parameter ${\stackrel{}{{W}_{0}}}^{IBVD}$, from Equation (8), a BGS height anomalies dataset is required. Subsequently, after estimating the parameter, it must be evaluated, as well as the ${\stackrel{}{{W}_{0}}}^{IBVD\_SS}$ parameter calculated by Sánchez and Sideris (2017Sánchez, L. and Sideris, M.G. 2017. Vertical datum unification for the International Height Reference System (IHRS). Geophysical Journal International, pp.ggx025. https://doi.org/10.1093/gji/ggx025.
https://doi.org/10.1093/gji/ggx025...
).
Therefore, from the 102 RN/GPS stations in the study area, 68 stations have been separated to estimate the parameters, being called reference points (r); and 34 stations have been separated to evaluation of the parameter, being called check points (t). The ratio of 2/3 and 1/3 has been defined empirically. The designation of the points in each onethird has been made looking for a homogeneous spatial distribution throughout the study region, but without a quantitative criterion.
According to the three onethirds procedure, an initial onethird configuration of points was called Alpha (A). Subsequently, two other adopted onethirds were named Bravo (B) and Charlie (C), using exchanges between the reference and check points of the onethird part A. In this way, all stations were, in some of the three onethirds, considered reference and check points. Thus, the question to test the use of RN/GPS stations for both purposes, allowing transformation parameters in each onethird to be obtained, validated from all stations and compared.
The evaluation of the parameter ${\stackrel{}{{W}_{0}}}^{IBVD\_SS}$ in the national approach transformation, as well as the local parameter ${\stackrel{}{{W}_{0}}}^{IBVD}$, estimated in each onethird, was based on the RMSE calculation of the discrepancy’s dataset between normal height values:
where ${H}_{{Q}_{i}}^{Ntrf\_IBVD}$ is the ith normal height obtained from GNSS/GGM/RTM referred to IBVD at a given Q check point; ${H}_{{Q}_{i}}^{NIBVD}$ is the ith reference BGS normal height (RN) at a given Q check point; and t is the number of check points.
However, in the context of the local modelling, before calculating it in each onethird, the presence of outliers in the dataset of discrepancies values ${\zeta}_{P}^{BGS}{\zeta}_{P}^{GGM/RTM}$ was verified. First, in the detection step, a JarqueBera (JB) normality test has been performed (Miot 2017Miot, H.A. 2017. Avaliação da normalidade dos dados em estudos clínicos e experimentais ENT#091;Assessment of data normality in clinical and experimental studiesENT#093;. Jornal Vascular Brasileiro, 16 (2), pp.8891. https://doi.org/10.1590/16775449.041117.
https://doi.org/10.1590/16775449.041117...
), with a 95% confidence level.
After the test, with the identification of outlier points, in the adaptation process, the hypothesis of spatial correlation was verified, that is, if the outliers points were spatially grouped in a subregion. Therefore, a spatial data clustering analysis technique has been used, based on the standard deviation value. This technique considered that spatially correlated outliers can express a different local bias from that verified in the rest of the data in the sample set, rather than simply the presence of gross error. Thus, for the study region, there was the possibility of carrying out of linkage studies for one or more areas, called study subregions, with different local parameters.
In the case of linkage to the IHRS, there is no ${\stackrel{}{{W}_{0}}}^{IBVD\_SS}$ parameter to be estimated. To calculate the discrepancies were used only the check points. Furthermore, in the case of linkage to the future IHRS, there is no possibility of calculating the accuracy since there are no values to be taken as a reference.
3. Results and Discussion
3.1 Definition of study subregions and onethird parts for the case of local modeling
For the case of local modeling, after obtaining the discrepancies between the height anomalies ${\zeta}_{P}^{BGS}$ and ${\zeta}_{P}^{GGM}$ , calculated for the 102 GPS/RN connection stations in the study region, a histogram was generated, which is shown SP in Figure 2.
As can be seen, most of the discrepancies were located between 0 and 0.4 m, with a mean equal to 0.259 m and standard deviation equal to 0.343 m. However, occurrences of values above 0.7 m can also be verified, moving the histogram away from a normal distribution.
When analyzing this question, the JB normality test was applied to the dataset of the study region. The test result showed an index higher than the defined critical value for 95% confidence level.
The presence of outliers was verified using the Data Snooping method. After the process, the points considered outliers were those with the greatest discrepancy between the calculated height anomalies, located in the right part of the histogram. After applying the spatial data clustering analysis, it was verified that the outlier points were all clustered in the southeast area of the study region, indicating a spatial correlation. This corroborates the result seen in the right part of the histogram. Since the subregion of the outlier points is in a coastal and mountainous area, the difference in the systematic component may be related GGM/RTM inaccuracies. The outlier points subregion was named SP2 and the other SP1, without excluding any of the points in this classification.
After ranking the points in each subregion, the normality of the associated datasets was verified. The JB normality test was applied to the dataset of the SP1 study subregion and the ShapiroWilk normality test was applied to the dataset of the SP2 study subregion, as it is a sample of less than 50 observations. In both cases, the datasets were within the confidence level and can be considered normally distributed, although in SP2 the sample is small for conclusions and the test results may not be significant.
The histograms of discrepancies between the height anomalies ${\zeta}_{P}^{BGS}$ and ${\zeta}_{P}^{GGM/RTM}$ in the SP1 and SP2 study subregions are shown in Figure 2. Through these results, the idea of estimating different local parameters in each subregion has been adopted.
Histogram of discrepancies (m) between height anomalies ${\zeta}_{P}^{BGS}$ and ${\zeta}_{P}^{GGM}$ for region SP and subregions SP1 and SP2.
Regarding the question of evaluating the estimated parameters, the procedure of onethirds was applied with the stations available in each subregion, alternately Alpha, Bravo and Charlie, as mentioned in item 2.3. The criterion used was the search for three groups with points distributed evenly in space and in discrepancies. In Figure 3, the check points in each onethird of the subregion are presented.
3.2 Analysis of the linkage results
In the case of local modeling, using the reference points defined for each onethird in each subregion, it was possible to determine the ${\stackrel{}{{W}_{0}}}^{IBVD}$ local transformation parameters, according to Equation (8). The calculated values are shown in Table 1, as well as their respective standard deviations. For comparison purposes, the value of ${\stackrel{}{{W}_{0}}}^{IBVD\_SS}$ obtained by Sánchez and Sideris (2017Sánchez, L. and Sideris, M.G. 2017. Vertical datum unification for the International Height Reference System (IHRS). Geophysical Journal International, pp.ggx025. https://doi.org/10.1093/gji/ggx025.
https://doi.org/10.1093/gji/ggx025...
) is also presented.
${\stackrel{}{{W}_{0}}}^{IBVD\_SS}$ and ${\stackrel{}{{W}_{0}}}^{IBVD}$ local transformation parameters and standard deviations.
The estimated parameters values indicate satisfactory agreement between the defined onethirds. While in subregion SP1 the discrepancy values are in the order of 10^{2} m^{2}/s^{2} or ~10^{3} m, in subregion SP2, in the worst case (Alpha/Charlie) the value is 0.14 m. This indicates that the local parameters can be estimated regardless of the location of the points, especially in the SP1 subregion, which has the greatest number of points. Regarding the difference for the value of ${\stackrel{}{{W}_{0}}}^{IBVD\_SS}$, the mean for subregion SP1 was 9.409 m^{2}/s^{2} or ~1 m and for subregion SP2 it was 3.410 m^{2}/s^{2} or ~ 0.35 m, indicating differences that can be significant for some applications, mainly in the SP1 subregion. The smaller difference between the values of ${\stackrel{}{{W}_{0}}}^{IBVD}$ and ${\stackrel{}{{W}_{0}}}^{IBVD\_SS}$ in SP2 indicates greater adherence and better applicability of the national parameter in this subregion.
After estimating the transformation parameters, using the r reference points in each onethird of the study subregion, the determination of normal heights ${H}_{Q}^{Ntrf\_IBDV}$ for the t check points have been carried out, according to Equation (9), for the local case. Then, the sets of differences ${\delta}_{t}^{trf\_IBVD}$ between the heights ${H}_{Q}^{Ntrf\_IBDV}$ and ${H}_{Q}^{NIBVD}$ provided by BGS have been obtained.
Complementarily, the differences ${\delta}_{t}^{trf\_IBVD\_SS}$ have been also calculated for the heights ${H}_{Q}^{Ntrf\_IBVD\_SS}$ derived from ${\stackrel{}{{W}_{0}}}^{IBVD\_SS}$, and the differences ${\delta}_{t}^{trf\_IBVD}$ for the heights ${\delta}_{t}^{trf\_IBVD}$ derived from ${\stackrel{}{{W}_{0}}}^{IBVD\_SS}$. Thus, for standardization purposes, for each one of the same t check points of each onethird it was calculated:
In Figures S1 and S2 (supplementary files) the differences for all points are presented, by onethirds, in the two study subregions.
For the linkage to IBVD, as can be seen, for both subregions, in all onethirds, the differences ${\delta}_{}^{trf\_IBVD}$ were smaller, with occurrence balanced of positive and negative values, leading the mean of the discrepancies to close to zero. This result indicates the possible presence of only random errors or random errors together with systematic errors of low magnitude. On the other hand, the differences ${\delta}_{}^{trf\_IBVD\_SS}$ were higher, with only positive values, and mean value around 0.95 m for the SP1 subregion and 0.35 m for the SP2 subregion. This result indicates a certain systematic bias in the results, as expected from the results in Table 1 with equal displacement values. Consequently, it is possible to verify that there was a better adjust of the data in the case of local modeling. In general, it is also possible to verify that the choice of onethirds did not change the results. The differences in this topic did not show any spatial correlation with the topography.
In the case of linkage to the future IHRS, in subregion SP1 the difference values ${\delta}_{}^{trf\_IHRS}$ were all positive with a mean value around 0.55 m. On the other hand, in subregion SP2 with occurrence balanced of positive and negative values, bringing the mean of the differences to close to zero. In this case, while in subregion SP1 there is a certain systematic bias, in subregion SP2 there is a possible presence of random errors or random errors with systematic errors of low magnitude.
In order to statistically verify whether systematic bias are significant or not, at a confidence level of 95%, statistical tests with Student’s t distribution (t £ 30) were performed in the two subregions, using the three onethirds, with the following hypotheses:
where ${\mu}_{\delta}$ is the mean of a given dataset of differences for a given onethird part/subregion. Theoretical statistics for subregions SP1 and SP2 are ±2.05 and ±3.18, respectively. The values of the calculated t student statistics are shown in Table 2.
As can be seen, the differences ${\delta}_{}^{trf\_IBVD\_SS}$, in all onethirds, in both subregions, cannot be considered free of significant systematic bias at a 95% confidence level. On the other hand, differences ${\delta}_{}^{trf\_IHRS}$ cannot be considered free of significant systematic biases at a 95% confidence level only in the SP1 subregion.
As presented in item 2.3, to obtain the accuracy of the transformation in each onethird of each subregion, the RMSE of each dataset of differences ${\delta}_{}^{trf\_IBVD}$ and ${\delta}_{}^{trf\_IBVD\_SS}$ has been calculated. Additionally, the mean RMSE have been also calculated. The results are shown in Table 3. As mentioned, in the case of linkage to the future IHRS, it is not possible to calculate the accuracy, since there are no height values to be taken as reference.
${\mathit{R}\mathit{M}\mathit{S}\mathit{E}}_{\u2206{\mathit{H}}^{\mathit{N}}}$ of discrepancies ${\delta}_{}^{trf\_IBVD}$ e ${\delta}_{}^{trf\_IBVD\_SS}$ (m).
As can be seen, the agreement in the RMSE values between the onethirds is predominantly at the level of mm, with a maximum difference of the order of cm.
In subregion SP1, the accuracies were practically ten times better when using the local modeling approach, regardless of the onethird part of reference points used. The mean difference was 0.90 m. The values obtained with the local parameter were approximately 0.10 m.
On the other hand, in subregion SP2, the differences in accuracy between the local and national modeling approach were smaller, in the order of 0.22 m, with the local modeling presenting more accurate values. In general, the values were less accurate in relation to the SP1 subregion. The more accurate values with the national modeling approach in subregion SP2 in relation to subregion SP1 is a direct consequence of the greater proximity between the values of ${\stackrel{}{{W}_{0}}}^{IBVD}$ and ${\stackrel{}{{W}_{0}}}^{IBVD\_SS}$ in this sub region, as indicated above with the data in Table 1.
4. Conclusions and outlook
The use of GNSS positioning integrated with GGM data, refined by DTMs, to determine normal heights in Brazil, referred to a LVD such as IBVD or IHRS, configures an alternative to the traditional spirit leveling method. However, there are two issues to be considered, the lower accuracy and the difference in terms of vertical datum. While the first question is solved by limiting the use of such methodology to certain types of applications, in the second, the solution is to carry out the transformation between HRFs.
Considering the second question, the objective of this work was to evaluate the use of GNSS ellipsoidal heights, in conjunction with height anomalies provided by GGM, refined by RTM technique, to obtain normal heights in Brazil, referred to the IBVD and future IHRS. In the case of linkage to IBVD, local and national modelling approach have been analyzed using the parameter obtained by Sánchez and Sideris (2017Sánchez, L. and Sideris, M.G. 2017. Vertical datum unification for the International Height Reference System (IHRS). Geophysical Journal International, pp.ggx025. https://doi.org/10.1093/gji/ggx025.
https://doi.org/10.1093/gji/ggx025...
). The study region was the state of São Paulo.
For the linkage, the methodology indicated in Ihde et al. (2017Ihde, J., Sánchez, L., Barzaghi, R., Drewes, H., Foerste, C., Gruber, T., Liebsch, G., Marti, U., Pail, R., and Sideris, M. 2017. Definition and proposed realization of the International Height Reference System (IHRS). Surveys in Geophysics, 38 (3), pp.549570. https://doi.org/10.1007/s1071201794093.
https://doi.org/10.1007/s107120179409...
) was adapted. The adaptation was made to directly provide the normal height values already transformed. In the case of linkage to IBVD with local modeling, the presence of outliers in the dataset of differences between the BGS and GGM/RTM height anomalies was firstly investigated. After the spatial clustering analysis, the spatial correlation in the outliers was detected and two subregions were defined for modeling.
Thus, in the linkage to IBVD, it was found that the differences between the local and national parameters approach for the SP1 and SP2 subregions were ~1 m and ~0.35 m, which may be significant for some applications, especially in the SP1. This, together with the histogram analysis of the differences, corroborates the idea of using a parameter better adjusted to the local case.
In the parameters validation step, for the case of linkage to IBVD, in both subregions, regardless of the configuration of reference and check points used, more accurate results have been obtained with local modeling approach than with national modeling approach. Also in this context, a significant systematic bias (95% confidence level) in the differences of normal heights referred to IBVD has been verified using the national modeling approach. For the case of linkage to the future IHRS, only in subregion SP1 was there a significant systematic bias (95% confidence level) in differences from normal heights.
With this research, the perspective is broadened of determining normal heights of points of interest in a study region in Brazil, with a fieldwork of only a GNSS survey of the point. This survey integrated with a postprocessing, using the transformation parameter, allows obtaining the normal height of a given point referred to IBVD and the future IHRS with a decimetric precision. Although the calculated accuracies are below those desired by the IAG/GGOS, the local modeling approach is more accurate than a national modeling approach, with possibilities for adjustments in procedures and part of the improvement process that is sought at this level of research.
The methodology comprises a modeling with open data with wide access, provided by government agencies and institutions and widely accepted in the scientific community. Thus, the same procedure performed here can be extended to other states and regions of Brazil. Additionally, it is recommended analysis of:

other GGM/RTM solutions;

improvements in the calculation of the RTM effect, using tesseroids and harmonic correction;

studies of other linkage strategies, such as the one based on the solution of the GBVP (Rummel and Teunissen 1988Rummel, R. and Teunissen, P. 1988. Height datum definition, height datum connection and the role of the geodetic boundary value problem. Bulletin Géodésique , 62 (4), pp.477498. https://doi.org/10.1007/BF02520239.
https://doi.org/10.1007/BF02520239... ) and; 
more rigorous analysis of clusters to obtain smaller subregions, to improvements of accuracy.
Finally, it is worth noting that since the geodetic stations have vertical movement, because of temporal variations of physical heights, the parameters calculated in the research are valid for a given time and region (Ferreira et al. 2019Ferreira, V.G., Montecino, H.D., Ndehedehe, C.E., del Rio, R.A., Cuevas, A., and de Freitas, S.R.C. 2019. Determining seasonal displacements of Earth’s crust in South America using observations from spaceborne geodetic sensors and surfaceloading models. Earth, Planets and Space, 71 (1), pp.84. https://doi.org/10.1186/s4062301910622.
https://doi.org/10.1186/s406230191062...
)suggesting that South America’s hydrological cycle is changing, which impacts its water availability and, consequently, the Earth’s surface due to its elastic response to the surface mass loading/unloading. Therefore, we analyzed 3 to 15 years of vertical crustal displacements (VCDs. Due to the lack of an estimate of vertical velocity and acceleration, due to the difficulties of modeling it is causative factors, vertical coordinates may change over time at the centimetric level (Blewitt, Hammond, and Kreemer 2018Blewitt, G., Hammond, W., and Kreemer, C. 2018. Harnessing the GPS Data Explosion for Interdisciplinary Science. Eos, 99. https://doi.org/10.1029/2018EO104623.
https://doi.org/10.1029/2018EO104623...
).Thus, consideration is given to the need to update the parameters, within what is recommended by the methodology, whenever there is a change in epoch or region.
REFERENCES
 Amante, C and Eakins, B W 2009. ETOPO1 Global Relief Model converted to PanMap layer format. application/zip. https://doi.org/10.1594/PANGAEA.769615
» https://doi.org/10.1594/PANGAEA.769615  Blewitt, G., Hammond, W., and Kreemer, C. 2018. Harnessing the GPS Data Explosion for Interdisciplinary Science. Eos, 99. https://doi.org/10.1029/2018EO104623
» https://doi.org/10.1029/2018EO104623  Castro, H.D.M., Jamur, K.P., and Freitas, S.R.C. de 2012. Exploration of the residual topography effect on some functionals of the gravity field. Revista Brasileira de Geofísica, 30 (4), pp.431. https://doi.org/10.22564/rbgf.v30i4.230
» https://doi.org/10.22564/rbgf.v30i4.230  Drewes, H. 2006. The changing objectives in geodetic research. Zeitschrift für Geodäsie, GeoInformation und Landmanagement, 131 (5), pp.17.
 Ferreira, V.G., Montecino, H.D., Ndehedehe, C.E., del Rio, R.A., Cuevas, A., and de Freitas, S.R.C. 2019. Determining seasonal displacements of Earth’s crust in South America using observations from spaceborne geodetic sensors and surfaceloading models. Earth, Planets and Space, 71 (1), pp.84. https://doi.org/10.1186/s4062301910622
» https://doi.org/10.1186/s4062301910622  de Freitas, S.R.C. 2015. SIRGASWGIII activities for unifying height systems in Latin America. Revista Cartográfica, (91). DOI:10.35424/rcarto.i91.452
» https://doi.org/10.35424/rcarto.i91.452  Heck, B. and Seitz, K. 2007. A comparison of the tesseroid, prism and pointmass approaches for mass reductions in gravity field modelling. Journal of Geodesy, 81 (2), pp.121136. https://doi.org/10.1007/s0019000600940
» https://doi.org/10.1007/s0019000600940  Hirt, C., Featherstone, W.E., and Marti, U. 2010b. Combining EGM2008 and SRTM/DTM2006.0 residual terrain model data to improve quasigeoid computations in mountainous areas devoid of gravity data. Journal of Geodesy , 84 (9), pp.557567. https://doi.org/10.1007/s0019001003951
» https://doi.org/10.1007/s0019001003951  IBGE 2019. Reajustamento da Rede Altimétrica com Números Geopotenciais REALT2018 ENT#091;Readjustment of the Altimetric Network with Geopotential Numbers REALT2018ENT#093; ENT#091;onlineENT#093;. Available from: https://biblioteca.ibge.gov.br/visualizacao/livros/liv101666.pdf
» https://biblioteca.ibge.gov.br/visualizacao/livros/liv101666.pdf  Ihde, J., Sánchez, L., Barzaghi, R., Drewes, H., Foerste, C., Gruber, T., Liebsch, G., Marti, U., Pail, R., and Sideris, M. 2017. Definition and proposed realization of the International Height Reference System (IHRS). Surveys in Geophysics, 38 (3), pp.549570. https://doi.org/10.1007/s1071201794093
» https://doi.org/10.1007/s1071201794093  Ince, E.S., Barthelmes, F., Reißland, S., Elger, K., Förste, C., Flechtner, F., and Schuh, H. 2019. ICGEM  15 years of successful collection and distribution of global gravitational models, associated services, and future plans. Earth System Science Data, 11 (2), pp.647674. https://doi.org/10.5194/essd116472019
» https://doi.org/10.5194/essd116472019  MacMillan WD (1930): Theoretical Mechanics, vol 2: the Theory of the potential. McGrawHill, New York (reprinted by Dover Publications, New York 1958).
 Miot, H.A. 2017. Avaliação da normalidade dos dados em estudos clínicos e experimentais ENT#091;Assessment of data normality in clinical and experimental studiesENT#093;. Jornal Vascular Brasileiro, 16 (2), pp.8891. https://doi.org/10.1590/16775449.041117
» https://doi.org/10.1590/16775449.041117  Moritz, H. 1980. Geodetic reference system 1980. Bulletin Géodésique, 54 (3), pp.395405. https://doi.org/10.1007/BF02521480
» https://doi.org/10.1007/BF02521480  Nicacio, E., Dalazoana, R., and Freitas, S.R.C. de 2018. Evaluation of recent combined global geopotential models in Brazil. Journal of Geodetic Science, 8 (1), pp.7282. https://doi.org/10.1515/jogs20180008
» https://doi.org/10.1515/jogs20180008  Pail, R., Fecher, T., Barnes, D., Factor, J.F., Holmes, S.A., Gruber, T., and Zingerle, P. 2018. Short note: the experimental geopotential model XGM2016. Journal of Geodesy , 92 (4), pp.443451. https://doi.org/10.1007/s0019001710706
» https://doi.org/10.1007/s0019001710706  Rapp, R.H. 1989. The treatment of permanent tidal effects in the analysis of satellite altimeter data for sea surface topography. Manuscripta Geodaetica
 Rummel, R. and Teunissen, P. 1988. Height datum definition, height datum connection and the role of the geodetic boundary value problem. Bulletin Géodésique , 62 (4), pp.477498. https://doi.org/10.1007/BF02520239
» https://doi.org/10.1007/BF02520239  Sánchez, L. and Sideris, M.G. 2017. Vertical datum unification for the International Height Reference System (IHRS). Geophysical Journal International, pp.ggx025. https://doi.org/10.1093/gji/ggx025
» https://doi.org/10.1093/gji/ggx025  Tozer, B., Sandwell, D.T., Smith, W.H.F., Olson, C., Beale, J.R., and Wessel, P. 2019. Global Bathymetry and Topography at 15 Arc Sec: SRTM15+. Earth and Space Science, 6 (10), pp.18471864. https://doi.org/10.1029/2019EA000658
» https://doi.org/10.1029/2019EA000658  Zingerle, P., Pail, R., Gruber, T., and Oikonomidou, X. 2020. The combined global gravity field model XGM2019e. Journal of Geodesy , 94 (7), pp.66. https://doi.org/10.1007/s00190020013980
» https://doi.org/10.1007/s00190020013980
Publication Dates

Publication in this collection
24 June 2022 
Date of issue
2022
History

Received
13 Aug 2021 
Accepted
18 May 2022