LOCAL ORBITAL FRAME OBSERVER FOR LEO DRAG-FREE SATELLITES

LOCAL ORBITAL FRAME OBSERVER FOR LEO DRAG-FREE SATELLITES Enrico S. Canuto 1, Luca Massotti2, Matilde Santos3 1

Politecnico di Torino, Dipartimento di Automatica e Informatica, Corso Duca degli Abruzzi 24, 10129 Torino, Italy [email protected] 2 European Space Agency,ESA-ESTEC/EOP-SFP, 2200 AG Noordwijk ZH, The Netherlands [email protected] 3 Dpto. Arquitectura de Comp. y Automatica,Universidad Complutense de Madrid 20040 Madrid, Spain [email protected] Abstract: This paper is concerned with on-board real-time CoM position and rate determination as input data to local orbital frame (LORF) for Low-Earth Orbit drag-free satellites. Study and simulation results are justified by Drag-Free and Attitude Control of the GOCE satellite (Gravity field and steady-state Ocean Circulation Explorer), the LORF being the reference for satellite attitude and scientific data. The paper focuses on modelling issues to obtain accurate orbit dynamics at lower frequencies for reducing integration error, because of the narrow-band filter requested by wide-band measurement errors. A further remedy in this sense is the addition of a second-order disturbance dynamics bounding unexplained noise components. Simulated results are presented with reference to GOCE mission. Copyright © 2004 IFAC. Keywords: Orbital frame, drag-free, satellite, state observer. 1. INTRODUCTION This paper is concerned with on-board real-time CoM position and rate determination as input data to local orbital frame (LORF) determination for LEO (Low-Earth Orbit) drag-free satellites. Methods and results can be extended to other orbits and satellites. Study and simulation results are justified by DragFree and Attitude Control (DFAC) of the GOCE satellite (Gravity field and steady-state Ocean Circulation Explorer), the LORF being the instantaneous reference for satellite attitude and scientific data (Canuto, 2003, 2007b, 2007c). Accuracy requirements directly descend from attitude and ask for microradian real-time errors in the socalled mission Measurement Bandwidth (MBW) from 1 mH to 0.1 Hz to be compared with an overall sub-milliradian error range, because of relaxed lowfrequency attitude below 1 mHz. Such a discrepancy depends on the gravity field measurement principle combining a posteriori precise orbit determination (POD) from on-board GPS data and gravity gradient

measures from on-board ultra-precise gradiometer. POD will provide low-frequency (hence long spatial wavelength) gravity components to be corrected by the fine spatial structure made available by gradiometer. The LORF provides the orientation of the osculating instantaneous orbit with respect to an inertial frame as the Equatorial Earth centered at some date (usually Julian 2000) or with respect to a mean Earth-fixed circular orbit only determined by mean Earth gravity, being the satellite drag-free. In the latter case, Euler angles of the LORF-to-mean orbit transformation can be shown to be combination of the ratios between CoM position and rate perturbations with respect to mean orbit radius r 6.6 × 106 m and absolute orbital speed v 7.8 × 103 m/s , respectively, the latter values pertaining to GOCE altitude around h = 250 km . Direct LORF determination in the order of microradian would imply CoM position and rate measurement errors below few meters and 1 cm/s , respectively, at a sampling rate greater than 0.1 Hz. Although space-borne GPS receivers approach this

limits at a sampling rate of 1 Hz, position and rate measurements need be filtered to cope with nonstationary GPS errors and to extrapolate measurements during sampling time as requested by higher control rates > 1Hz, as in GOCE mission. The objective is to dispose of a robust algorithm, free of measurement statistics, and tuned to quite conservative GPS errors (Kaplan, 1996) as shown in Table 1 . Table 1

GPS timing and errors.

Parameter Sampling time Position error Rate error Bias Delay

Unit s m m/s

Value 1 <30 (1σ) 0.03 (1σ) negligible negligible

