Low-degree coefficients of the conventional geopotential model
Within the process of deriving satellite accelerations from GPS observations, the reduced-dynamic Precise Orbit Determination (POD) approach is one of the most complete and accurate strategies. On the one hand, high-precision GPS solutions are independent of the user’s satellite dynamics but their solutions are more sensitive to geometrical GPS factors ( and ). On the other hand, a dynamical method estimates the satellite position and velocity at a single epoch for which the resulting model trajectory best fits the tracking observables. These dynamics involve the double integration and linearization of the Newton-Euler’s equation of motion, the force model, the parameters to be estimated and the orbital arc length. The combination of GPS observables with the dynamic approach can counterbalance the disadvantages of the dynamic miss-modeling and the GPS measurement noises (e.g., ) with the high GPS precision, the robustness of the dynamic approach and the use of additional force measurements and models (e.g., accelerometer measurements, time-varying gravity, irradiative accelerations, atmospheric drag, thruster firings, magnetic field, etc.). Since the POD products do not use to give information about the acceleration of the user’s satellite, accurate interpolation and subsequent numerical differentiation must be performed to the POD solution. For this purpose, a feasible methodology with the use of the arc-to-chord threshold is performed in this scheme of deriving satellite accelerations.
Non-gravitational accelerations are obtained by measuring the force to keep a proof mass exactly at the spacecraft’s center of mass, where the gravity is exactly compensated by the centrifugal force. Plus and minus drive voltages are applied to electrodes with respect to opposite sides of the proof mass, whose electrical potential is maintained at a dc biasing voltage. Unfortunately, this dc level is the source of bias and bias fluctuations of the most electrostatic space accelerometers. The non-gravitational forces acting on LEO satellites comprise the atmospheric drag, irradiative accelerations, thruster firings, the relative proof mass offset and the Lorentz force.
Several methodologies have been developed and applied for comparing the computed non-conservative accelerations with the accelerometer measurements. The non-conservative force models augmented with estimated empirical accelerations, for instance, have shown a good agreement with the accelerometer data  and were replaced with the accelerometer measurements in the reduced-dynamic POD, using a method with strong constraints in the cross-track and radial directions . In , the non-gravitational accelerations were calculated as piece-wise constant empirical accelerations via the reduced-dynamic POD approach with a standard Bayesian weighted least-squares estimator. Here, the regularization was applied to stabilize an ill-posed solution and only the longer wavelengths were recovered, at best in the along-track direction, with a bias in the cross-track direction. It was concluded that no meaningful solution could be obtained in the radial direction.
Concerning the acceleration approach, where accelerations are derived from a numerical differentiation along precise orbits, the second derivative of the Gregory–Newton interpolation scheme was used in , a 16-order polynomial in , and a 6-order-9-point-scheme in . After accelerations are calculated, the least-squares fitting can be applied, and the possible correlations between calibration parameters removed by using the (fitted) autoregressive-moving-average (ARMA) models .
However, results of the above mentioned studies were limited not only by the bias caused by applying a numerical differentiation, but also by the consequences of applying the least-squares estimation techniques. For example, the need of setting strong constraints or regularizations was probably caused by correlations from systematic errors. For this purpose, a new approach to accurately determine the GPS-based accelerations of GRACE mission was developed and applied in . Along with the benefit of using high accuracy and precise methodology and models, these systematic errors were modeled and removed from the computed non-gravitational accelerations of GRACE. In the following sections, this new approach is reviewed in detail and the calibration of GRACE accelerometers is given as an example of its application for LEO satellite measurements.
2. Methods and models
2.1. Accelerations from GPS-based precise orbit determination
First, a brief description of the GPS-based POD strategy is given to understand what kind of solutions is dealt with. Both GPS satellite (superscript s) and GPS receiver (subscript r) experience a clock offset (δt), causing their respective internal time delays in the overall GPS system time t as:
where is the signal travelling time. The number of carrier cycles between the reception and transmission of the signal is the phase difference plus an integer number of cycles
The components of the phase difference can be written in function of the track-continuity reference time t0 and frequency wave f as
The number of carrier cycles changes to
with the real valued parameter , which is constant over a continuous tracking.
The pseudo-range phase observation is obtained when multiplying the number of carrier cycles by the signal wavelength λ=c/f. Additionally, atmospheric effects , systematic errors and the thermal measurement noise must be taken into account.
and the circumflex mark indicates the value calculated in each iteration. The matrix adjustment v=[H] Δy for the linearized equation of observation is given for each row as
Note that is constant over a continuous tracking.
The solution must be solved in iterative process and null values can be used for Δy in the first iteration.
Through combinations of different frequency observations, several errors and parameters can be eliminated:
Ionosphere free linear combination.
Geometry-free carrier phase linear combination (to detect cycle slips).
Wide-lane and narrow-lane linear combinations (for ambiguity resolution applications).
Detailed methodology can be found, e.g., in  or . Differenced GPS observations has the advantage of eliminating or reducing several common error sources, such GPS satellite clock offsets and common biases from hardware delays. Distinguished by the carrier phase differencing level, high precision data processing procedures are basically:
The zero-differences or Precise Point Positioning (PPP) processing, where un-differenced observations from a single GPS receiver are supplemented by pre-computed IGS ephemeris and clock products (e.g., [12, 11]).
The single-differences processing, where the use of two receivers results in the cancelation of the satellite clock bias and common systematic errors .
Detailed developments and algorithmic implementation can be found, e.g., in . Regarding now the satellite dynamics, the first-order time-differential equations are
where the solution of a state vector y(t) is given by the position and velocity in function of time.
and p accounts for any force-model parameters to estimate. The residual vector is described by the difference between the measurements z and the modeled solution h, both as a function of time.
Reminding that the circumflex mark indicates the value calculated in each iteration, linearizing all quantities around a reference state
H gives the partial derivatives of the modeled observations with respect to the state vector at the reference epoch t0. The orbit determination problem is now reduced to the linear least-squares problem of finding Δy(t0). In first iteration, null values can be used, and the system solution is the correction to these values.
Note that H should be updated in each iteration, but can be replaced by an initial constant value, with the consequent iteration increasing and computation decreasing. Detailed algorithmic implementation and strategies can be found, e.g., in [22, 23, 11, 12, 24] or .
Since the POD products do not use to give information about the satellite’s acceleration, interpolation and subsequent numerical differentiation must be performed to the POD solution. In order to keep the error of interpolation small enough, a low-degree polynomial is not sufficient, high-degree polynomials introduce undesired oscillations, and the FFT approach is not considered in presence of data gaps and outliers . The best alternative, as seen in , is to use a piece-wise interpolation, such as splines or Hermite polynomials. Different algorithms for interpolation were compared in , interpolating odd from even original samples, and the committed error evaluated by a simple difference between interpolated and original data. Finally, 8-data point piece-wise Lagrange interpolation was chosen, which provided a white noise error of standard deviation of ~10nm/s, from evaluating the error committed at 10 s sampling (odd from even original samples). Similar results were obtained when the piece-wise cubic Hermite interpolation was tested.
When calculating total accelerations by simple differentiation of velocities, the approximations to numerical derivatives have been found to produce large bias . To avoid this error of arc-to-chord approximation , interpolated velocities must be differentiated by an increment of time (Δt), which minimizes the error committed at a given threshold. The modified three-point formula upgraded with the arc-to-chord threshold (Δt) is given by  in the form
where simple and double quotation mark denotes the simple or double arc to chord approximation respectively, and each ti is equispaced by Δt.
2.2. Time-varying gravity field effects
The conventional gravity model based on the EGM2008 , describes with Stokes’ coefficients the static part of the gravitational field and the underlying background for the secular variations of some coefficients. In addition, when computing the gravitational forces acting on the satellite, other time-varying effects must be taken into account. These include the third body tide caused by the Moon and Sun, the solid Earth tides, the ocean tides, the solid Earth pole tide, the ocean pole tide and the relativistic terms.
The geopotential field V at the point (r, φ, λ) is expanded in spherical harmonics with up to degree N as:
where t0 is the epoch J2000.0 and the values of , and rates of are given in Table 1
|Coefficient||Value at 2000.0||Rate (yr -1)|
To provide a mean figure axis coincident with the mean pole and consistent with the Terrestrial Reference Frame, the values for the coefficients and are given by:
Recent values of , and are adequate for 10−14 accuracy, e.g. the values of the present conventional model (−0.48416948 10−3, 2.4393836 10−6 and −1.4002737 10−6 respectively). The variables and (in radian) represent the IERS conventional mean pole:
where t0 is the year 2000.0 and the coefficients and (mas yr -i) are given in Table 2
|Before 2010.0||After 2010.0|
The gravitational acceleration of a third body on the satellite  can be described as a difference between the accelerations of the satellite and the Earth caused by a third body B.
where rsat and rB are the geocentric coordinates of the satellite and of a third body of mass MB, respectively. Since accelerations on near-Earth satellites from other planets actions are relatively small (< 0.1nm/s2), only Luni-solar accelerations are calculated. Moon and Soon coordinates can be interpolated from the solar and planetary ephemerides (DE-421) provided by the Jet Propulsion Laboratory (JPL) in the form of Chebyshev approximations. The evaluation of these polynomials yields Cartesian coordinates in the ICRS for the Earth-Moon barycenter b
where µ* denotes the ratio of the Earth's and the Moon's masses. Since the changes induced by the Earth’s solid tides  due to its rotation under effects of ellipticity and Coriolis force, can be described in terms of the Love numbers, variations in the low-degree Stokes’ coefficients can be easily computed. Dependent and independent frequency corrections are calculated using lunar and solar ephemerides, Doodson’s fundamental arguments, the nominal values of the Earth’s solid tide external potential Love numbers and the in-phase and out-of-phase amplitudes of the corrections for frequency-dependent Love values. First, changes induced by the tide generating potential in the normalized geopotential coefficients for both n = 2 and n = 3 for all m are given by the frequency-independent corrections in the form:
knm is the nominal Love number for degree n and order m,
rB is the distance from geocenter to Moon or Sun,
φB is the body-fixed geocentric latitude of Moon or Sun,
λB is the body-fixed east longitude (from Greenwich) of Moon or Sun.
Here, anelasticity of the mantle causes knm and k(+)nm to acquire small imaginary parts (Table 3).
|Elastic Earth||Anelastic Earth|
|n||m||knm||k(+) nm||Re(knm)||Im(knm)||k(+) nm|
To calculate rB, φB and λB, geocentric Moon and Sun Cartesian coordinates must be rotated from the ICRS to the ITRS and transformed to the spherical coordinates as usually.
Frequency dependent corrections are computed as the sum of contributions from a number of tidal constituents belonging to the respective bands.
is the six-vector of Doodson’s fundamental Lunisolar arguments (τ, s, h, p, N’, ps);
is the six-vector of multipliers of the fundamental Lunisolar arguments (j=a, b, c);
Aj is the In-phase (ip) and out-of-phase (op) amplitudes (j=a, b, c);
Doodson’s variables are related to Delaunay’s by
and Delaunay’s fundamental Lunisolar arguments can be computed as
Here, t is measured in Julian centuries of Barycentric Dynamical Time (TDB), but the Terrestrial Time (TT) can be used in practice (assuming a difference in the CIP location smaller than 0.01 µas).
Accounting for the dynamical effects of ocean tides, the periodic variations in the normalized Stokes’ coefficients are calculated based on the most recent ocean tide model EOT11a . The potential coefficients for the mass redistribution effect of ocean tides are available in the form of cnmCos, snmCos and cnmSin, snmSin (including the loading potential and the Doodson-Warburg phase corrections) up to maximum degree and order 120 for each tide s, and obtained by
 also provided the influences of additional minor tide constituents that are not included in the tide model EOT11a and should not be neglected in Low Earth Orbiters. This function evaluates the contribution of altogether 256 tides.
