FR2933513A1 - Procede et dispositif d'analyse frequentielle de donnees - Google Patents

Abstract

Le procédé d'analyse fréquentielle de données comporte : - une étape (310) d'entrées de signaux issus d'un premier capteur, - une étape (315) d'entrée de signaux issus d'au moins un deuxième capteur, chaque deuxième capteur étant positionné à proximité du premier capteur pour que les signaux issus de chaque deuxième capteur soient fortement corrélés avec les signaux issus du premier capteur, - une étape d'estimation, pour chaque capteur, d'une fonction de transfert ou modèle faite à partir de l'ensemble des signaux du premier et de chaque deuxième capteur et - une étape (320) d'extraction des propriétés structurelles du système à partir de chacun des modèles estimés.

Description

10 La présente invention concerne un procédé et un dispositif d'analyse fréquentielle de données. Elle s'applique, en particulier, à l'analyse de données d'essais d'ouverture de domaine de vol des avions. La présente invention s'applique notamment, dans le domaine aéronautique, aux commandes de vols, par exemple l'analyse et le contrôle en 15 vol des modes vibratoires de la structure, dans le domaine automobile, aux études et contrôles des vibrations des véhicules, en électrodynamique (contrôle des machines productrices d'électricité), notamment dans le domaine nucléaire, à la surveillance vibratoire du coeur du réacteur, en mécanique (étude et contrôle de pièces mobiles), en sismique (étude des signaux utilisés en 20 prospection pétrolière) et en zoologie (étude des sons émis par les animaux). Le but de la présente invention est d'estimer, en cours de test (en cours de vol, dans le cas d'un avion) les caractéristiques du véhicule et, en particulier, les fréquences de résonnance et les caractéristiques spectrales. En d'autres termes, de la très grande quantité d'informations qui proviennent des 25 capteurs installés sur le véhicule, il s'agit d'extraire les signatures pertinentes très rapidement, voire en temps réel. L'ensemble des caractéristiques du système ainsi déterminées permettent aux concepteurs d'améliorer la structure du système en vue d'augmenter son confort, son domaine de vol, sa consommation, ... 30 Les signaux que l'on souhaite analyser sont composés d'un bruit auquel s'ajoutent un ou des signaux sinusoïdaux dont les fréquences et les 1 amplitudes sont susceptibles de varier au cours du temps. Il s'agit d'estimer ces fréquences et ces amplitudes en temps réel. Dans le domaine de l'aéroélasticité, discipline étudiant les interactions entre l'aérodynamique, les forces inertielles et élastiques, le phénomène de flottement est une instabilité oscillatoire de la structure de l'avion (voilure, fuselage, empennage ...) très dangereuse puisqu'elle est susceptible d'affecter l'intégrité de cette structure en l'endommageant jusqu'à rupture. C'est une combinaison de deux ou plusieurs mouvements de la structure de l'avion de natures différentes, qui, avec les déphasages appropriés, permet aux forces aérodynamiques d'apporter de l'énergie au système. Le phénomène stable est alors transformé en un phénomène instable pour lequel l'énergie n'est plus dissipée : c'est un couplage de modes. L'exemple le plus couramment cité est le flottement de voilure résultant d'un couplage entre les modes de flexion et de torsion qui, en quadrature de phase, conduisent à des forces aérodynamiques de portance de même sens que le déplacement et, ainsi, à des oscillations divergentes. Plusieurs paramètres sont influents dans la caractérisation du flottement : la masse, la raideur, la forme de la structure ainsi que les conditions opérationnelles telle que la vitesse. Afin de se prémunir de ce phénomène, les avionneurs doivent l'étudier et démontrer, s'il existe, que son seuil d'apparition est situé au-delà de la vitesse maximale de fonctionnement (augmentée de 15%). Des essais en soufflerie sont tout d'abord réalisés, complétés par des essais de vibration au sol de la structure de l'avion. Des études théoriques permettent alors de définir un domaine exempt de flottement à partir duquel l'ouverture totale du domaine de vol pourra être effectuée pas à pas en excitant la structure de l'avion. Des méthodes d'identification des paramètres modaux sont utilisées pour extraire, en quasi-temps réel, les valeurs de fréquence et d'amortissement et étudier leur évolution dans le domaine de vol. L'analyse des données temporelles issues des essais de flottement est complexe : les données sont entachées de bruit et nécessitent une mise en forme par traitement du signal (notamment filtrage et sous-échantillonnage). Plusieurs capteurs sont aujourd'hui utilisés pour extraire globalement et de manière automatique les paramètres modaux de la structure aéro-élastique. On connaît le document WO03005229 qui décrit un système d'analyse fréquentielle de signaux issus d'un capteur. Cependant la résolution de cette analyse est limitée. La présente invention vise à remédier à ces inconvénients. A cet effet, la présente invention vise, selon un premier aspect, un procédé d'analyse fréquentielle d'un système, caractérisé en ce qu'il comporte : - une étape d'entrées de signaux issus d'un premier capteur, - une étape d'entrée de signaux issus d'au moins un deuxième capteur, chaque deuxième capteur étant positionné à proximité du premier capteur pour que les signaux issus de chaque deuxième capteur soient fortement corrélés avec les signaux issus du premier capteur, - une étape d'estimation, pour chaque capteur, d'une fonction de transfert ou modèle faite à partir de l'ensemble des signaux du premier et de chaque deuxième capteur et - une étape d'extraction des propriétés structurelles du système à partir de chacun des modèles estimés. Un modèle considère le signal issu d'un capteur comme la sortie d'un filtre excité par un bruit blanc. Les propriétés structurelles comportent, par exemple, les propriétés spectrales, les fréquences, amplitudes, phases à l'origine, amortissements, modes. Ainsi le modèle représentatif des modes structuraux est considéré comme linéaire. Grâce à la mise en oeuvre de la présente invention, on effectue un traitement en temps réel en réalisant un suivi en ligne des couples fréquence/amortissement. La présente invention permet de s'assurer en temps réel que le comportement du système, par exemple de l'avion, est satisfaisant puisque l'on dispose en temps réel des propriétés structurelles du système. On améliore ainsi les méthodes d'analyse utilisées tout en répondant aux contraintes grandissantes de gain de temps et donc de réduction des coûts. Selon des caractéristiques particulières, au cours de l'étape d'estimation, on considère les signaux issus des capteurs comme des polynômes. Grâce à ces dispositions, la représentation des signaux est compacte, du fait que le nombre de coefficients des polynômes est très inférieur au nombre d'échantillons de signaux mis en oeuvre. Selon des caractéristiques particulières, au cours de l'étape d'estimation, on résout une équation récurrente linéaire à coefficients variant lentement dans le temps et, dans l'espace, entre les capteurs. Selon des caractéristiques particulières, l'étape d'estimation comporte : - une étape de modélisation adaptative récursive sur le temps, l'ordre et l'espace des capteurs et - une étape d'estimation des modes pour chaque ordre en fonction du résultat de l'étape de modélisation adaptative. Selon des caractéristiques particulières, chaque étape d'entrée de signaux issus de capteurs comporte une étape de réduction en temps réel du niveau de bruit des signaux issus de capteurs précédant l'étape de modélisation adaptative. Selon des caractéristiques particulières, l'étape d'estimation des modes comporte une étape d'extraction de paramètres du modèle en fonction du résultat de l'étape de modélisation adaptative.