2. REFERENCE FRAMES AND REQUIREMENTS The inertial frame RJ = {O, i J , jJ , k J } centered in the Earth COM O , is the equatorial frame at the date J2000. The LORF RO = {C , i O , jO , k O } , centered in the spacecraft CoM C , is defined by the instantaneous orbit orientation, and specifically by the motion direction v / v of the CoM, v being the inertial velocity, and by the orbital plane orthogonal to the angular momentum h = ms r × v , ms denoting spacecraft mass and r CoM position (Fig. 1). The LORF is the reference frame for science measurements and attitude control. The LORF axes are defined by: i O = v / v , jO = r × v / r × v , k O = i O × jO . (1) The axes from i O to k O will be referred to as alongtrack, out-of-plane and radial, respectively. v

iO kJ

kO i

k

CoM

jO

Assuming a quasi-circular circular orbit, eO can be shown to be related to inertial position and velocity estimation errors Δr and Δ v through − Δr ⋅ jO 1 eO ≅ Δr ⋅ i O − Δ v ⋅ k O / ωO , r = r , v = ωO r = v , r Δ v ⋅ jO / ωO

where ωO = μ / r 3 1.2 mrad/s is the mean orbital rate at GOCE altitude. Assume now the LORF matrix to be directly measured from GPS measurements y r (t j ) = r (t j ) + vr (t j ) , (4) y v ( t j ) = v ( t j ) + v v ( t j ) , t j = jTg

which are collected from GPS or similar satellite-tosatellite receiver at uniform sampling rate f g = 1/ Tg = 1 Hz . Properties of GPS measurement errors v r and v v are not recalled here (Kaplan, 1996), but they are assumed to be discrete-time (DT) white and Gaussian noise, and their components to be statistically independent. Assume position and velocity scalar errors to have different dispersions bounded as in Table 1 , denoted with σ r and σ v respectively. Then, assuming equality Δr ( t j ) = v r ( t j ) , Δ v ( t j ) = v v ( t j ) (5) in (3) yields the following white noise PSD SOx SOy ( f ) ≅ SOz

O

Fig. 1.

LEO polar orbit and reference frames.

The matrix R O = [ i O jO k O ] accomplishes the LORF-to-inertial coordinate transformation, and defines the reference attitude to be tracked by the spacecraft during the mission. The orientation error ˆ , can be eO , caused by the on-line estimate R O defined as ˆ , E ≅ I + ΔE R O EO = R O O O 0 ΔEO = eOz −eOy

−eOz 0 eOx

eOy eOx . −eOx , eO = eOy 0 eOz

σ + (σ v / ω O ) 2 r

σ v / ωO

2

6.5 μ rad ≤ 8.5 Hz (6) 5.5

Bounds to LORF errors and attitude.

LORF error Attitude LEO polar orbit

r

σr

for each LORF component error. Notice PSD to mean the root of the unilateral power spectral density. Requirements to LORF estimation errors come from spacecraft attitude as the latter is defined by the body-to-LORF transformation. LORF errors should be a fraction of the residual attitude as shown in Table 2 , derived from GOCE requirements.

Variable

iJ

2Tg

f ≤ 0.5 f g = 0.5Hz

Table 2 jJ

(3)

Overall RMS [μrad] 200 370

MBW

PSD

⎡μrad/ Hz ⎤ ⎣ ⎦ 4.3 7.9

The most stringent requirements in Table 2 occur in the so-called MBW from 1 mHz to 0.1 Hz. The direct LORF measurement, being not compliant according to (6), justifies a real-time filter guaranteeing requirements with some margin. The filter will be designed in the form a state observer around a stylized DT dynamics called Embedded Model, as suggested by (Canuto 2004, 2007a). 3. THE EMBEDDED MODEL

(2)

First, orbit dynamics will be derived in continuoustime. The DT model will be derived for each inertial

coordinate by exploiting their weak interaction due to quasi-circular orbit and quasi-spherical gravity field. 3.1 Orbit dynamics

and to longitude λ and latitude θ , through the Earth’s rotation, such that z = r sin θ . Denote the Earth’s gravitational potential with U ( ρ ) , which is usually expanded into complex spherical harmonics Ynm (θ , λ ) of degree n and order m, scaled by their complex spectrum K nm , n +1 n μ ∞ ⎛ R⎞ U ( r ,θ , λ ) = ∑ ⎜ ⎟ ∑ K nmYnm (θ , λ ) (8) R n = 0 ⎝ r ⎠ m =− n where R is the mean Earth radius, μ is the Earth gravitational constant, and r = R + h , h being the CoM altitude. The fine structure of fluctuations according to Kaula (Bertotti and Farinella, 1990) mainly depends on the degree n and is approximated by the so-called Kaula’s rule 10−5 n −2 , n >> 1, K n =

