Research

Wei Ling, Jinheng Li,b, Xinyu Xu,b,*, Shengjun Zhng, Yongqi Zho

a School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China

b Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan University, Wuhan 430079, China

c School of Resources and Civil Engineering, Northeastern University, Shenyang 110004, China

Keywords:Gravity field model GOCE GRACE Satellite altimetry Block-diagonal least-squares

ABSTRACT This paper focuses on estimating a new high-resolution Earth’s gravity field model named SGG-UGM-2 from satellite gravimetry, satellite altimetry, and Earth Gravitational Model 2008 (EGM2008)-derived gravity data based on the theory of the ellipsoidal harmonic analysis and coefficient transformation(EHA-CT). We first derive the related formulas of the EHA-CT method, which is used for computing the spherical harmonic coefficients from grid area-mean and point gravity anomalies on the ellipsoid. The derived formulas are successfully evaluated based on numerical experiments.Then,based on the derived least-squares formulas of the EHA-CT method,we develop the new model SGG-UGM-2 up to degree 2190 and order 2159 by combining the observations of the Gravity Field and Steady-State Ocean Circulation Explorer (GOCE),the normal equation of the Gravity Recovery and Climate Experiment(GRACE), marine gravity data derived from satellite altimetry data, and EGM2008-derived continental gravity data. The coefficients of degrees 251-2159 are estimated by solving the block-diagonal form normal equations of surface gravity anomalies (including the marine gravity data). The coefficients of degrees 2-250 are determined by combining the normal equations of satellite observations and surface gravity anomalies.The variance component estimation technique is used to estimate the relative weights of different observations. Finally, global positioning system (GPS)/leveling data in the mainland of China and the United States are used to validate SGG-UGM-2 together with other models, such as European improved gravity model of the earth by new techniques (EIGEN)-6C4, GECO, EGM2008, and SGG-UGM-1 (the predecessor of SGG-UGM-2). Compared to other models, the model SGG-UGM-2 shows a promising performance in the GPS/leveling validation. All GOCE-related models have similar performances both in the mainland of China and the United States, and better performances than that of EGM2008 in the mainland of China. Due to the contribution of GRACE data and the new marine gravity anomalies, SGG-UGM-2 is slightly better than SGG-UGM-1 both in the mainland of China and the United States.

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 Steady-State 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.

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 har-monic 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.

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

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.1. Ellipsoidal harmonic analysis and coefficient transformation method

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 rΔg 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.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:

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 ofcoefficients,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 qin-2kare 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-17throughout 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. 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.

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-9around degree 5. The degree error RMS values greater than 10-11are 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 × 10-5m·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).

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

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 kindsof 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.

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

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

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-41 mHz 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 Δt = 5 s [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

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].

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.

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

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.

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 theoriginal 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.

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

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. 3. The schematic diagram of combining the NEQs of GRACE satellite, GOCE satellite, and surface gravity anomalies.

4.4. Computation of the SGG-UGM-2 model

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

Fig. 5. The scheme of computing 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

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.

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 EIGEN-6C4 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. 6. The degree RMS of the coefficient differences between the three models(SGG-UGM-2, EIGEN-6C4, and GECO) and EGM2008.

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.

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 westeast 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.

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

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

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

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.

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].

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.

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.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-UGM-2, and (e) GECO.

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.

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 Statistical results of comparison with GPS/leveling data in Qingdao (152 points) (unit: m).

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

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.

6. Conclusions

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. 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.

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

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

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

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