Selon des caractéristiques particulières, l'étape d'extraction de paramètres du modèle comporte une étape d'inversion d'une matrice de polynômes d'ordre N et de dimension égale au nombre de capteurs. Selon des caractéristiques particulières, l'étape d'estimation des modes est adaptée à fournir les paramètres de chacun des modèles constituant un ensemble d'informations redondantes qui permet de réduire la variance des modes estimés. Selon des caractéristiques particulières, l'étape de modélisation adaptative réalise une modélisation de type ARMA ( Autorégressive à moyenne ajustée ).

Selon des caractéristiques particulière ladite modélisation de type ARMA est effectuée à chaque instant, pour chaque capteur et pour tous les ordres considérés.

Selon des caractéristiques particulières, l'étape d'estimation comporte une étape de d'inversion d'une matrice de polynôme qui est une matrice inter-spectrale symétrique représentant, sur sa diagonale principale, la densité spectrale de puissance de chacun des capteurs et, dans les autres termes, les inter-spectres. Selon des caractéristiques particulières, l'étape de modélisation adaptative comporte une récursion sur le temps, sur l'ordre du modèle pour N = [1, 2, ..., Nmax] et le nombre de capteurs des étapes suivantes : - calcul des vecteurs des erreurs de prédiction linéaire arrière et avant, - calcul des matrices de corrélation partielle avant et arrière, - calcul des matrices de covariance des erreurs de prédiction linéaire arrière et avant, - calcul de la matrice de la puissance de l'erreur de prédiction linéaire avant et arrière, - calcul direct de vecteurs gain / n y1 ea,n_ xnk-kt-kn \ k=0 avec alpha un scalaire, lambda un facteur d'oubli et 20 - calcul du vecteur 8_ n par récurrence à partir de la connaissance de ban et - calcul des matrices Ak représentant les modèles, pour k = 1 à N. Selon des caractéristiques particulières, le procédé objet de la présente invention, tel que succinctement exposé ci-dessus, comporte une 25 étape de classification des modes issus de l'étape d'estimation des modes en mettant en oeuvre l'une des deux contraintes suivantes : - un seul mode par classe issu d'un même modèle et - les estimations ont toutes le même poids indépendamment de la provenance de l'estimation. 30 Selon des caractéristiques particulières, les dits signaux sont représentatifs d'accélérations de la structure d'un avion. N N în-1 Selon un deuxième aspect, la présente invention vise un programme d'ordinateur chargeable dans un système informatique, ledit programme contenant des instructions permettant la mise en oeuvre du procédé objet de la présente invention, tel que succinctement exposé ci-dessus.

