The Refined Gravity Field Models for Height System Unification in China

Abstract

A unified height datum is essential for global geographic information resource construction, ecological environment protection, and scientific research. The goal of this paper is to derive the geopotential value for the Chinese height datum (CNHD) in order to realize the height datum unification in China. The estimation of height datum geopotential value usually depends on high-precision global gravity field models (GFMs). The satellite gravity missions of the Gravity Recovery and Climate Experiment (GRACE) and Gravity field and steady-state Ocean Circulation Exploration (GOCE) provide high-accuracy, medium–long-wavelength gravity field spectra, but satellite-only GFMs are limited to medium–long wavelengths, which will involve omission errors. To compensate for the omission errors in satellite-only GFMs, a spectral expansion approach is used to obtain the refined gravity field models using the EGM2008 (Earth Gravitational Model 2008) and residual terrain model (RTM) technique. The refined GFMs are evaluated by using high-quality GNSS/leveling data, the results show that the quasi-geoid accuracy of the refined DIR_R6_EGM2008_RTM model in China has optimal accuracy and, compared with the EGM2008 model and the DIR_R6 model, this refined model in China is improved by 9.6 cm and 21.8 cm, and the improvement ranges are 35.7% and 55.8%, respectively. Finally, the geopotential value of the Chinese height datum is estimated to be equal to 62,636,853.29 m²s⁻² with respect to the global reference level defined by W₀ = 62,636,853.4 m²s⁻² by utilizing the refined DIR_R6_EGM2008_RTM model and 1908 high-quality GNSS/leveling datapoints.

1. Introduction

With the emergence of Global Navigation Satellite System (GNSS), users can obtain consistent ellipsoidal height at global scale. The ellipsoid height relative to a given geocentric ellipsoid can be obtained quickly and accurately by using GNSS. However, the ellipsoidal height is not related to the Earth’s gravity field. The height related to the Earth’s gravity field usually refers to the orthometric or normal height, which is strictly based on the geopotential number C_P ($C_P=W_0^{LVD}-W_P$, where $W_0^{LVD}$ is the geopotential value for the local vertical datum, and W_P indicates the gravity potential for point P). The local vertical datum refers to the geoid (or quasi-geoid), which is assumed to be coincident with the local mean sea level (MSL). Importantly, even though all local height datums are related to the MSL, the vertical offsets between them may be up to 2 m at global scale [1]. This is due to the fact that the MSL presents geographical and time-dependent variations, but they are a consequence of the natural sea dynamics. Therefore, the vertical datums vary among countries or regions, affecting and restricting the sharing and exchange of global geospatial information, and the unification of the height datums has become a key point in the field of geodesy.

According to the Global Geodetic Observing System (GGOS) of the International Association of Geodesy (IAG), a global vertical datum related to Earth’s gravity field should be established. One of the main works of GGOS is to support global geometric and physical heights with centimeter accuracy within a global framework, and to unify all existing physical height systems [2,3]. A global reference system defines constants, conventions, models, and parameters as the necessary basis for the mathematical representation of geometric and physical quantities [4]. According to the IAG resolution No 1, 2015, the gravity potential value for the International Height Reference System (IHRS) is released equal to 62,636,853.4 m²s⁻² [5,6,7]. According to this resolution, the existing local vertical datum systems can be integrated into the IHRS, which will ensure the consistency of the global height datum systems. The fundamental approach for height datum unification is the geodetic boundary value problem (GBVP) method [8,9,10,11,12,13]. This method has been widely applied in areas with good coverage of surface gravity data. However, in areas where surface gravity data are poor or restricted, a feasible option for height datum unification is to use high-accuracy GFMs. The GFMs provide expected accuracy at scales from centimeter to decimeter [11]. With the development of high-precision and high-resolution global gravity flied models, the GFMs have become feasible for the unification of vertical datums [14,15,16,17,18,19,20,21,22].

The GFMs have many advantages for determining the geopotential values of vertical datums, but the results depend largely on the accuracy of the utilized GFMs. The satellite gravity missions such as GRACE [23] and GOCE [24] provide unprecedented information for the medium-wavelength gravity field [25,26], which can improve the accuracy of medium- and long-wavelength geoid. The primary goal of the GOCE satellite is to derive about 1–2 cm geoid accuracy to a target resolution of about 100 km [27]. The GRACE/GOCE-based GFMs have the medium–long wavelength information with high precision; however, these models have certain omission errors. The omission error for geoid (quasi-geoid) height reaches about 32 cm for a GFM up to degree 200, according to Kaula’s rule and the variance model of Tscherning and Rapp [28,29,30,31]. This omission error cannot be ignored for realizing the height datum unification. In addition to the omission errors of the GFMs, the spectral accuracy of the selected GFM, the uncertainties of GNSS-derived ellipsoidal heights, and the accumulation leveling errors must be considered to determine the vertical datum geopotential with a high accuracy.

