Gravity field recovery from orbit information using the Lagrange formalism

Gravity field reconstruction via the analysis of kinematic orbit positions has been proven to provide essential information for Earth system research purposes. For this aim, various approaches have been developed and applied to exploit kinematic orbits. In addition to those existing methods, in this paper we present a new technique. By means of a series of simulation studies we demonstrate that the novel method is comparable with the hitherto proposed techniques. As the main difference with existing methods, our approach is based on the so-called Lagrange coefficients, i.e., a semi-analytical description of the satellite motion. For this reason, we denote the technique to as the Lagrange formalism. The low sensitivity to the priori information about the gravity field, and less influence of the polar gap are of its characteristics. The investigations demonstrate that the idea of the Lagrange method in determining the Earth’s gravity field could represent comparable results in term of quality with other approaches.


Introduction
In the recent decades, a strong scientific interest has emerged aiming to better understand the physics of the Earth system.This interest implied (and partially still implies) the need for a tremendous improvement of our knowledge about the Earth's gravity field -with respect to both accuracy and (spatial) resolution.Today, thanks to the realization of a series of space gravimetry missions, we are in the fortunate situation having highly valuable data sets available, which allow us to stimulate progress in the Earth system research.Branches of science which make use of space-borne gravity field information include geophysical and geotechnical research [e.g., Tassara et al. 2007, Migliaccio et al. 2008], geodynamic research and seismology [e.g., Mantovani et al. 2001, Jordan and Watts 2005, Swain and Kirby 2006, Chen et al. 2007], oceanography [e.g.Bingham et al. 2011, Knudsen et al. 2010], hydrology [e.g., Tapley et al. 2004a, Schmidt et al. 2008, Krogh et al. 2010], cryospheric research [e.g., Chen et al. 2006, Velicogna 2009, Harig and Simons 2012, Jacob et al. 2012], and sea-level research [e.g., Riva et al. 2010, Baur et al. 2013, Jensen et al. 2013].
Space gravimetry is intrinsically connected to the satellite missions CHAllenging Minisatellite Payload [CHAMP, Reigber et al. 2002], Gravity Recovery And Climate Experiment [GRACE, Tapley et al. 2004b], and Gravity field and steady-state Ocean Circulation Explorer [GOCE, ESA 1999, Pail et al. 2011].According to the different observation principles underlying these projects, the collected measurements can be used to infer different ranges of the gravitational spectrum.In this contribution, we have a closer look on High-Low Satellite-to-Satellite Tracking (HL-SST), which is common to all three aforementioned missions.HL-SST means that a low-Earth orbiter is tracked by high-altitude Global Navigation Satellite System (GNSS) satellites such as the Global Positioning System (GPS).The technique allows to recover long-to medium-wavelength gravitational features up to spherical harmonic degree and order (d/o) of about 100 [e.g., Flechtner et al. 2010, Prange, 2011].
A variety of methods has been proposed to infer gravity field information by analyzing kinematic orbits derived from HL-SST observations.These methods include the energy balance approach [Han et al. 2002, Visser et al. 2003], the acceleration approach [Reubelt et al. 2003, Ditmar andvan Eck van der Sluijs 2004], the short-arc approach [Mayer-Gürr 2006], and the celestial mechanics approach [Prange et al. 2009].For more information about these methods, we refer the reader to the work by Baur et al. [2014] and the references therein.The study by Baur et al. [2014] showed that -apart from energy balance -the mentioned approaches provide equivalent results against the background of GOCE real data analysis; this finding is in agreement with what has been expected from a theoretical point of view.As pointed by Ditmar and van Eck van der Sluijs [2004], and later confirmed by others, the energy balance approach is characterized by its shortcoming due to the fewer number of the observation equations.However, it is valuable to notice that the results are obtained from the analysis of the GPS-derived kinematic orbit only.Then, it has not meant that the energy balance approach should be disregarded for any application, e.g., for GRACE LL-SST data analysis or computing the combined SST and SGG solution for GOCE DATA.
In this paper, we present a new approach for gravity field recovery from HL-SST derived orbits.It is based on the Lagrange formalism for satellite motion analysis [Curtis 2005, Sharifi andSeif 2011].The core feature of this new approach is that the unknown gravity field parameters (spherical harmonic coefficients) are computed from a semi-analytical description of the satellite motion.We compare the results obtained from the Lagrange method with those based on the (point-wise) acceleration approach [Reubelt 2009, Baur et al. 2012].In this context, the acceleration approach is representative for any of the hitherto proposed methods for kinematic orbit analysis -which, as outlined before, excluding energy balance have been shown to be equivalent in performance.
We aim at the ability and specific characteristics of the Lagrange formalism in the gravity field modeling in a closed-loop simulation procedure.In this approach, the derivative of the satellite orbit with respect to the force parameters is directly computed using the semi-analytical formulation instead of the variational equations.Unlike the variational equation, all derivations (i.e, design matrix entries) have been computed for all parameters altogether using the proposed formulation.Then, there is no need to solve a system of six differential equations or equivalently six definite integrals by numerical quadrature for each parameter as described in Beutler [2004] and Beutler et al. [2010].It could be very useful especially for modeling the gravity field up to high degree to avoid high computation costs.In this paper, any observable can be defined for modeling the Earth's gravity field based on the orbit perturbation in situ point-wise measurement, whereas Kaula's perturbation theory analyzes the accumulated orbit perturbation [Mayer-Gürr et al. 2005].There is an ability to write the observation equation by means of all elements of the state vector together.This characteristic could be an advantage if the satellite velocity would be directly observed in future missions for example using the Doppler shift measurements.In addition, by the help of the semi-analytical formulation of the Lagrange formalism, the orbit perturbation during a time sub-interval could be interpreted as mass distribution.This ability enables the proposed approach to model the regional gravity field.

