Trajectographie passive sonar par estimation particulaire à maximum de vraisemblance

Dix-septième colloque GRETSI, Vannes, 13-17 septembre 1999 Trajectographie passive sonar par estimation particulaire à maximum de vraisemblance André Monin , Gérard Salut , Vincent Teulière Laboratoire d'Analyse et d'Architecture des Systèmes. Centre National de la Recherche Scientique. 7, avenue du Colonel Roche 31077 Toulouse Cedex 4 France [email protected],[email protected],[email protected] Résumé  Cet article introduit la technique particulaire à redistributions locales adaptatives, pour la résolution des problèmes d'estimation à maximum de vraisemblance. La performance de la technique est illustrée sur le problème de trajectographie passive SONAR Abstract  This paper introduces the particle method with local adapted redistributions, for the solution of maximum likelyhood estimation problems. The performance of this technique is illustrated by the SONAR target motion analysis problem. 1 Introduction. les approches de type Kalman étendu. Contrairement aux ltres exprimés dans la base cartésienne, le ltre exprimé La trajectographie passive SONAR consiste à estimer les dans la base polaire modié ne présente pas de problèmes éléments cinématiques d'une cible (distance, vitesse et cap) de stabilité et l'estimation est asymptotiquement non bià partir de mesures bruitées de ses émissions acoustiques. aisée. Dans l'approche particulaire, le repère choisi inue Le problème le plus classique est celui de la trajectogra- aussi sur la qualité de l'estimation. La modélisation pophie à partir de la mesure seule de l'azimut de la cible laire modiée permet de mieux répartir l'information sur ([1],[3],[4],[6]). Le porteur doit alors eectuer des manoeu- chacune des variables d'état et donc de réduire signicavres an d'assurer l'observabilité. L'adjonction des vari- tivement le nombre de particules nécessaire à une bonne ations Doppler au vecteur de mesure assure l'observabilité estimation soit du temps de calcul. De plus, l'algorithme du problème même lorsque le porteur reste immobile. Deux utilise une métrique locale pour ses redistributions qui approches sont utilisables pour traiter ce problème : le s'avère dénie par un calcul de type Kalman étendu. Le ltrage à minimum de variance ([5] et ses références) et choix du modèle polaire modié s'avère bénéciaire à ce l'estimation à maximum de vraisemblance ([2],[12]). Ces titre également. deux approches sont usuellement mises en oeuvre par apSoit un observateur (O) immobile et un mobile cible (C) proximations locales (ltre de Kalman étendu, algorithme observé dans le plan (voir gure 1 décrivant une trajectoire de type Newton) et posent généralement des problèmes rectiligne uniforme. d'initialisation et de comportement global. Nous introduisons ici l'outil de résolution global qu'est la technique particulaire déjà abordée avec succès [7] au sens du ltrage à minimum de variance et brevetée par ailleurs [9]. Cette technique s'avère ici supérieure en performance et en robustesse, au sens du maximum de vraisemblance, lorsqu'on l'assortit de redistributions probabilistes adaptées à la métrique de la vraisemblance. Elle se compare également favorablement aux approximations locales précitées, bien que l'application ne pose pas de dicultés multimodales. Nord K y Cible V Vy D Vx x A Observateur Ouest 2 Modélisation. 2.1 Modèle polaire modié Le modèle choisi est le modèle polaire modié [1]. En eet, cette modélisation a largement prouvé son ecacité dans Fig. 1: Localisation de la cible Le vecteur d'état du modèle polaire modié est consti- 705 Dix-septième colloque GRETSI, Vannes, 13-17 septembre 1999 tué des quatre composantes suivantes : 2 3 2 3 _ x1(t) A(t) 6 7 6 _ 7 7 = 6 D(t)=D(t) 7 X(t) = 64 xx2(t) (1) 5 4 A(t) 5 3 (t) x4(t) 1=D(t) Le modèle polaire modié se déduit du modèle cartésien (x; y; Vx ; Vy ). Il est possible d'exprimer l'état X(t) en fonction de l'état initial X(0) par : 2 3 x1 (0) 6 (tx1(0))2 + (1 + tx2 (0))2 77 6 X(t) = 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 tx21 (0) + (1 + tx2(0))2 (tx1(0))2 + (1 + tx2 (0))2 1 (0) x3(t0 ) + arctan 1 +txtx 2 (0) 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 (2) 2 2.2 Vecteur d'observation Le vecteur d'observation est constitué de l'azimut et de la fréquence perçue sous eet Doppler. Il s'exprime donc dans la base polaire modiée par : ( yA (t) = x3 (t) + vA (t) (3) yf (t) = x5(t)(1 , cxx24((tt)) ) + vf (t) où c est la célérité du son dans l'eau supposée constante, vA et vf sont deux bruits gaussiens indépendants, de moyenne nulle et de variances respectives RA et Rf . Y (t) = 3 Estimation particulaire à maximum de vraisemblance. 3.1 Principe de la technique Considérons le problème d'estimation paramétrique à maximum de vraisemblance non linéaire où la dynamique du système est décrite par :  xt = (x0; t) (4) yt = h(xt ) + vt où x est le vecteur d'état, y l'observation et v le bruit d'observation additif Gaussien de variance R. L'estimateur maximise la vraisemblance de x0 conditionnellement aux observations Y0T = fyk gk=0;:::;T T x0 = max c (5) x0 log pxjY (x0 jY0 ) = max x0 VT (x0) où VT (x0) = T (yk , h(xk ))T R,1 (yk , h(xk )) X k=0 PxjY (x0)jY0T )  N X i=1 (6) i x0 (7) i où x0 est la mesure de Dirac de support xi0 et i le poids de la particule xi0 i N X i = exp VTi = k=1 x4 (0) (tx1(0)) + (1 + tx2 (0))2 Une cinquième composante doit être rajoutée à l'état, c'est la pulsation propre acoustique de la cible y5 (t) = y5 (0) = F0. p L'approche particulaire (originellement esquissée dans [8]) consiste à approcher la mesure de probabilité conditionnelle pxjY (x0 jY0T ) par une somme nie et pondérée de mesures de Dirac (voir [10] et ses références) : exp VTk (8) où Vtk est la vraisemblance de la particule k conditionnelle aux mesures sur l'intervalle temporel [0; T] dénie en (6). L'application simple de cette approche bayésienne entraîne que le poids des particules dégénère avec le temps. Le poids des particules éloignées de la condition initiale à estimer devient très faible. Ces particules contribuent faiblement à la construction de l'estimateur : la loi des grands nombres est sous-utilisée. Pour un horizon de temps important tous le poids du support particulaire se concentre sur une seule particule. La densité de probabilité conditionnelle n'est alors plus représentée que par une seule particule. Pour pallier cette diculté, la méthode employée en ltrage particulaire consiste à redistribuer les particules évanescentes selon la densité de probabilité conditionnelle acquise dans tout le domaine. Chacune des positions xi0 du support reçoit un nombre de particules proportionnel à son poids i . Cette procédure de régularisation a pour eet de regrouper les particules dans les zones les plus probables de la densité de probabilité conditionnelle. 3.2 Redistribution diuse adaptative Dans la convergence asymptotique de l'estimation paramétrique, il n'existe pas de perturbations dynamiques. Au cours des redistributions successives, la particule la plus vraisemblable attire toutes les autres particules. L'ensemble du support nit par se concentrer sur la position la plus proche et l'estimation qui en découle hérite du biais initial. La gure 2 illustre ce problème pour un mobile en mouvement de rotation uniforme observé en cartésien. 1 vitesse angulaire vitesse estimee support particulaire 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0 10 20 30 iteration 40 50 60 Fig. 2: Estimation paramétrique - redistribution classique Il faudrait donc une distribution du support fxi0gi=1;:::;N inniment dense pour déterminer le maximum de la den706 Dix-septième colloque GRETSI, Vannes, 13-17 septembre 1999 N X probabilité conditionnelle pxjY (x0 jY0T ) avec exac4. Calcul de l'estimée c x0 = i xi0 Le nombre de particules étant nécessairement limi=1 sité de titude. ité, ce ne peut pas être le cas. La méthode de redistribution diuse adaptative consiste à redistribuer les particules non pas sur une seule et même position mais dans un voisinage autour d'une position donnée. La métrique de ce voisinage est obtenue naturellement à partir de la matrice d'information de Fisher F (xi0). L'inverse de F (xi0 ) se propage selon les mêmes équations que la matrice de covariance d'un lisseur de Kalman local autour de la trajectoire nominale [11] de chaque particule (la linéarité est conséquence de la localisation dans un voisinage). Soit une particule xi dont la dynamique est décrite par le système (4). Les équations de prédiction et de ltrage donnent : où 5. Évolution du support  xik+1 = f(xik )  k=k+1 6. Test de redistribution si le nombre de particules N aller en de poids inférieur à 101N est supérieur à 10 7, sinon retour en 2. 7. Redistribution  Calcul du nombre ni de particules à redistribuer autour de chaque position xi0 soit la partie entière de Ni .  Calcul pour chacune des positions retenues (de poids supérieur à 21N ) de P0ijk.  Redistribution de ni particules autour de la position xi0 par tirage selon N (xi0; P0ijk)  Réinitialisation de la fenêtre k = 0  Réinitialisation des vraisemblances Vki = 0 8. Retour en 2. Ptijt,1 = Fti,1Pti,1jt,1Fti,1T (9) T T Kti = Ptijt,1Hti [HtiPtijt,1Hti + R],1 (10) Ptijt = (I , Kti Hti )Ptijt,1 (11) t) i = @h(xt ) Fti = @(x et H (12) t @x(t) jx @xt x La covariance de la tâche de redistribution sur la position initiale P0ijt = F ,1 (xi0) est obtenue par l'équation de lissage : Pti,1 = Fti,1Pti Fti,1,T (13) i i initialisée par Pt = Ptjt. i t i t Poids initial 1/N dx dy Fig. 3: Principe de redistribution Le principe de redistribution est illustré par la gure 3. Autour de chacune des positions est redistribué un nombre de particules lles unitaires proportionnel au poids acquis par la particule mère. L'algorithme d'estimation particulaire se résume donc à : 1. Initialisation  Tirage du support fxi0gi=1;:::;N selon la loi a priori p(x0).  Initialisation de la fenêtre temporelle k = 0 2. Calcul de la vraisemblance Vki = Vki,1 , (yk , h(xik ))T R,1(yk , h(xik )) 3. Calcul des poids N X i i = exp Vk = exp Vkj j =1 3.3 Lisseur de Kalman conditionnel Lorsque une des composantes de l'état est observée linéairement conditionnellement aux autres composantes de l'état (c'est le cas de la fréquence propre du mobile F0 ), il est possible d'utiliser un estimateur exact sur cette composante conditionnellement aux autres . Le système étudié est donc de la forme : 8 < xt = (x0 ; t) t = Gt,1 + !t (14) : yt = h(xt )t + vt où  et h sont des fonctions non linéaires, G une fonction linéaire, !t et vt deux bruits gaussiens indépendants de moyenne nulle et de variance respectives Qt et Rt. Il sut donc de remarquer que conditionnellement à x, la partie dynamique/observation de  est linéaire,  est donc estimable de manière optimale par un lisseur de Kalman. Les équations du lisseur sont classiques. L'expression de la vraisemblance devient : T X (yt ,h(xit ))T (h(xit )Ptijt,1h(xit)T R),1(yt ,h(xit )) t=0 (15) où Ptijt,1 est la covariance de l'erreur de prédiction. L'estimée b est donnée par : VTi (x0) = d 0jt = N X i=1 i b0i jt (16) Cette méthode permet d'économiser un nombre signicatif de particules. La fréquence propre du mobile est estimée par cette méthode. 707 Dix-septième colloque GRETSI, Vannes, 13-17 septembre 1999 4 Application simulée mouvement rectiligne uniforme. Il est clair que ces résultats sont adaptables au problème ltrage avec bruits de Les seules informationsinitiales sont les mesures de l'azimut dynamique. Ainsi, il est possible d'appliquer cette techet de la fréquence. Aucune autre informationn'est disponible nique particulaire pour estimer la trajectoire d'un navire sur le cap, la vitesse et la distance. L'azimut et la fréquence cible eectuant des manoeuvres. sont donc distribués uniformément autour de leurs mesures initiales. La vitesse maximale est xée à 30 noeuds, cette valeur respecte la réalité physique des divers navires ex- Références istants. La distance est supposée comprise entre 0:5 et [1] Aidala and Hammel.  Utilisation of modied polar 128km. coordinates for bearing only-tracking . IEE TransA des ns comparatives, nous avons testé l'algorithme actions on automatic controls, 28(3):283293, 1983. sur la trajectoire dénie dans l'article [5]. La position initiale est (x(0) = 20km; y(0) = ,1:8km), la vitesse initiale [2] Jaufret and Bar-Shalom.  Track formation with bearing and frequency measuremnts in clutter est (Vx (0) = 0; Vy (0) = 9m=s), les écarts types de mesures . IEEE Transactions on aerospace and electronics pour l'azimut et la fréquence sont respectivement de 1:0 sytems, 26(6):9991010, Novembre 1990. degré et de 0:5 Hertz et la fréquence propre émise est de 300 Hertz. L'ellipsoïde de conance (99%) déduite de la [3] Kronham.  Bearing-only target motion analysis borne de Cramer Rao ainsi que le résultat de l'estimation based on a multihypothesis Kalman lter and adapparticulaire pour 50 trajectoires de bruits diérentes sont tative ownship motion control . IEE Proceedings représentés sur la gure 4. pour trois instants diérents on Radar Sonar and Navigation, 145(4):247252, 08 (5, 10 et 15mn). L'erreur d'estimation moyenne (R =R) 1998. passe au dessous du seuil de 10% au bout de 9mn pour [4] Nardone, Lindgren, and Gong.  Fundamental atteindre 5% en n de simulation. La gure 4 représente properties and performance of conventional bearingune trajectoire typique d'estimation. only target Target Motion Analysis . IEEE Transactions on Automatic control, pages 775,787, Sept 1984. [5] Passerieux, Pillon, Blanc-Benon, and Jaufret.  Target motion analysis with bearings and frequencies measurements . In Proceedings of the 22nd Asilomar Conference, Pacic Grove,CA, Nov 1988. [6] Peach.  Bearing-only tracking using a set of rangeparametrised extented Kalman lters . IEE Proceedings on Control theory Applications, 142(1):7380, 01 1995. [7] Rigal.  Filtrage non-linéaire,résolution particulaire et applications au traitement du signal . PhD thesis, Thèse de l'Université Paul Sabatier - LAAS/CNRS, Fig. 4: Estimation de la trajectoire de la cible. 1993. [8] G. SALUT.  Le thème non-linéaire en automatique et traitement du signal . Journées nationales du GdR traitement du signal et images du CNRS, septembre 1989. [9] G. Salut.  Procédé et système pour l'estimation optimale non-linéaire des processus dynamiques en temps réel . Brevet n INPI 94/07274, Europe extension n 95 95256.5-2206 USA n 8/750.313, juin 1994. [10] G. Salut.  Développement et applications du ltrage particulaire . Ecole des techniques avancées en Signal-Image-Parole de Grenoble Signal et non linéaire, août-septembre 1998. [11] James H. Taylor.  The Cramér-Rao estimation erFig. 5: Trajectoire d'un estimateur. ror lower bound computation for deterministic nonlinear systems . IEEE Transactions on automatic control, AC-24(2):343,345, avril 1979. 5 Conclusion [12] De Vlieger and Gmelig Meylig.  Maximum likelihood estimation for long-range target tracking using La nouvelle technique de redistribution diuse adaptative passive sonar measurements . IEEE Transactions a permise d'obtenir d'excellent résultats sur le problème on signal processing, 40(5):12161225, Mai 1992. de trajectographie SONAR où la cible est supposée en 8000 7000 t = 15mn 6000 5000 y (m) 4000 t = 10 mn 3000 2000 1000 t = 5mn 0 PORTEUR −1000 t=0 −2000 −0.5 0 0.5 1 1.5 2 x (m) 2.5 3 3.5 4 4.5 4 x 10 Position Vitesse 20 40 18 16 35 14 v (m/s) x (km) 12 30 10 8 25 6 4 20 2 15 −2 −1 0 1 2 3 y (km) 4 5 6 7 0 0 5 10 15 t (mn) 708