Estimation of Vertical Datum Parameters Using the GBVP Approach Based on the Combined Global Geopotential Models

Unification of the global vertical datum has been a key problem to be solved for geodesy over a long period, and the main challenge for a unified vertical datum system is to determine the vertical offset between the local vertical datum and the global vertical datum. For this purpose, the geodetic boundary value problem (GBVP) approach based on the remove-compute-restore (RCR) technique is used to determine the vertical datum parameters in this paper. In the RCR technique, a global geopotential model (GGM) is required to remove and restore the long wavelengths of the gravity field. The satellite missions of the GRACE (Gravity Recovery and Climate Experiment) and GOCE (Gravity field and steady-state Ocean Circulation Exploration) offer high accuracy medium–long gravity filed information, but GRACE/GOCE-based GGMs are restricted to medium–long wavelengths because the maximum degree of their spherical harmonic representation is limited, which is known as an omission error. To compensate for the omission error of GRACE/GOCE-based GGM, a weighting method is used to determine the combined GGM by combining the high-resolution EGM2008 model (Earth Gravitational Model 2008) and GRACE/GOCE-based GGM to effectively bridge the spectral gap between satellite and terrestrial data. An additional consideration for the high-frequency gravity signals is induced by the topography, and the residual terrain model (RTM) is used to recover the omission errors effect of the combined GGM. In addition, to facilitate practical implementation of the GBVP approach, the effects of the indirect bias term, the spectral accuracy of the GGM, and the systematic levelling errors and distortions in estimations of the vertical datum parameters are investigated in this study. Finally, as a result of the GBVP solution based on the combined DIR_R6/EGM2008 model, RTM, and residual gravity, the geopotential values of the North American Vertical Datum of 1988 (NAVD88), the Australian Height Datum (AHD), and the Hong Kong Principal Datum (HKPD) are estimated to be equal to 62636861.31 ± 0.96, 62653852.60 ± 0.95 and 62636860.55 ± 0.29 m2s−2, respectively. The vertical offsets of NAVD88, AHD, and HKPD with respect to the global geoid are estimated as −0.809 ± 0.090, 0.082 ± 0.093, and −0.731 ± 0.030 m, respectively.


Introduction
The vertical datum system is a reference used for determining the physical heights of a country or region. A local vertical datum (LVD) is usually determined by the local mean sea level (MSL) observed by single or multiple tide gauge stations. However, the local mean sea level does not exactly coincide with the global geoid, which makes the vertical datums of many countries and regions contain offsets [1,2]. It is essential to realize a unified vertical datum, and there are many applications that require a unified vertical datum system, such as engineering construction, environmental monitoring and scientific research.
Unification of the existing vertical datum systems around the world is one of the major tasks of the Global Geodetic Observing System (GGOS) of the International Association of Geodesy (IAG). The IAG released a resolution for the definition and realization of an International Height Reference System (IHRS) during the 2015 International Union of Geodesy and Geophysics (IUGG) General Assembly, where the global equipotential reference surface was fixed by the geopotential value of W0 = 6,2636,853.4 m 2 s −2 [3,4]. According to this resolution, the existing local vertical datum systems can be integrated into the IHRS, which will ensure the consistency of the global vertical datum systems. At present, three strategies have been extensively discussed for vertical datum unification: the geodetic levelling approach, the oceanographic approach, and the gravity field approach. The geodetic levelling approach involves directly measuring the potential differences between two or more LVDs via spirit levelling in combination with gravity measurements. This approach provides sub-millimeter relative accuracy over short distances, but the absolute accuracy (relative to the global vertical datum) is about ±2 m [5]. If the two LVDs cannot be connected by spirit levelling, the approach is restricted. The oceanographic approach is used to determine the datum offsets between different LVDs through the mean dynamic topography (MDT), whose accuracy depends on the quality and spatial resolution of the ocean model [6,7]. However, the ocean model has low quality around tide gauge stations. The gravity field approach estimates the vertical datum offset by comparing the discrepancies between the biased geometric geoid height determined by GNSS/levelling (Global Navigation Satellite System/levelling) and the unbiased gravimetric geoid height determined by the gravity field data. This approach mainly includes two strategies: (1) GNSS/levelling combined with the global geopotential model (GGM) method [8][9][10] and (2) the geodetic boundary value problem (GBVP) method [11][12][13][14][15]. With the development of satellite-based GGMs, GNSS/levelling combined with the GGM method provides absolute accuracy from centimeters to decimeters [5,16], but the GGMs do not completely represent the signals of the Earth's gravity field because their accuracy and resolution are limited. The GBVP approach is a rigorous method for LVD unification, which uses the gravity anomaly data of the local vertical datum to determine the gravimetric geoid by solving the geodetic boundary value problem. Compared with other approaches, the GBVP approach has obvious advantages.
The key to determine the vertical datum offsets using the GBVP approach is to determine the gravimetric geoid based on the remove-compute-restore (RCR) technique [17,18]. In this technique, the GGM and digital terrain model (DTM) provide the long and short wavelength parts of the gravity field spectrum, and the residual gravity field spectrum is provided by the residual gravity. Although the RCR technique has many advantages, its results depend largely on the accuracy of the GGM. With the implementation of gravity satellite missions like GRACE (Gravity Recovery and Climate Experiment) and GOCE (Gravity field and steady-state Ocean Circulation Exploration), these satellite missions offer unprecedented high-accuracy information for the medium-frequency of the gravity field spectrum; moreover, the accuracy of the medium-long wavelength in the static gravity field has been improved by two to three orders of magnitude [19][20][21]. However, the maximum expansion degree of GRACE/GOCE-based GGMs is limited, and these models have certain omission errors due to the influence of residual gravity field signals. Sánchez [22] demonstrated that the global mean omission error is about ±45 cm for GRACE/GOCE-based GGMs with a degree of 200. In addition to the omission errors of the GGMs, the factors that affect the estimation of vertical datum offsets Remote Sens. 2020, 12, 4137 3 of 23 also include the spectral accuracy of the GGM, the indirect bias term, and the systematic errors and distortions in the vertical control networks. Therefore, to determine high-accuracy vertical datum offsets, the aforementioned factors must be considered.
The objective of this study is to determine the vertical datum parameters of the USA, Australia, and Hong Kong based on the GBVP approach. For practical implementation of the GBVP approach to unify the vertical datum, the omission errors of the GOCE/GRACE-based GGMs, the spectral accuracy of the GGM, the indirect bias term, and the systematic errors and distortions in the vertical networks are investigated in this paper. To reduce the geoid omission error effect of the GOCE/GRACE-based GGMs, we combine GOCE/GRACE-based GGMs and the EGM2008 model to determine the combined GGMs via the weighting method. The combined GGMs include the high-accuracy medium and long wavelength information of the GRACE/GOCE-based GGMs and retain the short wavelength information of EGM2008. However, the maximum degree and order of the combined GGMs is only 2190, which is equivalent to a spatial resolution of 5 × 5 . The combined GGMs, however, fail to represent a high-frequency geoid signal with a spatial resolution greater than 5 × 5 . Considering the influence of the high-frequency gravity field signal above a degree 2190, the residual terrain model technique [23][24][25][26][27] is used to recover the high-frequency gravity field signal. Finally, the GBVP approach based on the combined GGM, residual terrain model (RTM), and residual gravity is used to estimate the vertical datum parameters of the North America Vertical Datum of 1988 (NAVD88), the Australia Height Datum (AHD), and the Hong Kong Principal Datum (HKPD).

