0% ont trouvé ce document utile (0 vote)
22 vues31 pages

Introduction au filtrage de Kalman

Le document présente le filtrage de Kalman, une méthode d'estimation de paramètres utilisée pour la fusion de données provenant de divers capteurs, essentielle dans les systèmes de navigation modernes. Il vise à illustrer le principe du filtrage de Kalman à travers des exemples simples, en établissant les équations nécessaires et en justifiant leur utilisation par la méthode des moindres carrés. Le texte aborde également des applications pratiques dans le domaine du positionnement dynamique, permettant au lecteur de développer ses propres applications.

Transféré par

Benjamin Jourdain
Copyright
© All Rights Reserved
Nous prenons très au sérieux les droits relatifs au contenu. Si vous pensez qu’il s’agit de votre contenu, signalez une atteinte au droit d’auteur ici.
Formats disponibles
Téléchargez aux formats PDF, TXT ou lisez en ligne sur Scribd
0% ont trouvé ce document utile (0 vote)
22 vues31 pages

Introduction au filtrage de Kalman

Le document présente le filtrage de Kalman, une méthode d'estimation de paramètres utilisée pour la fusion de données provenant de divers capteurs, essentielle dans les systèmes de navigation modernes. Il vise à illustrer le principe du filtrage de Kalman à travers des exemples simples, en établissant les équations nécessaires et en justifiant leur utilisation par la méthode des moindres carrés. Le texte aborde également des applications pratiques dans le domaine du positionnement dynamique, permettant au lecteur de développer ses propres applications.

Transféré par

Benjamin Jourdain
Copyright
© All Rights Reserved
Nous prenons très au sérieux les droits relatifs au contenu. Si vous pensez qu’il s’agit de votre contenu, signalez une atteinte au droit d’auteur ici.
Formats disponibles
Téléchargez aux formats PDF, TXT ou lisez en ligne sur Scribd

IG5

Sciences de la Mesure
 Positionnement dynamique 

Introduction au ltrage de Kalman

Jérôme Verdun

Rudolf Emil Kalman∗

* Photo réalisée par KuoShen Choong, chercheur au département d'ingénierie aéronautique de l'université de Sydney (Australie).
3

Avant-propos
Le ltrage de Kalman est une méthode d'estimation de paramètres dont l'un des principaux atouts
1

est de pouvoir prendre en compte, en temps réel, toute nouvelle observation. Cette méthode est lar-
gement répandue dans toutes les applications qui nécessitent la fusion de données en provenance de
divers capteurs - telle la robotique - et constitue l'outil de base de tous les systèmes de navigation
modernes qui équipent les satellites articiels, les missiles balistiques, les avions, les drones, les navires,
les submersibles habités et autonomes ainsi que les voitures  haut de gamme .
Le principe de base du ltrage de Kalman repose sur des idées simples, trop souvent occultées par la
lourdeur du formalisme mathématique. Ainsi, le présent document a-t-il pour objectif premier d'illus-
trer ce principe sur un exemple simple an que le lecteur en acquière une vision intuitive. Puis, les
équations du ltrage de Kalman sont établies et justiées à partir de la méthode des moindres car-
rés selon une approche développée par [Sage]. Parmi les nombreuses voies qui conduisent au ltre de
Kalman, cette approche progressive privilégie l'interprétation physique des équations qui conduit plus
aisément à la mise en ÷uvre du ltrage de Kalman dans les applications. La lecture assidue du présent
document devrait donc vous permettre de comprendre les utilisations du ltre de Kalman notamment
dans le domaine du positionnement dynamique, puis de développer assez rapidement vos propres ap-
plications.

Bonne lecture !

1. Rudolf Emil Kalman, ingénieur hongrois né en 1930 qui a travaillé aux États-Unis et en Suisse.
Table des matières

1 Navigation astronomique 7
1.1 Posons le problème . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.2 Ajoutons une mesure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3 De la dynamique en plus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.4 Cas limites . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2 Petit retour sur les moindres carrés 15


2.1 Moindres carrés sans pondération . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.2 Propriétés associées aux moindres carrés . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.2.1 Erreur d'estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.2.2 Résidu après estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.2.3 Espérance de l'erreur d'estimation . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2.4 Matrice de covariance de l'erreur d'estimation . . . . . . . . . . . . . . . . . . . 17
2.2.5 Matrice de covariance des résidus après estimation . . . . . . . . . . . . . . . . . 17
2.3 Moindres carrés avec pondération . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3 Mise à jour des paramètres 19


3.1 Cas statique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.1.1 Formulation générale . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.1.2 Cas limites . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.2 Cas dynamique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.2.1 Modèle d'évolution stochastique des paramètres . . . . . . . . . . . . . . . . . . 23
3.2.2 Solution à temps continu et discrétisation . . . . . . . . . . . . . . . . . . . . . 25
3.2.3 Meilleure estimation du vecteur des paramètres . . . . . . . . . . . . . . . . . . 27
3.2.4 Formulation du ltre de Kalman discret . . . . . . . . . . . . . . . . . . . . . . 28
3.2.5 Propriétés fondamentales du ltre de Kalman . . . . . . . . . . . . . . . . . . . 28
3.3 Récapitulation des équations du ltre de Kalman . . . . . . . . . . . . . . . . . . . . . 29
3.3.1 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.3.2 Remarques sur cette formulation . . . . . . . . . . . . . . . . . . . . . . . . . . 30

Bibliographie 31

5
Chapitre 1
Navigation astronomique
Cet exemple est librement inspiré de l'ouvrage de [Maybeck].

1.1 Posons le problème


An de comprendre le principe de fonctionnement du ltre de Kalman, imaginez que vous soyez un
navigateur perdu en pleine mer de nuit sans autre moyen de positionnement que la visée des étoiles
lointaines (gure 1.1). Pour simplier le problème, nous supposerons que votre navire se déplace dans
une dimension ce qui ne compromet en rien la généralité du problème. À l'instant t1 , vous estimez
que votre position est z1 . Pour rendre compte des erreurs aléatoires qui aectent votre mesure, vous
estimez également l'écart-type de votre mesure, soit σz1 . Soit x(t1 ) la position eective de votre navire
à l'instant t1 . La valeur réelle de x(t1 ) n'est pas connue, tout au plus les mesures peuvent-elles vous
en donner une estimation. Les valeurs de la position x(t1 ) sont donc distribuées de façon aléatoire,
autrement dit x(t1 ) est une variable aléatoire. Dans l'hypothèse où la distribution des erreurs aléatoires
suit une loi normale, la courbe représentant la densité de probabilité conditionnelle px(t1 )|z1 (x) de x(t1 )
sachant qu'une mesure z1 de cette dernière est disponible en fonction des valeurs possibles x de la
position est représentée sur la gure 1.2.

Figure 1.1  Perdu en pleine mer de nuit, comment retrouver votre route ?

7
8 CHAPITRE 1. NAVIGATION ASTRONOMIQUE

Cette densité présente un maximum en x = z1 ; la probabilité que la valeur de x(t1 ) se situe dans
l'intervalle [z1 − σz1 , z1 + σz1 ] (intervalle à 1σ ) est égale à 68, 3%. Vous pouvez donc admettre que la
meilleure estimation x̂(t1 ) de x(t1 ) compte tenu des informations disponibles à l'instant t1 (une seule
mesure) est donnée par
x̂(t1 ) = z1 , (1.1)
avec un écart-type associé σx (t1 ) donné par

σx (t1 ) = σz1 . (1.2)

px(t1)|z1 (x)
0.45

0.40

0.35

0.30

0.25

0.20

0.15

0.10

0.05 σz1
0.00
−2 0 2 4 z1 6 8 10 12
x

Figure 1.2  Représentation de la densité de probabilité conditionnelle px(t1 )|z1 (x) en fonction des
valeurs possibles x de la position. L'aire hachurée correspond à la probabilité que la variable aléatoire
x(t1 ) se situe dans l'intervalle [z1 − σz1 , z1 + σz1 ], soit 68, 3% pour une distribution normale.

1.2 Ajoutons une mesure


