InTechOpen uses cookies to offer you the best online experience. By continuing to use our site, you agree to our Privacy Policy.

Engineering » Aerospace Engineering » "Satellite Positioning - Methods, Models and Applications", book edited by Shuanggen Jin, ISBN 978-953-51-1738-4, Published: March 11, 2015 under CC BY 3.0 license. © The Author(s).

Chapter 3

GPS-based Non-Gravitational Accelerations and Accelerometer Calibration

By Andres Calabia and Shuanggen Jin
DOI: 10.5772/59975

Article top


Representation of the star-camera quaternion angles.
Figure 1. Representation of the star-camera quaternion angles.
Measured non-gravitational accelerations calibrated with recommended a-priory biases of [30] in solid line, and the GPS-based non-gravitational accelerations from this study in dashed line. GRACE-A on July 15th 2006. Plots are not equally scaled.
Figure 2. Measured non-gravitational accelerations calibrated with recommended a-priory biases of [30] in solid line, and the GPS-based non-gravitational accelerations from this study in dashed line. GRACE-A on July 15th 2006. Plots are not equally scaled.
Differences with respect to the solution of [30]. Small dots represent the standard deviation to the accelerometer measurements and solid line is the parameterized solution. Plots are not equally scaled.
Figure 3. Differences with respect to the solution of [30]. Small dots represent the standard deviation to the accelerometer measurements and solid line is the parameterized solution. Plots are not equally scaled.

GPS-based Non-Gravitational Accelerations and Accelerometer Calibration

Andres Calabia1, 2 and Shuanggen Jin1

1. Introduction

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 ([1] and [2]). 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., [3]) 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 [4] 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 [5]. In [6], 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 [7], a 16-order polynomial in [8], and a 6-order-9-point-scheme in [9]. 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 [9].

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 [10]. 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 τrs(t) 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  Nrs


The components of the phase difference can be written in function of the track-continuity reference time t0 and frequency wave f as

ϕr(t)=ϕr(t0)+(tt0) f+ (δtr(t) δtr(t0)) fϕs(tτrs(t))=ϕs(t0)+(tτrs(t)t0) f+ (δts(tτrs(t)) δts(t0)) f

The number of carrier cycles changes to

ϕrs(t)=τrs(t) f+[δtr(t) δts(tτrs(t))] f+Ars(t0)

with the real valued parameter  Ars, which is constant over a continuous tracking.

Ars(t0)=Nrs+ϕr(t0) δtr(t0) fϕs(t0)+δts(t0) f

The pseudo-range phase observation Lrs(t) is obtained when multiplying the number of carrier cycles by the signal wavelength λ=c/f. Additionally, atmospheric effects  Irs(t,f), systematic errors  Mrs(t) and the thermal measurement noise  εrs(t) must be taken into account.


where, ρrs(t)=c τrs(t)=  ǁrs(t- τrs (t))-rr(t)ǁ is defined as the geometric range between the GPS satellite and the receiver. The linearized carrier phase observation equation is written as:



(e^rxs,e^rys,e^rzs)=(x^r(t)xs(tτrs(t))ρ^rs(t), y^r(t)ys(tτrs(t))ρ^rs(t), z^r(t)zs(tτrs(t))ρ^rs(t))

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  Ars(t0)  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.

y=y^+Δy=y^+([H]t[H] )1[H]tv

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).

  • Multipath combinations.

Detailed methodology can be found, e.g., in [11] or [12]. 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  Mrs(t).

  • The double-differences processing, with the elimination of the receiver clock biases (e.g., [13, 14] and [15, 16, 17]).

  • The triple-differences processing, which eliminates the time-independent errors ([18, 19, 20]).

Detailed developments and algorithmic implementation can be found, e.g., in [21]. Regarding now the satellite dynamics, the first-order time-differential equations are

ddty(t)=f(t,y(t),p)=(r˙(t)r ¨(t,r,r˙,p))

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 ŷ(t0), the residual vector is approximately given by

vz(ti) h(ti,y^(t0))hy(t0)(y(t0)y^(t0))=Δ z(ti)HΔy(t0)

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 [25].

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 [26]. The best alternative, as seen in [7], is to use a piece-wise interpolation, such as splines or Hermite polynomials. Different algorithms for interpolation were compared in [10], 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 [9]. To avoid this error of arc-to-chord approximation [10], 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 [10] in the form

