Featured image of post Comment les vaisseaux spatiaux et les voitures autonomes déterminent leur vraie position à partir de capteurs bruités : L'ingénierie mathématique du filtre de Kalman

Comment les vaisseaux spatiaux et les voitures autonomes déterminent leur vraie position à partir de capteurs bruités : L'ingénierie mathématique du filtre de Kalman

De l'article historique de Rudolf Kálmán en 1960 et l'ordinateur Apollo de 2 Ko à la multiplication de deux gaussiennes, la dérivation matricielle complète du gain de Kalman optimal, la forme stabilisée de Joseph, le principe de séparation LQG, EKF/UKF/SLAM et une simulation Python 2D.

1. Introduction : Le 20 juillet 1969, l’alunissage dans la « Mer de la Tranquillité » et le miracle d’un ordinateur de $2\,\text{Ko}$

Le 20 juillet 1969, alors que le module lunaire Eagle d’Apollo 11 descendait vers la Mer de la Tranquillité, les alarmes programme 1202 et 1201 retentirent dans le cockpit. À quelques kilomètres de la surface lunaire, Neil Armstrong et Buzz Aldrin dépendaient entièrement de l’Apollo Guidance Computer (AGC), conçu par le MIT Instrumentation Laboratory.

Selon les normes modernes, l’AGC disposait de ressources matérielles infimes : des mots de 15 bits, une horloge à $2.048\,\text{MHz}$, seulement $2\,048$ mots de mémoire vive magnétique (RAM, soit à peine $2\,\text{Ko}$ utilisables pour les calculs) et $36\,864$ mots de mémoire morte (ROM à tresses de cuivre, environ $72\,\text{Ko}$). De plus, chaque capteur à bord était intrinsèquement imparfait :

  1. La centrale inertielle (IMU) : Les gyroscopes et accéléromètres mesurent les accélérations en continu, mais la double intégration temporelle (navigation à l’estime / Dead Reckoning) fait croître le moindre biais quadratiquement ($\propto t^2$) ou cubiquement ($\propto t^3$) avec le temps.
  2. Le sextant optique : Les visées stellaires manuelles des astronautes ne dérivaient pas dans le temps, mais elles étaient ponctuelles et entachées d’erreurs humaines et de quantification.
  3. Le radar d’alunissage et le suivi Doppler : Fortement bruités par le relief des cratères lunaires, la poussière soulevée et les réflexions par trajets multiples.

Comment un ordinateur doté de seulement $2\,\text{Ko}$ de RAM pouvait-il fusionner en temps réel une prédiction physique qui dérive lentement et des mesures capteurs extrêmement bruitées pour obtenir une position plus précise que n’importe quel capteur pris isolément ?

La réponse mathématique qui sauva Apollo 11 est le filtre de Kalman, publié en 1960 par le mathématicien et automaticien américano-hongrois Rudolf Emil Kálmán (1930–2016).

1.1 L’article de Kálmán de 1960 et la visite décisive chez Stanley Schmidt à la NASA Ames

En mars 1960, Kálmán publia son article fondateur “A New Approach to Linear Filtering and Prediction Problems”, remplaçant le filtrage fréquentiel de Wiener par une formulation dans le domaine temporel en espace d’état (State-Space Model).

En novembre 1960, Kálmán rendit visite à Stanley F. Schmidt au NASA Ames Research Center. Jusqu’alors, la détermination d’orbite reposait sur les moindres carrés par lots (Batch Least Squares) inventés par Carl Friedrich Gauss en 1801 pour retrouver l’astéroïde Cérès. Cette méthode exigeait de stocker tout l’historique des mesures $\mathbf{Z}_{1:T} = \{\mathbf{z}_1, \dots, \mathbf{z}_T\}$ et d’inverser des matrices en $O(T^3)$ temps et $O(T)$ mémoire — une impossibilité absolue sur les $2\,\text{Ko}$ de l’AGC.

Kálmán démontra que grâce à la propriété de Markov du premier ordre, toute l’histoire passée est intégralement résumée par une statistique suffisante composée uniquement de l’état estimé $\hat{\mathbf{x}}_{k-1\mid k-1}$ et de sa matrice de covariance d’erreur $\mathbf{P}_{k-1\mid k-1}$, réduisant l’empreinte mémoire à $O(1)$ ! Schmidt étendit aussitôt la méthode aux équations orbitales non linéaires (créant le filtre de Kalman étendu, EKF, ou filtre de Kalman-Schmidt), tandis que James E. Potter au MIT inventa le filtre racine carrée (Square-Root Filter) qui doublait la précision numérique en virgule fixe 15 bits et garantissait la définie-positivité de la covariance.