Essentials of the Lagrange formalism
The Lagrange coefficients (f and g series) have initially been introduced for the solution of the twobody problem [Curtis 2005].These coefficients are according to a Taylor series expansion of the position and velocity vectors.The solution of Newton's equation of motion expressed by means of the Lagrange matrix L (s 0 , q) reads (1) In Equation (1), s 0 and s(t) are the initial state vector (position and velocity vectors at t 0 ) and the solution of the equation of motion (state at arbitrary point in time t), respectively; q represent force model parameters (for instance, gravity field spherical harmonic coefficients).In this paper, the position and velocity vectors in time t are respectively denoted by r(t) and ṙ(t).The solution of Equation ( 1) is equivalent to the determination of the Lagrange matrix.This matrix has the structure (2) where the arrays F and G contain the Lagrange coefficients for the orbital motion in a (non-Keplerian) gravitational field.Here and below in this research, dot and double-dot placed over variables denote first and second time-derivatives, respectively.These arrays and their time-derivatives read [Sharifi and Seif 2011].
Where coefficients F n (t) and G n (t) satisfy the relation as (4) The superscript (n) indicates the n-th time-derivative term, and subscript n denotes n-th term of the n-th term of the Taylor series in Equation ( 3).Generally, as it was shown in [Sharifi and Seif 2011], the n-th order derivative of the position vector has the ability to be decomposed into the position and velocity vectors.Then, r (n) (t) could be formulated as a linear combination of matrices multiplied by r(t) and ṙ(t).These matrices are denoted by F n (t) and G n (t) in Equation (4).
By differentiating Equation ( 4) with respect to time, we have based on the Equation ( 4) and rearranging yield a reformulated form of Equation ( 5) as ( 6) As an alternative, the (n + 1)-order term could be easily obtained by advancing the sequence in Equation ( 4) as (7) The comparison between Equation ( 7) and Equation ( 6), give the recursive coefficients definitions (8a) (8b) Equation ( 8) is a fundamental relation to compute F n (t) and G n (t).The initial terms could be easily computed based on the position and velocity vectors, F 0 (t) = I and G 0 (t) = 0; F 1 (t) = 0 and G 1 (t) = I, (0 is zero matrix and I represents the identity matrix).To find next terms of the sequence (n ≥ 2), we need some information about the gravitational acceleration acting on the satellite.
As it can be concluded from Equation ( 8), F 2 (t) and G 2 (t) coefficients are required for the next terms.These coefficients could be resulted from the acceleration vector decomposition into the velocity and position vectors according to Equation (4).Based on the relation represented by Sharifi and Seif [2011), the Earth gravitational acceleration could be formulated as: Then, it could be found from Equations ( 9) and ( 4) that F 2 (t) = (a − b) I + A and G 2 (t) = 0.In Equation ( 9), A is a diagonal matrix with entries (10) The unknown quantities (a, b, e, h, l) can be expressed as a function of the satellite position vector coordinates (X(t), Y(t), Z(t)) and the partial derivatives of the gravitational potential V with respect to the polar spherical coordinates radial distance r, latitude ϕ, and longitude λ: It is valuable to notice that as the scalars are functions of the satellite position time series, they are functions of the time too.
The coefficient matrices of the Taylor series Equation (3) are computed based on recursive formulation Equation ( 8), needed substitution have been considered For more details about the Lagrange formalism we refer the reader to the work by Sharifi and Seif [2011].This formulation has been used for the Earth's gravity fi eld modeling in the next section.

