SEIZIÈME COLLOQUE GRETSI — 15-19 SEPTEMBRE 1997 — GRENOBLE 1121

Filtrage non linéaire en poursuite. Une approche par mélanges

gaussiens

Serge Prosperi
Département ES2I
École Supérieure d’Ingénieurs de Luminy
Luminy case 925,
13288 Marseille cedex 09

RÉSUMÉ ABSTRACT

Nous proposons de modéliser les bruits du modèle de filtrage par des We propose here to use a gaussian mixture model for the plant and
mélanges gaussiens, ce qui a été étudié dans [6], où des the measurement filter noise, which has already been considered in
approximations par mélanges gaussiens des lois de probabilité [6], where approximations with gaussian mixtures of the common
scalaires usuelles sont présentées. Dans le cas d’équations linéaires - scalar probability laws are presented. If the equations are linear, or
ou linéarisables - on en déduit des variantes des filtres de poursuite can be linearized, one obtain modified counterparts of the PDA
(du type PDA) intégrant les composantes du mélange. Nous nous Filter which accounts the components parameters of the mixture. We
intéressons ensuite à la densité exacte du vecteur d’état dont on then consider the exact conditional probability density of the state
observe qu’elle est elle - même un mélange gaussien, avec un calcul vector which is itself a gaussian mixture, with a recursive
récursif de ses composantes, et une réduction de leur nombre computation of its components, and a reduction of their number
conduisant à un algorithme approché. Ces résultats sont appliqués au which leads to a rough estimate algorithm. These results are applied
probléme de croisement de sources en sonar passif. Nous abordons to the crossing of passive sonar sources. We end the paper with an
pour finir le cas d’équations non linéaires, où les hypothèses des insight to the nonlinear equations case, for which model the classical
modèles classiques ne sont pas respectées. assumptions are not true.

d’origine ce n’est plus un modèle de filtrage linéaire (c’est un

problème exponentiel) :
1 Cadre du problème de poursuite
Xˆ = E ( X / Z , u ≤ t ) =
t t u ∑P ( H ) E ( X / Z ) (1.5)
t
H
Dans la théorie du filtrage linéaire avec variables d’état, où l’hypothèse H désigne une branche du graphe
les équations du modèle sont (X désignant le vecteur d’état et d’affectation. Les solutions (approchées) usuelles calculent
Y le vecteur de mesures) : l’état estimé de manière récursive - sans mémorisation des
X' (t) = A(t ) X(t ) + Γ(t) W(t) (1.1) mesures, dans les contextes mono ou multicible, soit en
Y (t ) = C(t ) X(t) + V(t ) (1.2) associant le plus proche voisin [5], soit en estimant le vecteur
V et W sont des bruits blancs gaussiens centrés, et les d’état comme une moyenne pondérée des estimations faites à
matrices A, Γ, C, ainsi que les covariances Q et R des partir des différentes mesures (Probabilistic Data Association,
processus W et V sont connues. La solution optimale [1]). Il s’agit dans tous les cas de méthodes sous optimales. Il
(l’estimateur bayésien du vecteur d’état) est fournie par le faut en outre disposer au préalable dans un contexte
filtre de Kalman. On peut le généraliser, par linéarisation des multicible d’une procédure d’initialisation des algorithmes de
équations, au modèle non linéaire : filtrage (détection - estimation des sources), fondée sur un test
X' (t) = f [t, X(t )] + Γ(t) W (t ) (1.3) de vraisemblance généralisée [3].
Y (t ) = h [t, X(t)] + V(t ) (1.4)
Pour obtenir le Filtre de Kalman Etendu. Les bruits sont 2 Cas de bruits non gaussiens
toujours dans ce cas supposés gaussiens, d’où une loi
conditionnelle gaussienne, caractérisée par ses deux premiers Dans le pistage de sources passives (sonar), un modèle
moments. En poursuite multipistes avec fausse alarme, On mélange gaussien permet de modéliser le phénomène de
observe simultanément plusieurs processus, X 1 ... X p, et croisement de deux sources, avec la présence de deux ou trois
plusieurs mesuresY 1 ... Y q où les valeurs Y j peuvent être des
modes. C’est l’objet des essais présentés au paragraphe 5.
mesures effectives effectuées sur les processus observés - on
ignore a priori lesquels, ou des fausses alarmes. L’observation De fait, toute loi continue peut être approchée avec une
est un vecteur Z ( t) = ( Y 1, ..., Y q) dont les composantes sont précision arbitraire par un mélange gaussien [2]. Cette
des candidats de mesures des processus en t . Du fait de la propriété n’a d’intérêt que si le nombre de termes est restreint,
non connaissance de l’affectation des mesures aux processus ce qui est souvent le cas en pratique, rendant l’approximation
1122
très intéressante par rapport aux autres approximations, issues
de développements de Taylor (Développements de Edgeworth
Qk = cov ( Vk ) = ∑ pk ,l cov ( Vk ,l ) (3.4)
k ,l

∑
et Cram Charlier). t
Si W(t) suit une loi mélange gaussienne ayant pour densité Pk / k = Gk Λk Gk + γ k ,l, r,0 Pk / k −1 +
f
W (t )
(x) = ∑ p f ( x) , ft , l étant une densité
t, l t,l
l,r
l
gaussienne. On indexe ( L) les applications t → l (t) .
(1 − ∑ γ k ,l, r,0 ) ( I − Gk Ck ) Pk / k −1 (3.5)
l,r