The Chinese vertical datum is determined by the MSL observed by the Qingdao tide gauge station during the 27 years from 1952 to 1979. The latest first-order leveling network in China, which is based on the 1985 vertical datum and was completed and put into use in 2017, is the most accurate and practical national modern vertical control network to date [32]. It is used for transmitting normal height at national scale. The aim of the paper is to determine the vertical datum geopotential of China based on the GFM by utilizing the latest normal height results. Thus, the 1985 vertical datum in China can be connected to the IHRS, which can provide a unified height datum for the construction of global geographic information resources, ecological environmental protection, and scientific research. In order to derive more accurate height datum geopotential of China based on GFM, we utilize the spectral expansion method by augmenting the GOCE/GRACE-based GFMs with the EGM2008 model and residual terrain model (RTM) technique [33,34,35,36,37], which can reduce the omission error of the satellite-only GFMs. In addition, considering the large east–west span in China, the systematic accumulated error may occur in the long-distance leveling, and GFM have certain systematic errors between the east and west of China, which has influence on determining the geopotential value of the Chinese height datum. This paper will further analyze the systematic errors influence. Finally, we will combine the refined GFMs and GNSS/leveling to preliminarily derive the geopotential value of the Chinese height datum (CNHD) towards height system unification in China.

The structure of the manuscript is as follows. The materials and method for estimating the Chinese height datum geopotential value are introduced in Section 2. Section 3 provides the results of the Chinese height datum geopotential value and specifically focuses on (a) spectral accuracy of GFMs, (b) the omission errors of the satellite-only GFMs in China, and (c) the determination of the refined GFMs. Finally, discussion and conclusions are provided in Section 4 and Section 5, respectively.

2. Materials and Methods

2.1. Materials

This section briefly introduces the materials used in this study. The main data required are (1) GNSS/leveling data; (2) global gravity field models; and (3) topographic data.

2.1.1. GNSS/LEVELING Data

The leveling networks in China contain the first-, second-, third-, and fourth-order leveling networks. The latest first-order leveling network observations, which were completed and placed into use in 2017, are used. This first-order leveling network consists of 148 loops with a length of over 125,746.5 km. Considering the large east–west span in China, the systematic accumulated error may occur in the long-distance leveling; therefore, the leveling data is handled by the least square adjustment in which the length correction of leveling rod, the non-parallel correction of level surface, and the gravity reduction are considered [32]. The maximum error in leveling is only about 3.57 cm (about 6185 km from Qingdao leveling origin) after adjustment. Finally, 1908 high-accuracy and evenly distributed GNSS/leveling datapoints are made available by the National Geomatics Center of China to determine the vertical datum geopotential of China. The GNSS ellipsoidal coordinates are based on ITRF2014 [38], and the GNSS coordinate accuracy reaches the millimeter level; in particular, the accuracy of GNSS ellipsoidal heights is about 5 mm. The distribution for the GNSS/leveling datapoints in China is shown in Figure 1.

Because the GNSS ellipsoidal height is based on a tide-free system, to ensure that the GNSS/leveling-based quasi-geoid height represents a tide-free system, the normal heights are transformed to the tide-free system by the following equation [39]:

$$H^{TF}=H^{MF}+0.68\left(0.099-0.296{\sin }^2\phi \right)$$

(1)

where H^TF represents the tide-free normal height, H^MF is the mean tide normal height, and $\phi$ is the latitude of the GNSS/leveling points.

2.1.2. Global Gravity Field Models (GFMs)

Table 1 shows the GFMs used in this paper. The direct approach is employed to derive the DIR_R5 and DIR_R6 models using GOCE satellite gravity observations, GRACE satellite gravity observations, and satellite laser ranging (SLR). The DIR_R5 model utilizes GOCE data from November 2009–October 2013 and GRACE data from the ten-year period (2003–2012). The DIR_R6 model utilizes GOCE data from October 2009–October 2013 and GRACE data from January 2007–November 2014. The TIM_R5 and TIM_R6 models are derived by the time-wise approach using GOCE observations from November 2009–October 2013 and September 2009–October 2013, respectively. The EGM2008 model combines multi-source gravity observation data to derive a maximum degree of 2190; however, this model is complete to d/o 2159. The gravity data utilized in EGM2008 model mainly refers to global gravity database of 5′ × 5′. The surface gravity data covering China is composed of two-divisions: one is the gravity data of 5′ for construction of the EGM2008, without any restrictions, the main coverage areas include eastern, southern, and Central China; another gravity data resource is permitted to use at a resolution of 15 arc-minute, these data cover other regions except for eastern, southern, and Central China.

For the consistency of tide system, the ${\overline{C}}_{20}$ coefficient of GFM is transformed to tide-free system using the following formula [40]:

$${\overline{C}}_{20}^{TF}={\overline{C}}_{20}^{ZT}-k_{20}\cdot \langle \Delta {\overline{C}}_{20}\rangle$$

(2)

where ${\overline{C}}_{20}^{TF}$ and ${\overline{C}}_{20}^{ZT}$ are the spherical harmonic coefficient under the tide-free system and the zero-tide system, respectively, $k_{20}=0.30190$ is loading Love number, and $\langle \Delta {\overline{C}}_{20}\rangle =-1.391412\cdot {10}^{-8}$ represents the value of tidal correction.

Table 1. The details for global gravity field models. The letter “S” in the third column represents satellite data, “G” represents ground observations, “A” represents altimetry observations, and “d/o” represents degree/order.