r ¨t0=limΔt0r ¨t0"=limΔt0r˙t1'r˙t(1)'Δt=limΔt0r ¨t22rt0r ¨t(2)(Δt)2

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 [27], 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:

V(r, φ, λ)=GMrn=0N(aer)nm=0n[C¯nmcos(mλ)+S¯nmsin(mλ)]P¯nm(sinφ)

The GM and ae values (398600.4415km3/s2 and 6378136.3m respectively) from the EGM2008 should be used as scaling parameters with its gravitational potential coefficients. In order to use the conventional static gravitational field properly and project it in time, the secular low degree C¯20 (zero-tide), C¯30 and C¯40 rates must be accounted for. These instantaneous values are given by:


where t0 is the epoch J2000.0 and the values of  C¯n0(t0), and rates of dC¯n0(t0)/dt  are given in Table 1

Coefficient Value at 2000.0 Rate (yr -1)
C¯20  (zero tide)-0.48416948·10−3 11.6·10−12
C¯30 0.9571612·10−6 4.9·10−12
C¯40 0.5399659·10−6 4.7·10−12

Table 1.

Low-degree coefficients of the conventional geopotential model

To provide a mean figure axis coincident with the mean pole and consistent with the Terrestrial Reference Frame, the values for the coefficients C¯21 and S¯21 are given by:

C¯21(t)=3 x¯p(t) C¯20x¯p(t) C¯22+y¯p(t) S¯22S¯21(t)=3 y¯p(t) C¯20y¯p(t) C¯22x¯p(t) S¯22

Recent values of  C¯20, C¯22 and S¯22 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 x¯p and y¯p (in radian) represent the IERS conventional mean pole:


where t0 is the year 2000.0 and the coefficients  x¯pi and y¯pi (mas yr -i) are given in Table 2

Before 2010.0 After 2010.0
Degree i x¯pi y¯pi x¯pi y¯pi

Table 2.

Coefficients of the IERS 2010 mean pole model

The gravitational acceleration of a third body on the satellite [22] can be described as a difference between the accelerations of the satellite and the Earth caused by a third body B.

Δr ¨sat=GMB(rBrsat|rBrsat|3rB|rB|3)

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 and the Sun b with respect to the barycenter of the solar system, while Moon positions r are given with respect to the center of the Earth. The geocentric position of the Sun can be computed as


where µ* denotes the ratio of the Earth's and the Moon's masses. Since the changes induced by the Earth’s solid tides [27] 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:

ΔC¯nmiΔS¯nm=knm2n+1B=,  GMBGM(aerB)n+1P¯nm(sinφB) eimλB   


ΔC¯4miΔS¯4m=k2m(+)5B=,  GMBGM(aerB)3P¯2m(sinφB) eimλB     m={0, 1, 2}


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

Table 3.

Nominal values of solid Earth tide external potential Love numbers

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);

n¯j 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);

j=a, b, c correspond to the parameters from Tables 6.5a, 6.5b and 6.5c given in [27].

Doodson’s variables are related to Delaunay’s by

τ= γ ss=F + Ωh=s  Dp=s  lN= Ωps=s  D  l

and Delaunay’s fundamental Lunisolar arguments can be computed as

γ= (67310.54841+(876600×3600+ 8640184.812866)t +0.093104t26.2×106t3)15 + 648000l= 0.00024470t4 + 0.051635t3 + 31.8792t2 + 1717915923.2178t + 485868.249036l= 0.00001149t4  0.000136t3  0.5532t2 + 129596581.0481t + 1287104.79305F= 0.00000417t4  0.001037t3  12.7512t2 + 1739527262.8478t + 335779.526232D= 0.00003169t4 + 0.006593t3  6.3706t2 + 1602961601.2090t + 1072260.70369Ω= 0.00005939t4 + 0.007702t3+ 7.4722t2  6962890.2665t + 450160.398036

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).

t= (JDTT 2451545)/36525JDTT=JDUTC+TAI  UTC+32.184

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 [28]. 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

Δc¯nm,s=cnmCos cos(n¯sβ¯)+cnmSin sin(n¯sβ¯)Δs¯nm,s=snmCos cos(n¯sβ¯)+snmSin sin(n¯sβ¯)

[28] 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 [27], can be readily computed in function of the wobble variables and calculated under sub-daily polar motion variations as