Changes in the geopotential value due to the centrifugal effect of pole motion, known as the Earth’s solid pole tides , can be readily computed in function of the wobble variables and calculated under sub-daily polar motion variations as
where m1 and m2 in seconds of arc are obtained from the difference between the polar motion and the IERS conventional mean pole (above defined) as
The standard pole coordinates of the parameters xp and yp are from the IERS (http://hpiers.obspm.fr/iers/eop/eopc04/) with additional components to account for the effect of ocean tides (∆xp, ∆yp)ocean and forced terms (∆xp, ∆yp)libration with periods of less than two days in space. These sub-daily variations are not part of the polar motion values published by the IERS and are therefore to be added after interpolation.
is the six-vector of Delaunay’s fundamental arguments (γ, l, l’, F, D, Ω);
is the six-vector of multipliers of the fundamental arguments (j=o, l);
are the amplitudes in xp for cosinus and sinus respectively (j=o, l);
are the amplitudes in yp for cosinus and sinus respectively (j=o, l);
j=o, l makes reference to oceanic and libration parameters, respectively.
Oceanic parameters are retrieved from Tables 8.2a and 8.2b, and libration parameters (n=2) from Table 5.1a of . The ocean pole tide, generated by the centrifugal effect of pole motion on the oceans, is calculated as a function of sub-daily wobble variables from the coefficients (and ) of the self-consistent equilibrium model . These perturbations to the normalized geopotential coefficients are given by
Only the main relativistic effects (described by the Schwarzschild field of the Earth itself ~1nm/s2), might be calculated, since the effects of the Lense-Thirring precession (frame-dragging) and the geodesic (de Sitter) precession are two orders of magnitude smaller at a near-Earth satellite orbit .
Finally, time-varying Stokes’ coefficients with up to a degree and order 120 can be computed under an increment of time small enough to desensitize from discontinuities (~3600 s) and the gravity for every satellite position calculated by using the first derivative of the gravitational potential in Cartesian coordinates. With the substitution of and , the first derivative of the gravitational potential of the Earth, in spherical coordinates, is calculated from the computed time-varying Stokes’ coefficients as
And the gravity
where the partial derivatives are
The derivative of the normalized associated Legendre function can be computed as
where the scale factor is
2.3. Accelerometer calibration
The approach here given (acceleration approach) is based on the fact that the GPS-derived accelerations can provide exact values to estimate the accelerometer calibration parameters. In order to compute the differences between the GPS-based non-gravitational accelerations and the accelerometer measurements, several transformations between reference systems are required (POD solutions are usually given in the ITRS and accelerometer measurements in the SBS).
Lunisolar direct tides and the relativistic effect are calculated in the ICRS, the others, in the form of Stokes’ coefficients. The computation must be done under an increment of time small enough to desensitize from temporal variations (~3600 s). Gravitational accelerations must be rotated to the ICRS as they were positions, because the Stokes’ coefficients have already included the contribution of the Earth’s rotation. The attitude of the satellites is based on celestial body observations, so the rotation to the SBS implies to pass through the ICRS. First, the gravitational accelerations must be subtracted from the GPS-based accelerations to obtain the non-gravitational accelerations. Then, computed non-gravitational accelerations can be rotated to the SBS, and the differences to the accelerometer measurements are evaluated by means of simple difference of median averages or by Least Squares adjustment.
For the following, the basic rotations in the tri-axial reference system are defined as
The rotation ICRS to SBS is derived from the star-camera quaternion (Fig. 1)
The transition from the ITRS to the ICRS is realized through a sequence of rotations that account for precession [P], nutation [N] and Earth rotation [S], including polar motion [PM] .
is the nominal mean Earth’s angular velocity corrected from Length Of Day (LOD). Similarly to polar motion, additional components should be added to the values from the IERS for UT1 and LOD, to account for the sub-daily effects of ocean tides.
where the effects of libration have been neglected due to its small size and
is the six-vector of Delaunay’s fundamental arguments (γ, l, l’, F, D, Ω);
is the six-vector of multipliers of the fundamental arguments;
are the oceanic amplitudes in UT1 for cosinus and sinus respectively;
are the oceanic amplitudes in LOD for cosinus and sinus respectively;
And the Earth rotation matrix [S] has the form
Polar motion consists of two quasi-periodic components and a gradual drift. The two main periodic parts are the Chandler wobble and the seasonal motions. The longer term variation is less well understood (motions in the Earth's core and mantle, water mass redistribution, and the isostatic rebound). The rotation matrix for polar motion is given by:
The standard pole coordinates xp and yp are from the IERS. The additional sub-daily variations of pole coordinates (Δxp and Δyp) and the parameter t (~Julian centuries of TT) have been formulated before.
The combined transformation for precession [P] and nutation [N] is defined in terms of co-declination , right ascension , and the instantaneous right ascension of the instantaneous pole (parameter ).
For convenience, the coordinates are redefined as
And the matrix [Q] becomes
The IAU 2006/2000A developments for the parameters X, Y and s are
is the 14-vector of multipliers of the fundamental arguments of the nutation theory;
is the 14-vector of fundamental arguments of the nutation theory, of which Lunisolar ones (l, l’, F, D and Ω) have been defined in the previous section, and the planetary ones (in radian) are:
The parameter t (~Julian centuries of TT) has been defined in the previous section. Amplitudes and the multipliers of the fundamental arguments are available at the IERS Conventions 2010 website (ftp://tai.bipm.org/iers/conv2010/chapter5/). After parameters X and Y are calculated, the IERS celestial pole offsets ΔX and ΔY can be added. These are provided by the IERS’ EOP.
3. Results and analysis
3.1. GPS-based non-gravitational accelerations from GRACE
The GNI_1A files provide the most precise position and velocity at 5 s interval, including formal error (dynamic-reduced approach computed by the GPS Inferred Positioning System software of JPL). GNV_1B files are practically unchanged from the previous processing level (GNI_1A files), since only smoothing on day boundaries is applied.
Time conversion between Level 1B files and the UTC is defined as
where TFILE is the time in the Level 1 data file and UTC0 is the January 1st of 2000 at 12h 00m 00s. Since satellite accelerations are not provided in the GNV_1B files, these are derived from the precise velocities by means of interpolation and subsequent numerical differentiation. 8-data point piece-wise Lagrange interpolation is used for interpolation as recommended by  or . In order to obtain a value to define the arc-to-chord interpolation threshold (Δt), total accelerations are calculated for several Δt and a value of Δt = 0.05 s should be chosen (error < 1nm/s2 in the arc-to-chord approximation). Then, total accelerations can be numerically obtained as the first derivative of the GNV_1B velocities.
After subtracting the time-varying gravity from the GPS-based accelerations, an unknown periodic error (~1.6 hour) of amplitude maxima of ~3µm/s2 can be identified in the YSBS axis of both GRACE satellites. ZSBS axes also reveal a slightly systematic behaviour, but two orders of magnitude smaller. The underlying signal is recovered by subtracting a sinusoidal function fitted on the envelope of the modulated amplitude. Since only sinusoidal functions are removed, the resulting solution remains unchanged from mean values . Figure 2 shows the results contrasted to accelerometer measurements calibrated with recommended a-priory biases of . When analyzing the solution for XSBS axes, aside the excellent agreement with the accelerometer precision, several local mismatches can be identified and probably attributed to the non-modelled local time-varying gravity, such as post glacial rebound, hydrologic cycle, etc., due to lack of accuracy of ocean tides models or possible external sources to the Earth’s gravity, e.g. Figure 2 at 19:30 h. Concerning the YSBS and ZSBS axes solutions, aside the bias detected with respect to , it can be clearly seen the high correlation with the accelerometer measurements, and consequently the good reference for accelerometer calibration.
3.2. Calibration of GRACE accelerometers
The twin satellites of the GRACE mission are equipped with three-axis capacitive SuperSTAR accelerometers to measure the non-gravitational forces. The precision of the XSBS and ZSBS axes is specified to be 0.1nm/s2, while 1nm/s2 for the YSBS axes. According to the ONERA Super Star accelerometers, the accelerations due to the relative proof mass offset and the Lorentz force can be neglected. Measurements at a second interval are included in ACC_1B files.
In order to evaluate the bias and bias fluctuations of the GRACE electrostatic accelerometers, the computed GPS-based non-gravitational accelerations are compared to the accelerometer measurements. Figure 3 shows the biases evaluated by simple difference of daily median averages for the 1st day and the 15th day of every month from 2003 to 2012 (240 days).
|Time span in MJD(UTC)|
|Time span in MJD(UTC)|
|Time span in MJD(UTC)|
Biases and fitted approximations are plotted in Figure 3, with respect the a-priori biases as recommended by . Since the nature of circular orbits implies a constant behaviour of the arc-to-chord error, the real magnitude (bigger) of radial accelerations (ZSBS axes) seems to cause a constant difference (~20 nm) in the solutions of . Furthermore, it is interesting to see that since electrostatic accelerometers are sensible to temperature changes, the correlation between YSBS biases and the β’ angle (angle between the Earth-Sun line and the orbit plane) is clearly recognized. Note that the β’ angle is defined such that it is zero when the Sun is within the orbit plane and, consequently, the perturbation of YSBS biases is minimized. The opposite situation happens when maximum β’ angle, in which the solar radiation has the same direction as the YSBS axes and maximizes its bias perturbation. These variations are disregarded in the polynomial fitting given in table 4, being this solution a more close approximation to values of β’ angle zero than the real β’ angle value.
In this chapter, a detailed theory and methodology to derive non-gravitational accelerations from GPS measurements has been review, developed and validated with GRACE satellites as an example. The results show good agreement with accelerometer measurements and demonstrate that this new approach is a good reference for accelerometer calibration. This methodology is based on the use of the arc-to-chord threshold for data interpolation, and the robust fitting of sinusoidal functions in case of finding systematic errors within the accurately derived non-gravitational accelerations.
With the development of space accelerometers, the calibration of current bias-rejection devices is not anymore required. Nevertheless, the non-gravitational accelerations can be determined accurately from the precise orbit ephemeris and, so far, it is, among others, a precious source of information for atmospheric and force model studies. With the increased accuracy of GPS positioning and dynamic measurements and models, the methodologies for GPS-based POD are becoming more and more accurate.
In order to guarantee an unbiased solution in accelerometer measurements, calibration parameters have been finally calculated without using any kind of regularization or constraint, and by using the GPS-based POD solution as a reference. Since the POD products do not use to give information about accelerations, the first derivatives of the precise-orbit velocities are calculated under the arc-to-chord interpolation-threshold. On the other hand, the modeled time-varying forces of gravitational origin and reference-system rotations are computed according to current IERS 2010 conventions (including sub-daily tide variations). In a case study from GRACE was shown that after subtracting the modeled time-varying gravity from the GPS-based accelerations, cross-track axes of both GRACE satellites are affected by a periodic error of unknown source. With the finality of extracting the underlying information from the resulting data, the systematic error is modeled and subtracted successfully by applying a robust fitting by sinusoidal functions. The resulting accelerations can serve as a reliable reference for the accelerometer calibration. In the future, this systematic error should have further considerations in POD software development.
ti; Instant i of time.
Δt; Increment of time.
c; Speed of light
U; Gravitational potential.
GM; Product of gravitational constant by mass.
ae ; Equatorial radius of the Earth.
Normalized associated Legendre functions of degree and order .
Normalized Stokes’ coefficient of degree and order for cosinus.
Normalized Stokes’ coefficient of degree and order for sinus.
Rotation vector of the Earth.
rSatellite position vector.
Satellite acceleration vector.
GPS; Global Positioning System.
POD; Precise Orbit Determination.
LEO; Low Earth Orbit
GRACE Gravity Recovery And Climate Experiment.
GA, GB; Satellite identifier (GRACE-A; GRACE-B).
EOP; Earth Orientation Parameters.
MJD; Modified Julian Date.
UTC; Coordinated Universal Time.
TAI; International Atomic Time.
ICRS; International Celestial Reference System.
ITRS; International Terrestrial Reference System.
SBS; Satellite Body System.
ORS; Orbit Reference System.
[P]; Rotation matrix for precession.
[N]; Rotation matrix for nutation.
[S]; Rotation matrix for sidereal time.
[PM]; Rotation matrix for polar motion.