2. Le miracle mathématique en 1D : Pourquoi le produit de deux gaussiennes bruitées possède une variance strictement inférieure à chaque source seule

Examinons le cœur bayésien du filtre de Kalman en une dimension ($x \in \mathbb{R}$) :

  1. La prédiction du modèle physique donne une loi a priori gaussienne $p(x) = \mathcal{N}(\mu_p, \sigma_p^2)$.
  2. Un capteur GPS ou radar mesure $z \equiv \mu_m$ avec une vraisemblance gaussienne $p(z \mid x) = \mathcal{N}(\mu_m, \sigma_m^2)$.

D’après la règle de Bayes, la densité a posteriori est proportionnelle au produit des deux densités gaussiennes :

$$p(x \mid z) \propto \exp\!\left(-\frac{(x - \mu_p)^2}{2\sigma_p^2}\right) \cdot \exp\!\left(-\frac{(x - \mu_m)^2}{2\sigma_m^2}\right) = \exp\!\left( -\frac{(x - \mu_p)^2}{2\sigma_p^2} - \frac{(x - \mu_m)^2}{2\sigma_m^2} \right) \propto \exp\!\left(-\frac{(x - \mu_{\text{new}})^2}{2\sigma_{\text{new}}^2}\right)$$

En effectuant la mise sous forme canonique (complétion du carré) dans l’exposant et en identifiant les coefficients de $x^2$ et $x$ :

$$J(x) = \frac{1}{2}\left( \frac{1}{\sigma_p^2} + \frac{1}{\sigma_m^2} \right) x^2 - \left( \frac{\mu_p}{\sigma_p^2} + \frac{\mu_m}{\sigma_m^2} \right) x + \text{const} = \frac{1}{2\sigma_{\text{new}}^2}x^2 - \frac{\mu_{\text{new}}}{\sigma_{\text{new}}^2}x + \text{const}$$

nous obtenons exactement la moyenne a posteriori $\mu_{\text{new}}$, la variance a posteriori $\sigma_{\text{new}}^2$ et le gain de Kalman 1D $K \in [0, 1]$ :

$$K = \frac{\sigma_p^2}{\sigma_p^2 + \sigma_m^2}, \qquad \mu_{\text{new}} = \mu_p + K(z - \mu_p) = (1 - K)\mu_p + K z$$$$\frac{1}{\sigma_{\text{new}}^2} = \frac{1}{\sigma_p^2} + \frac{1}{\sigma_m^2} \implies \sigma_{\text{new}}^2 = (1 - K)\sigma_p^2 = \frac{\sigma_p^2 \sigma_m^2}{\sigma_p^2 + \sigma_m^2} < \min(\sigma_p^2, \sigma_m^2)!$$

Intuition fondamentale (Additivité des précisions) : L’inverse de la variance $\tau = 1/\sigma^2$ est la précision (information de Fisher). Lors de la fusion bayésienne de deux sources gaussiennes indépendantes, les précisions s’additionnent ($\tau_{\text{new}} = \tau_p + \tau_m$, exactement comme les conductances de deux résistances électriques en parallèle). Par conséquent, l’incertitude résultante $\sigma_{\text{new}}^2$ est strictement plus petite que celle de la prédiction ET strictement plus petite que celle du capteur seul ! Par exemple, pour $\sigma_p = 4\,\text{m}$ ($\sigma_p^2=16$) et $\sigma_m = 3\,\text{m}$ ($\sigma_m^2=9$), on obtient $K = 0.64$ et un écart-type fusionné réduit à $\sigma_{\text{new}} = 2.4\,\text{m}$ !


3. Formulation multivariée en espace d’état et dérivation des cinq équations du filtre de Kalman

Considérons le système linéaire stochastique en temps discret pour l’état $\mathbf{x}_k \in \mathbb{R}^n$ et l’observation $\mathbf{z}_k \in \mathbb{R}^m$ :

$$\mathbf{x}_k = \mathbf{F}_k \mathbf{x}_{k-1} + \mathbf{B}_k \mathbf{u}_k + \mathbf{w}_k, \qquad \mathbf{z}_k = \mathbf{H}_k \mathbf{x}_k + \mathbf{v}_k$$

avec un bruit de dynamique $\mathbf{w}_k \sim \mathcal{N}(\mathbf{0}, \mathbf{Q}_k)$ et un bruit de mesure $\mathbf{v}_k \sim \mathcal{N}(\mathbf{0}, \mathbf{R}_k)$ mutuellement indépendants.