∂Ω ( r ) = 3 (1 − r / r ) I + γ 0 ( z / r ) I − Γ 1 .(12) 2

15 3 2 2 J2 ( R / r ) , Γ 1 = J2 ( R / r ) Γ 2 2 From (12), the perturbation term ∂Ω ( r ) of the orbit rate is clearly bounded by (13) ∂Ω ( r ) ≤ 3 ( ε + 2 J 2 ) 0.015 ,

γ0 =

Consider a LEO satellite where the CoM inertial position r and rate v = r can be measured by a GPS receiver at a uniform sampling rate f g . Assume the satellite is free-falling, i.e. a drag-free control cancels non-gravitational forces below a certain threshold. To this end, denote the residual non gravitational CoM acceleration, ideally kept to zero, with a . T The inertial position rT = [ x y z ] T = [ r1 r2 r3 ] , is related to Earth fixed coordinates ⎡cos λ cosθ ⎤ ρ = r ⎢⎢ sin λ cos θ ⎥⎥ , r = r , (7) ⎢⎣ sin θ ⎥⎦

Kn

g ( r ( t ) ) = −ωO2 ( I + ∂Ω ( r ) ) r ( t ) + δ g ( r )

n

( 2n + 1) ∑

m =− n

K nm

2

. (9)

By keeping the spherical and 2nd order terms due to Earth flatness in (8) and by confining higher order terms into a single anomaly δ U , one obtains 2 2 μ ⎛ 3 ⎛ R ⎞ ⎛⎛ z ⎞ 1 ⎞⎞ U ( r , z ) = ⎜1 − J 2 ⎜ ⎟ ⎜ ⎜ ⎟ − ⎟ ⎟ + δ U , (10) r ⎜⎝ 2 ⎝ r ⎠ ⎜⎝ ⎝ r ⎠ 3 ⎟⎠ ⎟⎠ where J 2 0.001 and the explicit term is now independent on the Earth’s rotation. The gravity acceleration g = ∇U is derived from (10) and the explicit component can be expressed into inertial coordinates as follows 2 2 ⎞⎞ μ ⎛ 3J ⎛ R ⎞ ⎛ ⎛ z ⎞ g = − 3 ⎜ I − 2 ⎜ ⎟ ⎜ 5⎜ ⎟ I − Γ ⎟⎟ r + δ g ⎟⎟ 2 ⎝ r ⎠ ⎜⎝ ⎝ r ⎠ r ⎜⎝ , (11) ⎠⎠ Γ = diag {1,1,3} and the anomaly δ g ( r ) corresponds to δ U . A further simplification comes by assuming a quasi circular orbit as in GOCE, i.e. having a low eccentricity ε < 0.003 . The explicit terms in (11) are expanded around r with the care of keeping perturbation terms of the same order of magnitude as J 2 and of confining residuals into δ g ( r ) . That amounts to 1st order expanding the spherical term and to zero-th order expanding the flatness terms, which gives the final expression

and it will be assumed to be periodic in time and to possess a long-term average. Eccentricity and flatness contributing by the same order of magnitude justify (12). Higher order terms would contribute to less than 0.0001. Orbit dynamics can then be written as 0 I ⎤ ⎡r ⎤ ⎡ ⎡r ⎤ ⎡0 ⎤ ⎢ v ⎥ ( t ) = ⎢ −ω 2 ( I + ∂Ω ( r ) ) 0 ⎥ ⎢ v ⎥ ( t ) + ⎢ I ⎥ a d (t ) , ⎣ ⎦ ⎣ ⎦ ⎣ O ⎦⎣ ⎦ ad (t ) = a (t ) + δ g ( t )

(14)