Les avantages, buts et caractéristiques particulières de ce programme étant similaires à ceux du procédé objet de la présente invention, tel que succinctement exposé ci-dessus, ils ne sont pas rappelés ici. D'autres avantages, buts et caractéristiques de la présente invention ressortiront de la description qui va suivre faite, dans un but explicatif et nullement limitatif, en regard des dessins annexés dans lesquels : - la figure 1 représente, schématiquement, un avion comportant un dispositif apte à implémenter le procédé objet de la présente invention, - la figure 2 représente des signaux issus de deux capteurs du dispositif illustré en figure 1, - la figure 3 représente, sous forme d'un logigramme, des étapes mises en oeuvre dans un premier mode de réalisation du procédé objet de la présente invention, - la figure 4 représente, sous forme d'un logigramme, des étapes mises en oeuvre dans un deuxième mode de réalisation du procédé objet de la présente invention, - la figure 5 représente un agencement de filtres mis en oeuvre au cours de l'une des étapes illustrées en figure 4, - la figure 6 représente, schématiquement, les fonctions successives mises en oeuvre dans un mode de réalisation d'un système de réduction de bruit, - la figure 7 représente, schématiquement, à chaque instant, les échantillons issus des capteurs constituant les entrées de l'algorithme du mode de réalisation illustré en figure 4, - la figure 8 représente, schématiquement des récursions mises en oeuvre dans le deuxième mode de réalisation illustré en figure 4, - la figure 9 représente, schématiquement, une évolution de classes d'une méthode de classification non-supervisée de type nuées dynamiques et - la figure 10 donne une illustration d'une fenêtre de validation mise en oeuvre dans le deuxième mode de réalisation illustré en figure 4. On observe, en figure 1, un avion 105 muni de deux capteurs proches entre eux 110 et 115 à l'avant de la voilure 120 et de deux capteurs proches entre eux 125 et 130, à l'arrière de la voilure 120. Dans un but explicatif, seuls deux couples de capteurs proches sont représentés en figure 1. On note, cependant, que, dans une implémentation réelle de la présente invention, plus de deux couples sont mis en oeuvre. Le terme de proche fait ici référence à des capteurs qui reçoivent des signaux fortement corrélés entre eux. Les capteurs proches reçoivent sensiblement les mêmes vibrations, décalées dans le temps et amorties différemment mais selon une fonction de transfert sensiblement linéaire. Les capteurs en question sont, par exemple, des accéléromètres. On observe, dans la figure 2, que le signal 205 issu d'un premier capteur d'un couple de capteurs comporte du bruit 210, et deux pics 215 et 220 et que le signal 255 issu du deuxième capteur du même couple de capteurs comporte du bruit 260 et deux pics 265 et 270. Le pic 265 correspond au pic 215 amorti et décalé dans le temps. Le pic 220 correspond au pic 270 amorti et décalé dans le temps. Comme on le comprend aisément, en mettant en oeuvre des ensembles (ici des couples) d'au moins deux capteurs proches, la présente invention permet une analyse des propriétés structurelle de l'avion. Ces propriétés structurelles comportent, par exemple, les propriétés spectrales, les fréquences, amplitudes, phases à l'origine, amortissements, modes. On observe, en figure 3, dans un premier mode de réalisation, le procédé d'analyse fréquentielle objet de la présente invention comporte, d'abord, une étape 305 de positionnement, sur la structure d'un système mécanique sujet à des vibrations, de groupes d'une pluralité de capteurs. Dans chaque groupe de capteur, au moins un capteur dit deuxième est positionné à proximité d'un capteur dit premier . Au cours du fonctionnement du système mécanique, on effectue une étape 310 d'entrées de signaux issus d'un premier capteur d'un dit groupe de capteurs et une étape 315 d'entrée de signaux issus d'au moins un deuxième capteur du même groupe de capteurs. Chaque étape d'entrée de signaux issus d'un capteur comporte une étape de réduction du niveau de bruit des signaux issus du capteur. Cette réduction de bruit peut être effectuée capteur par capteur, de manière connue ou sur un vecteur comportant, pour chacune de ses coordonnées, un signal provenant d'un capteur. Préférentiellement, cette fonction de débruitage est assurée par une décomposition sur une base d'ondelettes (Algorithme de Stéphane Mallat). Au cours des étapes 310 et 315, par exemple, pour une structure dont on recherche les fréquences propres inférieurs à 16 Hz, on échantillonne les signaux issus des capteurs à une fréquence bien supérieure, par exemple de 256 Hz en figure 6. On note que l'utilisation des ondelettes permet un traitement simple et rapide. Au cours d'une étape 320, on effectue l'extraction des fonctions de transfert par traitement des signaux issus du premier capteur et de chaque deuxième capteur. Ainsi le modèle représentatif des modes structuraux est considéré comme linéaire et on réalise un suivi en ligne des couples fréquence/amortissement sans tenir compte de l'entrée, c'est-à-dire de l'excitation injectée. Au cours de l'étape 320, on considère les signaux issus des capteurs comme des polynômes et on extrait les informations des signaux issus des capteurs, après prise en compte des modifications de phase et de valeur des dites informations entre les signaux issus des différents capteurs des différents groupes de capteurs. Au cours de l'étape 320, on réalise une étape 325 de résolution d'une équation récurrente linéaire à coefficients variant lentement dans le temps pour pouvoir estimer les modèles sur un intervalle de temps suffisamment stable. L'étape 325 comporte : - une étape 330 de modélisation adaptative récursive de type ARMA sur le temps (à chaque instant), l'ordre (pour chaque ordre considéré) et l'espace des capteurs (pour chaque capteur) et - une étape 350 d'estimation des modes pour chaque ordre en fonction du résultat de l'étape de modélisation adaptative. L'étape 330 réalise une modélisation de type paramétrique, de type ARMA (autorégresive à moyenne ajustée). Les échantillons issus des capteurs sont assemblés en un vecteur dont le nombre de composants est le nombre de capteurs considéré (voir figure 7). Par exemple, s'il y a quatre capteurs, les vecteurs considérés sont de dimension quatre. Plus généralement, dans la suite de la description, on appelle p le nombre de capteurs d'un groupe de capteurs proches. L'étape 330 effectue une modélisation dite récursive dans le temps car elle utilise les dernières estimations obtenues pour mettre à jour ses paramètres. Ainsi, dans ce mode de réalisation, la présente invention met en oeuvre des relations entre deux instants consécutifs car ils sont considérés comme corrélés et cohérents. Dans des modes de réalisation, ce sont les deux instants précédents qui servent à la récursion temporelle. L'ordre N optimal de la modélisation n'est pas nécessairement connu. Préférentiellement, on ne détermine pas cet ordre N optimal, de manière à obtenir un résultat en temps réel. Au contraire, on effectue la récursion jusqu'à un ordre Nmax suffisamment grand. On obtient ainsi un ensemble de modèles suffisamment important pour que l'information recherchée (les modes structuraux) y soit représentée.