Xˆ = E (X / Z , u ≤ t)
t t u t −1
G = Pk / k − 1 Ck Sk est le gain du filtre, avec
= ∑ PL P (H / L) E ( Xt / Z, H, L) (2.1) k

H, L S
k
= C k Pk / k −1
C
k
t
+ ∑ q
k, r
R
k,r
(3.6)
r

3 Filtrage Λk = ∑ t
γ k , l, r, j (ν k, l , r, j ) (ν k, l , r, j ) − ν k ν k
t
l ,r , j
Nous proposons deux approches du problème du filtrage,
où νk, l, r, désignent les innovations et νk les innovations
dans le cadre d’un modèle discret, l’une de type récursif, j
dérivée du modèle de filtrage linéaire classique, la seconde pondérées.
fondée sur le calcul direct de la densité de probabilité du
vecteur d’état. 3.2 Calcul direct de la loi de l’estimateur

3.1 Filtrage récursif Cette approche consiste à traiter directement la densité de

probabilité conditionnelle :
On peut donner une version séquentielle de l’algorithme de
poursuite (Analogue du PDA, que nous notons MPDA), sous
p (X / Z) = ∑ PL P (H / L) p ( X / Z, H, L) (3.7)
H, L
la forme d’un filtre non linéaire qui consiste à approcher la loi
a posteriori (loi mélange) du vecteur d’état par une loi Cette écriture peut là encore être faite de manière récursive
gaussienne de même espérance et de même covariance. On dans le cas d’un modèle (linéaire) discret avec les équations
peut alors proposer une écriture récursive : du filtrage non linéaire.
Si on suppose que la densité conditionnelle de l’état à l’étape

∑
k - 1 est un mélange gaussien :
Xˆ (t ) = β l, j Xˆ l, j (t )
∑
p (3.1) k −1
t, l
p ( xk −1
/Z ) = PL φ L ( x k − 1 ) (3.8)
l, j
L
β l , j = P ( H j / l, Z) , probabilité d’association de la mesure k
à l’étape suivante la densité conditionnelle p ( xk / Z ) s’écrit
j dans l’hypothèse d’un bruit gaussien de densité f est
∑ ∑
t,l
calculé sans difficulté par la règle de Bayes en supposant la α P
L
p
k, l
q
k,r
β k, j gk , r (yk , j − Ck x k )
probabilité de détection des sources connue et que la L l, r, j
répartition des fausses alarmes suit une loi de poisson [1].
Xˆ l , j (t) est l’estimateur obtenu par les équations de Kalman ∫ f k,l ( xk − Ak x k −1 ) φ L ( x k −1 ) dxk −1

