The earth’s gravity field recovery using the third invariant of the gravity gradient tensor from GOCE

Due to the independence of the gradiometer instrument’s orientation in space, the second invariant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_2$$\end{document}I2 of gravity gradients in combination with individual gravity gradients are demonstrated to be valid for gravity field determination. In this contribution, we develop a novel gravity field model named I3GG, which is built mainly based on three novel elements: (1) proposing to utilize the third invariant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_3$$\end{document}I3 of the gravity field and steady-state ocean circulation explorer (GOCE) gravity gradient tensor, instead of using the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_2$$\end{document}I2, similar to the previous studies; (2) applying an alternative two-dimensional fast fourier transform (2D FFT) method; (3) showing the advantages of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_3$$\end{document}I3 over \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_2$$\end{document}I2 in the effect of measurement noise from the theoretical and practical computations. For the purpose of implementing the linearization of the third invariant, this study employs the theory of boundary value problems with sphere approximation at an accuracy level of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$O(J_2^2\cdot T_{ij})$$\end{document}O(J22·Tij). In order to efficiently solve the boundary value problems, we proposed an alternative method of 2D FFT, which uses the coherent sampling theory to obtain the relationship between the 2D FFT and the third invariant measurements and uses the pseudo-inverse via QR factorization to transform the 2D Fourier coefficients to spherical harmonic ones. Based on the GOCE gravity gradient data of the nominal mission phase, a novel global gravity field model (I3GG) is derived up to maximum degree/order 240, corresponding to a spatial resolution of 83 km at the equator. Moreover, in order to investigate the differences of gravity field determination between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_3$$\end{document}I3 with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I_2$$\end{document}I2, we applied the same processing strategy on the second invariant measurements of the GOCE mission and we obtained another gravity field model (I2GG) with a maximum degree of 220, which is 20 degrees lower than that of I3GG. The root-mean-square (RMS) values of geoid differences indicates that the effects of measurement noise of I3GG is about 20% lower than that on I2GG when compared to the gravity field model EGM2008 (Earth Gravitational Model 2008) or EIGEN-5C (EIGEN: European Improved Gravity model of the Earth by New techniques). Then the accuracy of I3GG is evaluated independently by comparison the RMS differences between Global Navigation Satellite System (GNSS)/leveling data and the model-derived geoid heights. Meanwhile, the re-calibrated GOCE data released in 2018 is also dealt with and the corresponding result also shows the similar characteristics.