Models	d/o	Data	Tide System	Reference
EGM2008	2160	S(Grace), G, A	Tide-free	[41]
GO_CONS_GCF_2_TIM_R5	280	S(Goce)	Tide-free	[42]
GO_CONS_GCF_2_TIM_R6	300	S(Goce)	Zero-Tide	[43]
GO_CONS_GCF_2_DIR_R5	300	S(Goce, Grace, Lageos)	Tide-free	[44]
GO_CONS_GCF_2_DIR_R6	300	S(Goce, Grace, Lageos)	Tide-free	[45]

Table 1. The details for global gravity field models. The letter “S” in the third column represents satellite data, “G” represents ground observations, “A” represents altimetry observations, and “d/o” represents degree/order.

Models	d/o	Data	Tide System	Reference
EGM2008	2160	S(Grace), G, A	Tide-free	[41]
GO_CONS_GCF_2_TIM_R5	280	S(Goce)	Tide-free	[42]
GO_CONS_GCF_2_TIM_R6	300	S(Goce)	Zero-Tide	[43]
GO_CONS_GCF_2_DIR_R5	300	S(Goce, Grace, Lageos)	Tide-free	[44]
GO_CONS_GCF_2_DIR_R6	300	S(Goce, Grace, Lageos)	Tide-free	[45]

2.1.3. Topographic Data

Topographic data are used to calculate the RTM effect and recover the high-frequency gravity signals missing in the GFMs. The RTM represents the difference between the topographic surface and the long-wavelength reference terrain. The Shuttle Radar Topographic Mission (SRTM) V4.1 data [46] with a spatial resolution of 7.5″ × 7.5″ are used to represent the land over China. For the sea area, the SRTM15_PLUS V2 topographic data [47] are used, with a spatial resolution of 15″ × 15″. The rock-equivalent topography model RET2012 [48] works as the reference topography, and the reference terrain elevations can be computed by Equation (2) in [49]. In the coastal zone, the mass density of seawater is different from that of the standard topographic mass density. To avoid the necessity of distinguishing density changes in the computational process, the rock-equivalent topography (RET) method can be used to compress the water depth to the equivalent rock height [50].

$$H^{\ast }=H\left(1-{\rho }_w/\rho \right)$$

(3)

where $H^{\ast }$ and H are the compressed water depth and the original water depth, respectively, ${\rho }_w$ denotes the density of seawater, and $\rho$ represents the standard topographic mass density.

The SRTM and SRTM15_PLUS data can be merged to obtain the unified topographic data on land and sea. The merging process is mainly divided into two steps: (1) the bicubic interpolation method is employed to interpolate the sea topographic data into a spatial resolution of 7.5″ × 7.5″, and Equation (3) is employed to process the sea water depths; (2) the land topographic data are combined with the water depth data obtained by step (1). Figure 2a shows the merged topography. The residual topographic masses (RTM elevations) are the difference between the SRTM merged data (Figure 2a) and the reference terrain. Figure 2b shows the residual terrain model elevations.

2.2. Methods for Determining the Height Datum Geopotential Value

The vertical offset $\delta H$ between the local height datum and the global geoid can be expressed as

$$\delta H=h-H-\zeta$$

(4)

where H represents the normal heights, h is the ellipsoidal heights, and $\zeta$ is the height anomalies from GFM (see Figure 3).

The height anomaly $\zeta$ can be calculated by using the satellite-only GFMs, but the resulting omission errors cannot be ignored. Therefore, the EGM2008 model can be used to extend the satellite-only GFMs to obtain the refined GFMs via the spectral expansion approach. The height anomaly ${\zeta }_{GGM}$ determined by the refined GFMs can be expressed as

$${\zeta }_{GGM}=\left(\frac{GM-G{M}_0}{r\cdot \overline{\gamma }}-\frac{W_0-U_0}{\overline{\gamma }}\right)+{\zeta }_{GOCE/GRACE}|{}_2^l+{\zeta }_{EGM2008}|{}_{l+1}^{2160}$$

(5)

where ${\zeta }_{GOCE/GRACE}|{}_2^l$ is the height anomaly determined by the satellite-only GFMs truncated to the degree l, and ${\zeta }_{EGM2008}|{}_{l+1}^{2160}$ is the height anomaly represented from degrees 201 to 2160 of EGM2008. GM₀ = 3.986005000 × 10¹⁴ m³s⁻² and U₀ = 62,636,860.8500 m²s⁻² are constants of the gravitational constant and the normal potential value of the GRS80 ellipsoid [51], respectively; r is the geocentric radial for computation point; GM represents geocentric gravitational constant used in the GFM; $\overline{\gamma }$ is the mean normal gravity; and W₀ = 62,636,853.4 m²s⁻² is the geopotential value of the global geoid [52].

The RTM technology is used to recover the short-scale signal beyond degree 2160 in Equation (5). The RTM represents the difference (residual terrains) between the topographic surface and the reference surface. The gravitational potentials of residual terrains can be expressed as follows [53]:

$$T=G\rho {\displaystyle \underset{\mathsf{\Omega }}{\int }\frac{1}{r}}d\mathsf{\Omega }$$

(6)

where T is the gravitational effect for the residual terrains, G denotes the gravitational constant, r is the distance between the attraction mass and the computation point, $\mathsf{\Omega }$ is the volume for the residual terrains, and $\rho$ is the standard topographic mass density.

