Logiciel de restitution de trajectoire de satellite-application au cas de satellites d'observation

Acta Astroaautica Vol. 10, No. 3, pp. 117-124, 1983 Printed in Great Britain.

0094-5765/83 $3.00+ .00 Pergamon Press Ltd.

LOGICIEL DE RESTITUTION DE TRAJECTOIRE DE SATELLITE-APPLICATION AU CAS DE SATELLITES D'OBSERVATIONt M. GRIMAROet A. M. MArNGUY Office National d'Etudes et de Recherches A6rospatiales(ONERA), 92320, CMtillon, France (Received 23 December 1982) R6sum6---Cet article pr6sente un logiciel de trajectographie de satellite bas6 sur un algorithme de filtrage de Kalman et de lissage de Rauch. Une adaptation originale de cet algorithme permet d'obtenir des estimations statistiques a posteriori des bruits d'6volution et d'observation et de leur 6cart type. L'analyse des r6sultats concemant ces r6sidus permet de r6gler de fafon optimale les niveaux a priori des 6carts types sur les bruits et 6ventuellement de d6tecter des erreurs de mod61es. Une telle utilisation est pr6sent6e, pour des satellites d'observation, en simulationet en traitement de mesures r6elles (satellite TIROS N). Alntraet--This paper presents an orbit determination program based on two algorithms: Kaiman filtering and Ranch smoothing.Novel adaptation of these algorithms lets obtain a posteriori statistical estimations of evolution and measurement noises and of their standard deviation. Analysis of results about these noise estimations lets regulate optimally the a priori levels of noise standard deviations and detect eventually model errors. Such use is presented for Earth observation satellites, at first with simulated data, then with real data (TIROS N satellite). 1. INTRODUCTION L'un des probl~mes auquel se heurte l'utilisation des m6thodes du type filtrage de Kalman, pour restituer des trajectoires de satellites, est la mod61isation des bruits charg6s de repr6senter les incertitudes sur les 6quations du mouvement du satellite. Or le module utilis6 pour les bruits et le niveau a priori de leur 6cart type conditionnent la pr6cision de localisation du satellite en i'absence de mesures, c'est-~-dire entre les p6riodes de visibilit6 au-dessus des stations. Cela a une importance primordiale lorsque le nombre de stations est r6duit (en pr6sence d'une settle station, la dur6e de la pr6diction peut atteindre 12 h). Pour traiter des probl6mes analogues (trajectoires d'avions et de missiles, exp6rience spatiale BIRAMIS) un logiciel de traitement de mesures a 6t6 mis au point/~ I'ONERA. I1 utilise des m6thodes 6prouv6es, le •trage de Kalman pour le traitement en temps r6el et le lissage de Rauch pour le traitement en temps diff6r6. Une extension particuli~rement utile de ces techniques permet, au lissage, d'obtenir des estimations a posteriori des bruits d'6volution et d'observation et de leur 6cart type. L'analyse de ces valeurs permet en effet de v6rifier la bonne coh6rence des hypotheses initiales. Des erreurs de mod61e peuvent ~tre 6ventuellement d6tect6es grgce ~t l'6volution de ces r6sidus au cours du temps. Cet article pr6sente l'utilisation de ce logiciel pour la restitution des trajectoires de satellites. Les m6thodes utilis6es sont rappel6es et des applications plus particufi6res aux satellites d'observation sont ensuite donn6es: tout d'abord des simulations montrant la mise au point d'un module simplifi6 de train6e atmosph6rique; puis des traitements de mesures r6elles sur le satellite

tPaper presented at the 33rd Congress of the International Astronautical Federation, Paris, France, 27 September-2 October 1982.

TIROS N. I1 est montr6 comment les r6sultats concernant les r6sidus d'6volution et d'observation permettent de r6gler de fa~on optimale les niveaux a priori des 6carts types des bruits et d'am61iorer les modules. 2. ALGORITHMESEl" LOGICIELDETRAJECTOGRAPHIE 2.1 Filtrage de Kalman L'algorithme utilis6 applique le filtrage de Kalman classique[1] /t des syst~mes non lin6aires, /t 6volution continue et observations discr/~tes [4]. L'6volution du syst~me dynamique est r6gie par l'6quation vectorielle: ::(t) = f(x(t), t) + G(t)v(t),

(l)

oil: x(t) est le vecteur d'6tat, de dimension n; [ e s t le vecteur des d6riv6es non bruit6es de l'6tat, de dimension n; v(t) est un vecteur de bruit, de dimension m; G(t) est la matrice de pond6ration du bruit d'6volution, de dimension (n x m). Le bruit d'6volution v(t) est un bruit blanc, continu, gaussien, de moyenne nulle et de matrice de covariance Q(t) telle que: E[v(t)v(s) r] = Q(t)8(t - s) (~(t - s) = fonction de Dirac) L'6tat est estim6 ~ l'aide de mesures, rang6es ~ un instant tk dans un vecteur Yk, de dimension p, obtenu partir de i'6tat ~ l'aide des 6quations d'observation: Yk = h(xk, tk)~+ wk,

(2)

of1 Yk est ie vecteur de mesures, de dimension p; h est le vecteur de l'observation non bruit6e, de dimension p; wk est un vecteur de bruit, de dimension p. Le bruit d'observation wk est un bruit blanc, discret, gaussien, de moyenne nulle et de matrice covariance Rk 117

118

M. GRIMARDet A. M. MAINGUY

• une r6estimation est faite s'il y a une mesure ~ tk+,. L'~tat estim~ ~+11~+~ est donn6 par:

teUe que: E(w~w~f] = Rk'Sk~ (Sk~= symbole de Kronecker).

)~k+llk+ I = £k+llk +