Gravity fi eld recovery
The motion of a satellite depends on the forces acting on it.For low-Earth orbiters, the dominating force can be attributed to the Earth's gravitational attraction.Expressing the gravitational potential by means of a spherical harmonic series with the degree-l and order-m series coeffi cients (also referred to as Stokes coeffi cients) C ¯ lm and S ¯ lm [e.g., Heiskanen and Moritz 1967], we fi nd q = {C ¯ lm , S ¯ lm } for the vector of force model parameters introduced in Equation ( 1).Given the a priori (or approximate) force parameters q 0 , we seek to fi nd corrections ∆q to q 0 in a way that the difference between the propagated orbit and the observed orbit becomes minimal in a least-squares sense.Denoting the estimated corrections by ∆q ˆ , the recovered gravity fi eld parameters obey (13) The optimization problem could not be optimally applied for entire time span because of the propagation error.Then, the long-arc orbit is subdivided into equidistant sub-intervals and this procedure is formulated for each sub-interval as demonstrated in Based on Equation ( 15), the observable ∆s ˆ (t n ) is written based on the difference between the observed and computed orbits for each sub-intervals.In this equation, the observed satellite state vector is denoted by s obs (t n ), and s com (t n , q) is the computed state vector obtained from the Lagrange method.
Since the sub-interval is usually larger than the defi ned step-size (h) for solving equations of motion, the equations of motion have to be solved using a step-bystep strategy with M = I n /h steps on the time interval [t n−1 , t n ].Then, the computed orbit at t n in a gravity fi eld with parameter with force model parameter q, denoted by s com (t n , q), reads (16) where s(t n − h) is the state vector computed at the previous step.Based on the strategy illustrated in Figure (1b), the s obs (t n −1 ) the sub-interval I n .A priori orbit is propagated using a priori force parameters q 0 from t n −1 to t n .The priori parameters are updated during the iterative process using the non-linear least square approach and the computed orbit becomes closer and closer to the observed one until the solution converges to least square solution.Linearization of Equation ( 15) yields ( 17) Inserting the linearized form of the state vector into Equation ( 15) yields (18) or equivalently (19) where ∆l is difference between the observed and computed orbit.As seen in Equation ( 19), we seek for values ∆q ˆ in such a way that the sum of squared inconsistencies (contained in the misfi t vector d ˆ ) fi nds its minimum.The solution to this least-squares problem reads [e.g., Koch 1999] (20) Now, let's have a closer look at the design matrix H.The design matrix consists of the derivatives of the satellite state vector with respect to the force model parameters.
The derivative of the i-th element of the satellite state vector at epoch t n with respect to q can be obtained by taking the derivatives of Equation ( 16) with respect to q (for summary t n − h is replaced by t´): The superscript "com" has been neglected for simplicity.After some reformulation and simplifications we find that  24) is the backbone for gravity field recovery using the Lagrange formalism.The derivatives of the Lagrange coefficients with respect to the Stokes coefficients are given in the Appendix.