Dans le mode de réalisation illustré en figure 3, l'étape 330 comporte une étape 335 comporte une étape de traitement d'une matrice inter-spectrale symétrique représentant, sur sa diagonale principale, la densité spectrale de puissance de chacun des capteurs et, dans les termes extra diagonaux, les inter-spectres.

L'étape 335 comporte une étape 340 de récursion sur le temps avec initialisation à l'instant n pour l'ordre N=0 comportant une récursion sur l'ordre du modèle N = [1, 2, ..., Nmax], des étapes suivantes : - calcul des vecteurs des erreurs de prédiction linéaire arrière et avant, - calcul des matrices de corrélation partielle avant et arrière, - calcul des matrices de covariance des erreurs de prédiction linéaire arrière et avant, - calcul de la matrice de la puissance de l'erreur de prédiction linéaire avant, - calcul direct de vecteurs gain, dont la dimension est le nombre de ( n Capteurs : ea, N _ n k N t N N N n ù \~ kkn an avec alpha un scalaire, lambda un facteur d'oubli et N N N - calcul du vecteurn =an-1éa,n e b n par récurrence à partir de la connaissance de a n et - calcul des matrices Ak représentant les modèles pour k = 1 à N.

A partir de l'ensemble de modèles obtenu à la fin de l'étape 330, on effectue un tri pour extraire les modes structuraux. Les modèles sont considérés comme des rapports de polynômes (modèles dits ARMA , acronyme de AutoRegressive Moving Average pour moyenne mobile autorégressive ou autorégressive à moyenne ajustée). Pour extraire les valeurs de paramètres qui représentent les modes structuraux, on minimise les matrices de covariance Ea et Eb pour chaque instant et pour chaque ordre. Pour cela, on calcule les vecteurs ea et eb qui représentent, respectivement, le produit de l'inverse de la matrice Ea par un vecteur ea et le produit de l'inverse de la matrice Eb par un vecteur Lb.

Pour réaliser ces calculs en temps réel, on calcule directement le vecteur ea,nN pour toutes les valeurs de 1 à Nmax. Puis, on calcule eb,nN par récurrence connaissant ea,nN. Le nombre de matrices Ak (carrées de dimension p), matrices coefficients des modèles à estimer à l'instant n pour tous les ordres allant de 1 30 à Nmax est égal à Nmax x (Nmax - 1)/2, soit 105 pour N = 15. Le principe de la classification des modes structuraux est de considérer que, si un paramètre à un ordre donné est pertinent, on le retrouvera à un ordre supérieur. La classification des paramètres est non supervisée et consiste à rechercher des objets semblables issus de modèles différents. Comme exposé plus loin, on met en oeuvre, pour cette classification, des constructions de trajectoires.

L'étape 350 fournit des paramètres de chacun des modèles constituant un ensemble d'informations redondantes qui permet de réduire la variance des modes estimés. L'étape 350 comporte une étape 355 d'extraction de paramètres du modèle en fonction du résultat de l'étape de modélisation adaptative. L'étape 355 comporte une étape 360 d'inversion d'une matrice de polynômes d'ordre N et de dimension égale au nombre de capteurs. L'étape 360 comporte une étape 365 de décomposition de Cholesky. Dans son deuxième mode de réalisation particulier, le procédé objet de la présente invention suit une procédure temps réel, décrite en regard de la figure 4, dédiée à l'analyse des données d'essais d'ouverture de domaine de vol. Ce mode de réalisation permet de traiter chaque information avant l'apparition de la suivante, sans tenir compte de l'excitation injectée à la structure. Les caractéristiques de la méthode répondent en totalité aux contraintes de sécurité et de réduction des coûts énoncées en préambule et permettent d'améliorer la procédure d'ouverture du domaine de vol en fournissant une analyse modale plus performante. Les signaux à analyser sont des mesures d'accélération effectuées sur la structure primaire de l'avion. Chacune des opérations de la procédure d'analyse temps réel de p signaux de type accélération, décrite en regard de la figure 4, est détaillée dans les paragraphes suivants. La première étape 405 du procédé d'analyse consiste à effectuer un débruitage des p signaux issus des p capteurs d'un même groupe de capteurs proches, par exemple à l'aide de l'algorithme pyramidal proposé par S.

