Exécuter un filtre de Kalman en Python : bruit, données manquantes et valeurs aberrantes

Le réglage de Kalman exprime l'incertitude d'un modèle de mouvement et de ses observations. Cette expérience estime la position et la vitesse, puis compare les prédictions en cas de données manquantes, les réactions à une valeur aberrante et le rejet dû à une surconfiance. Toutes les entrées sont synthétiques ; les valeurs obtenues ne constituent pas des garanties de précision matérielle.

Modèle et unités

L'état comprend la position [m] et la vitesse [m/s]. Seule la position est observée, toutes les 0,1 s :

L'écart-type des observations est de 0,7 m, ce qui donne R = 0,49 m². La matrice de covariance du processus est q ggᵀ avec g = (dt²/2, dt), représentant une perturbation d'accélération indépendante à chaque étape. Ici, q = 0,4 a pour unités (m/s²)² ; il ne s'agit pas d'une densité spectrale de bruit blanc continu. Ne transposez pas cette valeur à une autre discrétisation sans en vérifier la signification. Consultez l'introduction au filtre de Kalman pour les calculs de prédiction et de mise à jour.

Pour la théorie, voir Simo Särkkä (2013), Filtrage et lissage bayésiens, section 4.3. La trajectoire, les conditions de porte et de comparaison utilisées ici sont issues de notre propre conception expérimentale.

Quatre configurations, entrée identique

Enregistrez le script d'expérience dans un répertoire de sortie vide et exécutez :

python3 engineering_labs.py

L'exécution a été vérifiée avec Python 3.12.3, NumPy 1.26.4 et Matplotlib 3.6.3. Le fichier génère également les autres expériences d'ingénierie liées. Il utilise une graine de 42 et 300 échantillons. La vitesse réelle passe de 1 à 1,4 m/s à 12 s. Trente observations, comprises entre 7,0 et 9,9 s, sont manquantes, et 12 m sont ajoutés à l'observation à 16 s. La position et la vitesse initiales sont respectivement de 0 m et 1 m/s ; la matrice de covariance initiale a pour éléments diagonaux 1 m² et 1 (m/s)². Chaque configuration reçoit exactement les mêmes observations.

Normalisation du résidu

Soit e le résidu d'observation avant correction et S sa variance prédite. Le seuil rejette une observation lorsque :

Ce seuil s'applique à ce résidu normalisé scalaire, et non arbitrairement à des observations multidimensionnelles. Les mesures manquantes ou rejetées déclenchent une prédiction sans correction. L'implémentation utilise la mise à jour de la matrice de covariance de Joseph et vérifie la semi-définition positive. L'incertitude augmente en cas de données manquantes, ce qui diffère d'une distance de rejet fixe en mètres.

Les résultats révèlent un excès de confiance

Diagram 1 · Use the button to switch views

La plage verticale est limitée à −3 à 4 m à des fins de comparaison. La valeur aberrante de 12 m et les erreurs importantes avec un faible R s'étendent au-delà. Examinez les fichiers CSV filtré, CSV à petit R et résumé complet pour vérifier la présence de valeurs non rognées. Les observations manquantes dans le CSV sont représentées par des valeurs NaN et ne sont pas comptabilisées comme des rejets.

Le filtre rejette deux observations, bien qu'une seule valeur aberrante artificielle ait été injectée. Un seuil ne permet pas d'identifier uniquement les échantillons délibérément corrompus. Une sous-estimation de la variance des observations d'un facteur 100 entraîne 237 rejets : la dispersion normale devient incompatible avec le niveau de confiance du filtre. Les corrections perdues empêchent alors de suivre la variation de vitesse, ce qui augmente encore les résidus. L'ajout d'un filtre seul ne garantit pas la robustesse.

Analyse de vos propres journaux

Vérifiez d'abord les horodatages, les unités et les images. Estimez la variance des observations sur les intervalles stationnaires, puis examinez les erreurs du modèle de mouvement et les résidus pendant le mouvement. Introduisez le filtrage en dernier, en indiquant le taux de données manquantes, le taux de rejet et l'incertitude, ainsi que l'erreur quadratique moyenne (RMSE). Le matériel réel nécessite une référence indépendante pour calculer des erreurs comparables.

Une faible valeur de Q peut lisser une courbe, mais augmente le décalage. Augmenter Q n'est pas systématiquement préférable. Ce classement s'applique uniquement à cette trajectoire, à cette initialisation, à ce bruit et à cette fenêtre d'évaluation.

Vérification complémentaire

La colonne sigma_m représente l'écart type estimé de la position. Tracez l'estimation plus ou moins deux fois sigma et comparez l'intervalle de données manquantes avec la variation de vitesse. Cela permet de distinguer l'erreur réelle de l'auto-évaluation de l'estimateur. Poursuivez avec le débogage de la fusion de capteurs pour diagnostiquer la synchronisation et la fonction de transfert avant de modifier la covariance.

Exécuter le même calcul localement

Extrayez le fichier ZIP source exécutable complet ou enregistrez le point d'entrée de l'expérience, les noyaux générés et le code source PID original dans un seul dossier. PID utilise la bibliothèque standard Python ; Kalman et ICP nécessitent NumPy. Les commandes originales de génération de figures restent disponibles.

python3 -m pip install numpy==1.26.4
python3 experiment.py kalman --parameters '{"q":0.4,"r":0.49,"gate":false}' > result.json

Activez/désactivez la porte pour comparer la réponse des valeurs aberrantes, puis définissez q à 0,002 ou r à 0,0049 et comparez l'erreur quadratique moyenne (RMSE) et le nombre de rejets. Les observations sont fixes d'une exécution à l'autre à partir du fichier CSV enregistré.

python3 reproduce.py --experiment kalman