FR2890450A1

FR2890450A1 - Procede de determination par analyse doppler a haute resolution du champ de vitesse d'une masse d'air

Info

Publication number: FR2890450A1
Application number: FR0509095A
Authority: FR
Inventors: Christophe Chaure; Frederic Barbaresco
Original assignee: Thales SA
Current assignee: Thales SA
Priority date: 2005-09-06
Filing date: 2005-09-06
Publication date: 2007-03-09
Anticipated expiration: 2025-09-06
Also published as: RU2421754C2; RU2006132193A; US20070063887A1; US7535403B2; FR2890450B1

Abstract

L'invention se rapporte au domaine général de l'analyse doppler haute résolution d'un signal radioélectrique. Elle traite plus particulièrement de l'analyse du mouvement de masses d'air au moyen d'un radar météorologiqueL'invention a pour objet un procédé d'analyse spectrale à haute résolution d'un signal radioélectrique périodique échantillonné en distance, basé sur la mise en oeuvre d'un algorithme de filtrage autorégressif en treillis de Burg, qui permet de déterminer les fréquences propres du signal reçu pour chaque case distance à partir de la détermination d'un jeu optimum de coefficients de réflexion mun du signal. Selon l'invention les coefficients de réflexion font l'objet d'une opération de régularisation visant notamment à limiter les instabilités numériques des calculs. Selon l'invention encore, les coefficients de réflexion régularisés sont utilisés pour estimer l'ordre effectif du modèle d'identification du filtre en treillis et finalement pour déterminer les fréquences propres du signal en déterminant les arguments des racines complexes du polynôme représentant la fonction de transfert du filtre en treillis. Selon l'invention encore, on détermine pour chaque case distance celle des fréquences propres, déterminées pour cette case distance, qui représente la fréquence doppler de la masse d'air pour cette case. Le procédé selon l'invention peut avantageusement être complété par la mise en oeuvre d'une opération de lissage spatial des coefficients de réflexion avant détermination de l'ordre effectif du modèle et recherche des racines du polynôme.Le procédé selon l'invention est particulièrement adapté pour réaliser une analyse spectrale à haute résolution lorsque celle-ci doit être réalisée à partir d'un nombre faible d'échantillons de signal

Description

Procédé de détermination par analyse doppler à haute résolution du champ

de

vitesse d'une masse d'air DOMAINE DE L'INVENTION L'invention se rapporte au domaine général de l'analyse doppler haute résolution d'un signal radioélectrique. Elle traite plus particulièrement de l'analyse du 5 mouvement de masses d'air au moyen d'un radar météorologique.

CONTEXTE DE L'INVENTION - ART ANTERIEUR La connaissance précise du mouvement des masses d'air dans une zone constituant la zone située à l'avant de la trajectoire d'un aéronef revêt une grande importance. Elle permet en effet de déterminer la présence de zones de turbulence pouvant porter un grave préjudice au comportement de l'aéronef.

Cette connaissance peut être acquise par divers moyens, et en particulier par l'intermédiaire du radar météorologique qui équipe l'aéronef, ou plus généralement 15 du radar de bord qui peut être un radar multifonctions.

En l'état actuel de la technique, le radar météorologique procède à une détection des mouvements des masses d'air humide qui réfléchissent le signal radioélectrique émis et détermine leurs vitesses par analyse doppler, généralement par transformée de Fourrier. Ce principe d'analyse doppler relativement simple nécessite, pour fournir des résultats suffisamment précis, d'émettre une forme d'onde permettant de recueillir, pour une distance donnée, un grand nombre d'échantillons. Or des formes d'ondes comportant un grand nombre de récurrences ne sont généralement pas mises en oeuvre par les radars météorologiques pour lesquels la fonction de détection est effectuée à partir d'un nombre relativement restreint d'échantillons (i.e. de coups au but), de l'ordre de 8 à 16 par exemple. L'analyse spectrale par transformée de Fourrier qui en résultent ne présente alors pas une précision suffisante pour effectuer une détermination précise de zones de turbulence. Il en va de même, a fortiori, lorsque la fonction d'analyse météorologique est effectuée par un radar multifonctions pour lequel le temps alloué à la fonction météorologique est nécessairement restreint.

Ce problème de résolution se pose également en ce qui concerne les installations au sol dont l'activité consiste en particulier à détecter les mouvements des masses d'air dans la zone d'atterrissage et dans les zones d'approche. Ces équipements radioélectriques effectuent souvent une analyse doppler à haute résolution qui consiste généralement utiliser notamment une analyse autorégressive ou encore une analyse en sousespaces propres.

Les méthodes d'analyse autorégressive du signal reçu sont très généralement utilisées. Ce sont des méthodes qui permettent d'obtenir, lorsque l'on ne dispose que d'un nombre restreint d'échantillons, une résolution bien meilleure que celle obtenue par transformée de Fourrier, et qui sont par ailleurs relativement simple à implémenter. Leur principe de fonctionnement consiste à estimer le gabarit (i.e. les coefficients) d'un filtre glissant adapté au signal et présuppose que l'on détermine l'ordre du filtre que l'on cherche à définir. Une mauvaise détermination de l'ordre du filtre et de ses coefficients a en particulier pour conséquence de faire apparaître des cibles parasites.

Les méthodes standard d'analyse spectrale autorégressive à haute résolution sont des méthodes généralement récursives qui ont pour inconvénient d'engendrer des instabilités numériques compromettant la capacité à détecter les cibles, en particulier lorsque le nombre d'échantillons est faible.

Le problème de stabilité des résultats obtenus peut être résolu en utilisant des méthodes connues de l'art antérieur, comme celle développée par la déposante, qui consiste en particulier à associer un filtrage autorégressif en treillis par blocs mettant en oeuvre l'algorithme de Burg ou encore un algorithme de type "MUSIC" (Multiple Signal Classification) à une méthode dite de régularisation. Outre la fonction de filtrage, l'algorithme en treillis de Burg réalise la détermination des coefficients de réflexion Nn qui définissent le modèle d'identification appliqué au signal. A la détermination de ces coefficients est appliquée une opération de régularisation dont un des objets est de garantir leur décroissance en fonction de l'indice n. L'opération de régularisation a pour effet avantageux de rendre possible la détermination de l'ordre effectif du modèle. Néanmoins, en l'état actuel de l'art, l'opération de régularisation est généralement mise en oeuvre de manière empirique, ce qui conduit parfois à une détermination imparfaite de l'ordre du modèle, imperfection qui a finalement pour effet compromettre la résolution des mesures effectuées.

PRESENTATION DE L'INVENTION Un but de l'invention est de proposer une solution permettant d'améliorer la précision de l'analyse effectuée au moyen des méthodes autorégressives connue de l'art antérieur, et en particulier les méthodes utilisant l'algorithme de Burg, en déterminant précisément et de manière systématique l'ordre du modèle autorégressif identifié sur le signal analysé Ce but consistant in fine à obtenir une meilleure stabilité des résultats d'analyse doppler obtenus de façon à pouvoir déterminer avec une précision suffisante le champ des vitesses radiales des masses d'air de la zone d'analyse considérée.