having combined the residual non gravitational acceleration and gravity anomalies into a d . CoM coordinates in (14) are each other connected through the weak perturbing term ∂Ω ( r ) . 3.2 Discrete-time dynamics Embedded model according to Canuto (2004, 2007a) is made by two interconnected parts: the controllable dynamics and the disturbance dynamics. The orbit being drag-free the input to controllable part is assumed to be quasi-zero except for gravity acceleration. The latter is partitioned into a position feedback to be part of the controllable dynamics (spherical and J2 terms) and residual components treated as an unknown disturbance driven by arbitrary signals. Consider a time unit T = Tg / ng , with ng ≥ 1 and assume T << TO = 2π / ωO , which implies (14) to be LTI during T and the discrete-time version to become time varying. The generic discrete time will denoted by ti = iT . In the following we shall limit to a single generic CoM coordinate rk , k = 1, 2,3 with increment vk [m], collected into the state xTk = [ rk vk ] . To this end, denote the current angular increment associated to rk with α ki = ω kiT , where

{

}

(

)

diag α12i , α 22i , α 32i = ω O2 T 2 I + ∂Ω ( r ( ti ) ) .(15)

The discrete-time equation holds sin α ki ⎤ ⎡ cosα ki ⎡ rk ⎤ ⎡r ⎤ ⎢ α ki ⎥⎥ ⎢ k ⎥ ( i ) + ⎢v ⎥ ( i + 1) = ⎢ v ⎣ k⎦ ⎢⎣ −α ki sin α ki cosα ki ⎥⎦ ⎣ k ⎦ (16) ⎡ sin (ω ki ( ti +1 − τ ) ) ⎤ ti +1 ⎢ ⎥ +∫ ⎢ ω ki ⎥ adk (τ ) dτ ti ⎢T cos (ω ( t − τ ) ) ⎥ ki i +1 ⎣ ⎦

where adx is a coordinate of a d . By exploiting T << TO = 2π / ωO and guaranteeing unit circle eigenvalues, (16) can be further simplified for real—time computation by replacing

trigonometric functions with polynomial expansions up to 2nd order terms in α ki , which provides ⎡1 − α ki2 / 2 1 − α ki2 / 4 ⎤ ⎡ rk ⎤ ⎡ rk ⎤ i + = 1 ( ) ⎢ ⎥ ⎢ ⎥ (i ) + ⎢ ⎥ 2 1 − α ki2 / 2 ⎦ ⎣vk ⎦ ⎣ vk ⎦ ⎣ −α ki . (17) ti +1 ⎡ ( t i +1 − τ ) I ⎤ +∫ ⎢ ⎥ adk (τ ) dτ ti T ⎣ ⎦ Integration in (17) is solved (i) by defining a disturbance increment d k (i ) , in length units [m], holding d k (i ) = T ∫

ti +1

ti

adk (τ ) dτ ,

(18)

and (ii) by dropping the direct effect of adk on as the former just vk (i ) = rk ( i + 1) − rk (i ) , corresponds to time integration of the unique disturbance source d k (i ) . The result is the DT version of a time-varying oscillator having very long period with respect to time unit T and subject to an acceleration disturbance: x k ( i + 1) = Ack ( i ) x xk ( i ) + Bc d k ( i ) ⎡1 − α ki2 / 2 1 − α ki2 / 4 ⎤ ⎡ 0 ⎤ .(19) Ack ( i ) = ⎢ ⎥ , Bc = ⎢ ⎥ 2 2 1 − α ki / 2 ⎦ ⎣1 ⎦ ⎣ −α ki Note that Ack (i ) is slowly time-varying with the orbit angular increment α ki . The LORF observer (Section 4) will be designed to have steady eigenvalues through time-varying feedback gains. Equation (19) must be completed with disturbance dynamics as below.

⎡1 Ad = ⎢ ⎣0 ⎡0 Hc = ⎢ ⎣1

1⎤ ⎡0 , Gd = ⎢ ⎥ 1⎦ ⎣0 0⎤ ⎡0 , Gc = ⎢ ⎥ 0⎦ ⎣1

1 0⎤ 0 1 ⎥⎦ . 0 0⎤ 0 0 ⎥⎦

(21)

