HUST-GOGRA 2018 s : A new gravity field model derived from the combination of GRACE and GOCE data

A new satellite-only gravity field model entitled HUST-GOGRA2018s is developed by the combination of GRACE and GOCE data in this study. The modified dynamic approach is applied for GRACE data processing, while the space-wise least square method with a cascade filter is utilized for GOCE data processing. The GRACE-only model HUST-Grace2016s and GOCE-only model HUST-GOCE2018s are then computed, respectively. Our new developed GRACE-only model HUSTGrace2016s performs better than AIUB-GRACE03S, GGM05S, Tongji-GRACE01S at higher degrees, and the quality of our GOCE-only model HUST-GOCE2018s is also better than that of GO_CONS_GCF_2_TIM_R2 and GO_CONS_GCF_2_SPW_ R2. The combination is subsequently implemented by the superposition of GRACE and GOCE full normal equations. During the combination, the optimal weight is determined by the least-squares combined adjustment method with parametric covariance approach (LS-PCA) and the spectral combination method, respectively. The comparison result demonstrates that LS-PCA is more proper for the combination. As a result, the final HUST-GOGRA2018s model is developed. Validated by external gravity field models, the results demonstrate that the HUST-GOGRA2018s is dominated by GRACE data for the spherical harmonic coefficients lower than degree 60 and GOCE data for the spherical harmonic coefficients higher than degree 150, and its performance is better than that of GOCO01S. Article history: Received 11 January 2018 Revised 15 October 2018 Accepted 2 November 2018