ΔC¯21=1.333 109(m1+0.0115 m2)ΔS¯21=1.333 109(m20.0115 m1)

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 ( 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.

Δxp=o[ xocoscos(n¯oF¯)+xosin sin(n¯oF¯)]+l[xlcos cos(n¯lF¯)+xlsin sin(n¯lF¯)]Δyp=o[yocos cos(n¯oF¯)+yosin sin(n¯oF¯)]+l[ylcos cos(n¯lF¯)+ylsin sin(n¯lF¯)]


F¯ is the six-vector of Delaunay’s fundamental arguments (γ, l, l’, F, D, Ω);

n¯j is the six-vector of multipliers of the fundamental arguments (j=o, l);

xjcos,xjsin are the amplitudes in xp for cosinus and sinus respectively (j=o, l);

yjcos,yjsin 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 [27]. 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 (A¯nm and B¯nm) of the self-consistent equilibrium model [29]. These perturbations to the normalized geopotential coefficients are given by



Rn=ω2 ae2GM4πGρwge(1+k'n2n+1)
ω=7.292115·10-5rad-1;ρw=1025kg/m3;G=6.67428·10-11m3/(kg·s2);ge=9.7803278m/s2;k’2=-0.3075, k’3=-0.195, k’4=-0.132, k’5=-0.1032, k’6=-0.0892;γ2R=0.6870andγ2I=0.0036.

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 [27].

Δr ¨sat=GMc2 rsat3{[4GMrsat(r˙satr˙sat)]rsat+4(rsatr˙sat)r˙sat}

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 P¯nm=P¯nm(sinφ) and P¯'nm=P¯nm(sinφ)/φ, 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)

rSBS=[Rib] rICRS=[q02+q12q22q322(q1q2+q0q3)2(q1q3q0q2)2(q1q2q0q3)q02q12+q22q322(q2q3+q0q1)2(q1q3+q0q2)2(q2q3q0q1)q02q12q22+q32]rICRS


q0=cos(ψ/2)     q1=sin(ψ/2)cosαq2=sin(ψ/2)cosβq3=sin(ψ/2)cosγ

Figure 1.

Representation of the star-camera quaternion angles.

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] [27].

rICRS=[P][N][S][PM]rITRSr˙ICRS=[P][N][S]{[PM]r˙ITRS+ω[PM]rITRS}r ¨ICRS=[P][N][S]{[PM]r ¨ICRS+ω(ω[PM]rITRS)+2ω[PM]r ˙ITRS}



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.

ΔUT1=oOUT1cos cos(n¯oF¯)+OUT1sin sin(n¯oF¯)ΔLOD=oOLODcos cos(n¯oF¯)+ OLODsinsin(n¯oF¯)

where the effects of libration have been neglected due to its small size and

F¯  is the six-vector of Delaunay’s fundamental arguments (γ, l, l’, F, D, Ω);

n¯o  is the six-vector of multipliers of the fundamental arguments;

OUT1cos,OUT1sin are the oceanic amplitudes in UT1 for cosinus and sinus respectively;

OLODcos,OLODsin are the oceanic amplitudes in LOD for cosinus and sinus respectively;

These oceanic parameters are retrieved from Tables 8.3a and 8.3b of [27].

And the Earth rotation matrix [S] has the form



ERA= 2π(0.7790572732640+1.00273781191135448JDUT);JDUT=JDUT12451545;UT1=UTC+ (UT1UTC) +ΔUT1;

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:

[PM]=R3(s) R2(xp+Δxp) R1(yp+Δyp)



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 d, right ascension  E, and the instantaneous right ascension of the instantaneous pole (parameter  s).

[Q]=[P][N]=R3(E) R2(d) R3(E) R3(s)

For convenience, the coordinates are redefined as


And the matrix [Q] becomes



a =0.5 + (X2+Y2)/8

The IAU 2006/2000A developments for the parameters X, Y and s are

X=106{16617''+2004191898'' t429782.9''t2198618.34'' t3+7.578''t4++5.9285''t5+j=04i[(ajsin)itjsin(A¯n¯) +(ajcos)itjcos(A¯n¯) ]}
Y=106{6951''25896''t22407274.7''t2+1900.59'' t3+1112.526'' t4++0.1358''t5+j=04i[(bjsin)itjsin(A¯n¯)+(bjcos)itjcos(A¯n¯) ]}
s=XY2+106{94.0''+3808.65''t122.68''t272574.11'' t3+27.98'' t4+15.62''t5+j=04i[(cjsin)itjsin(A¯n¯)+(cjcos)itjcos(A¯n¯) ]}


n¯ is the 14-vector of multipliers of the fundamental arguments of the nutation theory;

A¯ 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 (ajsin)i, (ajcos)i, (bjsin)i,(bjcos)i,(cjsin)i,(cjcos)i  and the multipliers of the fundamental arguments are available at the IERS Conventions 2010 website ( 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 [10] or [7]. 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 [10]. Figure 2 shows the results contrasted to accelerometer measurements calibrated with recommended a-priory biases of [30][30]. 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 [30], it can be clearly seen the high correlation with the accelerometer measurements, and consequently the good reference for accelerometer calibration.


Figure 2.

Measured non-gravitational accelerations calibrated with recommended a-priory biases of [30] in solid line, and the GPS-based non-gravitational accelerations from this study in dashed line. GRACE-A on July 15th 2006. Plots are not equally scaled.

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)
Axis 52720-53720 53720-55390 55390-55670 55670-56276
XSBS GAa1.8661E-107.6744E-112.6610E-09-2.1745E-09
XSBS GBa1.3180E-10-3.6244E-12-3.1007E-107.3652E-10
YSBS GAa-7.3899E-09-8.6972E-10-3.6298E-082.1178E-08
YSBS GBa-1.2166E-08-1.9175E-092.4402E-082.7737E-08
Time span in MJD(UTC)
Axis 52720-53005 53005-55166 55166-55562 55562-56276
ZSBS GAa2.5641E-094.1747E-116.9776E-10-7.2715E-12
Time span in MJD(UTC)
Axis 52720-53005 53005-55287 55287-55562 55562-56276
ZSBS GBa3.9394E-09-5.8487E-11-1.9218E-097.5564E-10

Table 4.

Fitted parameters for bias calibration of GRACE accelerometers

[i] - Equation: bias = ax2+bx+c, where x = (MJD(UTC)-55555)/100.

Biases and fitted approximations are plotted in Figure 3, with respect the a-priori biases as recommended by [30]. 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 [30]. 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.


Figure 3.

Differences with respect to the solution of [30]. Small dots represent the standard deviation to the accelerometer measurements and solid line is the parameterized solution. Plots are not equally scaled.

4. Summary

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.

5. Nomenclature

δjk;  Kronecker’s delta.

ti; Instant i of time.

Δt; Increment of time.

c; Speed of light

U; Gravitational potential.

; Earth.

⊙; Sun.

⊘; Moon.

GM; Product of gravitational constant by mass.

ae ; Equatorial radius of the Earth.

P¯lm;  Normalized associated Legendre functions of degree l and order m.

C¯lm;  Normalized Stokes’ coefficient of degree l and order m for cosinus.

S¯lm;  Normalized Stokes’ coefficient of degree l and order m for sinus.

ω;  Rotation vector of the Earth.

r;  Satellite position vector.

;  Satellite velocity vector.

r ¨; Satellite acceleration vector.

(φ,λ);  Latitude, longitude.

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.


1 - Jin S G, Luo O F and Gleason S. Characterization of diurnal cycles in ZTD from a decade of global GPS observations. J. Geodesy 2009; 83(6) 537-545.
2 - Jin S G, Dam T van and Wdowinski S. Observing and understanding the Earth system variations from space geodesy. J. Geodyn. 2013; 72 1-10.
3 - Jin S G, Zhang L J and Tapley B D. The understanding of length-of-day variations from satellite gravity and laser ranging measurements. Geophys. J. Int. 2011; 184(2) 651-660
4 - Helleputte T van, Visser P. GPS based orbit determination using accelerometer data. Aerospace Science and Technology 2008; 12 (6) 478–484.
5 - Helleputte T van, Doornbos E, Visser P. CHAMP and GRACE accelerometer calibration by GPS-based orbit determination. Advances in Space Research 2009; 43 1890–1896.
6 - Ijssel J van den, Visser P. Performance of GPS-based accelerometry: CHAMP and GRACE. Advances in Space Research 2007; 39 1597–1603.
7 - Reubelt T, Götzelmann M, Grafarend E W. Harmonic analysis of the Earth’s gravitational field from kinematic CHAMP orbits based on numerically derived satellite accelerations. In: Observation of the Earth System from Space. Berlin Springer, 2006. Part I 27-42.
8 - Ditmar P, Kuznetsov V, Sluijs A A V E van der, Schrama E, Klees R. DEOS CHAMP-01C70’: a model of the Earth’s gravity field computed from accelerations of the CHAMP satellite. Journal of Geodesy 2006; 79 586-601.
9 - Bezděk A. Calibration of accelerometers aboard GRACE satellites by comparison with POD-based nongravitational accelerations. Journal of Geodynamics 2010; 50(5) 410-423.
10 - Calabia A, Jin SG, and Tenzer R (2015) A new GPS-based calibration of GRACE accelerometers using the arc-to-chord threshold uncovered sinusoidal disturbing signal, Aerospace Sci. Technol. doi: 10.1016/j.ast.2015.05.013.
11 - Kroes R. Precise relative positioning of formation flying spacecraft using GPS. Dissertation, Delft University, Publications on Geodesy 61, Delft. 2006.
12 - Shabanloui A. A New approach for a Kinematic-Dynamic Determination of Low Satellite Orbits Based on GNSS Observations. Ph.D Thesis at the University of Bonn, Germany. 2008.
13 - Wu Y, Jin S G, Wang Z and Liu J. Cycle slip detection using multi-frequency GPS carrier phase observations: A simulation study. Adv. Space Res 2010; 46(2) 144-149.
14 - Švehla D and Rothacher M. Kinematic Orbit Determination of LEOs Based on Zero or Double Difference Algorithms Using Simulated and Real SST Data. Vistas for Geodesy in the New Millenium Springer Verlag 125 322–328. 2002.
15 - Švehla D and Rothacher M. CHAMP double difference kinematic POD with ambiguity resolution. In First CHAMP Science Meeting, Potsdam, Germany. 2003.
16 - Seeber G. Satellite geodesy, 2nd edition, Walter de Gruyter Berlin New York. 2003.
17 - Habrich H, Beutler G, Gurtner W, Rothacher M. Double difference ambiguity resolution for GLONASS/GPS carrier phase. Proceedings of ION GPS-99, 12th International Technical Meeting of the Satellite Division of the Institute of Navigation, Nashville, Tennessee, September 14–17. 1999.
18 - Reubelt T, Austen G and Grafarend E W. Harmonic analysis of the Earth’s gravitational field by means of semi-continuous ephemerides of a low Earth orbiting GPS-tracked satellite. Case study: CHAMP. Journal of Geodesy 2003; 77(5-6) 257–278.
19 - Byun S H. Satellite orbit determination using triple-differenced GPS carrier phase in pure kinematic mode. Journal of Geodesy 2003; 76 (9-10) 569–585.
20 - Grejner-Brzezinska D A, Shum C K and Kwon J H. A practical algorithm for LEO orbit determination. Journal of Navigation 2002; 49(3) 127–135.
21 - Hofmann-Wellennhof B, Lichtenegger H and Wasle E. Global Navigation Satellite Systems (GNSS). Springer-Verlag Wien New York. 2008.
22 - Montenbruck, O and Gill, E. Satellite orbits: Models, methods and applications. Berlin: Springer. 2013.
23 - Tapley B D, Schutz B E and Born G H. Statistical Orbit Determination. Elsevier Academic Press, Burlington, San Diego, London. 2004.
24 - Deo MN, Zhang K, Roberts C, Talbot NC. An investigation of GPS precise point positioning methods. Paper presented at SatNav 2003, 6th International Symposium on Satellite Navigation Technology Including Mobile Positioning & Location Services, Melbourne, Australia, July 22–25. 2003.
25 - Bisnath, S B. Precise Orbit Determination of Low Earth Orbiters with a Single GPS Receiver-Based, Geometric Strategy. Dissertation, Department of Geodesy and Geomatics Engineering, Technical ReportNo. 220, University of New Brunswick, Fredericton, New Brunswick, Canada. 2004.
26 - Weigelt M, Sneeuw N. Numerical Velocity determination and Calibration Methods for CHAMP Using the Energy Balance Approach. International Association of Geodesy Symposia 129 54-59. 2005.
27 - Petit G, Luzum B. IERS conventions (2010). IERS technical note 36, Frankfurt am Main. 2010.
28 - Rieser D et al. The ocean tide model EOT11a in spherical harmonics representation. 2012
29 - Desai S D. Observing the pole tide with satellite altimetry J. Geophys. Res. 2002; 107(C11) 3186
30 - Bettadpur S. Recommendation for a-priori Bias & Scale Parameters for Level-1B ACC Data 2009; GRACE technical note no. 2, V.2.