Results
The gravity field determination from the third invariant depends on two kinds of data, i.e. the GOCE gravity gradient measurements and the synthetic gravity gradients derived from a priori gravity field. The GOCE gravity gradient measurements of the components V xx , V yy , V zz and V xz in the GRF are taken from the Level-2 product EGG_NOM_2 (ESA) from November 1, 2009 to August 1, 2012. The synthetic gravity gradients of the components V xy and V yz are derived from the global gravity field model EIGEN-5C. Then a new gravity field model (I3GG) is obtained by applying the theory of boundary value problems for the third invariant with sphere approximation and the modified 2D FFT method proposed based on the coherent sampling theory and the pseudo-inverse via QR factorization. Details are given in the "Methods" section.
In order to evaluate its performance, the EGM2008 is taken as the reference and then we compare the error degree amplitudes between I3GG and other models, including the widely-used GOCE models (i.e. GO_ CONS_GCF_2_DIR_R2 (DIR_R2), GO_CONS_GCF_2_TIM_R5 (TIM_R5) and GO_CONS_GCF_2_SPW_R5 (SPW_R5)), EIGEN-5C and IGGT_R1, as shown in Fig. 1. It should be pointed out that models DIR_R2, TIM_ R5, SPW_R5 and IGGT_R1 are developed only from GOCE measurements and the regularization is applied to the coefficients of TIM_R5, SPW_R5 and IGGT_R1 at degrees above 200 while DIR is recovered without regularization 16,19 , which causes a larger maximum degree of recovery in the spectrum. From Fig. 1, it is seen that the maximum recovery degree of I3GG is 240, corresponding to a spatial resolution of 83 km at the equator. Since the signals of the recovered gravity model below the bandwidth of the band-pass filter mainly come from the reference model, results show that the I3GG and EIGEN-5C are relatively consistent at degrees below 27 ( ≈ 5 mHz). The degree amplitudes of GOCE-derived models become higher than EIGEN-5C in the transition band (degrees 28-69). For degrees , all models except EIGEN-5C are very close to each other because that the signals in this bandwidth are mainly provided by GOCE gravity gradients measurements. The bumps of these curves are related to the improvements from GOCE mission to the inclusion of low-accuracy terrestrial data in certain regions in the EGM2008 model 4 . The amplitudes of the EIGEN-5C model are slightly higher in this bandwidth due to the contribution of GRACE data rather than GOCE data. Briefly, thanks to the advantages of the gradiometer measurements, the accuracy of the GOCE models adopted in this paper outperforms the EGM2008 and EIGEN-5C models in this spectral range. Between the models without regularization, the model www.nature.com/scientificreports/ I3GG is lower than the DIR_R2 for higher degrees (> 200). But in the same range of degrees, both of them are higher than the models which benefit from the regularization, i.e. the TIM_R5, SPW_R5 and IGGT_R1, or from the terrestrial gravity field data i.e. the EIGEN-5C. Figure 2 shows the I3GG coefficients deviation to EGM2008 and EIGEN-5C in logarithmic scale. The spherical harmonic coefficients of I3GG and EGM2008 shown in Fig. 2a fits well to each other for degree/order 0 to 70 and the coefficients deviation becomes larger for higher degrees/orders, which means that GOCE gravity field models contribute new information in high-frequency parts. The similar situation happens between the coefficients of I3GG and EIGEN-5C except for the fitted range that is up to degrees/orders 100, as shown in Fig. 2b. It can be also concluded that spherical harmonic coefficients of EIGEN-5C are more accurate than that of EGM2008 at degrees between 70 and 100, if we use the I3GG as a standard.
In addition to the spectral comparisons displayed in Figs. 1 and 2, the differences between these models are investigated in the spatial domain. First, we computed the cumulative geoid height deviation between the I3GG and EGM2008 (Fig. 3a) and EIGEN-5C (Fig. 3b) up to degree/order 210. From Fig. 3, it is obvious that at regions where no or poor quality terrestrial data are available (Himalaya, Africa, Amazonas and Antarctica) new gravity field information is added, which is due to the advantages of the gravity gradient measurements of the GOCE mission in the medium and short wavelength ranges.
Then we compute the RMS of the differences between these models from 1 • × 1 • geoid height grids between −80 • and 80 • latitude (i.e. without the polar caps). The RMS values are computed for the common wholefrequency range (degrees 0-210) and for the medium-frequency range (degree 70-150), where these models have a similar behaviour, as given in Table 1. Moreover, in order to inspect the differences between the invariant I 2 and I 3 in the gravity field determination, the same processing strategy is applied on the second invariant I 2 measurements of GOCE mission and we obtained another gravity field model named I2GG. Its RMS values of the geoid differences between other models are also given in Table 1 and accordingly its difference degree amplitudes from the I3GG are shown in Fig. 4.
For the common whole-frequency range (0-210), Table 1 shows that the I3GG solution is closer to the invariant solutions I2GG and the IGGT_R1 (RMS = 0.076-0.091 m) than to the other GOCE models DIR_R2, TIM_R5, SPW_R5 (RMS = 0.119-0.137 m). In this case, the I3GG and I2GG solutions show a stronger consistency with the EIGEN-5C than the other GOCE models. When replacing the standard model by the EGM2008, the I3GG and IGGT_R1 solutions have lower RMS values of geoid differences than the ones except the regularized, i.e. the TIM_R5 and SPW_R5. For the medium-frequency range (70-150) it is visible that the difference among GOCE modes is much smaller (RMS = 0.014-0.078 m) than the EGM2008 (RMS = 0.081-0.120 m). Since I3GG, I2GG and IGGT_R1 have taken the EIGEN-5C as a priori model, they are closer to the EIGEN-5C than the DIR_R2, TIM_R5 and SPW_R5, which holds true both for the medium-frequency and the common whole-frequency range.
Meanwhile, the re-calibrated GOCE data released in 2018 is also dealt with and the corresponding gravity field model named I3GG_R2 and I2GG_R2 from the invariant I 2 and I 3 , respectively, as shown in Fig. 1. It indicates that the I3GG_R2 and I2GG_R2 also agree well with other models for degree 70 to 150, which is the most sensitive bandwidth for the signals detected by GOCE mission. For lower and higher degrees they are slightly lower than I3GG and I2GG, respectively. The same phenomenon appears in the comparison between TIM_R5 and TIM_R6 (GO_CONS_GCF_2_TIM_R6). It should be pointed out that the Kaula-regularization applied to coefficients of degrees/orders 201-300 of both TIM_R5 and TIM_R6 so that they are almost the same in the range. Table 2 shows the RMS valus of geoid differences of TIM_R5 versus TIM_R6, I2GG versus I2GG_R2 and I3GG versus I3GG_R2 are 0.011, 0.008 and 0.039 for the medium-frequency range,and 0.038, 0.025 and 0.082 for the whole-frequency range, respectively. It indicates the re-calibration leads more improvements to the lower  www.nature.com/scientificreports/ and higher degrees than the medium-frequency range. Additionally, if TIM_R6 is used as a standard, Table 2 shows that for the lower and higher degrees the rms values of geoid differences of the models from the invariant I 3 , i.e. the I3GG and I3GG_R2, are lower than the ones from I 2 , i.e. I2GG and I2GG_R2, respectively. Certainly, because that I2GG_R2 and I3GG_R2 are also derived from the re-calibrated GOCE data, their the rms values of geoid differences with TIM_R6 are lower than the ones of I2GG and I3GG, respectively. Under the condition that the linearization error is negligible, the measurement noise has the main influence on the gravity field determination. As mentioned in the "Methods" section, the theoretical computation indicates that the effect of measurement noise on I 3 is lower than that on I 2 for GOCE mission, which also holds true for the realistic data. In order to avoid the deviations from methodology and data processing, the comparison between the I3GG and I2GG is investigated in this section. From Table 1, it is shown that the difference RMS values of I3GG are smaller than that of I2GG when compared to the EGM2008 or EIGEN-5C. Specifically, for the common whole-frequency range (0-210), the ratio of the difference RMS values of the I3GG to I2GG is 87.3% (0.131/0.150) when the EGM2008 is used as the reference, while the ratio is 77.2% (0.088/0.114) when EIGEN-5C is used. The two ratios are bounced around the theoretical value 81.6% [ √ 4/3/ √ 2 from Eq. (14)], and their relative deviation to the theoretical value are about 7% and 5%, respectively. From Fig. 4, it is also illustrated that the maximum degree of I2GG is 220, which is lower than that of I3GG (degree 240). However, they are close to each other in the medium-frequency range, which is the signal sensitive band of GOCE mission and has the highest signal to noise ratio. The same rerult can be obtained by the comparison between I3GG_R2 and I2GG_R2. Such comparisons above show that both the invariant I 2 and I 3 can recover the gravity field effectively from the GOCE measurements under the sphere approximation. Additionally, compared to the linearization error, the comparison indicates that the measurement noise is the major factor which, affects the accuracy and the resolution of the gravity field determination.
An independent comparison with external data is made using geoid heights determined by GNSS positioning and leveling (GNSS/leveling). The details of the data processing procedures can be found in the study [20][21][22] . In order to verify and validate the method for the invariants, the models from the same GOCE data without  Considering the spatial resolution and the transition region effect, a Gaussian filter with a filter width of 47 km was applied to both the models and the GNSS/leveling data. Table 3 shows the results for the I3GG in comparison to other GOCE gravity field models using GNSS/leveling points of the European countries, including Norway, UK, Belgium, France, Germany and Netherlands. The RMS difference results show that I3GG has a relative comparative advantage over I2GG in Belgium and France.  . Difference degree amplitudes of the gravity field models I3GG and I2GG compared to the EGM2008. www.nature.com/scientificreports/ Meanwhile, in the rest countries they have alomst the same performance and their RMS differences are smaller than 1 cm . It should be pointed that the advantage of I3GG in the higher degrees is suppressed by the low-pass filter. Specifically, I3GG and I2GG fits better in comparison to the other models in UK. We also see that I3GG has an advantage over another second-invariant-derived model IGGT_R1 in all countries except Netherlands even though the IGGT_R1 applied regularization with the a-priori model EIGEN-5C. I3GG performs better in Norway, UK, Belgium, France and Germany of 0.3, 2.1, 2.9, 0.4 and 0.1 cm when compared with IGGT_R1. When compared with other GOCE-derived models DIR_R2, TIM_R5, SPW_R5, I3GG is in the middle level of the list in these countries. Briefly, it shows that the medium-to short-wavelength parts of the Earth's gravity field can be obtained effectively by using I 3 from GOCE.