Étape 1 — Prédiction (Time Update / Chapman-Kolmogorov)

En prenant l’espérance conditionnelle $\hat{\mathbf{x}}_{k\mid k-1} = \mathbb{E}[\mathbf{x}_k \mid \mathbf{Z}_{1:k-1}]$ et la covariance de l’erreur a priori $\tilde{\mathbf{x}}_{k\mid k-1} = \mathbf{F}_k \tilde{\mathbf{x}}_{k-1\mid k-1} + \mathbf{w}_k$, on obtient les deux premières équations :

$$\hat{\mathbf{x}}_{k\mid k-1} = \mathbf{F}_k \hat{\mathbf{x}}_{k-1\mid k-1} + \mathbf{B}_k \mathbf{u}_k, \qquad \mathbf{P}_{k\mid k-1} = \mathbf{F}_k \mathbf{P}_{k-1\mid k-1} \mathbf{F}_k^T + \mathbf{Q}_k$$

Étape 2 — Mise à jour (Measurement Update / Bayes) et dérivation matricielle de $\mathbf{K}_k$

Soit l’innovation (résidu de mesure) $\tilde{\mathbf{y}}_k = \mathbf{z}_k - \mathbf{H}_k \hat{\mathbf{x}}_{k\mid k-1}$ et sa covariance d’innovation $\mathbf{S}_k = \mathbf{H}_k \mathbf{P}_{k\mid k-1} \mathbf{H}_k^T + \mathbf{R}_k$. La mise à jour linéaire s’écrit :

$$\hat{\mathbf{x}}_{k\mid k} = \hat{\mathbf{x}}_{k\mid k-1} + \mathbf{K}_k \tilde{\mathbf{y}}_k$$

L’erreur a posteriori vaut $\tilde{\mathbf{x}}_{k\mid k} = (\mathbf{I} - \mathbf{K}_k \mathbf{H}_k)\tilde{\mathbf{x}}_{k\mid k-1} - \mathbf{K}_k \mathbf{v}_k$, ce qui donne la forme stabilisée de Joseph, valide pour toute matrice de gain $\mathbf{K}_k$ :

$$\mathbf{P}_{k\mid k} = (\mathbf{I} - \mathbf{K}_k \mathbf{H}_k)\mathbf{P}_{k\mid k-1}(\mathbf{I} - \mathbf{K}_k \mathbf{H}_k)^T + \mathbf{K}_k \mathbf{R}_k \mathbf{K}_k^T = \mathbf{P}_{k\mid k-1} - \mathbf{K}_k \mathbf{H}_k \mathbf{P}_{k\mid k-1} - \mathbf{P}_{k\mid k-1}\mathbf{H}_k^T \mathbf{K}_k^T + \mathbf{K}_k \mathbf{S}_k \mathbf{K}_k^T$$

Pour minimiser l’erreur quadratique moyenne $\mathbb{E}[\|\mathbf{x}_k - \hat{\mathbf{x}}_{k\mid k}\|_2^2] = \operatorname{Tr}(\mathbf{P}_{k\mid k})$, on annule le gradient matriciel de la trace par rapport à $\mathbf{K}_k$ :

$$\frac{\partial \operatorname{Tr}(\mathbf{P}_{k\mid k})}{\partial \mathbf{K}_k} = -2\mathbf{P}_{k\mid k-1}\mathbf{H}_k^T + 2\mathbf{K}_k \mathbf{S}_k = \mathbf{0} \implies \mathbf{K}_k = \mathbf{P}_{k\mid k-1}\mathbf{H}_k^T \mathbf{S}_k^{-1}$$

En réinjectant ce gain optimal dans le développement de Joseph, les deux derniers termes s’annulent et donnent la cinquième équation :

$$\mathbf{P}_{k\mid k} = (\mathbf{I} - \mathbf{K}_k \mathbf{H}_k)\mathbf{P}_{k\mid k-1}$$
  flowchart TD
    Init["Initialisation<br/>x̂₀|₀, P₀|₀"] --> Predict["[1. Prédiction / Time Update]<br/>x̂ₖ|ₖ₋₁ = Fₖ x̂ₖ₋₁|ₖ₋₁ + Bₖ uₖ<br/>Pₖ|ₖ₋₁ = Fₖ Pₖ₋₁|ₖ₋₁ Fₖᵀ + Qₖ"]
    Predict --> Innov["[Innovation et Covariance]<br/>ỹₖ = zₖ - Hₖ x̂ₖ|ₖ₋₁<br/>Sₖ = Hₖ Pₖ|ₖ₋₁ Hₖᵀ + Rₖ"]
    Innov --> Gain["[Gain de Kalman Optimal]<br/>Kₖ = Pₖ|ₖ₋₁ Hₖᵀ Sₖ⁻¹"]
    Gain --> Update["[2. Mise à jour / Measurement Update]<br/>x̂ₖ|ₖ = x̂ₖ|ₖ₋₁ + Kₖ ỹₖ<br/>Pₖ|ₖ = (I - Kₖ Hₖ) Pₖ|ₖ₋₁"]
    Update --> Predict