toujours dans l’hypothèse où Y = Y j et d’un bruit d’état où α est une constante de normalisation. Il s’agit bien encore
gaussien de densité f . d’un mélange gaussien. On peut calculer les densités de
t,l
probabilité qui apparaissent à l’intérieur de cette somme.
On peut traiter de manière analogue le cas où le bruit de
Mais le résultat de ce calcul est connu, puisqu’il s’agit des
mesure est un mélange gaussien, de proportions q et de
t,r densités gaussiennes dont les moments sont fournis par le
densités gt , r . L’estimateur a la même forme que filtre de Kalman.
Cette technique conduit donc ici encore à mettre en parallèle
précédemment, dans les expressions (2.1) et (3.1) :
un banc de filtres. Comme dans tout problème de filtrage non

∑
linéaire, on est amené à estimer la densité de probabilité, c’est
Xˆ (t ) = p
t, l
p
t, r
β l, r, j Xˆ l, r, j (t ) (3.2) à dire une fonction multivariable, mais représentée ici par un
l, r , j nombre fini de paramètres. L’inconvénient est que ce nombre
La covariance du filtre MPDA peut être calculée croît - rapidement - si on ne limite pas le nombre de termes du
immédiatement à partir de (3.2), dans le cas discret : mélange. Cela paraît difficile à mettre en oeuvre en dimension
On pose γ k ,l, r, j = pk ,l qk ,r β k,l ,r, j Le produit des élevée, ou avec des bruits nécessitant un grand nombre de
composantes pour le mélange gaussien, mais peut être
probabilités des composantes pour les bruits d’état et de
proposée en écoute passive (par exemple en sonar) avec
mesure et de la probabilité d’association.
mesures angulaires.
tt
Dans tous les cas de figure, il faut néanmoins réduire en
Pk / k −1 = Ak Pk −1 / k −1 Ak + Γk cov (Vk ) Γk (3.3) continu le nombre de composantes du mélange.
 1123