Discussion
This proposal for the gravity field recovery from the invariant in satellite gradiometry is to avoid the errors in the reconstitution procedure of the satellite's attitude. In this paper, it is found that the third invariant I 3 is more competitive than I 2 on the gravity field determination GOCE mission, although I 2 is easier to be dealt with and get more attention from researchers. Accordingly, we use the third invariant I 3 of the GOCE gravity gradient tensor to obtain the I3GG. The theory of boundary value problems is adopted to implement linearization of the third invariant gravity gradients with an accuracy of O(J 2 2 · T ij ) . In order to solve the boundary value problems efficiently, we propose an alternative method of 2D FFT, which uses the coherent sampling theory to obtain the relationship between the 2D FFT and the third invariant I 3 measurements and applies the pseudo-inverse via QR factorization to transform the 2D Fourier coefficients to spherical harmonic ones. The complexity of this algorithm reaches a level of O(l 4 max ) without the precomputation of the integrals by recursions or numerical integrations. Based on the theoretical bases mentioned above, the GOCE gravity gradient data of the nominal mission phase and the lower degree information from EIGEN-5C were combined in the construction of a satellite-only gravity field model to a maximum degree of 240, corresponding to a spatial resolution of 83 km at the equator. The results indicate that the I3GG agrees well with the other GOCE gravity field models in the medium-frequencies range (degree 70-150), which is the signal sensitive band of the GOCE mission. For the lower frequencies range, I3GG is close to the EGM2008 and EIGEN-5C and their fitted ranges are up to degree/ order 70 and 100, respectively. From the point of view of the spatial domain, at regions where no or poor quality terrestrial data are available (Himalaya, Africa, Amazonas and Antarctica) new gravity field information is added. Moreover, for comparison purposes, we applied the same processing strategy on the second invariant I 2 measurements of the GOCE mission and obtained I2GG with a maximum degree of 220, which is 20 degrees lower than that of I3GG. Their RMS values of the geoid differences indicate that the effects of measurement noise of I3GG is about 20% lower than that of the I2GG model when compared to the EGM2008 or EIGEN-5C. The same characteristics can be also concluded when we dealt with the re-calibrated GOCE data. Briefly, the results show that the approach presented in this study is an effective way to obtain the gravity field models with high accuracy and spatial resolution from the third invariant of gravity gradient tensor. It is believed that it provides a viable alternative tool and viewpoint for studies on satellite gravity gradients.