The residual terrains are decomposed into a set of rectangular-prism mass. Figure 4 shows the coordinate system definition of a rectangular-prism mass body. The computation point P is the origin of this coordinate system; this means that the coordinates defining have to be transformed by a shift (see Equation (1) in [54]). To obtain the gravitational effects of a residual mass distribution, the result of Equation (6) for a single rectangular-prism can be calculated using the following equation:

$$\begin{array}{l}T\left(P\right)=G\rho {\displaystyle \underset{x_1}{\overset{x_2}{\int }}{\displaystyle \underset{y_1}{\overset{y_2}{\int }}{\displaystyle \underset{z_1}{\overset{z_2}{\int }}\frac{1}{r}}}}dxdydz=G\rho |||xy\ln \left(z+r\right)+yz\ln \left(x+r\right)+zx\ln \left(y+r\right)-\\ \begin{array}{ccc}& & \end{array}\hspace{1em}\hspace{1em}\frac{x^2}{2}{\tan }^{-1}\left(\frac{yz}{xr}\right)-\frac{y^2}{2}{\tan }^{-1}\left(\frac{xz}{yr}\right)-\frac{z^2}{2}{\tan }^{-1}\left(\frac{xz}{zr}\right)|{}_{x_1}^{x_2}|{}_{y_1}^{y_2}|{}_{z_1}^{z_2}\end{array}$$

(7)

where x₁, x₂, y₁, y₂, z₁, and z₂ describe the corner coordinates of the prism faces in Figure 4.

The total gravitational effects of the residual terrains for computation point are derived by a summation of gravitational effects in Equation (7):

$$T\left(P\right)={\displaystyle \sum _{i=1}^{M}T{\left(P\right)}_i}$$

(8)

where M represents the number for prism mass elements.

The height anomaly ${\zeta }_{RTM}$ caused by the residual mass distribution can be expressed by

$${\zeta }_{RTM}=\frac{T\left(P\right)}{{\gamma }_P}$$

(9)

where ${\gamma }_P$ is normal gravity for computation point P.

${\zeta }_{RTM}$ represents the quasi-geoid height signals of the GFM beyond degree 2160. Then, the height anomaly $\zeta$ determined by the refined GFMs can be further expressed as follows:

$$\zeta ={\zeta }_{GGM}+{\zeta }_{RTM}=\left(\frac{GM-G{M}_0}{r\overline{\gamma }}-\frac{W_0-U_0}{\overline{\gamma }}\right)+{\zeta }_{GOCE/GRACE}|{}_2^l+{\zeta }_{EGM2008}|{}_{l+1}^{2160}+{\zeta }_{RTM}$$

(10)

The global or local vertical datum is a gravitational equipotential surface. Therefore, the vertical datum offset determined by Equation (4) should theoretically be a fixed constant. However, vertical offsets contain certain discrepancies due to the influences of systematic errors and random errors. Systematic errors may contain possible systematic errors in GFM, accumulated leveling errors and the GNSS height ellipsoidal errors [55,56,57]. To reduce these random and systematic errors, the vertical datum offset is estimated by applying a planar correction parametric model. For each GNSS/leveling point P, the observation equation can be formulated as follows:

$$l_P=h_P-H_P-{\zeta }_P=\delta \overline{H}+a_1\left(B_P-B_0\right)+a_2\left(L_P-L_0\right)\cos B_P$$

(11)

where $\delta \overline{H}$ is the unknown height datum offset, $\left(B_P,L_P\right)$ are the geodetic coordinates for the computation point, $\left(B_0,L_0\right)$ are the geodetic coordinates for leveling origin in the local vertical datum zone, and a₁ and a₂ are the north–south tilt and east–west tilt, respectively. If there are n GNSS/leveling benchmarks in the local vertical datum zone, according to Equation (11), the function model can be expressed as follows:

$$\underset{l}{\underbrace{\left[\begin{array}{c}{h}_1-H_1-{\zeta }_1\\ h_2-H_2-{\zeta }_2\\ ⋮\\ h_n-H_n-{\zeta }_n\end{array}\right]}}=\underset{A}{\underbrace{\left[\begin{array}{c}1\begin{array}{cc}& \left(B_1-B_0\right)\begin{array}{cc}& \left(L_1-L_0\right)\cos B_1\end{array}\end{array}\\ 1\begin{array}{cc}& \left(B_2-B_0\right)\begin{array}{cc}& \left(L_2-L_0\right)\cos B_2\end{array}\end{array}\\ ⋮\\ 1\begin{array}{cc}& \left(B_n-B_0\right)\begin{array}{cc}& \left(L_n-L_0\right)\cos B_n\end{array}\end{array}\end{array}\right]}}\cdot \underset{x}{\underbrace{\left[\begin{array}{c}\delta \overline{H}\\ a_1\\ a_2\end{array}\right]}}$$

(12)

The unknown parameter x in Equation (12) can be further determined according to the least square adjustment of the system, denoted by

$$\mathit{x}={\left({\mathit{A}}^T\mathit{P}\mathit{A}\right)}^{-1}{\mathit{A}}^T\mathit{P}\mathit{l}$$

(13)

where $\mathit{P}={\mathit{D}}_{\mathit{l}\mathit{l}}^{-1}$ is the weight matrix, ${\mathit{D}}_{\mathit{l}\mathit{l}}$ represents the error variance–covariance matrix for the observed values. Assuming the involved terms in Equation (11) are uncorrelated to each other, the ${\mathit{D}}_{\mathit{l}\mathit{l}}$ can be specified by