Mallat, en utilisant des bases orthonormées d'ondelettes. Cet algorithme trouve son origine dans les travaux de Burt et Adelson datant de 1983 qui portaient sur la vision et la compression d'images. Cet algorithme, d'une grande simplicité de mise en oeuvre, a une charge de calcul proportionnelle au nombre d'échantillons à traiter.

Le caractère non linéaire du traitement complique la mise en oeuvre, d'autant plus que le banc de filtre n'est pas causal. Les filtres sont agencés de la façon illustrée en figure 5, en mettant en oeuvre des ondelettes.

L'algorithme de Mallat est absolument temps réel, il est généralisé au traitement simultané des p échantillons à chaque instant. Le principe retenu repose sur le sur-échantillonnage du signal afin de permettre le débruitage par analyse et synthèse.

Ainsi, au cours de l'étape 405, on effectue, d'abord, une étape 406 de décomposition, en sous-bandes, des signaux. Au cours de cette analyse, ou décomposition, la suite formée par l'échantillonnage du signal continu est considérée, dans un premier temps, comme étant l'approximation de ce signal à une certaine échelle liée à la discrétisation (l'échantillonnage correspond de fait à l'échelle la plus fine).

Par convention, cette échelle d'approximation correspond à j= 0. Nous partons donc de données qui appartiennent au sous-espace vo : la suite d'échantillons zconstitue alors l'ensemble des données que l'on souhaite analyser. La relation entre les sous-espaces d'approximation est v;_1= o w; . Il suffit donc de décomposer le signal discret, c'est-à-dire la suite d'échantillons, sur les deux sous-espaces vlet wl pour avoir la suite d'échantillons à la résolution 2-1 soit (f yrl,k) et(f 0l,k) . Les relations suivantes donnent la récurrence pour deux résolutions successives : (.Î 0 j,k) _ L hn-2k (f 0 j-1,n) et (f j,k) _ L gn-2k (f 0 j-1,n ) neZ neZ Ces équations représentent respectivement les produits de convolution g_ ,, *(f O1_l,n) et h_n *(.f suivis d'une décimation par deux. Ainsi les coefficients d'approximation (f oj,k) et de détails(f yrj,k) à l'échelle i sont calculés à partir de ceux obtenus à l'échel le j - i par une simple opération de filtrage par les filtres Ô(f)etf(f)suivie d'une décimation.

Puis, on effectue une étape 407 un seuillage des coefficients de la décomposition, au cours de laquelle, on ne garde que les coefficients de la première sous-bande. On effectue ensuite une étape 408 de synthèse, ou reconstruction, au cours de laquelle la reconstruction du signal à partir de la connaissance des projections sur les sous-espaces d'approximation est faite de la façon suivante : (f jù1,n jù1,n = L (f j,k J,k + L (f j,k I T j,k keZ keZ keZ (f jù1,n) E (f j,k j,k jù1,n) + L (f j,k )( j,k jù1,n) keZ keZ hnù2k gnù2k soit : (f 0 jù1,n) ù L hnù2k (f 0 j,k) + L gnù2k (f j,k ) keZ keZ

La reconstruction est une opération duale de la précédente. Elle est obtenue par filtrage numérique précédé d'une interpolation sur les coefficients d'approximation et de détails résultants de la décomposition. La structure générale du débruiteur est donc celle illustrée en figure 6. On note que la fréquence maximale considérée est la moitié de la fréquence f d'échantillonnage. Dans le signal brut, le signal recherché se situe dans les fréquences les plus basses. Après décomposition du signal sur des bandes de fréquences allant de 0 à f/16, de f/16 à f/8, de f/8 à f/4 et de f/4 à f/2, on effectue un seuillage par mise à zéro des coefficients des ondelettes. Puis on effectue une synthèse pour fournir un signal débruité comportant le signal recherché. Puis, au cours d'une étape 410, on effectue une modélisation adaptative de type ARMA ( AutoRegressive Moving Average ). Cette étape de modélisation du signal est de type paramétrique, type qui permet d'obtenir une analyse spectrale du signal étudié, et récursive en temps, en ordre et sur l'espace des capteurs. L'étape 410 comporte une étape 411 de détermination d'un vecteur de prédiction linéaire avant et d'un vecteur de prédiction linéaire arrière. Etendue au cas vectoriel, chaque détermination d'un vecteur de prédiction consiste à exprimer xn en une combinaison linéaire des N derniers vecteurs d'échantillons, avec la représentation illustrée en figure 7 qui représente la prédiction linéaire avant en forme vectorielle sur l'espace des p capteurs. Il s'agit donc de modéliser le signal directement dans l'espace du temps. Si xn représente la suite des échantillons, le modèle est alors :

N xn ùLAkxn-k k=1 On rappelle que le vecteurxn a pour composantes les p échantillons courants des p capteurs. Dans cette expression, lesAk sont des matrices de dimension p qui correspondent au nombre de capteur, N étant l'ordre de la modélisation.

La généralisation des méthodes scalaires donne une nouvelle écriture de l'erreur et de sa prédiction : N en = Xn ù xn en = Xn ù LAk xn-k (1) k=1 ( N Soit : 2(z) = 1ù L Akz-k x(z) _ cl)(z)x(z) (2) k=1 D'un point de vue externe, en est la sortie d'un filtre RIF excité par la suite vectorielle des échantillonsxn . La propriété de linéarité permet d'inverser le processus : xn apparaît alors comme la sortie d'un filtre excité par en. Ce filtre est obtenu en inversant la matrice polynômiale 43(z), il est stable et de type RII (acronyme de à réponse impulsionnelle infinie ). Pour la détermination d'une prédiction linéaire arrière l'estimation xn_N de Xn_N est exprimée comme une combinaison linéaire des vecteursLn-N+k}k=1,...,N . Les propriétés sont analogues à celles obtenues pour la prédiction linéaire avant. Puis, au cours d'une étape 412, on effectue le calcul des matrices de covariance et des matrices de corrélation partielle avant et arrière et des matrices des puissances des erreurs de prédiction avant et arrière. Puis, on minimise la matrice de covariance des erreurs de prédiction linéaire avant et arrière, à priori et à posteriori, construite à partir des notions de prédiction linéaire avant et arrière définies auparavant. La matrice de covariance de l'erreur de prédiction linéaire avant à posteriori est définie par: Ean =Lxn kektek =LXn-kLk t An k1//tLk t An k1 ) k=1 k=1 avec 0 1 est appelé facteur d'oubli ou d'adaptation.

Au cours d'une étape 415, on obtient des matrices coefficients du modèle d'ordre N pour toutes les valeurs de N. A cet effet, lors d'une étape 416, le vecteur4optimum est obtenu lorsque la matrice de covarianceEan est minimum, soit : aEa,n aAn Au cours d'une étape 417, un raisonnement semblable, tenu à partir de l'estimation de la matrice de covariance de l'erreur de prédiction linéaire arrière a priori donne : n

Eb,n = L Xn k 1Xk-N t Bn 1~k //t k N t Bn 1~k // k=1 Un algorithme intitulé ESA effectue, dans le cadre de l'extension multi-capteurs, le calcul des deux vecteurs définis ci-après : N r N 1-1 N 0 b,n ù 1-b,n J a,n N (EN 1-1 N e a,n+1 ù Ea,n+1 J a,n+1 Le vecteur Dan, calculé directement grâce à un l'algorithme ESA , est le produit d'une matrice par un vecteur. La matrice est la somme pondérée jusqu'à l'instant n+1du produit dyadique du vecteur eN et le vecteur est ce ùa,n même vecteur à l'instant n +1 soit N ùa,n+1 Le vecteur8 bn+1 est obtenu par récurrence à partir de 0 an . =0 Puis, au cours d'une étape 420, on effectue la résolution du problème inverse, c'est-à-dire l'extraction des paramètres des modèles. A partir des équations du modèle étendues au cas multi-capteurs (1) et (2) décrites ci-dessus : N et, = xt, ùx, et, = xt, - L, Ak xt,-k (1) k=1 / N Soit : e(z) = 1ù L Akz-k x(z) _ cl)(z)x(z) (2) k=1 N Avec cl^(z) = I - LAkz-k , matrice de polynômes en z d'ordre N = k=1

[1 ,Nmax], de dimension p x p, p étant le nombre de capteurs analysés.

e(z) est un vecteur d'erreur inconnu, par contre sa matrice de covariance est estimée. Ainsi pour la faire apparaître nous calculerons : e(z)t'e(z) = c(z)x(z)tx(z)t cl)(z) Soit la matrice inter-spectrale : x(z)tx(z)t= c(z)-1e(z)t'e(z)cI(z) Cette matrice est symétrique, sur sa diagonale principale apparaît la densité spectrale de puissance de chacun des capteurs. Les termes extra diagonaux sont les inter-spectres. L'extraction des paramètres des modèles consiste à effectuer à chaque instant et pour chaque ordre l'inversion de la matrice de polynômes d'ordre N et de dimension p : c(z) . Pour cela, la décomposition de Cholesky de la matrice e(z)t'e(z) =LtL , donne : x(z) tx(z) = .1)(z)-1 L(z)tL(z)t .1)(z)-1 =4(Z)-1 L(Z)) \~(Z) 1 L(z)) =(A(Z)-1 U(z)L(z))t (A(Z)-1 U(z)L(z))