Un marin plus expérimenté que vous réalise ses propres visées à un instant t2 ≃ t1 de sorte que l'on
peut admettre que la position du navire n'a quasiment pas changé. Son estimation de la position est z2
avec un écart-type σz2 < σz1 . La densité de probabilité conditionnelle px(t2 )|z2 (x) est représentée sur la
gure 1.3. Comparée à la gure 1.2, cette courbe apparaît plus étroite autour du maximum x = z2 ; les
valeurs possibles de x(t2 ) sont donc moins dispersées autour de z2 que ne le sont celles de x(t1 ) autour
de z1 . L'estimation x̂(t2 ) = z2 est donc plus précise que x̂(t1 ) = z1 .
1.2. AJOUTONS UNE MESURE 9

px(t2)|z2 (x)
0.45

0.40

0.35

0.30
σz2
0.25

0.20

0.15

0.10

0.05

0.00
−2 0 2 4 6 z2 8 10 12 x
Figure 1.3  Représentation de la densité de probabilité conditionnelle px(t2 )|z2 (x) en fonction des
valeurs possibles x de la position. La mesure z2 étant plus précise que z1 , la courbe se resserre autour
du maximum x = z2 (trait continu) par rapport à la courbe précédente (trait pointillé).

La question cruciale à ce niveau est de savoir comment combiner les deux mesures z1 et z2 disponibles à
l'instant t2 pour déterminer une meilleure estimation de x(t2 ) ≃ x(t1 ) que celle obtenue avec seulement
l'une des deux mesures. D'après les résultats théoriques que nous établirons plus loin (cf chapitre 3),
l'estimation optimale de x(t2 ) compte tenu des mesures z1 et z2 disponibles correspond à la simple
moyenne pondérée µ suivante :
σz22 σz12
x̂(t2 ) = µ = z1 + z2 , (1.3)
σz12 + σz22 σz12 + σz22

avec un écart-type σ qui vérie :


1 1 1
2
= + . (1.4)
σ σz12 σz22
La densité de probabilité conditionnelle px(t2 )|z1 ,z2 (x) est alors représentée sur la gure 1.4. Elle cor-
respond toujours à une distribution normale avec pour paramètres la moyenne pondérée µ et l'écart-
type σ . De par sa dénition, l'écart-type σ est toujours plus petit que σz1 et σz2 (σ ≤ min(σz1 , σz2 )).
L'estimation de x(t2 ) qui tient compte des deux mesures est donc toujours plus précise que celles basées
sur seulement l'une des deux, et ce même si l'une des mesures est de très mauvaise qualité (σz1 ≫ σz2 ).

Si les deux mesures sont de même qualité alors σz1 = σz2 et les deux poids des mesures z1 et z2 sont
égaux à 12 dans l'équation 1.3 :
σz22 σz12 1
= = .
σz12 + σz22 σz12 + σz22 2
10 CHAPITRE 1. NAVIGATION ASTRONOMIQUE

px(t2)|z1,z2 (x)
0.45

0.40

0.35

0.30
σ
0.25

0.20

0.15

0.10

0.05

0.00
−2 0 2 4 6 µ 8 10 12 x
Figure 1.4  Représentation de la densité de probabilité conditionnelle px(t2 )|z1 ,z2 (x) en fonction des
valeurs possibles x de la position. Cette courbe est plus étroite autour du maximum x = µ que ne le
sont celles (trait pointillé) représentant la même distribution avec seulement une mesure disponible
(z1 ou z2 ).

L'estimation x̂(t1 ) correspond alors à la moyenne arithmétique des deux mesures :


z1 + z2
x̂(t1 ) = ,
2
σ σ
avec l'écart-type σx (t1 ) = √z1 = √z2 (cf équation 1.4).
2 2
En remarquant que le poids de la mesure z2 dans l'équation 1.3 peut également s'écrire :
σz12 σz2 σz22
= 1 ,
σz12 + σz22 σz22 σz12 + σz22
il est clair que
σz12 σz22
> ,
σz12 + σz22 σz12 + σz22
σz2
dès que 1
σz2
> 1 ⇒ σz1 > σz2 .
2
C'est donc la mesure la plus précise z2 (resp. la moins précise z1 ) qui intervient dans l'estimation
x̂(t1 ) avec le poids le plus fort (resp. le plus faible). L'estimation x̂(t2 ) = µ possède donc la propriété
attendue d'une combinaison pertinente : en tirant partie de la redondance, elle présente une précision
plus grande que n'importe quelle estimation n'impliquant que l'une des mesures de la combinaison

À partir de l'expression de x̂(t2 ), il est possible de donner une autre interprétation de la combinaison
1.3. Il vient très simplement :
σz12
x̂(t2 ) = z1 + (z2 − z1 ),
σz12 + σz22
1.3. DE LA DYNAMIQUE EN PLUS 11

soit, en tenant compte l'équation 1.1,


x̂(t2 ) = x̂(t1 ) + K(t2 )(z2 − x̂(t1 )), (1.5)
où l'on a posé
σz12
K(t2 ) = . (1.6)
σz12 + σz22
Ainsi, la meilleure estimation x̂(t2 ) qui prend en compte les deux mesures z1 et z2 est égale la somme
de la meilleure estimation x̂(t1 ) ne prenant en compte que la mesure z1 et d'un terme correctif pro-
portionnel à la diérence entre la nouvelle mesure z2 et l'estimation x̂(t1 ).
Le facteur de proportionnalité K(t2 ), appelé également gain de Kalman, dépend des incertitudes res-
pectives sur x̂(t1 ) et z2 . En eet, d'après les relations 1.2 et 1.6, il vient :
σx2 (t1 )
K(t2 ) = . (1.7)
σx2 (t1 ) + σz22

De même, l'évolution de la variance sur l'estimation peut être obtenue en exprimant la variance σx2 (t2 )
en fonction de σx2 (t1 ) à partir des relations 1.4 et 1.6 ; il vient successivement :
1 σz21 + σz22 1
2
= 2 2
=
σ σz1 σz2 K(t2 )σz2
 
1 1 1
= − ,
K(t2 ) σ 2 σz21

d'où, il résulte :  
1 1 1
−1 2
= ,
K(t2 ) σ K(t2 )σz21
soit
σ 2 = (1 − K(t2 ))σz21 .
Sachant que σx2 (t2 ) = σ 2 et σx2 (t1 ) = σz21 , il vient nalement :
σx2 (t2 ) = σx2 (t1 ) − K(t2 )σx2 (t1 ). (1.8)
Puisque 1 − K(t2 ) ≤ 1, on en déduit que σx2 (t2 ) ≤ σx2 (t1 ) ⇒ σx (t2 ) ≤ σx (t1 ), et ce d'autant plus que le
gain K(t2 ) est proche de 1. Il est intéressant de noter à ce stade que toute l'information que nous avons
sur notre position à l'instant t2 est donnée par les valeurs de l'estimation x̂(t2 ) et la variance σx2 (t2 ) ;
elles permettent en eet de xer complètement la densité de probabilité conditionnelle px(t2 )|z1 ,z2 (x)
représentée sur la gure 1.4.

1.3 De la dynamique en plus


Supposons à présent que vous naviguiez continûment sans chercher à modier votre vitesse. Le bon
modèle pour traduire l'évolution de votre position x(t) à l'instant t > t2 est alors donné par :
dx
= u + w, (1.9)
dt
où u représente la vitesse nominale du navire et w un bruit qui rend compte de l'incertitude sur la
connaissance de la vitesse réelle du navire étant donné les perturbations dues à la houle. Nous suppo-
serons que ce bruit de modèle a une distribution normale de moyenne nulle et de variance σw2 .
12 CHAPITRE 1. NAVIGATION ASTRONOMIQUE