$${\mathit{D}}_{\mathit{l}\mathit{l}}={\mathit{D}}_{\mathit{h}\mathit{h}}+{\mathit{D}}_{\mathit{H}\mathit{H}}+{\mathit{D}}_{\zeta \zeta }$$

(14)

where ${\mathit{D}}_{\mathit{h}\mathit{h}}$ and ${\mathit{D}}_{\mathit{H}\mathit{H}}$ are the error variance–covariance matrix of the ellipsoidal and normal heights, respectively, and ${\mathit{D}}_{\zeta \zeta }^{}$ represents the error variance–covariance matrix for quasi-geoid heights.

${\mathit{D}}_{\zeta \zeta }^{}$ might be determined from the errors of the GFMs and RTM quasi-geoid heights. Voigt and Denker [58] and Grombein et al. [17] showed that the uncertainty of the topography-implied gravity signals is at the sub-mm level; therefore, the errors of RTM quasi-geoid heights can be considered negligible. However, the error variance–covariance matrix of the GFM might generally not be available [17]. In addition, the uncertainties for ellipsoidal heights and normal heights are usually not available. Therefore, it is not possible to obtain ${\mathit{D}}_{\mathit{l}\mathit{l}}$ in practical cases. Based on the above reasons, we assume P = I in this paper, where I is the identity weight matrix.

After removing the errors effects, we will obtain vertical offsets with considering corrections by

$$\delta H_P=\left(h_P-H_P-{\zeta }_P\right)-a_1\left(B_P-B_0\right)-a_2\left(L_P-L_0\right)\cos B_P$$

(15)

The geopotential value $W_{P_0}^{LVD}$ for point P can be expressed as follows:

$$W_{P_0}^{LVD}=W_0-\delta H_P\cdot {\overline{\gamma }}_P$$

(16)

where ${\overline{\gamma }}_P$ denotes the mean normal gravity, which is computed by Equations (4)–(60) in [59].

Finally, we can derive the geopotential value $W_0^{LVD}$ of the local height datum by

$$W_0^{LVD}=W_0-\frac{{\displaystyle \sum _{P=1}^m\delta H_P\cdot {\overline{\gamma }}_P}}{m}=W_0-\delta W_0^{LVD}$$

(17)

where m is the number of GNSS/leveling benchmarks.

3. Results

3.1. Spectral Accuracy Evaluation for GFMs

According to the spherical harmonic coefficients and their formal errors, the spectral accuracy of a GFM is evaluated by using the commission errors (the degree error, cumulative degree error, difference degree error, and cumulative difference degree error) and signal-to-noise ratio (SNR). The geoid degree error can be expressed as the square root of the error degree variances, and the cumulative geoid degree error is represented as the square root of the cumulative geoid error degree variances. A relative comparison of satellite-based GFM and EGM2008 can be estimated by the difference degree error and the cumulative difference degree error. These specific calculation formulas can be found in [60].

Figure 5a shows the geoid degree errors for GFMs. We can see from the figure that the geoid degree errors of the satellite-only GFMs before about degree 200 are lower than the degree error of EGM2008. Figure 5b shows the geoid cumulative degree errors for the GFMs. The geoid cumulative degree errors for DIR_R5 and DIR_R6 models are lower than that of the EGM2008 model in the whole spectral domain. The geoid cumulative degree errors for TIM_R5 and TIM_R6 models are lower than that of the EGM2008 model before approximately degree 260. The above analysis shows that the satellite-only GFMs have higher medium–long wavelength accuracy than the EGM2008 model, and the satellite-based GFMs are characterized by noise as increases in degree.

Figure 6a indicates the SNRs for the GFMs. The SNRs of satellite-only GFMs are better than that of the EGM2008 model before about degree 200. The results show that the satellite-only GFMs have a strong geoid signal in the medium–long wavelength band, and the EGM2008 model has a strong geoid signal in the short wavelength band. Figure 6b shows the difference degree error and cumulative difference degree error of the satellite-only GFMs, taking the EGM2008 model as a reference. We can see from the figure that the geoid difference degree errors for DIR_R5 and DIR_R6 are lower than geoid difference degree errors of the TIM_R5 and TIM_R6 before about degree 100, and cumulative geoid difference degree errors of the DIR_R5 and DIR_R6 are lower than cumulative difference degree errors of TIM_R5 and TIM_R6 before about degree 150. This is because the DIR_R5 and DIR_R6 GFM take advantage of GRACE gravity data and LAGEOS satellite laser ranging (SLR) data.

3.2. The Omission Errors for Satellite-Only GFMs

To quantify the magnitudes of these omission errors of satellite-only GFMs, we assume that 200 is the general expansion degree of the satellite-only GFM. The magnitudes of the height anomaly ${\zeta }_{201\text{-}2160}^{EGM2008}$ using the degrees 201 to 2160 of EGM2008 and the height anomaly ${\zeta }_{RTM}$ inferred from the RTM can indicate the omission errors of the satellite-only GFM. The specific calculation process for ${\zeta }_{RTM}$ can be found in Section 3.3.