2nd order dynamics (20) is in agreement with Kaula’s rule since Fourier frequency f [Hz] is related to n by f = nωO / ( 2π ) = nfO 0.19n , (22) which implies the time profiles of the Earth’s gravity potential and accelerations to show a PSD rolling off at -40 dB/dec just after resonances at the 1st and 3rd harmonics of fO due to J2 terms. Resonances are proven by restricting (12) to J2 terms and assuming a circular polar orbit, corresponding to z ( t ) = r sin (ωO t ) in (12) and θ ( t ) = sin (ωO t ) in (7) less a phase ϕO = 0 set to zero for simplicity, which leads to ⎡cos λ cos ωOt ⎤ g J 2 ( t ) = −ωO2 r γ 0 I sin 2 (ωOt ) − Γ 1 ⎢⎢ sin λ cos ωOt ⎥⎥ ⎣⎢ sin ωOt ⎦⎥

(

)

(23).

Expansion of trigonometric terms shows 1st and 3rd harmonics alone to appear as expected. Fig. 2 shows the PSD of the gravity acceleration anomalies in J2000 coordinates from measured coefficients up to n = 36 (to limit simulation time over several days) and extrapolation from Kaula’s rule. Fig. 2 fully confirms the adopted EM.

3.3 Disturbance dynamics Disturbance dynamics is synthesized starting from the experimental PSD of the CoM disturbances (Canuto 2007a, 2007c). Residual non-gravitational acceleration should be retained a wide-band noise within the performance BW, except at lower frequencies, where lateral components might remain uncompensated to save propellant (Canuto, 2007b). The noise will be modelled as DT white noise w0k ; drifts can be included in the gravity anomaly dynamics as follows. Gravity anomalies and the modeling errors of spherical and J2 terms resulting from approximations leading to (19) are modelled from the spectral density of the gravity anomalies rolling off at -40dB/dec beyond 3rd order harmonics as in Fig. 2, from a LEO polar orbit at h = 250 km . To this end, the total disturbance d k is decomposed into the sum of white noise w0 x and the combination of random drifts zTk = [ z1k z2 k ] , driven by a pair of white noise w1k and w2k , which are collected together with w0k into wk . The following 2nd order dynamics, to be replicated for each coordinate, follows z k (i + 1) = Ad z k (i ) + Gd w k (i ) , (20) Bc d k (i ) = H c z k (i ) + Gc w k (i ) and the relevant matrices hold

Fig. 2.

PSD of the gravity acceleration anomalies.

3.4 Embedded model and model error The complete EM of a single coordinate is the combination of (19) and (20), and is rewritten below by dropping the subscript k and adding the subscript m to indicate EM: ⎡ Ac ( i ) H c ⎤ ⎡ x m ⎤ ⎡ xm ⎤ ⎡ Gc ⎤ ⎥ ⎢ ⎥ (i ) + ⎢ ⎥ w ( i ) ⎢ z ⎥ (i + 1) = ⎢ Ad ⎦ ⎣ z m ⎦ , (24) ⎣ m⎦ ⎣Gd ⎦ ⎣ 0 y m ( i ) = x m (i )

where Ac (i ) = Ac ( rm ( i ) ) depends on the state through the term α i2 , rewritten from (15) as

α i2 = α 02 (1 + ∂Ω ( rm ( i ) ) ) .

Equation size is denoted by

(25)

n y = dim y m = 2, nw = dim w = 3

.

(26)

nc = dim x = 2, nd = dim z = 2 Denote with y ( i ) the sampled position and rate as measured by the GPS receiver. They are related to y m by the model error e = y − y m , which, according to Canuto (2007a), is the composition of neglected dynamics ∂P ( y m ,...) in fractional form and residual noise v including GPS measurement noise. The model error relation holds y ( i ) = y m (i ) + ∂P ( y m ,...) + v ( i ) . (27) Expression of ∂P ( ⋅) is not derived for lack of space.

4. NOISE ESTIMATOR AND STATE OBSERVER Only fundamentals of the theory are presented, missing error convergence and eigenvalue design. Error convergence can be proven in presence of modeling errors following the error loop concept in Canuto (2007a). Eigenvalue design may exploit highfrequency harmonic analysis as in Canuto (2007b, 2007c). 4.1 State prediction