La gure 1.5 montre l'évolution de la densité de probabilité conditionnelle px(t)|z1 ,z2 (x) d'après le
modèle d'évolution temporelle 1.9, depuis l'instant t2 de la dernière mesure jusqu'à l'instant t− 3 qui
précède l'acquisition d'une troisième mesure à l'instant t3 . L'estimation de la position est conforme
à un déplacement du navire à la vitesse constante u ; la relation entre les estimations x̂(t−
3 ) et x̂(t2 )
s'écrit donc :
x̂(t−
3 ) = x̂(t2 ) + u(t3 − t2 ). (1.10)
Dans le même temps, l'incertitude sur cette estimation croît de par la présence du bruit sur le modèle
et puisqu'aucune nouvelle information autre que ce modèle d'évolution ne vient la contraindre. La
3 ) et σx (t2 ) s'écrit alors :
relation qui lie les variances σx2 (t− 2

σx2 (t− 2 2 2
3 ) = σx (t2 ) + σw (t3 − t2 ) . (1.11)

px(t)|z1,z2 (x)
0.40

0.35

0.30

0.25

0.20
t2
0.15
t > t2
σ(t2) t−
3 > t
0.10

0.05 σ(t)
σ(t−
3)
0.00
−10 0 10 20 30 40 50 60 70 80 x

x̂(t2) x̂(t)
x̂(t−
3)

Figure 1.5  Représentation de l'évolution de la densité de probabilité conditionnelle px(t)|z1 ,z2 (x) en
fonction des valeurs possibles x de la position entre les instants t2 et t−
3 . La courbe s'élargit progressi-
vement autour du maximum x = x̂(t) ce qui indique que l'incertitude sur l'estimation de la précision
s'accroît de par la présence du bruit de modèle w sur la vitesse u.
1.4. CAS LIMITES 13

Ainsi, avant la mesure prévue à l'instant t3 , l'estimation de la position est-elle x̂(t3 ) avec l'écart-type
σx (t3 ).
Soit z3 la nouvelle mesure disponible à l'instant t3 . En appliquant rigoureusement les relations 1.5, 1.8
et 1.7, l'estimation de la position à l'instant t3 peut s'écrire :

x̂(t3 ) = x̂(t− −
3 ) + K(t3 )(z3 − x̂(t3 )), (1.12)

avec pour variance :


σx2 (t3 ) = σx2 (t− 2 −
3 ) − K(t3 )σx (t3 ), (1.13)
où le gain de Kalman K(t3 ) est donné par :
σx2 (t−
3)
K(t3 ) = − . (1.14)
σx (t3 ) + σz23
2

L'équation 1.12 rend compte de la modication à apporter à l'estimation x̂(t− 3 ) donnée par le modèle
d'évolution lors de l'apparition d'une nouvelle mesure ; elle permet donc d'apporter une correction
3 ) à l'estimation prédite σx (t3 ). L'eet de l'apport d'une nouvelle mesure est donc une
K(t3 )σx2 (t− 2 −

diminution de la variance proportionnelle à σx2 (t− 3 ), égale à K(t3 )σx (t3 ). Le rôle fondamental du gain de
2 −

Kalman apparaît donc clairement puisque ce dernier intervient à la fois sur l'amplitude de la correction
et celle de la diminution de la variance (et donc l'augmentation de la précision) de l'estimation.

1.4 Cas limites


L'examen de cas limites à la lumière des relations 1.12, 1.13 et 1.14 est très instructive sur le compor-
tement du ltre de Kalman.

Supposons tout d'abord que la mesure z3 soit de mauvaise qualité par rapport à l'estimation dis-
ponible x̂(t−
3 ), ce qui se traduit par σz3 ≫ σx (t3 ). La relation 1.14 indique alors que K(t3 ) ≪ 1. Dans

le cas limite où σz23 → ∞, K(t3 ) → 0 et, d'après 1.12, x̂(t3 ) = x̂(t−


3 ) avec σx (t3 ) = σx (t3 ).

La nouvelle mesure n'a alors pas été prise en compte dans l'estimation x̂(t3 ) qui reste égale à celle
prédite x̂(t−
3 ). Le modèle d'évolution est ici prépondérant pour la détermination de l'estimation x̂(t3 ).

Supposons à présent que ce soit le modèle d'évolution  vitesse constante  qui soit à son tour peu
able, ce qui se traduit par σw2 ≫ σx (t2 ). Dans ce cas, σx2 (t−
3 ) ≫ 1 (relation 1.10) et K(t3 ) ≫ 1. À la
limite, il vient
σw2 → ∞ ⇒ σx2 (t−
3 ) → ∞ ⇒ K(t3 ) → 1.

Il en résulte, d'après 1.12 :


x̂(t3 ) = x̂(t− −
3 ) + (z3 − x̂(t3 )) = z3 .

Le mesure z3 est donc ici prépondérante puisque c'est elle qui xe la valeur de l'estimation x̂(t3 ) avec
une incertitude nulle σx (t3 ) = 0 (équation 1.13). Le modèle d'évolution est donc ici complètement
occulté.

Enn, si une estimation de la position devient quasiment certaine juste avant la mesure z3 , autre-
ment dit σx (t−
3 ) = 0, alors le gain de Kalman est nul K(t3 ) = 0 et plus aucune mesure supplémentaire
ne pourra changer la valeur de l'estimation x̂(t3 ) = x̂(t−
3 ).

Bien que nous n'ayons encore rien démontré, l'examen de ce petit exemple permet d'illustrer les aspects
principaux du ltrage de Kalman : une méthode d'estimation qui permet d'estimer un paramètre qui
14 CHAPITRE 1. NAVIGATION ASTRONOMIQUE

varie dans le temps en combinant de façon optimale des mesures disponibles et un modèle d'évolution
temporelle tous deux aectés de leurs propres incertitudes. Les équations 1.10 à 1.14 respectent la
structure des phases dites de prédiction et de correction de l'algorithme du ltre de Kalman qui sont
présentes dans toutes les formulations de ce ltre.

Fin du chapitre 1
Chapitre 2
Petit retour sur les moindres carrés
Les développements mathématiques présentés dans la suite sont tirés de [Rogers].

An de comprendre la notion de combinaison optimale évoquée dans l'exemple du chapitre 1, il nous
faut revenir sur la méthode d'estimation par moindres carrés. Dans ce qui suit, nous supposerons que
les paramètres estimés par a méthode ne varient pas dans le temps.

2.1 Moindres carrés sans pondération


Soit à estimer les n composantes inconnues (x1 , x2 , . . . , xn ) d'un vecteur x ∈ Rn à partir des p mesures
(z1 , z2 , . . . , zp ) qui forment le vecteur z ∈ Rp , sachant que le vecteur des paramètres x est lié au vecteur
des observations par une relation de la forme :

z = Hx + v, (2.1)
où H est une matrice de dimension p × n et v ∈ Rp le vecteur des résidus dont les composantes sont
les diérences entre les observations z et leur modèle linéaire Hx.
Nous supposerons, en outre, que les mesures sont redondantes, c'est-à-dire p ≫ n et que la matrice H
est de rang égal à n.
La meilleure estimation x̂ du vecteur des paramètres x au sens des moindres carrés est celle qui, si elle
existe, minimise la fonction coût JLS dénie par :
JLS (x) = (z − Hx)T (z − Hx) = ∥z − Hx∥22 , (2.2)
où (z − Hx)T désigne la transposée de la matrice z − Hx et ∥ . . . ∥2 la norme 2 (norme euclidienne)
sur Rp .

Le taux de variation de cette fonction coût entre les points x et x + ∆x s'exprime par :
JLS (x + ∆x) − JLS (x) = (z − H(x + ∆x))T (z − H(x + ∆x)) − (z − Hx)T (z − Hx)
= (z − Hx)T − (H∆x)T ((z − Hx) − (H∆x)) − (z − Hx)T (z − Hx)


= −(z − Hx)T H∆x − (H∆x)T (z − Hx) + (H∆x)T H∆x