Figure 7a indicates the omission error of satellite-only GFM represented by ${\zeta }_{201\text{-}2160}^{EGM2008}$, Figure 7b shows the omission errors of satellite-only GFM represented by the ${\zeta }_{RTM}$, and Figure 7c shows the total omission errors obtained from the sum of ${\zeta }_{201\text{-}2160}^{EGM2008}$ and ${\zeta }_{RTM}$.

Table 2. Statistics of the omission errors for satellite-only GFMs in China. Unit: m.

Omission Errors (m)	Max	Min	Mean	Std
${\zeta }_{201\text{-}2160}^{EGM2008}$	2.885	−3.529	−0.215	0.601
${\zeta }_{RTM}$	0.264	−0.248	−0.006	0.027
${\zeta }_{201\text{-}2160}^{EGM2008}$$+{\zeta }_{RTM}$	2.892	−3.548	−0. 221	0.602

show statistics for the omission errors of satellite-only GFM in China. We can see from the figure that the omission errors in the rugged regions are larger (mid-west of mainland China) than in other regions. From Table 2, we can see that the omission errors of the satellite-only GFM in mainland China reach the decimeter level, and the largest amplitude is about 350 cm. The omission error signals represented by the RTM are centimeter level in mainland China. Therefore, the omission errors for the satellite-only GFMs must be considered for unification of China vertical datum.

3.3. The Refined GFMs Obtained by the Spectral Expansion Approach

A spectral expansion approach is used to obtain the refined gravity field models using the EGM2008 (Earth Gravitational Model 2008) and residual terrain model (RTM) technique. Because the degree errors of the satellite-only GFMs increase with the increase in degree and order, the noise starts to dominate the signals at high degree and order. Therefore, it is not possible to derive an optimal GFM by combining the maximum degree of satellite-based GFMs and EGM2008. It is not a reasonable strategy to obtain a refined GFM by combining two models directly, and the choice of the optimal degree for the combination is very important. The specific steps used to determine the optimal combination degrees are as follows: (a) the satellite-only GFMs are truncated to degree l (l = 10, 20, 30……N, where N is the maximum degree for satellite-only GFMs) as the 0~l degree of the refined GFMs, and the l~2160 degree of the refined GMFs is supplemented by the corresponding degree of the EGM2008 model; (b) the GNSS/leveling-based height anomaly is used to check the height anomaly determined by the refined GFMs; (c) the combination degree presenting the best accuracy is chosen as the optimal combination degree of the refined GFMs. Figure 8 indicates the standard deviations of the height anomalies differences between GNSS/leveling and the refined GFMs with varying combinations of degree. As seen in Figure 8, the accuracy of the combined GFMs is basically consistent with that of the EGM2008 model before about degree 90. From degree 90 onwards, there are obvious differences in accuracy. We select the optimal combination degree when the standard deviations minimized. Figure 8 shows that 230, 240, 220, and 240 are the optimal combination degrees for obtaining the refined GFMs by combining the TIM_R5, TIM_R6, DIR_R5, and DIR_R6 models, respectively, with the EGM2008 model.

The RTM technology can further be utilized to obtain higher frequency gravity field signals of the refined GFMs. The RTM gravitational potential is not harmonic when the computation point is below the reference elevation surface. In order to solve the problem of non-harmonic, a harmonic correction can usually be performed by downward continuation through a Bouguer plate [33,61]. Thus, the harmonic correction for quasi-geoid (or geoid) can be considered as zero; however, harmonic correction for gravity can be considered as $4\pi G\rho \mathsf{\Delta }H$ ($\mathsf{\Delta }H$ represents the height difference between the topography and reference terrain). The choice of the RTM integration radius is crucial. For RTM quasi-geoid height, an integration radius of ~200 km is suitable [49]. Therefore, an integration radius of ~200 km is used herein to determine the RTM quasi-geoid heights. According to unified topographic data (in Figure 2), we can obtain RTM quasi-geoid height by using Equations (7)–(9). Figure 9 shows the quasi-geoid contribution of the RTM with a resolution of 7.5” × 7.5”, with maximum, minimum, mean, and standard deviation values of 0.301 m, −0.280 m, 0.001 m, and 0.020 m, respectively.

Table 3. Statistics of the height anomaly differences between GNSS/leveling and the GFMs. Unit: m.

Models	Max	Min	Mean	STD
EGM2008	1.162	−2.079	0.061	0.269
DIR_R6	2.624	−2.475	−0.114	0.391
DIR_R5	2.536	−2.507	−0.112	0.396
TIM_R6	2.578	−2.457	−0.107	0.391
TIM_R5	2.549	−2.432	−0.114	0.400
TIM_R5_EGM2008	1.044	−1.804	0.045	0.185
TIM_R6_EGM2008	1.099	−1.806	0.050	0.184
DIR_R5_EGM2008	1.007	−1.689	0.048	0.183
DIR_R6_EGM2008	1.072	−1.787	0.044	0.181
EGM2008_RTM	1.256	−1.981	0.068	0.261
TIM_R5_EGM2008_RTM	1.042	−1.803	0.053	0.178
TIM_R6_EGM2008_RTM	1.097	−1.806	0.058	0.177
DIR_R5_EGM2008_RTM	1.005	−1.688	0.056	0.176
DIR_R6_EGM2008_RTM	1.069	−1.787	0.052	0.173