Avec

2 Fp1 2 x(z) f x(z) = A(z)-1 tn(z)-1 Fp1 2 Fpp 2 Les polynômes F représentent les numérateurs des différentes fonctions de transfert, les dénominateurs (en fait l'unique dénominateur, d'après le théorème de superposition) sont les valeurs propres. N N n-N Posons FN(z) = ~az q F(z)tF,N(z) = L Lan-manz -m q=0 n=0 m=n aoaN alaN a2aN aNaN 0 aoaN F."(Z)`F.'Y(z)= aoaN-2 aoao aoao a,a1 aoaN z N Les coefficients du numérateur sont les racines carrés des N+1 éléments de la N+1 ième colonne de la matrice inter spectrale (sans le signe). Le traitement de la modélisation adaptative peut être représenté comme en figure 8. Dans ce schéma du traitement de la modélisation mufti- capteurs, on voit l'imbriquation des trois récursions 805, 810 et 815. Puis, au cours d'une étape 425, on effectue une classification des modes. A chaque instant, N modèles de dimension {1:p:N*p} avec p le nombre de capteurs, sont estimés. L'extraction des paramètres de chacun des modèles (valeurs de fréquence et amortissement) est un ensemble d'informations redondantes qui permet de réduire la variance de l'estimation. Nous cherchons maintenant à classer ces données issues de la modélisation au sens d'un certain critère sous les contraintes suivantes : - un seul mode par classe issu d'un même modèle, - les estimations ont toutes le même poids indépendamment de la 20 provenance de l'estimation. L'absence de connaissance a priori (la densité en particulier) a amené à privilégier une méthode de classification non-supervisée de type nuées dynamiques qui consiste à trouver les classes naturelles (implicites) pour rassembler des données non étiquetées. 25 Cette méthode, illustrée ci-après figure 9, répond de façon satisfaisante à l'ensemble des exigences selon l'algorithme suivant : - au cours d'une étape 426, on effectue un choix d'une partition initiale en K classes, - au cours d'une étape 427, on effectue une recherche systématique pour chaque donnée de la meilleure classe ; calcul de la distance des données aux barycentres et affecter l'élément à la classe dont le centre lui est le plus proche (en utilisant par exemple une distance euclidienne ou de Kullback-Leibler) et - au cours d'une étape 428, on effectue une mise à jour des barycentres de la classe émettrice et de la classe réceptrice .