Une technique de lin6arisation 6tendue permet d'obtenir des 6quations lin6aires discr6tes d'6volution et d'observation. Entre les instants t~ et t~+,, le mod61e d'6volution non bruit6 x(t) = l ( i ( t ) , t),

t ~ [tu, t~+d,

(3)

{

- ~ (t, s) = F(t)~(t, s),

(4)

¢(s, s) = Id, F(t) = ~ (£(t), t).

(5)

Le module d'6volution lin6aire discret a pour forme: x~÷, = F~xk + Vk,

(6)

Pk+tlk+~ = Pk+Hk - Kk+tHk*~Pk+tlk.

h(-~+llk, t~+,) Kk+~ est foumi par:

Zk+l = Yk+l --

Kk+I

(11)

est l'innovation et le gain

T T -1 = Pk+tI~Hk+,(R~+, + Hk+,Pk+,lkHk+,) .

(12) 2.2 Lissage de Rauch Le lissage permet de corriger les estimations fournies par le filtrage en tenant compte, sur une p6riode de temps donn6e, non seulement des observations ant6rieures mais aussi des observations ult6rieures ~ l'instant tk. Cette m6thode en temps diff6r6 fair appel ~ des 6quations r6trogrades dont la forme g6n6rale est: XalN = :¢~1~- P~Ik FkTP k l + l l k ( ' ~ / + l l k - - " ~ k + l [ N ) , T -I -I PklN = Pklk -- PklkFk P k+~lk(Pk+~lk -- P k + ~ m ) P k+HkFgP~lk.

(13)

avec F~ = ¢P(tk+l, t~)

v~ = f[+~ ¢(tk+l. $)G(s)v(s) ds. Le bruit Vk est un bruit blanc, discret, gaussien, de moyenne nulle et de matrice de covariance Qk:

~

tk+l

Q~ =

(10)

La matrice de covariance de l'erreur d'estimation Pk+~l~+l a pour expression:

du filtre est lin6aris6 autour de '2(h,)=.~klk=E[xdy~ . . . . . yk], puis discr6tis6, en utilisant la matrice de transition q~(t, s), d6finie par:

K~+,z~+,

¢(t~+,, s)G(s)Q(s)Gr(s)~bT(tk+,, s) ds.

Cette formulation prgsente l'inconvgnient d'utiliser I'inverse de la matrice Pk+llk (dimension n x n), dont le calcul est cofiteux et imprgcis pour des dimensions d'6tat n 61ev6es. La formulation de Rauch[2-3], 6tendue par Aumasson[4], permet d'6viter cette inversion en introduisant un vecteur t0klN[2-3] et une matrice TklN[4] auxiliaires, d6finis par:

k

~ tokls = P k +1l l k ( X ^k + l l N

(7)

t Tkm Le module d'observation est lin6aris6, ~ l'instant tk+,, autour de la meilleure 6valuation de xk+,, ~ savoir ~Ck+llk = E[xk+dyl . . . . . y~]: Yk+l = Hk+,Xk+l + Wk+,,

p~t+llk -

-- 3 ( k + l l k ) ,

P -1 k+HkPk+IINP -1 k+Hk.

(14)

t0klN et Tkm sont calculgs par r6currence, en partant de oJr~m= 0 et TwIN =0, ~ l'aide des 6quations r6trogrades:

(8) [ T k - l m = F- k T k l s F- kT + H k TR -k- I Hk

avec:

(15)

avec:

Oh . tk+3. Ha+I = tg-~(Xk+llk'

Fk = [ I d - H r k ( R k + H k P k l k - l H kr)- i I'IkPklk-I]Fk, T

La construction r6cursive de la meiUeure estimation de l'6tat xk+t h l'instant tk+t, compte tenu des mesures ant6rieures A tk+l et de la mesure 6ventuelle fi tk+l, et de sa matrice de covariance P = E[(£ - x)(£ - x)T], proc~de du sch6ma suivant: • une pr6diction ~k+llk de l'6tat h l'instant tk+l est obtenue en int6grant l'6quation d'6volution non bruit6e (3), h partir de £(tk)= iklk. La matfice de covariance pr6dite Pk+llk est donn6e par: Pk+.k =

FkPktkFkT + Qk.

(9)

Rk = Rk + HkPklk-,H ~.

L'estimation liss6e ~kl~ et la matrice de covariance PklN sont alors obtenues par: { -~kIN = Xklk-I + Pklk-ltOk-HN PklN = Pklk-, - Pklk-t Tk-IINPklk-l. •

%.

,

(16)

Cette formulation utihse l'mverse de la matfice/~k, qui a d6j~ 6t6 calcul6e au filtrage pour l'obtention de la matrice de gain Kk et qui n'est que de dimension p × p.

119

Logiciel de restitution de trajectoire de satellite application au cas de satellites d'observation 2.3 Extensions du lissage: r~sidus d'~volution et d' observation Le lissage, sous la form¢ pr6sent6e ci-dessus, permet d'obtenir des r6sultats statistiques int6ressants: les estimations a posteriori des bruits d'6volution et d'observation, appel6es r6sidus[4]. L'estimation a posteriori du bruit d'6volution discr6tis6 VklN-- E[vk/y~ . . . . . yN] est donn6e par: VklN = ~k+,N - Fk~klN.

(17)

On d6montre que: 1)kin = Qkt0klN

E[ VklNVkrlr~]= QkTkmQk"

(18)

v~ls, ayant m~.me dimension n que le vecteur d'6tat, est di~cilement comparable au bruit d'6volution continu v(t), de dimension m, qui est la v6ritable entr6e du filtre. I1 est donc int6ressant de d6finir un bruit d'6volution "pseudo-continu" ~, de m6me dimension que v(t). On utilise la pseudo-inverse G ~ = ( G kTG k ) - , Okr de la matrice de pond6ration G ~ = G ( h ) , introduite darts l'6quation (1): ~5~= 1

G~'vk

moyenne nulle et coh6rents avec les 6carts types calcul6s. 2.4 Presentation du Iogiciel de trajectographie de satellite: param$tres d'~tat et modiies d'&olution Le logiciel d'orbitographie de I'ONERA utilise les m6thodes pr6sent6es ci-dessus pour estimer la trajectoire d'un satellite h partir de mesures d61ivr6es par une ou plusieurs stations de poursuite. L'6volution du satellite est repr6sent6e ~ l'aide d'un syst6me de param6tres "orbitaux vectoriels" (Fig. 1): a = demi-~and axe, e = vecteur excentricit6, dirig6 du centre de la Terre vers le p6rig6e, de module 6gal /~ rexcentricit6 e, i = vecteur inclinaison, unitaire sulvant le moment cin6tique, ~ -- 1"1+ co + v = angle bris6 rep6rant ia position du satellite le long de l'orbite, par rapport fi un axe inertiel (fl, ascension droite du noeud ascendant; ~o, argument du p6rig6e; v, anomalie vraie). Remarque: Ces variables sont au nombre de 8, mais les relations entre les vecteurs e et i (e.l = 0 et I~il]-- 1) les ram~nent en f a i t h 6 variables ind6pendantes. A la diff6rence des param6tres orbitaux classiques, les param~tres orbitaux vectoriels sont d6finis pour tousles types d'orbites, y compris les orbites circulaires (e -- 0) et 6quatoriales (i = z).

(19)

zI P~lenordmoye~/ Satellite // // ..dP~rig(~e / i I i

0k est un bruit blanc, discret, de matrice de covariance: G~'T -~ Q(tk).

E[~kfik r ] = G~ ~

Le r6sidu d'6volution pseudo-continu f)klN e t matrice de covariance sont donn6s par:

1

T

(20) sa

Point

G k Q~okl., (21)

¥S0

Fig. 1. Param/~tres orbitaux vectoriels.

E[VktNV kl~] = tk+t -- t--------~ G k QkYk~QkG k . L'estimation a posteriori du bruit d'observation Wk,N = E[wdy~ . . . . . ys] est fournie par: Wkls = Yk- Hk~klN.

(22)

On d6montre que:

WklN= Zk --HkPklk-I(Ok-IIN, E[wkmw kTl~]

=

Rk - HkPklNHkT.

(23)

Le calcul de ces r6sidus est peu cofiteux car il utilise des termes calcul6s au filtrage (Zk, Hk, P~l~-', Rk, Qk) et au lissage (
Cette repr6sentation a l'avantage de ne faire intervenir que les acc616rations perturbatrices (autres que racc616ration k6pl6rienne) clans les 6quations d'6volution, ce qul,/L m6thodes d'int6gration identiques, assure une meilleure pr6cision du calcul que pour des param~tres du type r, V oi~ les acc616rations perturbatrices sont explicitement ajout6es /L l'acc616ration k6pl6rienne. Les 6quations d'6volution correspondant /L ces variables d'6tat a, e, i, ~ sont: d = 2a: V.%,

1 = ~'0tp ^ l l + v ^ (r ^ ~/.)). i = ~1 (r ^ V D - i.(r h'Ye) 1, h

~=~-g

ri~sin~+i, cos~,

l+iz

(24) ,,

t~/r,.

120

M. GRIMARDet

avec: r, V: rayon vecteur centre Terre-satellite et vitesse absolue du satellite (calcul6s ~ partir de a, e, i, , ) , r = Ilrll

h = r ^ V: moment cin6tique, h = Ilhll ~/,: acc616ration perturbatrice /~: constante de gravitation terrestre i , i~, i~: coordonnres de i dans O~y~. L'accrlrration perturbatrice 3'~ peut ~tre 6crite sous la forme: ~'o = 3%, + ~r + 'YP~ + "gLS (25) avec: %0, accrlrration perturbatrice de gravitation terrestre, ~'T, accrlrration de trainre atmosphrrique, 3'~, accrlrration de pression de radiation, 3'LS, accrlrration perturbatrice de gravitation lunisolaire. Les travaux prrsentrs ici ont plus particuli&ement port6 sur la modrlisation des erreurs faites lors du calcul des accrlrrations de gravitation terrestre et de trainre atmosphrrique. Ces erreurs ont en effet une influence prrpondrrante sur la prrcision de localisation d'un satellite h basse altitude. Ces incertitudes sur le modrle d'rvolution sont introduites sous la forme de "bruits colorrs" qui sont des processus markoviens du premier ordre. Ceci prrsente l'intrrrt d'avoir des 6carts types stationnaires rrglables (fonctions des niveaux des bruits blancs et des temps de corrrlation) pour ces termes correcteurs. Ainsi, l'accrlrration perturbatrice de gravitation terrestre a pour expression: Tpo, = grad U ~ + A%o, .

1

a~/tpOt = -- 7 a ~ / p Ot + ~pot Tpot

(26)

avec: UM, potentiel perturbateur de gravitation terrestre donn6 par le module M, choisi compte tenu de l'application envisa#e, A~/~ot, vecteur d'erreurs de potentiel, ~'~o,, temps de corrrlation, g,,ot, vecteur de bruits blancs continus gaussiens. De fagon/i d6coupler les probl~mes, les incertitudes sur la trafnre atmosphrrique ont 6t6 introduites de faqon ind6pendante en agissant sur la masse volumique:

(27)

avec:

Co, S, m, coefficient de trainre, maitre couple et masse du satellite, p, masse volumique, 0,, masse volumique calculre,

A. M. MAINGUY Ap, terme correcteur de masse volumique, rp, temps de corrrlation, Gp, bruit blanc continu gaussien. Les termes correcteurs A'ypot et Ap sont des variables d'rtat dont la variance est quasi-stationnaire. Les niveaux respectifs de A~pot (erreurs de gravitation terrestre) et de A~/T induite par Ap (erreur de trafnre) drpendent de l'altitude. Des param~tres supplgmentaires peuvent 6tre introduits dans le vecteur d'gtat: ce sont des constantes real connues intervenant dans le modrle d'rvolution et dans le module d'observation (position de station, frgquence 6raise . . . . ). Leur 6quation d'rvolution est bien stir du type ~ = 0. Le logiciel de trajectographie de I'ONERA permet donc d'estimer, ~ l'aide des mrthodes de filtrage de Kalman et de lissage de Rauch, un vecteur d'rtat compos6 comme suit: • paramrtres de trajectoire: a, e, i, tp • termes correcteurs de l'rvolution: A%ot, Ap • biais d'rvolution et d'observation. 3. RESULTATS DE SIMULATIONS

3.1 Presentation Des 6tudes de simulation ont 6t6 effectures pour un satellite hrliosynchrone situ6 sur une orbite quasi-circulaire d'altitude moyenne 500km. Elles ont permis d'examiner l'influence des erreurs de trafnre atmosphrrique, importantes ~t cette altitude, sur la prrcision de Iocalisation et de srlectionner un module simple de masse volumique en restitution d'orbite (eqn (27)). Des mesures de distance, fournies par une seule station de poursuite (Toulouse) avec une prrcision de 15 m, sont simulres ~ partir d'une trajectoire intrgrre numrriquement. Les "prriodes de visibilitY" sont composres de 2 (ou 3) passages consrcutifs au-dessus de la station/~ environ 1 h 30 mm d'intervalle, chaque # r i o d e de visibilit6 6tant srparre de la suivante par une "#riode de prrdiction", sans mesures, d'environ 12 h. Les forces perturbatrices utilisres pour cette intrgration sont rrduites au potentiel de gravitation terrestre et ~ la trafnre atmospMrique. Un modrle d'atmosph~re complexe et rraliste (DTM [5]) permet de calculer la masse volumique atmosphrrique p en prenant en compte les effets actuellement connus agissant sur l'atmosph~re (flux solaire, indice #omagnrtique, effets diurne et annuel...). A partir de ces mesures, une trajectoire du satellite est estimre en utilisant le logiciel prrsent6 ci-dessus. Les forces perturbatrices prises en compte dans le module d'rvolution (rqn (1) et (24)) sont la gravitation terrestre et la trainre. Cette 6tude en simulation ayant pour but d'examiner i'influence des erreurs de traln6e seules sur la prrcision de localisation, le m~me potentiel de gravitation terrestre que pour la simulation des mesures est employr. Par contre, il est fait appei ~t un modrle simplifi6 bruit6 de masse volumique (rqn (27)) pour le calcul de la trainre atmosphrdque. En effet, la drtermination a priori de la masse volumique h l'aide du modrle complexe (et cotiteux) DTM [5] est hasardeuse, compte tenu

Logiciel de restitution de trajectoire de satellite application au cas de satellites d'observation de la mauvaise connaissance a priori des donntes d'entrde. Du far de la quasi-circularit6 de rorbite, il n'est pas ntcessaire d'introduire un module exponentiel en fonction de raltitude. Trois formulations de la masse volumique simplifite p, sont successivement prtsenttes: --masse volumique constante: p~ = p,,, /~,, = 0,

121

~PxlO t~(kg/mJ) 2,S2

-

1,5-. I

TM ffi P s i m u l a tion

I

,

B

l

;=Pro

-

(28)

. . . . . .

Ps~,Psj

~,

w+ v

~o

--masse volumique avec effet diurne simplifi6 autour de la masse volumique moyenne, constante, p,~:

H I =Igh

I~o

2'70

(.) ,_

3'6o

H I =2h

H I =Igh

Fig. 2. Modtles de masse volumique:6volution sur une orbite. ~r - H~ m~)]J Ps2--Pm 1 + k cos~(H~

=0

(29) d'tvolution continu est de dimension 1, et pour ie modtle P,3, il est de dimension 2. Les rtsultats prtsentts montrent comment ranalyse des rtsidus d'tvolution permet de rtgler les niveaux a priori des 6carts types des bruits et d'amtliorer les modtles eux-mtmes.

ob H, est rheure locale, et HIm~ = 14 h. --masse volumique avcc effet diurne simplifi6 autour d'une masse volumique moyenne, bruitte:

Ap,. = - ~

1

Ap,, + ~,,,.

3.2 R~sultats et analyse des r~sidus d'~volution: amelioration des modules de masse volumique Quatre types de traitement, dtfinis darts le Tableau 1, sont prtsentts. Les prtcisions sur la position longitudinale, h ia fin d'une ptdode de prtdiction de 12h, faisant suite h deux jours de traitement des mesures, leur sont assocites dans le Tableau 2. C'est en effet cette composante qui est la plus sensible aux erreurs de trainte. La comparaison des cas 1 et 2 montre les possibilitds de rtajustement des niveaux des 6carts types des bruits a priori. Les Figs. 3 (cas 1) et 4 (cas 2), reprtsentent le rtsidu d'tvolution pour le bruit ~o et son 6cart type. Rappelons que, d'aprts les hypothbses de ralgorithme, ce rdsidu doit 8tre de moyenne nuUe et les dispersions autour de cette moyenne cohtrentes avec rtcart type (Section 2.3). La Fig. 3 fait apparaitre un choix trop optimiste du niveau a priori de rdcart type du bruit, alors que celui-ci est plus judicieux stir la Fig. 4. Ceci est confirm6 par les rtsuitats sur la prtcision d'estimation (Tableau 2). Dans le premier cas, le niveau trop faible du bruit conduit effectivement ~, un 6cart type 6galement peu 61ev6 sur la composante longitudinale, qui n'est pas en accord avec rtcart constatt; Dans le deuxi~me cas, par contre, rensemble est plus cohtrent.

(30)

Les 6volutions sur une orbite des modtles p,~, ps2, ps3 sont compartes h la masse volumique p de ia simulation sur la Fig. 2. La masse volumique moyenne Pm est calculde a priori, en fonction de l'altitude et des activitds solaire et gtomagnttique prtvisibles, et ajustte par le filtre h l'aide des mesures, en tant que param~tre d'ttat. Dans tous les cas, le temps de corrdlation ~ du bruit color6 Ap (tqn (27)) est pris 6gal au quart de la ptriode orbitale, ce qui doit permettre de compenser les erreurs it court terme comme l'absence ou la mauvaise modtlisation de reffet diurne (Fig. 2). Par contre, pour le module P~3, le temps de corrtlation Tp~ du bruit color6 Ap,, (tqn (30)) est pris 6gal ~ 6h, ce terme dcvant prendre en compte des erreurs sur la masse volumique moyenne, dues essentieilement aux variations d'indices gdomagnttiques et de flux solaire journalier. Le vecteur d'ttat identifi6 par l'algorithme de restitution est composd des six variables de trajectoire et des termes correcteurs de trainte p,., Ap dans le cas des modules p,, et P~2, puis Ap,, pour le troisi~me modble. Pour les deux premiers modtles, le vecteur de bruit

Tableau 1. Simulationavec erreurs de trainte--Description des cas traitts (1' "unitt" de bruit blanc 6 est le kg/m3/sI/2)

Cas

Mod61e de p,

Bruit a priori sur Ap (¢~ = Torbl4) Ecart type du Niveau stabruit/~p tionnaire de
p~t: constante

1.8 x 10-~4

p~t: constante

9 x 10-t4

10-13 2.4 x 10-12

Bruit a priori sur Ap,. (%, = 6 h) Ecart type du Niveau stationbruit/~p, naire de %p (kg/m3) "

4.8 x

P,2: effet diurne p~3: effet

1.8x 10-~4

4 . 8 x 10-13

diurne + effet moyen bruit6

1.1 × 10-t'

3 X 10-13

m

9.6x 10-16

10-13

122

M. GRIMARDet A. M. MAINGUY Tableau 2. Simulationavec erreurs de trainre--Prrcision en longitudinal: prrdiction A 12h (6 = valeur estimrevaleur de rrfrrence, en simulation seulement) Cas

Ecart type ~r (m)

1 2 3 4

Ecart (m)

1100 4000 I100 650

( • ' Ek /)N = (ubbl )) (ubbl

lO~s

2-

,. 0

,,"~j';,, /x

......

- bv

r t(h)_ -

-2.

2100 1500 1200 500

Fig. 5. Rrsidu d'rvolution du bruit sur Ap. Cas 3:--modrle p,., (effet diurne simpliflr);----rcarttype a priori du bruit: 1.8x 10-1~ubbl.

'~.'~p) k / N = lO 's

8-

1.6,

6o

A 4 A A

A / /VEEE\

4" 20

I .........

..,-s-,:---,.,--,

tr/,~

,--, t~..-

, (ubbl)

a~

^,

0

;,.y'

L ~ '-, - t Xla l

t(hJ_

-O.t.

Fig. 3. Rrsidu d'rvolution du bruit sur Ap. Cas l:--mod~le p,~ (constante);---rcart type a priori du bruit: 1.8x 10-~ubbl ("unitr" de bruit blanc = ubbl = kg/m~/s~n).

-0,8

#,'" ,;

-'

-/,2

Fig. 6. R~sidu d'~volution du bruit sur Ap. Cas 4:--modUle 9~ (effet diurne simplili6et masse volumique moyerme bruitEe);-6cart type a priori du bruit: 1.1 × 10-a~ubbl. 12$.

40 ¸

-8'

- 12-~ -

/ V v,v:v ;'

, .....

i • *

;, 1 L~

9

•

~l"," i

Ii~

16

Fig. 4. REsidu d'6volution du bruit sur Ap. Cas 2:--modUle P~ (constante);--4cart type a priori du bruit: 9 × l0-~4ubbl. Sur la Fig. 4, l'rvolution du bruit restitu6 est #riodique. I1 existe donc une erreur de module drterministe et une amrlioration de celui-ci doit ~,tre possible. La #riode du phrnom~ne 6tant 6gale h la #riode orbitale, rerreur de module a 6t6 attribure ~ la non prise en compte de l'effet d'heure locale sur la masse volumique. Le module p,2 (rqn (29)) permet de s'en rapprocher avec une formulation peu coi~teuse, mais laisse subsister des erreurs (Fig. 2). Cette amrlioration permet cependant une diminution substantielle du niveau du bruit ~ao (Fig. 5) et entraine une meilleure prrcision de localisation (cas 3 du Tableau 2). II appara~t alors tr~s nettement une drrive sur le rrsidu d'~volution (Fig. 5), qui peut ~tre imputre ~ une variation du niveau de la masse volumique moyenne. Sur de longues prriodes de temps, elle est due aux variations d'activitrs solaire et #omagnrtique. L'utilisation du module p~3 (rqn (30)), prenant en compte cette erreur par l'intermrdiaire d'un bruit sur la masse volumique moyenne pro, permet d'amrliorer les rrsultats (Fig. 6 et cas 4 du Tableau 2) et de diviser par deux la valeur des 6carts et 6carts types sur la composante longitudinale.

Cet exemple montre bien l'intrrrt des rrsidus d'rvolution. L'analyse de la cohrrence des rrsidus avec leur 6cart type permet de rrgler de fa~on optimale les niveaux a priori des 6carts types des bruits. L'rtude de la forme des rrsidus a conduit ~ amrliorer le module d'rvolution lui-m~me, ce qui a entrain6 6videmment une meilleure prrcision de localisation. Enfm, il appara~t que, dans le cas d'un satellite ~ basse altitude, soumis ~ une atmosphere dense, un module simplifi6 de trainre atmosphrrique peut suffire pour obtenir une bonne qualit6 de localisation. 4. T R M T E M E N T DES M E S U R E S R E E L L E S SUR TIRO~; N

Afin de valider ce logiciel de trajectographie de satellite, il 6tait indispensable de traiter des mesures rreUes et de comparer les estimations de trajectoire avec celles obtenues avec un logiciel bien 6prouv6 4.1 PrEsentation Dans le cadre du syst~me ARGOS de localisation par satellites, mis au point par le CNES, des mesures de type Doppler sont rrguli~rement faites sur le satellite TIROS N (hrliosynchrone ~ une altitude de 850 km), permettant d'effectuer la trajectographie de rrfrrence utilisre pour le syst~me. Ces mesures, foumies par le CNES [7], ont 6t6 traltres ~ ralde du logiciel de I'ONERA et les rrsultats obtenus comparrs ~ ceux du CNES, 6tablis avec une mrthode de moindres carrrs. Apr~s analyse, il a sembl6 suliisant, pour la prriode de restitution consid&re (3 jours), de ne prendre en compte que les effets du potentiel de gravitation terrestre et de la trafnre atmosphrrique. Le potentiel terrestre est reprrsent6 par le module GEMIO-A de la NASA[6] et la trainee atmosphrrique par le modrle simplili6 du para-

Logiciel de restitution de trajectoire de satellite application au cas de satellites d'observation graphe pr6c6dent, avec la masse volumique p~2. Les termes correcteurs A%ot pour le potentiel et Ap pour la trafn~e sont utilis~s, avec les bruits blancs ~ot et ~ao. Pour la # d o d e 6tudi6e (7 au 9 mai 1980), dix stations bien r6parties sur le globe terrestre 6taient actives[7]. Les fr6quences Fo~ des balises 6mettrices n'6tant 6gales ~t la valeur nominale #o = 401650 000 Hz qu'~ 10 Hz pros, elles sont introduites dans le vecteur d'6tat en rant que biais d'observation, en plus des six variables de trajectoire et des quatre termes correcteurs des modules d'~volution. La pr6cision des mesures est prise 6gale A 1 Hz, ce qui englobe les erreurs de mesures proprement dites et les erreurs dues A la non-correction des effets de r6fraction ionosph6rique et troposph6dque. Divers r6sultats sont pr6sent6s, concernant la pr6cision de restitution et l'analyse des r6sidus d'observation et d'6volution. 4.2 Precision de restitution avec dix stations La trajectoire estim6e en temps diff6r6 (filtrage et lissage), sur trois jours de mesure, est compar6e avec la trajectoire du CNES, obtenue ~ l'aide d'un traitement par moindres carr6s sur trois jours. En prenant en compte les mesures fournies par les dix stations actives, la pr6cision obtenue peut 8tre ~valu6e l'aide des 6carts # par rapport ~t I'estimation du CNES et des 6carts types tr maximaux atteints sur ies trois jours (Tableau 3). Les 6carts maximaux sont inf6deurs ~ 3tr et Tableau 3. Traitement des mesures r&:lles sur TIROS N avec 10 stations. Ecarts types maximaux et &:arts maximaux par rapport la trajectoire CNES Altitude a ~ . (m)

10

I~lmax

30

Longitudinal 40 100

Transversal 25 60

(m)

les valeurs affich6es sont en bon accord avec la pr6cision de restitution du CNES. Un exemple de l'6volution de ces 6carts et 6carts types est donn6 sur la Fig. 7, pour la composante Iongitudinale, ~ la fin de la troisi~me

journ6e. La qualit~ de la restitution est aussi tr6s bonne en ce qui conceme ridentification des fr6quences des balises. L'~cart avec les fr~quences obtenues par le CNES est inf~rieur ~ 0.18 Hz, pour des ~carts types allant de 0.07 0.14Hz. Par exemple, pour la station de Kourou, la fr~quence estim6e par le Iogiciel ONERA est 401650009.09 Hz avec un ~cart type de 0.14 Hz et l'~cart avec la fr~quence estim~e par le CNES est de 0.09 Hz. 4.3 Analyse des r~sidus d'observation Les r6sidus d'observation concernent les fr6quences Doppler mesur6es. Le bruit de mesure a priori est de 1 Hz. La Fig. 8 pr6sente le r6sidu d'observation, obtenu au lissage, pour un passage au-dessus de la station de tLe CNES a confirm~ que cette station a effectivement un mode de fonctionnement d~gradE.

60. 40-

123

Erreur Iongitudinale (m)

~.

20" 0

A

-20. - 40- 60-

Nes/ures

~.

Fig. 7. Traitement de mesures r6elles: filtrage et lissage sur 3 jours, avec 10 stations--Pr&:ision longitudinale. (. . . . ): &:art type; ( - - ) : &:art (estim&:-valeurCNES).

3 ~ Wk/N (Hz)

3.] Wk/N (Hz)

2-

2-

t.

.---.p--

"":: t ( ~ ) . o

t~

i.." Yo

-I-

. . . . .

4'0

' ,-."

30' *~. . . .

-1

_e_

-2-

4'o

-2~

-3- ~-

Ib

o

-31 (a)

-4-

(b)

Fig. 8. REsidu d'observation pour un passage au-dessus de Kourou. (a) avec toutes les mesures. (b) avec les mesures de site > 10°. Kourou. Le r6sidu de la Fig. 8a a 6t6 obtenu en utilisant toutes les mesures du passage en visibilit6. Il apparaft nettement que ies mesures en d6but et en fin de visibilit6, c'est-~-dire ~ site tr~s bas, sont mauvaises. Physiquement, cela vient de la non correction des erreurs de r6fraction ionosph6rique et troposph6rique. Celles-ci sont en effet maximales pour les mesures ~ site bas, le trajet des ondes 6tant le plus long. Ceci est confirm6 par la Fig. 8b qui montre que, lorsque les mesures de site inf6rieur ~ 10° sont rejet6es, le r6sidu est beaucoup plus coh6rent avec son 6cart type. Ce rejet des mesures ~ site bas, d6cid6 ~ la suite de cette analyse du r6sidu, a permis d'am61iorer sensiblement la pr6cision des estimations de la trajectoire. Malgr6 ce rejet de mesures, le r6sidu d'observation pour la station d'Ancon reste assez 61ev6 et en dehors de son ~cart type (Fig. 9).t Vu le niveau des r6sidus, un 6cart type a priori de 2 Hz serait mieux adapt6 pour cette station. d

Wk/N(Hz)

2. t. 0

-

Ib

2~

YO

t(mn ] 4'0

.

t1

-2,

le

Fig. 9. R6sidu d'observation pour un passage au-dessus de Ancon. Seules les mesures de site > 10° sont prises.

124

M. GRIMARDet A. M. MAINGUY

4.4 Influence du nombre de stations sur les bruits d'Evolution et sur la precision de restitution La valeur des 6carts types, obtenus ~ des instants o/~ aucune mesure n'est faite, est d'autant plus li6e au niveau a priori des 6carts types des bruits d'6volution que l'instant consid6r6 est "loin" d'un instant de mesure. Le r6glage de ces niveaux est donc primordial, surtout si les p6riodes sans mesures sont longues, comme c'est le cas lorsqu'une seule station de poursuite est utilis6e. L'analyse des r6sidus d'6volution permet d'effectuer correctement ces r6glages. Ainsi, pour TIROS N, un r6glage de bruits (vecteur yoo~ et ~a,,)'satisfaisant a 6t6 obtenu, pour la restitution avec dix stations de mesure, comme le montre la Fig. 10, pour la composante radiale ~a~ de %o~: le r6sidu est bien l'int6rieur de son 6cart type (l'unit6 de bruit blanc pour ~a~ est le m/s~/~). On constate d'autre part que le r6sidu et son 6cart type augmentent en dehors des p6riodes de mesures, lh or) les erreurs de module d'6volution ont le plus d'importance. Le m6me r6glage de bruits 6tant utilis6 pour une restitution avec une seule station de poursuite (Toulouse), l'analyse des r6sidus a montr6 que les niveaux a priori des 6carts types des bruits 6taient alors insuilisants: pour la m6me p6riode que pour la Fig. 10, le r6sidu sur ~a~ pour une seule station est donn6 sur la Fig. l l. Le r6sidu est nettement en dehors de l'6cart type. Pour avoir des r6sidus coh6rents avec leur 6cart type, dans ce cas de la restitution avec une seule station, il a donc fallu augmenter les niveaux a priori des 6carts types des bruits, par rapport au cas de la restitution avec dix stations.

6¸

'~&lr ) k/N x10l (ubb2)

2.

i!

' ?--/ Mesures

5,+' Nesures

Fig. 10. Traitement de mesures rrelles avec 10 stations. Rrsidu d'rvolution du bruit sur la composante radiale de A'Yoot,6cart type a priori du bruit: 7.4 × 10 ~ubb2 ("unitr" de bruit blanc = ubb2 = tulsa/").

l(~"~lr) 2~

-

k/N*'

I0"~

(ubb2)

3-i

I Fig. 11. Traitement de mesures r~elles avec I station. R~sidu d'~volution du bruit sur la com~osante radiale de ~pot, ~cart type a priori du bruit: 7.4 x 10- ubb2 (p~riode sans mesures).

Tableau 4. Traitement des mesures rrelles TIROS N avec une station. Ecarts types maximaux et 6carts maximaux par rapport la trajectoire CNES Altitude C~max(m) I~1o,~ (m)

50 160

Longitudinal 250 500

Transversal 110

Ce'ci s'explique aisrment: un module d'rvolution n'est jamais parfait; la prise en compte de mesures aprrs une prriode de prrdiction permet de drtecter l'amplitude des erreurs dues ~ celui-ci. Pour un module donnr, cette amplitude est d'autant plus grande que la prriode de prrdiction est plus longue. Rappeions que ces erreurs sont prises en compte par le bruit d'rvolution. D'autre part, la disparition d'un certain nombre de mesures bien rrparties conduit a une moins bonne observabilit6 des param~tres. Cette augmentation des bruits, pour la restitution avec une station, ainsi que le moins grand nombre de mesures entrainent une drgradation tr~s nette de la prrcision par rapport au traitement avec dix stations, comme le montre la comparaison des Tableaux 3 et 4. 5. CONCLUSION Rappelons que la m~thode de filtrage de Kalman associre au lissage de Rauch permet de disposer en permanence (avec ou sans mesures) des 6carts types des param~tres de localisation. Les algorithmes d'rvaluation des bruits d'observation et surtout des bruits d'rvolution permettent ensuite de contrrler la validit6 de la prrcision de localisation en fonction des diverses hypotheses (modrles d'rvolution plus ou moins complexes, modrles d'observation, etc...). Les rrsultats du traitement des mesures rrelles sur TIROS N prouvent le bon fonctionnement du logiciel de trajectographie de satellite mis au point h I'ONERA. Ce logiciel est actuellement un outil de recherche mais peut 6ventuellement servir de base ~ l'61aboration d'un logiciel oprrationnel employant les m~mes mrthodes de traitement. REFERENCES

1. R. E. Kalman, A new approach to linear filtering and prediction problems. J. Basic Engng, Trans. ASME 82, 35--45 (1960). 2. H. E. Rauch, F. Tung and C. T. Striebel, Maximum likelihood estimates of linear dynamic systems. AIAA 3. 3, 1 445-1 450 0965). 3. H. E. Rauch, Linear smoothing techniques. AGARD ograph 139: Theory and Applications of Kalman Filtering, Chap. 8. 4. C. Aumasson and T. Cotillon, Presentation d'un code num~rique de filtrage et de lissage de Kalman-Rauch pour des syst~mes non lin~aires, ~ ~volution continue et ~ observation discrete. Note Technique ONERA No. 1981-5 (1981). 5. F. Barlier, C. Berger, J. L. Falin, G. Kockarts et G. Thuillier, A thermospheric model based on a satellite drag data. Ann. G~ophysique 34(1). 9-24 (1978). 6. F. J. Lerch, S. M. Klosko et R. E. Laubscher, Gravity model improvement using GEOS3 (GEM9 and 10). Geodynamics Branch-Goddard Space Flight Center (1977). 7. A. Cariou et M. Winterholer, ARGOS-Doun~es disponibles sur bandes magn~tiques, concernant les mesures d'orbitographic, les param~tres d'orbitographie, les rrsultats d'orbitographic, les r~sultats de localisation. Rapport CNESCTIEMT/MTIMSI263 (1979).