On peut également utiliser un estimateur Markovien à un
ordre s quelconque, de manière analogue à ce que propose 5 Croisement de sources (sonar passif)
s
Singer dans [8], ce qui requiert le calcul de m filtres, si l’état
estimé est un mélange comportant m composantes
gaussiennes. Cela peut être vu soit comme une limitation de 5.1 Le problème
la propagation des équations (3.8) à un nombre fini s de pas,
Le contexte est celui de deux sources donnant lieu à des
ou comme une écriture du filtre MPDA à l’ordre s.
mesures angulaires (azimut) qui correspondent à des pics
L’hypothèse faite ici est dans tous les cas (contrairement à
d’énergie dans le domaine spatial, représentées dans le cas de
[8]) que la loi a priori de l’état est une loi mélange.
mouvements lents par une équation d’état approchée. On
suppose que l’accélération est nulle, ce qui est le cas en
4 Modélisation dehors du point de contact entre les sources de bruit et
l’observateur. Lorsque les sources se rapprochent l’une de
l’autre, les mesures sont entachées d’un biais, jusqu’à un
4.1 Choix des bruits de modèle point où elles ne sont plus résolues - on peut noter que ce
problème de capture a déjà été abordé (JPDAM, [1]), mais
Le choix du mélange (nombre de composantes et paramètres)
dans le contexte d’échos radar, avec une discrétisation en
peut être fait par apprentissage (identification de mélange
cases distance. Le biais moyen est estimé par une phase
gaussien) si on dispose d’un nombre suffisant de mesures
d’apprentissage (il dépend des niveaux de puissance
pour chaque source, ce qui n’est pas toujours faisable pour
rayonnée par les sources). Ici, le bruit de mesure n’est plus
des systèmes dynamiques non stationnaires. Si les bruits sont
gaussien. En outre, dans la phase de transition entre les zones
stationnaires, cette loi peut être évaluée paramétriquement, en
résolue (avec présence de deux mesures) en non résolue (une
fonction des caractéristiques des processus, dans la mesure où
seule mesure), chaque mesure peut être considérée comme
on peut les identifier (c’est à dire s’affranchir à cette étape du
suivant un mélange à deux composantes.
problème d’affectation). Cela ne posera pas de gros
problèmes en présence de forts rapports signal sur bruit.
Dans le cas contraire, ce problème est analogue au problème
5.2 Résultats
classique de normalisation. Les essais ont été réalisés jusqu’ici uniquement sur signaux
La question de l’identification de lois mélange a reçu un synthétiques, en générant des sorties de voies, avec un critère
grand intérêt depuis K. Pearson, et a instigué de nombreuses de détection par seuil qui fournit des mesures correspondant à
méthodes (méthode des moments, approximation stochastique
une probabilité de détection PD (estimée) ainsi que des fausses
ou encore techniques issues des problèmes de classification
alarmes. Nous avons pris pour vecteur d’état un vecteur à 4
comme les nuées dynamiques, algorithme EM) car les
composantes qui comprend les positions et vitesses des deux
équations de vraisemblance sont trop complexes pour être
sources, et identifié la loi du bruit en azimut dans la phase de
résolues directement. Il s’agit d’un problème difficile, en
croisement en fonction de la distance estimée des sources et
particulier pour l’estimation du nombre de composantes. Si la
de leur puissance relative, par un mélange gaussien à 4
dimension des données est faible, et le nombre de modes
composantes, puis filtré les pistes par un filtre de type MPDA,
connu (ou borné), le choix de la méthode n’est pas
et comparé avec l’utilisation d’une technique NNA et d’un
fondamental, et une approche de type vraisemblance est ici
filtrage PDA conjoint des deux sources (JPDA). Les filtres
bien adaptée.
sont initialisés avec les valeurs exactes.
Le scénario retenu est celui de deux sources de même
4.2 Réduction du nombre de composantes puissance, en l’absence de manoeuvre pendant le croisement.
Pour la détermination des composantes retenues dans la Contrairement à [3] nous ne tenons pas compte de l’énergie
propagation du calcul de la densité du vecteur d’état, on peut des sources, et le réglage du seuil fournit simplement le taux
chercher la meilleure approximation de cette loi sous la forme de fausse alarme (un rapport signal sur bruit de 10 dB assure
d’un nouveau mélange ayant un nombre N de composantes une PD supérieure à 0.9 dans tous les cas). L’écart - type des
déterminé. mesures hors croisement est de 1 degré.
Mais ici on connaît une expression analytique de la loi Nous représentons à la figure 5.1 l’erreur d’estimation en
(mélange gaussien) et on veut réduire le nombre de fonction de la densité µ de fausses alarmes exprimée en
composantes. On peut pour cela proposer une méthode directe nombre moyen d’alarmes par degré, observée sur toute la
d’estimation des paramètres de la nouvelle loi, en minimisant phase de croisement avant capture. le modèle exact de bruit
une distance entre les deux fonctionnelles, comme l’écart est supposé connu.
quadratique - ce qui revient à la même complexité que le La figure 5.2 présente également l’écart type de l’erreur, mais
problème d’identification, ou utiliser une méthode ad hoc (de en l’absence de fausses alarmes, en fonction de la distance d
type classification) ce que nous avons mis en oeuvre. des deux pistes exprimée en degrés, ce qui permet de juger de
1. On néglige les proportions inférieures à un seuil, fonction manière intrinsèque la discrimination des sources. Le modèle
du nombre de modes retenu. pour le bruit de mesure a été déterminé une seule fois pour la
2. On établit des classes de composantes sur un critère de génération des mesures et repris par la suite, et est ici estimé à
distance normalisée des moyennes, et on traite ces classes de chaque tirage.
façon indépendante. Ces valeurs ont été obtenues en simulation de type Monte
Carlo sur 100 tirages. Nous n’avons pas comparé les résultats
à la borne d’estimation en poursuite [4] car il était difficile de posteriori - comme le caractère gaussien est perdu par
tenir compte dans la matrice de Fisher du modèle exact pour transformation non linéaire. Il faut dans ce cas pouvoir
les lois mélange, mais cette comparaison est certainement linéariser les équations du filtrage. De fait, on peut donner
souhaitable pour avoir un étalon absolu des performances de une version “étendue“ de l’estimation de la densité de
la poursuite. probabilité du vecteur d’état conditionnellement aux mesures
dans le cas d’un modèle linéarisable, avec des versions à
l’ordre 1 ou 2 du filtre de Kalman Etendu, aussi bien pour le
filtrage classique que pour le problème de poursuite. Dans le
7
cas le plus général (équations 6.1 à 6.5) on n’a pas d’autres
6
NNA ressources, sauf dans l’hypothèse de petites perturbations, que
5
4 de discrétiser l’espace des états, ce qui n’est possible que si le
JPDA
3 vecteur d’état est de faible dimension. Nous reviendrons
P
D MPDA ultérieurement sur ces deux points, qui font l’objet d’une
2
1 étude spécifique.
0
0 0,2 0,4 0,6 0,8 1
7 Conclusion
. Fig. 5.1 Ecart - type σ = f (µ) Lors du croisement de sources en sonar passif, l’hypothèse
de bruits gaussiens n’est pas respectée, et nous proposons ici
de modéliser les bruits du modèle de filtrage par des
mélanges gaussiens. Dans le cas d’un modèle linéaire, on peut
2 écrire de manière analytique les équations du filtre optimal,
c’est à dire la densité de probabilité conditionnelle de l’état,
1,5
qui est un mélange gaussien, et utiliser un estimateur
1 NNA markovien d’ordre r (pour r = 1, on retrouve les équations du
JPDA PDA modifiées), ou un calcul récursif de la densité en
0,5 limitant le nombre de termes du mélange. Nous avons testé
MPDA
l’utilisation du filtre d’ordre 1 (MPDA) sur des signaux
0
0 0,5 1 1,5 2 2,5
synthétiques pour un croisement de sources, et vérifié l’intérêt
de cette approche, en comparant avec l’utilisation des filtres
classiques en poursuite. Les performances ne semblent pas
. Fig. 5.2 Ecart - type σ = f (d) être très sensibles à la qualité de l’estimation des paramètres
du mélange. Compte tenu de l’évolution des puissances de
calcul depuis 25 ans, l’approche par modèles de lois mélanges
6 Modèles non linéaires déjà proposée dans [6] peut retrouver aujourd’hui un regain
d’intérêt.
Dans le modèle général de filtrage non linéaire, les équations
de transition d’état et de mesure peuvent s’écrire:
X’ (t) = Φ [ X (t) ] (6.1)
8 Bibliographie
Y (t) = ψ [ X (t) ] (6.2) [1] Y. Bar Shalom, T. Fortmann : “Tracking and Data
Φ et ψ étant des fonctions aléatoires. L’équation (6.1) n’a de Association” Academic Press, 1988.
sens que si on peut intégrer le second membre (intégrale [2] W. Feller : “An Introduction to Probability Theory and its
stochastique), par exemple si X k+1 = Φ [k, Xk, W k+1], où Φ Applications”, T2, John Wiley, 1971.
est une fonction inversible par rapport à sa troisième variable, [3] S. Prosperi : “New Passive Tracking Algorithms using
ce qui permet la récursion de la loi p (Xk / Y1,..,Yk). Power Information with application to Passive Sonar”,
On peut également formuler le problème général de poursuite ICASSP, Mai 1991, Toronto.
avec un modèle non linéaire : [4] S. Prosperi : “Bornes d’estimation en poursuite de
Xj’ (t) = Φ [ Xj (t) ] (6.3) pistes”, GRETSI , Septembre 1993, Juan les pins.
Yj (t) = ψ [ Xj (t) ] (6.4) [5] R. A. Singer, J. J. St
Z (t) = ζ [Y1 (t), ..., Yp (t)] (6.5)
La fonction aléatoire ζ = γr ° χs ° σ est la composée d’une
permutation, d’une troncation d’ordre s, et d’une adjonction
de r valeurs aléatoires. Si les opérateurs Φ et ψ sont
linéaires, le problème est linéaire (mais il ne s’agit pas ici
d’un problème de filtrage linéaire).
Dans le cadre général du filtrage non linéaire on n’a pas
d’autre choix que d’estimer directement la densité de
probabilité conditionnelle de l’état.
On peut encore proposer de modéliser les bruits par des lois
mélange, mais ce caractère n’est plus respecté pour la loi a