Puis on retourne aux étapes 427 et 428 jusqu'à la convergence. Au cours d'une étape 430, on effectue une construction des trajectoires des modes. La problématique consiste à suivre en temps réel les trajectoires d'un ensemble de cibles correspondant aux fréquences des modes dont le nombre évolue dans le temps. La structure de l'algorithme est construite autour d'un filtre de Kalman par piste suivie. Lorsqu'un ensemble de mesures est fourni par les modèles, on cherche à les corréler aux pistes existantes. L'objectif est ici de sélectionner, parmi les mesures reçues, celles qui sont susceptibles de provenir de la cible à partir de laquelle la mesure est prédite. Un principe souvent utilisé consiste à définir une fenêtre, communément appelée gating , autour de la prédiction faite. L'architecture générale de traitement est basée sur les principes suivants : - prédiction de l'état à l'instant n+1 sachant n à partir des trajectoires 25 connues, - l'association mesures-cibles consiste à comparer les mesures avec celles prédites à partir des trajectoires connues. Ce traitement doit permettre, non seulement, d'entretenir les pistes déjà existantes mais aussi d'en initialiser de nouvelles, et éventuellement d'éliminer celles qui correspondent à des cibles 30 ayant quitté l'espace des observations. C'est de la qualité de ces fonctions que dépendent les performances du suiveur de pistes. - le filtrage avec la mise à jour de l'état, du gain de Kalman ainsi que des matrices de covariance. Le filtre de Kalman permet, dans cette procédure, de poursuivre plusieurs cibles où la prédiction tient un rôle fondamental. Il fournit, pour chaque cible, une estimation filtrée de l'état au sens du minimum de variance, prédit l'état et permet le calcul du gating . La solution est constituée par l'ensemble de deux systèmes d'équations de prédiction et de filtrage, soit : Equations de prédiction Calcul de l'état prédit z(n+ln) = Aî(nn) Calcul de la matrice de covariance de P+lln = Al nInAT +Q,, l'état prédit Calcul de la mesure prédite z(n+ln)=Hz(n+ln) Equations de filtrage Calcul de l'innovation s(n+ln) = 2(n+1)ù2(n+ ln) Calcul de la matrice de covariance de T~ T D l'innovation Sn+1 = Hl n+1I nH + 1\n Calcul du gain du filtre (gain de Kalman) Kn+~ = P+1InHT Sn+l Estimation de l'état à l'instant n+1 z(n + 1 n + 1) = z(n + 1 n) + K _ n+1z Calcul de matrice de covariance P+ln+~ = (I ûKn+~H)P+ln associée La fenêtre de validation permet, pour chaque cible, de sélectionner les mesures susceptibles d'appartenir à la cible. Le principe est de définir une zone, un volume dans l'espace des observations, autour de la mesure prédite. La taille de cette zone est définie grâce aux propriétés statistiques de la mesure prédite (Gaussienne en l'occurrence). D'une manière générale, la dimension de ce volume doit être judicieusement choisie. De lui, en effet, dépend le tri des mesures et la probabilité pour que la mesure en provenance de la cible soit à l'intérieur de la surface délimitant ce volume. La figure 10 donne une illustration de la fenêtre de validation. La technique d'association mesure-cible est la partie centrale de la procédure de suivi de cible. De nombreuses techniques existent parmi lesquelles certaines ne gèrent pas l'apparition et la disparition des pistes. Il faut donc prévoir un mécanisme supplémentaire pour réaliser cette gestion. Une façon simple consiste à adopter les règles suivantes : - règle 1 : Toute mesure qui n'est associée à aucune piste existante est considérée comme l'initialisation d'une nouvelle piste. - règle 2: Une piste est confirmée (détectée) si au moins Nd mesures consécutives lui ont été associées. - règle 3 : Une piste est considérée comme disparue si au moins NI mesures consécutives ne lui ont pas été associées.

La méthode hongroise permet de résoudre le problème d'affectation estimation-mesure par la recherche d'un coût minimal grâce à la méthode de résolution particulière suivante: Soient m ressources à affecter à m taches et soit C la matrice des coûts d'affectation. Une affectation quelconque est définie par m couples notés (1,x), (2,y), ..(k, t), ...(m, u) avec (x, y, ...u) : permutation de {1,2,...,m}. A une affectation particulière correspond une dépense ou coût total : F(x, y, ..t, ..,u) = C1,x +C2,y + ... + Ck,t +...+Cm,u Le problème consiste alors à déterminer (x, y, ..t, ..,u) de façon à rendre F minimale.