shows the quasi-geoid accuracy statistics of different GFMs in China. From Table 3, without considering the influence of the RTM on the quasi-geoid height, the quasi-geoid accuracies of the refined GFMs in mainland China are better than 18.5 cm. Among them, the DIR_R6_EGM2008 model has the best accuracy, at 18.1 cm. Compared with the EGM2008 model, the quasi-geoid accuracy of the DIR_R6_EGM2008 model is improved by 8.8 cm. On the other hand, these results also indicate that the quasi-geoid medium-long wavelength accuracy of EGM2008 model in China is poor. Considering the influence of the RTM quasi-geoid, the mainland China quasi-geoid accuracies determined by the refined GFMs are better than 17.8 cm, and the DIR_R6_EGM2008_RTM model has the optimal quasi-geoid accuracy with a standard deviation of 17.3 cm. Compared with the EGM2008 model, the quasi-geoid accuracy of the DIR_R6_EGM2008_RTM model is improved by 9.6 cm. Compared with the TIM_R5, TIM_R6, DIR_R5, and DIR_R6 models, the quasi-geoid accuracy of the DIR_R6_EGM2008_RTM model is improved by 0.400 m, 0.391 m, 0.396 m, and 0.391 m to 0.178 m, 0.177 m, 0.176 m, and 0.173 m, respectively.

To validate our results of the refined GFMs, the EIGEN-6C4 [62], GECO [63], SGG-UGM-1 [64], SGG-UGM-2 [65], XGM2016 [66], and XGM2019 [67] models are used. These six higher-degree GFMs further improve the medium–long wavelength by adding GOCE data.

Table 4. Statistics of the height anomaly differences between GNSS/leveling and six higher-degree GFMs. Unit: m.

Models	Max	Min	Mean	STD
EIGEN-6C4	1.007	−1.696	0.048	0.187
GECO	1.579	−1.703	0.041	0.223
SGG-UGM−1	1.003	−1.671	0.052	0.194
SGG-UGM−2	1.003	−1.704	0.051	0.191
XGM2019	1.705	−1.737	0.081	0.213
XGM2016	1.016	−1.757	−0.020	0.214

summarizes the statistics for height anomaly differences between GNSS/leveling and six higher-degree GFMs. These models in Table 4 perform better than the EGM2008 model in China, which can assume the largest impact from the contribution of GOCE solution. The EIGEN-6C4 outperforms GECO, SGG-UGM-1, and SGG-UGM-2 models in China, which is mainly due to the differences in use of GOCE data. Comparing Table 3 and Table 4, we can find that the refined GFMs outperform EIGEN-6C4, GECO, SGG-UGM-1, XGM2019, XGM2016, and SGG-UGM-2 in mainland China; the major improvement of these models can be attributed to the GOCE data and topography signals.

3.4. Determination for the Geopotential Value of Chinese Height Datum

The geopotential W₀ = 62,636,853.4 m²s⁻² adopted as reference level for the IHRS is used herein. Based on the GFM, the vertical datum geopotential value for China is determined. Thus, the vertical datum in China is connected into the IHRS.

In the analysis presented in Section 3.3, we conclude that the DIR_R6_EGM2008_RTM model has optimal quasi-geoid accuracy in China. Therefore, we choose this model to derive the Chinese height datum geopotential value. According to Equation (4), we can obtain vertical offset values of each GNSS/leveling point. Figure 10a represents the height anomaly differences between the GNSS/leveling and the DIR_R6_EGM2008_RTM model in China, thus providing the spatial distribution for the vertical datum offsets. However, it can be found from Figure 10a that the discrepancies exist in mainland China. The vertical offsets have certain discrepancies due to the influence of various factors. Among them, systematic error is a major contributor. Therefore, the planar corrections surface is used further through Equation (11). Then, we can derive vertical offsets after applying the correction surfaces by Equation (15). Figure 10b shows correction surface at GNSS-leveling points. Figure 10c shows vertical offsets after applying the correction surfaces. It can be seen from Figure 10b that there are certain east–west systematic errors in China, and the maximum error is about 12.2 cm. Systematic errors mainly contain errors in the GFM and leveling, The maximum first-order leveling error in China is only about 3.57 cm after height networks adjustment [32]. Therefore, we can conclude that the major systematic effect contributor to the estimation of the height datum geopotential values of China mainly comes from GFM error, and the systematic errors from GFM are more obvious in western China.

The systematic errors influence on the determination of the Chinese height datum geopotential values is evaluated. The results with and without consideration of systematic error effect are compared, as shown in

Table 5. Statistics of the without-planar corrections and with-planar corrections scenarios for the height datum geopotential value in China based on DIR_R6_EGM2008_RTM model. Unit: m²s⁻².

Model	Scenarios	Max	Min	Mean
DIR_R6_EGM2008_RTM	Without-planar corrections	62,636,870.86	62,636,842.93	62,636,852.89 ± 1.75
	With-planar corrections	62,636,870.88	62,636,843.31	62,636,853.29 ± 1.69

. Table 5 presents the numerical results of both scenarios based on the DIR_R6_EGM2008_RTM model. As shown in Table 5, we can see that there is a minor improvement of 0.06 m²s⁻² for DIR_R6_EGM2008_RTM model when considering the planar corrections. Finally, the geopotential value for the Chinese height datum is derived to be equal to 62,636,853.29 m²s⁻² based on the DIR_R6_EGM2008_RTM model, considering planar corrections.