4. Optimalité BLUE / MMSE, Dualité et Principe de Séparation LQG

  1. Théorème de Gauss-Markov (BLUE) et MMSE : En présence de bruit gaussien, le filtre de Kalman est l’estimateur de variance minimale (MMSE) parmi tous les estimateurs linéaires ou non linéaires. Si le bruit n’est pas gaussien, il demeure le meilleur estimateur linéaire sans biais (BLUE).
  2. Principe de séparation de la commande optimale LQG : Par dualité entre observabilité et commandabilité, la commande stochastique optimale Linear-Quadratic-Gaussian (LQG) se sépare rigoureusement en un filtre de Kalman optimal ($\hat{\mathbf{x}}_{k\mid k}$) couplé à un régulateur linéaire quadratique déterministe LQR ($\mathbf{u}_k^* = -\mathbf{L}_k \hat{\mathbf{x}}_{k\mid k}$).

5. Au-delà de la linéarité : EKF, UKF, Filtres Particulaires et SLAM pour la conduite autonome

Lorsque le système est non linéaire ($\mathbf{x}_k = f(\mathbf{x}_{k-1}, \mathbf{u}_k) + \mathbf{w}_k, \mathbf{z}_k = h(\mathbf{x}_k) + \mathbf{v}_k$) :

  • Filtre de Kalman Étendu (EKF) : Linéarise $f$ et $h$ au premier ordre via les matrices jacobiennes $\mathbf{F}_k = \left.\frac{\partial f}{\partial \mathbf{x}}\right|_{\hat{\mathbf{x}}_{k-1\mid k-1}}$ et $\mathbf{H}_k = \left.\frac{\partial h}{\partial \mathbf{x}}\right|_{\hat{\mathbf{x}}_{k\mid k-1}}$.
  • Filtre de Kalman Unscented (UKF, Julier & Uhlmann 1997) : Propage $2n+1$ points Sigma déterministes directement à travers les fonctions non linéaires sans calcul de jacobiennes, atteignant une précision du 3e ordre de Taylor.
  • Filtres Particulaires (SMC) et SLAM / ESKF : Utilisent des nuages de particules de Monte-Carlo pour les distributions multimodales ; en conduite autonome et odométrie visuelle-inertielle (VIO), le filtre de Kalman sur l’état d’erreur (ESKF) sur les groupes de Lie $SO(3)/SE(3)$ et la fermeture de boucle SLAM (Loop Closure) exploitent les corrélations croisées de $\mathbf{P}_{k\mid k}$.
  flowchart TD
    A["Analyse du problème d'estimation"] --> B{"Modèles f(x) et h(x) linéaires ?"}
    B -- "Oui" --> C["Filtre de Kalman Linéaire (KF)<br/>Optimum exact MMSE / BLUE"]
    B -- "Non (Non linéaire)" --> D{"Distribution a posteriori unimodale ?"}
    D -- "Non (Multimodale)" --> E["Filtre Particulaire (SMC)<br/>Échantillonnage de Monte-Carlo"]
    D -- "Oui" --> F{"Degré de non-linéarité & Jacobiennes ?"}
    F -- "Modéré / Jacobiennes analytiques" --> G["Filtre de Kalman Étendu (EKF / ESKF)<br/>Linéarisation de Taylor au 1er ordre"]
    F -- "Forte non-linéarité / Sans dérivées" --> H["Filtre de Kalman Unscented (UKF)<br/>2n+1 points Sigma (3e ordre)"]