Determination of the Vertical Datum Parameters Based on the Geodetic Boundary Value Problem (GBVP) Approach
According to the generalized Burns' formula [28], for a point P, the geoid height N P can be expressed as: where ∆W 0 is the difference between the geopotential W 0 of the global geoid and the normal geopotential U 0 of the reference ellipsoid, δW LVD 0 is the potential difference between the geopotential W 0 of the global geoid and the geopotential W LVD 0 of the local vertical datum (see Figure 1), γ is the normal gravity on the reference ellipsoid, and T P is the disturbing potential. In this study, the geopotential W 0 = 62636853.4 m 2 s −2 of the global geoid released by the IHRS is used. The ellipsoid parameters used in this paper are based on the GRS80 ellipsoid [29].
Remote Sens. 2020, 12, x 5 of 23 theoretically, be a fixed constant. However, vertical offsets have certain differences due to the influence of systematic errors such as systematic errors and distortions in the leveling network and random errors in computation of the ellipsoidal, as well as geoidal heights. To reduce the random and systematic errors, the vertical datum offset will be estimated by applying a parametric model [32,33]. For each GNSS/leveling point P, the observation equation can be formulated as: where H δ is mean vertical offset, ( ) can be expressed as: where γ is the mean value of the normal gravity of all the GNSS/levelling benchmarks.
Consequently, we can further determine the geopotential value 0 LVD W of the local vertical datum as: Figure 1. Vertical datum parameters (left of red dotted line) and relations between the ellipsoidal heights, levelling height and geometric geoid height at benchmark P (right of the red dotted line). Note that P and P' are the same.

Determination of the Combined Global Geopotential Models (GGMs) by the Weighting Method
The GGM N in Equation (10) can be determined by the GRACE/GOCE-based GGMs. However, the omission errors of the GRACE/GOCE-based GGMs must be considered. EGM2008 makes full use The unification of the vertical datum requires the gravity anomalies ∆g based on the global datum, but we can only obtain the local vertical datum gravity anomalies ∆g l in the practice case. The relationship between the ∆g and the ∆g l can be expressed as [30]: where R is the mean Earth radius. The solution of the GBVP in terms of its disturbing potential is [14,31] where δGM is the difference between the gravitational constant of the Earth and the gravitational constant of the reference ellipsoid, and S(ψ) is the stokes kernel function. By applying Equation (3) to (1) and using N 0 = (δGM/Rγ − ∆W 0 /γ), N Stokes = (R/4πγ) Ω S(ψ)∆g l dσ, and δH = δW LVD 0 /γ, Equation (1) can be written as follows: where N 0 is the zero-degree term of the geoid height, and δH is the vertical offset between the LVD and the global vertical datum. Using the RCR technique, N P Stokes in Equation (4) can be computed from the GGM and stokes integration with the residual gravity anomalies. To obtain high-frequency geoid information, the RTM is used to recover the short-wavelength geoid in this paper. Then, Equation (4) can be written as follows: 2π Ω δH·S res (ψ)dσ (5) where N P GGM and N P RTM are the GGM geoid height and the RTM geoid height, respectively, ∆g res = ∆g l − ∆g GGM − ∆g RTM is the residual gravity anomalies, and the fifth term on the right-hand side of Equation (5) is the residual geoid N P res . The modified Stokes kernel function S res (ψ) can then be expressed as: where N is the maximum degree of the GGM, and P n (cos ψ) is the Legendre polynomials. The sixth term in the right-hand side of Equation (5) is the indirect bias term. If there are k vertical datum zones, and the vertical offset δH is constant in each datum zone, the indirect bias term N ind can be expressed as: where ψ 1 is the spherical distance between the computation point and the vertical datum zone Ω i , ψ 2 = ψ 1 + ∆ψ i , ∆ψ i is the width of the datum zone, and δH i represent the average of all datum offsets integrated over azimuth for vertical datum zone i. The integration element in Equation (7) is in terms of polar coordinates and carried out the integration over azimuth. Then, they are left with integration over spherical distance. Gerlach and Rummel [30] and Amjadiparvar et al. [14] showed that the indirect bias term can be neglected when a GGM with a higher degree and order is used. The indirect bias term is also evaluated in Section 4.3 in this paper. The result shows that the effect of the indirect bias term is equal to 1 mm when higher degree GGMs (300 degree and order) are used, as the maximum degree of the GGM increases, the indirect bias effect becomes smaller. Therefore, the effect of the indirect bias term can be ignored. The geoid height N P is also computed from GNSS/levelling as: where h P is the ellipsoidal height, and H P is the levelling height (see Figure 1). By combining Equations (5) and (8), the vertical offset δH P for point P can be expressed as: The datum surface corresponding to the global or local vertical datum is a gravitational equipotential surface. Therefore, the vertical datum offset δH P determined by Equation (9) should, theoretically, be a fixed constant. However, vertical offsets have certain differences due to the influence of systematic errors such as systematic errors and distortions in the leveling network and random errors in computation of the ellipsoidal, as well as geoidal heights. To reduce the random and systematic errors, the vertical datum offset will be estimated by applying a parametric model [32,33]. For each GNSS/leveling point P, the observation equation can be formulated as: where δH is mean vertical offset, (B P , L P ) is the geodetic coordinates, (B 0 , L 0 ) is the mean value of the geodetic coordinates of all GNSS/levelling benchmarks in the local vertical datum zone, and a 1 and a 2 are the north-south tilt and east-west tilt, respectively. The vertical datum parameters can be expressed in the form of the vertical offset δH, the potential difference δW LVD 0 , and geopotential W LVD 0 of the local vertical datum (see Figure 1). The δH can be estimated by means of least square adjustment (LSA) of the system in Equation (10). According to the relationship between the vertical offset and potential difference, the potential difference δW LVD 0 can be expressed as: where γ is the mean value of the normal gravity of all the GNSS/levelling benchmarks. Consequently, we can further determine the geopotential value W LVD 0 of the local vertical datum as:

Determination of the Combined Global Geopotential Models (GGMs) by the Weighting Method
The N GGM in Equation (10) can be determined by the GRACE/GOCE-based GGMs. However, the omission errors of the GRACE/GOCE-based GGMs must be considered. EGM2008 makes full use of terrestrial, airborne and satellite altimetry gravity data [34], so this model can effectively represent the short wavelength gravity field signal. Consequently, the EGM2008 model can be used to extend GOCE/GRACE-based GGMs to compensate for omission errors. We combine GOCE/GRACE-based GGMs and the EGM2008 model to determine the combined GGMs via the weighting method [35][36][37]. The degree 0-N of the combined GGM (N is the maximum degree of the GOCE/GRACE-based GGMs) is determined by spectrum combination, and the degree from N+1 to 2190 is supplemented by the corresponding degree and order of the EGM2008 model. The combined GGMs can have the high-accuracy medium and long wavelength information of the GRACE/GOCE GGMs and retain the Remote Sens. 2020, 12, 4137 6 of 23 short wavelength information of EGM2008. The spherical harmonic coefficients combined using by the weighting method can be expressed as: where C EGM08 nm and S EGM08 nm are the EGM2008 spherical harmonic coefficients of degree n and order m, C GOCE/GRACE nm and S GOCE/GRACE nm are the GOCE/GRACE-based GGMs' spherical harmonic coefficients of degree n and order m, and p C nm and p S nm represents the spectral weight of the EGM2008 model, and they can be calculated as: where σ EGM08 C and σ EGM08 S are the coefficient standard deviations of the EGM2008 model, and σ GOCE/GRACE C and σ GOCE/GRACE S are the coefficient standard deviations of GOCE/GRACE-based GGMs. According to the error propagation law, the coefficient standard deviations of the spectrum combination can be calculated as:

Data Sets
This section will briefly introduce the data sources used. The main data required are: (1) global geopotential models; (2) GNSS/levelling data; (3) gravity data; and (4) high resolution topographic data.

Global Geopotential Models (GGMs)
There are seven GGMs used in this study: DIR_R6 [38], TIM_R6 [39], GOSG01S [40], IfE_GOCE05s [41], IGGT_R1 [42], SPW_5 [43], and EGM2008 [34], as shown in Table 1. Except for TIM_R6, which is based on the zero-tide system, the other models are based on a tide-free system. For the spherical harmonic coefficients, only C 20 is affected by the tide system. Hence, the C 20 coefficient of TIM_R6 is transformed from zero-tide to tide-free system according to Sánchez et al. [44]: where C TF 20 is the spherical harmonic coefficient under the tide-free system, and C ZF 20 is spherical harmonic coefficient under the zero-tide system.

GNSS/Levelling Data
For mainland America, a GNSS/levelling data set consisting of 23961 benchmarks was made available by the National Geodetic Survey (NGS). These data include the ellipsoidal heights in NAD83 and levelling heights with respect to the North American Vertical Datum of 1988 (NAVD88). The Australian State and Territory geodetic agencies provide a set of 7545 GNSS/levelling points, consisting of GNSS-based ellipsoidal heights in ITRF14 and levelling heights with respect to the Australian Height Datum (AHD). We used 122 GNSS/levelling benchmarks evenly distributed in Hong Kong, the ellipsoidal heights are given in ITRF96, and the levelling heights determined by precise geometric leveling or trigonometric leveling are given with respect to the HKPD. The spatial distribution of the GNSS/levelling benchmarks for Mainland America, Mainland Australia, and Hong Kong are shown in Figure 2.

GNSS/Levelling Data
For mainland America, a GNSS/levelling data set consisting of 23961 benchmarks was made available by the National Geodetic Survey (NGS). These data include the ellipsoidal heights in NAD83 and levelling heights with respect to the North American Vertical Datum of 1988 (NAVD88). The Australian State and Territory geodetic agencies provide a set of 7545 GNSS/levelling points, consisting of GNSS-based ellipsoidal heights in ITRF14 and levelling heights with respect to the Australian Height Datum (AHD). We used 122 GNSS/levelling benchmarks evenly distributed in Hong Kong, the ellipsoidal heights are given in ITRF96, and the levelling heights determined by precise geometric leveling or trigonometric leveling are given with respect to the HKPD. The spatial distribution of the GNSS/levelling benchmarks for Mainland America, Mainland Australia, and Hong Kong are shown in Figure 2.
To ensure consistent calculations, the geometric coordinates in mainland America and Hong Kong are transformed into ITRF2014. In this study, the tide-free system will be consistently used for all quantities. The levelling heights are transformed from mean tide to the tide-free system based on the formula below [45]: where H TF is the tide-free levelling height, H MF is the mean tide levelling height, and ϕ is the latitude of the GNSS/levelling benchmarks. To ensure consistent calculations, the geometric coordinates in mainland America and Hong Kong are transformed into ITRF2014. In this study, the tide-free system will be consistently used for all quantities. The levelling heights are transformed from mean tide to the tide-free system based on the formula below [45]: where H TF is the tide-free levelling height, H MF is the mean tide levelling height, and ϕ is the latitude of the GNSS/levelling benchmarks.

Gravity Data
In total, 1,633,499 gravity data in the USA were made available by the National Oceanic and Atmospheric Administration (NOAA). These data mainly include land gravity observations in Mainland America, Alaska, and Hawaii along with airborne gravity data near the coastline, as well as a small number of gravity observation data of Canada and Mexico. After removing gravity observation data from Alaska, Hawaii, and airborne gravity data, we selected 822,301 gravity observation data distributed across the land. For Australia, we used 1,835,358 gravity observation data from the Geoscience Australia's national gravity database; we also used 610 gravity points in Hong Kong. The spatial distribution of the gravity station points for the USA, Australia, and Hong Kong are shown in Figure 3a,c,e.
Mainland America, Alaska, and Hawaii along with airborne gravity data near the coastline, as well as a small number of gravity observation data of Canada and Mexico. After removing gravity observation data from Alaska, Hawaii, and airborne gravity data, we selected 822,301 gravity observation data distributed across the land. For Australia, we used 1,835,358 gravity observation data from the Geoscience Australia's national gravity database; we also used 610 gravity points in Hong Kong. The spatial distribution of the gravity station points for the USA, Australia, and Hong Kong are shown in Figure 3a,c,e.
We calculated normal gravity via the closed form of the 1980 international gravity formula [29] and applied second-order free air and atmospheric corrections [28] to the point gravity values to obtain free-air gravity anomalies. The bicubic interpolation method was employed to interpolate the point gravity anomalies data into 1 1 ′ ′ × gravity anomaly grid for the USA and Australia and 3 3