P (γ ) = γ 4 + c3γ 3 + c2γ 2 + c1γ + c0 c3 = lov + α i2

(

)(

)

(

)

c2 = l1v + α i2 + l0 r 1 − α i2 / 4 + α i2 α i2 / 2 + lov / 2 ,

(

) / 2 + (1 − α

(32)

c1 = 1 − α i2 / 4 l1r + α i2l1v / 2 + l2v c0 = α l

2 i 2v

2 i

)

/ 4 l2 r

where the polynomial coefficients, depending on the closed-loop eigenvalues λh , h = 1,… , n, are constant and show multiple solutions to exist. To find a unique solution, let us simplify (32) by setting α i2 = 0 , which leads to c3 lov , c2 l1v + l0 r (33) c1 l2 v + l1r , c0 l2 r and implies four alternatives to exist if useless coefficients are set to zero as in Table 3 . Table 3

Case l1v l2v 0 0 0 1 0 X 2 X 0 3 X X X means not zero

Noise estimator alternatives.

l0r X X 0 0

l1r X 0 X 0

Higher noise Intermediate Intermediate Lower noise

State prediction is strictly related to noise design ending into equation (20). The model error e = y − y m has to real-time estimate noise w . Different than classical observers where driving noise is acting on each state variable, EM noise channels are respected implying no driving noise to directly perturb the CoM rate v as discussed in Section 3.3. As explained in Canuto (2007a), this constraint might require a dynamic noise estimator when (28) n = nc + nd > nw × n y ,

Table 3 and (33) show (i) rate measurement being capable of estimating all noise components as in case 3, (ii) position measurement being necessary to guarantee closed-loop stability as implied by last equation in (33). As position measurement is much noisier than rate as Table 3 shows, case 3 may be selected as it favours Kalman filter guideline of reducing to a minimum measurement noise passing through feedback channels. Formal proof can be developed.

which is not the present case, as position and rate measurements are available from GPS. A static noise estimator applies and holds w ( i ) = L ( i ) e (i ), e (i ) = y ( i ) − yˆ m ( i ) , (29)

4.2 Predictor-corrector

where (i) bar and hat account for noise estimator inaccuracies, (ii) the gain matrix L(i ) , sized nw × n y = 6 > n = 4 is time-varying to force closedloop dynamics (EM+noise estimator) to be LTI. As the matrix L(i ) is oversized with respect to EM dynamics, some entries could be set to zero, while respecting stability. Denote entries as follows ⎡l0 r l0 v ⎤ L(i ) = ⎢⎢ l1r l1v ⎥⎥ (i ) , (30) ⎢⎣l2 r l2 v ⎥⎦

Actually prediction equation (31) can only be implemented with L(i ) = 0 because of to lower measurement rate. This implies conversion to predictor-corrector scheme. Denote with t j = jTg measurement times, such that ti j = t j for i j = ng j . At each time t j the prediction correction is invoked ⎡ Kc ( j ) ⎤ ⎡ xm ⎤ ⎡ xˆ m ⎤ ⎢ z ⎥ ( j ) = ⎢ zˆ ⎥ (i j ) + ⎢ K j ⎥ e ( i j ) , (34) ⎣ m⎦ ⎣ m⎦ ⎣ d ( )⎦ where the variable gains depend on the predictor gain L ( i j ) through the following equation

⎡ K c ( j ) ⎤ ⎡ Ac ( i j ) H c ⎤ ⎡ Gc ⎤ ⎥ ⎢ ⎥ L(i j ) . (35) ⎢ ⎥=⎢ Ad ⎦⎥ ⎣Gd ⎦ ⎣ K d ( j ) ⎦ ⎣⎢ 0 and compute the characteristic polynomial P (γ ) of At each step i the predictor (31) is invoked with the state observer L(i ) = 0 and initial conditions ⎡ Ac ( i ) − Gc L ( i ) H c ⎤ ⎡ xˆ m ⎤ ⎡ xˆ m ⎤ xˆ m (i j ) = xm ( j ), zˆ m (i j ) = zm ( j ) . (36) ⎥ ⎢ ˆ ⎥ (i ) + ⎢ zˆ ⎥ ( i + 1) = ⎢ −G L i z A d ( ) d ⎦⎣ m⎦ ⎣ m⎦ Computing burden due to variable gains is reduced ⎣ .(31) by decomposing gains into steady and variable parts. ⎡G ⎤ + ⎢ c ⎥ L(i )y (i ), yˆ m (i ) = xˆ m ( i ) ⎣Gd ⎦ 5. SIMULATED RESULTS Straightforward computation, upon definition of γ = 1 − λ , λ denoting a generic eigenvalue, yields −1

