An alternative method to improve gravity field models by incorporating GOCE gradient data

The aim of this paper is to present an alternative method that can be used to improve existing gravity field models via the application of gradient data from Gravity field and Ocean Circulation Explorer (GOCE). First, the proposed algorithm used to construct the observation equation is presented. Then methods for noise processing in both time and space domains aimed at reducing noises are introduced. As an example, the European Improved Gravity model of the Earth by New techniques (EIGEN5C) is modified with gradient observations over the whole lifetime of the GOCE, leading to a new gravity field model named as EGMGOCE (Earth Gravitational Model of GOCE). The results show that the cumulative geoid difference between EGMGOCE and EGM08 is reduced by 4 centimeters compared with that between EIGEN5C and Earth Gravitational Model 2008 (EGM08) up to 200 degrees. The large geoid differences between EGMGOCE and EIGEN5C mainly exist in Africa, South America, Antarctica and Himalaya, which indicates the contribution from GOCE. Compared to the newest GOCE gravity field model resolved by direct method from European Space Agency (ESA), the cumulative geoid difference is reduced by 7 centimeters up to 200 degrees. ABSTRACT


Introduction
The Gravity field and steady-state Ocean Circulation Explorer (GOCE) satellite (European Space Agency, 1999) was launched on March 17 th , 2009 and was terminated on November 11 th , 2013. By combining the gradient and orbit data from GOCE, three groups of gravity field models with increasingly high accuracies have been derived by the ESA (European Space Agency) GOCE level 2 processing facility. Indeed, in addition to GOCE data, a number of other kinds of observations can also be used to derive gravity field models, including observations from altimeter satellites (Marchenko, 2003a;Flechtner et al., 2006), the CHAllenging Minisatellite Payload (CHAMP) (Reigber et al., 1996) and the Gravity Recovery and Climate Experiment (GRACE) (Tapley et al., 2004), as well as some ground observations (Marchenko, 2003a;Pavlis et al., 2012). The use of such data has enabled a range of institutes to resolve more than 100 gravity field models, released by the International Centre for Global Earth Models operated by the German Research Centre for Geosciences. In comparison to other kinds of gravity information, GOCE data are of higher accuracy in some special frequency bands, and the spatial consistency of the accuracy is superior to that of ground observations (Marchenko et al., 2016). Nevertheless, deriving an accurate gravity field model using GOCE gravity gradient observations in combination with other data remains an important issue (Gilardoni et al., 2013).
One possible approach to this question is to initially construct a series of normal equations based on all the observations (Pail et al., 2010) and then resolve them using least square estimation method (LSM). However, as the large quantity of data makes the computational cost of this process prohibitive, an alternative method is to utilize GOCE gravity gradient data to improve existing gravity field models that have been derived from other kinds of observations. If this approach is adopted, there is no need to process the other (non-GOCE) observation data over again.
In this study, we adopt the latter method in which the existing gravity field models are selected as reference models to be improved. To resolve the gravity field model, two approaches can be used, i.e., LSM and spherical harmonic analysis. Currently, the former is more commonly applied for processing GOCE gradient data for gravity field recovery, for example, the resolving of time-wise solution (Pail et al., 2011) and direct solution of the GOCE gravity field models (Bruinsma et al., 2010) by ESA, DGM-1S (Hashemi et al., 2013), etc. One of the main advantages of the LSM approach is the ease with which a gravity field model can be resolved using gradient and orbit data in combination. Considering that the existing gravity field models are improved mainly by gradient data in this study, we adopt spherical harmonic analysis to recover the gravity field model.
Resolving of a gravity field model via the use of spherical harmonics can be treated as a Boundary Value Problem (BVP) (Rummel, et al., 1993), the first step of which is to construct a so-called boundary value condition based on the observation equation. A great deal of research has been conducted in this area. Rummel and van Gelderen (1992) analyzed the spectral characteristics of different gravity gradient components and generated a series of two-dimensional Fourier expressions. van Gelderen and Rummel (2001) summarized the solution of a general geodetic BVP under gravity gradient boundary conditions. Martinec (2003) discussed Green function to gradiometric BVP. In addition, invariants of gradient tensors (Rummel, 1986, Sacerdote & Sanso, 1989, Vermeer, 1990, Marchenko, 2003b, Baur et al., 2008 were discussed to recover gravity field models by LSM. Holota (1989) and Yu and Zhao (2010) constructed the boundary conditions using invariants of a gravity gradient tensor. Those conditions contain few attitude errors but create linear errors. Marchenko et al.(2016) investigated a special version of the space-wise method based on the second degree radial derivatives of GOCE missions. Their results show that the second order derivatives of EGG_TRF_2 data can improve the middle and long waves of the geo-potential (Marchenko et al., 2016).
The radial gravity gradients are utilized in this study to construct the boundary condition. It is noteworthy that Fuchs et al. (2013) proposed a superior and highly accurate method for vertical gravity gradient computation, which they successfully applied to a discussion of the Tokyo earthquake. By this method, Bouman et al.(2016) calculated satellite gravity gradient grids for geophysics. Meng et al.(2016) proposed an alternative combination of V xx , V yy and V zz to compute a new V zz which then is used for computing the radial gravity gradient. Since accurate radial gravity gradients can be obtained, we believe that they can be used as a boundary condition to improve some existing gravity field model and to obtain new ones. This is the starting point of this paper.

Radial gravity gradient computation
GOCE measurements include the V xx , V yy , and V zz components of the gravity gradient tensors (Stummer et al., 2012). Thus, if no noise is observed, the Laplace equation can be applied, as follows: This leads us to obtain two vertical gravity gradient values ( i.e., V zz and V zz ' ) that are, however, unequal to one another because of noise. Fuchs et al. (2013) addressed this by using mean values of V zz and V zz ' to improve the accuracy of vertical gravity gradients. Meng et al. (2016) proposed the following combination to compute new vertical gravity gradients, Compared to V zz , noises of V zz c , are reduced about 42.3% (Meng et al., 2016). However, V zz c , are the vertical components of gravity gradient tensors in Gradiometer Reference Frame (GRF) and not the radial gravity gradients, hence the coordinate transformation is needed. Although Fuchs et al. (2013) and Meng et al. (2016) proposed new methods for computing new vertical gravity gradients, horizontal gravity gradiens (ie., V xx and V yy ) are also needed for the data rotation. However, no literature discussed how to compute high accurate V xx and V yy . Indeed, according to the Laplace equation, we can construct the following equations, This leads us to obtain two x-components (i.e., V xx and V xx ' ) and two y-components ( i.e., V yy and V yy ' ) of a gravity gradient tensor. We provide combinations for these values by solving the following equations: Thus, as long as  ii jj i x y z , , , 2 = ( ) are known, it is possible to obtain optimal factors from equation (6) and (7). Indeed, if the components of the gradient tensor are independent, equation (6) and (7) The error distribution of the GOCE gradient tensor applied by Fuchs et al. (2013) is used here, as follows: Solutions for a and b , c and d can therefore be derived from equation (8) and equation (9), as follows: » 0 91 0 . The noises of horizontal gravity gradients are reduced by about 8.9%. In this paper, besides V xx c , and V yy c , expressed by equation (13), V zz c , expressed by equation (2) are used to compute the radial gravity gradients by the following coordinate transformation, where R denotes the transformation matrix from the GRF coordinate system to the Local North Oriented Frame (LNOF)  coordinate system. In LNOF the Z direction is same as the radial direction of the Earth. It should be noted that, owing to the poor accuracy of V xy and V yz , we replace them with U xy and U yz respectively, that were calculated using the reference model (i.e., the model to be modified). Since accuracies V xx c , , V yy c , , and V zz c , are higher than those of V xx , V yy , and V zz , this means that the accuracy of V rr can be improved via the use of V xx c , , V yy c , , and V zz c , .
Residual gradient data were used in this study for data rotation as follows: EGGT R   T  T  T   T  T  T   T  T  T   LONF   T  which are obtained from the reference model. T xy and T xz are both equal to zero because V xy and V yz are both replaced by U xy and U yz , respectively. The use of this rotation approach enables us to obtain T rr .
The variable T rr , obtained from equation (15), was then used to construct the following boundary value problem ( Wan & Yu, 2013), as follows: where f denotes T rr , and S refers to the mean surface of GOCE orbits. Thus, the solution of the boundary value problem (16) where GM denotes the gravitational constant of the Earth, a is the semimajor axis of the Earth, r, ,

Noise processing of GOCE gradient data
Before using GOCE gradient observations to recover a gravity field model, the observed noise should be processed first (Schuh, 2003;Reguzzoni & Tselfes, 2009;Bruinsma et al., 2010;Pail et al., 2011;Gatti et al., 2014). Three types of noise should be considered. The first one is low frequency noise. As implied by the name, this type of noise is located in the low frequency band (i.e., lower than 0.005 Hz in this study) and constitutes the main type of noise in GOCE gradient data. This low frequency noise is large in magnitude, and within a known frequency band. The second type of noise is random noise. This noise is smaller in magnitude and is distributed throughout the full frequency band. The combinations presented in equations (2) and (13) are used to reduce this random noise. Hence, the techniques of white noise processing are not discussed here. The third one is gross error. Their magnitude is very large and its distribution is scattered. Gross error is found in all forms of observational data used in the gravity field recovery, including gradient, orbit, and attitude data. Gross error can also be created by a kind of computing, e.g., the warm effect of filtering.

a. Time domain processing
Forward and backward Finite Impulse Response (FIR) filtering was initially applied to the residual gravity gradient data, which are the difference values of the actual GOCE gradient data and the simulated gradient data from the reference gravity field model (i.e., the reference model to be modified). As noted, EIGEN5C model was selected here as the reference model, which was derived by GFZ and CNES ( National Centre for Space Studies) using five years of GRACE data as well as some ground observations. We only selected 0-300 degrees / orders of EIGEN5C as the reference model, considering the limits of measurement band wise (MBW) and the accuracy of the gradient observations.
To show the filtering effect, we processed the gradient data provided by EGG_NOM2 from ESA over the short time interval between 0:0:1.6 on October 17th , 2010 and 23:59:59.6 on October 26th , 2010, using a band-pass FIR filter. The pass band frequency is 0.005~0.1 Hz. Fig. 1 shows the Power Spectral Density (PSD) of the data after filtering. It is clearly that signals in the pass band remained unchanged by this filtering processing, while those beyond the pass band were largely removed and only occupy a very small proportion in the final data obtained from the filtering. This filtering step therefore enabled us to remove low frequency noise while saving the signals during the pass band. It is also noteworthy that no phase drift occurred (e.g., Yu & Wan, 2011); thus, the data after filtering can still be considered as gravity gradient signals.
In order to verify the effects of FIR processing, we perform a statistical analysis to the absolute values of Laplace Equation. Table 1 gives the results. Smaller values of this statistic indicate a higher accuracy; thus, according to Table 1, the accuracy was improved by a factor of 159172 by FIR filtering in terms of the mean value.

b. Space domain processing
Although the dominant noises have been removed in the time domain via FIR filtering, we still need to remove the residual gross errors. These errors were processed in the space domain. When the radial gradients of the disturbing gravity were obtained using equation (15), we were able to perform grid processing on the satellite orbit surface. We first divided the surface into girds, such as 20' × 20' and then conducted statistical processing. Mean values, variances, and standard errors of the data in each grid were then obtained, and we removed the data that deviated from the mean value by more than twice the standard deviation. This processing led to the removal of approximately 3.84% of the data.
Processing was carried out in the space domain because gross errors are sometimes stable over short time periods, perhaps because of an abnormal satellite situation, making these errors difficult to remove in the time domain. However, even though the data in the same grid are from different points in time, their values should be the same or close to each other. If some abnormal values exist, there must be gross errors that can be isolated.

Results and Discussion
All the gradient data obtained by GOCE throughout the lifetime of the project are processed in this paper. As described above, we first processed the noises in both time and space domains. Radial disturbing gravity gradients were then calculated using equation (15). In this section, the EIGEN5C is taken as an example to be improved with the proposed method. The final solution is derived by integration based on the orthogonality of spherical harmonic functions given by equation (16) and equation (17). The new model is named as the Earth Gravitational Model of Gravity field and Ocean Circulation Explorer (EGMGOCE).

Comparisons in spherical domain
The degree variances of EGMGOCE and EGM08 (Pavlis et al., 2012), are compared in Fig. 2 and the geoid heights are compared in Fig. 3.
In Fig. 2, Variance_EIGEN5C denotes the degree variance differences between EGM08 and EIGEN5C, and Variance_EGMGOCE denotes the differences between EGM08 and EGMGOCE. In Fig. 3, EIGEN5C-EGM08 and EGMGOCE-EGM08 indicate the cumulative geoid differences between Figure 1. Comparison of PSD before and after FIR filtering EIGEN5C, EGMGOCE and EGM08 respectively. These two figures illustrate that the initial model was mainly improved in the degrees from 100 to 200. If EGM08 is considered as the true model, the geoid accuracy of the new model is higher than that of the initial model by 4 cm up to 200 degree.
For a comparison with GOCE models, EGMD5 from ESA is selected, which is the newest GOCE gravity field model derived by the direct method. The degree variance differences and cumulative geoid differences between EGMGOCE, EIGEN5C and EGMD5 are shown in Fig. 4 and Fig. 5, respectively. The new model is clearly much closer to EGMD5 than EIGEN5C, particularly at degrees from 80 to 230. If EGMD5 is considered as the true model, the geoid accuracy of the new model is higher than that of the initial model by 7 cm up to 200 degree. These results indicated that the extraction of the GOCE gravity gradient signals and modification of EIGEN5C were successful.

Comparisons in space domain
The differences between EGMGOCE and EIGEN5C can be seen as the contribution from GOCE gradient observation. In order to present GOCE contribution in space domain, we selected 0-200 degrees / orders of EIGEN5C, EGMGOCE, EGM08 and EGMD5 to compute geoid values at mean surface of the Earth. The data were calculated at grids of 1˚×1˚. Statistics of their differences are given in Table 2 and Table 3. According to Table 2, if EGM08 is the true model, the accuracy of EIGEN5C is improved by 4.6 cm because of the use of GOCE gradient data. According to Table 3, if EGMD5 is the true model, the accuracy of EIGEN5C is improved by 7 cm. In order to show the differences of the geoids between EIGEN5C and EGMGOCE, Fig. 6 and Fig. 7 show the geoids calculated by EIGEN5C and EGM08 respectively, and Fig. 8 shows their differences. Table 4 gives the statistics of absolute values of the differences. According to this table, mean value of the difference is about 5 cm and the standard deviation is about 9 cm. According to Fig. 8, large differences mainly occur in Africa, South America, Antarctica and Himalaya (Yi, 2001;Rummel & Gruber, 2012). This indicates that the data used to derive EIGEN5C has a little poorer accuracy in these areas, since accuracy of GOCE gradient data has a high space consistence, so that accuracy of EIGEN5C can be improved by GOCE gradient observations.

Conclusion
This paper discusses the method for modifying existing gravity field models using GOCE gradient observations. We used combinations of V xx , V yy and V zz to compute new horizontal and vertical gravity gradients, which can reduce white noises of them by 8.9%, 8.9% and 42.3% respectively.
Compared with the conventional method, accuracies of the radial gravity gradients can be improved correspondingly.
Strategies for processing noise within GOCE gradient observation data in both time and space domains are also discussed in this study. Band pass FIR filtering is adopted in the time domain and grid processing based on the statistical properties of each grid is implemented in the space domain. These strategies enable the removal of low frequency noises and gross errors. The combinations given by equations (2) and (13) can reduce white noises.
The EIGEN5C model was improved as an example in this study using all available GOCE gradient data, and a new gravity field model named as EGMGOCE was derived. Compared with EGM08, modifications using GOCE gradient data are mainly reflected in the degrees from 100 to 200. Up to the 200th degree, the cumulative geoid difference is reduced by 4cm. Compared with EGMD5, modifications from GOCE data are mainly reflected in degrees from 80 to 230, and up to the 200th degree, the cumulative geoid difference is reduced by 7 cm. These results showed that the GOCE gradient observations can contribute meaningfully to the development of EGMGOCE model and demonstrate the effectiveness of the methods adopted in this paper. Other gravity field models can also be processed similarly to derive new models. Figure 2. Differences in the degree of variance between EIGEN5C, EGMGOCE and EGM08

Figure 3
Cumulative geoid height differences. Red line: differences between EIGEN5C and EGM08; Green line: differences between EGMGOCE and EGM08.    Since the radial gravity gradients as well as the horizontal components of gravity gradients can be derived with high accuracy, they could not only be utilized to recover gravity field models, but also be used to extract and interpret geophysical signals directly, such as Moho geometry (Ye et al., 2016 ), earthquake (Fuchs, et al., 2013), etc. Hence, the methods proposed in this paper can also be powerful tools for researches in geophysics.