T
= − (z − Hx)T H∆x + ((z − Hx)T H∆x + ∆xT H T H∆x.
T
En remarquant que (z − Hx)T H∆x est un scalaire, il vient (z − Hx)T H∆x = (z − Hx)T H∆x .
La fonction coût peut donc se développer pour tout x ∈ Rn et tout ∆x ∈ Rn par :
JLS (x + ∆x) = JLS (x) − 2(z − Hx)T H∆x + ∆xT H T H∆x. (2.3)

15
16 CHAPITRE 2. PETIT RETOUR SUR LES MOINDRES CARRÉS

Les deuxième et troisième terme de cette relation correspondent donc respectivement aux diérentielles
première et seconde de la fonction coût. Les extremums de la fonction coût satisfont donc l'équation
matricielle (z − Hx)T H = 0TRn où 0Rn désigne le vecteur nul de Rn . En prenant la transposée membre
à membre de cette dernière égalité, on obtient le système normal d'équations :
H T Hx = H T z. (2.4)
Si le rang de H est égal à n, alors la matrice carrée de dimension n, H T H , est inversible. Le système
normal d'équations admet alors une solution unique x̂, extremum de la fonction coût, qui s'exprime
par :
x̂ = (H T H)−1 H T z. (2.5)
De plus, la matrice H H est également symétrique, dénie et positive. Ainsi, pour tout ∆x ∈ Rn ,
T

∆xT H T H∆x ≥ 0 et ∆xT H T H∆x = 0 uniquement lorsque ∆x = 0Rn . On en déduit que x̂ correspond
au minimum de la fonction coût et constitue l'estimation du vecteur des paramètres x par moindres
carrés.
Exercice 1 (Théorème fondamental) : Montrer que si H est une matrice de dimension p × n de rang
n alors H T H est une matrice inversible de dimension n, symétrique, dénie et positive.

2.2 Propriétés associées aux moindres carrés


2.2.1 Erreur d'estimation
L'erreur d'estimation ε est dénie par ε = x − x̂.
Il vient alors successivement :
ε = x − x̂
= x − (H T H)−1 H T z
= x − (H T H)−1 H T (Hx + v)
= x − (H T H)−1 H T H x − (H T H)−1 H T v
| {z }
In

où In désigne la matrice identité de dimension n.


Il vient nalement :
ε = −(H T H)−1 H T v (2.6)
.

2.2.2 Résidu après estimation


Le vecteur des résidus après estimation v̂ est dénie par v̂ = z − H x̂.
Il vient successivement :
v̂ = Hx + v − H x̂
= H(x − x̂) + v
= Hε + v
= −H(H T H)−1 H T v + v,
soit nalement,
v̂ = (Ip − H(H T H)−1 H T )v (2.7)
où Ip désigne la matrice identité de dimension p.
2.3. MOINDRES CARRÉS AVEC PONDÉRATION 17

2.2.3 Espérance de l'erreur d'estimation


Le vecteur des résidus est supposé être d'espérance nulle, soit E(v) = 0.
La matrice de covariance Σv du vecteur des résidus v est alors donnée par Σv = E(vv T ). La diagonale
de cette matrice est constituée par les variances sur les observations σzi , i = 1, 2, . . . , p xées a priori
compte tenu de ce que nous connaissons de la qualité des observations. Elle peut être complétée en ses
termes non diagonaux cov(zi , zj ) pour (i, j) ∈ {1, 2, . . . , p}2 , i ̸= j , lorsque ces dernières sont connues.
En utilisant la linéarité de l'espérance, il vient successivement, d'après 2.6 :
E(x − x̂) = −E (H T H)−1 H T v
 

= −(H T H)−1 H T E(v)


= 0Rn .
L'erreur d'estimation est d'espérance nulle. L'estimateur des moindres carrés est dit sans biais.

2.2.4 Matrice de covariance de l'erreur d'estimation


La matrice de covariance Σε de l'erreur d'estimation ε est dénie par :
Σε = E (x − x̂)(x − x̂)T .
 

La diagonale de cette matrice contient toutes les variances des paramètres estimés σx̂2i , i = 1, 2, . . . , n.
Les termes non diagonaux correspondent aux covariances entre les paramètres estimés cov(x̂i , x̂j ) pour
(i, j) ∈ {1, 2, . . . , n}2 , i ̸= j .
Il vient successivement d'après 2.6 :
Σε = E ((H T H)−1 H T v)(v T H(H T H)−1 )
 

= (H T H)−1 H T E(vv T )H(H T H)−1 ,


soit nalement,
Σε = (H T H)−1 H T Σv H(H T H)−1 . (2.8)

2.2.5 Matrice de covariance des résidus après estimation


La matrice de covariance des résidus après estimation est dénie par :
Σv̂ = E (z − H x̂)(z − H x̂)T ,
 

puisque, d'après 2.7, E(v̂) = (Ip − H(H T H)−1 H T )E(v) = 0Rp .

Il vient alors :
Σv̂ = (Ip − H(H T H)−1 H T )E(vv T )(Ip − H(H T H)−1 H T )T ,
d'où nalement,
Σv̂ = (Ip − H(H T H)−1 H T )Σv (Ip − H(H T H)−1 H T )T . (2.9)

2.3 Moindres carrés avec pondération


La recherche de l'estimation des moindres carrés peut être modiée pour tenir compte de la qualité
des mesures grâce à une matrice de poids, carrée de dimension p, symétrique, inversible, utilisée pour
exprimer la nouvelle fonction coût dénie par :
JW LS (x) = (z − Hx)T W (z − Hx). (2.10)
18 CHAPITRE 2. PETIT RETOUR SUR LES MOINDRES CARRÉS

Exercice 2 (Solution des moindres carrés pondérés) : Démontrer que l'estimation des moindres carrés
pondérés est donnée par :
x̂ = (H T W H)−1 H T W z, (2.11)
et que la matrice de covariance de l'erreur d'estimation Σε s'exprime par :
Σε = (H T W H)−1 H T W Σv W H(H T W H)−1 .

Dans le cas très fréquent en pratique où la matrice de poids est prise égale à l'inverse de la matrice de
covariance des résidus W = Σ−1v , il vient :

Σε = (H T W H)−1 H T W
| W
−1 T −1 T −1 T T −1
{z } W H(H W H) = |(H W H){z H W H}(H W H) ,
Ip In

d'où il résulte :
Σε = (H T W H)−1 . (2.12)
Exercice 3 (Formule de propagation de la variance) : Soit X un vecteur aléatoire et ΣX sa matrice
de covariance. On suppose que le vecteur Y est lié au vecteur X par la relation matricielle Y = AX
où A est une matrice complètement déterministe.
h i
Exprimer la matrice de covariance ΣY de Y à partir de la relation ΣY = E (Y − E(Y )) (Y − E(Y ))T
en fonction de ΣX . La relation obtenue est appelée formule de propagation de la variance.
Exercice 4 : Démontrer les relations 2.8 et 2.9 en utilisant la formule de propagation de la variance.

Fin du chapitre 2
Chapitre 3
Mise à jour des paramètres
Nous allons à présent étudier le problème de la mise à jour de l'estimation de paramètres par moindres
carrés lorsque de nouvelles observations sont disponibles. Nous envisagerons successivement deux cas,
le premier dans lequel les paramètres seront supposés indépendants du temps (cas statique), le second
dans lequel les paramètres seront susceptibles de varier dans le temps en suivant une loi d'évolution
connue (cas dynamique). Pour alléger les notations, nous noterons P la matrice de covariance de
l'erreur d'estimation Σε et R la matrice de covariance des résidus a priori Σv .

3.1 Cas statique


3.1.1 Formulation générale
Soit x̂m l'estimation des moindres carrés obtenue à partir des m observations qui entre dans la relation
z = Hx + v dénie au chapitre 2. Cette dernière peut s'écrire sous la forme matricielle suivante :
      
z1 h11 h12 · · · h1n x1 v1
 z2   h21 h22 · · · h2n   x2   v2 
 ..  =  .. .. .. ..   ..  +  ..  .
      
 .   . . . .  .   . 
zm hm1 hm2 · · · hmn xn vm

Supposons à présent qu'une nouvelle (m + 1)ième mesure zm+1 soit disponible de sorte qu'il faille
modier la relation z = Hx + v en
     
z H v
= x+ ,
zm+1 hT vm+1

où l'on a posé hT = [hm+1,1 hm+1,2 · · · hm+1,n ].

On se propose de déterminer la nouvelle estimation par moindres carrés x̂m+1 sous la forme :

x̂m+1 = x̂m + ∆x, (3.1)

où ∆x correspond à une correction de l'estimation obtenue avec seulement m mesures.

La fonction coût associée à ce problème peut s'écrire sous la forme :


    T  −1      
z H R 0m,1 z H
JRW LS (x) = − x − x , (3.2)
zm+1 hT 0Tm,1 r−1 zm+1 hT

19
20 CHAPITRE 3. MISE À JOUR DES PARAMÈTRES

où R désigne la matrice de covariance sur les résidus des m observations zi , i = 1, 2, . . . , m, r la va-


riance sur la mesure m + 1, zm+1 et 0m,1 = [0, 0, . . . , 0]T .
| {z }
m fois  
R 0m,1
La forme de la nouvelle matrice de covariance des résidus d'observation R = indique
0Tm,1 r
que la nouvelle mesure zm+1 est supposée complètement décorrélée des m autres.

En notant temporairement    
z H
Z= et H = ,
zm+1 hT
et en utilisant la relation 2.11, il vient :

x̂m+1 = (HT R−1 H)−1 HT R−1 Z.

Sous forme matricielle, cette relation s'écrit :


T  !−1  T 
R−1 0m,1 R−1 0m,1
   
H H H z
x̂m+1 = .
hT 0Tm,1 r−1 hT hT 0Tm,1 r−1 zm+1

En développant cette dernière expression, il vient successivement :


−1
R−1 H R−1 z
   
 T
  T 
x̂m+1 = H h H h
r−1 hT r−1 zm+1
−1
H T R−1 H + hr−1 hT H T R−1 z + hr−1 zm+1 ,

=

soit nalement :
x̂m+1 = Pm+1 H T R−1 z + hr−1 zm+1 , (3.3)


où Pm+1 correspond à la matrice de covariance de l'erreur d'estimation associée à x̂m+1 .

Cette dernière s'exprime par :


−1 −1
Pm+1 = HT R−1 H = H T R−1 H + hr−1 hT = (Pm−1 + hr−1 hT )−1 ,

où la quantité H T R−1 H correspond à l'inverse de la matrice de covariance de l'erreur d'estimation


associée à x̂m , c'est-à-dire Pm−1 .

Il vient donc nalement :


−1
Pm+1 = Pm−1 + hr−1 hT . (3.4)
Avant de poursuivre la formulation de l'estimation x̂m+1 , il convient de démontrer un lemme qui
permettra de transformer la relation 3.4.
Exercice 5 (Lemme) : Soient A, B et D, trois matrices inversibles telles que :

A−1 = B −1 + C T D−1 C

où C est une matrice telle que D + CBC T soit inversible.


Alors
A = B − BC T (D + CBC T )−1 CB.
3.1. CAS STATIQUE 21

En vériant les hypothèses du lemme pour A = Pm+1 , B = Pm , C = hT et D = r (relation 3.4), on


constate tout d'abord que :
D + CBC T = r + hT Pm h.
Cette dernière quantité est un scalaire non nul comme somme d'un réel positif non nul r (variance
de zp+1 ) et d'une forme quadratique positive hT Pm h ; l'inverse de r + hT Pm h existe et la matrice
D + CBC T est inversible.