′′ ′′ ×
gravity anomaly grid for Hong Kong. Because the point gravity anomalies are mainly distributed across the land, for offshore gravity data, we chose the latest grav.img.28.1 (ftp://topex.ucsd.edu/pub/global_grav_1min/) model derived from multi-mission satellite altimetry. For USA and Australia, the land gravity grid data were combined with offshore gravity data to obtain unified resolution (1 arc-minute) gravity anomaly grid data for land and sea. For Hong Kong, the bicubic interpolation method was employed to interpolate the 1 1 ′ ′ × offshore gravity data into a 3 3 ′′ ′′ × grid, then, the land gravity grid data were combined with the 3 3

′′ ′′ ×
offshore gravity data to obtain unified resolution (3 arc-second) gravity anomalies grid data on land and sea. Figure 3b,d,f shows the unified gravity anomalies grid in the USA, Australia, and Hong Kong, respectively.

Topographic Data
The topographic data were used to calculate the RTM effect. SRTM (Shuttle Radar Topographic Mission) data [46] with a spatial resolution of 7.5 7.5 ′′ ′′ × were used for USA and Australia over land, and SRTM with a spatial resolution of 3 3

′′ ′′ ×
was used for Hong Kong over land. In the sea area, the latest SRTM15_PLUS V2 topographic data [47] were used, with a spatial resolution of 15″×15″. We calculated normal gravity via the closed form of the 1980 international gravity formula [29] and applied second-order free air and atmospheric corrections [28] to the point gravity values to obtain free-air gravity anomalies. The bicubic interpolation method was employed to interpolate the point gravity anomalies data into 1 × 1 gravity anomaly grid for the USA and Australia and 3 × 3 gravity anomaly grid for Hong Kong. Because the point gravity anomalies are mainly distributed across the land, for offshore gravity data, we chose the latest grav.img.28.1 (ftp://topex.ucsd.edu/pub/global_grav_1min/) model derived from multi-mission satellite altimetry. For USA and Australia, the land gravity grid data were combined with offshore gravity data to obtain unified resolution (1 arc-minute) gravity anomaly grid data for land and sea. For Hong Kong, the bicubic interpolation method was employed to interpolate the 1 × 1 offshore gravity data into a 3 × 3 grid, then, the land gravity grid data were combined with the 3 × 3 offshore gravity data to obtain unified resolution (3 arc-second) gravity anomalies grid data on land and sea. Figure 3b,d,f shows the unified gravity anomalies grid in the USA, Australia, and Hong Kong, respectively.

Topographic Data
The topographic data were used to calculate the RTM effect. SRTM (Shuttle Radar Topographic Mission) data [46] with a spatial resolution of 7.5 × 7.5 were used for USA and Australia over land, and SRTM with a spatial resolution of 3 × 3 was used for Hong Kong over land. In the sea area, the latest SRTM15_PLUS V2 topographic data [47] were used, with a spatial resolution of 15" × 15". RET2012 was used as a topographic reference model. The RET2012 reference terrain elevations are synthesized using [48]: (HC nm cos mλ + HS nm sin mλ)P nm (cos θ) (18) where N = 2160 is the maximum degree, λ and θ are the longitude and the geocentric co-latitude of the computation points, respectively, and P nm (cos θ) is the associated fully-normalized spherical Legendre functions of degree n and order m.
In coastal areas, the influence of water depth on RTM cannot be ignored. Because the mass-density of sea water is inconsistent with the standard topographic mass density, to avoid the need to distinguish density changes in the calculation process, the rock-equivalent topography (RET) method was used to compress the water depth into the equivalent rock height [49,50], then, the SRTM and SRTM15_PLUS data were merged to determine the unified land and sea topographic data. This merging was accomplished in a two-step procedure. Firstly, the RET method was employed to process the sea water depths, and the bicubic interpolation method was employed to interpolate the data into a spatial resolution the same as that of the land topographic data; then, the land data were combined with the water depth data to obtain unified resolution topographic data on land and sea. Figure 4a,c,e shows the topography merged for the USA, Australia, and Hong Kong, respectively. The RTM is the difference between the topography merged and reference terrain. Figure 4b,d,f shows the RTM elevation for the USA, Australia, and Hong Kong, respectively.

Spectral Accuracy Evaluation of the GGMs
In this study, the spectral accuracy of each GGM is evaluated according to its degree error and signal-to-noise ratio (SNR). The degree error is computed as the square root of the error degree variances, and the SNR represents the ratio of the geoid signal to the degree error. The calculation formulas for degree error and SNR can be found in Ustun et al. [51]. Figure 5a shows the degree errors of the applied GGMs in terms of geoid height. It can be seen from the figure that the degree error of the EGM2008 model gradually becomes lower than that of the GRACE/GOCE-based GGMs after about degree 220, indicating that the spectral accuracy of the EGM2008 model is better than that of the GRACE/GOCE-based GGMs with a higher degree and order. The degree errors of the DIR_R6, TIM_R6, GOSG01S, IfE_GOCE05s, and IGGT_R1 models before about degrees 220, 203, 202, 209, 175, respectively, are lower than the degree errors of EGM2008 model. From about degrees 65 to 197, the degree error of the SPW_5 model is smaller than that of the

Spectral Accuracy Evaluation of the GGMs
In this study, the spectral accuracy of each GGM is evaluated according to its degree error and signal-to-noise ratio (SNR). The degree error is computed as the square root of the error degree variances, and the SNR represents the ratio of the geoid signal to the degree error. The calculation formulas for degree error and SNR can be found in Ustun et al. [51]. Figure 5a shows the degree errors of the applied GGMs in terms of geoid height. It can be seen from the figure that the degree error of the EGM2008 model gradually becomes lower than that of the GRACE/GOCE-based GGMs after about degree 220, indicating that the spectral accuracy of the EGM2008 model is better than that of the GRACE/GOCE-based GGMs with a higher degree and order. The degree errors of the DIR_R6, TIM_R6, GOSG01S, IfE_GOCE05s, and IGGT_R1 models before about degrees 220, 203, 202, 209, 175, respectively, are lower than the degree errors of EGM2008 model. From about degrees 65 to 197, the degree error of the SPW_5 model is smaller than that of the EGM2008. The results indicate that the DIR_R6, TIM_R6, GOSG01S, IfE_GOCE05s, IGGT_R1, and SPW_5 models have better spectral accuracy in medium and long wavelengths.
The geoid signal power of GGMs at different wavelengths could be expressed by SNR. Figure  5b shows the SNR of the GGMs. The SNRs of the DIR_R6, TIM_R6, GOSG01S, IfE_GOCE05s, and IGGT_R1 models are better than those of the EGM2008 model before about degrees 220, 204, 203, 209 and 175, respectively. The SNR of the SPW_5 model is better than that of the EGM2008 model from degrees 65 to 198. The results show that the GRACE/GOCE-based GGMs have a stronger geoid signal in the medium and long wavelength, and the EGM2008 model has a stronger geoid signal in the short wavelength. From the perspective of the SNR and degree error, the spectral accuracy of the DIR_R6 model is significantly better than that of the other GRACE/GOCE GGMs, because the DIR_R6 model takes advantage of GRACE gravity data and LAGEOS Satellite Laser Ranging (SLR) data.

Omission Errors of the GGMs
The GRACE/GOCE-based GGMs have high medium and long wavelength accuracy, but the maximum degree of these models is limited, as there are certain omission errors affected by the roughness of the remaining (residual) gravity field signals. To quantify the magnitude of these omission errors, we assume that degree 300 is the general expansion degree of the GRACE/GOCEbased GGMs and that the geoid represented from degree 301 to 2190 in the EGM2008 model can be viewed as the omission errors of the GRACE/GOCE-based GGMs. Figure 6 shows the spatial distribution of the omission errors for the USA, Australia, and Hong Kong, the corresponding statistics are presented Table 2. It can be seen from the Table that the effect of the omission errors reaches the dm level in the USA and Australia, while the impact of omission errors in Hong Kong is at the cm level. The omission errors of the largest amplitude with -159 cm and 88 cm are reached in USA and Australia, respectively, and the mean omission errors in the USA and Australia reach -2.2 and -2.4 cm, respectively. In Hong Kong, the mean omission error is about 1 cm. In Figure 5, we can also see that the omission errors in the rugged regions (for example, the western region of the USA) are larger. Therefore, the omission errors of the GRACE/GOCE-based GGMs must be considered for the purpose of the vertical datum offset estimation. The geoid signal power of GGMs at different wavelengths could be expressed by SNR. Figure 5b shows the SNR of the GGMs. The SNRs of the DIR_R6, TIM_R6, GOSG01S, IfE_GOCE05s, and IGGT_R1 models are better than those of the EGM2008 model before about degrees 220, 204, 203, 209 and 175, respectively. The SNR of the SPW_5 model is better than that of the EGM2008 model from degrees 65 to 198. The results show that the GRACE/GOCE-based GGMs have a stronger geoid signal in the medium and long wavelength, and the EGM2008 model has a stronger geoid signal in the short wavelength. From the perspective of the SNR and degree error, the spectral accuracy of the DIR_R6 model is significantly better than that of the other GRACE/GOCE GGMs, because the DIR_R6 model takes advantage of GRACE gravity data and LAGEOS Satellite Laser Ranging (SLR) data.

Omission Errors of the GGMs
The GRACE/GOCE-based GGMs have high medium and long wavelength accuracy, but the maximum degree of these models is limited, as there are certain omission errors affected by the roughness of the remaining (residual) gravity field signals. To quantify the magnitude of these omission errors, we assume that degree 300 is the general expansion degree of the GRACE/GOCE-based GGMs and that the geoid represented from degree 301 to 2190 in the EGM2008 model can be viewed as the omission errors of the GRACE/GOCE-based GGMs. Figure 6 shows the spatial distribution of the omission errors for the USA, Australia, and Hong Kong, the corresponding statistics are presented Table 2. It can be seen from the Table that the effect of the omission errors reaches the dm level in the USA and Australia, while the impact of omission errors in Hong Kong is at the cm level. The omission errors of the largest amplitude with −159 cm and 88 cm are reached in USA and Australia, respectively, and the mean omission errors in the USA and Australia reach −2.2 and −2.4 cm, respectively. In Hong Kong, the mean omission error is about 1 cm. In Figure 5, we can also see that the omission errors in the rugged regions (for example, the western region of the USA) are larger. Therefore, the omission errors of the GRACE/GOCE-based GGMs must be considered for the purpose of the vertical datum offset estimation.

Effect of the Indirect Bias Term
The effect of the indirect bias term is assessed by means of Equation (7). The indirect bias term is not only closely related to the vertical offset H δ and integral function, but it corresponds to δ H being integrated over the respective datum zone (weighted average using stokes function as distance depending weights). The local mean sea level does not coincide with the global geoid, and the resulting discrepancies between the different local vertical datum zones cause global height datum offsets of about 1~2 m. To evaluate the indirect bias term, we assume 1 H δ = m. The datum zones are divided into near zone and far zone. The datum zone where the calculation point is located is called the near zone-that is, 1 =0 ψ ; the datum zone outside calculation point is called the far zonethat is, 1 0 ψ ≠ . Then, the indirect bias term of the calculation point is divided into two aspects: one is the indirect bias effect of the near zone, and the other is the indirect bias effect of the far zone. Figure 7a shows the indirect bias effect as a function of the size ψ

Effect of the Indirect Bias Term
The effect of the indirect bias term is assessed by means of Equation (7). The indirect bias term is not only closely related to the vertical offset δH and integral function, but it corresponds to δH being integrated over the respective datum zone (weighted average using stokes function as distance depending weights). The local mean sea level does not coincide with the global geoid, and the resulting discrepancies between the different local vertical datum zones cause global height datum offsets of about 1~2 m. To evaluate the indirect bias term, we assume δH = 1 m. The datum zones are divided into near zone and far zone. The datum zone where the calculation point is located is called the near zone-that is, ψ 1 = 0; the datum zone outside calculation point is called the far zone-that is, ψ 1 0. Then, the indirect bias term of the calculation point is divided into two aspects: one is the indirect bias effect of the near zone, and the other is the indirect bias effect of the far zone. Figure 7a shows the indirect bias effect as a function of the size ∆ψ of the datum zone for different maximum degrees N. This figure reflects the indirect bias term effect of the near zone to the calculation point. The indirect bias effect reaches a maximum of about 2 cm for N = 50. As the degree of the GGM increases, the indirect bias effect gradually becomes smaller, especially for N = 300, where the maximum indirect bias effect is only 0.1 cm. point gradually becomes smaller with an increase of 1 ψ . Therefore, the corresponding 1 ψ when the indirect bias term effect is the largest can be selected to evaluate the total indirect bias term effect of the far zone and the near zone to the calculation point. The total indirect bias term effect for different maximum degrees N is shown in Figure 7c, where the maximum total indirect bias effect is only 0.1 cm for N = 300. Therefore, the indirect bias term effect can be ignored by using GGMs of a high degree and order.

Determination the Combined GGM
The omission errors of the GRACE/GOCE-based GGMs must be considered for the purpose of vertical datum offset estimation. The DIR_R6 model remains better spectral accuracy than the TIM_R6, GOSG01S, IfE_GOCE05s, IGGT_R1, and SPW_R5 models (see Section 4.1). Therefore, we built a combined model DIR_R6/EGM2008 by combining the DIR_R6 model and the EGM2008 model using the weighting method (Section 2.2). Figure 8 shows the geoid height signals and geoid errors of the DIR_R6, EGM2008, and DIR_R6/EGM2008 models. The geoid height signal of DIR_R6/EGM2008 is so close to that of the EGM2008 model that the two curves are identical after degree 300. The geoid height signal provided by the DIR_R6 model shows no difference with DIR_R6/EGM2008 before about degree 220, and after degree 220, it is lower than that of the DIR_R6/EGM2008 model. The geoid degree error of the  It can be seen from Figure 7b that the indirect bias term effect of the far zone to the calculation point gradually becomes smaller with an increase of ψ 1 . Therefore, the corresponding ψ 1 when the indirect bias term effect is the largest can be selected to evaluate the total indirect bias term effect of the far zone and the near zone to the calculation point. The total indirect bias term effect for different maximum degrees N is shown in Figure 7c, where the maximum total indirect bias effect is only 0.1 cm for N = 300. Therefore, the indirect bias term effect can be ignored by using GGMs of a high degree and order.

Determination the Combined GGM
The omission errors of the GRACE/GOCE-based GGMs must be considered for the purpose of vertical datum offset estimation. The DIR_R6 model remains better spectral accuracy than the TIM_R6, GOSG01S, IfE_GOCE05s, IGGT_R1, and SPW_R5 models (see Section 4.1). Therefore, we built a combined model DIR_R6/EGM2008 by combining the DIR_R6 model and the EGM2008 model using the weighting method (Section 2.2). Figure 8 shows the geoid height signals and geoid errors of the DIR_R6, EGM2008, and DIR_R6/EGM2008 models. The geoid height signal of DIR_R6/EGM2008 is so close to that of the EGM2008 model that the two curves are identical after degree 300. The geoid height signal provided by the DIR_R6 model shows no difference with DIR_R6/EGM2008 before about degree 220, and after degree 220, it is lower than that of the DIR_R6/EGM2008 model. The geoid degree error of the DIR_R6/EGM2008 is approximately the same as that of the DIR_R6 model before degree 200, after which, the geoid degree error of DIR_R6/EGM2008 is significantly lower than that of the DIR_R6 model. The degree error of the DIR_R6/EGM2008 model is lower than that of the EGM2008 model before about degree 250. The cumulative degree error of the DIR_R6/EGM2008 model is consistent with the DIR_R6 model before degree 150, and from degree 150 onwards, it becomes significantly lower than the cumulative degree error of the DIR_R6 model. Compared with the EGM2008 model, the cumulative degree error of DIR_R6/EGM2008 is significantly lower than that of the EGM2008 model over the entire spectral domain.
Remote Sens. 2020, 12, x 14 of 23 DIR_R6/EGM2008 is approximately the same as that of the DIR_R6 model before degree 200, after which, the geoid degree error of DIR_R6/EGM2008 is significantly lower than that of the DIR_R6 model. The degree error of the DIR_R6/EGM2008 model is lower than that of the EGM2008 model before about degree 250. The cumulative degree error of the DIR_R6/EGM2008 model is consistent with the DIR_R6 model before degree 150, and from degree 150 onwards, it becomes significantly lower than the cumulative degree error of the DIR_R6 model. Compared with the EGM2008 model, the cumulative degree error of DIR_R6/EGM2008 is significantly lower than that of the EGM2008 model over the entire spectral domain. Therefore, the best combination of the EGM2008 model and DIR_R6 model can be achieved by using the weighting method. The DIR_R6/EGM2008 model offers stronger geoid signal and higher accuracy compared to the EGM2008 and DIR_R6 models. The DIR_R6/EGM2008 model provides the high-accuracy medium and long wavelength information of the DIR_R6 model and retains the short wavelength information of EGM2008.

Residual Gravity Anomalies
The residual gravity anomalies =  [52] proposed the lumped coefficient approach (LCA), but this approach cannot be applied to irregular grid surfaces. The LCA approach regards the calculated grid surface as a regular surface (a sphere or ellipsoid) and does not consider the influence of the actual terrain surface. Hirt [53] presented a gradient approach. According to this approach, the gravity field functions and its derivatives of chosen order can be computed on the sphere or the ellipsoid employing the LCA, and subsequently continued to an irregular surface. Therefore, we use the gradient approach to determine the GGM g Δ in this study.
The RTM gravity anomalies are estimated via prism integration [54] based on planar approximation. Because the spatial resolution of the reference surface is equivalent to that of the DIR_R6/EGM2008 model, the RTM gravity anomalies can represent the high-frequency gravity field signal with a spatial resolution of less than 5 arc-minutes. Figure 9 shows the residual gravity anomalies res g Δ and the corresponding frequency statistics histograms. The resolution of residual gravity anomalies in the USA and Australia is 1 1 ′ ′ × , and the resolution of residual gravity anomalies in Hong Kong is 3 3 ′′ ′′ × . It can be seen from the histograms that the residual gravity anomalies in the USA, Australia, and Hong Kong maintain a normal Therefore, the best combination of the EGM2008 model and DIR_R6 model can be achieved by using the weighting method. The DIR_R6/EGM2008 model offers stronger geoid signal and higher accuracy compared to the EGM2008 and DIR_R6 models. The DIR_R6/EGM2008 model provides the high-accuracy medium and long wavelength information of the DIR_R6 model and retains the short wavelength information of EGM2008.

Residual Gravity Anomalies
The residual gravity anomalies ∆g res = ∆g l − ∆g GGM − ∆g RTM are used to calculate the residual geoid N res in Equation (9). The ∆g GGM is determined by the DIR_R6/EGM2008 model. However, calculating the high-resolution ∆g GGM point by point will reduce computational efficiency. To resolve the problem of numerical computational efficiency, Colombo et al. [52] proposed the lumped coefficient approach (LCA), but this approach cannot be applied to irregular grid surfaces. The LCA approach regards the calculated grid surface as a regular surface (a sphere or ellipsoid) and does not consider the influence of the actual terrain surface. Hirt [53] presented a gradient approach. According to this approach, the gravity field functions and its derivatives of chosen order can be computed on the sphere or the ellipsoid employing the LCA, and subsequently continued to an irregular surface. Therefore, we use the gradient approach to determine the ∆g GGM in this study. The RTM gravity anomalies are estimated via prism integration [54] based on planar approximation. Because the spatial resolution of the reference surface is equivalent to that of the DIR_R6/EGM2008 model, the RTM gravity anomalies can represent the high-frequency gravity field signal with a spatial resolution of less than 5 arc-minutes. Figure 9 shows the residual gravity anomalies ∆g res and the corresponding frequency statistics histograms. The resolution of residual gravity anomalies in the USA and Australia is 1 × 1 , and the  Table 3. Statistics of the gravity anomalies and each reduction term as well as the final residuals for the three study areas of the USA, Australia, and Hong Kong. Unit: mGal (1mGal = 10 −5 m/s 2 ).

Estimation of Vertical Datum Parameters
The geopotential W 0 = 62,636,853.4 m 2 s −2 of the global vertical datum released by the IHRS is used here. Based on the GBVP approach, the vertical datum parameters of the USA, Australia, and Hong Kong are determined. The key to determine the vertical datum parameters via the GPVP approach is to use the RCR technique to determine the gravimetric geoid. The N GGM is determined by the DIR_R6/EGM2008 model, the RTM geoid height N RTM is estimated by prism integration [25], and the residual geoid height N res is calculated by stokes integration using 1D-FFT (One Dimensional Fast Fourier Transform) approach. The stokes integration radius is as same as the computation area. The RTM integration radius depends on the gravity functional, as well as on the spectral energy of the RTM elevations, for RTM geoid height (the spatial scales shorter than 5 arc-min), an integration radius of~200 km is suitable [50]. Figure 10a shows the geoid contribution of residual gravity and RTM in USA, with a maximum, minimum, mean, and standard deviation of 0.466 m, −0.430 m, 0.001 m, and 0.02 m, respectively. Figure 10c shows the geoid contribution of residual gravity and RTM in Australia, with a maximum, minimum, mean, and standard deviation of 0.539 m, −0.372 m, 0.012 m, and 0.030 m, respectively. Figure 10e shows the geoid contribution of residual gravity and RTM in Hong Kong, with a maximum, minimum, mean, and standard deviation of 0.073 m, −0.039 m, −0.005 m, and 0.010 m, respectively. The geoid contribution of residual gravity and RTM contains most its power in the higher frequencies of the gravity field, and larger geoid signals are strongly correlated with the topography. Figure 10b,d,f represents the difference between the geometric geoid heights obtained by GNSS/levelling and the gravimetric geoid determined by the GBVP approach based on DIR_R6/EGM2008, RTM, and residual gravity, thus providing the spatial distribution of the vertical datum offsets. Figure 10b,d,f shows that the vertical datum offsets have east-west or north-south discrepancies that are mainly caused by systematic tilt errors and distortions in the vertical networks.
The influence of systematic tilt errors and distortions on the estimation of vertical datum offsets is evaluated according to Equation (10). The results with and without consideration of tilts and distortions are compared in the following. For readability, these two cases will be called "non-tilts" and "with-tilts" scenarios, respectively. For the three study areas, Table 4 presents the numerical results of both scenarios under the different approaches. As shown in the Table 4, considering the influence of the systematic tilt and distortions, the accuracy of the vertical offset for the USA determined by the GBVP approach was improved by 69.1% (from 29.1 cm to 9.0 cm) in terms of the standard deviation. Moreover, the accuracy of the vertical offset for Australia determined by the GBVP approach was improved by 44.6% (or from 16.8 cm to 9.3 cm) in terms of the standard deviation, and the accuracy of the vertical offset for Hong Kong determined by the GBVP approach was improved by 18.9% (from 3.7 cm to 3.0 cm) in terms of the standard deviation. Based on the above analysis, there are significant accuracy improvements on the estimation of vertical datum offsets when considering the systematic tilts and distortions. power in the higher frequencies of the gravity field, and larger geoid signals are strongly correlated with the topography. Figure 10b,d,f represents the difference between the geometric geoid heights obtained by GNSS/levelling and the gravimetric geoid determined by the GBVP approach based on DIR_R6/EGM2008, RTM, and residual gravity, thus providing the spatial distribution of the vertical datum offsets. Figure 10b,d,f shows that the vertical datum offsets have east-west or north-south discrepancies that are mainly caused by systematic tilt errors and distortions in the vertical networks. The influence of systematic tilt errors and distortions on the estimation of vertical datum offsets is evaluated according to Equation (10). The results with and without consideration of tilts and distortions are compared in the following. For readability, these two cases will be called "non-tilts" and "with-tilts" scenarios, respectively. For the three study areas, Table 4 presents the numerical results of both scenarios under the different approaches. As shown in the Table 4, considering the influence of the systematic tilt and distortions, the accuracy of the vertical offset for the USA determined by the GBVP approach was improved by 69.1% (from 29.1 cm to 9.0 cm) in terms of the standard deviation. Moreover, the accuracy of the vertical offset for Australia determined by the GBVP approach was improved by 44.6% (or from 16.8 cm to 9.3 cm) in terms of the standard deviation, and the accuracy of the vertical offset for Hong Kong determined by the GBVP approach was improved by 18.9% (from 3.7 cm to 3.0 cm) in terms of the standard deviation. Based on the above analysis, there are significant accuracy improvements on the estimation of vertical datum offsets when considering the systematic tilts and distortions.
From Table 4, considering the influence of the systematic tilts and distortions, the GBVP solution features an improvement of 8.2% compared to the DIR_R6/EGM2008 solution in terms of the standard deviation for determining the vertical offset in the USA. For Australia, the GBVP solution provides an 11.4% improvement compared to the DIR_R6/EGM2008 solution in terms of the standard  Table 4, considering the influence of the systematic tilts and distortions, the GBVP solution features an improvement of 8.2% compared to the DIR_R6/EGM2008 solution in terms of the standard deviation for determining the vertical offset in the USA. For Australia, the GBVP solution provides an 11.4% improvement compared to the DIR_R6/EGM2008 solution in terms of the standard deviation. In Hong Kong, the GBVP solution experiences a 14.3% improvement compared the DIR_R6/EGM2008 solution in terms of the standard deviation. Therefore, the accuracy of the GBVP solution is better than that using only DIR_R6/EGM2008 solution. this is because the RTM and residual gravity can compensate for the omission errors of the DIR_R6/EGM2008 model. Considering the influence of the systematic tilts and distortions, the vertical datum parameters are estimated according to Equations (10) to (12). Table 5 shows the vertical datum parameters for the USA, Australia, and Hong Kong. The geopotential value of the NAVD88 is estimated as equal to 62636861.31 ± 0.96m 2 s −2 , and the vertical offset with respect to geoid is −0.809 ± 0.090 m, which means that the NAVD88 is about 80.9 cm below the geoid. The geopotential value of the AHD is 62653852.60 ± 0.95 m 2 s −2 , and the vertical offset with respect to geoid is 0.082 ± 0.093 m, which means that the AHD is about 8.2 cm above the geoid. The geopotential value of the HKPD is 62636860.55 ± 0.29 m 2 s −2 , and the vertical offset with respect to geoid is −0.731 ± 0.030 m, which means that the HKPD is about 73.1 cm below the geoid.  [4] used the GBVP approach to determine the vertical datum parameters of South America. We chose Brazil and Argentina as representative in South America. The potential differences between the Brazil and Argentina vertical datum and global geoid are 3.79 ± 0.18 m 2 s −2 and 6.51 ± 0.49 m 2 s 2 , respectively. North America can use the same vertical datum of NAVD88. Figure 11 shows the potential differences between the North America, Australia, Brazil, Argentina, and Hong Kong vertical datums and the global geoid. By determining the potential differences between the global existing LVDs and the global geoid (as shown in Figure 11), the LVDs can be incorporated into the IHRS to ensure the global vertical datum system's consistency.
Remote Sens. 2020, 12, x 19 of 23 existing LVDs and the global geoid (as shown in Figure 11), the LVDs can be incorporated into the IHRS to ensure the global vertical datum system's consistency.

Conclusions
The unification of the global vertical datum is a key problem to be solved for geodesy for a long time, the main goal of which is to determine the vertical datum parameters for local vertical datum. In this study, we evaluate the vertical datum parameters of the USA, Australia, and Hong Kong using the GBVP approach. Meanwhile, the geoid omission errors of the GGM, the spectral accuracy of the GGM, the indirect bias term, and the systematic errors and distortions on the estimation of the vertical datum parameters are analyzed. These factors should be considered for the implementation of the GBVP approach to unify the vertical datum.
The effect of the indirect bias term is equal to 1 mm when higher degree GGMs (300 degree and order) are used. Therefore, this effect can be neglected safely when implementing the GBVP approach. To analyze the influence of the systematic tilts and distortions on the estimation of the vertical datum parameters, the parametric model in Equation (10) is used to absorb the existing systematic errors and distortions. The results show that there are significant accuracy improvements on the estimation of the vertical datum offsets when considering the systematic tilts and distortions in vertical networks. The GRACE/GOCE-based GGMs have lower degree errors and stronger geoid signal in the medium-long wavelength, but the maximum expansion degree of the GRACE/GOCEbased GGMs is limited. These models have certain omission errors. The effect of omission errors reaches the decimeter level in the USA and Australia, while the impact of the omission error is at the centimeter level in Hong Kong.
The omission errors of GOCE/GRACE-based GGMs cannot be ignored when unifying the vertical datum. To compensate for the omission errors of GOCE/GRACE-based GGMs, the EGM2008 model is used to extend the GOCE/GRACE-based GGMs to determine the combined GGMs. Because the spectral accuracy of the DIR_R6 model is significantly better than other GRACE/GOCE-based GGMs, we built a combined model, DIR_R6/EGM2008, by combining the DIR_R6 model and the EGM2008 model using the weighting method. The combined DIR_R6/EGM2008 model provides a stronger geoid signal and a lower geoid degree errors than the EGM2008 and DIR_R6 models. To obtain high-frequency geoid information caused by the Earth's topographic masses, the RTM is used to recover the short-wavelength geoid in this paper. The reference terrain model with a resolution of 5 5 ′ ′ × and a modified stokes kernel function are used to allow the contributions of the RTM and

Conclusions
The unification of the global vertical datum is a key problem to be solved for geodesy for a long time, the main goal of which is to determine the vertical datum parameters for local vertical datum. In this study, we evaluate the vertical datum parameters of the USA, Australia, and Hong Kong using the GBVP approach. Meanwhile, the geoid omission errors of the GGM, the spectral accuracy of the GGM, the indirect bias term, and the systematic errors and distortions on the estimation of the vertical datum parameters are analyzed. These factors should be considered for the implementation of the GBVP approach to unify the vertical datum.
The effect of the indirect bias term is equal to 1 mm when higher degree GGMs (300 degree and order) are used. Therefore, this effect can be neglected safely when implementing the GBVP approach. To analyze the influence of the systematic tilts and distortions on the estimation of the vertical datum parameters, the parametric model in Equation (10) is used to absorb the existing systematic errors and distortions. The results show that there are significant accuracy improvements on the estimation of the vertical datum offsets when considering the systematic tilts and distortions in vertical networks. The GRACE/GOCE-based GGMs have lower degree errors and stronger geoid signal in the medium-long wavelength, but the maximum expansion degree of the GRACE/GOCE-based GGMs is limited. These models have certain omission errors. The effect of omission errors reaches the decimeter level in the USA and Australia, while the impact of the omission error is at the centimeter level in Hong Kong.
The omission errors of GOCE/GRACE-based GGMs cannot be ignored when unifying the vertical datum. To compensate for the omission errors of GOCE/GRACE-based GGMs, the EGM2008 model is used to extend the GOCE/GRACE-based GGMs to determine the combined GGMs. Because the spectral accuracy of the DIR_R6 model is significantly better than other GRACE/GOCE-based GGMs, we built a combined model, DIR_R6/EGM2008, by combining the DIR_R6 model and the EGM2008 model using the weighting method. The combined DIR_R6/EGM2008 model provides a stronger geoid signal and a lower geoid degree errors than the EGM2008 and DIR_R6 models. To obtain high-frequency geoid information caused by the Earth's topographic masses, the RTM is used to recover the short-wavelength geoid in this paper. The reference terrain model with a resolution of 5 × 5 and a modified stokes kernel function are used to allow the contributions of the RTM and residual gravity to the geoid further compensate for the omission errors of the DIR_R6/EGM2008 model. Finally, the vertical datum parameters of the USA, Australia and Hong Kong are determined using the GBVP approach based on DIR_R6/EGM2008, RTM and residual gravity. The geopotential value of the NAVD88 is estimated as equal to 62,636,861.31 ± 0.96 m 2 s −2 , and the vertical offset with respect to the global geoid is −0.809 ± 0.090 m, which means that the NAVD88 is about 80.9 cm below the geoid. The geopotential value of the AHD is 62,653,852.60 ± 0.95 m 2 s −2 , and the vertical offset with respect to the global geoid is 0.082 ± 0.093m, which means that the AHD is about 8.2 cm above the geoid. The geopotential value of the HKPD is 62,636,860.55 ± 0.29 m 2 s −2 , and the vertical offset with respect to the global geoid is −0.731 ± 0.030 m, which means that the HKPD is about 73.1 cm below the geoid. With the vertical datum parameters, we can thus connect the local height system with the global height system to ensure the global vertical datum system's consistency.