4. Discussion

The satellite-only GFMs have high spectral accuracy and strong geoid signals. How-ever, the maximum expansion degree of these satellite-based models is limited and there are certain omission errors as a result of the gravity field attenuation at the height of satellite orbit. From Figure 7 and Table 2, we can see that the omission errors of the satellite-only GFM in mainland China reach the decimeter level, and the largest amplitude is about 350 cm. The omission error signals represented by the RTM are centimeter level in mainland China. Therefore, the omission errors for the satellite-only GFMs must be considered for unification of China vertical datum. The satellite-only GFMs have higher spectral accuracy and stronger geoid signals at medium–long wavelength, and EGM2008 model has stronger geoid signals at short wavelength. Therefore, combining the GRACE/GOCE-based GFMs and EGM2008 model to obtain refined GFMs is a feasible strategy.

We combine the GRACE/GOCE-based GFMs and EGM2008 model to derived refined GFMs. The accuracy trend of the refined GFMs in Figure 8 mainly depends on the EGM008 model, the satellite-only GFMs, and GNSS/leveling resources. The GNSS/leveling have a good accuracy and quality; especially, the systematic errors in leveling have been greatly weakened by height network adjustment. Therefore, the results presented in Figure 8 mainly reflect the accuracy of the combined GFMs in China. The accuracy of the combined GFMs relies heavily on the improvement of the satellite-only models. Gruber and Willberg [56] demonstrated that 80% accuracy improvement for high-resolution GFMs compared with EGM2008 reveal the contribution of GOCE solution to medium wavelengths. Therefore, high-quality GNSS/leveling resources and the satellite-only GFMs increase the quality and reliability of the combined GFMs in this study. In the combination process, the satellite-only GFM using GOCE, GRACE, and LAGEOS laser ranging data should be used, which can better meet the quasi-geoid medium–long wavelength signal.

The RTM technology can further be utilized to obtain higher-frequency gravity field signals of the refined GFMs. In the spectral expansion process, the spatial resolution of the refined GFMs obtained by combining the satellite-only GFMs and EGM2008 model is 5′. Because the RET2012 reference topographic model has the same resolution as the refined GFMs, the reference model serves as a high-pass filter that can filter out low-frequency features from the SRTM data. As a result, the residual terrain height can imply gravity field signals at shorter scales than the spatial resolution of the refined GFMs and further compensate for the omission error of the refined GFMs. It can be seen from Figure 9 that larger quasi-geoid signals are strongly correlated with topography. However, it should be noted that the RTM quasi-geoid signals may contain some possible implications due to density anomaly and the uncertainties in the harmonic correction [48].

Finally, we use the refined GFMs to preliminarily derive the geopotential value of the Chinese height datum by utilizing the latest normal height results. The refined GGMs provide obvious accuracy improvement advantage and provide guidance for developing a high-accuracy quasi-geoid in China. Prior to this work, He et al. [20] determined the geopotential value of China’s height datum as 62,636,853.47 m²s⁻², which fits well with our value of 62,636,853.29 m²s⁻², considering the standard deviation. The geopotential value of Chinese height datum is determined by utilizing the latest national first-order leveling network and GNSS points in this paper, which are representative to some extent.

We can expect that the refined GFMs can provide a guidance for determining the quasi-geoid or the geopotential value of the Chinese height datum. However, the combination of the satellite-only GFM with the EGM2008 in this study is by a pure complementation of the spherical harmonic coefficients at a specific degree; however, such a procedure might cause a spectral gap between both models, which is different from rigorous combination that is carried out on the basis of the normal equations and covariance by a least-squares. In the next step, the rigorous combination will be considered to derive the refined GFMs. In addition, the well-distributed surface gravity data for determining the geopotential value of the Chinese height datum also play a crucial role and need to be considered further.

5. Conclusions

In this paper, we used a spectral expansion to derive the refined GFMs by combining the EGM2008 and the satellite-only GFMs. The results show that 230, 240, 220, and 240 are the optimal combination degrees for determining the refined GFMs by combining the TIM_R5, TIM _R6, DIR_R5, and DIR_R6 models with the EGM2008 model, respectively. To consider the influence of higher-frequency gravity field signals caused by topography, the RTM is utilized to further compensate for the omission errors in the refined GFMs. The mainland China quasi-geoid accuracies determined by the refined GFMs are better than 17.8 cm, and the DIR_R6_EGM2008_RTM model has an optimal quasi-geoid accuracy with of 17.3 cm. To validate our results of the refined GFMs, the EIGEN-6C4, GECO, SGG-UGM-1, and SGG-UGM-2, XGM2019, and XGM2016 models were used. The results show that the refined GFMs outperform EIGEN-6C4, GECO, SGG-UGM-1, XGM2019, XGM2016, and SGG-UGM-2 in mainland China; the major improvement of these models can be attributed to the GOCE data and topography signals.

The systematic error effects for determining the geopotential value of the Chinese height datum are considered. The results show that the refined GFMs show minor improvements when considering the planar corrections. The major systematic effect contributor for determining the Chinese height datum geopotential values mainly comes from GFM error. Finally, the Chinese height datum geopotential value is derived to be equal to 62,636,853.29 m²s⁻² based on the refined DIR_R6_EGM2008_RTM model when considering planar corrections.