Le lemme permet donc d'établir l'identité suivante :


Pm+1 = Pm − Pm h(r + hT Pm h)−1 hT Pm . (3.5)
En utilisant l'équation précédente, il est possible d'exprimer diéremment la relation 3.3 donnant x̂m+1
par
x̂m+1 = Pm+1 H T R−1 z + hr−1 zm+1


Pm − Pm h(r + hT Pm h)−1 hT Pm H T R−1 z + hr−1 zm+1


 
=
= Pm H T R−1 z +Pm hr−1 zm+1 − Pm h(hT Pm h + r)−1 hT Pm H T R−1 z
| {z }
x̂m

−Pm h(hT Pm h + r)−1 hT Pm r−1 zm+1 .

En factorisant Pm dans tous les termes en facteur de hr−1 zm+1 , il vient :


x̂m+1 = x̂m + Pm In − h(hT Pm h + r)−1 hT Pm hr−1 zm+1
 

−Pm h(hT Pm h + r)−1 hT x̂m .

Le deuxième terme de cette dernière relation peut s'écrire en distribuant h à gauche :


Pm In − h(hT Pm h + r)−1 hT Pm hr−1 zm+1 = Pm h − h(hT Pm h + r)−1 hT Pm h r−1 zm+1 ,
   

où l'expression entre crochets correspond au vecteur h de dimension n multiplié par le scalaire


hT Pm h r
1 − (hT Pm h + r)−1 hT Pm h = 1 − T
= T .
h Pm h + r h Pm h + r
Il vient alors
x̂m+1 = x̂m + Pm h(hT Pm h + r)−1 zm+1 − Pm h(hT Pm h + r)−1 hT x̂m .

Les équations de mise à jour de l'estimation x̂m+1 peuvent nalement s'écrire :


x̂m+1 = x̂m + km (zm+1 − hT x̂m ), (3.6)
avec km = Pm h(hT Pm h + r)−1 ,
et Pm+1 = Pm − km hT Pm = (In − km hT )Pm ,

où In désigne la matrice identité de dimension n.

Ce sont là les équations fondamentales qui permettent de calculer, à partir d'une estimation établie
pour un ensemble de mesures, la nouvelle estimation lorsqu'une nouvelle mesure est disponible.
Remarque 1 : Le vecteur ∆x dénie par la relation 3.1 est nalement donné par :

∆x = km (zm+1 − hT x̂m ) avec km = Pm h(hT Pm h + r)−1 .


22 CHAPITRE 3. MISE À JOUR DES PARAMÈTRES

Remarque 2 : Le vecteur des paramètres estimés x̂ est appelé vecteur d'état. Les quantités x̂m et
x̂m+1 correspondent donc à deux états qui dièrent par leur nombre de mesures (une de plus pour
x̂m+1 ).

Remarque 3 : La quantité km correspond au . Ce dernier dépend de l'incertitude


gain de Kalman

sur la nouvelle mesure par l'intermédiaire de la variance r présente dans le terme (hT Pm h + r)−1 .
Le terme restant hT Pm h correspond à la matrice de covariance du vecteur hT x̂m selon la formule de
propagation de la variance. Ce vecteur fournit une estimation de la nouvelle observation zm+1 à partir
des paramètres x̂m estimés avec les m mesures zi , i = 1, 2, . . . , m ; la matrice hT Pm h correspond donc
à l'incertitude sur cette observation.
Remarque 4 : La matrice de covariance de l'erreur d'estimation Pm voit parfois sa dénomination
raccourcie en covariance de l'erreur . Il ne faut néanmoins moins jamais oublier qu'il s'agit bien d'une
matrice.
Exercice 6 : À l'aide des relations 3.6, retrouver les résultats de l'exemple de la navigation astrono-
mique traitée au chapitre 1.
Exercice 7 : Vous disposez de m mesures zi , i = 1, 2, . . . , m d'une grandeur x réputée constante.
Quelle est la meilleure estimation de x au sens des moindres carrés en supposant que toutes les me-
sures sont décorrélés, de variance respectives σi2 , i = 1, 2, . . . , m ? Comment est modiée cette estima-
tion lorsque l'on dispose d'une mesure zm+1 supplémentaire, décorrélée des autres et de variances σm+1 2
?

On précisera dans chacun des cas la matrice de covariance de l'erreur d'estimation.


Exercice 8 : Montrer que le gain de Kalman km peut également s'exprimer par km = Pm+1 hr−1 .

La démarche présentée ici peut aisément se généraliser au cas où tout un vecteur de nouvelles mesures
est disponible.
Exercice 9 : On suppose à présent qu'un nouvel ensemble de mesures est disponible représenté par
le vecteur zs de dimension s. La fonction coût de ce problème d'estimation peut alors s'écrire sous la
forme matricielle suivante :
    T  −1      
z H R 0m,s z H
JRW LS (x) = − x − x ,
zs Hs 0Tm,s Rs−1 zs Hs

où Hs est une matrice de dimension s × n et Rs la matrice de covariance des résidus a priori sur les s
nouvelles mesures.

On appelle respectivement x̂ et x̂s les estimations obtenues avant et après l'ajout des nouvelles mesures
et P et Ps les matrices de covariance de l'erreur associées à ces estimations. En utillisant la même
démarche que celle employée dans le paragraphe 3.1.1, démontrer les relations suivantes :
x̂s = x̂ + ks (zs − Hs x̂), (3.7)
avec ks = P HsT (Hs P HsT + Rs )−1 ,
et Ps = (In − ks Hs )P.

3.1.2 Cas limites


Envisageons les deux cas limites correspondant à une nouvelle mesure zm+1 de très bonne qualité
r ≪ hT Pm h et de très mauvaise qualité r ≫ hT Pm h.
3.2. CAS DYNAMIQUE 23