Table 4

Axis x y z

LORF error RMS [μrad ] .

Total (target) 0.48 (200) 4.0 (200) 1.6 (200)

MBW 0.08 0.8 0.04

HF 0.001 0.03 0.03

Total RMS in Table 4 is largely lower than target in Table 2 , showing filter efficiency. Fig. 3 shows time history of the along-track error: the MBW component is much lower than the whole error as expected, being progressively attenuated by the filter narrow BW.

Fig. 4 shows the observer BW to approach 1 mHz, close to orbit frequency fO 0.2 mHz , because of measurement errors. Only accurate orbit and disturbance dynamics can cope with such a narrow BW: Fig. 5 overlaps the pitch attitude with the corresponding LORF error. Two mission phases are shown: in the former one ending at about 80 ks, attitude is just measured by star trackers affected by noise and bias; in the latter one, data fusion between payload accelerometers and star trackers allow cancelling high frequency noise (bias of course remains). 100 Attitude - Pitch (y) LORF error

50 Angle [microrad]

Simulated results concern GOCE mission during 42 hours (150 ks) corresponding to about 30 orbits. Table 4 shows the a-posteriori error statistics (RMS) for the different LORF axes: total RMS is shown together with MBW and higher frequency components. Being a real-time estimate, low frequency dominates due to relevant components of the measurement errors which are integrated during the filter time constants, as harmonic analysis, not treated here, would have predicted.

0 -50 -100 -150

2 LORF x error data2

-200 0

1.5

5

10 Time [ks]

Angular position [μrad]

1

Fig. 5.

Residual attitude and LORF error.

6. CONCLUSIONS

0

-1

-1.5 0

50

100

150

Time [ks]

Fig. 3. LORF along-track error and MBW component (narrow strip).

Considerations from Fig. 3 are confirmed by the error PSD shown in Fig. 4, having a large margin with respect to target bound. Out-of-plane and radial (lateral) components overlap, while along-track error looks much lower. The reason is due to residual non gravitational accelerations, ideally zero due to dragfree control, but actually non zero at lower frequencies along the lateral orbit movements for propellant saving (Canuto 2007b, 2007c). 10

5

LORF error x LORF error y LORF error z Target bound [μrad/√Hz

7

0.5

-0.5

10

0

10

Fig. 4.

15 x 10

-6

10

-4

10 Frequency [Hz]

LORF error PSD.

-2

10

0

A real-time observer for LORF estimation has been presented designed as a predictor-corrector around the EM of orbit dynamics made by controllable and disturbance dynamics. Simulated results are compatible with GOCE target bounds with a good margin. REFERENCES Bertotti B. and P. Farinella (1990) Physics of the Earth and the Solar System. Kluwer Academic Pu.; Dordrecht. Canuto E., P. Martella and G. Sechi (2003) “Attitude and drag control: an application to the GOCE satellite”. Space Science Reviews. 108 (1-2): 357-366. Canuto E. and A. Rolino (2004) “Multi-input digital frequency stabilization of monolithic lasers”, Automatica, 40 (12), 2139-2147. Canuto E. (2007a) “Embedded Model Control: outline of the theory”, ISA Trans. 46 (3), to appear. Canuto E. (2007b) “Drag-free and attitude control for the GOCE satellite”, submitted to Automatica. Canuto E. (2007c) “Drag-free control of the GOCE satellite: noise and observer design”, submitted to IEEE Trans. on Control Systems Technology. Kaplan E.D. ed. (1996) Understanding GPS: Principles and Applications. Artech House Pu., Boston..