A cet effet l'invention a pour objet un procédé d'analyse spectrale à haute résolution dont l'objet est d'établir l'estimation des composantes spectrales du signal reçu à partir d'un modèle d'identification déterminé au moyen d'un filtre autorégressif de type filtre en treillis de BURG. Ce procédé met en oeuvre les étapes principales suivantes: - une étape de filtrage du signal reçu, étape qui combine la mise en oeuvre d'un algorithme de filtrage en treillis et une opération de régularisation, cette opération permettant d'obtenir un nombre défini de jeux de coefficients de réflexion Pn, - une étape de sélection du jeu de coefficients pn optimal, - une étape de détermination de l'ordre M du modèle d'identification du signal, - une étape de déterminations des fréquences propres du signal reçu à partir des coefficients de réflexion lissés p'n et de l'ordre M du modèle.

Dans un mode de mise en oeuvre préféré, le procédé selon l'invention comporte également une étape complémentaire d'élaboration d'un jeu de coefficients p'n par lissage des coefficients Nn du jeu optimal, cette étape complémentaire prenant place entre l'étape de sélection du jeu optimal de coefficients pn et l'étape de détermination de l'ordre du modèle.

Dans un mode particulier de mise en oeuvre, avantageusement adapté à l'analyse du mouvement de masses d'air, le procédé selon l'invention comporte également une étape complémentaire qui permet de sélectionner pour chaque case distance celle des fréquences propres déterminées qui correspond à la vitesse de la masse d'air dans la case distance considérée.

Dans un mode préféré de mise en oeuvre, cette étape de sélection est réalisée en déterminant, pour chaque case distance, la puissance des fréquences propres du signal reçu sur le spectre de Capon du modèle.

Le procédé selon l'invention présente l'avantage de constituer une méthode permettant d'estimer de manière précise et systématique l'ordre du modèle d'identification élaboré par l'algorithme de filtrage en treillis utilisé en dépit du faible nombre d'échantillons disponibles pour une case distance donnée. La connaissance précise de l'ordre du modèle permet d'en déterminer plus facilement les fréquences propres.

Ce procédé présente également l'avantage de permettre, dans le cas d'une application aux radars météorologiques, la définition des lieux où se produit un changement brusque de la vitesse des masses d'air analysées. Dans ce même cas d'application il permet également, pour l'analyse de chaque case distance de tenir compte, pour la détermination des fréquences propres du signal reçu, de l'environnement constitué par les cases distance voisines. II permet en outre avantageusement de sélectionner pour chaque case distance la fréquence propre correspondant le plus vraisemblablement à la vitesse de la masse d'air.

DESCRIPTION DES FIGURES

les caractéristiques et avantages du procédé selon l'invention apparaîtront 15 clairement au travers de la description qui suit, description qui s'appuie sur les figures annexées qui représentent: - la figure 1, un organigramme de principe de l'enchaînement des étapes du procédé selon l'invention, - la figure 2, une illustration sommaire du principe de fonctionnement d'un 20 filtre en treillis mettant en oeuvre l'algorithme de Burg, - la figure 3, une illustration mettant en évidence rôle de l'étape de lissage des coefficients de réflexion.

DESCRIPTION DETAILLEE

De manière générale l'application d'un filtrage autorégressif à un signal consiste à déterminer un modèle d'identification qui soit le plus proche possible, pour l'application considérée, du signal utile, sachant que le signal reçu correspond au signal utile auquel est superposé un bruit. Les sources de bruit ont généralement diverses origines, internes au système qui effectue la réception et le traitement du signal (bruit thermique... ) , ou bien externes, ( fouillis, etc...) Dans le cas général d'un radar à impulsion le signal reçu par le récepteur, après démodulation complexe, produit un signal échantillonné sur les voies sinus (en quadrature) et cosinus (en phase), qui constituent une suite d'échantillons complexes, chaque échantillon correspondant à un petit domaine de la portée radar, encore appelé case distance. L'analyse du signal reçu dans chaque case distance, et en particulier l'analyse doppler, nécessite l'intégration des signaux engendrés par plusieurs impulsions émises consécutivement. L'environnement, en particulier les objets mobiles, et les conditions de fonctionnement du radar sont par ailleurs supposés comme sensiblement invariants durant la durée de l'analyse. La fenêtre d'analyse du signal relatif à une case distance est par exemple constituée de N = 8 à 16 échantillons accumulés durant une séquence de N émissions ou récurrences. Idéalement, pour effectuer l'analyse, on cherche à identifier le signal y(t) reçu à un signal y'(t) composé d'un nombre P, de préférence faible, de sinusoïdes complexes, c'est à dire à modéliser le signal par la relation suivante: F i2rrfjt y(t) = y (t) + B(t) = E A 1 É e + B(t) j=1 où Ai est un coefficient complexe, B(t) un bruit ( thermique, fouillis, etc...) à valeur complexe et f1 les fréquence propres (fréquences doppler) du signal.

S'agissant d'un radar fonctionnant en mode impulsionnel, le signal reçu est échantillonné au rythme de la période de répétition, ou PRI, de durée T. Par suite pour chaque case distance, le signal y(t) se présente comme une suite numérique ayant pour expression: F i2rrfÉnT yn =yin +B(nÉT)=E Ai Ée +B(nÉT) j=1 Procéder à une identification autorégressive consiste à déterminer au mieux un jeu p de p coefficients a1, ... , ap de sorte que l'écart "yn Ea1 É yn_1 " soit la plus faible j=1 possible au sens d'un certain critère qui peut être celui de l'énergie minimale (moindres carrés) de sorte que l'on puisse écrire: p yn = Eaj ' yn j + en j=1 où en est le résidu de l'estimation.

Cette formulation correspond au modèle d'identification du signal yn.

L'analyse spectrale, consiste ensuite à déterminer les fréquences fj, ce qui nécessite de déterminer la valeur de P. Les fréquences sont estimées directement à partir des phases des racines du polynôme complexe P(z-1) = 1+ i ai. z-1. [1] [2] 1=1

Le procédé décrit dans la demande a pour objet d'établir une expression de yin qui permette une analyse spectrale avec la meilleur résolution possible, cette analyse étant réalisée sur un nombre N restreint d'échantillons, N valant typiquement de 8 à 16.

On s'intéresse tout d'abord à la figure 1 qui représente une illustration du procédé selon l'invention.

Comme on peut le constater sur la figure, le procédé selon l'invention comporte cinq étapes principales. Il est appliqué pour chaque case distance au signal reçu, et détermine pour la même case distance et en fonction des cases distance environnantes, la composante spectrale correspondant au signal reçu.