Dans le premier cas :


km → Pm h(hT Pm h)−1 ⇒ hT km → hT Pm h(hT Pm h)−1 = 1.

À partir des relations 3.6, en multipliant par hT à gauche l'équation donnant x̂m+1 , il vient :
hT k →1
hT x̂m+1 = hT x̂m + hT km (zm+1 − hT x̂m ) = hT x̂m + hT km zm+1 − hT km hT x̂m −→
m
zm+1 ,

puis en multipliant la matrice Pm+1 par hT à gauche et par h à droite, il vient :


hT k →1
hT Pm+1 h = (hT − hT km hT )Pm h −→
m
0.

Ce résultat indique que la mesure zm+1 est obtenue sans erreur par hT x̂m+1 avec une variance hT Pm+1 h
nulle. L'estimation x̂m+1 est en quelque sorte contrainte de restituer sans erreur la mesure réputée de
très grande qualité.

Dans le second cas :


km → Pm hr−1 .
Il vient donc :
x̂m+1 → x̂m + Pm hr−1 (zm+1 − hT x̂m ) = (In − r−1 Pm hhT )x̂m + r−1 Pm hZm+1 ≈ In x̂m = x̂m ,

et
Pm+1 → (In − r−1 Pm hhT )Pm ≈ In Pm = Pm .
Ce résultat indique qu'une nouvelle mesure zm+1 de mauvaise qualité n'est que peu prise en compte
dans l'estimation x̂m+1 puisque à la limite, cette dernière égale x̂m , estimation obtenue sans la mesure
zm+1 , avec la même incertitude (Pm+1 = Pm ).

Ainsi, les deux cas limites de prépondérance de la mesure et du modèle d'évolution illustrés dans
le chapitre 1, trouvent-ils là leur justication. Le réglage d'un ltre de Kalman consiste à mettre en
balance les contributions de la mesure et du modèle d'évolution à un instant donné en jouant sur
les niveaux de bruit accordés respectivement aux mesures (bruit de mesure ) et au modèle d'évolution
(bruit de modèle ).
Exercice 10 : Retrouver les résultats de l'étude des cas limites de l'exemple du navire (cf chapitre 1)
à partir des résultats généraux précédemment établis.

3.2 Cas dynamique


3.2.1 Modèle d'évolution stochastique des paramètres
Nous supposons à présent que les paramètres à estimer varient dans le temps. Le vecteur des paramètres
x ∈ Rn devient alors la fonction vectorielle suivante :
x : R → Rn
t 7→ x(t) = [x1 (t), x2 (t), . . . , xn (t)]T .

Un modèle d'évolution temporelle de x(t) est alors constitué par un système d'équations diérentielles
impliquant les paramètres à estimer.
24 CHAPITRE 3. MISE À JOUR DES PARAMÈTRES

Dans l'exemple du chapitre 1, la trajectoire du navire est supposée parcourue  à vitesse constante 
ce qui a été traduit par le modèle d'évolution :
dx
= u,
dt
où u représente la vitesse du navire.

Il s'agit là d'une équation diérentielle linéaire du premier ordre, qui est ici complètement déter-
ministe ; en eet, la connaissance d'une condition initiale x(t0 ) = x0 où t0 est l'instant où la position

du navire est x0 permet de déterminer, de façon exacte, l'évolution de la position du navire à un instant
ultérieur ; en eet, la postion à l'instant t est donnée par :

x(t) = x0 + u(t − t0 ).

L'incertitude sur la pertinence de ce modèle a été prise en compte en introduisant le bruit w sous la
forme :
dx
= u + w,
dt
ce qui rend alors le modèle d'évolution non plus déterministe, mais stochastique.

Dans la suite, nous supposerons que le vecteur des paramètres x(t) à l'instant t satisfait le système
diérentiel du premier ordre suivant :
dx
= −Ax(t) + Bu(t) + Cw(t), (3.8)
dt
où A est une matrice carrée de dimension n, B , C deux matrices de dimensions respectives n × q et
n × l, u(t), un fonction vectorielle déterministe de R → Rq appelée commande, et w(t) un vecteur de
bruit blanc à l composantes.
Remarque 5 : Le choix d'un système diérentiel linéaire du premier ordre n'est pas réducteur
puisque :
 tout système non linéaire peut se ramener à un système linéaire au voisinage d'une solution
approchée ;
 toute équation diérentielle linéaire d'ordre n peut se ramener à la résolution d'un système
diérentiel de n équations à n inconnues.
Exercice 11 : Vérier la première assertion de la remarque précédente en donnant le système dié-
rentiel correspondant à l'équation diérentielle suivante :
d2⃗r
= −k⃗r,
dt2
où k est une constante positive et ⃗r le rayon vecteur de coordonnées cartésiennes ⃗r = [x(t), y(t)]T .
dx dy
On utilisera le vecteur d'état X = [x(t), ẋ(t), y(t), ẏ(t)]T où ẋ(t) = et ẏ(t) = .
dt dt
Exercice 12 : Reprendre l'exercice précédent avec l'équation diérentielle suivante :

d2⃗r
= −k∥⃗r∥2⃗r.
dt2
3.2. CAS DYNAMIQUE 25

Remarque 6 : Le bruit blanc w(t) est un processus aléatoire qui satisfait les deux propriétés suivantes :

∀t ∈ R, E[w(t)] = 0 ;
∀(t, t ) ∈ R , E[w(t)wT (t′ )] = Q(t)δ(t − t′ ),
′ 2

où Q(t) est une matrice de covariance de dimension l et δ la distribution de Dirac.

Ce bruit est souvent appelé bruit de modèle .

3.2.2 Solution à temps continu et discrétisation


La solution du système 3.8 peut s'exprimer à partir de la connaissance de la condition initiale x(t0 ) =
x0 ∈ Rn sous la forme suivante :
Z t Z t
x(t) = e −(t−t0 )A
x(t0 ) + e −(t−τ )A
Bu(τ )dτ + e−(t−τ )A Cw(τ )dτ. (3.9)
t0 t0

Remarque 7 : Dans l'expression e−(t−t0 )A x(t0 ), le terme e−(t−t0 )A représente l'exponentielle de matrice
dénie par :
1 1
e−(t−t0 )A = In + (t − t0 )A + ((t − t0 )A)2 + . . . + ((t − t0 )A)i + . . . , i ∈ N.
2! i!
Il s'agit donc d'une matrice carrée de même dimension que A. Il conviendra de ne pas inverser e−(t−t0 )A
et x(t0 ) dans le produit e−(t−t0 )A x(t0 ).
Exercice 13 : Dans le cas où les matrices A, B et C qui interviennent dans la relation 3.8 dépendent
explicitement du temps, montrer que la solution du système diérentiel peut s'écrire :
Rt Z t Rt
Z t Rt

x(t) = e t0 A(τ )dτ
x(t0 ) + e− τ A(τ ′ )dτ ′
B(τ )u(τ )dτ + e− τ A(τ ′ )dτ ′
C(τ )w(τ )dτ. (3.10)
t0 t0

An de préciser la meilleure estimation de x(t) selon ce modèle d'évolution temporelle, il convient de
poser l'hypothèse suivante :

(3.11)
 
∀τ ∈ R, E (x(t0 ) − E[x(t0 )])wT (τ ) = 0n,l ,

où 0n,l est la matrice nulle de dimension n × l.

De plus, puisque E[x(t0 )] n'est pas aléatoire,

E (x(t0 ) − E[x(t0 )])wT (τ ) = E x(t0 )wT (τ ) − E(x(t0 )) E(wT (τ )),


   
| {z }
0Tl
R

donc :
∀τ ∈ R, E x(t0 )wT (τ ) = 0n,l .


Cette hypothèse traduit la complète décorrélation entre le bruit w(t) avec la condition initiale x(t0 ) et
donc avec tout les vecteurs de la forme e−(t−t0 )A x(t0 ) puisque pour tout τ ∈ R :

E e−(t−t0 )A x(t0 )wT (τ ) = e−(t−t0 )A E x(t0 )wT (τ ) = 0n,l .


   
| {z }
0n,l
26 CHAPITRE 3. MISE À JOUR DES PARAMÈTRES