Implementation of the proposed approach
In addition to the data sampling (dt), there are two time spacings needed to be set during the method implementation, step-size and time sub-interval.The orbit propagation step-size (h) in the Lagrange formalism is set to 5 s in this paper.The time sub-interval denoted by I n is needed to be defined in minimizing the problem based on Equation (15).Using the larger sub-interval I n in the proposed method leads to more difference between the computed and observed orbits, as illustrated in Figure (1a).It means more power of the observable signal, ∆s(t n ), that will be useful for modeling the gravity field especially higher-order harmonics.On the other hand, the larger sub-interval, the higher level of the prop-agation error in the Lagrange formalism.1 min is known for the sub-interval to be a proper trade-off between these factors.As described in Equation ( 16), we need the full satellite state vector at the start of each sub-interval to propagate the satellite orbit for sub-interval I n .The required velocity is numerically derived from the time series of the satellite position vectors by using 9-point Newton approach in case of the kinematic orbit.Although in the previous section the observation equations are formed for satellite state vector, only the first three terms of the state vector should be used for the gravity field modeling using the kinematic orbit.
For each sub-interval, three observation equations are written based on three position components.In order to form the design matrix, the partial derivatives of each observable equation with respect to all parameters are generated by using the semi-analytical formulation of the Lagrange approach.A set of observables is defined as time series of observables collected at time interval over the entire time span, I 0 , I 1 , ..., I K .The observation equations for sub-interval In are written by the help of the only observations at the beginning and the end of I n , s obs (t n -1) and s obs (t n ), without containing other ones within I n .Other observed orbit information would be imported into gravity field modeling problem by defining new sets of the observables.In order to involve all data, we need N = I n /dt observable sets.If the first sub-intervals of the first set starts at t 0 , one of the i-th set starts at t 0 + (i − 1)dt (i = 1, ..., N).The sub-intervals of the second set (shown with red color in Figure ( 1a)) is defined by dt time-shifting.As an example, for 1 min time interval, 2 sets are defined to involve all observables data with 30 sec sampling.

