Brevet FR2949882A1 — publié le 20110311

METHODE DE TRAITEMENT DU SIGNAL COMPORTANT UNE ETAPE D'ESTIMATION DE LA MATRICE MEDIANE D'UN ECHANTILLON DE MATRICES La présente invention concerne une méthode de traitement du signal comportant une étape d'estimation d'une matrice médiane d'un échantillon de matrices. Les méthodes classiques de traitement du signal radar sont basées sur des techniques appelées TFAC pour "Taux de Fausse Alarme Constant". De manière générale, une méthode TFAC consiste à soustraire au signal une moyenne de l'ambiance du signal estimée sur une fenêtre glissante. Concrètement, la moyenne de l'ambiance est soustraite au logarithme de l'amplitude du signal en sortie de filtres FFT ("Fast Fourier Transform") disposés en banc, puis la différence est comparée à un seuil sur chaque voie FFT. Un OU logique permet de décider une détection si le seuil est dépassé sur l'une des voies FFT. Malheureusement, pour des formes d'onde Doppler avec des rafales courtes de moins de 16 impulsions par exemple, l'utilisation d'un banc de filtres FFT souffre plusieurs inconvénients, dont notamment une résolution faible. Tout d'abord, si le Doppler d'une cible est entre deux filtres, alors sa puissance est également répartie entre les deux filtres. De plus, un fouillis intense comme le fouillis de sol pollue l'ensemble des filtres à cause des lobes secondaires importants. Ainsi, avec un échantillon de 8 impulsions et en présence d'un fouillis de sol puissant, qui ne devrait théoriquement perturber que le filtre FFT correspondant à la vitesse nulle, l'ensemble des 8 filtres FFT peuvent être perturbés et la résolution Doppler de chaque filtre du banc est faible pour discriminer deux fréquences proches. L'invention se propose d'éviter certains des inconvénients précités des méthodes TFAC, qui trouvent essentiellement leur origine dans l'utilisation des filtres FFT. Elle a notamment pour but de maintenir une bonne résolution Doppler même pour des formes d'onde avec des rafales courtes. Pour cela, l'invention propose de moyenner de façon anisotrope l'information Doppler du signal. A cet effet, l'invention a pour objet une méthode de traitement du signal comportant une étape d'estimation d'une matrice X médiane d'un échantillon de N matrices {Bk }N., où N est un entier strictement supérieur à 1. Dans un mode de réalisation préféré, l'étape d'estimation de la 5 matrice médiane X peut comporter une étape d'estimation de la solution du problème d'optimisation: N N X = arg min Id(X, Bk) - 81 log(d (X, Bk) + 8) k=1 k=1 où d est une distance matricielle et 8 est un réel strictement positif, Avantageusement, la méthode peut permettre de détecter des 10 cibles avec un radar, les N matrices {Bk}k , pouvant être les matrices de covariance de séries d'un signal reçu par le radar. Par exemple, les séries du signal reçu par le radar peuvent être échantillonnées spatialement, de manière à permettre des traitements d'antenne, ou elles peuvent être échantillonnées spatio-temporellement, de 15 manière à permettre des traitements de type STAR Les séries du signal reçu par le radar peuvent également être échantillonnées temporellement, de sorte que la matrice médiane X fournit une estimation du spectre Doppler de l'ambiance. Avantageusement, si la distance matricielle géodésique entre la 20 matrice de covariance représentant le signal à tester et la matrice médiane X est supérieure à un seuil S, où S est un réel positif, alors une détection peut être décidée. Dans un autre mode de réalisation, l'étape d'estimation de la matrice médiane X peut comporter une étape d'estimation de la solution 25 d'une pluralité de problèmes d'optimisation définis itérativement par: Xm+1 = arg min , d 2 (X' Bk ) X k=1 d(Xm,Bk) où X0 est une matrice quelconque et d est une distance matricielle, la suite (X,,,)ä1EN convergeant vers X. Dans encore un autre mode de réalisation, chaque matrice Bk 30 étant une matrice carrée à n lignes et n colonnes où n est un entier strictement supérieur à 1, l'étape d'estimation de la matrice médiane X peut comporter une étape de paramétrisation des N matrices {Bk }k, conformément à un modèle autorégressif complexe, de sorte que à chaque matrice Bk correspond une paramétrisation e (^> = [Po ... j où po est la puissance du signal et p, à sont les coefficients de réflexion du modèle autorégressif complexe, une paramétrisation 0.e, médiane des N paramétrisations k("~ Ik , étant estimée par résolution du problème d'optimisation: \ 0med =argmin Ed(0,0k) où d est la distance riemannienne découlant de la métrique Kâhlérienne définie, dans le système de coordonnées o' =[po _Jr, par: 2 n , 2 dsn = JJ n.1 1 141,1 I\ Po _~ ~1ù~r12~ L'invention a encore pour principaux avantages qu'elle améliore la détection sur des fouillis inhomogène présentant des discontinuités abruptes du spectre Doppler, comme des cibles basses altitudes sur un fouillis de sol 15 ou comme des petites cibles sur un fouillis de mer. Elle améliore également la détection en présence de cibles proches ou de signaux parasites dans la zone d'estimation de l'ambiance fouillis. D'autres caractéristiques et avantages de l'invention apparaîtront 20 à l'aide de la description qui suit faite en regard de dessins annexés qui représentent : la figure 1, par un diagramme, un premier exemple de chaîne de traitement selon l'invention; les figures 2a, 2b, 2c et 2d, par un diagramme et trois graphes 25 respectivement, un deuxième exemple de chaîne de traitement selon l'invention; la figure 3, par un graphe, les performances obtenues en termes de détections et de fausses alarmes grâce à l'invention. 30 Les méthodes classiques dans la littérature permettent de calculer une matrice moyenne X au sens du barycentre d'un échantillon de matrices {Bk}k I. Il s'agit alors de résoudre le problème d'optimisation (1) qui suit, où d désigne la distance géodésiques entre deux matrices:10 N X =argminEd2(X,Bk) (1) k=1 En général, moyenne et médiane sont deux indicateurs statistiques différents. La moyenne reste fortement sensible aux données aberrantes, alors que la médiane est un estimateur statistique beaucoup plus robuste. Pour calculer une matrice médiane X d'un échantillon de matrices {Bk }k 1 , il s'agit de résoudre le problème d'optimisation (2) qui suit: N X = argminEd(X,Bk) (2) k=1 Une méthode itérative pour estimer la matrice X solution de (2) par descente de gradient consiste, à partir d'une matrice initiale Xo quelconque, à construire par itérations successives une suite de matrices (Xn),,EN qui 15 devrait converger vers X de la manière suivante: EN 10g(x;'"Bk Xn-ll2) Xn+1 = X1/2e k°11i1ogUGI /2Bkxnii2)l 1/2 Xn M avec lllog(xn-"2Bkxn-"211 = E log2 (2, ) V i=1 et det(Xn ù 2Bk ) = 0 où s est une valeur scalaire choisie arbitrairement faible pour assurer la 20 convergence, où 2411,..., /In } sont les valeurs propres de la matrice Xn 12BkXn 1 2 , et où la distance géodésique entre Xn et Bk est d(Xn,Bk) =I hog(Xn"ZBkXn-- I F avec 11.IIF la norme de Frobenius définie par VAL = JTrace(A.A+). 25 Malheureusement, pour des raisons de non différentiabilité aux points de données {Bk}k 1, la suite des itérés (Xn)nEN diverge aux voisinages des points de données {Bk }k 1. Donc, en pratique, la convergence n'est assurée que si les {Bk }k l "entourent" la matrice médiane. II s'agit là d'un inconvénient que l'invention se propose de surmonter, en s'inspirant notamment de l'algorithme de Weiszfeld et de sa généralisation proposée par Thomas Fletcher. En effet, l'algorithme de Weiszfeld permet de construire la médiane dans l'espace euclidien Rn et une généralisation proposée par Thomas Fletcher permet de l'appliquer à toute variété V de courbure sectionnelle positive ou nulle munie d'une mesure de probabilité. 10 Malheureusement, les éventuelles propriétés de convergence de l'algorithme de Weiszfeld dans les variétés générales sont inconnues et sujettes à caution. Qui plus est, pour assurer la convergence de la méthode de Fletcher, il est nécessaire de tester à chaque itération si l'on atteint un point du support de la mesure de probabilité. C'est pourquoi, dans la présente 15 demande, la méthode selon l'invention est appliquée à des variétés riemanniennes, de sorte qu'elle converge quel que soit le point initial et qu'il n'est pas nécessaire de tester l'appartenance des itérés au support de la mesure de probabilité. 20 Dans un mode de réalisation, l'invention propose avantageusement de remplacer le problème d'optimisation (2) qui précède par le problème d'optimisation régularisé (3) qui suit: N N X = arg min E d(X , Bk ) ù SE log(d(X, Bk ) + (5) (3) X k-1 où 8 est un réel positif strictement. L'invention propose alors la méthode itérative de base qui suit comportant deux boucles 1 et 2 et dans laquelle, quelle que soit une matrice 30 initiale X0, la suite de matrices (Xn )neN converge vers la matrice X solution de (3): k=1 25 Boucle 1: E \ Il N lOg(Xn1/2Bkx 1/2 _ I/2L' k-1 V Iog(Xn1 /2BkX;1/1 1/2 Xn+l = Xn ' Xn avec >5Q tant que log(tYn+12Xnxn+12) F Boucle 2: CN log(x, 2Bkx-u2) 1q 1/2 k-i10g (X"U2BkX I12 ~~ +Ô,, 1/2 Xn+1 = Xn e Xn avec Sn = n tant que llog(Xn+l2XnXn+12'11 F récemment par le module 1. Par la suite, cette matrice médiane sera appelée 30 "matrice médiane glissante", car elle évolue en permanence avec les nouvelles matrices de covariance fournies par le module 1. Un module 3 compare la matrice médiane glissante fournie par le module 2 à la matrice de covariance courante fournie par le module 1, qui représente le signal à tester à l'instant courant. Si la distance matricielle entre la matrice de covariance 6 courante Xcellule_._test et la matrice médiane Xmedlan est supérieure à un seuil S comme suit : d (Xmedian , X celluleùsousùtest) = g( 'iz '" 10 `Xmec/mnXcelluleùsousùtestXmedian t alors une détection est décidée. Ainsi, l'invention permet d'estimer une valeur approchée de la solution au problème d'optimisation régularisé (3) aussi précise que souhaité. Une variante de l'invention peut être de remplacer le problème (3) d'optimisation régularisé en norme LI qui précède par le problème (3bis) d'optimisation régularisé en norme L2 qui suit: N X min d2(X'Bk) m+l arg X k=l d(Xm,Bk) Dans le problème (3bis), la valeur obtenue à une itération permet de redéfinir une nouvelle fonction de coût en norme L2 qui, par descente de gradient, fournit la valeur optimale suivante. On construit ainsi une série de fonctions de coût en norme L2, chaque étape permettant de définir une fonction de coût par la valeur qui minimise la fonction de coût de l'étape précédente. On peut montrer que l'on tend ainsi vers la solution minimisant la fonction de coût en norme L1, c'est-à-dire la médiane. Il faut noter que cette solution n'est prouvée que dans le cas des variétés de courbures sectionnelles négatives ou nulles. (3bis)25 L'invention propose la méthode itérative qui suit dans laquelle, quelle que soit une matrice initiale X0, la suite de matrices (Xm)m{N converge vers la matrice X solution du problème (3bis): Initialisation: calcul de Xo,o Déterminer ko tel que f (Bko) = min{f(Bk0) :1 k 0 tel que: 10g~Bk 1/2~k Bk 1l2~ 0 D /2 2 X0,0 = Bko2e A I Iog(Bk~l BA8ko Ik Bko12 N venfie f (Xo,o) converge k=1 11og(X,,l/2BkXmli2 L A l'étape m: Faire Arrêt: Xm+i = Xm+l,n si _.ge N k=1 1 11og(Xm+11 2BkXm+12 )l }kN=., fournies successivement par le module 21, une paramétrisation médiane glissante emed . Un module 23 peut comparer la paramétrisation médiane glissante emed fournie par le module 22 à la paramétrisation courante fournie par le module 21, qui représente le signal à tester. Si la distance ARC entre la paramétrisation courante et la paramétrisation médiane glissante emed est supérieure à un seuil S, alors une détection est décidée. Dans ce qui suit, il est explicité comment le module 22 calcule la paramétrisation médiane glissante emed 15 selon l'invention. Pour définir une métrique dans cette nouvelle paramétrisation 0(0 = [P, ... unir, on utilise le biais de la géométrie affine duale introduite en géométrie de l'information, pour laquelle on définit une géométrie kâhlérienne via la fonction g qui suit, en prenant comme fonction potentielle 20 (Rn) l'entropie d'un modèle gaussien multivarié de moyenne nulle: a2,-ig" aHiaH; OU =ûlog(detR,, )ûnlog(ir.e)=E 0-1 (nûk).ln[lû1,-4121+ n.ln[z.e.Po] k=1 et H =ûR 25 En considérant le système de coordonnées [Po ,u, a métrique Kâhlérienne est alors la suivante: 2= Cd~J2 n -1 dsn n. ù +E(nûi Po La variété R+ x Dn-' est de courbure sectionnelle négative ou nulles, bornée par -1 (c'est la courbure sectionnelle constante sur D). La médiane s'obtient alors par résolution du problème d'optimisation (4) suivant: 30 N B,ä,d =argmiinEd(B,Bk) (4) k=1 où d est la distance riemannienne découlant de la métrique Kâhlérienne 5 définie précédemment. Supposant de plus que les points de données ne sont pas inclus dans une unique géodésique, la médiane est unique et s'obtient en appliquant à ce cas particulier la méthode itérative de base à deux boucles décrite 10 précédemment. Pour la simplicité des calculs, on résout le problème dans R+ x H"-' , où H désigne le demi-plan de Poincaré. Pour cela, on ramène les coefficients de réflexion dans H grâce à c-' défini par : C' :Dû*H 1+,u ,uHZ=i 1-,u La transformation C est connue sous le nom de transformation de Mobius. On utilise ici la transformation inverse C. En effet, les géodésiques sont alors des cercles centrés sur l'axe des abscisses comme illustré par la figure 2b, ce qui permet d'exprimer simplement la carte exponentielle. 25 Ainsi, les données étant (81,...9N) avec Bk = (P(k> äu, k), . äunkj )) _ (P(k) ,u(k) ) une méthode itérative pour résoudre le problème (4) selon l'invention peut alors comporter les 5 étapes suivantes, dans laquelle (z[iki)kEN converge vers ,uJmed pour tout 1 compris entre 1 et n-1: 30 Etape 1: Faire Pomed = exp(mediane(log(P(') ), log(P(2) ),..., log(Pr))) ; Etape 2: Initialiser p à 0; Etape 3: Pour k entre 1 et N Faire z(k) = C-' (p(k') Etape 4: Faire M = maxlz(l 15 20 10 11 Etape 5: Initialiser Z(°) = Ezl(0)...zn_,(0)r Pour r entre 1 et n-1 faire Début tant que Faire p-p+i Initialiser tzp) = 0 Pour k entre 1 et N faire ù Z(k)I2 (P)2 C(k) 1 1 ` 2 Re(z; k) _ Z1(P) ) z(p) ùC(k) (P) (k)\ 1 tz,k,l =ixslgne(Zl ùzl l(P)-Ci1+1/p t(ZP)l=t(ZP)1+t Z,k,l Fin pour Faire Im t(P) cz(p) =Re(Zl(P))+ (z,l) Im(Zl(P)) Re(tzP)) (P) = - arg(Zr(P) ù cZ (p> ) 15 Z1(P+') - = cz,p, (P) Zl ù cz,p) 1 çbl(P) ù 2Mk signe(Re(Z,(P) ù Cz )) exp(i 7 Z(P) ùcz ) Zl(P) ù Cz (p) +1/p 20 25 Tant que tzp) > 1/p et p La figure 2d illustre les coefficients de réflexion du modèle autorégressif complexe. Les ronds représentent les coefficients de réflexion bruts, les croix représentent les coefficients de réflexion moyens et les losanges représentent les coefficients de réflexion médians. La figure 3 illustre par un graphe l'amélioration des performances qui peut être obtenue grâce à l'invention. II s'agit d'un graphe qui permet de représenter, en échelle logarithmique, la probabilité de détection par rapport à la probabilité de fausse alarme. Une courbe CINV représente la courbe 10 obtenue en utilisant une chaîne de traitement selon l'invention. Une courbe CFFT représente la courbe obtenue en utilisant une chaîne de traitement classique à base de FFT. Les performances de détection sont nettement améliorées sur 15 fouillis inhomogènes présentant des discontinuités abruptes du spectre Doppler, comme un fouillis de mer ou un fouillis de sol, grâce au non-lissage des zones de discontinuités par l'estimateur d'ambiance, qui est une propriété du médian et du lissage anisotrope. De plus, la détection devient insensible à la présence de cibles 20 proches ou de signaux parasites dans la zone d'estimation de l'ambiance. L'invention décrite précédemment a encore pour principaux avantages qu'elle permet l'utilisation de formes d'onde doppler régulières non-uniformes, qui peuvent être entrelacées pulse à pulse. 25 Par ailleurs, il apparaît clairement à l'homme du métier que l'exemple d'implémentation décrit ci-dessus dans le cas particulier du traitement Doppler à partir des matrices de covariance des séries temporelles du signal reçu par le radar se transpose immédiatement et sans aucune difficulté au traitement d'antenne radar classique, à partir de matrices 30 de covariance de séries spatiales du signal reçu par le radar, ou encore aux traitements adaptatifs spatio-temporelles de type STAP ("Space-Time Adaptive Processing") à partir de matrices de covariance de séries spatio-temporelles du signal reçu par le radar ou encore aux techniques de retournement de type DORT ("Decomposition of the time reversai operator"). 35 II apparaîtra également clairement à l'homme du métier que la méthode décrite ci-dessus pour estimer la matrice médiane d'un échantillon de matrices peut trouver des applications dans de nombreux autres domaines, comme dans le domaine de l'imagerie médicale pour estimer la matrice de covariance de matrices de rotation, mais également dans le domaine du "data mining", de l'analyse ordinale des données ou encore de la cryptographie quantique. REVENDICATIONS1. Procédé de traitement d'un signal reçu par un radar pour détecter des cibles, caractérisé en ce qu'il comporte une étape d'estimation d'une matrice X médiane d'un échantillon de N matrices de covariance {Bk}k l de séries du signal reçu par le radar où N est un entier strictement 5 supérieur à 1, l'étape d'estimation de la matrice médiane X comportant une étape d'estimation de la solution du problème d'optimisation: N N X = argmin1d(X, Bk ) ù gllog(d(X, Bk ) + 5) k=1 k=1 où d est une distance matricielle et Ô est un réel strictement positif. 10 2. Procédé selon la revendication 1, caractérisé en ce que les séries du signal reçu par le radar sont échantillonnées spatialement, de manière à permettre des traitements d'antenne, ou en ce que les séries du signal reçu par le radar sont échantillonnées spatio-temporellement, de manière à permettre des traitements de type STAP. 15 3. Procédé selon la revendication 1, caractérisé en ce que les séries du signal reçu par le radar sont échantillonnées temporellement, de sorte que la matrice médiane X fournit une estimation du spectre Doppler de l'ambiance. 4. Procédé selon la revendication 2 ou 3, caractérisé en ce que, si la distance matricielle géodésique entre la matrice de covariance représentant le signal à tester et la matrice médiane X est supérieure à un seuil S, où S est un réel positif, alors une détection est décidée. 25 5. Procédé selon la revendication 1, caractérisé en ce que l'étape d'estimation de la matrice médiane X comporte une étape d'estimation de la solution d'une pluralité de problèmes d'optimisation définis itérativement par: N d2(X,Bk) 30 X m+1 = arg min 1 k=1 d(Xm,Bk) 20où X0 est une matrice quelconque et d est une distance matricielle, la suite (X,n)ä,N convergeant vers X. 6. Procédé selon la revendication 1, caractérisé en ce que, chaque matrice 5 Bk étant une matrice carrée à n lignes et n colonnes où n est un entier strictement supérieur à 1, l'étape d'estimation de la matrice médiane X comporte une étape de paramétrisation des N matrices {Bk}k 1 conformément à un modèle autorégressif complexe, de sorte que à chaque matrice Bk correspond une paramétrisation e(n) _ [po u, ~n ,lT 10 où po est la puissance du signal et à sont les coefficients de réflexion du modèle autorégressif complexe, une paramétrisation e,,,ed médiane des N paramétrisations {Bk(�}k I étant estimée par résolution du problème d'optimisation: emed = arg min d (6, ek k=1 15 où d est la distance riemannienne découlant de la métrique Kâhlérienne définie, dans le système de coordonnées 0(n) = [po , par: (dg" 141,1 2 dsï =n. ù° +1(nùi) z Po / ,=, 1ùICI 7. Dispositif de traitement d'un signal reçu par un radar pour détecter des 20 cibles, caractérisé en ce qu'il comporte: - un module (1) de construction de matrices de covariance B de séries du signal reçu par le radar; - un module (2) d'estimation d'une matrice médiane X des N matrices de covariance successives {Bk}k 1, où N est un entier 25 strictement supérieur à 1, fournies le plus récemment par le module (1) de construction, par estimation de la solution du problème d'optimisation: N N X =argmin1d(X,Bk)ùgElog(d(X,Bk)+ k=1 k=1 30 où d est une distance matricielle et est un réel strictement positif;un module (3) de comparaison de la matrice médiane X fournie par le module d'estimation (2) à la matrice de covariance courante fournie par le module de contruction (1). 8. Dispositif selon la revendication 7, caractérisée en ce que les séries du signal reçu par le radar sont échantillonnées spatialement, de manière à permettre des traitements d'antenne, ou en ce que les séries du signal reçu par le radar sont échantillonnées spatio-temporellement, de manière à permettre des traitements de type STAP. 9. Dispositif selon la revendication 7, caractérisée en ce que les séries du signal reçu par le radar sont échantillonnées temporellement, de sorte que la matrice médiane X fournit une estimation du spectre Doppler de l'ambiance. 10. Dispositif selon la revendication 8 ou 9, caractérisée en ce que, si la distance matricielle géodésique entre la matrice de covariance représentant le signal à tester et la matrice médiane X est supérieure à un seuil S, où S est un réel positif, alors une détection est décidée. 20 11. Dispositif selon la revendication 7, caractérisée en ce que l'estimation de la matrice médiane X inclut d'estimer la solution d'une pluralité de problèmes d'optimisation définis itérativement par: L d2(X,Bk) Xm+1ù arg min k=1 d(Xm,Bk) 25 où Xo est une matrice quelconque et d est une distance matricielle, la suite (X�,),,,EN convergeant vers X. 12. Dispositif selon la revendication 7, caractérisée en ce que, chaque matrice Bk étant une matrice carrée à n lignes et n colonnes où n est un 30 entier strictement supérieur à 1, l'estimation de la matrice médiane X inclut une paramétrisation des N matrices {Bk}k l conformément à un modèle autorégressif complexe, de sorte que à chaque matrice Bk correspond une paramétrisation g('o =[po ,1, .. /jn_,y où po est la puissance154 du signal et à Dry, sont les coefficients de réflexion du modèle autorégressif complexe, une paramétrisation emed médiane des N paramétrisations 1Bk(n)Jk i étant estimée par résolution du problème d'optimisation: B,ned = arg min E d (B, Bk ) 9 k_1 où d est la distance riemannienne découlant de la métrique Kâhlérienne définie, dans le système de coordonnées g(^) =[po un_1y, par: ds dP ll2 n_1 Id `Iz ,1, =n1--°) +Ë(nùi)(1ùf4z o =1