De plus, le bruit w(t) est naturellement décorrélé avec le terme de commande de l'équation 3.9 puisque
pour tout t′ ∈ R :
Z t   Z t
−(t−τ )A ′
E e Bu(τ )dτ T
w (t ) = e−(t−τ )A Bu(τ ) E(wT (t′ )) dτ = 0n,1 .
t0 t0 | {z }
0T
l,1

Enn, l'étude de la corrélation avec le terme de bruit de l'équation 3.9 conduit, pour tout t′ ∈ R, à :
Z t   Z t
−(t−τ )A ′
E e Cw(τ )dτ T
w (t ) = e−(t−τ )A CE(w(τ )wT (t′ ))dτ
t0 t
Z 0t
= e−(t−τ )A CQ(t′ )δ(t′ − τ )dτ.
t0

An de garantir également la décorrélation du bruit de modèle w(t) et le bruit de mesure représenté
par le vecteur v(t) dans l'équation d'observation z(t) = Hx(t) + v(t) où z(t) est le vecteur des m
observations disponibles à l'instant t et H une matrice de dimension m × n, la condition suivante est
également imposée :
∀(t, t′ ) ∈ R2 , E w(t)v T (t′ ) = 0. (3.12)
 

À ce stade, on peut utiliser l'équation d'évolution 3.9 pour déterminer le vecteur des paramètres à
l'instant tk+1 , x(tk+1 ), en fonction de fonction de celui obtenu à l'instant tk , x(tk ) où tk et tk+1 , k ∈ N,
sont deux instants consécutifs de la suite croissante d'instants (tn )n∈N . Nous verrons que cette suite
correspond aux instants de disponibilité des mesures. Il vient alors :
Z tk+1 Z tk+1
−(tk+1 −tk )A −(tk+1 −τ )A
x(tk+1 ) = e x(tk ) + e Bu(τ )dτ + e−(tk+1 −τ )A Cw(τ )dτ.
tk tk

Posons xk+1 = x(tk+1 ) et xk = x(tk ) vecteurs d'état respectivement aux instants tk et tk+1 . On dénit
également la matrice de tansition Φk:k+1 entre les états k et k + 1 par la relation suivante :
Φk:k+1 = e−(tk+1 −tk )A . (3.13)
Enn, étant donné la proximité des instants tk et tk+1 en pratique, il est légitime de réaliser les
approximations ci-après :
Z tk+1
e−(tk+1 −τ )A Bu(τ )dτ ≈ (tk+1 − tk )Bu(tk+1 ) = Bk+1 uk+1 ;
tk

Z tk+1
e−(tk+1 −τ )A Cw(τ )dτ ≈ (tk+1 − tk )Cw(tk+1 ) = Ck+1 wk+1 ,
tk

où l'on a posé :
Bk+1 = (tk+1 − tk )B et Ck+1 = (tk+1 − tk )C, (3.14)
puis,
uk+1 = u(tk+1 ) et wk+1 = w(tk+1 ). (3.15)
En tenant compte des équations 3.13, 3.14 et 3.15, l'équation d'évolution peut s'écrire sous la forme :
xk+1 = Φk:k+1 xk + Bk+1 uk+1 + Ck+1 wk+1 . (3.16)
Cette dernière équation correspond à une forme discrétisée de l'équation d'évolution 3.9. Cette mo-
délisation discrète est absolument nécessaire pour tenir compte du fait que les mesures ne sont pas
continûment disponibles dans le temps, mais seulement à des instants discrets.
3.2. CAS DYNAMIQUE 27

3.2.3 Meilleure estimation du vecteur des paramètres


La meilleure estimation x̂k+1 du vecteur des paramètres xk+1 étant donné le modèle d'évolution 3.16
s'obtient en calculant l'espérance du vecteur xk+1 , soit x̂k+1 = E(xk+1 ).

En appliquant l'opérateur linéaire espérance uniquement sur les termes aléatoires, il vient :
x̂k+1 = E (Φk:k+1 xk + Bk+1 uk+1 + Ck+1 wk+1 ) = Φk:k+1 E(xk ) + Bk+1 uk+1 + Ck+1 E(wk+1 ),
et puisque E(xk ) = x̂k et E(wk+1 ) = 0, il vient nalement :
x̂k+1 = Φk:k+1 x̂k + Bk+1 uk+1 . (3.17)
Déterminons à présent la matrice de covariance de l'erreur d'estimation Pk+1 dénie par :
Pk+1 = E (xk+1 − x̂k+1 )(xk+1 − x̂k+1 )T .
 

Avec les équations 3.16 et 3.17, l'erreur d'estimation dans l'état k + 1, εk+1 est donnée par :
εk+1 = xk+1 − x̂k+1 = Φk:k+1 xk + Bk+1 uk+1 + Ck+1 wk+1 − (Φk:k+1 x̂k + Bk+1 uk+1 )
= Φk:k+1 (xk − x̂k ) + Ck+1 wk+1 = Φk:k+1 εk + Ck+1 wk+1 ,
où l'on a posé εk = xk − x̂k .

La détermination de la matrice Pk+1 peut donc se poursuivre selon les calculs suivants :
Pk+1 = E(εk+1 εTk+1 ) = E (Φk:k+1 εk + Ck+1 wk+1 )(Φk:k+1 εk + Ck+1 wk+1 )T
 

= Φk:k+1 E(εk εTk )ΦTk:k+1 + Φk:k+1 E(εk wk+1


T T
)Ck+1 +
T T T T
Ck+1 E(wk+1 εk )Bk+1 + Ck+1 E(wk+1 wk+1 )C .
On remarque tout d'abord que E(εk εTk ) correspond à la matrice de covariance de l'erreur de l'état k ,
soit Pk et que E(wk+1 wk+1
T
) n'est rien d'autre que la matrice de covariance du bruit blanc w(t), soit
Qk+1 .
T
D'autre part, explicitons le terme E(εk wk+1
T
) = E wk+1 εTk :
T T T T
E(εk wk+1 ) = E(xk wk+1 ) − x̂k E(wk+1 ) = E(xk wk+1 ).
| {z }
0Tl
R

Puisque le bruit w est décorrélé de x(t) à tout instant, il vient : E(xk wk+1
T
) = 0n,l où 0n,l est la matrice
nulle de dimension n × l ; on en déduit nalement que :
T
) = 0n,l et E wk+1 εTk = 0Tn,l .