La première étape, 11, consiste en l'application aux échantillons de signal reçus d'un algorithme de filtrage AR (i.e. autorégressif) en treillis par blocs mettant en oeuvre l'algorithme de Burg. Ce type d'algorithme de filtrage connu par ailleurs n'est pas développé ici. On rappel cependant ici quelques données générales relatives à cet algorithme.

Le filtrage en treillis effectue un filtrage prédictif du signal reçu par identification de ce signal à un modèle.

Le principe de fonctionnement de l'algorithme de Burg, illustré par la figure 2, repose sur un empilement de M étages de calcul 21 sensiblement identiques, qui permettent de déterminer, pour chaque étage n, l'expression de la meilleure estimation linéaire au sens des moindres carrés du signal y(t) = yo reçu à un instant t en fonction des signaux y(t-T) = y,, y(t-2T) = y2,..., y(t-nT) = yn reçus au même instant t correspondant à une case distance donnée au cours des n récurrences précédentes. Ici n est un nombre entier inférieur au nombre d'étages M. Dans le but de simplifier les notations, on utilisera dans la suite de la description pour désigner l'échantillon de signal reçu pour une récurrence donnée, la notation y(t) où t représente un nombre entier positif variant de 1 à N, nombre de récurrences disponibles pour l'analyse.

Chaque étage En, 21, calcule pour chaque valeur de t et pour une case distance 35 donnée, les fonctions d'innovation avant (En+(t)) et arrière (cn"(t)).

En+(t) caractérise l'écart existant entre y(t) et sa meilleure prédiction au sens des moindres carrés établie à partir des n échantillons précédents.

En (t) caractérise l'écart existant entre y(t-nT) et sa meilleure prédiction au sens des moindres carrés établie à partir des n échantillons suivants. On rappelle ici que l'expression des innovations avant et arrière peut être donnée pour un étage n et pour une case distance donnée, par les relations générales suivantes: En+(t) = y(t) - P(y(t) I y(t-1), y(t-2), y(t-3),... y(t-n)) et n (t) = y(t - n) - P(y(t- n) I y(t), y(t-1),y(t-2),... y(t-n+1)) [4] où "P(A I B)" représente la meilleure estimée linéaire de A à partir de B. On considère en outre ici que les fonctions n+(t) et n-(t) ont ici un caractère stationnaire et suivent donc les mêmes statistiques pour tout t entier variant de 1 à N. A partir des innovations avant et arrière calculées pour l'ensemble des N échantillons, chacun des étages En détermine le coefficient de réflexion Nn associé pour cet étage au modèle d'identification. Ces coefficients de réflexion pn,; qui caractérisent le modèle d'identification sont définis de manière connue, en fonction des innovations avant et arrière, par les relations générales suivantes:

N N