Pour la mise en oeuvre du procédé objet de la présente invention, dans un mode préféré de réalisation de la présente invention, on prévoit un ordinateur à usage général muni d'un programme d'ordinateur chargeable dans cet ordinateur, ledit programme contenant des instructions implémentant les étapes et algorithmes détaillés ci-dessus.

Claims (15)

1. REVENDICATIONS1. Procédé d'analyse fréquentielle de données, caractérisé en ce qu'il comporte : - une étape (310) d'entrées de signaux issus d'un premier capteur, - une étape (315) d'entrée de signaux issus d'au moins un deuxième capteur, chaque deuxième capteur étant positionné à proximité du premier capteur pour que les signaux issus de chaque deuxième capteur soient fortement corrélés avec les signaux issus du premier capteur, - une étape d'estimation, pour chaque capteur, d'une fonction de transfert ou modèle faite à partir de l'ensemble des signaux du premier et de chaque deuxième capteur et - une étape (320) d'extraction des propriétés structurelles du système à partir de chacun des modèles estimés.
2. 2. Procédé selon la revendication 1, caractérisé en ce que, au cours de l'étape (320) d'extraction, on considère les signaux issus des capteurs comme des polynômes.
3. 3. Procédé selon l'une quelconque des revendications 1 ou 2, caractérisé en ce que, au cours de l'étape (320) d'extraction, on résout une équation récurrente linéaire à coefficients variant lentement dans le temps et, dans l'espace, entre les capteurs.
4. 4. Procédé selon l'une quelconque des revendications 1 caractérisé en ce que l'étape (320) d'extraction comporte : - une étape de modélisation adaptative récursive sur le temps, l'ordre et l'espace des capteurs et - une étape d'estimation des modes pour chaque ordre en fonction du résultat de l'étape de modélisation adaptative.
5. 5. Procédé selon la revendication 4, caractérisé en ce que chaque étape d'entrée de signaux issus de capteurs comporte une étape de réduction en temps réel du niveau de bruit des signaux issus de capteurs précédant l'étape de modélisation adaptative.
6. 6. Procédé selon l'une quelconque des revendications 4 ou 5, caractérisé en ce que l'étape d'estimation des modes comporte une étape d'extraction de paramètres du modèle en fonction du résultat de l'étape de modélisation adaptative.
7. 7. Procédé selon la revendication 6, caractérisé en ce que, l'étape d'extraction de paramètres du modèle comporte une étape d'inversion d'une matrice de polynômes d'ordre N et de dimension égale au nombre de capteurs.
8. 8. Procédé selon l'une quelconque des revendications 4 à 7, caractérisé en ce que l'étape d'estimation des modes est adaptée à fournir les paramètres de chacun des modèles constituant un ensemble d'informations redondantes qui permet de réduire la variance des modes estimés.
9. 9. Procédé selon l'une quelconque des revendications 4 à 8, caractérisé en ce que l'étape de modélisation adaptative réalise une modélisation de type paramétrique.
10. 10. Procédé selon l'une quelconque des revendications 4 à 9, caractérisé en ce que l'étape de modélisation adaptative réalise une modélisation de type ARMA ( Autorégressive à moyenne ajustée ).
11. 11. Procédé selon la revendication 10, caractérisé en ce que ladite modélisation de type ARMA est effectuée à chaque instant, pour chaque capteur et pour tous les ordres considérés.
12. 12. Procédé selon l'une quelconque des revendications 4 à 11, caractérisé en ce qu'il comporte une étape d'inversion d'une matrice de polynôme qui est une matrice inter-spectrale symétrique représentant, sur sa diagonale principale, la densité spectrale de puissance de chacun des capteurs et, dans les autres termes, les inter-spectres.
13. 13. Procédé selon l'une quelconque des revendications 4 à 12, caractérisé en ce que l'étape de modélisation adaptative comporte une récursion sur le temps avec initialisation à l'instant n pour l'ordre N=0 comportant une récursion sur l'ordre du modèle N=[1, 2, ..., Nmax] des étapes suivantes :- calcul des vecteurs des erreurs de prédiction linéaire arrière et avant, - calcul des matrices de corrélation partielle avant et arrière, - calcul des matrices de covariance des erreurs de prédiction linéaire arrière et avant, - calcul de la matrice de la puissance de l'erreur de prédiction linéaire avant, - calcul direct de vecteurs gain n ON ù z-k YN t YN N ùa,n ù Sùk Cùk n k=0 / N N N avec alpha un scalaire, lambda un facteur d'oubli et e-ä = aä_, ea ,, - calcul du vecteur ebn par récurrence à partir de la connaissance de ean et - calcul des matrices Ak représentant les modèles, pour k = 1 à N. 15
14. 14. Procédé selon l'une quelconque des revendications 1 à 13, caractérisé en ce qu'il comporte une étape de classification des modes issus de l'étape d'estimation des modes en mettant en oeuvre l'une des deux contraintes suivantes : - un seul mode par classe issu d'un même modèle et 20 - les estimations ont toutes le même poids indépendamment de la provenance de l'estimation.
15. 15. Programme d'ordinateur chargeable dans un système informatique, ledit programme contenant des instructions permettant la mise en oeuvre du procédé selon l'une quelconque des revendications 1 à 14. 25