E(εk wk+1
Finalement, l'équation de propagation de la variance entre les états k et k + 1 est donnée par :
Pk+1 = Φk:k+1 Pk ΦTk:k+1 + Ck+1 Qk+1 Ck+1
T
. (3.18)
Les équations 3.17 et 3.18 constituent les équations de propagation de l'estimation des paramètres au
cours du temps, sachant que ces derniers obéissent à un modèle stochastique d'évolution temporelle.
Exercice 14 : En utilisant les équations 3.17 et 3.18, démontrer les relations ci-après utilisées dans
l'exemple du navire au chapitre 1 :
x̂(t−
3 ) = x̂(t2 ) + u(t3 − t2 ),

avec
σx2 (t− 2 2 2
3 ) = σx (t2 ) + σw (t3 − t2 ) .
28 CHAPITRE 3. MISE À JOUR DES PARAMÈTRES

3.2.4 Formulation du ltre de Kalman discret


Nous sommes à présent en mesure de donner les équations complètes du ltre de Kalman discret.
Soit x̂k le vecteur des paramètres estimés à l'instant tk (vecteur d'état dans l'état k ) et Pk la matrice
de covariance de l'erreur d'estimation associée à x̂k . L'évolution du vecteur des paramètres entre les
instants tk et t− k+1 est supposé suivre le modèle d'évolution temporelle 3.16, de sorte que si xk+1 et Pk+1
− −

désigne respectivement le vecteur des paramètres et la matrice de covariance de l'erreur à l'instant


k+1 , il vient :
t−

x̂−
k+1 = Φk:k+1 x̂k + Bk+1 uk+1 ,

Pk+1 = Φk:k+1 Pk ΦTk:k+1 + Ck+1 Qk+1 Ck+1
T
.

À l'instant tk+1 , les m mesures disponibles sont représentées par le vecteur


zk+1 = [zk+1,1 , zk+1,2 , . . . , zk+1,m ]T ,

de matrice de covariance Rk+1 ; ce dernier est supposé lié au vecteur des paramètres xk+1 par :
zk+1 = Hk+1 xk+1 + vk ,

où Hk+1 est une matrice de dimension m × n et vk le vecteur des m résidus a priori à l'instant tk+1 .

Le vecteur zk+1 vient s'ajouter au k vecteurs de mesures déjà acquis z1 , z2 , . . . , zk où zi ∈ Rm , i =


1, 2, . . . , k , jusqu'à l'instant t−
k . D'après les relations 3.7, la meilleure estimation x̂k+1 des paramètres
étant donné le nouveau vecteur de mesures à l'instant tk+1 est donnée par :
x̂k+1 = x̂− −

k+1 + Kk+1 zk+1 − Hk+1 x̂k+1 ,

où l'on a posé le gain de Kalman Kk+1 égal à :


− −
Kk+1 = Pk+1 T
Hk+1 (Hk+1 Pk+1 T
Hk+1 + Rk+1 )−1 .

L'équation de propagation de la matrice de covariance entre les instants t−


k+1 et tk s'écrit alors :


Pk+1 = (In − Kk+1 Hk+1 )Pk+1 .

3.2.5 Propriétés fondamentales du ltre de Kalman


Le ltre de Kalman possède des propriétés fondamentales qui portent sur l'erreur d'estimation
[Simon]. Soient xk le vecteur des paramètres recherchés et x̂k son estimation réalisée avec le ltre de
Kalman. L'erreur d'estimation ek s'exprime par :
ek = xk − x̂k . (3.19)
Puisqu'il existe une incertitude sur le modèle d'évolution, le vecteur xk est aléatoire. De même, le
vecteur estimé x̂k est aléatoire puisque dépendant des mesures aectées par une erreur aléatoire. Le
vecteur d'erreur ek est donc aléatoire en tant que diérence de deux vecteurs aléatoire.

Supposons que l'on cherche, à chaque instant k , l'estimateur x̂k de sorte à minimiser l'espérance
E eTk Sk ek où Sk est une matrice symétrique, dénie, positive. Remarquons que si Sk est une matrice
diagonale, l'expression de E eTk Sk ek s'écrit :


E eTk Sk ek = Sk (1)E e2k (1) + Sk (2)E e2k (2) + . . . + Sk (n)E e2k (n) ,
   
3.3. RÉCAPITULATION DES ÉQUATIONS DU FILTRE DE KALMAN 29

où Sk (1), Sk (2), . . . , Sk (n) sont les éléments de la diagonale de Sk et ek (1), ek (2), . . . , ek (n) sont les
composantes de ek .

Les propositions ci-après constituent les réponses au problème de minimisation ci-dessus.


1. Si les bruits de modèle {wk } et de mesure {vk } sont blancs gaussiens, à moyenne nulle, et non
corrélés, alors le ltre de Kalman est la solution au problème posé.
2. Si les bruits de modèle {wk } et de mesure {vk } sont blancs, à moyenne nulle, et non corrélés, alors
le ltre de Kalman est le meilleur ltre linéaire, solution au problème posé. Il peut cependant
exister des ltres non linéaires qui donnent une meilleure solution.
3. Si les bruits de modèle {wk } et de mesure {vk } sont corrélés ou colorés, le ltre de Kalman doit
être modié pour répondre au problème posé.
4. Si les équations d'évolution et/ou de mesures sont non linéaires, il faut utiliser les formulations
non linéaires du ltre de Kalman.
Dans l'équation de mise à jour x̂k+1 = x̂− −
, la quantité −
 
k+1 + K k+1 zk+1 − H k+1 x̂ k+1 zk+1 − Hk+1 x̂ k+1
est appelée innovation. C'est eectivement la part des mesures qui contient de nouvelles informations
sur les paramètres estimés. L'innovation est un processus aléatoire  qui correspond à un bruit blanc
centré dont la matrice de covariance est Hk+1 Pk+1 Hk + Rk+1 . En fait, le ltre de Kalman agit
− T

comme un algorithme qui  blanchit  les mesures et donc en extrait le maximum d'informations.
L'analyse statistique des innovations permet de vérier la conformité de la loi de probabilité empirique
de l'innovation avec un bruit blanc centré avec la bonne matrice de variance-covariance. Cette analyse
permet en particulier de vérier le bon fonctionnement d'un ltre de Kalman.

3.3 Récapitulation des équations du ltre de Kalman


Les diérentes étapes du ltrage de Kalman sont ici retranscrites sous la forme d'un algorithme en vue
de sa programmation.

3.3.1 Formulation
Hypothèses fondamentales
Les bruits de modèle wk et de mesure vk sont des bruits blancs décorrélés entre eux :

∀(k, k ′ ) ∈ N2 , E(wk ) = 0Rl et E(wk ) = 0Rm ,


E(wk wkT′ ) = 0l,l pour k ̸= k ′ et E(wk wkT ) = Qk ,
E(vk vkT′ ) = 0m,m pour k ̸= k ′ et E(vk vkT ) = Rk ,
E(wk vkT′ ) = 0l,l .

Si x̂0 = E(x(t0 )) désigne la condition initiale alors ∀k ∈ N, E(wk x̂T0 ) = 0l,n .

Initialisation
Le premier vecteur d'état est xé par la condition initiale : x̂0 = E(x(t0 )) ; la matrice de covariance de
l'erreur d'estimation est P0 .

Prédiction

x̂−
k+1 = Φk:k+1 x̂k + Bk+1 uk+1 ,

Pk+1 = Φk:k+1 Pk ΦTk:k+1 + Ck+1 Qk+1 Ck+1
T
.
30 CHAPITRE 3. MISE À JOUR DES PARAMÈTRES

Correction

x̂k+1 = x̂− −

k+1 + Kk+1 zk+1 − Hk+1 x̂k+1 ,
− −
Kk+1 = Pk+1 T
Hk+1 (Hk+1 Pk+1 T
Hk+1 + Rk+1 )−1 ,

Pk+1 = (In − Kk+1 Hk+1 )Pk+1 .

3.3.2 Remarques sur cette formulation


L'expression de la matrice de variance-covariance donnée par :

Pk+1 = (In − Kk+1 Hk+1 )Pk+1

bien que très simple à programmer, ne garantit pas que la matrice obtenue soit symétrique ou dénie,
positive. Aussi Peter Joseph a proposé dans les années 60 une formulation dite stabilisée qui s'exprime
par [Simon] :

Pk+1 = (In − Kk+1 Hk+1 )Pk+1 (In − Kk+1 Hk+1 )T + Kk+1 Rk+1 Kk+1
T
. (3.20)

La formulation donnée par l'équation (3.20) est celle à privilégier dans toute programmation de l'al-
gorithme du ltre de Kalman.

La formulation du ltre de Kalman indique que Pk+1 −


, Kk+1 et Pk+1 ne dépendent pas des mesures
représentées par le vecteur zk+1 . La calcul du gain Kk+1 peut donc être calculé avant l'arrivée des
mesures ce qui constitue un gain de temps essentiel notamment pour les applications temps réel. En
outre, la connaissance de la matrice de variance-covariance Pk+1 renseigne sur les performances du
ltre en termes de barres d'erreur sur les paramètres estimés par le ltre de Kalman. La possibilité de
calculer cette matrice sans avoir réalisé eectivement les mesures autorise la réalisation de simulation.

Fin du chapitre 3
Bibliographie
[Sage] ,
Sage, A. P. et Melsa, J. L. Estimation Theory with Applications to Communications and
Control, Éd. MacGraw-Hill, New York, 1971.

[Maybeck] Maybeck Peter S., Stochastic Models, Estimation, and Control, Mathematics in
Science and Engineering, Volume 141-1, Éd. Richard Bellman, Université de Californie Sud,
1979 ; republié par Navtech Book and Software Store, 1994.
Premier chapitre disponible à l'adresse : [Link]
[Rogers] Rogers Robert M., Applied Mathematics in Integrated Navigation Systems, Éd. AIAA
Education Series, Joseph A. Schetz, Reston, Virginie, 2003.

[Simon] ,
Simon Dan Optimal State Estimation , Éd. John Wiley & Sons, 2006.

**********

31

Vous aimerez peut-être aussi