Methods
There are two basic theoretical aspects for the gravity field determination from I 3 : the first is to formulate the relationship between the gravity potential coefficients and I 3 , and the second aspect is the 2D FFT for dealing with the boundary value problems. The former relates to the linearization of invariant and its error level, while the latter relates to the data processing flow and efficiency, which will be discussed in the following two subsections, respectively.

Invariant theory and the related boundary value problems.
In this subsection the theory of the gravitational gradient tensor's invariant and the related boundary conditions is introduced. The Earth's gravity potential V satisfies the Laplace equation and leads the gravitational gradient tensor to the symmetric and tracefree, and therefore the following representation for invariant system I 1 , I 2 , I 3 is adopted 15 : Both I 2 and I 3 , rather than I 1 , are generally considered to be suitable to recover the gravity field models since that the first invariant is the trace of the gravitational tensor and its value is zero. In the next step, the invariant are linearized by the perturbations calculation relative to a priori known reference solution subject to (1a) are the priori invariant from reference gravity gradients U ij , and O(T 2 ij ) and O(T 3 ij ) indicate the influences of terms which are non-linear in invariants. The disturbance gravity gradients T ij are defined as T ij = V ij − U ij . Eqs. (2) and (3) are the linearized equations for invariant with an accuracy of O(T 2 ij ) and O(T 3 ij ) , respectively 15 . Further researches focus on the simplicity of the linearization and reduces the computation to an efficient level with the sphere approximation 11,17 . The reference gravity potential U in the Eqs. (2) and (3) is approximated with Ũ = GM/r under the sphere approximation, and its second derivatives can be represented as where GM is gravitational parameter of the Earth (G is the universal gravitational constant and M is the mass of the Earth). This approximation makes δI 2 and δI 3 with the accuracy of O(J 2 · T ij ) and O(J 2 2 · T ij ) , respectively. Then one obtains the boundary value problems as follows 11 In this study the following second-order radial derivative for the Earth's disturbance potential T rr is adopted 23 : where (r, θ, ) are the geocentric spherical coordinates (radius, co-latitude, longitude), R is reference radius. l and m are degree and order of spherical harmonic, P lm (cos θ ) are the fully normalized Legendre functions, C lm and S lm are fully normalized spherical harmonic coefficients. The C lm and S lm are the unknown parameters that should be estimated from observations as the solution of the gravity field recovery problem. By combining the Eqs. (6) and (7) we have access to the gravity field determination from the third invariant.
In addition, we compare the effects of measurement noises of the gravity gradients on the invariant I 2 and I 3 under the sphere approximation. The effects of measurement noise on them are also different since that different combinations of the invariant I 2 and I 3 . Substituting Eq. (4) into (2) and (3) and ignoring the linearization errors O(J 2 · T ij ) and O(J 2 2 · T ij ) , one obtains: where T rr, I 2 and T rr, I 3 are second-order radial derivatives derived from the invariants I 2 and I 3 , respectively. Under the assumption that the measurement noises are uncorrelated 24 , error propagation then gives For the purpose of exploiting the influence of the measurement noises, here three situations are employed on the invariants: on or out of the sphere surface C lm cos m +S lm sin m P lm (cos θ ) www.nature.com/scientificreports/ 1. When the noise levels of T 11 , T 22 and T 33 are the same, it holds that σ := σ 11 = σ 22 = σ 33 , which results in 2. When the noise levels of T 11 and T 22 are twice that of T 33 , it holds that σ := 1 2 σ 11 = 1 2 σ 22 = σ 33 , which results in 3. When the noise level of T 33 is twice that of T 11 and T 22 , it holds that σ := σ 11 = σ 22 = 1 2 σ 33 , which results in From the above results, it is concluded that when the noise level of T 33 are not worse than that of T 11 and T 22 , the effect of measurement noise on I 2 is better than that on I 3 . On the contrary, when the noise level of T 33 is higher than that of T 11 and T 22 , the effect on I 3 is better than that on I 2 . In reality, the measurement noises of GOCE mission corresponds to the third situation 2 , i.e. σ 11 = σ 22 = 1 2 σ 33 = 10 mE/Hz 1/2 , which makes I 3 more competitive than I 2 on the gravity field determination from this point of view. This issue has been discussed from the results of realistic data in the Results section.
Principle of the alternative 2D FFT method. The 2D FFT method can be employed for the boundary value problems to determine the gravity field from the third invariant efficiently. As mentioned above, this method has two major steps, i.e. obtaining the 2D Fourier spectrum and transformation into spherical harmonic coefficients 18 . In this study an alternative 2D FFT method is developed to obtain the spherical harmonics from the third invariant in practice, which makes modifications to both steps.
First of all, we present a concise relation between the 2D FFT and the spatial signal based on the coherent sampling. The relation is developed from an explicit formula for reconstructing exactly a one-dimensional (1D) signal from the magnitude and phase of its FFT under the condition of coherent sampling, which is to avoid the picket fence and spectral leakage effects. In the case of a time series with a sampling rate of f s and N samples, while the value of the k-th sample of its 1D FFT is (a + ib) , the corresponding signal in time t domain s k (t) can be obtained as follows: where f k = (k − 1)f s /N is the frequency, PH k = arctan (b/a) is the phase and A k = √ a 2 +b 2