INTRODUCTION
It is significant for different disciplines, such as geophysics, geodesy and seismology, to construct global static gravity field models with high precision and spatial resolution, which are vital strategic data to understand structure of the Earth interior and solve the problems of resource and disaster (Kao et al. 2017;Tanaka et al. 2019).With the implementation of the Gravity Recovery and Climate Experiment (GRACE) on 17 March 2002 (Tapley et al. 2004) and the Gravity field and steady-state Ocean Circulation Explorer (GOCE) on 17 March 2009 (Floberghagen et al. 2011), various GRACE-only gravity field models including AIUB-GRACE03S (Jäggi et al. 2012), GGM05S (Tapley et al. 2013), ITG-Grace2010s (Mayer-Gürr et al. 2010), ITSG-Grace2014s (Mayer-Gürr et al. 2014), ITU-GRACE16 (Guo et al. 2015), and Tongji-GRACE01S (Chen et al. 2015), and independent GOCE models including GO_CONS_GCF_2_ DIR_R1, GO_CONS_GCF_2_TIM_R1 and GO_CONS_ GCF_2 _SPW_R1 (Pail et al. 2011), have been developed in recent years.In addition, many combination models purely determined by GRACE and GOCE data, GOCO01S (Pail et al. 2010), GOGRA04S (Yi et al. 2013), DGM-1S (Farahani et al. 2013), GGM05G (Bettadpur et al. 2015), and ITU_ GGC16 (Akyilmaz et al. 2016), have been also released.More details about these models can be found in the International Center for Global Gravity Field Models (ICGEM, http://icgem.gfz-potsdam.de/home).Thus, constructing Earth gravity field model is a hot topic and has been widely researched by scientists from all over the world.
The K-Band ranging (KBR) system in GRACE can Terr. Atmos. Ocean. Sci., Vol. 30, No. 1, 97-109, February 2019 be regarded as a gradiometer with an approximate 220 km baseline, which can measure the gravity field signal at large scales.In contrast, the distance between two high precision accelerometers located in the same axis of satellite gravity gradiometer (SGG) in GOCE is only 50 cm, which can detect the gravity signal at small scale.Due to the special north-south tracking pattern in GRACE mission, the sectorial and near-sectorial spherical harmonic coefficients are determined with poor quality in GRACE-only gravity field models (Wang et al. 2012;Zhou et al. 2016).Fortunately, these errors can be reduced by GOCE mission.However, due to the color noise in gravity gradients (GGs) observations, it cannot precisely model the gravity field models at low degrees with GOCE mission.Besides, the polar holes in GOCE mission also affect the modelling precision of zonal and near-zonal spherical harmonic coefficients.All these problems can be improved by GRACE observations.Hence, it is necessary to implement the combination of GRACE and GOCE for improvement of the precision and spatial resolution of gravity field models.Therefore, a new combination model is purely determined by the GRACE and GOCE data in this study.
The orbits of satellite mission provide important location information for various of payloads, and they can be also used for gravity field model determination.For instance, Jäggi et al. (2011) determined the GOCE-only gravity field model with kinematic orbits using the celestial mechanics approach.Baur et al. (2014) utilized GOCE orbits to compare the solutions derived from celestial mechanics approach, short arc approach, point-wise acceleration approach and averaged acceleration approach.Visser et al. (2014) and Jäggi et al. (2015) discussed the possibility of determining temporal gravity field model solely with kinematic orbits.However, in these works, the observations in different directions are taken as an equal weight.Thus, the influence of the weights for the observations at different directions is discussed in this study.
During GRACE data processing, the different range rate processing strategies would seriously affect the accuracy of the final solutions.For instance, using the Hill's equation, Kim (2000) demonstrated that any force components at a certain frequency would cause the perturbations of range rates at 1-cpr (circle-per-revolution). Visser (2005) also proved the necessity of introducing kinematic empirical parameters during gravity field model determination with range rates.As for these kinematic empirical parameters, Zhao et al. (2011) have noted that separate estimated the gravity field model and kinematic empirical parameters would lead to a reduction of low frequency noise.These results demonstrate the importance of processing range rates for GRACE-only gravity field mode determination.To obtain a model with higher accuracy, several efforts are conducted in this work.
The purpose of this work is to determine a new gravity field model with the combination of GRACE and GOCE data.Section 2 provides the details about the data sets as well as data processing strategies.In section 3, three new models, GRACE-only, GOCE-only, and their combined models, are developed, respectively.The performance of our new models is assessed by some representative GRACE-only models, GOCE-only models, and satellite-only GRACE/GOCE models.The conclusions are presented in section 4.

Data Sets
The Release 02 of GRACE Level 1B observations, spanning from January 2003 to April 2015, are introduced to determine GRACE-only gravity field model (Case et al. 2004).The input Level 1B data sets published by Jet Propulsion Laboratory (JPL) include KBR1B, ACC1B, and SCA1B, which denote the observations of K-band ranging system, accelerometers and star cameras, respectively.As for orbits, we make use of the pure kinematic orbits computed by Graz University of Technology, Institute of Geodesy (former ITSG).
There are many data gaps in these data sets.Considering the observations are presented as smooth curves, the cubic spline interpolation is used.The interpolation errors for various data sets with different data gaps are shown in Table 1.The interpolation error for range rate observations is larger than 2.0 × 10 -7 m s -1 , which is comparable to the nominal accuracy of range rate.To avoid such a large interpolation error, none interpolation is applied to range rate.In other words, for the epochs without range rate observations, they are not introduced for the gravity field determination.The same processing strategy is also applied to the kinematic orbits with more than 2 gaps.For the data gaps of ACC1B and SCA1B, the cubic spline interpolation is used if the missing observations are no more than 10.The GOCE Level 1b and Level 2 observations, spanning from 1 November 2009 to 11 January 2010, are utilized to determine GOCE-only gravity field model (Gruber et al. 2010;Frommknecht et al. 2011).The reduced dynamic orbits in SST_PRD_2 are exploited to fill the gaps in pure kinematic orbits SST_PKI_2.The common mode accelerations EGG_CCD_1i in EGG_NOM_1b are used to model the non-conservative forces.Considering these accelerations are in the gradiometer frame, the quaternions in EGG_ IAQ_1i are applied to rotate them to the inertial frame.The GGs are in EGG_NOM_2 in the gradiometer frame.To rotate these observations to the Earth fixed frame, the quaternions in EGG_IAQ_1i and the rotation matrix formed by IERS2010 (Petit and Luzum 2010) are applied.

SST Data Processing
Satellite-to-satellite tracking (SST) data include the kinematic orbits from GPS receivers, the ranges and their differentials from K-band ranging system.These data are reduced to the corresponding residuals by subtracting the integration counterparts.During orbit integration, the 14 orders of Gauss-Jackson integrator with median correction is applied (Berry and Healy 2004;Zhou et al. 2016).The arc step and arc length are set to 5 seconds and 24 hours, respectively.The following prior force models are introduced: (1) Ocean tide model EOT11a complete to degree and order 180 according to Savcenko and Bosch (2012).
(2) N-body perturbations including direct and indirect terms are removed based on DE421 (Folkner et al. 2013).
(3) Solid earth tides include the impacts of frequency dependent and independent terms, and the permanent tide is removed according to Petit and Luzum (2010).(4) Both solid earth pole tides and ocean pole tides are considered.Solid earth pole tides are calculated according to Petit and Luzum (2010), while ocean pole tides are determined using Desai (2002) model.( 5) Atmosphere and Oceanic variability is removed by AOD1B RL05 model complete to degree and order 100 (Flechtner and Dobslaw 2013).( 6) General Relativistic effect is estimated according to Petit and Luzum ( 2010).( 7) Non-conservative forces are determined by accelerometer observations, and they are calibrated by determining scale and bias factors per arc.
In the traditional dynamic approach (Reigber 1989), the observations at three orthogonal axes are equally weighted.The solution of gravity field model X can be summarized as where N and W are the coefficient matrix and adjoint matrix of normal equation, respectively.D is the covariance matrix and v W is the variance of unit weight.V represents the post- fit residual, and P is the weight matrix.n and s are the number of observations and unknown parameters, respectively.If we take the observations of three different directions as an equal weight, the equally weighted solutions can be solved.Due to the rotation invariance features, the equally weighted solutions computed in different frames, including Earth fixed frame, inertial frame and local north-oriented frame, have a good agreement with each other.Here, we can also use the solutions X i derived from the ith direction of observations and its corresponding unit weight variance where N i is the coefficient matrix of normal equation which is solely determined with the ith direction of observations.N c is the combined coefficient matrix of normal equation.
As seen in the equation, during the combination, the contributions of multi-direction observations are considered.As a result, using the resolution matrix R i , the final weighted solution X c is determined.As shown in Fig. 1, when the contributions of multi-direction observations are taken into consideration, the weighted solution outperforms the solutions with equal weights.Hence, the weighted dynamic approach is introduced during the gravity field model determination with GRACE and GOCE kinematic orbits.
Since the low-frequency noise in range rates seriously affects the accuracy of gravity field model determination, it is necessary to introduce empirical parameters to fit range rate residuals (e.g., Tapley et al. 2005;Visser 2005;Liu et al. 2010;Tangdamrongsub and Hwang 2016).Zhao et al. (2011) have discussed two different low-frequency noise processing strategies, and they concluded that the spherical harmonic coefficients and empirical parameters should be simultaneously determined.Based on this study, we deduced the mathematical calculation formulas of these two different low-frequency noise processing strategies in Zhou et al. (2017a), and a new strategy was created to simultaneously filter the design matrix and observation vector of observation equation.Using this new strategy, we have developed a new time series of monthly gravity field models HUST-Grace2016 (Zhou et al. 2017b), which have good agreement with CSR Release05, JPL Release05, and GFZ Release05.In this study, we also make use of this new processing strategy to remove the low-frequency noise in range rate residuals.The post-fit range rate residuals in 1 May 2005 are shown in Fig. 2, and the corresponding root mean square error (RMS) is 0.207 μm s -1 , which is close to the results in Meyer et al. (2012).More details about GRACE data processing can refer to Luo et al. (2015Luo et al. ( , 2016) ) and Zhou et al. (2017a).

SGG Data Processing
Although the GGs in EGG_NOM_2 have been calibrated and corrected by Gruber et al. (2010), they are still contaminated by colored noise and systemic error (shown in Fig. 3a).Thus, we design a cascade filter which is the combination of moving-average (MA) and autoregressive (AR) models to reduce the systemic error and colored noise, respectively.The differential equation of ARMA model can be summarized as where c(n) is a time series of color noise derived from the white noise w(n) with ARMA model.a i (i = 1, …, p) and b j (j = 1, …, q) are the coefficients of AR and MA components, respectively.To determine the coefficients, the Levinson-Durbin recursive algorithm is used here (Franke 1985;Pail et al. 2010).
The GGs are filtered by ARMA models.As shown in Fig. 3b, if the GGs are filtered only by AR(800) filter, the characteristic of the SQRT-PSD is still shown as colored noise.But after fitting the GGs with MA(60) + AR(800) cascade filter, the SQRT-PSD for the three main diagonal components of GOCE GGs V xx , V yy , and V zz remains almost the same in the measurement bandwidth, see in Fig. 3c.
Using the processed GGs, we calculated a gravity field model.The spherical harmonic coefficient differences with respect to EGM2008 are shown in Fig. 4a.Due to the 6.7° polar holes, the normal equation derived from GOCE GGs are seriously ill-conditioned, and the zonal and near zonal spherical harmonic coefficients are not well determined.The Kaula regularization is employed to mitigate the effects of polar holes.The optimal regularization parameter is selected by minimizing the geoid degree height relative to EGM2008.As it is shown in Fig. 4b, the precision of spherical harmonic coefficients improves significantly.It should be noted that the low degree portion in Fig. 4b are still not well determined, which need to be improved by GOCE orbits or GRACE observations.

HUST-Grace2016s
Using the data spanning from January 2003 to April 2015, a new GRACE-only gravity field model is determined.The new model entitled HUST-Grace2016s is complete to spherical harmonic degree and order 160.To assess the  accuracy of our new model, we compare the geoid degree heights of several representative GRACE-only gravity field models relative to EIGEN-6C4 (Förste et al. 2014).The latest EIGEN-6C4 model is developed by the combination of GOCE, GRACE, Laser Geodynamics Satellite (Lageos), gravity data and altimetry data, which perform better than a GRACE-only static gravity field model.To assess our GRACE-only model, the recent GRACE models AIUB-GRACE03S, GGM05S, Tongji-GRACE01, ITG-Grac-e2010s, and ITSG-Grace2014s are also introduced for comparison.As shown in Fig. 5, in terms of the degrees higher than 60, the ITG-Grace2010s and ITSG-Grace2014s, which are determined by short arc approach, are closest to EIGEN-6C4.In addition, since about 13 years of observations are applied, our model outperforms better than AIUB-GRACE01S, GGM05S, and Tongji-GRACE01.In terms of the degrees lower than 60, the performances of our HUT-Grace2016s model is slightly lower than that of ITSG-Grace2014s, but it is better than that of the other GRACE-only models.
The RMSs of the global gravity anomalies between GRACE-only models and EIGEN-6C4 are presented in Table 2.The RMS values for the degrees lower than 60 are smaller than 0.004 mGal for all of the GRACE-only models, which demonstrates that these models are in extremely good consistency with EIGEN-6C4 model at low degrees.For the model up to degree 140, the RMSs values are approximately 1 mGal.It indicates that the precision of the static gravity field model derived from GRACE has achieved up to approximately 1 mGal in terms of the gravity anomalies with the spatial resolution smaller than 150 km.The statistic result in Table 2 also shows that HUST-Grace2016s model performs better than AIUB-GRACE03S, GGM05S, and Tongji-GRACE01 model.The spherical harmonic coefficient differences of recent GRACE-only gravity field models are presented in Fig. 6, and our model is in good accordance with other models.The high precision of spherical harmonic coefficients in degrees lower than 80 is well determined, which can be regarded as the contribution of K-band ranging system.Due to the gravity signal attenuation with the raise of orbital altitude, the spherical harmonic coefficient differences at high degrees increase rapidly, hence the GOCE data need to be introduced.

HUST-GOCE2018s
Using the processed GOCE data mentioned in section 2.3, a new GOCE-only gravity field model entitled HUST-GOCE2018s is determined.Since the computation task grows exponentially with the increase of the truncated degree and order of solved gravity field model, our new model is only determined to degree and order 210.
One of the key processes of GOCE gravity field model determination is the selection of optimal power between GOCE kinematic orbits and GGs.Here we compared two weight determination methods, i.e., the spectral combination method and the least-squares combined adjustment method with parametric covariance approach (LS-PCA).The spectral combination method is easy to implement and the calculation speed is very fast (Zhong et al. 2012), while the LS-PCA needs iterations and it is very time consuming.Using these two weight determination methods, two GOCE-only models named HUST-GOCE2018s (Spectral) and HUST-GOCE2018s (LS-PCA) are calculated in the present study.As the geoid degree height shown in Fig. 7, these two models have compatible performance at degrees higher than 80, but HUST-GOCE2018s (LS-PCA) performs better than HUST-GOCE2018s (Spectral) at lower degrees.Therefore, LS-PCA is selected as the optimal approach to combine GOCE kinematic orbits and GGs here, and the corresponding model is named as HUST-GOCE2018s.
Figure 7 shows the comparison between our new model and official ESA models in terms of geoid degree height relative to EGM2008.GO_CONS_GCF_2_DIR_R2, GO_ CONS_GCF_2_SPW_R2, and GO_CONS_GCF_2_TIM_ R2, which are respective complete to degree and order 240, 240, and 250, are the second release of official ESA models.They are computed by direct approach, time-wise approach and space-wise approach, respectively (Pail et al. 2011).Because of the increasing noise at low-frequency of GGs, the lower degrees of GOCE-only models are dominated by errors in kinematic orbits.The performance of our model for the degrees lower than 90 is better than GO_CONS_GCF_2_ SPW_R2 and GO_CONS_GCF_2_TIM_R2, which can   be regarded as the improvement of our modified dynamic method compared with the energy integral method applied in time-wise approach and space-wise approach.This feature is also clear in Fig. 8, in which the errors in zonal and near zonal spherical harmonic coefficients of HUST-GOC-E2018s are smaller than those determined by time-wise approach and space-wise approach.In addition, the lower degrees of GO_CONS_GCF_2_DIR_R2 are dominated by the prior information of reference model, which results in higher accuracy than other GOCE-only models.The errors in terms of geoid degree height are close to those degrees between 90 and 180 for different GOCE models, while the spherical harmonic coefficient differences of HUST-GOC-E2018s shown in Fig. 8 are better than those of other three ESA official models, which are probably the result of different regularization methods.In the frequency band higher than degree 180, our models also match the EGM2008 better, which may be due to the effect of different regularization methods and different truncated degrees.

HUST-GOGRA2018s
HUST-GOGRA2018s complete to degree and order 210 is a satellite-only gravity field model derived from the combination of GRACE and GOCE data.The input data are the same as those for the computation of HUST-Grace2016s and HUST-GOCE2018s, and the combination is implemented by adding the corresponding full normal equations built by GRACE and GOCE data.The Kaula regularization is introduced as a constraint for zonal and near zonal spherical harmonic coefficients.The optimal weight factor of each individual data set is estimated by parametric covariance approach.
Figure 9 shows the percentage contribution of GRACE and GOCE solution to the combined model HUST-GO-GRA2018s.The percentage contribution of GRACE solution for the lower degree portion (up to degree 60) is close to 1.0, indicating the dominate contribution of inter-satellite range rates to the combined model.It is a mixture domain between degree 60 and degree 150, in which the contribution of GRACE solutions decreases while that of GOCE solutions increases.Comparing the spherical harmonic coefficient differences of HUST-GOGRA2018s in Fig. 10a and HUST-GOCE2018s in Fig. 8a, we can find that the errors in the zonal and near zonal spherical harmonic coefficients of the former solution are smaller than those of latter solution.It can be better understood by Fig. 9a, which demonstrates that the zonal and near zonal harmonic coefficients lower than degree 120 is exclusively based on GRACE-only solution.Large values can be found for the degrees higher than 150 in GOCE solution, indicating the significant contribution of GOCE data at this part.
Figures 10 and 11 presents the comparison result in terms of spherical harmonic coefficient differences and geoid degree height, respectively.Due to the high sensitivity of the long wavelength part of gravity spectrum, the low degree errors of GRACE solutions (HUST-Grace2016s and ITG-Grace2010s) and combined solutions (HUST-GO-GRA2018s and GOCO01s), in terms of both geoid degree height and cumulative geoid error, are smaller than those of GOCE solutions (HUST-GOCE2018s and GO_CONS_ GCF_2_DIR_R2).In addition, the contribution of GOCE solution for the short wavelength part of gravity spectrum can be clearly seen at the high degrees, in which the errors of GOCE solutions and combined solutions are smaller than GRACE solutions.For example, the geoid degree heights of HUST-GOCE2018s and HUST-GOGRA2018s are gradually smaller than HUST-Grace2016s and ITG-Grace2010s after degree 120.Moreover, the spherical harmonic coefficient differences of our combined solution are slightly smaller than those of GOCO01s, which was derived from almost the same time span of observations as HUST-GO-GRA2018s.It demonstrates the good performance of our new combine model.
To compare the models derived from different observations, the cumulative geoid errors for different models are also presented in Table 3.The key points from Table 3 can be summarized as: (1) Since GRACE data is sensitive to the long wavelength part of gravity spectrum, the GRACEonly models (HUST-Grace2016s and ITG-Grace2010s) and combined models (HUST-GOGRA2018s and GOCO01s) outperforms better than the GOCE-only models (HUST-GOCE2018s and GO_CONS_GCF_2_DIR_R2) at the degrees lower than 60; (2) Our HUST-GOGRA2018s model performs the best between degree 90 and degree 160, which is the mixture part for the contribution of both GRACE and GOCE data; (3) For the spherical harmonic coefficients higher than degree 160, our HUST-GOGRA2018s model also performs the best, which is mainly derived from the contribution of GOCE data.In terms of the cumulative geoid heights, our HUST-GOGRA2018s model performs better than GOCO01s.

CONCLUSIONS
In this paper, a new GRACE-only gravity field model named HUST-Grace2016s is developed by a modified dynamic approach.In terms of data processing, the interpolation errors for different types of observations are analyzed.For the orbit and range rate observations, none of interpolation is recommended.The contribution of the observations in different directions is also considered.The result demonstrates that the weighted solution performs better than the equally weighted solution.As for the GGs processing, the cascade filter is applied.The PSD of the filtered GGs demonstrates the good performance of this processing strategy.In addition, before calculating the combined solution, the different approaches for determining the optimal weight    are also compared.The comparison result indicates that LS-PCA performs better than the spectral combination method.
Base on the improved data processing strategies, the new gravity field models are developed.Compared with other representative GRACE-only model, the performance of our model is better than AIUB-GRACE03S, GGM05s, and Tongji-GRACE01 in terms of geoid degree height and spherical harmonic coefficient differences.The GRACEonly model has achieved up to approximately 1 mGal in terms of the gravity anomalies with the spatial resolution smaller than 150 km.The new GOCE-only gravity field model named HUST-GOCE2018s also performs better than GO_CONS_GCF_2_SPW_R2 and GO_CONS_GCF_2_ TIM_R2.In addition, for the spherical harmonic coefficients higher than degree 180, the performance of HUST-GOC-E2018s is also better than GO_CONS_GCF_2_DIR_R2. Finally, a new combined gravity field model HUST-GO-GRA2018s is computed from the GRACE and GOCE data by LS-PCA.In terms of spherical harmonic coefficient differences, geoid degree height and cumulative geoid error, HUST-GOGRA2018s performs better than GOCO01s.

Acknowledgements
The authors are very grateful to two anonymous reviewers for their valuable comments, which helped to improve the quality of the manuscript significantly.The GRACE and GOCE data in this paper are provided by JPL and ESA, respectively.This work is mainly supported by the National Natural Science foundation of China (Grant No. 41704012, No. 41504015, No. 41474019).It is also partly sponsored by the Project funded by China Postdoctoral Science Foundation (Grant No. 2016M592337) and the Open Research Fund Program of the State Key Laboratory of Geodesy and Earth's Dynamics (Grant No. SKLGED2017-1-2-E, No. SKLGED2018-1-3-E).

Fig. 1 .
Fig. 1.Geoid degree height of solutions derived from GOCE kinematic orbits in December 2009, the equal weight (red line) and optimal weight (blue line) of the observations at three orthogonal axes are added respectively.

Fig. 3 .Fig. 4 .
Fig. 3. Square-root of error power spectral density (SQRT-PSD) for the three main diagonal components of GOCE GGs V xx , V yy , and V zz : (a) before filtering with 1 second sampling; (b) after filtering only by AR(800) filter; (c) after filtering by MA(60) + AR(800) cascade filter.The dot line is the boundary of measurement bandwidth.

Table 2 .
RMSs of the global anomaly differences between GRACE-only static gravity field models and EIGEN-6C4 (mGal).

Table 3 .
Cumulative geoid errors of different models relative to EGM2008 (cm).