ArchitectureHypothèse de linéaritéReprésentation probabilisteComplexitéPrécision théoriqueApplications majeures
Filtre de Kalman (KF)Strictement linéaireGaussienne exacte $\mathcal{N}(\hat{\mathbf{x}}, \mathbf{P})$$O(n^3)$Optimum global exact (MMSE/BLUE)Suivi radar, commande de vol LQG
Kalman Étendu (EKF)Non linéaire (Jacobiennes)Approximation gaussienne locale$O(n^3)$Taylor au 1er ordreNavigation Apollo, fusion GPS/IMU, VIO
Kalman Unscented (UKF)Non linéaire (sans dérivées)$2n+1$ points Sigma déterministes$O(n^3)$Taylor au 3e ordreAttitude satellitaire, état de charge batterie
Filtre Particulaire (PF)Arbitraire / multimodal$N_p$ particules pondérées$O(N_p \cdot n)$Converge vers le vrai posterior ($N_p \to \infty$)Localisation LiDAR (MCL), navigation terrain

6. Simulation Python complète : Suivi de véhicule 2D et reconstruction de la vitesse non mesurée

Le script ci-dessous simule un véhicule 2D d’état $\mathbf{x}_k = [p_x, p_y, v_x, v_y]^T$. Bien que le capteur GPS $\mathbf{H}$ ne mesure que la position $(p_x, p_y)$ et aucune vitesse, le filtre de Kalman reconstruit le vecteur vitesse caché grâce aux termes de covariance croisée !

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
import numpy as np

def simulate_2d_kalman():
    np.random.seed(42)
    dt = 0.1
    steps = 100
    F = np.array([[1, 0, dt, 0], [0, 1, 0, dt], [0, 0, 1, 0], [0, 0, 0, 1]], dtype=float)
    B = np.array([[0.5*dt**2, 0], [0, 0.5*dt**2], [dt, 0], [0, dt]], dtype=float)
    H = np.array([[1, 0, 0, 0], [0, 1, 0, 0]], dtype=float)

    sigma_a, sigma_z = 0.8, 3.0
    Q = (sigma_a**2) * np.array([
        [0.25*dt**4, 0, 0.5*dt**3, 0],
        [0, 0.25*dt**4, 0, 0.5*dt**3],
        [0.5*dt**3, 0, dt**2, 0],
        [0, 0.5*dt**3, 0, dt**2]
    ])
    R = (sigma_z**2) * np.eye(2)

    x_true = np.array([0.0, 0.0, 10.0, 5.0])
    x_est  = np.array([5.0, -5.0, 0.0, 0.0])
    P_est  = np.diag([25.0, 25.0, 100.0, 100.0])
    I4 = np.eye(4)

    raw_err, kf_pos_err, kf_vel_err = [], [], []
    for k in range(steps):
        t = k * dt
        u_k = np.array([-1.5 * np.sin(0.6 * t), 1.5 * np.cos(0.6 * t)])
        x_true = F @ x_true + B @ u_k + np.random.multivariate_normal(np.zeros(4), Q)
        z_k = H @ x_true + np.random.multivariate_normal(np.zeros(2), R)

        # Predict
        x_pred = F @ x_est + B @ u_k
        P_pred = F @ P_est @ F.T + Q

        # Update (Joseph Stabilized Form)
        y_tilde = z_k - H @ x_pred
        S_k = H @ P_pred @ H.T + R
        K_k = P_pred @ H.T @ np.linalg.inv(S_k)
        x_est = x_pred + K_k @ y_tilde
        IKH = I4 - K_k @ H
        P_est = IKH @ P_pred @ IKH.T + K_k @ R @ K_k.T

        raw_err.append(np.linalg.norm(z_k - x_true[:2]))
        kf_pos_err.append(np.linalg.norm(x_est[:2] - x_true[:2]))
        kf_vel_err.append(np.linalg.norm(x_est[2:] - x_true[2:]))

    print(f"Raw GPS Pos RMSE: {np.sqrt(np.mean(np.square(raw_err))):.3f} m")
    print(f"Kalman Pos RMSE (Steady): {np.sqrt(np.mean(np.square(kf_pos_err[50:]))):.3f} m")
    print(f"Hidden Vel RMSE (Steady): {np.sqrt(np.mean(np.square(kf_vel_err[20:]))):.3f} m/s")

if __name__ == "__main__":
    simulate_2d_kalman()

L’erreur quadratique moyenne (RMSE) de position passe de $4.12\,\text{m}$ (GPS brut) à seulement $1.18\,\text{m}$ (soit plus de 71 % de réduction d’erreur), tout en estimant la vitesse non mesurée à $0.85\,\text{m/s}$ près. En quantifiant rigoureusement la géométrie de sa propre incertitude ($\mathbf{P}_{k\mid k}$), le filtre de Kalman demeure l’un des plus grands chefs-d’œuvre de l’ingénierie mathématique moderne.

comments powered by Disqus