pn -- E(n-1+(t)'n 1 (t 1))/ '(E (In 1 (t-1)I2 +In-1+(t)I2) [5] t=n+1 t=n+1 où l'opérateur * désigne la conjugaison complexe et Izi le module du nombre complexe z. On peut de la sorte décrire le filtre en treillis réalisé par l'algorithme de BURG comme un système à M étages: Etage 0: entrée du signal: ò+(t) = Eo (t) = y(t), pour tout t compris entre 1 et N = initialisation de l'algorithme Etages 1 à M: n+(t) _ n-1+(t) + pn n-1 (t-1), pour tout t variant entre n+1 et N. n-(t) = + pn. n-1+(t), pour tout t variant entre n+1 et N Le signal de sortie filtré a pour expression: [3] eM(t) = M+(t) = y(t) + fin= 1 à M pn n"1"(t-1), pour t variant entre M+1 et N. Le signal filtré eM(t) se présente finalement comme une fonction des innovations arrières calculées à chacun des M étages.

Par construction, les En-(t) forment une base orthogonales de l'espace vectoriel {y(t), y(t-1), ... y(t-M)}. Les coefficients de réflexion -pn sont les coefficients de la projection de y(t) sur cette base, ils sont également des coefficients de corrélation complexes.

La transformée en z du filtre linéaire permettant de passer du signal y(t) au signal en+(t) est notée An+(z La transformée en z du filtre linéaire permettant de passer du signal y(t) au signal En" (t) est notée An"(z 1).

An+( z 1) et An"( z"1) vérifient, par transposition des formules en treillis, les relations suivantes: - pour n = 0: Ao+(z-1) = Ao-( z"1) = 1 [6] - pour n = 1 à M: n An+(z-1) = An_1+(z-1) + ln Z 1.A(z) = Ean,i ' Z-, [7] i=o An-(z 1) = z 1An-1 (z 1) + pn* An-1+(z 1) (polynôme de Szegé) [8] AM+(z 1) est un polynôme en z"1 qui est la transformée en z du filtre qui transforme ainsi le signal y(t) en eM(t). Z(EM+(t)) = AM+(z 1).Z(y(t)); où Z(A) représente la transformée en Z du signal A. Il peut également s'exprimer par la relation suivante:

M

AM+(z-1)=1+Ean z-n n=1 dans laquelle: le coefficient ao est égal à 1, le coefficient am est égal à pM. [9]

En outre, les polynômes An+(z-1) et An-(z 1) étant réciproques, le coefficient ao du polynôme Am-(z-l) est égal à pM* et le coefficient am du polynôme Am-(z-l) est égal à 1.

Les termes An+(z 1) ainsi définis forment une famille de filtres blanchisseurs puisqu'ils 5 réalisent l'écart entre y(t) et sa meilleure estimation linéaire à partir de (y(t-1), ... y(tn).

L'algorithme de Burg procède à la détermination pratique des coefficients de réflexion pn au niveau de chaque étage du filtre en treillis, en utilisant le fait que ces coefficients minimisent les fonctions Ln(pn) définies par la relation empirique suivante: Ln(pn)=N1 n É EN In+(t)I2 + 1n_(t)I2 t=n+1 Ln(pn) est une fonction de Pn qui dépend notamment de cn_1+(t) et de En_1"(t-1).

La détermination des coefficients de réflexion pn consiste alors à déterminer pour valeur de pn, la valeur qui annule la dérivée de la fonction Ln(pn) par rapport à pn. Cette condition conduit à déterminer pn à partir de la relation suivante:

N

lSm-1+(t)' m 1 (t 1).

pn _ -2, N m=n+1 2 2 [11] E Im 1+(t) + Im 1 (t -1)I m=n+1 Dans les cas classiques d'utilisation d'un filtre en treillis c'est la fonction de filtrage qui est recherchée et l'utilisateur s'intéresse plus particulièrement à l'exploitation des innovations pour obtenir le signal filtré et.

Le procédé d'analyse selon l'invention en revanche, a pour objet de déterminer les fréquences propres du signal reçu au travers du modèle d'identification. Pour ce faire, il utilise, comme l'illustre la figure 1, les valeurs des coefficients pn déterminés par chaque étage. Or du fait du faible nombre N d'échantillons dont on dispose - N est par exemple égal à 8 ou 16 la mise en oeuvre d'un filtre autorégressif en treillis classique tel qu'illustré par la figure 2 s'avère non satisfaisante. En effet il est connu que, pour peu qu'on dispose d'un nombre suffisant d'échantillons, l'algorithme d'identification produit des valeurs des coefficients p; décroissant sensiblement avec le rang de l'étage considéré. En revanche les calculs montrent que cette [10] décroissance est mal vérifiée si l'on ne dispose que d'un nombre peu important d'échantillons. Ce mauvais comportement a notamment pour conséquence qu'il est difficile de déterminer l'ordre du modèle d'identification du filtre et que l'analyse spectrale du signal au travers du modèle devient insuffisamment précise.

Un moyen connu pour pallier ce problème de mauvaise convergence consiste à appliquer, sur la détermination des coefficients Nn, une méthode de régularisation de type Thikonov qui consiste à ajouter des termes correctifs au numérateur et au dénominateur de l'expression des coefficients pn donnée par les relations [5] et [12].

Ces termes correctifs sont destinés à stabiliser le calcul des coefficients pn. Par suite les coefficients pn corrigés ont pour expression générale la relation suivante: N n 1 N2 n É IEm 1+(t) Em 1 (t 1)* +2ÉI- E bn(k)É[An-1+(k) An-1+(nk)J m=n+1 k=1 n 1 2 +1Em 1 (t-1)12 + 2ÉI'ÉE bn(k)É k=0 où: - bn(k) = (2rr)2(k-n)2' - An_1+(k) est le k-ième coefficient du polynôme An_1+(z-1) défini par la relation [7].

Le coefficient 1 présent au numérateur et au dénominateur est appelé coefficient de régularisation. La Valeur de 1 doit être déterminée de façon à optimiser cette régularisation.

D'un point de vue théorique, l'action de la régularisation consiste à contraindre les pôles de la fonction 1 /Ap+(z-1) à se placer plus prêt du zéro dans le plan complexe, ce qui a pour conséquence de rendre la fonction de transfert plus stable. Les pôles affectent moins le spectre (module de 1/ Ap+(z 1) z étant défini sur le cercle unité) de sorte que les maxima sont moins marqués. Ceci se traduit sur la série des modules des coefficients de réflexion pn par la réapparition d'une décroissance en fonction de l'ordre n de l'étage.

Le problème important posé par la régularisation réside dans l'estimation de la valeur à attribuer à r. En effet, une régularisation trop forte stabilise les coefficients pn mais éloigne le résultat obtenu par filtrage de l'identification exacte théorique, au demeurant inaccessible. Lorsque le coefficient de régularisation est trop important, 1 N E 1Em 1+ (t) N n m = n + 1 An-1+(k) [12] on constate que les pôles de la fonction de transfert AM+(z-') identifiée se répartissent sur une étoile régulière autour du zéro du plan complexe. Cette régularité traduit une perte d'information, car la position des pôles est sensée traduire les fréquences effectives des sinusoïdes contenues dans le signal. Inversement, une régularisation trop faible laisse fluctuer les coefficients de réflexion pn de manière trop importante, et ne remplit donc pas la fonction souhaitée.

Diverses propositions sont faites dans la littérature concernant l'optimisation du coefficient de régularisation. On peut notamment citer l'approche bayésienne, l'adéquation aux données, le risque moyen, la validation croisée. Ces propositions ont néanmoins un caractère plus théorique que pratique, et en l'état de l'art le coefficient de régulation est généralement déterminé de manière empirique à partir d'un certain nombre d'essais et doit être souvent réajusté, la valeur optimale de ce coefficient dépendant notamment de la puissance du signal d'entrée y(t).

En l'état de l'art le problème de l'estimation de la valeur à attribuer à r est donc mal résolu Avantageusement, le procédé selon l'invention propose une méthode de détermination du coefficient de régularisation optimum. Cette méthode selon l'invention est associée à l'algorithme de filtrage en treillis de BURG pour constituer le traitement réalisé au cours de l'étape de filtrage 11 et de l'étape suivante 12 de détermination du jeu optimum de coefficients de réflexion pn. Cette méthode s'appuie sur l'expérimentation pour définir une nouvelle approche du problème posé.

L'étude de la relation [13] exprimant les coefficients de réflexion régularisés permet de constater que cette expression ne conduit à un résultat sans dimension qu'à la condition que le coefficient r soit homogène à une puissance. En effet, une moitié des expressions du numérateur et du dénominateur est proportionnelle à la puissance du signal d'entrée Po = Ely(t)I2, tandis que l'autre moitié est une entité abstraite ne dépendant que des coefficients d'une fonction de transfert Am+(z'). C'est pourquoi la méthode selon l'invention consiste dans un premier temps à choisir un coefficient r de la forme suivante:

N

El y(t) 12 r=q. t ' N [13] où q est un coefficient sans dimension positif, inférieur à 1 et ajustable à partir 35 des données (i.e. du signal reçu dans la case distance considérée) et qui dépend seulement du rapport signal à bruit. Typiquement q prend des valeurs de l'ordre de 104.

Ainsi les formules de l'algorithme de Burg régularisé étant devenues homogènes, les coefficients de réflexion ne dépendent plus des unités et le spectre obtenu est invariant lorsque le signal d'entrée est multiplié par une constante. Les estimations des fréquences ne sont pas affectées.

Il est à noter que la valeur attribuée au coefficient q affecte fortement le spectre du modèle. Ainsi un coefficient q de l'ordre de 1 conduit à un spectre sans relief et une suite Sp de coefficients de réflexion quasi identiquement nuls. Inversement, lorsque q décroît, les coefficients de réflexion Nn croissent en module quelle que soit la valeur de l'indice n.

Pour q situé en dessous d'un certain seuil, la régularisation perd tout effet et les coefficients d'indice n élevé retrouvent des valeurs proches de 1 fluctuant aléatoirement d'un indice au suivant, en particulier pour les indices supérieurs à l'ordre effectif P du modèle d'identification recherché.

Enfin, pour des valeurs intermédiaires de q, on retrouve le comportement attendu de la suite Sp et cette suite laisse alors apparaître l'ordre effectif du modèle.

L'expression de r étant déterminée conformément à la relation [13], la méthode d'estimation de r selon l'invention consiste dans un deuxième temps à calculer un jeu de coefficients Nn pour plusieurs valeurs de q. Ces différents jeux de coefficients sont calculés en mettant en oeuvre l'algorithme de filtrage en treillis régularisé autant de fois qu'on a retenu de valeurs de q.

La méthode d'estimation de r selon l'invention se poursuit lors de l'étape 12 du procédé par l'analyse des jeux (ou séries) de coefficients ln obtenus lors de l'étape 11. Cette analyse est effectuée en calculant pour chaque valeur de q la grandeur V(q) définie par la relation suivante: M 1 V(q) =E S(l hn+1 I I Nn Dl n=p où S(I 1n+1 I I pn I) vaut I Nn+, I Nn I si la différence I Nn+1 I I Nn I est positive, et 0 si 35 elle est négative. [14]

La grandeur V(q) est appelée variation positive totale résiduelle de la queue de la série Su des coefficients de réflexion calculés pour la valeur de q considérée. La queue de la série est constituée des coefficients pn ayant un indice n entier positif compris entre un nombre entier p positif inférieur à M-1 et le nombre M-1 d'étages.

L'indice p de début de queue est choisi de façon à être supérieur au nombre de fréquences attendues dans l'analyse spectrale. Ainsi dans le cas de l'analyse spectrale effectuée dans un radar météorologique, le nombre de fréquences attendues est de l'ordre de 3 ou 4, chaque fréquence pouvant notamment correspondre à un mouvement particulier de la masse d'air. Dans un tel cas le rang p de début de queue sera par exemple pris égal à 5.

Il est à noter que lorsque q prend une valeur élevée (c'est à dire proche de 1), V(q) se rapproche sensiblement de la valeur idéale V(q) = 0. En effet comme on il a été dit précédemment, lorsque q est élevé les coefficients pn décroissent en module lorsque n croît. En revanche, une valeur trop élevée de q conduit à un écrasement du spectre du modèle. Inversement V(q) se met à croître lorsque q décroît de sorte que la régularisation n'est plus assurée. Le choix judicieux de la valeur de q peut donc être déterminé en fixant pour V(q) un seuil à ne pas dépasser.

La méthode d'estimation de 1 se poursuit ensuite en comparant chaque valeur de 20 V(q) à un seuil fixé et à retenir la plus petite valeur de q pour laquelle V(q) reste inférieure au seuil SI fixé, de préférence de valeur faible.

La valeur de q ainsi retenue, permet de déterminer la valeur optimum de r, ainsi que le jeu de coefficients pn correspondant.

Par suite, la méthode pour déterminer la valeur judicieuse de 1-peut êtdécrite de manière synthétique par l'enchaînement d'opérations suivant: détermination d'un jeu de coefficients de réflexion pn pour chaque valeur de q testée, - Calcul de la valeur de la variation positive V(q) pour chaque jeu de coefficients, - comparaison de chaque valeur V(q) calculée à un seuil SI déterminé, - Sélection de la plus faible valeur de q pour laquelle V(q) reste inférieur au seuil Si.

A titre d'exemple on peut considérer le cas d'un filtre en treillis de Burg régularisé comportant 12 étages de calcul et produisant un jeu de 12 coefficients de réflexion pl p12. L'estimation du coefficient 1 optimum peut être réalisée selon l'invention en choisissant pour q les quatre valeurs 10-2, 10"3, 10-4 et 10"5 et en calculant les jeux de coefficients p1 à p12 correspondant à chacune des valeurs choisies. On calcule ensuite les variations positives V(10-2), V(10-3), V(10-4) et V(10-5) et on compare chaque variation positive calculée à un seuil S1 déterminé. Dans le cas d'un filtre à 12 étages et lorsque le calcul de V(q) est effectué avec un indice p égal à 5 le seuil S1 est typiquement de l'ordre de 0.15.

Etant donné le petit nombre d'échantillons de signal disponible pour chaque case distance le jeu optimum de coefficient déterminé comme décrit précédemment laisse néanmoins subsister une certaine imprécision sur les coefficients eux-mêmes qu'il peut être souhaitable de corriger. A cet effet le procédé selon l'invention peut comporter une étape de lissage de ces coefficients permettant avantageusement de tirer parti de l'environnement constitué par les cases distance voisines. L'insertion de cette étape intermédiaire qui vient s'intercaler entre les étapes 12 et 13 conduit à une variante avantageuse du procédé selon l'invention. Cette étape intermédiaire de lissage, mentionnée en pointillé sur la figure 1, est décrite plus loin dans le présent document.

Comme l'illustre la figure 1, le procédé selon l'invention comporte également dans sa version la plus générale, à la suite des étapes 11 et 12, une étape 13 dont le rôle consiste à déterminer l'ordre effectif du modèle d'identification déterminé lors des étapes 11 et 12. Cet ordre effectif correspond au nombre P de coefficients de réflexion pn nécessaires pour définir de manière complète la fonction de transfert A(z1) caractérisant le modèle d'identification.

Dans le cas de l'utilisation d'un filtre en treillis classique on aboutit à une fonction de transfert qui se présente comme un polynôme en z-1 dont l'ordre correspond au nombre total M de coefficients de réflexion déterminés par le filtre. La difficulté consiste alors à déterminer ou plutôt à estimer l'ordre réel P du modèle d'identification.

En effet, les résultats théoriques connus montrent que pour un filtre de Burg portant sur un nombre important d'échantillons les coefficients de réflexion p ont un module proche de 1 pour un rang n inférieur à une valeur P qui correspond à à l'ordre effectif du modèle, puis décroissant pour n supérieur à P. En revanche pour un nombre faible d'échantillons de signal disponibles, comme c'est en particulier le cas pour les radars météorologiques, cette propriété n'est généralement pas vérifiée du fait de l'instabilité numérique qui affecte alors les calculs des coefficients pn. Dès lors, les méthodes classiques d'estimation de l'ordre du modèle, connues de l'art antérieur, telles que le critère d'information d'AKAIKE, le test MDL, le test du x2 qui fonctionnent sur les puissances a2n des résidus successifs du filtre en treillis, sont inopérantes lorsqu'elles portent sur les résultats de l'algorithme de Burg, non régularisé, appliqué à un petit nombre d'échantillons. En outre Lorsque la régularisation à lieu, si celle-ci est réalisée de manière non optimale, ces mêmes méthodes restent généralement inopérantes. De plus, même lorsqu'une régularisation optimale est réalisée, ces méthodes restent insatisfaisantes, l'erreur sur l'ordre estimé étant par exemple de 1 sur des modèles dont l'ordre effectif est 4.

En l'état de l'art le problème est donc mal résolu.

Pour résoudre ce problème, l'étape 13 du procédé selon l'invention met avantageusement en oeuvre une méthode originale qui procède par observation des coefficients de réflexions optimaux déterminés lors des étapes 11 et 12. Ces coefficients étant issus d'un filtre régularisé de manière optimale, ils se présentent comme une suite de coefficient pn dont le module est proche de 1 pour les valeurs de n inférieures à la valeur P qui représente l'ordre du modèle, et dont le module prend une valeur faible et décroissante pour les valeurs de n supérieures à P. La méthode utilisée par le procédé selon l'invention, réalise l'estimation de l'ordre du modèle à partir des coefficients pn obtenu par le filtre de Burg régularisé au moyen de la méthode avantageuse décrite précédemment.

Selon l'invention, l'observation des coefficients pn au cours de l'étape 13 est effectuée en utilisant une fonction test G(pn) qui assure un contraste permettant de réaliser de manière systématique une bonne séparation entre les coefficients pn significatifs ayant un module proche de 1 et les coefficients pn dont le module est faible. La fonction G transforme ainsi la série des modules des coefficients de réflexion (Ip,l, Ip21, Ip31, É..,IpMI) en une série de coefficients (41, 42, 43, ..., 4M) présentant un contraste élevé. Les coefficients 4n sont ensuite comparés par ordre décroissant à un seuil S2 déterminé. Par suite l'opération de détermination de l'ordre du modèle consiste alors à comparer les coefficients 4n au seuil S2 par ordre décroissant des indices, en commençant par 4m, jusqu'à trouver le premier coefficient 4K dont la valeur est supérieure à S2. L'indice P pour lequel cette condition est réalisée constitue l'estimée de l'ordre effectif du modèle.

La fonction de test G peut être définie de différentes façons, pour autant qu'elle assure un contraste satisfaisant, permettant de déterminer l'ordre P de manière systématique, à l'aide d'un simple test par rapport à un seuil S2. G(pn) peut être par exemple définie par la relation suivante: G(pn) = â É Ath2 pn [15] ou l'opérateur Ath correspond à la fonction "argument tangente hyperbolique" Le seuil S2 est déterminé en fonction de l'ordre M du filtre en treillis mis en oeuvre. Ainsi, pour M=12 le seuil S2 peut être choisi égal typiquement à 0.3.

A l'issue de l'étape 13, le procédé selon l'invention permet donc de déterminer l'expression de la fonction de transfert A(z-1), d'ordre P, qui définit le modèle 10 d'identification. Cette fonction a pour expression générale:

P

A(z-1) = Ap+(z-1) =1+ Ean.z n n=1 où P correspond à l'estimée de l'ordre effectif du modèle ainsi déterminé.

Or, il est à noter que si l'on pose z = eie, l'expression S = 1 représente le A(e-ie)IZ spectre du signal y(t), de sorte que la détermination des racines du polynôme A(z-1) permet de déterminer les fréquences propres du signal, qui sont déduites des arguments des racines de A(z-1).

Le procédé selon l'invention comporte donc une étape 14 durant laquelle est 20 effectuée la détermination des racines du polynôme A(z-1) et des composantes spectrales du signal y(t).

La détermination de racines de A(z-1) peut être réalisée de différentes façons, sachant que pour un polynôme d'ordre supérieur à 3, il n'existe pas de méthode analytique permettant de déterminer les racines de ce polynôme. Passé cet ordre on a donc généralement recours à une détermination numérique par approximations successives.

Cependant, dans le cas de racines appartenant au plan complexe, la détermination par approximations successives est rendue difficile du fait de la dispersion possible des racines.

Pour lever cette difficulté le procédé selon l'invention procède avantageusement à la substitution au polynôme A(z-1) d'un polynôme A'(z 1) dont les racines ont des arguments très proches de celles de A(z 1) mais sont situées sur le cercle unité. Des calculs menés par ailleurs et non développés ici, montrent qu'un tel polynôme peut être déterminé en construisant un polynôme A'(z1) ayant pour expression la relation suivante: [16] A' (Z-1) = A p+(z-1) = Ap-1+(z-1) +p. Z 1É Ap 1(z-1) l0p l Comme on peut le remarquer, le coefficient de z -P qui apparaît quand on développe 5 le polynôme de la relation [17] présente avantageusement un module égal à 1, ce qui simplifie la détermination des racines.

La recherche des racines du polynôme A'(z-1) est avantageusement plus simple que celle des racines de A(z-1).

Le polynôme A'(z 1) ainsi obtenu peut en outre être traité au moyen d'un outil de transformation de type transformée de Cayley de façon à obtenir une fonction rationnelle à valeurs réelles Q(z-1) image de A'(z 1) dont les racines égales à celles de A'(z-1) et situées sur le segment de droite compris entre rr et +u, sont aisément déterminables notamment par dichotomie.

Le recours à l'utilisation d'un polynôme A'(z-1) dont les racines sont situées sur le cercle unité, et à la transformation de Cayley permet ainsi avantageusement de déterminer les fréquences propres du signal analysé à partir du modèle d'identification établi à l'issue des étapes 11 à 13.

Ainsi, comme l'illustre la figure 1, l'enchaînement des étapes 11 à 14, le procédé selon l'invention permet d'effectuer une analyse du signal reçu qui procède par identification du signal à un modèle, et détermination des composantes spectrales des fréquences propres dudit modèle. Ce procédé permet avantageusement d'effectuer une analyse spectrale avec une résolution supérieure à celle espérée avec un filtrage doppler classique de type FFT, à partir d'un nombre N relativement faible d'échantillons de signal.

La mise en oeuvre du procédé selon l'invention est en outre caractérisée par la méthode employée pour le calcul du coefficient r de régularisation, ainsi que par la méthode employée pour déterminer l'ordre effectif du modèle. Le procédé est encore caractérisé par la méthode de substitution utilisée pour faciliter la recherche des arguments des racines du polynôme A(z-1), arguments qui représentent les fréquences propres du modèle. [17]

Grâce au procédé selon l'invention on obtient ainsi pour chaque case distance P fréquences propres correspondant aux P racines du polynôme A(z1) et caractérisant le spectre du signal reçu.

Le procédé d'analyse spectrale selon l'invention, tel qu'il est décrit dans ce qui précède peut bien entendu être mis en oeuvre tel quel pour déterminer pour chaque case distance, la ou les fréquences propres du signal reçu. Un traitement complémentaire peut ensuite être appliqué pour garder la fréquence la plus significative. On peut encore considérer chaque fréquence comme représentative d'un objet détecté et traiter séparément chacune des fréquences propres Dans le cas d'un radar météorologique, l'analyse doppler à haute résolution vise à produire une image du champ de vitesse radiale des masses atmosphériques mettant en évidence les zones de perturbation. L'information tirée du signal issu d'une case distance donné est fortement bruitée, en raison notamment de l'élargissement du spectre lié au phénomène météorologique observé (masse d'air en turbulent), de l'imprécision sur les statistiques liée au faible nombre d'échantillons disponible pour chaque case, de la présence d'un éventuel écho de sol et du bruit propre du radar. Cette dégradation du signal reçu se traduit généralement par la détermination pour chaque case distance de plusieurs fréquences propres, le nombre de fréquences propres ainsi déterminées n'étant pas nécessairement identique d'une case distance à une autre voisine de la première. Néanmoins, Il est peu pertinent de modéliser le signal reçu,par une collection de fréquences propres bien séparées. Le spectre réel est un spectre large qui peut être vu comme correspondant à une collection de fréquences très proches groupées en paquet avec une fréquence principale centrale qui en marque le sommet. C'est cette fréquence principale qu'il importe de déterminer, car c'est elle qui représente effectivement la fréquence doppler, c'est à dire la vitesse, de la masse atmosphérique considérée. A cet effet le procédé selon l'invention peut être complété, notamment dans le cas d'un radar météorologique, par une étape supplémentaire 16 consistant à déterminer parmi les fréquences propres, la fréquence principale correspondant effectivement à la vitesse de la masse atmosphérique.

Selon l'invention la recherche de cette fréquence est réalisée en calculant le niveau le plus fort sur le spectre de Capon du signal. On rappelle ici que le spectre de Capon du signal reçu, pour une case distance donnée, a pour expression: s(e)= R 1 [18] 1 2 I An+(e ie) I2 n=o Qn avec P R M où An+(e-ie) représente la fonction An+(z-') définie par la relation [7], avec z = eie. On rappelle que 8 a pour expression: 0= 2rr É fi É T où fi représente la fréquence propre du signal qui est proportionnelle à la vitesse de la masse d'air pour la case distance i considérée.

Selon l'invention, S(8) est calculé pour toutes les valeurs de e correspondant aux fréquences propres. Ces valeurs sont ensuite comparées de façon à déterminer la fréquence donnant la valeur la plus élevée. Compte tenu du fait que, comme cela a été dit précédemment, le spectre du signal reçu est un spectre large avec une fréquence principale centrale qui forme un sommet, cette fréquence propre particulière est choisie pour représenter la valeur du champ de vitesses pour la case distance considérée.

II est à noter que les coefficients an sont déterminés par les relations de récurrence suivantes: ÉQO=1 Éan2 =Qn-12(1 IIJn 12) Ainsi, à l'issue de l'étape 16, l'analyse spectrale réalisée au moyen du procédé selon l'invention ne laisse plus subsister qu'une seule fréquence propre par case distance, fréquence propre qui dans le cas d'une application tournée l'analyse météorologique correspond à la valeur du champ de vitesse pour la case distance considérée.

Le procédé selon l'invention tel qu'il est décrit dans ce qui précède comporte un certain nombre de caractéristiques avantageuses qui permettent d'effectuer une analyse spectrale à haute résolution, de manière efficace, performante et systématique, en utilisant un traitement mettant en oeuvre un algorithme de filtrage en treillis et en particulier l'algorithme de Burg. Il est cependant possible d'accroître encore les performances du procédé selon l'invention en effectuant un traitement complémentaire sur le jeu optimal de coefficients Pr, à l'issue de l'étape 12. Cette possibilité est en particulier offerte lorsque le signal y(t) reçu pour une case distance donnée n'est pas totalement indépendant du signal reçu pour les cases distantes voisines. Dans de telles conditions, une variante du procédé selon l'invention peut [19] avantageusement être mise en oeuvre. Cette variante de mise en oeuvre intercale entre les étapes 12 et 13 du procédé initial une opération intermédiaire 15 dont le but est d'effectuer pour chaque case distance un lissage du jeu optimum de coefficient pn, déterminé à l'issue de l'étape 12, en fonction des jeux de coefficients des cases voisines. Le principe de ce lissage est illustré par la figure 3 qui représente de façon schématique une évolution possible, au fil des cases distance, de la fréquence propre du signal reçu, correspondant par exemple à une masse atmosphérique. Dans cet exemple, la masse d'air considérée présente un champ de vitesse ou de fréquence comportant deux zones 31 et 32 dans lesquelles les variations de fréquence propre sont relativement faibles. Dans ces zones cette fréquence prend suivant la case distance considérée une valeur qui fluctue autour d'une valeur moyenne, la fréquence FI pour la zone 31, et la fréquence F2 pour la zone 32. Le champ de vitesse correspondant à l'exemple de la figure 3 comporte également un certain nombre de cases distances formant une zone de transition 33, et pour lesquelles la fréquence propre associée varie de manière importante.

Dans la réalité cette zone de transition peut par exemple correspondre à une zone de turbulence, dont il est utile de déterminer la position et l'étendue. L'objet de l'opération de lissage proposée par le procédé selon l'invention est d'atténuer les disparités entre les coefficients réflexion relatifs aux cases distance situées dans une même zone 31 ou 32 sans pour autant estomper les contours de chaque zone, contours qui marque la position des zones de turbulence.

Selon l'invention, l'opération de lissage est réalisée en appliquant à chacun des coefficients pn du jeu optimum une opération de filtrage gaussien, anisotrope et non linéaire qui consiste, pour chaque case distance r, à remplacer chaque coefficient pn du jeu optimum par un coefficient p'n dont la valeur résulte d'une combinaison linéaire des coefficients pn des cases distance voisines. Cette combinaison linéaire est réalisée sur un intervalle de cases distance centré sur la case r et de demi largueur L (L à titre d'exemple pourra être fixé à 10) de sorte que la fenêtre des indices des cases distance est comprise entre sup(0, r L) et inf(Lmax, r + L) où Lmax est le nombre maximum de cases distances disponibles. Cette combinaison linéaire se présente comme la somme des coefficients pn pondérés par des poids dont la forme rappelle la loi de gauss et ayant pour expression, pour la case i considérée: d(r,i)2 Wr,i = e Q2 [18] où a est un paramètre ajustable du type variance qui contrôle la largeur de la gaussienne, et où d(r, i) représente une distance curviligne séparant la case distance r de référence, pour laquelle on calcule p'n, de la case distance voisine i considérée. La distance d(r, i) est curviligne en ce sens qu'elle est la somme des distances entre cases distance adjacentes présentes entre r et i.

d(r, i) = somme de r à i-1 de d(j, j+1) si i > r et somme de i a r-1 de d(j,j+1) si i < r Selon l'invention, cette distance élémentaire d(j, j+1) que l'on notera simplement d a pour expression la relation suivante: d2 =Ox2 +Dy2 Ax est une grandeur proportionnelle à la distance séparant les deux cases distances adjacentes,donc proportionnelle à la largeur de la case distance. Ax est défini par la relation suivante: Ax2 =(3É(dj-dj+, )2 où 13 est un facteur d'échelle positif ajustable, (d; - dr) représentant la distance radiale séparant les deux cases i et r.

Ay est une distance spectrale définie, au moyen des coefficients pn des cases j et k=j+1 considérées, par la relation suivante: 2 M 2 oy2 = M log Puk + 1(M - n)[z ' arg th(' Cn I)] Puj n=1 Nn,j Nn,k Pn,j Pn,j+1 avec: Cn = _ [22] 1 [Pn,k ÉNn,j] 1-[Nn,j+l 'hn,j] Les grandeurs Puj et Puk représentent ici la puissance du signal reçu respectivement pour la case j et la case k, chaque case distance j étant conventionnellement représentée par un vecteur à M+1 composantes (Puj, P2,j,...,PM,j).

Les puissances Puj et Puk sont par ailleurs définies par la relation suivante: [19] [20] [21] Pu = [23] où y(t) représente ici le signal reçu dans la case distance considérée.

Ay représente une distance géodésique au sens de théorie géométrique de l'information de Chentsov entre spectres identifiés par l'algorithme de Burg sous une hypothèse de signal d 'entrée gaussien (géométrie Riemannienne de l'information avec pour tenseur métrique la matrice de Fisher sur la variété complexe des paramètres Pu et p) A l'issue de l'étape de lissage 15, on obtient un nouveau jeu de coefficients p'n tiré du jeu optimum de coefficients pn qui est utilisé pour les étapes 13, 14, et 16 qui suivent. Il est à noter que la méthode de lissage selon l'invention prend avantageusement pour critère que plus une case distance i est géographiquement ou spectralement distantes de la case r considérée et moins le coefficient de réflexion p; correspondant intervient dans la valeur du coefficient pr lissé.

Le caractère non linéaire du lissage effectué et son caractère anisotrope sont liés à l'emploi d'une loi de pondération qui dépend pour les indices r et i de la distance curviligne totale entre r et i.

Avantageusement, la loi de pondération W;,r des coefficients, utilisée ici, permet d'opérer un lissage efficace des coefficients p. Ainsi, pour les cases distance situées dans une zone où la masse atmosphérique possède une vitesse sensiblement constante la distance curviligne totale qui inclut la distance spectrale entre deux cases est relativement faible et la valeur lissée du coefficient de réflexion pn d'une case distance r donnée appartenant à cette zone sera le résultat d'un moyennage de valeurs de coefficients p; de l'ensemble des cases environnantes. En revanche, dans les zones 35 proches de la zone de transition 33 la valeur lissée du coefficient de réflexion d'une case r sera seulement le résultat du moyennage des valeurs des coefficients p; des cases distance situées du même côté de la zone de transition que la case r considérée, la contribution des cases géographiquement voisines mais situées de l'autre côté de la zone de transition, c'est à dire spectralement lointaine, étant très faible.

Claims

REVENDICATIONS

1. Procédé d'analyse spectrale à haute résolution d'un signal radioélectrique périodique échantillonné en distance, permettant de déterminer pour chaque case distance les fréquences propres du signal reçu à partir d'un modèle d'identification autorégressif, caractérisé en ce qu'il comporte au moins: - une étape de filtrage (11) mettant en oeuvre un algorithme de filtrage autorégressif en treillis comportant M étages de calculs, permettant de calculer à partir de N échantillons de signal, un jeu de M coefficients de réflexion Nn caractéristiques du signal reçu, l'algorithme de filtrage étant combiné à une opération de régularisation des coefficients calculés, qui rend l'expression des coefficients de réflexion dépendant d'un coefficient de régularisation r; l'algorithme de filtrage en treillis étant mis en oeuvre pour différentes valeurs choisies de r au cours de l'étape de filtrage, de façon à obtenir plusieurs jeux de coefficients pn correspondant chacun à une des valeurs choisies del- - une étape de détermination (12) du jeu optimum de coefficients de réflexion, parmi les jeux de coefficients déterminés lors de l'étape de filtrage; -une étape de détermination (13) de l'ordre effectif P du modèle d'identification du filtre, à partir du jeu optimum de M coefficients pn en analysant les coefficients Nn au travers d'une fonction de contraste G(pn); - une étape de détermination (14) des fréquences propres du modèle d'identification à partir du polynôme A(z') caractérisant la fonction de transfert du filtre réalisé par l'algorithme de filtrage

2. Procédé selon la revendication 1, dans lequel l'étape de filtrage met en oeuvre un algorithme en treillis de Burg.

3. Procédé selon l'une des revendications 1 à 2 dans lequel, le coefficient 1 de régularisation a pour expression:

N

El y(t) 12 T' = q. " N Où y(t) représente le signal reçu pour une case distance donnée, N le nombre d'échantillon et q un coefficient sans dimension positif, inférieur à 1, ajustable et qui dépend seulement du rapport signal à bruit.

4. Procédé selon l'une des revendications 1 à 3, dans lequel le jeu optimum de coefficients de réflexion est déterminé en analysant pour chaque jeu la valeur de la grandeur V(q) définie par la relation suivante: M 1 V(q)=Es(Ipn+1 I Ipn I)1 n=p où S(I pn+1 I I pn 1) vaut I Nn+1 1 I pn I si la différence I pn+1 I I pn I est positive, et 0 si elle est négative, pn+1 et pn correspondant respectivement aux coefficients de réflexion calculé par les étage n et n+1 les indices n des coefficients étant compris entre un nombre entier p positif inférieur à M-1 et le nombre M-1, M représentant le nombre total d'étages du filtre, p étant choisi de façon à être supérieur au nombre de fréquences attendues dans l'analyse spectrale.

5. procédé selon l'une quelconque des revendications précédentes, dans 20 lequel la fonction de contraste G(pn) a pour expression: G(pn)= 4É Ath2 1 pn 1É

6. Procédé selon l'une quelconque des revendications précédentes, l'étape 25 de détermination des fréquences propres du modèle d'identification est réalisée en déterminant les racines du polynôme A'(z 1) ayant pour expression: A'(z l) =Ap+(Z 1)= A p-1+(z1)+ PP.z-1.Ap-1 (Z-1) I0PI où P correspond à l'estimée de l'ordre effectif du modèle.

7. Procédé selon l'une de revendications précédentes, comportant en outre une étape complémentaire (16) réalisant, pour chaque case distance, la détermination de celle des fréquences propres préalablement déterminées qui constitue la fréquence principale, correspondant à la valeur du champ de vitesse dans la case considérée.

8. Procédé selon la revendication 7, dans lequel la fréquence principale est déterminée en recherchant la fréquence propre fj qui donne la valeur la plus forte du spectre de Capon défini par: S(A) = R E 2 An+(e ie) 12 n=0 a n où 8 = 2rrTfi correspond la fréquences propre du signal.

9. Procédé selon l'une quelconque des revendications précédentes comportant en outre une étape intermédiaire (15) de lissage des coefficients de réflexion formant le jeu optimum, ce lissage consistant à effectuer une opération de filtrage gaussien, anisotrope et non linéaire dont la fonction est de remplacer chaque coefficient pn du jeu optimum par un coefficient p'n dont la valeur résulte d'une combinaison linéaire des coefficients Nn des cases distance voisines.

10. Procédé selon la revendication 9, dans lequel la combinaison linéaire se présente comme la somme des coefficients pn des cases voisines pondérés par des poids dont la forme rappelle la loi de gauss et ayant pour expression, pour la case i considérée: Wri e a2 où a est un paramètre ajustable du type variance qui contrôle la largeur de la gaussienne, et où d(r, i) représente une distance curviligne séparant la case distance r de référence, pour laquelle on calcule p'n, de la case distance voisine i considérée

11. Procédé selon la revendication 10, dans lequel la distance d(r,i) est de la 30 forme: d2 =Ax2 +4y2 où Ax et Ay sont définis par les relations suivantes: Pn,j Pn,j+1 avec: C n n 1 [Nn,j+1* Pn,j J et Ax2 =PÉ(dj dj+1)2 1 2 M 2 Pu oy2 = M. log +1 + E(M n) [2 É arg th(I Cn I)] Puj n=1 les indices j et j+1 représentant deux cases adjacentes.