Brevet FR2955950A1 — publié le 20110805

PROCEDE ET DISPOSITIF PERMETTANT DE DETERMINER UN PARAMETRE CENTRAL ASSOCIE A DES COEFFICIENTS D'IDENTIFICATION OBTENUS PAR UN ALGORITHME DE BURG La présente invention concerne un procédé exécuté notamment au sein d'un récepteur comprenant un ou plusieurs capteurs de réception d'un signal complexe, ledit procédé permettant de différencier au sein dudit signal complexe des informations utiles à un traitement de signal en amont du procédé d'informations considérées comme non utiles. 1 o Dans le cas d'une application radar à impulsion, le signal reçu par le récepteur après démodulation complexe produit un signal échantillonné sur les voies sinus et cosinus, qui constituent une suite d'échantillons complexes, chaque échantillon correspond à un petit domaine de la portée radar encore appelée case distance. 15 Dans le cas d'un procédé de lissage d'un signal, l'idée sera d'éliminer les informations utiles détectées du bruit majoritaire. Le procédé et le système associé selon l'invention s'appliquent par exemple pour des radars Doppler. Le procédé selon l'invention est applicable de manière plus générale 20 aux systèmes temps-réels embarqués du fait du peu de calcul nécessaire. Dans le cas d'une application radar, l'objectif premier d'un radar est de pouvoir détecter des cibles dans une plage locale étendue dans la dimension distance radiale radar-cible et dans la dimension angulaire. Cette plage est 25 découpée en cases distance qui contiennent chacune un nombre N de coups au but (qui peut être très limité), c'est à dire de mesures complexes (YI, Y2, YN) issues du récepteur. Dans la suite de la description, ces N mesures sont appelées échantillon. Ces mesures représentatives de l'écho contiennent à la fois 30 l'information relative aux cibles et celle relative au bruit d'environnement ou fouillis. On isole classiquement les cases contenant les cibles par un seuil adapté au fouillis local (TFAC). 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 de transformée de Fourier ou 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 peut être perturbé et la résolution Doppler de chaque filtre du banc est faible pour discriminer deux fréquences proches. Un autre objectif d'un radar, à des fins par exemple de cartographie, est de distinguer et de séparer des zones localement homogènes dans un fouillis 25 globalement inhomogène. C'est particulièrement le cas en radar doppler météo où il s'agit d'analyser des perturbations atmosphériques et de distinguer en leur sein des zones homogènes de caractéristiques spectrales différentes. C'est également le problème dans le cas d'application de surveillance 30 maritime où l'on est amené à séparer des zones de fouillis homogènes dans 3 2955950 un ensemble inhomogène et où l'on doit détecter des cibles faibles dans le fouillis local. Le problème général est donc de déterminer la caractéristique moyenne d'une collection de P échantillons constitués chacun de N mesures 5 complexes, soit pour isoler un échantillon remarquable, soit pour éliminer un élément perturbateur, soit pour classer l'ensemble en lots distincts. Ceci peut être appliqué dans d'autres domaines que le traitement radar, dans tous les cas où les mesures se présentent sous la forme d'un signal complexe, qui peut aussi être un signal réel. 10 Ce problème est abordé depuis déjà plusieurs années dans le monde du traitement du signal, tant pour les applications radars que pour toutes sortes d'applications similaires. L'une des solutions explorées à ce jour porte sur la recherche de médiane dans les espaces de matrices de corrélation c'est à dire les matrices complexes symétriques définies positives, placée 15 dans le cadre de la théorie géométrique de l'information. Les méthodes ainsi proposées font appel à la diagonalisation, à l'extraction de racines carrées, au produit et à l'exponentiation des matrices carrées avec les problèmes de conditionnement que cela pose, itérés sur tous les éléments de la collection d'échantillons dont on cherche la médiane et itérés sur le nombre de pas 20 nécessaire à la convergence de l'algorithme. De plus les matrices de corrélation empiriques sont toujours entachées d'une forte incertitude sur de petits échantillons de mesures. Ces algorithmes sont difficilement applicables à des systèmes temps-réel embarqués dans l'état actuel de la technologie des calculateurs. L'invention se propose notamment d'éviter certains des inconvénients précités des méthodes TFAC, qui trouvent essentiellement leur origine dans l'utilisation des filtres FFT. L'objet de la présente invention est donc de déterminer la caractéristique moyenne ou médiane d'une collection de P échantillons constitués chacun de N mesures complexes, soit pour isoler un échantillon 4 remarquable, soit pout éliminer un élément perturbateur, soit pour classer l'ensemble en lots distincts. La présente invention propose une méthode simple de calcul de la médiane ou de la moyenne, et de façon plus générale d'un élément central représentatif des échantillons mesurés, par un algorithme du gradient à partir d'ensembles de mesures radar Doppler de taille N, lorsque ces données font l'objet d'une identification auto-régressive (AR) par l'application notamment de l'algorithme de BURG (filtrage en treillis). La méthode s'applique également pour obtenir d'un algorithme de type Multiple Signal Classification (clusering). L'objet de la présente invention concerne un procédé permettant de déterminer un paramètre central, tel que la médiane ou la moyenne, sur un ensemble de résultats d'identification associé à un signal complexe reçu sur un ensemble de capteurs d'un récepteur, conduisant à p collections E1,...EP de M-uplets de nombres complexes qui correspondent à m mesures du signal complexe, lesdites mesures étant issues du récepteur Ek = LXkl,Xk 2,...Xk,M J k=1,..,P ledit procédé comprenant l'application d'un algorithme d'identification auto-régressive en treillis d'ordre N-1, sur ledit signal complexe, fournissant P points {W1, W2,....Wp} ou points d'ancrage qui sont invariants lors de la mise en oeuvre des étapes du procédé selon l'invention, caractérisé en ce qu'il comporte au moins les étapes suivantes : Chaque point Wk est une constitué d'une suite de N coefficients de la forme Wk = (Wk , Wk, , Wk, , Wk ') dans laquelle Wk =1og(1 E II X k•j II2) et (Wk ,..., Wk Wk-') sont les N-1 premiers coefficients de réflexion complexes provenant de ladite identification auto-régressive complexe par mise en oeuvre d'un l'algorithme d'identification en treillis appliquée à l'ensemble des nombres, (Xk,,,...,Xk.j,...,Xk.M) ledit procédé s'applique sur un ensemble général appelé H dans lequel se situent lesdits points d'ancrage et ledit ensemble général H est l'ensemble 5 2955950 des suites de la forme Wk û (Wk , Wk,..., Wk,.., WkN-') de N coefficients dont le premier coefficient est un nombre réel et les coefficients suivants des nombres complexes de module strictement inférieur à 1, ledit procédé comprend également l'application d'un algorithme de descente 5 du gradient qui opère sur une fonction définie sur l'ensemble général H, ladite fonction étant obtenue, soit comme la somme des distances des points d'ancrage à un point quelconque de l'ensemble général lorsque le procédé vise à calculer une médiane, soit comme la somme des carrés des distances des points d'ancrage à un point quelconque de l'ensemble général lorsque le 10 procédé vise à calculer une moyenne, et en ce que: ledit procédé est caractérisé en ce qu'il réalise une boucle visant à minimiser la fonction ci-dessus nommée par un algorithme de gradient comprenant un nombre fini de pas au cours desquels est calculé, à chaque 15 pas dudit algorithme, un vecteur gradient se dirigeant vers un point de convergence qui est le résultat visé par le procédé, et en ce que chaque pas de l'algorithme de gradient comporte au moins les étapes suivantes qui transforment un point z de l'ensemble général H provenant du pas précédent de l'algorithme de gradient en un nouveau point z de 20 l'ensemble général : Phase 0 : prendre le point z de l'ensemble général issu de l'algorithme à l'étape précédente j-1, pour la première étape on initialisera le point z de départ en un point d'ancrage quelconque, Phase 1, (13) : effectuer sur chacun des points d'ancrage une rotation 25 généralisée de Mebius directe de centre z pour obtenir un jeu de points P modifiés correspondant à l'étape j, W'm, Phase 2, (14) : pour chacun des points modifiés W'm, déterminer le vecteur Vm tangent à H en 0, Phase 3, (15a, 15b) : déterminer pour chacun des points modifiés, le vecteur 30 unitaire Um tangent à H en 0, 6 2955950 Phase 4, (16) : calculer la combinaison linéaire convexe des vecteurs unitaires pondérés par un coefficient am, Phase 5, (17) : multiplier la combinaison linéaire S par un nombre réel correspondant à l'incrément du gradient à l'étape j afin d'obtenir un vecteur 5 tangent à H en 0, Phase 6, (16) : enrouler le vecteur A sur la variété riemannienne H le long de la géodésique issue de 10 Z'= 0 pour obtenir le point Z' de H. A Do, 13.11 .th(1011),..., IQN 1l .th(~oN-lh Ou bien : Z' = [AQ,Al 03. (IAlI),...,ON l Y'OION_l1)] avec 0o(x) = th(x) six 2 so et th(e~ six 5 eo x sa 15 Phase 7, (19) : effectuer la rotation inverse de Môbius de centre z pour obtenir le nouveau candidat z à l'étape j+1, Phase 8, (20) : comparer la valeur de la norme de S à une valeur seuil et si cette valeur est inférieure à la valeur seuil alors arrêter l'itération et définir l'élément central (point médian ou moyenne) comme étant le dernier point 20 déterminé, sinon, faire j=j+1 et retourner à la phase 0, ledit procédé comportant une étape au cours de laquelle connaissant le point médian ou le point moyen selon l'objet du procédé, on calcule la distance entre un point de l'ensemble H et le point médian ou la moyenne issu des étapes précédentes puis on compare cette distance à une valeur 25 seuil définie par le besoin de I 'application, et si cette distance excède le seuil on décrète présence d'une cible. La valeur de seuil peut être définie par exemple de manière expérimentale. L'invention concerne aussi un dispositif d'émission et de réception de signaux permettant de déterminer un paramètre central, tel que la médiane ou la moyenne, sur un ensemble de résultats d'identification associé à un signal complexe reçu sur un ensemble de capteurs d'un module de réception caractérisé en ce qu'il comporte au moins un processeur adapté à exécuter les étapes du procédé exposées ci-dessus. L'utilisation de ce procédé en radar doppler est de plusieurs natures : • réaliser un TFAC pour séparer les cibles utiles du fouillis assimilé à la quantité médiane évaluée sur une zone localement homogène dans de fouillis globalement inhomogène ; • séparer et comparer entre elles des zones de fouillis localement homogènes, en estimant sur chacune une quantité médiane et en calculant la distance qui sépare ces dernières. Ce problème est en particulièrement celui du radar doppler météo qui doit séparer les diverses zones homogènes d'une perturbation 15 atmosphérique globalement inhomogène. Le procédé selon l'invention fournit un algorithme rapide pour déterminer un point médian ou au besoin une moyenne sur un ensemble de résultats d'identification auto-régressif. Une mesure de la distance entre la 20 médiane et chacun des échantillons ou entre plusieurs médianes est alors disponible. Le but étant, par exemple: de lisser le signal par une valeur moyenne ou plus pertinemment une valeur médiane, 25 d'écarter des éléments anormaux, d'extraire des éléments remarquables, de comparer et classer des sous-ensembles de caractéristiques différentes. 30 D'autres caractéristiques et avantages de la présente invention apparaîtront mieux à la lecture de la description qui suit d'exemples de différentes mises en oeuvre de l'invention données à titre illustratif et nullement limitatif annexé des figures qui représentent : • La figure 1A, un synoptique d'un exemple de système pouvant mettre en oeuvre le procédé selon l'invention dans le cas d'une application radar, la figure 1B la représentation temporelle des signaux émis et reçus par une antenne du système, la figure 1C un schéma des étapes selon l'art antérieur et la figure 1 D, un schéma des étapes selon l'invention, • La figure 2A, un synoptique des étapes du procédé mises en oeuvre pour une application radar, et la figure 2B un synoptique des étapes du procédé selon l'invention dans le cas à une dimension, • Les figures 3 et 4, le résultat du procédé selon l'invention dans le cas à une dimension, • La figure 5, le résultat de l'application du procédé selon l'invention dans le cas de gaussiennes, et • La figure 6, un synoptique du procédé selon l'invention appliqué pour des matrices de corrélation. Afin de mieux faire comprendre l'objet de la présente invention, la description qui suit va comporter des rappels sur les notions de métrique riemannienne et d'autres éléments utiles à la compréhension du problème. L'explication donnée dans le cadre d'une application radar dans un premier temps à titre illustratif et nullement limitatif sera généralisée à d'autres applications. En résumé, la présente invention propose une méthode simple de calcul de la médiane par un algorithme du gradient à partir d'ensembles de mesures radar doppler de taille N, lorsque ces données font l'objet d'une identification auto-régressive (AR) par l'application de l'algorithme de BURG (filtrage en treillis), par exemple ou de tout autre algorithme ayant une fonction similaire. Cette médiane est définie dans le cadre de la théorie géométrique de l'information qui adapte le langage de la géométrie riemannienne à la collection des coefficients de réflexion obtenus par l'algorithme de BURG. La présente invention propose aussi de calculer une moyenne, définie dans le même cadre, les étapes de ces deux variantes de mise en oeuvre vont être décrites ci-après. La présente invention est de portée générale en traitement du signal à valeur complexe, c'est-à-dire pour des signaux complexes mesurés par des capteurs dont certains exemples sont exposés ci-après. Nature et modélisation du siqnal dans le cas d'une application radar. Le siqnal utile En sortie du récepteur, un radar doppler analysant un signal provenant d'un ensemble de K cibles évoluant aux vitesses Vk, prend la forme d'une somme de sinusoïdes complexes bruitées et échantillonnées à la période T: K Sn = Ak.e i,2.e.fdk.n.T + Bru ta k=1 fdk = 22k est la fréquence Doppler n = 1,...,N(À longueur d'onde de la porteuse) Le nombre de N mesures disponibles (nombre de coups au but) sur un radar moderne peut être de l'ordre de 8 à 16. Le siqnal maioritaire Le signal majoritaire dans la plupart des cases distance se présente sous forme d'un bruit corrélé (fouillis). Le signal utile est rare, le signal parasite est majoritaire, il faut pouvoir les séparer. L'analyse L'analyse vise à identifier les fréquences propres du signal qui sont les fréquences Doppler fdk pour déterminer les vitesses Vk. L'identification auto-régressive (AR) est d'une part adaptée aux signaux périodiques, puisqu'en l'absence de bruit, une somme finie de signaux périodiques, vérifie une équation récurrente de type AR, et est d'autre part adaptée aux signaux de 10 2955950 bruit corrélé dont on souhaite également faire l'analyse spectrale. L'estimation du modèle AR adapté se fait dans le cadre de la présente invention par une analyse en treillis selon l'algorithme de BURG amélioré par une régularisation. 5 Cette analyse qui s'interprète comme un procédé d'orthogonalisation de type Hilbert-Schmidt, substitue au calcul des coefficients du modèle AR, le calcul de N-1 coefficients complexes appelés coefficients de réflexion, nécessairement de module strictement inférieur à 1, qui permettent par un calcul connu de l'homme de l'art de retrouver le spectre et les fréquences 10 propres du signal. Modélisation du signal radar Le signal majoritaire prend la forme d'un bruit corrélé. II est donc naturel de considérer que : • chacun des N-uplets de mesures est la réalisation d'un vecteur 15 aléatoire gaussien complexe de dimension N, dont la distribution est caractérisée par un vecteur moyen et par une matrice de corrélation (complexe et hermitienne). • le signal est stationnaire sur cet horizon de quelques mesures. Ce qui revient à dire que la matrice de corrélation est de forme Toeplitz 20 (constante sur chacune des sous-diagonales). On supposera dans la suite que le vecteur moyen, s'il n'est pas nul, est systématiquement retranché du signal. Densité du vecteur aléatoire complexe gaussien ZN de matrice de corrélation RN: 25 P(Z\T/ RA) --(ZN mnRAT-1 .(Ziv_mNy . RN .e ZN = ZOO Z1,..., Zk;..., ZN -1Y Zi E C T MN = MO, M1,..., Mk;..., MN -1 m E C 11 2955950 On suppose toujours mN = 0 La matrice de corrélation RN est complexe hermitienne de forme Toeplitz : CO Cl ... CN-1 5 RN = Cl CO CI CN-1 CN-2 CO Médiane et Moyenne de P variables aléatoires Gaussiennes complexes. Dans ce paragraphe, il s'agit de définir la notion de moyenne ou de médiane d'une collection de P variables aléatoires gaussiennes complexes 10 centrées N-dimensionnelles et d'établir un estimateur simple et efficace d'un point de vue algorithmique permettant le calcul effectif de cette médiane et de cette moyenne sur une réalisation. On est donc confronté à un espace de variables aléatoires complexes gaussiennes caractérisées par un jeu de paramètres complexes qui fait de 15 cet espace une variété abstraite de dimension N (variété = surface courbe au sens général). Il est naturel de se placer dans le cadre de la théorie géométrique de l'information, et pour cela, de définir une structure riemannienne pertinente sur cette variété (produit scalaire local défini sur chacun des espaces 20 vectoriels tangents dont la matrice dépend de façon différentiable des paramètres). Ce cadre permet de parler de distance entre 2 points qui sont ici chacun, rappelons-le, une distribution aléatoire gaussienne. On peut alors parler de point barycentre ou de point médian d'une collection de P points de 25 la variété. Barycentre Un point barycentre Mo sera un minimum local de la fonction définie comme la somme des carrés des distances du point Mo aux P points de la collection de départ appelés points d'ancrage. Ce minimum sera unique dans des conditions raisonnables de dispersion des P points sur la variété. A contrario, sur un cercle il n'existe pas de barycentre de N points répartis en étoile. Cependant ce type de problème ne se pose pas sur la variété considérée dans la présente demande de brevet. Domaine médian Un point médian Me sera un minimum local de la fonction définie comme la somme des distances du point Me aux P points d'ancrage. Les points médians sont rarement uniques, ils forment un domaine fermé à l'intérieur de la variété. A titre d'exemple, le domaine médian de 2 nombres réels est le segment qui joint les 2 nombres. Le gradient de la fonctionnelle à minimiser Soit 8(a, x) la distance géodésique entre un point d'ancrage a et le point x. Sur une variété riemannienne, la fonction b est dérivable en dehors du point « a ». Son gradient en x est un vecteur û(x) unitaire au sens de la métrique locale au point x, tangent à la géodésique joignant les 2 points et dirigé vers de a vers x. On peut alors atteindre la moyenne par un algorithme de descente du gradient visant à annuler le vecteur gradient de la fonction globale: P P N d(x) = E 8(xi, x)2 Vd(x) = E 8(xi, x).üi(x) l gk1(x).Uik.ull =1 i=1 1=1 k,1=1 De même un point du domaine médian sera obtenu par un algorithme de gradient visant à annuler le vecteur gradient de la fonction : P P N d(x) = E 8(xi, x) Vd(x) = E iii(x) I gk1(x).uik.u11=1 1=1 i=1 k,1=1 Ceci suppose de calculer les géodésiques ainsi que les vecteurs tangents. Préambule théorique Définition de la variété et paramétrage La présente invention repose sur l'analyse auto-régressive (AR) de chacun des vecteurs aléatoires, analyse qui se fait sous la forme de l'analyse 13 2955950 « en treillis » dont l'algorithme de BURG fournit un estimateur. Une variante connue de l'art antérieur permet d'introduire une « régularisation » dans l'algorithme de BURG réduisant la sensibilité au bruit des coefficients du modèle de rang le plus élevé du modèle. Cette analyse en treillis n'est rien 5 de plus qu'un procédé d'orthogonalisation de type Hilbert-Schmidt exploitant la stationnarité des N mesures aléatoires de l'échantillon. Ce sont les formules en treillis bien connues de l'Homme du métier et rappelées en annexe qui aboutissent à déterminer un jeu de N-1 paramètres complexes de module strictement inférieur à 1 appelés coefficients de 10 réflexion ,crl,..., UN _1 et à un Nème paramètre Po qui est la variance commune de chacune des variables aléatoires. = [Po, /11,..., f IN -1] = [el, 92,..., ON] PoER+ù{0} ED={zllzl La rotation de Môbius dans D. La recherche des géodésiques est considérablement simplifiée lorsque l'on remarque que dans D muni de la structure de variété riemannienne 2 définie par l'élément de longueur ds2 = IdzI , la rotation de Môbius de (1-IzI)centre « a » contenu dans D, qui au point z de D associe le point R[a)(z) avec R{a](z) = z û a , est une bijection de D dans D et que cette 1 transformation est également une isométrie pour la structure riemannienne. C'est à dire qu'elle laisse inchangées les distances. Elle transforme les géodésiques en géodésiques, et envoie le centre « a » de la rotation vers le point 0 origine du plan complexe. L'avantage qui résulte de l'utilisation de cette transformation est que : • la géodésique reliant 2 points quelconques a et z de D est transformée en une géodésique reliant le centre du plan complexe au point R[a](z) de D. • la géodésique ainsi transformée est le segment de droite dans le plan complexe reliant 0 à la transformée de z, c'est à dire qu'elle est portée par un rayon du disque. La rotation généralisée dans H. II est alors possible de définir la rotation généralisée T de centre A qui envoie H dans H: TA(0) = [q ù qo, R[al](zl), R[a2](z2)..., R[aN - 1](zN - 1)] 0 = [q, zl, z2,..., zN - A=[go,al,a2,...,aN I] a,ED i=1,...,N-1 avec R[a](z) = zùa qui transforme composante à composante les éléments de H et qui envoie le centre A de la rotation sur 0 = élément de H. La 1 ere composante qui est un nombre réel est translatée. Les composantes suivantes qui sont dans des nombres complexes de D subissent une rotation de Mébius élémentaire. Un vecteur tangent à la variété différentiable H au point 0 est envoyé composante à composante sur un vecteur tangent à H en TA(0)de manière isométrique sur chaque composante. Le ds2 partiel est conservé, comme le ds2 global est une somme pondérée par des coefficients constants des ds2 partiels, le ds2 global est lui-même conservé. T est donc une transformation isométrique de H. Les géodésiques issues de A sont transformées en géodésiques issues de 0 dans H. Cependant, contrairement à ce qui se passe dans D, les géodésiques dans H ne sont pas de simples droites. Seules les projections sur chacune des composantes complexes sont des segments de droites dans D. La métrique de H en 0. Précisons la notion de vecteur unitaire en 0. Dans l'espace vectoriel tangent à H en 0, la métrique se réduit à : ds2 = N.dq2 + EkNH'(N û k). dzkI2 Le plan tangent en 0 à H est alors un espace vectoriel euclidien muni du produit scalaire qui se déduit de l'élément de longueur ci-dessus. En effet la métrique est euclidienne aux pondérations près et l'espace vectoriel est un espace vectoriel réel de dimension 1 + 2*(N-1). Equation de la géodésique dans H entre 0 et un point W quelconque de H. 2 ds2 =ao.dro2k ~Idzkl2 ao=N ak=Nùk -IZkI k=1,..N-1 17 2955950 Passons en coordonnées polaires (r module de z et e sa phase) sur un chemin joignant le point 0 au point W = (Wo,W1,..Wk,..,WN -1)T dans H: ds 2 = ao.dro 2 + ~N-lak drk' +n2 .d 612 k-1 2 (1ûrk2) 2 =dxk2 +dyke =drk2 +rk2.dOk2 avecxk =cos(Bk) yk =sin(Ok) dzk dzk = dxk + i.dyk 5 L'intégrale qui doit être minimisée pour obtenir la géodésique s'écrit symboliquement : 2ù W 2 N-1 drk' +rk2.d8k2 I~~ 2 N-1 drk j = ~ao.dr0 + k-1 ak. 2)2 ao.dro + k-1 ak. 2 0 0 (1 ù rk o (1 ù rk ) La fonctionnelle est minorée par la valeur qu'elle prend lorsque tous les 0, sont constants, ce qui signifie que sur chaque disque de Poincaré D, la 15 projection de la géodésique est une portion de rayon du disque. Le problème se déplace donc dans H'= R x [0,1[N-' ds 2 = ao.dro Z + r l'ak.d (arg th(rk)) En faisant le changement de variable 4' H'-> RN [ro, ri, r2,..., rN -1] I >[ro, arg th(ri), arg th(r2),..., arg th(rN -1)] = [Ro, Ri, R2,..., RN -1] Dans cet espace la métrique définie par l'élément de longueur ds devient : ds2 = ao.dRo2 +~k l~ak.dRk2 10 20 25 Puisque les ak sont constants, il s'agit de géométrie euclidienne, les géodésiques sont des segments de droites de RN d'équation paramétrique en À 2 E [0,11: i(2) = [2.Ro max, .t.Rl max,..., /i..R(N - 1) max] Longueur de l'arc géodésique, La longueur d(0,W) de l'arc géodésique entre le point 0 de H et le point W de H, après la transformation P, peut se calculer dans l'espace euclidien RN puisque l'arc est en l'occurrence un segment de droite: Nù1 IIZ J (1) = ao.Ro .2+ Ek=1 ak.Rk max = d(0, W ) Rkmax=argth(rkmax)=argth(IWk) k= 1,..., N-1 Ro max = ro max = IWoI = q max puisque la première composante est inchangée dans toutes les transformations envisagées. La longueur de l'arc géodésique dans RN entre 0 et le point de paramètre À sera alors : 2.D(0,W) et cette longueur est inchangée au cours des différentes transformations 20 puisque ces transformations sont isométriques et conservent donc les longueurs. Equation paramétrique de la géodésique issue de 0 dans H: Revenons dans l'ensemble H par transformation inverse pour obtenir l'équation paramétrée par À de la géodésique: 25 9%x(il) = [2.Ro max, exp(i 81).th(2.Rl max),..., exp(i B(N 1).th(2.R(N - max)] 8k est la phase de Wk dans D. En effet, comme on l'a précédemment noté, la géodésique dans H, en restriction à l'une des composantes complexes, est 10 15 19 2955950 un segment de rayon du disque D. ek désigne simplement la phase dans D de ce rayon. exp(iO) = Wk IWkI indépendant de À. H(2)= 2.Wo, WI .th(2.argth(W1I),..., WN WN II .th(/t.argth(IWN-'I) IWII 5 et ITtH(1)II= (1)11=./ao.Womax2+1%T -1 ak.argth(Wk)2 =d(O,W) Vecteur tangent à H en 0. Un vecteur tangent à l'arc paramétré dans H s'obtient en dérivant H(2,I) par rapport à h: V(2) = [Ro max, exp(i ei).RI max .(1 th 2 (2,.R) max)),..., exp(i B(N - 1)).R(N -1) max .(1- th 2 (2.R(N -1) max))] En À =0, c'est à dire au point 0 de H on obtient : 15 V (0) = [Ro max, exp(l01).RI max,..., exp(i9(N - I)).R(N -1) max)] Mais : exp(ia) = Wk IWkI Rk max = argth(Wkl) k=1,...N-1. Le vecteur tangent V(0) (qui n'est pas encore unitaire) s'écrit donc : V(0)=[romax,IW W) WN -1 ll.argth(IWII),...,IWN I l •argth(WN-II)] Il reste à normaliser dans la métrique locale en 0 (voir paragraphe précédent) pour obtenir le vecteur U(0): unitaire : U(0) = v (o) IIV(0)II 20 25 5 10 20 25 avec : IIV(0)II= Jao.goX2+Ek I'ax.argth2(IWkI) =d(0,W) c'est à dire : IIV(0)II = IN.gomax2+Ek 1'(Nûk).argth2(WkJ) = d(0,W) Equation paramétrique de la géodésique en fonction d'un vecteur tangent à H en O. L'équation paramétrique de la géodésique s'exprime en fonction d'un vecteur quelconque tangent à H en 0 : V = (Vo,V1,..,.VN-1) Lorsque le maximum de [2.IVk ] est petit devant l'unité pour k = 0 à N-1, l'approximation suivante est valide: 9ÎH(2) [2.Vo, 2.Vi /Î.VN -1] Expression de la rotation généralisée inverse. TA(e) = [q û qo, R[al](zi), R[a2](z2)..., R[aN - 1](ZN - ,)] TA'(9) =[q+qo,R-' [al](zl), R-' [a2](z2)..., R-1[aN- ](zN-1)] e=[q,ZN-1] A = [go, a,, a2,..., aN - ,] R-1 [a](z) = z + a 1+â.z avec la condition que a et z sont dans D. 2.Vo, V1 .th(2.IVii),..., VN - 1 th(/1.IVN - 1 ) IV1I VN - Il - RH(2) = 21 2955950 La distance entre les points U et V de H s'écrit selon une formule connue de l'art antérieur : Uo = log(Puo) Description de l'invention Schéma synoptique La figure 1A schématise un exemple de système permettant la mise en oeuvre du procédé selon l'invention pour l'exemple de l'application radar. Sur cette figure sont représentés un dispositif 1 d'émission et de réception de signaux, comprenant un module émetteur 2 d'une onde émise Se allant se réfléchir sur une cible 3, un module de réception 4 du signal Sr réfléchi par la cible et un processeur 5 relié au module de réception et adapté à mettre en oeuvre les étapes du procédé. Le dispositif d'émission-réception 1 comporte une sortie 6 reliée par exemple à un dispositif de traitement 7 des informations issues du procédé. La figure 1B représente une antenne et le lobe d'émission. L'antenne effectue un balayage angulaire. La cible se trouve dans le lobe de l'antenne et réfléchi un signal. Les cases distances habituellement utilisées dans le domaine des radars sont représentées. Un diagramme temporel représente la répartition et l'enchainement des émissions et des échos cibles qui sont reçus au niveau du récepteur. Le temps aller-retour est égal à At=[2(d+vt)c]. T correspond à une période de répétition des impulsions. Pour une case distance k, où il y a M coups au but. Après réception sur les voies I et Q, soient X1, ..Xp les p mesures complexes du signal reçu. Avant démodulation X A e -12g (f°+MlT = + Bruit 1 l L(U,V)=1/N.(VoùUo)2+Ik 11(Nùk).argth2( Vk-Uk 1ù Vk.Uk ) Vo = log(Pro) Après démodulation les mesures s'expriment sous la forme : XI = Ale-'2' (f» T + Bruit Où fd est la fréquence Doppler pour la case k et fo la porteuse La figure 1C schématise le traitement des signaux après démodulation selon l'art antérieur. On remarque la mise en oeuvre d'un ensemble de filtres FFT sur chacune des cases distance pour déterminer une valeur de TFAC et une valeur qui sera comparée à une valeur seuil pour déterminer la présence ou non d'une cible. La figure 1 D représente l'enchainement des étapes du procédé selon l'invention qui détermine pour chaque échantillon reçu par une analyse AR et pour une case distance donnée une valeur des N coefficients de réflexion, 0 k pour j= 1, N. Les coefficients de réflexion sont ensuite transmis au module de calcul de la médiane TFAC pour déterminer la valeur des coefficients de réflexion médian 0jm. Les coefficients de réflexion médian sont ensuite utilisés pour effectuer le calcul de la valeur de la distance pour une case distance k donnée. Si cette valeur est supérieure à une valeur seuil, alors le procédé décrète la présence d'une cible dans la case distance, sinon, il ne détecte pas de cible. Les étapes du schéma de la figure ID relatives au procédé selon l'invention vont être détaillées ci-après. La figure 2A représente un synoptique des étapes du procédé mises en oeuvre pour une application radar, et la figure 2B un synoptique des étapes du procédé selon l'invention dans le cas à une dimension. L'alqorithme de gradient pour un point médian ou la moyenne û FIG.2A Initialisation : La donnée de départ est un ensemble de P échantillons {E1,... EP} de M mesures complexesM N, 10 23 2955950 L'algorithme de BURG (éventuellement régularisé), 11, appliqué à chacun de ces échantillons génère P points de l'ensemble H qui sont des N-uplets de coefficients, l'ensemble étant référencé Bp, 11, sur la figure 2A. La 1 ere composante du N-uplet Bp est le logarithme népérien de la 5 puissance moyenne estimée de l'échantillon qui sera traité différemment des autres coefficients. Les composantes suivantes sont les N-1 coefficients de réflexion complexes du modèle AR. II en résulte P points de H : {W~,...,WP} qui sont les points d'ancrage de 10 l'algorithme qui restent fixes tout au long de la procédure. L'initialisation (étape j = 0), démarre avec un point candidat Z qui peut être n'importe quel point, par exemple l'un des P points d'ancrage ou l'origine 0 de H. La boucle : 15 La boucle sur l'étape L: qui se décompose en 7 phases. A l'étape j, de l'algorithme, on dispose d'un point Z E H candidat issu de l'algorithme à l'étape précédente j-1 Z = (zo,zl,...zN -1) R1-1 phase 1, 13 Effectuer sur chacun des P points d'ancrage {W1,...,WP} la rotation 20 généralisée directe de Môbius de centre Z pour obtenir un jeu de P points modifiés dans H {W'1,...,W'P}propre à l'étape j. Wm = (Wm,O, Wm,1,...Wm, N -1) m=1,..., P est transformé dans cette phase en : 25 W, m = (-W m,O, W~ m,1,...W' m, N -1) = (W' m,O - Z0, R[ZI](Wm,l ),..., R[zn - 1](Wm, N - 1)) Où la première composante w'm,o qui correspond à la puissance du signal est translatée de z0, tandis que chacune des N-1 autres composantes subit une rotation de Môbius élémentaire 24 2955950 R[z](w w - z 1ù z.w R1-2 phase 2, 14 Pour chacun des points modifiés W'm, calculer le vecteur Vm (qui doit être interprété comme un vecteur tangent en 0 dans H) selon les formules. W'm,l W m,N-1 Vm = [w '.,o, arg th(W m,l ),..., arg th(W m, N - Il)] 1W'm,11 W'm,N-11 Pour échapper aux divisions par 0 une variante consiste à effectuer le calcul suivant : Vm = [ W m,o, w' m,1 .0V (w'71 1),..., W m,N-1 •0V (W ,, N_1 )] avec : w (x) =argth(ev) six e, et argth(x) six ev 6'v x ev est un nombre réel positif petit devant 1 par exemple 0.1 R1-3-a phase 3, 15a si l'objet de l'algorithme est le calcul d'un point médiane Pour chacun des points Vm, calculer le vecteur unitaire Um (qui doit être interprété comme un vecteur tangent en 0 dans H) selon les formules. Um = avec I vmll ~Vm = . z + ,k '(N ù k). arg th 2 (Vm,k Pour échapper aux divisions par 0 une variante consiste à exécuter le calcul 25 suivant : 5 10 15 20 Um = lVmll Si llVmll ~ SU et Sm si lVm _ 0 am nombre réel On prendra par exemple : am = 1 P 15 R1-5 phase 5, 17 Multiplier le vecteur S par le nombre réel y., > 0 incrément du gradient à l'étape j, pour obtenir un vecteur A qui est encore un vecteur tangent à H en O. A=y,..S 20 En développant selon les composantes de H, 0 s'écrit: A = [Ao, A1,..., Ak,..., ON - 1 ] R1-6 phase 6, 18 Calculer le vecteur Z' de H qui s'interprète comme I«( enroulement » du vecteur A sur la variété riemannienne H le long de la géodésique issue du 25 point 0. 26 2955950 r ° Z'= °o, ' th(I°,I) ,...° N .th(I°Nù1I) I I I - °1 °Na Pour échapper aux divisions par 0 la variante de calcul suivant sera utilisée: Z' = [° °1 çb° (I°, ),..., ° N_1 çb° (I°N 1 )j avec çb° (x) = thzx) six >_ e° et th(e) six e° e° est un nombre réel positif petit devant 1 par exemple 0.1 On pourra prendre également l'approximation. Z' ,-=, [°0,°1 °1] si toutes les composantes sont de module petit devant 1. Posons : Z'= [z'o,Z'l,...,Z'k,...,Z'Nù1] R1-7 phase 7, 19 Effectuer sur Z' la rotation inverse de Môbius généralisée de centre Z qui est le point candidat au pas j pour obtenir le nouveau point candidat Z au pas j+1. Znouveau = Tzù1 (Z') = [z'+zo, R 1 [zl ](z'1 ), Rù' [z2 ](z'2 )..., R 1 [ZN 1 ](z'N_1 )] i = 1,...N-1 R1-8 critère d'arrêt, 20 La boucle s'arrête lorsque la norme du vecteur S obtenu à la phase 4 descend en dessous d'un nombre réel positif proche de 0 qui est le seuil d'arrêt de la boucle. IISII = .\jN.So + Lk l'(N ù k). arg th 2 ( Sk ) seuil arrêt R-' [z1 ](z' z' +Z 1+z;.z'1 avec : 27 2955950 Le point médian ou la moyenne selon l'obiet de la procédure est le dernier point Z calculé, 21. R1-9-a, 22 rétablissement de l'échelle linéaire pour la puissance médiane. 5 Après arrêt de la boucle, comme la première composante du N-uplet considérée correspondait au logarithme népérien de la puissance moyenne de l'échantillon, la puissance médiane ou moyenne finale, (suivant que l'objet de la procédure est le calcul de la médiane et celui de la moyenne), s'obtient en rétablissant l'échelle linéaire de la première composante: 10 Puissance Médiane = exp(Zo) R1-9-b, 23 rétablissement de l'échelle linéaire pour la puissance moyenne. 15 Puissance Moyenne = exp(Zo) Une étape après avoir déterminé le point médian ou la moyenne peut être la suivante : connaissant le point médian ou le point moyen selon l'objet du procédé, on calcule la distance entre un point d'ancrage et le point issu du 20 procédé par l'application, par exemple, de la formule précédente : puis on compare cette distance à une valeur seuil définie par le besoin de I 'application, et si cette distance excède le seuil on décrète présence 25 d'une cible. L'exemple qui suit correspond au cas particulier où l'algorithme utilise un seul coefficient. L(U,V) = N.(Vo-Uo)2 +~k ~'(Nùk).argth2( VkùUk 1 - Vk.Uk ) 28 2955950 Lorsque le traitement radar envisagé n'utilise que le premier coefficient de réflexion qui est alors le coefficient de corrélation, l'algorithme se simplifie. Dans ce cas en particulier la puissance du signal est considérée comme égale à 1. Le logarithme de la puissance est donc égal à 0 et la 5 première composante de H est toujours nulle. L'espace H se réduit alors au disque D. C'est par exemple le cas dans un mode perturbation d'un radar doppler météo. 10 Algorithme du gradient pour un point médian ou moyen en dimension 1 ù FIG.2B Le cadre de travail de l'algorithme est alors le disque D de Poincaré D = {zl z E C et Izl W'm=R[z](Wm)- Wmùz 1ùz.Wm R2-2-a phase 2, 31a si l'objet de l'algorithme est le calcul d'un point médian (idem R1-2 et R1-3a) Calculer pour chacun des points modifiés le nombre complexe unitaire Um selon les formules. W' m m IW'1 Pour échapper aux divisions par 0 on préférera le calcul suivant : Ov(x)= 1 six ≤. ev et 1 si xev ev x avec Um = W'm •Y7 ('W'm I ) ev est un nombre réel positif petit devant 1 par exemple 0.1 R2-2-b phase 2, 31b si l'objet de l'algorithme est le calcul de la 15 moyenne (idem R1-2 et R1-3b) Ne rien faire, et renommer le nombre W'm en Um. R2-3 phase 3 (idem R1-4), 32 Calculer la combinaison linéaire convexe S des nombres Um pour m=1 à P. 20 P S=Eam.Um m=lp l an, =1 m=1 am >_ 0 25 réel et coefficient de « pondération » des nombres unitaires Um pour effectuer la combinaison linéaire convexe 30 2955950 On pourra prendre par exemple : R2-4 phase 4 (idem R1-5), 33 Multiplier S par le nombre réel y,. > 0 où y, est l'incrément du gradient à l'étape j, afin d'obtenir le vecteur A S.. A=y;.S R2-5 phase 5 (idem R1-6), 34 « Enrouler » le vecteur A sur la variété riemannienne D le long de la 1 o géodésique issue de 0 pour obtenir le point z' de D. z'= IQ .th(AI) Pour échapper aux divisions par 0 on préférera le calcul suivant : z'= A.Oo (DI) avec 0o (x) = th(x) si x >_ eo et . th(e°) six R2-6 phase 6 (idem R1-7), 35 Effectuer sur z' la rotation inverse de centre z pour obtenir le nouveau nombre candidat z à l'étape j+1. 5 25 Znouveau = Tz 1 (Z') = R ' [z]( Z' Z'+Z 1+ z.z' R2-7 critère d'arrêt (idem R1-8), 36 La boucle s'arrête lorsque le module du nombre complexe S de la phase 4 est inférieur à un nombre réel positif proche de 0 qui est le seuil d'arrêt de la 5 boucle : IS 5 seuil arrêt Une étape après avoir déterminé le point médian ou la moyenne peut être la suivante : connaissant le point médian ou le point moyen selon l'objet du procédé, on calcule la distance entre un point d'ancrage et le point issu du 10 procédé par l'application de la formule précédente : puis on compare cette distance à une valeur seuil définie par le besoin de I 'application, et si cette distance excède le seuil on décrète présence 15 d'une cible. Les fiqures 3 et 4 représentent le résultat de l'algorithme de médiane en dimension 1 Par exemple, sur la figure 4, 10 points complexes dans D schématisés par des (X) points notés Wi et correspondant aux points d'ancrage ou 20 coefficients de corrélation dans cet exemple de réalisation et plusieurs points + notés xi représentent les points candidats qui vont intervenir dans l'algorithme. La convergence de l'algorithme est notamment représentée par les différents points xi et la valeur médiane par le carré sur la figure. En prenant les étapes décrites précédemment, après application de la 25 rotation de Môbius de centre xl sur l'ensemble de la figure comprenant les différents points Wi et xl , le point xl se retrouve au centre du cercle unité C, donc en xO (voir la flèche FI sur le schéma). Les points Wi sont devenus des points modifiés Wi' (qui ne sont pas représentés sur la figure). Les arcs géodésiques à déterminer entre le point central et les différents points L(U,V)= N.(VoùUo)2+k~'(Nùk).argth'( Vk-Uk 1 - Vk.Uk ) 32 2955950 modifiés Wi' correspondent à segments de rayon du cercle unité. On calcule pour chacun des points modifiés Wi', le nombre complexe unitaire Um correspondant qui peut être interprété comme un vecteur tangent en 0 à la variété riemannienne D selon les formules R2-2-a dans le cas de la médiane, puis on effectue une combinaison linéaire convexe S des nombres Um pour m = 1 à P le nombre de points d'ancrage considérés. Lors de l'étape suivante, on va avancer dans la direction du vecteur S sur une longueur d'un pas y pour obtenir un nouveau point x2 de la variété différentiable. x2 est un nouveau nombre complexe. 1 o Puis on va effectuer une rotation de Môbius inverse de centre xl sur le point x2 pour obtenir le nouveau candidat à l'étape j+1. Les points modifiés Wi' n'ont pas à être inversés puisque cette opération aurait pour seul effet de reconstituer les points d'ancrage Wi qui sont de toute façon conservés en mémoire. 15 Ces étapes vont être réitérées tant que le critère d'arrêt n'est pas vérifié, pour obtenir la médiane représentée par le carré sur la figure. Alqorithme du gradient un point médian en dimension réduite Dans l'exposé explicatif initial, il a été supposé que la dimension N de l'espace H était égale à celle des échantillons de base sur lesquels 20 fonctionne l'algorithme de BURG régularisé. Cette hypothèse n'est nullement obligatoire. Disposant d'échantillons de taille M, l'algorithme de BURG régularisé peut être interrompu à un nombre d'étages N inférieur à M. L'espace H sur lequel l'invention fonctionne les algorithmes de médiane 25 et de moyenne prend alors la dimension N. Ceci est particulièrement utile dans les applications où les coefficients de réflexion d'ordre élevé sont peu significatifs, leur estimation consommant inutilement du temps calcul. Toutes les étapes qui ont été décrites pour le calcul en une dimension 30 en plusieurs dimensions s'appliquent. Application à la mise en groupe ou « clusterinq » de données Doppler pour la cartoqraphie de fouillis. Les techniques de regroupement ou « clustering » (cherchent à classifier des données qui sont déjà regroupées dans des classes décrites par des variables gaussiennes réelles non centrées et non réduites N(m,Q). Le regroupement de telles classes suppose de disposer d'une méthode permettant de reconnaître une entité centrale à tout un groupe de variables gaussiennes. La théorie géométrique de l'information est donc couramment ~o invoquée pour résoudre ce type de problème. II est évident qu'un algorithme déterminant la médiane d'un ensemble de gaussiennes est le bienvenu et de nombreux travaux vont dans ce sens. La présente invention s'applique immédiatement au problème de la détermination de la médiane d'un ensemble de P variables gaussiennes 15 réelles à 1 dimension. En effet dans le cadre de la théorie géométrique de l'information une variable aléatoire gaussienne réelle est vue comme l'élément d'une variété réelle à 2 dimensions paramétrée par le couple (m,a) avec m appartient à R ensemble des nombres réels et a appartient à R+* ensemble des nombres 20 réels strictement positifs. Cet ensemble ouvert est identifié au demi-plan complexe supérieur que nous noterons H. Cette variété est munie d'une structure riemannienne au moyen du tenseur métrique défini au point (m,a) par la matrice d'information Fisher de la densité gaussienne calculée en ce point. 25 H+ muni de sa structure riemannienne est le'/2 plan de Poincaré. La transformation de Cayley u = z ù l z = ùi. u +1 met H+ en bijection z+i u-1 bi-continue et bi-différentiable (difféomorphisme) avec le disque de Poincaré D. Il s'agit d'un changement de paramétrage (carte locale). Le tenseur métrique de H+ se transfert sur D pour obtenir le tenseur de Poincaré dans D. 5 34 2955950 Le tenseur de la géométrie de l'information sur les variables gaussiennes est égal à 8 fois le tenseur dans le disque de Poincaré. Tenseur de la géométrie de l'information Densité gaussienne : 1 1 1 x-m \ z f(xm,6) _ ù.exp(--. ) f2 -T, o 2\ 6 log-vraisemblance : 1 1 1 "xûm m,ff))=û2.log(21r)--.log(6)--. matrice d'information de Fisher : L(x,m,a)=1og(f(x 1 0 F(m, u) = ûE a 2 amati L ds2 =1E aL aLT dB;=ûEE a2 a2 - am2 L amati L a2 a62 L 2 6_ a2L de..de; ae;ae; 0 2 si l'on pose : x = y=6 z=x+i.y ( m\ z 2 dS2=I.ds2=1 dm2+2.do-21 dv.y +d6 -1 dx2+dy2 ldzlz 8 8 62 4 6-z 4 y2 4 (Im(z))2 Transformation de Cayley et passage de H+ dans D : z û a . u+1 u= qz= -a. z +i uû1 dS2 == 1 ldzlz - ldulz 4 (Im(z))2 (1û ul2 }2 35 2955950 L'algorithme de recherche d'un point médian/moyen parmi un ensemble de Gaussiennes L'algorithme de recherche d'un point médian ou moyen parmi un ensemble de P gaussiennes réelles consiste alors à appliquer la 5 transformation de Cayley T afin de transformer les P couples (mk, ok) en P points dans le disque D : (mk,ok)=Zk=-+1.0k >Uk=T(Zk)= Zk-1 k=1,...,P Zk+l 10 L'application de l'algorithme de médiane à une dimension dans D fournit un point médian ou moyen UD auquel est appliquée la transformation de Cayley inverse T-1 de manière à obtenir un couple (mD,QD) qui est le paramétrage de la gaussienne médiane ou moyenne. 15 MD . ZDT-'(UD)=û1 UD+I ZD= +1.OD(mD,Ore) UDù1 2 La figure 5 est le résultat de 10 gaussiennes (mk,ok) croix « x », la médiane étant représentée par le carré « ^ ». Algorithme de gradient moyenne/médiane pour matrices de corrélation de forme Toeplitz. Lorsque les blocs de données échantillonnées I et Q issus du récepteur sont présents sous la forme de matrices de corrélation, la méthode proposée par l'invention permet encore d'obtenir une matrice représentative de l'ambiante (TFAC). En effet des formules connues de l'homme de l'art permettent de passer des matrices de corrélation aux coefficients de réflexion et inversement. L'objet de ce paragraphe est de rappeler ces formules qui permettent l'extension de l'invention et traite donc du passage d'une matrice de corrélation complexe hermitienne définie positive de forme Toeplitz vers les coefficients de réflexion du modèle AR complexe équivalent et inversement. La matrice de corrélation RN est définie par la liste des coefficients de corrélation RN = Co C1 et Co CN_I CN-2 La liste des coefficients de réflexion est : {Poäu1,..,un,...,UN.} avec PO puissance du signal : Po = co Les coefficients de réflexions s'obtiennent partir des coefficients de 10 corrélation cn, de proche en proche, par l'algorithme itératif suivant sur l'indice n. Cette itération est déduite de l'algorithme de BURG régularisé et s'appuie sur des coefficients ak intermédiaires. L'itération s'écrit sous la forme de la boucle suivante : 15 ao =1 Pour n = 1 au pas de 1 jusqu'à N faire : { n-1 n-1 n-1 (a0 ,al , ,ak , n-1) /,in (a0 n, ,akn, ,an , an-1 ' ~ al > an ) Passage des coefficients ak-' du pas n-1 au coefficient de 20 réflexion ,un : n-1 n-1 Eak'.an (1+1).Ck-(1+1) +E Mn).ak-1 an k k,1=0 k=1 ~n ù n-1 n-1 an-1 an'.ck-1+ E Mn) jar' l2 k 1 k,1=0 k=0 Dans cette formule : ~kn) = 71.(2z)2.(k û n)2 Passage des coefficients ak-' et ,un aux coefficients ak du pas aô =1 ak = ak ' +,un an=k , k=1,...,n-1 n an = ,un } On prendra par exemple yi = co.10-4. A ce stade les coefficients de corrélations sont transformés en une suite de coefficients {1og(Po)äu1,..,un,..yuN.} qui constitue un point de l'ensemble H utilisable par l'algorithme de gradient ou de moyenne décrit aux chapitres précédents. La procédure précédente est évidemment appliquée à toutes les matrices de corrélation (de 1 à p) du jeu de données initial de matrices. Ceci fournit P points d'ancrage dans l'ensemble H auxquels est alors appliqué l'algorithme de gradient pour la médiane ou la moyenne en dimension N. Le résultat de l'algorithme de gradient est un élément de l'ensemble H qui est l'estimation du point médian ou du point moyen dans H. Le point de H ainsi obtenu est alors transformé en retour en une matrice de corrélation par l'application de l'itération suivante qui construit une suite de matrices carrées symétriques Rn de taille croissante. A chaque pas d'itération n, une matrice Rn de dimension (n+1)x(n+1) est obtenue à partir de la matrice Rn_1 de dimension (n x n) calculée à l'itération précédente en utilisant des coefficients complexes intermédiaires ak et un coefficient réel positif Enz : Au pas 0 : matrice(lxl) Ro = Po 602 = Po n: 5 38 2955950 Pour n= 1 au pas de 1 jusqu'à N faire : { I Le résultat de l'algorithme est RN, qui, suivant le type d'algorithme utilisé, est la matrice médiane ou moyenne recherchée. ag =1 avec ak ak-' +,un an-k' k=l,...,n-1 ann = /(n n-1 An = aQ an k Rnù1 ùRnù1.An ù Rn ù ].An Enz + A* .Rn ù I.An_ Rn = 39 2955950 Annexe 1 Exemple de programmation en dimension 1. Cet exemple montre la simplicité de l'algorithme en dimension 1, 5 z complexe ; w[m] m=1 ..., M points d'ancrage dont on cherche la médiane tous sont des nombres complexes de module inférieur strictement à 1. Arrêt = FAUX ; y = 1 ; z = 0; 10 j=0; tant que ((j Annexe 2 Identification autorégressive en treillis et coefficients de réflexion On dispose d'un ensemble de N variables aléatoires à valeurs complexes de moyenne nulle E = {Yo,...YN -1} Par identification auto-régressive, on entend la recherche d'un jeu de N-1 coefficients complexes a, tels que la quantité réelle: E N-1 YN-1 ûai.Y(N-1)-i) i=1 soit minimale. 2 Pour résoudre ce problème on se place dans l'espace des variables aléatoires complexes de carré intégrable sur lequel (X, Y) = E(X.Y) est un produit scalaire et l'écart type est une norme vectorielle. Cet espace est un espace de Hilbert. On peut alors adopter un langage géométrique et parler de projection orthogonale sur un sous-espace vectoriel. Notons P(X /Y) la projection orthogonale de la variable X sur l'espace linéaire Y. Pour 0_ e+k(t)=Y-P(Yt/Yt-1,...Yt-k) o(t) = eù o(t) = Yi e k(t)=Y-kûP(Ytk/Yt,...Y-k+1) Par stationnarité En+(t) et En (t) ont les mêmes statistiques pour tout t. En particulier : 41 Ile+k(t)II2 = IIE k(t»II2 = EIe+k(t)I2 = EI E k(t')I2 pour tout k, t et t' Par orthogonalité on a : E+k(t) = Yt ûP(Yt/Yt-1,...Y -k) = YtûP(Yt /Y 1,...Yt-k+1)ûP(Yt /E-k-1(tû1)) 5 _ E+k -1(t) +,U1k.E k -1(t -1) E k(t)=Y-kûP(Yt-kY,...Yt-k+1)=Yt-kûP(Y-kYt-1,...Yt-k+1)-P(Yt-kE+k-1(t)) = .E k -1(t -1) +,u2k.E+k -1(t) (e+k-1(t)£ k-1(t-1)) -1)) ,u1k(t) _ - 2 , IIE k-1(t-1)Il \ (E k-1(t-1)IE+k-1(t)) 10 ,ct2k(t) _ ù 2 t ,uk IIE+k -1(1 ,u2k=,ulk indépendants de t. Formules récurrentes en treillis : 15 e+o(t) = e-o(t) = Yt pour k = 1 jusqu'à N-1 : Jek(t) = E+k -1(t) +,uk.E k -1(t -1) E k(t) = e +k - 1(t - 1) +,Uk.E+k -1(0 Coefficients de réflexion dans leur formulation théorique: 20 E(E+k - 1(t).e k - 1(t -1)) ,uk = û EIE-k-1(tû 1)I2 ù E(E k - 1(t -1).E+k - 1(t)) ,uk = - 2 EI E+k - 1(t)I E(e+k-1(t).e k-1(t-1)) .JEI£+k-1(t) 2..\ Ele k-1(t-1) 2 => t[lk E D Les innovations arrières forment un système orthogonal (procédé d'orthogonalisation de Hilbert-Schmidt). En prenant t = N-1 : E-o(t)=Yt, e-1(t)=Yt-1-P(Yt-1/Yt), E 2(t)=Yt-2-P(Yt-1/Yt,Yt ) 6 k(t)=Yt-kûP(Yt-1/Yt,Yt-1...Y-k+1) Les coefficients de réflexion sont les coefficients de la projection orthogonale de Yt sur la base orthogonale des {£-k(t -1) /k = 0...Nû1} Y = -1 k=1 pk.& k - 1(t -1) + e+n(t) N-1 YN - 1 = [-k 11,L1k.£ k - 1((N -1) -1)] + £+N -1(N -1) = E ai.Y(N -1) - i + £+N -1(N -1) i=1 20 Les variances des innovations peuvent être exprimées en fonction des coefficients de réflexion et de la variance commune Po des Yt : 1le 2 = II£ ,(t) r = (1 - ~kl Z ).11£+k 1(t)IIZ = fl (1 - 1/4 2 ). YII2 = n (1 - /'I ».Po f-1 J=1 La matrice de passage de la base des Yt à la base orthogonale des 25 innovations arrières est semi-diagonale inférieure : 15 10 k-1 e-k(t) =Yt-kùP(Yt-1/Y,Y-1...Y-k+l)=Y-k+E%3k,J.Y-i J-0 La famille des e'- k(t) = s-k(t)2 est orthonormée. Os k(t>II rit = [Yt,..., Y - k,...Y - N + l ]T ~.t = [e". 0(0,...,60- k(t),...eN - 1(t)]T = L.IIt avec : E(IIN -1T.IIN -1) = RN matrice de corrélation du vecteur aléatoire. /30,0 0 0 0 0 0 0 0 L /k,1 /3k, k 0 0 0 /3Nù1,1 /Nù1,k ... /N-1,N-1 /3k, k = 1 IIS k(t)II E(EN.EN*) = IN = LN.E(IIN.HN").LN* = LN.RN.LN* Où E est l'espérance mathématique et RN matrice de corrélation de Ili , IN la matrice identité à N dimensions. N-1 N-1 = Ile k(N-1)IIz =PONn(1-I,uk2)N k k=1 k=1 1 2 det(RN) = 1 N-1 det(L)2 fl/3k,k k=0 15 N-1 ln(det(RN)) = N. ln(Po) + (N ù k).ln(1ù ,uk 2 ) k=1 44 2955950 Annexe 3 Entropie de la variable qaussienne et métrique Kâhlerienne dans H : 5 La densité gaussienne s'écrit : p(ZN / RN) _ IRN -I . exp(ùZN.RN ' .ZN*) = NN .1 RN -I . exp(ùtrace(RN ' .ZN.RN * )) 7z- 7r û ln(p(ZN / RN)) = N.ln(ir) + 1nJRNI + ZN.RN .ZN* = N.ln(r) + 1nIRNJ + trace(RN ' .ZN.RN* ) L'entropie (D de la variable gaussienne centrée complexe de dimension N s'écrit : (D(RN) = E(ûln(p(ZN/RN))) ù N. ln(TC) + 1nIRNI + trace(RN .E(ZN.ZN * )) = N.1n(z) + lnIRNI + N N-1 15 ch(RN) _ (D(Poäu1,...,uN -1) = N.ln(rz.e) + 1n RNI = N.1n(rz.e.Po) + E (N û k).ln(1û ,ukI2 ) k=I 19 2Po2~(PO, 11,...,uN -1) _ ù N Poe Nûk (1 _ II2 ) 2 Les dérivées partielles croisées sont nulles. Le Hessien de (D est donc 20 diagonal. D'où l'expression du produit scalaire dans l'espace vectoriel tangent à la variété H au point (Po, ui,...pN - 1) .10 Annexe 4 L'algorithme de BURG régularisé L'algorithme de BURG fournit sur un échantillon de N mesures z(k) k= 1,..,N un estimateur des coefficients de réflexion obtenu dans leur formulation théorique en annexe 2 Cet estimateur est régularisé et stabilise les coefficients d'ordre > 1. Initialisation : fo (k) = bo(k) = z(k) , k=1,...,N (N : nombre de mesures.) N P = ù.E 1z(k)r N k=l e') = 1 . Itération (n) : Pour n =1 à m 2 N n-1 (k).bn* (k û 1) + 2.E ) a(n 1) a ,,,n:k» n-1 -1 k k - - N _ n k=1 n-1 E If _1 (k)2 + (k -1)12 + 2.E e) ,kn-1)12 1 N û n k=n+l n k=O I a,,)n) = al') = akn 1) + pn.an,ï-kle , k=l,...,n-1 an) n û ,un f,, (k) = (k) + /in .bn_1 (k -1) bn (k) = bn_1 (k -1) + p.fn_1 (k) avec jain, = y, .(27t)2.(k n)2 on prendra yl = Po.coef avec coef = 10-3 ou 10-4. REVENDICATIONS1. Procédé permettant de déterminer un paramètre central, tel que la médiane ou la moyenne, sur un ensemble de résultats d'identification associé à un signal complexe issu d'un signal radar reçu sur un ensemble de capteurs d'un récepteur radar, conduisant à P collections E1,...Ep de M- uplets de nombres complexes qui correspondent à m mesures du signal complexe, lesdites mesures étant issues du récepteur radar Ek = LXk 1,Xk 2,...Xk M J k=1,..,P ledit procédé comprenant l'application d'un algorithme d'identification auto-régressive en treillis d'ordre N-1, sur ledit signal complexe issu d'un signal radar, fournissant P points {Wl, W2 Wp} ou points d'ancrage qui sont invariants lors de la mise en oeuvre des étapes du procédé selon l'invention, caractérisé en ce qu'il comporte au moins les étapes suivantes : chaque point Wk est une constitué d'une suite de N coefficients de la forme Wk = dans laquelle Wk =1og(M EII Xk j IIZ) et (Wk,..., W,.., Wk -') sont les N-1 premiers j=1 coefficients de réflexion complexes provenant de ladite identification auto-régressive complexe par mise en oeuvre d'un l'algorithme d'identification en 20 treillis appliquée à l'ensemble des nombres, (Xk,l,•••,Xk,;,•••,Xk,M) ledit procédé s'applique sur un ensemble général appelé H dans lequel se situent lesdits points d'ancrage et ledit ensemble général H est l'ensemble des suites de la forme Wk = (W,°, Wk,..., Wk,.., Wk -') de N coefficients dont le premier coefficient est un nombre réel et les coefficients suivants des 25 nombres complexes de module strictement inférieur à 1, ledit procédé comprend également l'application d'un algorithme de descente du gradient qui opère sur une fonction définie sur l'ensemble général H, ladite fonction étant obtenue, soit comme la somme des distances des points d'ancrage à un point quelconque de l'ensemble général lorsque le procédé 5 vise à calculer une médiane, soit comme la somme des carrés des distances des points d'ancrage à un point quelconque de l'ensemble général lorsque le procédé vise à calculer une moyenne, et en ce que: ledit procédé est caractérisé en ce qu'il réalise une boucle visant à minimiser ~o la fonction ci-dessus nommée par un algorithme de gradient comprenant un nombre fini de pas au cours desquels est calculé, à chaque pas dudit algorithme, un vecteur gradient se dirigeant vers un point de convergence qui est le résultat visé par le procédé, et en ce que chaque pas de l'algorithme de gradient comporte au moins les étapes 15 suivantes qui transforment un point z de l'ensemble général H provenant du pas précédent de l'algorithme de gradient en un nouveau point z de l'ensemble général : Phase 0 : prendre le point z de l'ensemble général issu de l'algorithme à l'étape précédente j-1, pour la première étape on initialisera le point z de 20 départ en un point d'ancrage quelconque, Phase 1, (13) : effectuer sur chacun des points d'ancrage une rotation généralisée de Môbius directe de centre z pour obtenir un jeu de points P modifiés correspondant à l'étape j, W'm, Phase 2, (14) : pour chacun des points modifiés W'm, déterminer le vecteur 25 Vm tangent à H en 0, Phase 3, (15a, 15b) : déterminer pour chacun des points modifiés, le vecteur unitaire Um tangent à H en 0, Phase 4, (16) : calculer la combinaison linéaire convexe des vecteurs unitaires pondérés par un coefficient ocrr,, Phase 5, (17) : multiplier la combinaison linéaire S par un nombre réel correspondant à l'incrément du gradient à l'étape j afin d'obtenir un vecteur tangent à H en 0, Phase 6, (16) : enrouler le vecteur A sur la variété riemannienne H le long de 5 la géodésique issue de 0 pour obtenir un point Z' de H, Phase 7, (19) : effectuer la rotation inverse de Môbius de centre z pour obtenir le nouveau candidat z à l'étape j+1, Phase 8, (20) : comparer la valeur de la norme de S à une valeur seuil et si cette valeur est inférieure à la valeur seuil alors arrêter l'itération et définir 10 l'élément central (point médian ou moyenne) comme étant le dernier point déterminé, sinon, faire j=j+1 et retourner à la phase 0, ledit procédé comportant une étape au cours de laquelle connaissant le point médian ou le point moyen selon l'objet du procédé, on calcule la distance entre un point de l'ensemble H et le point médian ou la moyenne 15 issu des étapes précédentes puis on compare cette distance à une valeur seuil définie par le besoin de I 'application, et si cette distance excède le seuil on décrète la présence d'une cible. 2. Procédé selon la revendication 1 caractérisé en ce que le signal est un signal issu d'un signal radar Doppler, voie I et Q échantillonné et analysé par bloc de longueur M N au moyen de l'algorithme de BURG, l'algorithme de BURG génère un ensemble de (N-1)uplets de coefficients qui sont les N-1 premiers coefficients de réflexion complexes de l'échantillon, ledit procédé nécessite P blocs issus du récepteur qui conduisent à p points de l'ensemble H par application de l'algorithme de BURG à chacun d'eux ; chacun des P points de l'ensemble H, est un N-uplet dont la première composante (indice 0) est un nombre réel qui est le logarithme naturel de la puissance Po de l'échantillon de signal, les N-1 composantes suivantes du point de H sont les N-1 coefficients de réflexion complexes, le procédé est caractérise en ce que : La phase 1, (13) génère un jeu de P points modifiés dans H {W'1,...,W'P}propre à l'étape j en appliquant une rotation de Môbius généralisée 5 Wm = (Wm,o, Wm,1,...Wm, N ù 1) W m = (W m,o, W' m, N ù 1) = (W' m,o ù Zo, R[ZI](Wm,l),..., R[Zn ù 1] (Wm, N ù 1)) R[z](w) = WùZ 1-z.w 10 Pour chacun des points modifiés, ledit procédé calcule le vecteur Vm selon les formules. W'm,1 W'm,N-1 Vm=[W'm,0, argth(1w'm,ll),..., .arg th(IW'm,Nùll)j 1W m,11 1 W m,N-11 15 si l'objet de l'algorithme est le calcul d'un point médiane, (15a) le procédé calcule pour chacun des points Vm le vecteur unitaire Um selon les formules. 11r/m11 = 'N.Vm 02 + Ek=11(N ù k). arg th 2 (I Vm k 1) 20 Um - II Vmll si l'objet de l'algorithme est le calcul de la moyenne, (15b) le procédé recopie le vecteur Vm dans le Um Um=Vm puis le procédé calcule la combinaison linéaire convexe S des vecteurs Um pour m=1 à P. 25 50 2955950 P SûEam.U, m=1P Fan, =1 m=1 am ? 0 am nombre réel 15 On prendra par exemple : a = m 1 P 5 L'étape suivante (14) est de multiplier S par le nombre réel yi > 0 incrément du gradient à l'étape j, pour obtenir A un vecteur. 0 = 7i.S A = [Ao, 01,..., Ok,..., AN ù 1] La phase 6 est exécutée de la manière suivante : D0, Ol .th(1A11),..., ANù1 .th(I ANù1I ) IA1 ON lI Ou bien : Z~= [AI) Al 0e (Io1 I),..., ONù1 0e (ON-1 I )~ avec (x) = th(x) six et th (E°) six x E° puis d'effectuer, (19) sur Z' la rotation inverse de centre Z pour obtenir le nouveau point candidat Z à l'étape j+1. Znouveau =Tzù1(Z,) = [z'+zo,R-1[z1] (z'1 ),R-'[z2](z'2 )...,Rù1[zNù1] (z'Nù1 )] 10 Z'= i = 1,...N-1 R-1[zt](z' z'.+z. 1+z;.z' 20 avec : réitérer lesdites étapes tant que la norme du vecteur S est supérieure ou égale à un nombre réel positif proche de 0 qui est le seuil d'arrêt de la boucle. 25 51 2955950 IISII = . JN.So2 + Ek 11(N ù k). arg the (I Sk) seuil arrêt et en ce qu'il comporte les étapes suivantes après arrêt de la boucle, exécuter dans le cas du calcul pour la puissance 5 médiane (22): Puissance Médiane = exp(Zo) ou dans le cas du calcul de la puissance moyenne (23) Puissance Moyenne = exp(Zo) 3 ù Procédé selon la revendication 2 caractérisé en ce que l'on détermine le vecteur tangent Vm en appliquant Vm = [W m,O, W m,l -0V (W m,l W m,N-1'OV (W m,N-1 )] avec : 0V (x) = arg th(£V) si x £V et arg th (x) six >_ £V £V x et le calcul du vecteur unitaire Um Um = IIVmII si Vm I % £U et EU sl 4 ù Procédé selon la revendication 3 caractérisé en ce que 25 Z'- [Ao,01 0] si toutes les composantes du vecteur A sont de module petit devant 1. 10 15 20 5 û Procédé selon la revendication 1 caractérisé en ce que le signal complexe est un signal radar et en ce que le cadre de travail de l'algorithme est alors le disque D de Poincaré D = {zl z e C et lzl 0 où yi est l'incrément du gradient à l'étape j, pour obtenir le nombre A qui correspond à un nouveau point positionné sur le parcours effectué suivant le vecteur gradient obtenu. A=.S enrouler A sur la variété riemannienne D le long de la géodésique issue de 0 pour obtenir le nombre z' de D. z'= ol.th(IAI) Effectuer, (35), sur z' la rotation inverse de centre z pour obtenir le nouveau nombre candidat z à l'étape j+1. z'+z Znouveau = Tz1(z') = R [z](Z') _ 1+ Z.Z' 20 critère d'arrêt (36) arrêter la boucle lorsque le module du nombre complexe S est inférieur à un nombre réel positif proche de 0 qui est le seuil d'arrêt de la boucle. 1 SI