《1. Introduction》

1. Introduction

High-resolution Earth’s gravity field model can be used for high-precision geoid determination [1,2], the unification of global height systems [3], the determination of dynamic sea surface topography [4], and exploring Earth’s interior structure [5]. With the advent of the new generation satellite gravity missions (Challenging Minisatellite Payload (CHAMP) [6], Gravity Recovery and Climate Experiment (GRACE) [7], and Gravity Field and SteadyState Ocean Circulation Explorer (GOCE) [8]), the accuracy of long to medium wavelength signals was improved greatly [9–16]. Meanwhile, surface gravity anomaly data, which is constructed from terrestrial gravity, satellite altimetry, airborne gravimetry, or fill-in gravity anomalies computed by residual terrain model (RTM) forward modeling, provides high precision short wavelength information [17]. Thus, high-resolution gravity field model can be obtained by combination of the gravity signals from satellite gravity data, satellite altimetry data, surface gravity data, airborne gravity data, shipborne gravity data, and terrain model.

The high-resolution gravity field models incorporating gravity data from the dedicated satellite missions published on the International Center for Global Earth Models (ICGEM) website (http://icgem.gfz-potsdam.de/tom_longtime) include Earth Gravitational Model (EGM)2008 [18], European improved gravity model of the earth by new techniques (EIGEN)-6C4 [19], GECO [20], gravity observation combination (GOCO) model GOCO05c [21], Experimental Gravity Field Model (XGM)2016 [22], and SGGUGM-1 [23]. Referring to the ICGEM website, the main attributes of these models are shown in Table 1. EIGEN-6C4 is a representative model of the EIGEN-series, which share almost the same calculation strategy and data sources. The new-generation satellite gravity missions contributed greatly to these models, and the accuracy of the long-to-medium-wavelength parts of the models has been improved substantially. EGM2008, which is currently the most frequently used gravity field model, is constructed with possibly the best global 5'×5' data set of gravity anomaly data from terrestrial observations, satellite altimetry, and fill-in gravity anomalies from RTM forward modeling and the GRACE normal equation (NEQ) of the Institute of Geodesy and Geoinformation of the University of Bonn (ITG)-GRACE03S satellite-only model. However, it does not contain any GOCE observations. A later model EIGEN-6C4 incorporates GOCE observations and Laser Geodynamics Satellite (LAGEOS) observations. However, the surface gravity anomalies on land areas contained in EIGEN-6C4 are taken from EGM2008. Compared to the block-diagonal NEQs of high-degree coefficients in EGM2008 and EIGEN-6C4 modeling, GOCO05c and XGM2016 are developed based on the combination of full NEQ systems up to full resolution, and therefore they use regionally varying weighting based on the varying quality of the terrestrial/ altimetry data. Moreover, compared to EIGEN-6C4, the gravity anomaly data used in GOCO05c and XGM2016 is independent of EGM2008. Both GECO and SGG-UGM-1 are calculated through the improvement of EGM2008 with GOCE data, and they are calculated using EGM2008-derived global gravity anomaly data and GOCE-only NEQs [22–24]. According to the validation results from Liang et al. [23], EGM2008, EIGEN-6C4, and SGG-UGM-1 have consistent accuracy in United States and the GOCE-related models (e.g., EIGEN-6C4 and SGG-UGM-1) have better performances in China. SGG-UGM-2 is different from SGG-UGM-1 in three main aspects: the use of ellipsoidal harmonic functions, the update of gravity anomaly data in marine areas, and the employment of the GRACE NEQ.

《Table 1》

Table 1 Main characteristics of the released high-resolution gravity field models.

a S is for satellite (e.g., GRACE, GOCE, and Laser Geodynamics Satellite (LAGEOS)), A is for altimetry, and G is for ground data (e.g., terrestrial, shipborne, and airborne measurements).

Among the above mentioned models, only EGM2008 is constructed based on the ellipsoidal harmonic functions. However, as the figure of the Earth can be closely approximated by an oblate ellipsoid of revolution, the errors induced by spherical approximation are bigger than those caused by ellipsoidal approximation. Thus, for modeling Ohio States University (OSU)91 [25], EGM96 [26], and EGM2008 [18], gravity anomalies were reduced onto the earth’s reference ellipsoid. In this situation, ellipsoidal harmonic analysis is more suitable than spherical harmonic analysis [27]. Hotine [28] and Jekeli [29,30] proposed the renormalized Legendre function of the second kind and derived the mutual transformation formulas between ellipsoidal and spherical harmonic coefficients. This transformation method was later numerically investigated by Gleason [31], and is called “Jekeli’s transformation” in this paper. In addition, Sebera et al. [32] extended the direct computation of the Legendre functions up to second derivatives and minimized the number of required recurrences by the hypergeometric transformation.

With respect to their applications in gravity field modeling, ellipsoidal harmonic functions are always employed in the following manner: The ellipsoidal harmonic coefficients are first calculated through ellipsoidal harmonic analysis with gravity anomalies on the reference ellipsoid, and then the ellipsoidal coefficients are transformed to spherical harmonic coefficients of disturbing potential (). This method of determining gravity field models is called the ellipsoidal harmonic analysis and coefficients transformation (EHA-CT) method for brevity in this paper. Implementations of the EHA-CT method can differ from each other, for example, the ellipsoidal harmonic analysis can be fulfilled by either numerical integration or the least-squares method. Rapp and Pavlis [33] derived a formula for the computation of 

directly from area-mean gravity anomalies on the ellipsoid based on Jekeli’s transformation. This formula was later used in computing gravity field models, such as OSU89 [33], OSU91 [25], EGM96 [26], IGG97LB [34], and MOD99 [35]. In computing the EGM2008 [18], at first, ellipsoidal harmonic coefficients () are calculated by ellipsoidal harmonic analysis with the least-squares method. Then Jekeli’s transformation is used to transform  to .

Although the EHA-CT method has been widely used in gravity field determination, there is still lack of thorough introduction and review of its different applications. In this paper, we will first give a short introduction to this method and we will then introduce our implementation and highlight the differences relative to other implementations. The weighting and sampling theory of Driscoll and Healy [36] is introduced in the ellipsoidal harmonic analysis. Moreover, its employment in other studies will also be reviewed. For example, we point out that 1/(n – 2k – 1) in the formula given by Rapp and Pavlis [33] should be 1/(n – 1). Although this might be a typo, it is still worth pointing out the problem to avoid its misuse. Recovering the gravity field model from gravity anomalies on the ellipsoid is also described as the ellipsoidal geodetic boundary value problem (GBVP). The EHA-CT method gives a solution to this problem; meanwhile, other alternative methods can be found in Refs. [36–43].

Gravity data from the vast ocean areas, which account for nearly 71% of the earth’s area, is necessary for modeling a highresolution gravity field model. Fortunately, radar altimeter data from more and more altimetry satellites can be used for recovering marine gravity anomalies. The released altimetry data include Geosat GM/ERM (17 d), ERS-1/GM (168 d), ERS/ERM (35 d), T/P/T/P Tandem (10 d), Jason-1/ERM (10 d), Envisat (35 d/30 d), Jason-2/ ERM (10 d), Jason-1/GM (406 d), CryoSat-2 (369 d), SARAL/AltiKa ERM (35 d), HY-2A (14 d), Jason-2/GM, and SARAL/AltiKa GM. The notation “d” in the brackets after the mission’s name means day, which indicates the repetition period for each altimeter mission. By combining these multiple sources of altimetry data, grid marine gravity anomalies in the latitude range of ±80.738 with a 1'×1' spatial resolution can be recovered based on either a numerical analysis method [44,45] or least squares collocation [46,47] with geoid height as the intermediate variable. The EGM2008- derived gravity anomalies were used to fill in the ocean area to determine SGG-UGM-1. In this paper, the selected altimetry data shown above are used to recover 1'×1' spatial resolution marine gravity anomaly data. By combination with EGM2008-derived data for the rest of the area, the global surface gravity anomaly data is formed for the development of the new model SGG-UGM-2. In addition, the GRACE satellite mission was in orbit for over 15 years and provided valuable data for recovering the long-wavelength part of the gravity field. Institute of Theoretical Geodesy and Satellite Geodesy (ITSG)-Grace2018 [48], which consists of constrained daily solutions, a high-resolution static field, and unconstrained monthly solutions, is the latest time series of the ITSG series model at the Institute of Geodesy in Graz University of Technical. The authors [48] provided the NEQ system of ITSG-Grace2018 on the ftp sever of their institute (ftp://ftp.tugraz.at/outgoing/ITSG/ GRACE/ITSG-Grace2018/). The NEQ of the static gravity field from GRACE is used in modeling SGG-UGM-2. Note that we intend to continuously develop SGG-UGM-series models and release them as alternatives to users on the ICGEM website. SGG-UGM-2 will also be available there.

The paper is divided into 6 sections. First, the principles of the EHA-CT method and the derivations of the formulas are given in Section 2. The derived discrete integral formulas and least-squares formulas of the EHA-CT method are evaluated in Section 3. The data processing strategies of forming the GOCE satellite NEQ, the determination of marine gravity anomalies, the combination of the NEQs of satellite data and gravity anomalies, and the scheme of determining SGG-UGM-2 are given in Section 4. The SGG-UGM-2 model is validated in Section 5. The conclusions are given in Section 6.

《2. Methodology》

2. Methodology

《2.1. Ellipsoidal harmonic analysis and coefficient transformation method》

2.1. Ellipsoidal harmonic analysis and coefficient transformation method

The arbitrary harmonic function f, which satisfies Laplace’s equation, can be represented by the ellipsoidal harmonic expansions in its harmonic domain as follows [30,32]:

where (u, δ, λ) are ellipsoidal-harmonic coordinates [1]; E is the linear eccentricity of the reference ellipsoid; b is the semiminor axis of the reference ellipsoid;  are the fully normalized surface spherical harmonic functions:

are the fully normalized associated Legendre functions of the first kind [1];  are the fully normalized ellipsoidal harmonic coefficients; n and m are the degree and order of the ellipsoidal harmonic functions and their relevant coefficients; and  is Jekeli’s renormalized function [30].

Thus, the harmonic function  is represented by the ellipsoidal harmonic expansions as follows [31,33]:

where r is the geocentric radius, is the gravity anomaly defined by Eqs. (2–287) in Ref. [1], is the product of r and , are the “ellipsoidal harmonic” coefficients of , and is the semimajor axis of the reference ellipsoid. The quotation mark on the term ellipsoidal harmonic here indicates that are not rigorous the ellipsoidal harmonic coefficients of , as there exists a scale factor in the expansion.

When gravity anomalies refer to the surface of the reference ellipsoid, Eq. (2) is simplified to the following equation:

where []E means that the data refer to the reference ellipsoid.

The ellipsoidal harmonic coefficients can be estimated if we have gravity anomaly data on the ellipsoid surface. Then there are two more steps for estimating the coefficients as follows [33]:

Based on the scheme in Eq. (4), first, are transformed to the spherical harmonic coefficients of [] E , . Then, are transformed to . The detailed transformation equations are given in the next section.

The approach described above is named as the EHA-CT method in this paper for abbreviation. Moreover, if the integration method is employed in the ellipsoidal harmonic analysis, the EHA-CT method is called the integral EHA-CT method. Otherwise, if the least-squares method is used instead, the method is called the least-squares EHA-CT method. In addition, the discrete observations are either grid area-mean values or point values. We will discuss the related formulas in detail in Sections 2.1.1 and 2.1.2.

2.1.1. Integral formulas of the EHA-CT method

According to Eq. (3), the integral formula for estimating the coefficients from gravity anomalies on the reference ellipsoid is obtained based on the orthogonality of the surface spherical harmonics as follows:

The coefficients are transformed to using Jekeli’s transformation as follows [30]:

For the explicit expression of , please refer to Jekeli [30]. The relation between coefficients and was given by Gleason [31] as follows:

where GM is the geocentric gravitational constant.

Thus, by combining Eqs. (5), (6), and (7), we obtain the following:

where = .

Even though Eq. (8) is rigorous, it is not possible to use it in gravity field determination as we only have discretized gravity anomalies in reality. Thus, the discretized form of Eq. (8) is needed. The discretization depends on the gravity anomaly values, in the form of point values or area means, and layout of the grid, such as an equal-angular grid.

For area-mean values in equal-angular grid, the discretization of Eq. (8) is shown as follows:

where i and j denote a cell residing on the i-th (i = 0, ..., N – 1) latitude belt ( "row”) and the j-th ( j = 0, ..., 2N – 1) meridional sector ( "column”) in a global equiangular grid composed of N rows by 2N columns on the reference ellipsoid; i_max and j_max are determined by i_max = N – 1 and j_max = 2N – 1, respectively; [ri ] E is the radius at the center of a block on the i-th latitude belt; are the area-mean gravity anomalies in a block; the area-mean value of [] E in the block is approximated as ; and is the integral of Legendre function. and  are computed as:  and . The detailed derivation of this equation can be seen in the Appendix A.

For the grid point gravity anomalies in and equal-angular grid, we have the following:

where and are the grid width in the latitudinal and meridional direction, respectively; is the grid point gravity anomaly in the block and, generally, is the value at the center of the block; and  is defined as = .

However, as the Legendre base functions lose their orthogonality in the discrete case, both Eq. (9) and Eq. (10) have discretization errors and will yield the approximate coefficients . Moreover, in the case of area-mean values, in the derivation of Eq. (9) there exists an assumption that the gravity anomalies are constant in a grid and equal to the area-mean value of the grid. However, in most cases, the observations will have fluctuations in each block and it will be less smooth. Therefore, the coefficients from the equation are smoothed values. As a refinement, the smoothing factor suggested by Colombo [49] is employed to de-smooth the coefficients. After applying the smoothing factor, the de-smoothed form of Eq. (9) is as follows:

Even though the smoothing factor is applied, Eq. (11) is still an approximated formula as the orthogonality problem of the discretized integral of Legendre functions still exists. To our limited knowledge, this problem has not been solved. Fortunately, in the case of point values, the orthogonality problem of discretized Legendre functions can be solved using certain weighting and sampling schemes. Two popular quadrature algorithms, the Gauss– Legendre (GL) quadrature method (or the second-order Neumann method) [50] and the quadrature following Driscoll/Healy (DH) method [36], were developed. In the DH method, the data should be distributed within a [N ×N ] (nmax = N/2 – 1) grid, where nmax is the maximum spherical harmonic degree of the coefficients, with a latitude parallel sampling of = 180°/N and a meridian sampling of = 360°/N, or on a [N ×2N] equiangular grid with = = 180°/N. The weight accounting for the orthogonality problem of the Legendre functions in the DH method is as follows:

By applying the weighting scheme in Eq. (12) to Eq. (10), we can derive the formula for estimating the coefficients from the grid point gravity anomalies on the ellipsoid as follows:

Similar to the DH method, the GL quadrature method can also be used to solve the orthogonality problem [50]. However, the grid in the GL method is an irregular grid with equidistant sampling along latitude parallels and variable sampling along meridians. On the meridians, the grid points are the zero-crossing of . A recent employment of the GL quadrature method was in the numerical evaluation of geodetic convolution integrals by Hirt et al. [51]. For more details, please refer to Refs. [50,51]. Considering that the data used for the gravity field determination are commonly provided in the equal angular grid, we prefer the DH method in this paper.

Eq. (11) and Eq. (13) can be used to compute coefficients using the grid area-mean and point gravity anomalies on the ellipsoid, respectively. To our limited knowledge, the DH method has not been employed in gravity field determination with gravity anomalies on the ellipsoid surface; therefore, we think that the derived formula Eq. (13) is new and helpful to readers. Note that Rapp and Pavlis [33] also derived the counterpart of Eq. (11), which is Eq. (20) in Ref. [33]. The basic principles of the derivation of these two formulas are exactly the same. A more detailed derivation of Eq. (11) is shown in the Appendix A. However, we find that 1/(n – 1) in Eq. (11) is different from 1/(n –2k – 1) in its counterpart in Ref. [33]. From our derivation procedure above, it is very clear that 1/(n – ) or 1/(n –2k – 1) is produced by the relation between the coefficients and (see Eq. (7)), which is only related to the degree n.

2.1.2. Least-squares formulas of the EHA-CT method

Based on Eq. (3), the surface ellipsoidal harmonic expansion for the grid area-mean and point gravity anomalies and are shown as follows:

and

According to Eq. (14), based on a standard Gauss–Markov model, the functional and statistical models for estimating the coefficients from the observations or on the reference ellipsoid are defined as follows:

where y is the vector of the observations or  ; A is the design matrix, xe is the vector of to be estimated, ε is the observation error vector; D{y } is the error variance–covariance matrix; P and Q are the weight matrix and its inverse, respectively; and is the unit weight variance.

Based on Eq. (16), the coefficients can be estimated by the least-squares technique. Then, according to the strategy in Eq. (4), the estimated coefficients are transformed to by Eqs. (6) and (7). However, when solving Eq. (16) by the least-squares method, the dimension of the normal matrix is extremely huge for determining a high-resolution gravity field. Thus, the computational work volume is extremely large, for example, in the determination of the gravity field model up to degree 2160 with the dimensions of the normal matrix of 4 669 917 × 4 669 917. The block-diagonal least-squares (BDLS) method can be used to reduce the computational tasks [18,19]. If the gravity anomalies referring to the ellipsoidal surface can satisfy the requirements for a blockdiagonal normal matrix [49], the normal matrix can be inversed by the BDLS method [18]. The grid equal angular gravity anomalies on the ellipsoidal surface can meet all the requirements except for the requirements for data weights. The data weights should be longitude independent and symmetric with respect to the equator [49]. However, considering real data conditions, it is always hard to meet such requirement for data weights. Therefore, some approximation is always employed. Liang et al. [52] analyzed the impact of the approximation of matrix P with unit matrix I, and the approximation error was moderate. Therefore, in this study, when forming the normal matrix, all the gravity anomalies have equal weight.

Meanwhile, by combining Eqs. (6), (7), and (14), we obtain the explicit relation between the grid area-mean and point gravity anomalies on the ellipsoid and the coefficients as follows:

Eq. (17) can be used to directly set up the observation equation with respect to the coefficients , as in Eq. (16) and to simulate the grid area-mean and point gravity anomalies on the reference ellipsoid. Note that ellipsoidal corrections are required if the EHA-CT method is used with raw terrestial gravity data.

The least-squares form of EHA-CT method was also employed in the determination of EGM2008 [18]. Its calculation scheme is summarized as follows:

Comparing Eq. (4) and Eq. (18), the processing strategy in this paper is different from the one in Ref. [18]. In Pavlis et al. [18], the observation coefficients are at first determined, and then are transformed to using Jekeli’s transformation. This data processing strategy is also rigorous and the corresponding formulas in Ref. [18] used to determine EGM2008 are correct.

《2.2. Combination of surface gravity anomalies and satellite gravity observations》

2.2. Combination of surface gravity anomalies and satellite gravity observations

Different data sets, such as gravity anomalies, GOCE observations and GRACE observations, can be combined using the leastsquares method in recovering a gravity field model. In the leastsquares method, when each data set is assumed to be uncorrelated with other data, the combined solution from multiple data sets is given in the following equation:

where Ai denotes the design matrices that establish the relationship between the unknown coefficients and the i data set; Pi is the weight matrix of the i data set; li is the observation vector; is the initial variance component with respect to the i data set; Ni and Ui are the normal matrix and vector with respect to the i data set; and denotes the relative weight of the i data set.

As the initial variance components are not accurate enough for the combination of the NEQs, the variance components are estimated iteratively by variance component estimation (VCE) method [53]. In the VCE method, the variance component is determined by the following [53]:

where k denotes the k-th iteration; represents the residual vector; and rk,j represents the redundancy. For more details, please refer to Ref. [53].

However, when using the VCE method on the NEQ level, on one hand, the weight sums of the residuals are not always provided when the NEQs are provided by others. On the other hand, if stochastic models of noise with respect to the normal matrices are different and not accurate enough, the classical VCE method will cause nonoptimal results [54]. On the contrary, the weighting scheme proposed in Ref. [55] is used here to estimate the variance components. In this weighting scheme, the relative weights of the NEQs are derived according to the noise of their corresponding solutions. The weights are derived by VCE on the solution level and then applied to the NEQs. The variance components and weights are given by the following [55]:

where is the combined solution in the k-th iteration; is the solution corresponding to the i NEQ; and Ncoef,j is the number of coefficients in the i NEQ.

In this way, the relative weights of the different data sets are estimated iteratively and the final combined solution is also determined in the iteration. In modeling the SGG-UGM-1 model in Liang et al. [23], since both of the unknown parameters of the observation equations for surface gravity anomalies and satellite gravity data are the spherical harmonic coefficients, the combination can be easily done based on the least squares approach.

However, in this study, the ellipsoidal harmonic base functions are employed when constructing the NEQ of surface gravity anomalies based on Eq. (16). Therefore, the unknown parameters in the NEQs are the ellipsoidal coefficients , while the unknowns in the satellite NEQs are always the spherical harmonic coefficients. Before combination, the unknown parameters should first be unified. We can transform the satellite observation equation with the spherical harmonic coefficients into the equation with the ellipsoidal harmonic coefficients .

Similar to Eq. (16), the observation equation of the satellite data for the spherical harmonic coefficients is written as follows:

where the subscript “s” indicates the quantities with respect to the satellite observations, and is the vector of to be estimated.

The transformation between the coefficients and is expressed in matrix form as follows:

where T se is a matrix computed by Eqs. (6) and (7); and is a vector of .

Upon substituting Eq. (23) into Eq. (22), we obtain the following:

Under the least squares criterion, the NEQ is given by the following:

Eq. (25) provides the relation of the NEQs with different parameters, thus we can transform the satellite observation equation with the spherical harmonic coefficients into the equation with the ellipsoidal harmonic coefficients . Then, different NEQs can be unified and combined. The combined estimated coefficients can finally also be converted to coefficients .

《3. Evaluation of the derived formulas of the EHA-CT method》

3. Evaluation of the derived formulas of the EHA-CT method

The main goal of this section is to evaluate the discrete integral formulas and least-squares formulas of EHA-CT derived in Section 2, and the discrete integral formula (Eq. (20)) in Rapp and Pavlis [33] using numerical experiments. On one hand, we want to ensure that the formulas used for the determination of SGGUGM-2 are correct. On the other hand, as shown in Section 2.1.1, the two discrete integral formulas used for estimating the spherical coefficnets with grid area mean gravity anomaly data on the ellipsoid, Eq. (11) in the paper and Eq. (20) in Rapp and Pavlis [33], are different in terms of the items 1/(n – 1) and 1/(n – 2k – 1), respectively. It is interesting to validate these items. The numerical test is based on a close-loop test. First, gravity anomaly data on the reference ellipsoid are simulated with an initial set of coefficients, and then the formulas are used to recover the input coefficients. The error of the estimated coefficients with respect to the input coefficients reflects the accuracy of the formulas.

Based on Eq. (17), we simulate two different datasets of gravity anomalies on the Geodetic Reference System 1980 (GRS80) reference ellipsoid [56] using the EGM2008 up to degree 2190 and order 2159. One dataset is the grid area-mean gravity anomalies, and another is the point gravity anomalies. The spatial resolution is 2'×2' .

First, we use Eq. (11) to recover the geopotential model (named Model1) from the simulated 2'×2' grid area-mean gravity anomalies. The degree error root mean square (RMS) of Model1 compared to EGM2008 is shown in Fig. 1 in red. As analyzed in Section 2.1.1, there are still discretization errors in Eq. (11), although the smoothing factors are employed. From the figure, the influences of the discretization errors in Eq. (11) on the coefficients of Model1 run up to the magnitude of 10–11 (except when otherwise specified, all the coefficient errors in this paper mean the absolute error compared to the "true” input EGM2008 model coefficients) for the area-mean gravity anomalies, which cannot be ignored if we want to fully recover the input coefficients.

Then, Model2 and Model3 are calculated using Eq. (13) with DH weights and Eq. (10) without DH weights, respectively from the simulated 2'×2' grid point gravity anomalies. The degree error RMS values of Model2 and Model3 are also shown in Fig. 1 in blue and dark gray, respectively. According to Fig. 1, the degree error RMS of Model2 is lower than the magnitude of 10–17 throughout the whole frequency band (2–2160), the level of which can be considered as the effect of computer truncation error. The degree error RMS of Model3 is very large and shows huge fluctuation in even and odd degrees. The accuracy of Model2 is far higher than those of the Model1 and Model3. Thus Eq. (13) is deemed to be accurate enough to recover the input coefficients, which also demonstrates the effect of the sampling weights used in it.

《Fig. 1》

Fig. 1. The degree error RMS of Model1 (red), Model2 (blue), Model3 (dark gray), Model4 (magenta), and Model5 (green) compared to EGM2008. (a) Degrees 2–2160; (b) degrees 2–200. The degree error RMS of the EGM2008 model coefficients is also shown here as the dashed black line.

Moreover, with the simulated 2'×2' grid point gravity anomalies, Model4 is estimated based on the least-squares method as shown in Section 2.1.2. The coefficients are first estimated by the least-squares technique with the observation equation implied by Eq. (14). The BDLS method and OpenMP [57] parallel computing technique are used to deal with the huge computing task, considering the coefficients to be estimated (e.g. degree and order up to 2160 corresponding to 4 669 917 parameters). Then, the estimated coefficients are transformed to by Eqs. (6) and (7). The degree error RMS of Model4 is also shown in Fig. 1 in magenta. The coefficient errors of Model4 are also small enough to be neglected, which shows that the least-squares formula of the EHA-CT method in this paper have enough accuracy to be used for gravity field determination. Meanwhile, compared to the integral method, the least squares method has advantages in its flexibility in concerning data covariance and its flexibility in the combination of different data types. Therefore, the least squares method is used for estimating SGG-UGM-2. Note that we have applied the methods used for Model3 and Model4 to the grid area-mean gravity anomalies. The results are very similar to those from grid point values.

Eq. (11) in this paper and Eq. (20) in Ref. [33] are both discretized formulas and differ in terms of two items 1/(n – 1) and 1/(n – 2k – 1). Therefore, both contain discretization error, which causes difficulty in distinguishing whether the error is caused by the formula differences or from discretization error. Therefore, it is hard to evaluate Eq. (11) and its counterpart in Ref. [33] using the equations themselves. However, as shown above, Eq. (13), which is a similar formula to Eq. (11) used for point values, can fully recover the coefficients from point gravity anomalies. The equations share the 1/(n –1). Analogously, if Eq. (20) in Rapp and Pavlis [33] is correct, we can get its corresponding formula for computing from point gravity anomalies by replacing the factor 1/(n – 1) in Eq. (13) with 1/(n –2k – 1):

Therefore, the errors caused by Eq. (26) reflect the influences of the 1/(n –2k – 1) item in Eq. (20) in Rapp and Pavlis [33]. This is the "trick” that we use to validate the equation. Based on Eq. (26), we recover a model up to degree and order 2160 (named Model5) from the simulated 2'×2' grid point gravity anomalies. The degree error RMS of Model5 is also shown in Fig. 1 in green. We can see that the degree error RMS of the long wavelength part of Model5 is relatively large and reaches a magnitude of 10–9 around degree 5. The degree error RMS values greater than 10–11 are mainly located at the low degrees (n < 50), which cannot be ignored. Thus it is inferred that 1/(n –2k – 1) in Eq. (20) in Ref. [33] is wrong although it might be a typo.

To further analyze the formulas derived in this paper and Rapp and Pavlis [33], the gravity anomalies and geoid errors on the reference ellipsoid of Model2 and Model5 are computed, and given in Table 2. The spatial distributions of the model-derived geoid errors of Model2 and Model5 are shown in Fig. 2. From Table 2, the RMS values of the gravity anomalies and geoid errors of Model5 up to degree and order 2160 are 0.18 mGal (1 mGal = 1×105 m·s–2 ) and 11 cm, respectively, which are far larger than those of Model2. For Model5, the maximum geoid error is 3.5 cm for the degree range 100–2160 and 1.9 cm for the degree range 200–2160. And, the maximum error of the gravity anomalies is 1.4 mGal for the degree range 100–2160 and 1.2 mGal for the degree range 200–2160. According to Fig. 2, Model5 shows large and systematic geoid errors, while the geoid errors of Model2 are far less than those of Model5. These results reflect the error level caused by the item 1/(n –2k – 1) in Eq. (26).

《Table 2》

Table 2 Statistics of the global gravity anomalies and geoid errors of Model2 and Model5 compared to EGM2008.

STD: standard deviation.

《Fig. 2》

Fig. 2. Spatial distribution of the model-derived geoid errors of (a) Model2 and (b) Model5 up to degree 2160.

《4. Computation of the high-resolution gravity field model SGGUGM-2》

4. Computation of the high-resolution gravity field model SGGUGM-2

In this paper, we combine altimetry data, satellite gravity data, and surface gravity anomalies to compute the high-resolution gravity field model SGG-UGM-2 up to degree and order 2160. The data processing strategies of different observations (satellite gravity, satellite altimetry data) will be discussed briefly in the following sections. Moreover, the strategy of combining the NEQs of the satellite observations and gravity anomalies is given.

《4.1. Forming the NEQs of GOCE and GRACE satellites》

4.1. Forming the NEQs of GOCE and GRACE satellites

To construct the NEQ of the GOCE satellite, the released GOCE’s EGG_NOM_2 and SST_PSO_2 products are used here [58]. The EGG_NOM_2 product mainly includes gravity gradient tensor (GGT) observations in gradiometer reference frame (GRF), the attitude quaternions EGG_IAQ_2 used for the transformation from inertial reference frame (IRF) to GRF, and the common-mode accelerations EGG_CCD_2C. The SST_PSO_2 product includes the kinematic orbits SST_PKI_2 (PKI orbits), the variance–covariance information SST_PCV_2 of the precise PKI orbits, reduceddynamic orbits SST_PRD_2, and the quaternions SST_PRM_2 used for the transformation from earth-fixed reference frame (EFRF) to IRF. The data period of the EGG_NOM_2 products is approximately 2.5 years starting from 1st of November, 2009. The data period of SST_PKI_2 product is approximately eight months starting from 1st of November, 2009. The sampling interval of all kinds of observations is 1 s. We only used the diagonal components (Vxx , Vyy , Vzz ) with high accuracy of the GGT to form the satellite observation NEQ [24]. The maximum degrees of the recovered model from satellite gravity gradient (SGG) and satellite-to-satellite in high-low mode (SST-hl) data (PKI orbits) are 220 and 130, respectively.

Based on the data described above, the key data processing strategies in forming the NEQ of the GOCE satellite are as follows:

(1) All SGG and SST-hl data are preprocessed, such as the data interpolation, outlier detection, coordinate system transformation, and epoch unification.

(2) The NEQ of SGG is formed independently based on the direct method [24]. A bandpass auto regressive moving-average (ARMA) filter with the pass-band of 5–41mHz is applied to both sides of the linear observation equation to deal with the colored noise in SGG data [59]. The maximum frequency fmax = 41 mHz of the pass-band approximately corresponds to the maximum degree of 220 of the geopotential model based on the formula fmax = Nmax=Tr , where Tr = 5383 s is one satellite orbital revolution [60].

(3) The NEQ of SST-hl is formed independently by the pointwise acceleration approach, and the observation residuals are computed [60–62]. The accelerations of satellite motion are derived from the kinematic satellite positions based on the extended differentiation filter (EDF5) technique with [62].

(4) The NEQs obtained from SGG observations and SST-hl observations are combined according to their variance components. For more details about the weighting strategies, please refer to Xu et al. [24].

This is a brief description of forming the NEQ of the GOCE satellite. The constructed GOCE NEQ is also the basis for determining the GOCE-only model GOSG01S and the high resolution gravity field model SGG-UGM-1. For more detailed description of the data processing in recovering this GOCE-only satellite gravity model, please refer to Xu et al. [24].

As mentioned in the introduction, the NEQ of ITSG-Grace2018 [48] is used as the NEQ of the static gravity field from GRACE in SGG-UGM-2, and the GRACE satellite observations are not used or processed here. The NEQ of ITSG-Grace2018 in SINEX [63] format is converted to the format defined in our software; thus, we can use it directly in the computation.

《4.2. Global marine gravity anomaly recovery》

4.2. Global marine gravity anomaly recovery

For recovering global marine gravity anomalies, multi-satellite altimeter datasets including Geosat, ERS-1, Envisat, T/P, Jason-1, CryoSat-2 and SARAL/AltiKa are collected and used [45]. The used satellite altimetry data sets and corresponding record numbers during the preprocessing procedure are collected in Table 3. The specific cycle number and time span are not investigated in constructing the global marine gravity model. It is well known that the geoid heights and vertical deflections derived from satellite altimeter measurements provide major input information to calculate marine gravity anomalies. In addition, the process of calculating the vertical deflection from sea surface heights can effectively restrain the radial orbit error and other long-wavelength corrections. The numerical-analytical method leads to reasonable skipping of the complicated crossover adjustment procedure, and yields a reliable accuracy according to previous numerical tests using the same altimeter measurements [64].

《Table 3》

Table 3 Data used in the computation of 1'×1' resolution marine gravity data.

Consequently, we first obtained the information on the vertical deflection from multi-satellite altimeter datasets through a series of joint processing procedures and recovered the desired marine gravity anomalies by the numerical-analytical method [45]. First, the raw waveforms from different altimeter missions were fitted and corrected using a two-pass waveform retracker [45] and resampled along profiles to a reasonable rate, aiming at enhancing both the accuracy and density of the available measurements. Second, the obtained measurements were transformed to sea surface heights using correction items provided in the standard products to constrain the corresponding effects of both path delay and the geophysical environment. Afterward, the along-track sea surface height gradients were calculated, while the along-track gradients of the EGM2008 were also interpolated for preliminary verification to detect outliers. Considering that the high frequency noise was amplified during the difference procedure, we used Parks– McClellan low-pass filters to obtain along-track filtered sea surface height gradients data. Then, the DOT2008A and EGM2008 models were selected respectively to interpolate and subtract from along-track observations to remove the effects of the sea surface topography and geoid height. The along-track residual vertical deflections were computed according to the velocity formulas of ground tracks. The relationship between the along-track residual vertical deflections and the two-dimensional components of residual deflections can be established as equations at each grid point.

Based on the above procedures, the directional components of the residual vertical deflection at gridding points were calculated. Then, the residual gravity anomalies were calculated according to the relationship formula between the gravity anomalies and vertical deflections. At last, a 1'×1' resolution marine gravity anomaly dataset was then computed after restoring the reference model. We compared this dataset with DTU10, DTU13, and SS V23.1 using three kinds of ship-measured data provided by the National Geophysical Data Center (NGDC) in Table 4 for three situations, which represent a shallow water area, non-shallow water area, and open sea, respectively. The results showed that our marine gravity anomaly dataset has a higher accuracy over non-shallow water area and open ocean areas compared to recently published models such as DTU10, DTU13, and SS V23.1, although the significant difference is quite close.

《Table 4》

Table 4 Validation information using NGDC shipboard gravity data over typical areas.

《4.3. Combination of the NEQs of GOCE and GRACE satellite observations and surface gravity anomalies》

4.3. Combination of the NEQs of GOCE and GRACE satellite observations and surface gravity anomalies

Since satellite and surface gravity anomaly data have different spectral sensitivities to gravity field, the method by which to properly make use of the gravity signal implied from satellite observations and surface gravity anomalies is very important for obtaining an optimal high-resolution gravity field model. In this study, we assume that the three kinds of observations, GRACE observations, GOCE observations, and surface gravity anomalies, are uncorrelated. Therefore, these observations can be easily combined on the NEQ level. We use the NEQ system of the ITSG-Grace2018 [48] model provided by the authors instead of processing the original GRACE observations. The combination method follows the degree-dependent NEQ combination technique that was used for the computation of EGM96 [26] and the EIGEN-series highresolution model [19]. The strategy of combining the NEQs of satellite gravity data and the gravity anomalies is shown in Fig. 3.

《Fig. 3》

Fig. 3. The schematic diagram of combining the NEQs of GRACE satellite, GOCE satellite, and surface gravity anomalies.

The maximum degrees of the NEQs of GRACE satellite, GOCE satellite, and surface gravity anomaly data are 200, 220, and 2159, respectively. As we have gravity anomalies on the ellipsoid, the block-diagonal form can be achieved according to the required rules [49], which is helpful for reducing the computation scale. As it is hard to directly combine the full satellite NEQs and the blockdiagonal ones from gravity anomalies, the combination of the NEQs is divided into two parts as shown in Fig. 3, namely, the low-degree part and the high-degree part. On one hand, in the high-degree part the coefficients of 251–2159 degrees are estimated by solving the block-diagonal form NEQs only from the surface gravity anomalies (including the marine gravity data derived by satellite altimetry). On the other hand, in the low-degree part, the coefficients of 2–250 degrees are determined by combining the NEQs of satellite observations and surface gravity anomalies. Degree 250 is called the "transition degree” here. The premise of such separation is that the two parts of the coefficients are uncorrelated; however, this is not true. In reality, incorporation via the satellite NEQ should not only change the coefficients in low degrees but also those in high degrees. This simplification will lead to some approximation. An appropriate selection of the "transition degree” should be well considered when combining the NEQs. On the positive side if this degree is larger, more coefficients are estimated with the combination of satellite observations and gravity anomalies. On the other side, when the degree is increased, the dimension of the combined normal matrix in the low-degree part will become larger, which can result in a much greater computational task. Considering the balance between the computational task and accuracy requirement, the "transition degree” is set as 250 in this paper. This can also help obtain a smooth transition from the full combined to the block-diagonal solution.

Considering that satellite observations are more sensitive to the long wavelength part of the gravity field compared to surface gravity anomalies, the signals of gravity anomalies corresponding to coefficients of less than degree 101 are removed. We select this special degree according to the geoid degree errors of the satellite solutions (ITSG-Grace2018 and GOSG01S) and gravity anomaly solution (EGM2008). As shown in Fig. 4, the geoid errors of EGM2008 reach maximum values around degree 100, while after degree 100 the errors decrease. Because our surface gravity anomalies on land are derived from EGM2008, it is reasonable to select the special degree band based on the performance of EGM2008. In addition, the residual gravity anomalies are used to form the NEQ of the coefficients of 101–250 degrees, which is combined with the NEQs of the GOCE and GRACE satellite observations.

《Fig. 4》

Fig. 4. Geoid degree errors of the EGM2008, GOSG01S, and ITSG-Grace2018 models.

《4.4. Computation of the SGG-UGM-2 model》

4.4. Computation of the SGG-UGM-2 model

The scheme of determining the SGG-UGM-2 is shown in Fig. 5. The detailed steps and data processing strategies of modeling SGGUGM-2 are as follows:

《Fig. 5》

Fig. 5. The scheme of computing the SGG-UGM-2 model.

(1) Using the processing strategy described in Section 4.1, the NEQ of the GOCE satellite from SGG and SST-hl observations with respect to coefficients up to degree and order 220 is constructed. Meanwhile, a GOCE-only model named GOSG01S is achieved after solving this NEQ. The NEQ of ITSG-Grace2018 from the authors is converted to the format defined in our software.

(2) The parameters of these two satellite NEQs are the spherical harmonic coefficients . As we intend to unify the NEQs for easy combination, the NEQs are transformed into NEQs with parameters of ellipsoidal harmonic coefficients using the transformation matrix shown in Section 2.2.

(3) Based on the method discussed in detail in Section 4.2, we estimate the 5'×5' marine grid gravity anomalies on the GRS80 reference ellipsoid using satellite altimetry data.

(4) The 5'×5' continental grid point gravity anomalies are computed from the EGM2008 on the reference ellipsoid using Eq. (27). The parameters GM and for determining SGG-UGM-2 are 3.986004415 ×1014 m3·s–2 and 6 378 136.3 m, respectively, which are same as those in the satellite NEQs. A global set of the gravity anomaly data set is formed by the combination of the recovered marine gravity anomalies in the previous step and EGM2008-derived gravity anomalies on land areas. Not only in this step but also throughout the whole study, the numerical problem of computing high degree and order normalized associated Legendre functions is solved using the method proposed by [65].

(5) With the constructed global gravity anomalies , the NEQ with parameters of the coefficients of 2–2159 degrees is formed in block-diagonal form. The BDLS method and the OpenMP parallel computing technique are employed here. Meanwhile, a set of coefficients of 2–2159 degrees is computed by solving this NEQ. To fulfill the requirement of the BDLS method, the unit matrix is selected as the weight matrix and the variance component applied to the NEQ is set as (5 mGal)2 , which is approximately derived from the accumulated gravity anomaly error of EGM2008 up to degree 2190.

(6) The residual gravity anomalies on the ellipsoid are computed based on the recovered coefficients in step (5) with Eq. (28); thus, only contain the signals of coefficients from degree 101 to 250. Afterwards, the block-diagonal form NEQ is constructed using as follows:

(7) The satellite NEQs with as the parameters in step (2) and the NEQ formed in step (6) are accumulated to get the combined NEQ system up to degree 250. All the relative weights among the NEQs are set to 1.0 originally. Then, the ellipsoidal harmonic coefficients of degrees 2–250 are estimated based on the strict least-squares method from the combined NEQ.

(8) After getting the combined solution in the low-degree part, the relative weight between the NEQs are updated using the VCE method at the solution level as shown in Section 2.2. Afterwards, the step (7) and step (8) are carried out using the updated weights until they converge. In this study, we only need four iterations, and the weights of the NEQs of surface gravity anomaly data, GOCE, and GRACE converge to 1.0, 1.9, and 1.3, respectively.

(9) Then, the final set of ellipsoidal harmonic coefficients in the degrees 2–250 is obtained in the iteration process and the remaining part, the coefficients in the degree range 251–2159 are those estimated in step (5).

(10) The ellipsoidal harmonic coefficients of 2–2159 degree determined in step (9) are transformed into the spherical harmonic coefficients of the degree 2–2190 and order 0–2159 by Eqs. (6) and (7). Finally, after adding the normal potential coefficients of GRS80 ellipsoid, the high resolution gravity field model SGG-UGM2 is obtained.

《5. Accuracy analysis of the SGG-UGM-2 model》

5. Accuracy analysis of the SGG-UGM-2 model

《5.1. Comparison with EGM2008 in the frequency and spatial domains》

5.1. Comparison with EGM2008 in the frequency and spatial domains

To analyze the accuracy of the SGG-UGM-2 model, we compute the degree RMS of the coefficient differences between our model and EGM2008, which are shown in Fig. 6 in red. The degree RMSs of the coefficient differences between the other two highresolution models (EIGEN-6C4 and GECO) and EGM2008 are also shown in Fig. 6. We also plot the spectra of the coefficient differences between the three models (SGG-UGM-2, EIGEN-6C4, and GECO) and EGM2008 in Fig. 7. According to Fig. 6 and Fig. 7, the three models are very close to each other below degree 160, especially SGG-UGM-2 and EIGEN-6C4, because all of them contain GOCE data. The signal differences of all three models begin to diverge from each other above degree 160, which is caused by the different surface gravity datasets and the different combination methods used in their modeling. We use the newly derived marine gravity anomalies, while GECO only uses the EGM2008-derived marine gravity anomalies and the gravity anomalies used for EIGEN-6C4 are very close to those used for EGM2008. This also results in the model coefficients of EIGEN-6C4 and GECO after degree 360 being more close to those of EGM2008 than those of SGG-UGM-2. Moreover, the coefficients of the GECO model are exactly same as those of EGM2008 after degree 360.

《Fig. 6》

Fig. 6. The degree RMS of the coefficient differences between the three models (SGG-UGM-2, EIGEN-6C4, and GECO) and EGM2008.

《Fig. 7》

Fig. 7. Spectra of the absolute value of the coefficient differences (represented by common logarithm lg x ) between the three models (a) SGG-UGM-2, (b) EIGEN-6C4, and (c) GECO and EGM2008.

The differences of the model-derived gravity anomalies between SGG-UGM-2 and EGM2008 are computed and shown in Fig. 8. Similarly, the differences between SGG-UGM-2 and EIGEN6C4 are also shown in Fig. 8. According to Fig. 8, the large differences between SGG-UGM-2 and EGM2008 are located at areas where there are no gravity data or only sparse gravity data used for compiling EGM2008, such as the Tibetan Plateau, South America, central Africa, and Antarctica. This indicates the contribution of the GOCE data to the SGG-UGM-2 and EIGEN-6C4 models. Moreover, there are also large differences around coast lines, which might reflect differences in the marine gravity anomaly data used between SGG-UGM-2 and EGM2008. The differences around coast lines between SGG-UGM-2 and EIGEN-6C4 have similar characteristics, because the marine gravity anomalies used for modeling EIGEN-6C4 are very close to EGM2008-derived gravity anomalies.

《Fig. 8》

Fig. 8. Spatial distribution of the model-derived gravity anomaly differences between SGG-UGM-2 and other models (a) EGM2008 and (b) EIGEN-6C4.

《5.2. Validation using global positioning system/Leveling data in the United States and the mainland of China》

5.2. Validation using global positioning system/Leveling data in the United States and the mainland of China

For analyzing the accuracy of the SGG-UGM-2 model, we first use 649 global positioning system (GPS)/Leveling points in the mainland of China [66] and 6169 GPS/leveling points in the United States [67] to validate the gravity field models. GPS/leveling data in the mainland of China refer to quasi-geoidal heights, while GPS/ leveling data in the United States refer to geoidal heights. Therefore, in the validation we use the models to compute the geoidal heights in the United States and the quasi-geoidal heights in China on the GPS/leveling points. Note that the GPS/leveling data in both China and the United States as well as the gravity field models to be validated use the tide-free system [67–69]. He et al. [69] showed that there is about a 70 cm tilt in the west–east direction in GPS/levelling datasets in the United States [69], while the west– east tilt of the data in the mainland of China is approximately 9 cm. The statistical results of the full differences between the quasigeoidal/geoidal heights of the SGG-UGM-2 model and the GPS/ leveling data in the United States and the mainland of China are given in Table 5 and Table 6. Note that the differences in Table 5 and Table 6 refer to full differences without removing any deterministic model. To compare these results with recently released high-resolution models, the validation results of EGM2008, EIGEN-6C4, SGG-UGM-1, and GECO are also given in the tables. Moreover, to validate these models, histograms of the differences with respect to the GPS/leveling data sets in the United States and the mainland of China are shown in Fig. 9 and Fig. 10.

《Table 5》

Table 5 Statistical results of comparison with GPS/leveling data in the United States (6169 points) (unit: m).

《Table 6》

Table 6 Statistical results of comparison with GPS/leveling data in the mainland of China (649 points) (unit: m).

According to Table 5, Table 6, Fig. 9, and Fig. 10, in the United States, the accuracies of EGM2008, EIGEN-6C4, SGG-UGM-1, SGGUGM-2, and GECO are very close to each other as suggested by the error STDs, and the error STDs of the models differ by less than 7 mm. The histograms corresponding to these models are very similar to each other, although the dispersion degree of the difference distribution is slightly high, which might indicate that the GPS/ leveling data in the United States have different quality levels in different regions. SGG-UGM-2 has the best performance and EGM2008 has the worst performance as suggested by their STDs, which are 0.277 and 0.284 m, respectively. However, in the mainland of China, the models behave inversely, and the error STDs of the models range from 0.157 to 0.240 m. The models (SGG-UGM-1, SGG-UGM-2, GECO, and EIGEN-6C4) including GOCE data have very similar accuracies, and EIGEN-6C4 performs the best with the error STD of 0.157 m. However, EGM2008 also behaves the worst in the mainland of China with the error STD of 0.240 m. Moreover, its difference relative to the best model is much bigger than the difference relative to the best model in the United States. The histograms corresponding to these models including GOCE data are also very similar to each other, as in the case of in the United States. Meanwhile, the dispersion degree of the difference distribution of EGM2008 is the highest, and the histogram shows obvious differences relative to other models. This may be caused by only low accuracy gravity data or no data in the mainland of China being available to EGM2008 developers. The other models overcome this problem by including GOCE data, which improves the accuracy of long wavelength signals in the mainland of China; therefore, their performances over China are better than that of EGM2008. In addition, both SGG-UGM-1 and SGG-UGM-2 have promising accuracies in the United States and the mainland of China, and can be regarded as improvements over EGM2008 because of combining GOCE and GRACE satellite observations and satellite altimetry data. However, due to the contribution of the new GRACE NEQ system and the new marine gravity anomalies, SGG-UGM-2 has a better performance than its predecessor SGG-UGM-1 in both the mainland of China and the United States. Because EIGEN-6C4 and SGG-UGM-2 share similar combination methods and input data, their performances in both the mainland of China and the United States are similar. EIGEN-6C4 uses more satellite gravity data (e.g., LAGEOS) than SGG-UGM-2. The relative weights of the surface gravity data and satellite data are determined with the modified VCE method in SGG-UGM-2, while they are empirically determined by the model validation result in EIGEN-6C4 [12,19].

《Fig. 9》

Fig. 9. The histograms of the differences with respect to the GPS/leveling data sets in the United States for (a) EGM2008, (b) EIGEN-6C4, (c) SGG-UGM-1, (d) SGG-UGM-2, and (e) GECO.

《Fig. 10》

Fig. 10. The histograms of the differences with respect to the GPS/leveling data sets in the mainland of China for (a) EGM2008, (b) EIGEN-6C4, (c) SGG-UGM-1, (d) SGG-UGM2, and (e) GECO.

Moreover, the models are validated through the variogram analysis of the differences with respect to the GPS/leveling data sets. Each variogram represents the variance of the differences between the model and the GPS control data set for pairs of points as a function of the lag distance. The computational method of the empirical variograms refers to Ref. [70]. The empirical variograms for the models in the mainland of China and the United States are shown in Fig. 11. Following Refs. [71,72], the term “gammavariance” representing the variance of the differences at a given lag distance is used here. The empirical variograms of EGM2008 and GECO show obvious differences from those of other three models in the mainland of China. EGM2008 has the higher gammavariance in both areas, while SGG-UGM-1, SGG-UGM-2, and EIGEN-6C4 have similar gammavariances in both areas. In both the mainland of China and the United States, the EIGEN-6C4 model has almost the lowest gammavariance, especially in cases of long distances, which indicates that long wave-length signal of EIGEN-6C4 performs the best.

《Fig. 11》

Fig. 11. The empirical variograms of the differences with respect to the GPS/leveling data sets in (a) the mainland of China and (b) United States for EGM2008, EIGEN-6C4, GECO, SGG-UGM-1, and SGG-UGM-2.

《5.3. Validation using GPS/leveling data in Qingdao and Taiwan》

5.3. Validation using GPS/leveling data in Qingdao and Taiwan

To validate the accuracy of SGG-UGM-2 in coastal regions and islands, we compare the model-derived quasi-geoidal/geoidal heights with GPS/leveling data in two coastal areas in China, Qingdao and Taiwan. The GPS/leveling data in Qingdao and Taiwan contain 152 points and 88 points, respectively. The statistical results of the full differences between the quasi-geoidal/geoidal heights of the SGG-UGM-2 model and the GPS/leveling data in Qingdao and Taiwan are shown in Table 7 and Table 8. Meanwhile, the histograms of the differences with respect to the GPS/leveling data sets in Qingdao and Taiwan are shown in Fig. 12 and Fig. 13. The empirical variograms for all the models in Qingdao and Taiwan are shown in Fig. 14, which represents the variance of the differences between the model and the GPS control data set for pairs of points as a function of the lag distance in Qingdao and Taiwan. As the height reference frames in Qingdao and Taiwan are the normal height and orthometric height respectively, the model-derived heights in Qingdao and Taiwan are quasi-geoidal heights and geoidal heights, respectively.

《Fig. 12》

Fig. 12. The histograms of the differences with respect to the GPS/leveling data sets in Qingdao for (a) EGM2008, (b) EIGEN-6C4, (c) SGG-UGM-1, (d) SGG-UGM-2, and (e) GECO.

《Fig. 13》

Fig. 13. The histograms of the differences with respect to the GPS/leveling data sets in Taiwan for (a) EGM2008, (b) EIGEN-6C4, (c) SGG-UGM-1, (d) SGG-UGM-2, and (e) GECO.

《Fig. 14》

Fig. 14. The empirical variograms of the differences with respect to the GPS/leveling data sets in (a) Qingdao and (b) Taiwan for EGM2008, EIGEN-6C4, GECO, SGG-UGM-1, and SGG-UGM-2.

According to Table 7 and Table 8, the error STD of SGG-UGM-2 is smaller than that of SGG-UGM-1 in Qingdao, but bigger than that of SGG-UGM-1 in Taiwan. On one hand, the statistics of the differences in Qingdao indicates that the newly included marine gravity anomaly data improved the high-resolution model in the coastal regions and islands. On the other hand, this is not true in Taiwan. The reason for this situation might be that surface gravity data with very good quality in Taiwan has been used for modeling EGM2008. Therefore, the newly included satellite data and marine gravity data in SGG-UGM-2 do not improve its accuracy in Taiwan. The same situation happened for SGG-UGM-1, EIGEN-6C4, and GECO. The error STDs of these models in Taiwan are also larger than those of EGM2008. In Qingdao, the histograms corresponding to all the models show very similar patterns, which is consistent with the statistical results in Table 7. However, in Taiwan, the histograms corresponding to all the models show obvious differences, which should be caused by the complex topography in Taiwan and its surroundings. According to Fig. 14, EGM2008 has the best performance at all distances, as indicated by the results shown in Table 7 and Table 8, especially at distances from 80 to 140 km in Qingdao, which approximately correspond to the degrees from 140 to 250. This frequency band can be greatly contributed by surface gravity anomalies, which has been proven by the situation in which GOCE-related models show no obvious improvement at areas with well covered surface gravity data, such as oceanic areas and the United States.

《Table 7》

Table 7 Statistical results of comparison with GPS/leveling data in Qingdao (152 points) (unit: m).

《Table 8》

Table 8 Statistical results of comparison with GPS/leveling data in Taiwan (88 points) (unit: m).

《6. Conclusions》

6. Conclusions

In this paper, we introduce the EHA-CT method and give its implementation strategies. The related formulas in the implementation strategies for computing the spherical harmonic coefficients from the grid area-mean and point gravity anomalies on the ellipsoid are derived. The DH weighting and sampling theory [36] is introduced for the ellipsoidal harmonic analysis. A review of the implementation of the EHA-CT method in Rapp and Pavlis [33] shows that the formula Eq. (20) in Ref. [33] contains a wrong item, which might be a typo. The simulation experimental results show that the formula Eq. (20) in Ref. [33] causes large errors in the long wavelength part of the gravity field model, while the corresponding formula derived in the paper is rigorous.

Moreover, based on the derived least-squares formulas of the EHA-CT method, we develop a new 5'×5' spatial resolution gravity field model SGG-UGM-2 up to degree 2190 and order 2159 by combining GOCE SGG and SST-hl observations, the ITSGGrace2018 NEQ system, marine gravity anomalies recovered from satellite altimetry data, and EGM2008-derived continental gravity data. The new SGG-UGM-2 model has a promising performance in the GPS/leveling validation and error analysis compared to EGM2008 in the frequency and spatial domains. The GPS/leveling data in China and the United States are used to validate the model SGG-UGM-2, together with EIGEN-6C4, SGG-UGM-1, GECO, and EGM2008. SGG-UGM-2 shows the best performance in the United States, as indicated by the statistics of the differences between model-derived quasi-geoidal/geoidal heights and GPS/leveling data, and their histograms and empirical variograms. Due to the contribution of the new GRACE NEQ and the new marine gravity anomalies, SGG-UGM-2 has a slightly better performance than that of its predecessor SGG-UGM-1 in both the mainland of China, the United States, and the coastal city Qingdao of China. This indicates that the methods used for developing SGG-UGM-2 are valid and can be used for developing future SGG-UGM series with available independent terrestrial gravity datasets (e.g. the mainland of China). In addition, the accuracy of the new model SGG-UGM-2 indicates that this model will provide an alternative for users.

《Acknowledgements》

Acknowledgements

We appreciate the help from Torsten Mayer-Gürr and Andreas Kvas for providing us the NEQ system of the ITSG-Grace2018 model. This research was financially supported by the National Natural Science Foundation of China (41574019 and 41774020), the German Academic Exchange Service (DAAD) Thematic Network Project (57421148), the Major Project of High-Resolution Earth Observation System, and Science Fund for Creative Research Groups of the National Natural Science Foundation of China (41721003), the Fundamental Research Funds for the Central Universities (N170103009). We also thank the editor and the anonymous reviewers for their constructive remarks that helped us to improve the quality of the manuscript.

《Compliance with ethics guidelines》

Compliance with ethics guidelines

Wei Liang, Jiancheng Li, Xinyu Xu, Shengjun Zhang, and Yongqi Zhao declare that they have no conflict of interest or financial conflicts to disclose.

《Appendix A. Supplementary data》

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.eng.2020.05.008.