Z(k) is the magnitude with
It is noted that the first point ( k = 1 ) represents the 0 Hz frequecy. There is a similar formula for reconstructing a 2D signal from the magnitude and phase of its FFT. In the case of a spatial signal with a sampling rate of f sx and N sampling lines in x direction, and a sampling rate of f sy and M sampling lines in y direction, while the value located in the h-th row and k-th column of its 2D FFT is (a + ib) , the corresponding signal in space domain s h,k (x, y) can be obtained as follows: where f x = (h − 1)f sx /N and f y = (k − 1)f sy /M are the frequencies in x direction and y direction, respectively. The information of magnitude and phase are and respectively. The coefficient Z(h, k) is defined as Next, owing the fact that the Legendre functions can be expanded into sums of cosine or sine series, the same principles mentioned above can be also applied to discuss the correspondence between 2D FFT and spherical harmonics 25 . An analytical and square integrable function f (θ, ) defined on the unit sphere (0 ≤ θ ≤ π, 0 ≤ ≤ 2π) can be expanded in a series of spherical harmonics 26 (12) σ rr,�I 2 = 2 3 σ , σ rr,�I 3 = σ www.nature.com/scientificreports/ For the sake of clarity, we first discuss the correspondence for a specific value of l and m: Using trigonometric product sum identities, we can obtain the coefficient u q lm from the normalized Legendre function P lm (cos θ ) as follows 18 for m is even, and for m is odd, with It is shown that the phase is zero in the latitude direction and the magnitude u q lm can be obtained by a prior computation from the decomposition of the Legendre functions. In the case of spherical harmonics with a sampling rate of f sθ and N sampling lines in θ direction, and a sampling rate of f s and M sampling lines in direction, while the value located in the h-th row and k-th column of its 2D FFT is (a + ib) , the Fourier coefficients qm and qm are obtained as follows: with The frequencies in the h-th row and k-th column are f θ = q 2π and f = m 2π , respectively. Considering that only one component with a specific value of l and m is involved, the spherical harmonic coefficients C lm and S lm of Eq. (20) can be directly obtained from qm and qm Eqs. (21) and (22) indicate that the coefficient products qm u q lm is essentially the magnitude of spatial signal cos (m ) cos qθ or cos (m ) sin qθ , and qm u q lm is the one of spatial signal sin (m ) cos qθ or sin (m ) sin qθ . If only a spherical harmonic signal with a specific value of l and m is sampled under the condition of coherent sampling, the spherical harmonic coefficients C lm and S lm are equal to qm and qm , respectively. When sampling the signal contains multiple spherical harmonics in Eq. (19), the coefficients qm and qm may contain the joint information of Fourier spectrum of these components since the signal aliasing appears in the spherical harmonics with the same order.
Then, in order to transform the 2D FFT of the third invariant into spherical harmonic coefficients, we present an alternative algorithm based on the pseudo-inverse via QR factorization in this step. Using trigonometric product sum identities, we can first obtain the normalized Legendre function P lm (cos θ ) as follows 26 : where h lmk is the decomposition factor of the associated Legendre functions. For each order m the decomposition factors can be defined as a factor matrix as follows: (19) f (θ, ) = l max l=0 l m=0 C lm cos (m ) +S lm sin (m ) P lm (cos θ) (20) f lm (θ, ) = C lm cos (m ) +S lm sin (m ) P lm (cos θ)  www.nature.com/scientificreports/ Therefore, we can easily obtain the spherical harmonic coefficients C lm and S lm from the Fourier coefficients qm and qm by solving linear equations as follows: where C coef is the vector of spherical harmonic coefficients C lm and S lm , D FFT is the vector of Fourier coefficients qm and qm . Since H is left-invertible, its columns are linearly independent and the QR factorization H = QR exists. We have We can compute the pseudo-inverse using the QR factorization to obtain the spherical harmonic coefficients. By taking advantages of the orthogonality between the even degree Legendre functions and the odd degree ones in the same order, we can reduce the computation based on the classification by parity of their degrees. Considering Eq. (6), one can see that the third invariant is derived from f (θ, ) by adding a multiplication factor F(l) for a specific degree l, where Finally, it leads to an explicit expression for computing the gravity field from I 3 measurements based on the Fourier analysis of spherical harmonics as follows where F is a diagonal matrix constructed by multiplication factors F(l). The time complexity of the method presented in this study is mainly determined by the 2D-FFT computation and the pseudo-inverse via QR factorization. When the maximum degree is l max , the time complexity of the 2D-FFT computation is O(l 2 max · log(2l max )) . Meanwhile, the pseudo-inverse has a running time of O(l 3 max ) for each order, so that the total complexity for all orders is O(l 4 max ).
Realistic data processing. The components V xx , V yy , V zz and V xz are provided in the gradiometer reference frame (GRF) with highest accuracy of 10 to 20 mE. However, the accuracy of the other gravity gradients is reduced by around two orders of magnitude because they are measured by the less sensitive axes of the accelerometers 2 . Therefore, the gravity field determination from the third invariant depends on two kinds of data, i.e. the GOCE gravity gradient measurements and the synthetic gravity gradients derived from a priori gravity field. The GOCE gravity gradient measurements of the components V xx , V yy , V zz and V xz in the GRF are taken from the Level-2 product EGG_NOM_2 (ESA) from November 1, 2009 to August 1, 2012. The synthetic gravity gradients of the components V xy and V yz are derived from the global gravity field model EIGEN-5C. The synthetic gravity gradient values are first computed in the local north oriented frame (LNOR) and then trans- www.nature.com/scientificreports/ formed to the GRF by the attitude quaternions products. Due to the fact that the components V xy and V yz are very small compared to the main diagonal components, the accuracy of the synthetic gravity gradients satisfies the requirement for the gravity field determination from invariant 15 . On the other hand, the gravity field model EIGEN-5C is employed as the reference model and used to fill the polar gap since the inclination of the orbit of GOCE mission is 96.7 •8 . Meanwhile, the disturbance gravity gradients T ij can be also obtained by the subtraction of the reference gravity gradients U ij of EIGEN-5C from the GOCE gravity gradients measurements. Then the disturbance gravity gradients T ij are filtered to 5-100 mHz by a finite impulse response (FIR) filter, which means that the signals of the recovered gravity model at lower frequencies (i.e. below 5 mHz) mainly come from the reference model. The gravity gradients also need a vertical reduction from the realistic orbit height to the mean orbital sphere. Considering that the eccentricity of GOCE is about 0.001 and the difference between the realistic orbit height and the mean one is less than 4 km, according to Eq. (5), we can reduce the third invariant in the radial direction by using the Taylor series expansion where δr is the radial distances between the average sphere S and realistic orbit S. The error caused by this procedure is less than 1 mE and meets the requirement of gravity field determination from GOCE mission 27 . Then, the measurements are gridded with a spacing of 0.5 • × 0.5 • at the mean orbital sphere. For the sake of clarity, the associated flow chart is displayed in Fig. 5.