Numerical studies
The accuracy and efficiency of the Lagrange formalism for gravity field recovery was numerically assessed in the framework of a series of closed-loop simulation studies.All results are based on synthetic CHAMP and GOCE orbit data sets, which have been an outcome of the IAG Special Commission 7 (SC7) activities [Ilk et al. 2003].This standard data set was presented to provide a simple platform for any scientist to validate their approach.The unique data set is very useful to compare the various approaches for the gravity field modeling.The simulated orbits contain gravity field information up to d/o 300 using the EGM96 gravity field model [Lemoine et al. 1998].The data is produced using numerical integration over one month by neglecting the non-gravitational accelerations.The noiseless data set contains position, velocity vectors with a sampling rate of 5 s.As the main target of this manuscript is the gravity field recovery from the satellite kinematic orbit, only the position vectors were used as input information.In order to investigate the impact of observation noise on the gravity field inference results, we occasionally contaminated the (noise-free) SC7 satellite orbits with white noise with zero mean and standard deviation 1 cm in each coordinate.The simulation setup for GOCE and CHAMP has been listed in Table 1.
Apart from the absolute performance of our new HL-SST analysis method, we were particularly interested in its performance relative to the more established approaches.As outlined in Sect.I, for this purpose the acceleration approach has been selected as representative for any of the previously proposed methods.The acceleration approach results presented in this paper are based on the implementation as described in Baur et al. [2012] and Baur et al. [2014].The necessary accelerations are derived from the satellite orbit by means of the 9-point Gregory-Newton interpolation as the numerical differentiator.The numerical differentiation leads to the temporal correlation between the accelerations derived from the noisy orbit, even Gaussian noise.This correlation has been considered in the acceleration inversion based on the scheme represented in Baur et al. [2012], by help of the segmentation of the total number of observation and applying the analytical covariance propagation.
Results are displayed in terms of degree-error root mean square (spectral domain representation) (25) and geoid hight differences (spatial domain representation).Both were computed relative to the EGM96 model, which for our simulations is the "truth" model.Accordingly, ( holds true.For the sake of completeness, it should be mentioned that EGM96 signal amplitudes were computed by replacing ∆C ¯ lm and ∆S ¯ lm with and respectively, in Equation (25).

Impact of a priori gravity field information
As described in Sect.2.2, gravity field recovery using the Lagrange formalism is a non-linear problem, and hence requires a priori information (in the sequel referred to as initial gravity field, IGF) of the Stokes coefficients for linearization.In order to demonstrate that our Lagrange results are independent of the adopted IGF, we conducted the following experiment based on analyzing the noisy CHAMP orbit with 30 s sampling.All gravity models are resolved up to d/o 70.
In a first computation, we considered the GGM02S model [Tapley et al. 2005] up to d/o 0 as the IGF.In that case three iterations were sufficient to get the final solution.The change insolution between last two iterations was about 10 -15 in term of RMSl for each degree.In a second computation, we used GGM02S up to d/o 2 as the IGF.This time, two iterations were needed to end up at the final solution.Importantly, the final solutions are quite similar (Figure 2), which demonstrates that the choice of the IGF is uncritical.For convenience, we made use of the GGM02S model up to d/o 2 as the IGF for all our further investigations.

CHAMP results
Figure 3 shows the degree-error variances of recovered gravity fields from both the Lagrange method and the acceleration approach exploiting the noise-free CHAMP orbit.All results were derived from 30 days of position information.The originally 5 s sampled simulated observations were used for the studies.
It can be stated from Figure 3 that both approaches provide comparable results, whereas the solutions based on the acceleration approach are superior, especially in the short-and middle-wavelengths.The increase of the degree variances near the maximum    resolution of the recovered gravity fields (l >≈ 115) can be attributed to the spectral aliasing (impact of the omission error).
Although the two methods behave slightly different, the main message at this stage is that the impact of "model errors" are well below the impact of observation noise.In order to demonstrate this, we repeated the computations (for 5 s and 30 s sampling) with the noise-contaminated CHAMP data set.The results are shown in Figure 4.The most striking issue turning out from the analysis of the noisy data set is that the estimate based on the Lagrange formalism is clearly comparable with the acceleration approach solution.The EGM96 signal is intersected at l ≈ 50 and l ≈ 60 by the error curve of the solutions derived from the acceleration approach using the down-sampled and original orbit information, respectively.This result is clearly supported by the Lagrange findings.
In addition to the presentation of the degree-error variances, Figure 5  Finally, we assessed the estimated gravity fields in the spatial domain, namely in terms of geoid height differences (Figure 6).As expected, the pattern of the Lagrange method solution exhibits equivalent noise level with the acceleration approach pattern.

GOCE results
In order to substantiate our findings from the previous section, we additionally conducted experiments using the simulated GOCE data.Similar to Figure 4 for the CHAMP case, Figure 7 present the results for the case of GOCE.Because of the lower GOCE satellite altitude (250 km opposed to 500 km), and hence the possibility to recover higher-degree short-scale features, the gravity field was reconstructed up to d/o 90.
The sun-synchronous orbit design of the GOCE satellite implies a data gap in the polar regions.Therefore, GOCE-derived gravity field information near the poles is rather poor.Transferred to the spectral domain this means that the (near-)zonal Stokes coefficients can only be recovered with considerably reduced accuracy.For this reason, in Figure 7 coefficients with m < 5 are excluded.
For the solutions with maximum d/o 90 in Figure 7, again the Lagrange method result is comparable with the estimate making use of the acceleration approach (this is in agreement with the findings from our CHAMP investigations).
The acceleration approach error curve crosses the EGM96 signal curve at degree l ≈ 85; for the Lagrange formalism, the intersection is at degree l ≈ 89.Inspecting the error of each individual coefficient in Figure 8 suggests that the influence of the polar data gap is less pronounced for the Lagrange method.
The quality assessment of the recovered gravity fields has been carried out in terms of geoid height differences in the northern near-polar region too in Figure 9.In addition, RMS of geoid height difference along parallels has been shown for both approaches

Conclusions
We showed the ability of the Lagrange formalism for gravity field recovery from HL-SST derived kinematic orbits.From our simulation studies we conclude that this new method is comparable with the acceleration approach (which can be seen as representative for the more "conventional" orbit analysis methods).Based on one month of synthetic CHAMP orbit data (contaminated by noise), using the Lagrange method we found the signal-to-noise ratio to be 1 at spherical harmonic degree l ≈ 60.Making use of the acceleration approach, the noise exceeds the signal from similar degree.This qualitative finding has been by analyzing a synthetic GOCE data set.The maximum resolvable degrees were found to l ≈ 70; moreover, it turned out that the polar gap is much less pronounced for the Lagrange formalism gravity field solutions compared to the results exploiting the acceleration approach.
Future work includes the extension of the formalism with regard to the application of the method on real HL-SST orbit data sets.
Figure (1a).The n-th sub-interval is denoted by I n in Figure (1b).For the sake of simplicity, based on Einstein summation convention we rewrite Equation (1) by setting t n−1 as the initial time and t n as the desired time for sub-interval [t n−1 , t n ] as (14) Then, the minimization problem has been formulated as (15) Seif et al. [2011] introduced the Lagrangian state transition

Figure 1 .
Figure 1.Illustration of the Lagrange method configuration and the principles of time spacing.(a) The configuration of the Lagrange approach observable, the first (blue) and second (red) sets of the observables.(b) Principle of the sub-interval and step-size.

Figure 2 .
Figure 2. Degree-error variances of gravity field solutions dependent on a priori information for linearization.Results computed via the Lagrange formalism, analyzing the noisy CHAMP orbit with 30 s sampling.

Figure 3 .
Figure 3. Degree-error variances of gravity field solutions analyzing the noise-free CHAMP orbit with 5 s sampling (orginal signal).

Figure 4 .
Figure 4. Degree-error variances of gravity field solutions analyzing the noisy CHAMP orbit (σ = 1 cm in each coordinate) with 5 s and 30 s sampling.

Figure 5 .
Figure 5. Actual Error (two top sub-figures) Formal Errors (two bottom sub-figures) (in logarithmic scale) of spherical harmonic coefficients analyzing the noisy CHAMP orbit with 30 s sampling; (a and c) Lagrange formalism, (b and d) acceleration proach.

Figure 6 .
Figure 6.Geoid height differences (starting at d/o 2 up to d/o 45) relative to EGM96 analyzing the noisy CHAMP orbit with 5 s sampling; (top) Lagrange formalism, (bottom) acceleration approach.Gaussian smoothing with a radius of 500 km was applied to suppress high-frequency noise.
displays the error of each individual coefficient (difference between the "true" and the recovered Stokes coefficients in logarithmic scale).Figure 5 displays the formal errors (left column) besides the empirical errors (right column) of the estimated spherical harmonic coefficients.The patterns of the actual and formal errors are in a good agreement for both approaches.

Figure 7 .
Figure 7. Degree-error variances of gravity field solutions analyzing the noisy GOCE orbit (σ = 1 cm in each coordinate) with 5 s and 30 s sampling due to the polar data gap, orders m < 5 were omitted in RMSl computation.

Figure 9 .
Figure 9. Geoid height differences (starting at d/o 2 up to d/o 70) relative to EGM96 analyzing the noisy GOCE orbit with 5 s sampling for the northern near-polar region; (left) Lagrange formalism, (right) acceleration approach.