Login| Sign Up| Help| Contact|

Patent Searching and Data


Title:
METHOD, DEVICE AND COMPUTER PROGRAM FOR MONITORING A ROTATING MACHINE OF AN AIRCRAFT
Document Type and Number:
WIPO Patent Application WO/2021/089936
Kind Code:
A1
Abstract:
The invention relates to a method (1) for monitoring a rotating machine (100) of an aircraft, wherein a measurement signal is acquired from the rotating machine. According to the invention, instantaneous frequencies (ƒĸ(t)) of sinusoidal components of the measurement signal are estimated, and, using a computing module (12), a plurality of successive iterations are carried out in each of which: • complex envelopes of the components are updated (C1), • parameters of a model of a noise of the signal are updated (C21) on the basis of the envelopes, • whether the model has converged from the preceding iteration to the present iteration is tested (C4), with a view to: o if not, carrying out a new iteration, o if so, performing a computation (D) of the complex envelopes on the basis of the iterations that have been carried out.

Inventors:
MARNISSI YOSRA (FR)
ABBOUD DANY (FR)
EL BADOUI MOHAMED (FR)
Application Number:
PCT/FR2020/051978
Publication Date:
May 14, 2021
Filing Date:
November 03, 2020
Export Citation:
Click for automatic bibliography generation   Help
Assignee:
SAFRAN (FR)
International Classes:
G01R23/16; G01H1/00; G01M1/08; G01M1/22; G01M13/045; G01M15/12; G01M15/14; G01R23/167; G01R31/00
Foreign References:
EP2085902A12009-08-05
CN107300467B2019-09-06
US20040125893A12004-07-01
US20150090036A12015-04-02
CN104316323B2017-01-11
FR3079976A12019-10-11
US3307408A1967-03-07
US3938394A1976-02-17
Other References:
PAN ET AL: "Adaptive Vold-Kalman filtering order tracking", MECHANICAL SYSTEMS AND SIGNAL PROCESSING, ELSEVIER, AMSTERDAM, NL, vol. 21, no. 8, 28 September 2007 (2007-09-28), pages 2957 - 2969, XP022277316, ISSN: 0888-3270, DOI: 10.1016/J.YMSSP.2007.06.002
MIN-CHUN PAN ET AL: "Adaptive Angular-Displacement Vold-Kalman Order Tracking", 2007 IEEE INTERNATIONAL CONFERENCE ON ACOUSTICS, SPEECH, AND SIGNAL PROCESSING 15-20 APRIL 2007 HONOLULU, HI, USA, IEEE, PISCATAWAY, NJ, USA, 15 April 2007 (2007-04-15), pages III - 1293, XP031463714, ISBN: 978-1-4244-0727-9
HONG JIANG ET AL: "Application of spectral kurtosis and vold-Kalman filter based order tracking in wind turbine gearbox fault diagnosis", 2017 9TH INTERNATIONAL CONFERENCE ON MODELLING, IDENTIFICATION AND CONTROL (ICMIC), IEEE, 10 July 2017 (2017-07-10), pages 49 - 53, XP033332362, DOI: 10.1109/ICMIC.2017.8321698
YU GUO ET AL: "Noise removal in Vold-Kalman order tracking based on independent component analysis", AUTOMATION AND LOGISTICS, 2008. ICAL 2008. IEEE INTERNATIONAL CONFERENCE ON, IEEE, PISCATAWAY, NJ, USA, 1 September 2008 (2008-09-01), pages 574 - 579, XP031329708, ISBN: 978-1-4244-2502-0
KANDUKURI SURYA TEJA ET AL: "Diagnostics of stator winding failures in wind turbine pitch motors using Vold-Kalman filter", IECON 2019 - 45TH ANNUAL CONFERENCE OF THE IEEE INDUSTRIAL ELECTRONICS SOCIETY, IEEE, vol. 1, 14 October 2019 (2019-10-14), pages 5992 - 5997, XP033669821, DOI: 10.1109/IECON.2019.8926983
Attorney, Agent or Firm:
REGIMBEAU (FR)
Download PDF:
Claims:
REVENDICATIONS

1. Procédé de surveillance d'au moins une machine tournante (100) d'un aéronef, procédé dans lequel

- on acquiert par un module (10) d' acquisition au moins un signal (y (t)) de mesure de la machine tournante, caractérisé en ce que

- on estime des fréquences (fk (t) ) instantanées de composantes sinusoïdales ( ak(n )) du signal (y(t)) de mesure par un module (11) d'estimation, - on exécute par un module (12) de calcul plusieurs itérations successives, dans chacune desquelles :

• on met à jour (Cl) des enveloppes complexes ( ak(n )) des composantes sinusoïdales,

• on met à jour (C21) des paramètres d'un modèle d'un bruit ( e(n ), l(n)) du signal (y(t)) de mesure à partir des enveloppes complexes ( ak(n )) mises à jour,

• on teste (C4) si le modèle converge de l'itération précédente à l'itération actuelle, pour :

¨ dans la négative, refaire une nouvelle itération, ¨ dans l'affirmative, effectuer un calcul (D) des enveloppes complexes ( ak(n )) à partir des itérations effectuées, à chaque itération, on met à jour (C22) des paramètres (yk) de régularisation du modèle par le module (12) de calcul à partir des enveloppes complexes ( ak ( n )) mises à jour, et, à chaque itération, on met à jour (C3) par le module (12) de calcul (12) une variable auxiliaire (v) de découplage des enveloppes complexes ( ak(n )) des composantes entre elles dans le modèle, à partir des enveloppes complexes ( ak(n )) mises à jour et des paramètres du modèle du bruit mis à jour.

2. Procédé suivant la revendication 1, caractérisé en ce que le modèle du bruit (e(n), l(n)) du signal (y(t)) de mesure comporte un bruit interfèrent (e(n)) structuré, dont la distribution de probabilité a priori de son spectre a une loi à queue lourde calculée par le module (12) de calcul à chaque itération.

3. Procédé suivant l'une quelconque des revendications précédentes, caractérisé en ce que le modèle du bruit (e(n), l(n)) du signal (y(t)) de mesure comporte un bruit interfèrent (l(n)) non structuré dont la distribution de probabilité a priori est donnée par une troisième distribution du type Gaussien calculée par le module (12) de calcul à chaque itération.

4. Procédé suivant l'une quelconque des revendications précédentes, caractérisé en ce que le modèle du bruit (e(n), l(n)) du signal (y(t)) de mesure résultant d'un modèle (e(n)) de bruit interfèrent structuré et d'un modèle de bruit interfèrent (l(n)) non structuré est la convolution d'une loi Gaussienne et d'une loi à queue lourde dont les paramètres sont calculés par le module (12) de calcul à chaque itération.

5. Procédé suivant l'une quelconque des revendications précédentes, caractérisé en ce qu'une loi de probabilité a priori conjuguée est affectée à des hyperparamétres des lois a priori des paramètres du modèle du bruit.

6. Procédé suivant l'une quelconque des revendications précédentes, caractérisé en ce qu'une loi de probabilité a priori de Gamma est affectée à au moins un des paramètres (yk) de régularisation et résulte d'une loi a posteriori de Gamma dont les paramètres sont calculés par le module (12) de calcul à chaque itération.

7. Procédé suivant l'une quelconque des revendications précédentes, caractérisé en ce que la variable auxiliaire (v) de découplage est définie par un modèle ayant une distribution de probabilité du type loi normale complexe circulairement symétrique qui est calculé par le module (12) de calcul à chaque itération.

8. Procédé suivant l'une quelconque des revendications précédentes, caractérisé en ce que le test de convergence comprend de tester à chaque itération actuelle par le module (12) de calcul que la différence entre la valeur absolue de l'enveloppe (ak(n)) calculée de l'itération actuelle et la valeur absolue de l'enveloppe (ak(n)) calculée de l'itération précédente est inférieure à un seuil prescrit.

9. Procédé suivant l'une quelconque des revendications précédentes, caractérisé en ce que l'on analyse par le module (12) de calcul les enveloppes complexes (ak(n)) ayant été calculées, et on déclenche par le module (12) de calcul la communication d'une alarme sur au moins une sortie physique (SP) en fonction des enveloppes complexes (ak(n)) ayant été analysées.

10. Dispositif (1) de surveillance d'au moins une machine tournante (100) d'un aéronef pour la mise en œuvre du procédé suivant l'une quelconque des revendications précédentes, le dispositif comportant

- un module (10) d'acquisition d'au moins un signal (y(t)) de mesure de la machine tournante, caractérisé en ce que le dispositif comporte - un module (11) d'estimation de fréquences (fk (t) ) instantanées de composantes sinusoïdales ( ak(n )) du signal (y(t)) de mesure,

- un module (12) de calcul configuré pour, à chacune de plusieurs itérations successives effectuées:

• mettre à jour (Cl) des enveloppes complexes ( ak(n )) des composantes sinusoïdales,

• mettre à jour (C21) des paramètres d'un modèle d'un bruit (e(n), l(n)) du signal (y(t)) de mesure à partir des enveloppes complexes ( ak(n )) mises à jour,

• tester (C4) si le modèle converge de l'itération précédente à l'itération actuelle, pour :

¨ dans la négative, refaire une nouvelle itération,

¨ dans l'affirmative, effectuer un calcul (D) des enveloppes complexes ( ak(n )) à partir des itérations effectuées, le module (12) de calcul étant configuré pour, à chacune des itérations successives effectuées:

- mettre à jour (C22) des paramètres (yk) de régularisation du modèle à partir des enveloppes complexes ( ak(n )) mises à jour, - mettre à jour (C3) une variable auxiliaire (v) de découplage des enveloppes complexes ( ak(n )) des composantes entre elles dans le modèle, à partir des enveloppes complexes ( ak(n )) mises à jour et des paramètres du modèle du bruit mis à jour.

11. Dispositif de surveillance suivant la revendication 10, caractérisé en ce que le module (12) de calcul est configuré pour analyser les enveloppes complexes (ak(n)) ayant été calculées, et est apte à déclencher la communication d'une alarme sur au moins une sortie physique (SP) du dispositif (1) en fonction des enveloppes complexes ( ak(n )) ayant été analysées.

12. Dispositif de surveillance suivant la revendication 10 ou 11, caractérisé en ce que la variable auxiliaire (v) de découplage est définie par un modèle ayant une distribution de probabilité du type loi normale complexe circulairement symétrique qui est calculé par le module (12) de calcul à chaque itération.

13. Programme d'ordinateur comportant des instructions de code pour la mise en œuvre du procédé de surveillance suivant l'une quelconque des revendications

1 à 9, lorsqu'il est exécuté sur un ordinateur.

14. Aéronef comportant au moins une machine tournante (100) et un dispositif (1) de surveillance de la au moins une machine tournante (100) suivant l'une quelconque des revendications 10 à 12.

Description:
DESCRIPTION

Titre : Procédé, dispositif et programme d'ordinateur de surveillance d'une machine tournante d'un aéronef

L'invention concerne un procédé, un dispositif et un programme d'ordinateur pour la surveillance d'au moins une machine tournante d'un aéronef. Notamment, il peut s'agir d'une machine tournante de propulsion de l'aéronef.

L'invention concerne le domaine de surveillance des machines tournantes par analyse des signaux vibratoires. Les éléments en rotation produisent des vibrations mécaniques pouvant résulter d'un balourd, de l'excentricité de masse, de roulements usés, d'arbres tordus, d'arbres mal alignés, ou autres. Les manifestations du déséquilibre ont pour conséquence de créer des vibrations à caractère cyclique dont la périodicité dépend du composant et du phénomène affecté. Par exemple, un défaut de balourd d'un rotor génère une composante purement sinusoïdale dont la fréquence instantanée est égale à celle de la rotation et d'amplitude proportionnelle au carré de la vitesse de rotation. Un autre exemple est le défaut d' alignement parallèle d'un arbre ayant comme conséquence de générer une vibration sinusoïdale à une fréquence double à celle de l'arbre. Des tels défauts, ainsi que d'autres, peuvent aboutir à la réduction de la durée de la vie des composants de la machine tournante ou un arrêt intempestif du système mécanique en cas d'endommagement (par exemple si l'arbre est cassé). Ces défauts peuvent également affecter d'autres éléments du système. Par exemple, un défaut d'alignement aboutit à une usure prématurée des paliers, joints, arbres et accouplements de la chaîne cinématique.

Pour ces raisons, il est souvent nécessaire de contrôler ces défauts en surveillant le signal vibratoire. L'opérateur est souvent intéressé à évaluer le niveau de balourd ou du mauvais alignement de l'arbre pendant une montée en régime (de la vitesse 0 à la vitesse de fonctionnement) ainsi qu'à la descente. Il est également souhaitable de savoir si le rotor passe par une ou plusieurs résonnances depuis son démarrage jusqu'à sa vitesse de fonctionnement normale. Un filtre suiveur (en anglais « tracking filter ») peut assurer cette fonction en calculant l'amplitude et la phase d'une composante sinusoïdale d'intérêt pendant une montée en régime (en anglais « run-up ») et une descente de régime (en anglais «shut-down »). Un filtre suiveur est équivalent à un filtre passe bande très étroite dont la fréquence centrale suit la fréquence de rotation du composant d'intérêt. Cette dernière est souvent proportionnelle à la fréquence de rotation de l'arbre de moteur.

Par exemple, si l'opérateur désire contrôler l'alignement d'un rotor d'une machine tournante, un filtre suiveur avec une fréquence centrale égale à deux fois la fréquence de rotation de l'arbre peut être appliqué pour extraire la composante d'intérêt. L'amplitude et la phase instantanée de cette composante sont des indicateurs précieux de diagnostic. L'amplitude peut indiquer le niveau et la sévérité du défaut d'alignement (dans certains cas, cela peut être tolérable) ou l'existence d'une résonnance. D'autre part, un changement brusque de la phase peut indiquer le passage par une résonnance.

Les premiers filtres suiveurs étaient implémentés par des circuits électriques analogiques, par exemple les brevets US3307408A et US3938394A. Les problèmes des filtres analogiques sont la complexité des circuits électriques qui nécessite souvent une consommation importante de puissance. Pour cela, les filtres numériques sont aujourd'hui préférés. Un filtre suiveur connu et largement utilisé dans le domaine des machines tournantes est le filtre de Vold-Kalman. La conception de ce filtre repose sur un modèle physique donné au signal vibratoire. En effet, les signaux vibratoires provenant d'organes mécaniques tournants peuvent être modélisés comme une somme de composantes harmoniques noyées dans un bruit additif. Chaque composante peut être représentée comme un signal modulé en phase et en amplitude par un signal basse fréquence. Les modulations lentes d'amplitude et de phase sont physiquement liées au fait que le changement de la vitesse de la machine a pour effet d'exciter différentes régions de la fonction du transfert du système. Cela résulte en un changement de gain et de phase accompagné avec le changement de la vitesse. L'expression analytique de ce modèle connu est donnée par:

La fréquence instantanée f k (n) décrit la variation du contenu fréquentiel au cours du temps. Elle est égale à la dérivée de la phase instantanée selon l'équation suivante : où Fs est la fréquence d'échantillonnage. Il est important de noter que les informations de l'amplitude et de la phase de la k-ième composante sont respectivement l'amplitude et la phase de l'enveloppe a k (n) associée.

La méthode de Vold-Kalman propose une estimation des enveloppes complexes a k (n), associés aux fréquences instantanées f k (n) à partir de deux équations. La première est appelée l'équation de données stipulant que chaque composante à suivre est égale au signal observé à une erreur additive près δ k (n) qui doit être la plus minimale que possible :

La méthode de Vold-Kalman nécessite aussi une information a priori sur la composante à suivre qui sera intégrée dans l'équation structurelle. L'enveloppe complexe à estimer a k (n) est un signal basse fréquence modulant la porteuse exp (jΦ k (n)), ce qui implique une structure régulière et lisse de l'enveloppe. L'idée de base derrière la méthode de Vold-Kalman est d'introduire des contraintes locales favorisant cette propriété. Cela emmène à une deuxième équation, appelée l'équation structurelle, qui tire avantage des différences discrètes pour agir comme un filtre passe bas: où V p est l'opérateur de différence d'ordre p (p désigne également l'ordre du filtre) et ∈ k (n) est un terme d'erreur qui doit être minimisé pour assurer la régularité de l'enveloppe. Généralement, les ordres 1 et 2 sont les plus utilisés et les équations structurelles dont elles dérivent, sont alors données respectivement par : [Math. 6]

(p = 1) : a k (n - 1) - 2 a k (n) + a k (n + 1) = ∈ k (n)

[Math. 7]

(p = 2): a k (n - 2) - 3 a k (n - 1) + 3 a k (n) + a k (n + 1) = ∈ k (n)

La sélectivité du filtre augmente pour des ordres croissants.

Compte tenu de l'équation structurelle et de celle de données, l'estimateur au sens de moindres carrés de l'enveloppe complexe [Math. 8] revient à minimiser la norme [Math. 9] des erreurs selon l'équation suivante :

[Math. 10] où

[Math. 11]

[Math. 12] et r k sont des poids positifs à fixer par l'opérateur pour contrôler le degré de régularité de l'enveloppe. La solution de ce problème est explicite et est donnée sous forme matricielle pour chaque composante k par l'équation suivante : [Math. 13] avec

[Math. 14]

[Math. 15] où A p est la matrice de différences discrètes d'ordre p et C k est la matrice diagonale formée par les termes [Math. 16] exp (jΦ k (1)) exp (jΦ k (N))

Une première possibilité pour l'estimation de plusieurs composantes [Math. 17] a 1 n), ..., a K (n) serait de procéder séquentiellement d'une manière indépendante à l'estimation de chaque composante c'est-à-dire :

[Math. 18]

Cette dernière est connue sous le nom d'estimation mono-composante. Une deuxième possibilité serait d'estimer conjointement ces composantes.

Dans ce deuxième cas, l'équation de données s'écrit :

[Math. 19] et l'estimateur au sens de moindres carrées du vecteur des enveloppes complexes

[Math. 20] est la solution du problème :

[Math. 21] qui est donnée sous forme matricielle par :

[Math. 22] avec

[Math. 23]

= [C 1 , ... , C K ]

, et D est la matrice block-diagonale formée par les sous-matrices D1; ... ; D k .

Cette estimation est alors multi-composante.

Toutefois, aucune de ces deux possibilités utilisant la méthode de Vold- Kalman connue ne donne satisfaction dans le domaine de l'aéronautique.

Il existe un premier problème non résolu par cette méthode de Vold-Kalman connue, dans le cas d'un bruit non-stationnaire et/ou non-gaussien. En effet, la méthode de Vold-Kalman est seulement optimale sous l'hypothèse que le bruit est Gaussien et stationnaire ce qui n'est pas vérifié dans l'aéronautique pour deux raisons :

Certes, le bruit aérodynamique peut être approximé par un bruit Gaussien mais il est loin d'être stationnaire surtout en vol.

Dans le moteur, il existe plusieurs organes tournants qui génèrent des vibrations périodiques et peuvent être modélisés par des modèles sinusoïdaux. Les composantes non prises en compte dans l'équation structurelle vont agir comme un bruit s'ajoutant au bruit aérodynamique et peuvent fausser l'estimation des enveloppes d'intérêt surtout en présence de couplage entre les composantes.

Pour ces raisons, le bruit (qui est ce qui reste dans le signal mesuré si on soustrait les composantes d'intérêt), est composé d'un bruit aléatoire de type Gaussien (stationnaire ou non) et d'un bruit sinusoïdal (composantes sinusoïdales résiduelles interférentes qui ne sont pas prises dans le modèle parmi les K composantes qu'on souhaite estimer). Ce bruit sinusoïdal est loin d'être Gaussien. En plus, en régime variable (quand le moteur tourne à une vitesse/charge non constante), ce bruit sinusoïdal est non stationnaire.

Afin de s'approcher de l'hypothèse du bruit Gaussien et d'atténuer les artéfacts aux instants de couplage, une estimation conjointe (multi-composante) du signal d'intérêt avec toutes les autres harmoniques qui le croisent ou qui sont très proches sur une période de temps est préférable. En effet, lorsque les différentes composantes interférentes sont estimées conjointement, l'équation de données garantit que l'énergie totale du signal sera distribuée entre ces composantes, et avec l'équation structurelle favorisant la régularité des solutions, l'estimation conjointe permet de mieux découpler les composantes dont les fréquences se croisent.

Ainsi, ce premier problème non résolu engendre plusieurs inconvénients :

L'inconvénient majeur de l'estimation multi-composante dans la méthode de Vold-Kalman réside dans un coût de calcul plus important par rapport à l'estimation mono-composante: la dimension de la matrice à inverser est K fois plus grande. En prenant en considération toutes les vibrations, K va être très grand.

Le plus souvent, l'opérateur n'est pas intéressé aux résultats d'estimation de composantes interférentes. Ces composantes sont juste estimées pour diminuer les artéfacts d'estimation des composantes d'intérêt.

Le nombre total de composantes sinusoïdales dans le signal mesuré par le capteur est le plus souvent inconnu.

Pour détecter la présence de composantes interférentes, l'opérateur doit faire une inspection du spectrogramme de signal multi-composant mesuré. De plus, il peut arriver qu'on ne connaisse pas les fréquences instantanées de ces harmoniques interférentes ce qui demande un travail supplémentaire pour estimer ces fréquences avant de procéder à l'estimation par Vold-Kalman, ce qui peut être un travail prohibitif au cas où il existe un grand nombre de composantes interférentes. Pour atténuer ces artéfacts, on peut utiliser une grande valeur de r k . Cela suppose toutefois qu'on connaisse les moments de couplage a priori ce qui n'est pas le cas en pratique. Mettre une grande valeur de r k sur tout le signal, n'est pas une solution fiable, ainsi qu'exposé ci- dessous en référence à la figure 4d, car cela conduit à une sous- estimation de l'amplitude de l'enveloppe. Il existe un deuxième problème non résolu par cette méthode de Vold-Kalman connue, qui concerne la difficulté de régler le paramètre r k . Les méthodes connues nécessitent le réglage manuel des paramètres r k , ou, d'une manière équivalente, de la largeur de bande associée à chaque sinusoïde. Le choix de ces paramètres r k affecte significativement les résultats d'estimation et, donc, un bon choix est requis. La qualité de l'estimation est très liée au choix du paramètre r k . Le paramètre r k peut s'interpréter statistiquement comme le rapport des variances des deux vecteurs aléatoires gaussiens centrés δ k et ∈ k respectivement. Il s'ensuit que plus r k est grand, plus petite sera la variance de ∈ k et donc la solution est de plus en plus influencée par l'équation structurelle résultante d'une enveloppe très lisse. Ceci revient aussi à dire que plus r k est grand, plus la largeur de bande du filtre associé à la k-ième fréquence est étroite et donc plus importante sera la sélectivité. La valeur optimale de r k dépend du niveau signal sur bruit. Plus le bruit est important, plus il faut utiliser une valeur grande de r k pour pouvoir supprimer le bruit. Toutefois, le réglage de ce paramètre r k reste difficile pour plusieurs raisons :

Si bruit est non-stationnaire, ce qui est le cas d'une machine tournante d'aéronef, notamment de propulsion de l'aéronef (le niveau de ce bruit peut changer en vol par exemple), l'utilisation d'une même valeur de r k n'est pas optimale.

L'utilisation d'une valeur très grande de r k peut causer des problèmes numériques. En effet, le produit matriciel donne une matrice semi-définie positive dont la valeur propre maximale augmente linéairement avec r k . Dans une estimation mono-composante, la matrice est une matrice symétrique définie positive (et donc inversible) dont la valeur propre minimale est égale à 1. Dans une estimation multi-composante, C T C est une matrice semi-définie positive, ainsi la valeur propre minimale de C T C + D T D peut être très faible. Ainsi, les valeurs des poids r k ne doivent pas dépasser une certaine limite, au-delà de laquelle la matrice D T D ) dans le cas d'une estimation multi-composante) peut devenir très mal-conditionnée, ce qui résulte d'un problème d'inversion et d'une amplification d'erreurs numériques. La figure 1 montre l'évolution du rapport signal sur bruit SNR de la composante estimée dans un signal mono-composant synthétique bruité en fonction de la valeur choisie du paramètre r = r k . On peut constater qu' au-delà d'une certaine limite du paramètre r k dépendant de l'ordre du filtre, ce rapport signal sur bruit SNR de la composante estimée chute fortement à cause d'un mauvais conditionnement de la matrice à inverser et que par conséquent les performances de la méthode sont dégradées.

Le problème de réglage du paramètre r k est plus difficile dans le cas ou des nombreuses harmoniques (c'est-à-dire que K est grand) sont à suivre. En effet, la valeur de r k ne dépend pas seulement du niveau du bruit mais aussi des statistiques des composantes à estimer. Ainsi considérer la même valeur pour tous les r k n'est pas généralement une solution fiable, vu que les différentes composantes à estimer peuvent avoir différentes plages d'énergie et différents taux de variation. Ce problème peut être rencontré par exemple si on souhaite estimer les composantes liées à l'engrènement de roues dentées rotatives dans la machine tournante et les modulations. En effet, les modulations ont généralement une faible énergie par rapport à la composante liée à l'engrènement.

Plusieurs défauts dans la machine tournante ont pour effet de générer des composantes sinusoïdales à des fréquences connues. Le problème revient à suivre cette fréquence et analyser son amplitude. Dans le cas où le défaut n'existe pas, l'amplitude devrait être faible, mais un mauvais réglage de la méthode Vold-Kalman peut produire une composante fictive dont l'amplitude est égale au bruit de fond lissé ce qui peut engendrer des erreurs dans le diagnostic. Une comparaison de l'amplitude estimée avec le bruit de fond peut résoudre la présence de ces faux positifs. Toutefois, cela suppose qu'on connaisse le bruit. Compte tenu de ces observations, il apparaît que le problème du réglage des paramètres r k n'est pas évident à résoudre et peut affecter la qualité du filtrage, ce qui compromet par conséquent l'efficacité de la méthode de Vold-Kalman.

L'invention vise à obtenir un procédé, un dispositif et un programme d'ordinateur pour la surveillance d'au moins une machine tournante d'un aéronef, qui permettent de résoudre les deux problèmes précités de couplage de composantes et du réglage du paramètre de la méthode Vold-Kalman.

A cet effet, un premier objet de l'invention est un procédé de surveillance d'au moins une machine tournante d'un aéronef, procédé dans lequel on acquiert par un module d'acquisition au moins un signal de mesure de la machine tournante, caractérisé en ce que on estime des fréquences instantanées de composantes sinusoïdales du signal de mesure par un module d'estimation, on exécute par un module de calcul plusieurs itérations successives, dans chacune desquelles : o on met à jour des enveloppes complexes des composantes sinusoïdales, o on met à jour des paramètres d'un modèle d'un bruit du signal de mesure à partir des enveloppes complexes mises à jour, o on teste si le modèle converge de l'itération précédente à l'itération actuelle, pour :

dans la négative, refaire une nouvelle itération,

dans l'affirmative, effectuer un calcul des enveloppes complexes à partir des itérations effectuées.

L'invention offre ainsi la possibilité de résoudre notamment le deuxième problème précité de réglage du paramètre de la méthode Vold-Kalman. L'invention permet l'automatisation du problème d'estimation, ainsi que de généraliser à des interférences non-stationnaires en particulier le couplage de composantes et la non- stationnarité du bruit aérodynamique. L'invention fournit ainsi un procédé et un dispositif d'extraction ou de suivi d'une ou de plusieurs composantes sinusoïdales vibratoires ayant une vitesse instantanée variable. En comparaison avec la méthode Vold-Kalman, l'invention se démarque par les avantages suivants :

Le problème est formulé d'une manière probabiliste. Les sorties de la méthode ne sont plus des estimateurs ponctuels mais des lois décrivant les solutions possibles, ce qui rend possible non seulement d'avoir des estimateurs ponctuels mais aussi une quantification de la confiance sur l'estimation.

L'invention permet l'automatisation de l'estimation en calculant, conjointement avec les enveloppes des estimateurs, des r k . Cela évite le réglage à la main ou heuristique du paramètre par les opérateurs, et maximise l'efficacité du filtrage.

Suivant un mode de réalisation de l'invention, toutes les variables inconnues du problème : les enveloppes complexes, les paramètres du modèle du bruit, les paramètres de régularisation sont estimés à partir de leurs lois a posteriori conditionnelles sachant les mesures. Les paramètres de ces lois a posteriori sont calculés par le module de calcul à chaque itération. Ainsi que cela est connu de l'homme du métier d'une manière générale, une loi a priori est définie par le fait qu'elle est choisie par l'utilisateur de façon à incorporer les caractéristiques désirées de la solution recherchée. Ainsi que cela est connu de l'homme du métier d'une manière générale, une loi a posteriori est définie par le fait qu'elle prend en compte les observations.

Suivant un mode de réalisation de l'invention, ledit calcul des enveloppes complexes dans ladite affirmative est effectué à partir de la moyenne de la loi estimée.

Suivant un mode de réalisation de l'invention, le modèle du bruit du signal de mesure comporte un bruit interfèrent structuré, dont la distribution de probabilité a priori de son spectre a une loi à queue lourde calculée par le module de calcul à chaque itération. Ainsi, suivant un mode de réalisation de l'invention, le modèle du bruit du signal de mesure comporte un bruit interfèrent structuré, dont la distribution de probabilité de sa transformée de Lourier appartient à la classe des distributions à queues plus lourdes que la loi normale. Cela favorise ainsi la parcimonie du bruit structuré sinusoïdal dans le domaine de Fourier.

Suivant un mode de réalisation, la loi à queue lourde du spectre du bruit de mesure est de type Bemoulli-Gaussien. Suivant un mode de réalisation, la loi à queue lourde est donnée par un barycentre entre une première distribution du type Gaussien calculée par le module de calcul à chaque itération et pondérée par une pondération, β ou β j , non nulle et calculée par le module de calcul à chaque itération, et une deuxième distribution de Dirac calculée par le module de calcul à chaque itération et pondérée par une pondération, 1 — β ou i — β j , non nulle et calculée par le module de calcul à chaque itération.

Suivant un mode de réalisation de l'invention, le modèle du bruit du signal de mesure comporte un bruit interfèrent non structuré Gaussien dont les paramètres sont calculés par le module de calcul à chaque itération. Ainsi, suivant un mode de réalisation de l'invention, le modèle du bruit du signal de mesure comporte un bruit interfèrent non structuré dont la distribution de probabilité a priori est donnée par une troisième distribution du type Gaussien calculée par le module de calcul à chaque itération.

Suivant un mode de réalisation de l'invention, le modèle du bruit du signal de mesure résultant d'un modèle de bruit interfèrent structuré et d'un modèle de bruit interfèrent non structuré est la convolution d'une loi Gaussienne et d'une loi à queue lourde dont les paramètres sont calculés par le module de calcul à chaque itération. Suivant un mode de réalisation, la loi à queue lourde du spectre du bruit de mesure est de type Bemoulli-Gaussien. Suivant un mode de réalisation, la convolution est une troisième distribution du type Gaussien calculée par le module de calcul à chaque itération et pondérée par une pondération, 1 - β ou 1 — β j , non nulle calculée par le module de calcul à chaque itération et d'une quatrième distribution du type Gaussien calculée par le module de calcul à chaque itération et pondérée par une pondération, β ou β j , de présence de bruit structuré, non nulle et calculée par le module de calcul à chaque itération, la troisième distribution et la quatrième distribution ayant des variances différentes l'une de l'autre et calculées par le module de calcul à chaque itération. Suivant un mode de réalisation de l'invention, les paramètres du modèle du bruit du signal de mesure comportent au moins la pondération β ou β j et/ou 1 — β ou 1 — β j .

Suivant un mode de réalisation de l'invention, les paramètres des lois du bruit structuré et du bruit non structuré sont considérés comme des variables aléatoires auxquelles sont affectées des lois de probabilités conjuguées avec le modèle du bruit. Suivant un mode de réalisation de l'invention, une loi de probabilité a priori conjuguée est affectée à des hyperparamétres des lois a priori des paramètres du modèle du bruit. Suivant un mode de réalisation de l'invention, le paramètre de précision de la loi Gaussienne du bruit non-structuré est modélisé par une loi Gamma. Suivant un mode de réalisation de l'invention, dans le cas d'une loi de Bemoulli- Gaussien pour le bruit structuré, la probabilité b de présence du bruit structuré est modélisée par une loi a priori uniforme et le paramètre de précision de la loi Gaussienne du bruit structuré est modélisé par une loi a priori Gamma.

Suivant un mode de réalisation de l'invention, une loi de probabilité Bêta a posteriori est suivie par la pondération β ou β j et/ou 1 — β ou 1 — β j et est calculée par le module de calcul à chaque itération.

Suivant un mode de réalisation de l'invention, une loi de probabilité de Gamma a priori est affectée à au moins un T j et une loi de probabilité de Gamma a posteriori est calculée par le module de calcul à chaque itération, où est la variance de la troisième distribution.

Suivant un mode de réalisation de l'invention, une loi de probabilité de Gamma a priori est affectée à au moins un h et une loi de probabilité de Gamma a posteriori est calculée par le module de calcul à chaque itération, où

[Math. 24] est la variance de la quatrième distribution, est la variance de la troisième distribution, est la variance de la première distribution du type Gaussien. Suivant un mode de réalisation de l'invention, les enveloppes complexes à estimer sont modélisées par une loi assurant la régularité des enveloppes et dépend d'un paramètre de régularisation y k qui joue le même rôle que le paramètre r k dans la méthode de Vold-Kalman.

[Math. 25] où

[Math. 26]

C a est une constante indépendante de [Math. 27]

Y 1 , ...,Y K et

[Math. 28] a

Suivant un mode de réalisation de l'invention, le paramètre y k est considéré comme une variable aléatoire et modélisé par une loi conjuguée Gamma. Suivant un mode de réalisation de l'invention, une variable auxiliaire est ajoutée au modèle pour découpler les différentes composantes à estimer. L'intérêt de cette variable auxiliaire n'est pas physique mais c'est seulement pour des raisons d'implémentations .

Suivant un mode de réalisation de l'invention, la variable auxiliaire de découplage est définie par un modèle ayant une distribution de probabilité du type loi normale complexe circulairement symétrique, qui est calculé par le module de calcul à chaque itération.

Suivant un mode de réalisation de l'invention, les premiers paramètres du modèle du bruit du signal de mesure comportent au moins une variance de la ou d'au moins l'une des distributions du type Gaussien. Suivant un mode de réalisation de l'invention, le module de calcul est configuré pour, à chacune des itérations successives effectuées: mettre à jour la loi des enveloppes complexes à partir des paramètres du modèle du bruit, la variable auxiliaire, les paramètres de régularisation et l'observation,

- mettre à jour les paramètres du modèle du bruit à partir des enveloppes mises à jour.

Suivant un mode de réalisation de l'invention, à chaque itération, on met à jour des paramètres de régularisation du modèle par le module de calcul à partir des enveloppes complexes mises à jour.

Suivant un mode de réalisation de l'invention, une loi de probabilité a priori de Gamma est affectée à au moins un des paramètres de régularisation et résulte d'une loi a posteriori de Gamma dont les paramètres sont calculés par le module de calcul à chaque itération.

Suivant un mode de réalisation de l'invention, à chaque itération, on met à jour par le module de calcul une variable auxiliaire de découplage des enveloppes complexes des composantes entre elles dans le modèle, à partir des enveloppes complexes mises à jour et des paramètres du modèle du bruit mis à jour.

Suivant un mode de réalisation de l'invention, la variable auxiliaire de découplage est définie par un modèle ayant une distribution de probabilité du type loi normale complexe circulairement symétrique, qui est calculé par le module de calcul à chaque itération.

Suivant un mode de réalisation de l'invention, le test de convergence comprend de tester à chaque itération actuelle par le module de calcul que la différence entre la valeur absolue de l'enveloppe calculée de l'itération actuelle et la valeur absolue de l'enveloppe calculée de l'itération précédente est inférieure à un seuil prescrit.

Suivant un mode de réalisation de l'invention, on analyse par le module de calcul les enveloppes complexes ayant été calculées, et on déclenche par le module de calcul la communication d'une alarme sur au moins une sortie physique en fonction des enveloppes complexes ayant été analysées. Un deuxième objet de l'invention est un dispositif de surveillance d'au moins une machine tournante d'un aéronef pour la mise en œuvre du procédé tel que décrit ci-dessus, le dispositif comportant

- un module d'acquisition d'au moins un signal de mesure de la machine tournante, caractérisé en ce que le dispositif comporte

- un module d'estimation de fréquences instantanées de composantes sinusoïdales du signal de mesure,

- un module de calcul configuré pour, à chacune de plusieurs itérations successives effectuées: o mettre à jour des enveloppes complexes des composantes sinusoïdales, o mettre à jour des paramètres d'un modèle d'un bruit du signal de mesure à partir des enveloppes complexes mises à jour, o tester si le modèle converge de l'itération précédente à l'itération actuelle, pour :

dans la négative, refaire une nouvelle itération,

dans l'affirmative, effectuer un calcul des enveloppes complexes à partir des itérations effectuées. Suivant un mode de réalisation de l'invention, on met à jour les enveloppes complexes des composantes sinusoïdales d'une façon indépendante.

Suivant un mode de réalisation de l'invention, le module de calcul est configuré pour, à chacune des itérations successives effectuées : mettre à jour des paramètres de régularisation du modèle à partir des enveloppes complexes mises à jour.

Suivant un mode de réalisation de l'invention, le module de calcul est configuré pour, à chacune des itérations successives effectuées :

- mettre à jour une variable auxiliaire de découplage des enveloppes complexes des composantes entre elles dans le modèle, à partir des enveloppes complexes mises à jour et des paramètres du modèle du bruit mis à jour. Suivant un mode de réalisation de l'invention, le module de calcul est configuré pour analyser les enveloppes complexes ayant été calculées, et est apte à déclencher la communication d'une alarme sur au moins une sortie physique du dispositif en fonction des enveloppes complexes ayant été analysées. Un deuxième objet de l'invention est un programme d'ordinateur comportant des instructions de code pour la mise en œuvre du procédé de surveillance tel que décrit ci-dessus, lorsqu'il est exécuté sur un ordinateur.

Un troisième objet de l'invention est un aéronef comportant au moins une machine tournante et un dispositif de surveillance de la au moins une machine tournante tel que décrit ci-dessus.

L'invention sera mieux comprise à la lecture de la description qui va suivre, donnée uniquement à titre d'exemple non limitatif en référence aux figures des dessins annexés, sur lesquelles :

[Fig. 1] représente schématiquement l'évolution du rapport signal sur bruit d'une composante estimée suivant une méthode de l'état de la technique,

[Fig. 2] représente un schéma explicatif du dispositif et du procédé de surveillance suivant un mode de réalisation de l'invention,

[Fig. 3] représente un organigramme explicatif d'une partie du procédé de surveillance suivant un mode de réalisation de l'invention,

[Fig. 4a] montre le spectrogramme initial d'un signal mesuré,

[Fig. 4b] montre un exemple d'enveloppe ayant été calculée par le dispositif et le procédé de surveillance suivant un mode de réalisation de l'invention,

[Fig. 4c] montre le spectrogramme d'un signal résiduel obtenu à partir des figures 4a et 4b par le dispositif et le procédé de surveillance suivant un mode de réalisation de l'invention,

[Fig. 4d] montre le spectrogramme d'un signal résiduel obtenu par un réglage d'une méthode de l'état de la technique,

[Fig. 4e] montre le spectrogramme d'un signal résiduel obtenu par un autre réglage d'une méthode de l'état de la technique, [Fig. 4f] montre le spectrogramme d'un signal résiduel obtenu par encore un autre réglage d'une méthode de l'état de la technique,

[Fig. 5] représente un synoptique modulaire du dispositif de surveillance suivant un mode de réalisation de l'invention.

Le dispositif 1 de surveillance d'une ou plusieurs machine(s) toumante(s) 100 embarquée(s) sur un aéronef est décrit ci-dessous en référence aux figures. La machine tournante 100 peut être une machine tournante de propulsion de l'aéronef, pouvant être une ou plusieurs turbomachine(s), comme par exemple un ou plusieurs turboréacteur(s) ou un ou plusieurs turbopropulseur (s). L'aéronef peut être un avion ou un hélicoptère. Le dispositif 1 de surveillance met en œuvre les étapes du procédé de surveillance de la machine tournante 100, décrit ci-dessous.

Aux figures 2 et 5, le dispositif 1 de surveillance comporte un module 10 d'acquisition d'un ou plusieurs signaux y (t) de mesure d'une ou plusieurs grandeurs sur la machine tournante 100 lors du fonctionnement de cette machine tournante 100 en fonction du temps t. Le signal y(t) de mesure peut être un signal y(t) de mesure de vibrations ou acoustique, par exemple d'un rotor et/ou d'un stator, le signal étant dû à la rotation du rotor par rapport au stator. Le signal y (t) de mesure de vibrations peut être mesuré sur un hélicoptère sur une génératrice de gaz et/ou sur une boîte de vitesse et/ou sur un engrenage et/ou sur un roulement. Le signal y(t) de mesure de vibrations peut être mesuré sur un avion, sur un stator d'une turbomachine. Le module 10 d'acquisition peut comporter un ou plusieurs capteurs, embarqués sur ou à proximité de la machine tournante 100 sur l'aéronef, pour mesurer le ou les signaux y(t), pouvant comprendre par exemple un ou plusieurs accéléromètres de mesure du signal vibratoire y(t) et/ou un ou plusieurs microphones de mesure du signal vibratoire y(t) acoustique, et/ou un ou plusieurs capteur (encodeur ou autre) de vitesse de rotation d'un arbre ou d'un rotor de la machine tournante 100, ou autres.

Suivant un mode de réalisation, le module 10 d'acquisition permet de mesurer et acquérir un signal vibratoire ou acoustique généré par la machine tournante 100 en fonctionnement. Le capteur permet d'obtenir un signal analogique représentatif du signal vibratoire ou acoustique qui est connecté à une chaîne 13 d'acquisition configurée pour fournir le signal numérique de mesure y (t) = y(nTe) avec T e la période d'échantillonnage non nulle, où T = nT e représentent différents instants d'échantillonnages, où n peuvent être par exemple des entiers naturels. Par conséquent, dans ce qui suit la variable n est une variable temporelle. Cette chaîne 13 d'acquisition comprend dans ce cas un conditionneur, un filtre analogique anti- repliement, un échantillonneur bloqueur et un convertisseur analogique - numérique.

Le signal y(t) de mesure comporte K composantes sinusoïdales d'intérêt x k = a k (n) exp (jΦ k (n) ) , dont la somme est égale à x(n) selon l'équation ci- dessous :

[Math. 29]

On note par I = {1, ... , K} l'ensemble des indices des composantes d'intérêt du signal y(t) de mesure.

Φ k. (h) désigne la phase de chaque composante sinusoïdale d'intérêt x k = a k (n)exp(jΦ k (n)). a k (n) désigne l'amplitude complexe (ou enveloppe complexe mentionnée ci- dessous) de chaque composante d'intérêt a k (n)exp (jΦ k (n)) ou la dérivée temporelle d'ordre quelconque de cette composante d'intérêt a k (n)exp (jΦ k (n))

Aux figures 2 et 5, le dispositif 1 de surveillance comporte un module 11 d'estimation de plusieurs fréquences f k (t) instantanées des composantes sinusoïdales x k = a k (n) exp (jΦ k (n)) du signal y(t) de mesure, où k est un entier naturel variable de 1 à K, et où K est un entier naturel prescrit supérieur à 1. Ces fréquences f k (t) instantanées sont les fréquences à surveiller dans la machine tournante 100 (par exemple le premier ordre du générateur de gaz, de la turbine libre, de l'arbre du réducteur, des fréquences d'engrènement de la boite de la vitesse et des engrenages dans le démarreur, ou autres).

Suivant un mode de réalisation de l'invention, une fréquence de rotation f ref (t) d'un arbre de référence de la machine tournante 100 (pouvant être par exemple l'arbre de son rotor) peut avoir été prescrite à l'avance dans le module 11 d'estimation. La fréquence de rotation f ref (t) peut être calculée à partir du signal de vitesse mesuré par un capteur de vitesse de rotation faisant partie du module 10 d'acquisition ou à partir du signal vibratoire ou acoustique y(t) fourni par le module 10 d'acquisition. Le module 11 d'estimation peut calculer la fréquence f k (t) instantanée de chaque composante sinusoïdale x k = a k (n) e x p (jΦ k (n)) comme étant proportionnelle à la fréquence de rotation f ref (t) selon l'équation f k (t) = c kfref (t), où chaque coefficient de proportionnalité c k de la fréquence f k (t) instantanée est réel. La connexion entre les différents composants de la chaîne cinématique de la machine tournante 100 étant rigide, les coefficients de proportionnalité c k peuvent avoir été prescrits à l'avance dans le module 11 d'estimation. A la figure 2, les éléments de la machine 100 qui sont surveillés sont précisés à l'avance et leur cinématique est prescrite dans le module 11 d'estimation, sans qu'il y ait besoin que l'utilisateur donne les fréquences f k (t) en entrées.

Par exemple dans ce qui suit, n est un entier naturel allant de 1 à N, où N est un entier naturel prescrit supérieur à 1.

On pose que le signal y (n) de mesure est égal à la somme des K composantes d'intérêt a k (n)exp(jΦ k (n)) à laquelle est ajoutée le signal (ou bruit) interfèrent l(n) non structuré de type Gaussien (également appelé à large bande) et le signal (ou bruit) interfèrent e(n) structuré (également appelé à bande étroite), selon l'équation suivante : y(n) = x(n) + e(n) + l(n)

Suivant un mode de réalisation de l'invention, le signal (ou bruit) interfèrent e (n) structuré est défini comme étant la partie du bruit interfèrent dont la transformée de Fourier est discrète, c'est-à-dire que ce bruit e(n) peut s'écrire comme la somme de signaux sinusoïdaux. En pratique, ce bruit e(n) est constitué des composantes sinusoïdales résiduelles non prises en compte dans le signal x(n). La transformée de Fourier du bruit interfèrent e(n) structuré est constituées de pics chacun liés à une sinusoïde. On dit alors que le spectre du bruit structuré e(n) est caractérisé par une allure parcimonieuse.

Suivant un mode de réalisation de l'invention, le signal (ou bruit) interfèrent non structuré l(n) est défini comme étant la partie du bruit interfèrent dont la transformée de Fourier n'est pas discrète Le bruit interfèrent non structuré l(n) est constitué des composantes aléatoires non-sinusoïdales dont le spectre n'est pas discret. Ce bruit interfèrent non structuré l(n) peut être par exemple le bruit aérodynamique et est modélisé par une loi Gaussienne Suivant l'invention, e(n) incorpore un signal aléatoire suivant une loi autre qu'une loi normale et doit favoriser un spectre parcimonieux.

Pour rappel et par définition d'une manière générale, une loi normale ou gaussienne complexe f(x) = NC(μ, σ 2 ) en une dimension réelle est définie comme la distribution de moyenne m et d'écart-type s ou variance s 2 ayant la densité de probabilité f(x) selon l'équation suivante :

[Math. 30] Le signal (ou bruit) interfèrent e (n) structuré représente le signal harmonique résiduel, et peut être noté sous la forme suivante :

[Math. 31] où

[Math. 32] est l'ensemble des indices des composantes sinusoïdales du signal interfèrent e(n) structuré, qui sont chacun distinct des indices de l'ensemble /.

Le signal interfèrent total ou bruit total du signal y (t) de mesure est égal à la somme b ( n ) du signal (ou bruit) interfèrent l( n ) non-structuré Gaussien et du signal (ou bruit) interfèrent e(n) structuré selon l'équation suivante : b(n ) = e(n) + l(n)

Le calcul des composantes sinusoïdales d'intérêt a k (n)exp(jΦ k (n)) est du type Bayésien. Ainsi, n, x(n), e(n) et l(n) sont les réalisations temporelles successives de variables aléatoires dont on va définir les lois ci-dessous. De plus, x(n), e(n) et l(n) sont indépendants l'un de l'autre.

Dans la suite, l'exposant R désigne la partie réelle et l'exposant I désigne la partie imaginaire (distinct de l'ensemble I précité).

On définit le vecteur y R comme étant le vecteur d'observation réel, formé des N valeurs réelles y(n) R de y(n) pour l'entier naturel n allant de 1 à N, selon l'équation suivante :

[Math. 33] où l'exposant T désigne la transposée.

On définit le vecteur y 1 comme étant le vecteur d'observation imaginaire, formé des N valeurs imaginaires y(n) 1 de y(n) pour n allant de 1 à N, selon l'équation suivante : [Math. 34]

Les vecteurs x R , x 1 , x k , x k , I R , I I , e R , e I sont également formé des N valeurs réelles ou imaginaires de x(n) ou l(n) et e(n), pour n allant de 1 à N. Modèle du bruit

Ci-dessous, on définit une loi d'observation du bruit b(n) ou modèle du bruit b(n). Ce modèle du bruit b(n) est calculé par un module de calcul 12 (ou unité 12 de surveillance) à la figure 2.

Le module de calcul 12 est numérique (par opposition à des circuits analogiques). Le module de calcul 12 est automatique dans le sens où il effectue automatiquement les étapes décrites ci-dessous. Ainsi que représenté à la figure 5, le module de calcul 12 peut être ou comprendre un ou plusieurs calculateurs, et/ou un ou plusieurs ordinateurs, et/ou une ou plusieurs machines de calcul, et/ou un ou plusieurs processeurs et/ou un ou plusieurs microprocesseurs. Le module de calcul 12 et/ou le module 10 d'acquisition peut comporter un programme d'ordinateur enregistré dans leur mémoire pour exécuter d'une manière automatique les étapes décrites ci-dessous. En outre, le dispositif 1 de surveillance peut comporter un ou plusieurs sorties physiques SP de présentation d'informations pour présenter à un utilisateur les enveloppes complexes a k ( n ) ayant été calculées par le module 12, cette ou ces sorties physiques SP pouvant comprendre un écran d'affichage pour afficher ces enveloppes complexes a k (n ) ayant été calculées (ou leur dérivée d'ordre quelconque (amplitude, phase, amplitude en fonction de vitesse, la ou les composantes estimées, etc.)). La mémoire est prévue pour stocker les sorties a k (n ) ayant été calculées par le module 12 de calcul.

Modèle du bruit e(n) structuré

Le modèle du bruit interfèrent e(n ) structuré peut être réalisé par l'un des modes de réalisation décrits ci-dessous.

Le modèle du bruit e (n) est défini comme étant un signal interfèrent structuré de la manière suivante. Suivant un mode de réalisation de l'invention, le spectre du bruit interfèrent e(n) structuré est défini par un modèle présentant un spectre parcimonieux, c'est-à-dire une distribution concentrée à zéro pour favorisant les valeurs nulles avec une grande probabilité et ayant des queues plus lourdes que la loi normale pour autoriser la présence des valeurs grandes. Par exemple, une bonne distribution de probabilité du spectre du bruit interfèrent e(n) structuré peut être définie par un modèle étant un mélange d'un nombre fïni/infïnie de distributions différentes pondérées.

Suivant un mode de réalisation de l'invention, le spectre du bruit interfèrent à bande étroite est défini par une somme pondérée d'une distribution de probabilité du type Gaussien et d'une autre distribution de probabilité du type non Gaussien.

Parmi les distributions de mélanges infinies, il y a par exemple la loi de Laplace et la loi de Student couramment utilisées pour modéliser les signaux parcimonieux et qui peuvent s'exprimer comme des mélanges infinis de Gaussiennes avec des paramètres de mélanges suivant des lois gamma et des lois inverse gamma respectivement. Parmi les distributions de mélanges finies, il y a par exemple le mélange de deux distributions : d'une distribution de Bernoulli et d'une distribution Gaussienne connue sous le nom de Bemoulli-Gaussien. Dans le cas du modèle de Bemouilli, le spectre du bruit interfèrent e(n) à structuré est défini par un modèle ayant une distribution de probabilité donnée par le barycentre entre une distribution du type Gaussien pondérée par une pondération (probabilité) b de présence du bruit interfèrent e(n) et une distribution de Dirac d 0 pondérée par une pondération 1 — β, par exemple selon l'équation suivante :

[Math. 35] où

F N désigne l'opérateur de transformée de Fourier discrète de taille N, cet opérateur étant divisé par [Math. 36] , d 0 est la distribution de Dirac positionnée à l'abscisse zéro,

0 ≤ β ≤ 1, est la probablité de présence du bruit structuré, désigne une distribution normale complexe de moyenne nulle, dont la variance est égale à Ainsi, la partie réelle de cette distribution normale complexe est une distribution gaussienne de moyenne nulle et de variance La partie imaginaire de cette distribution normale complexe est une distribution gaussienne de moyenne nulle et de variance Suivant un mode de réalisation, b est différent de zéro.

Ci-dessus, [F N e] n désigne l'opérateur de transformée de Fourier discrète de taille N, qui est divisé par [Math. 37] et est appliqué à e(n). Ainsi, on a

[Math. 38]

Et

[Math. 39] où Id désigné la matrice identité.

A la différence de la méthode de Vold-Kalman connue, où le bruit est supposé gaussien centré, ce mode de réalisation permet de prendre en compte facilement dans le modèle le bruit structuré afin d'améliorer l'efficacité de l'estimation au moment de couplage sans avoir besoin d'estimer les composantes interférentes. En régime variable, le spectre de e(n) peut avoir des pics plus élargis. On définit des transformées S j de Fourier, limitée à un segment temporel de taille N j pour e (n) comme suit :

[Math. 40] avec

[Math. 41]

S J = F Nj P j où P j est une matrice de sélection de taille [Math. 42]

N j x N N j est un entier naturel prescrit inférieur à N ( N j est la longueur de la fenêtre),

[Math. 43] avec

[Math. 44]

DFT Nj la transformée de Fourier discrète de taille N j , n est un entier naturel allant de 1 à N j ,

J est un entier naturel prescrit (nombre de fenêtres temporelles j), j est un entier naturel allant de 1 à J, désigne une distribution normale complexe de moyenne nulle, dont la variance est égale à β j est la pondération de présence du bruit structuré, affectée à

0 ≤ β j ≤ 1,

1 — β j est la pondération affectée à la distribution de Dirac positionnée à l'abscisse zéro. Les transformées S) de Fourier sont donc de plus courte durée que FN. Si en plus les matrices P j ont des supports disjoints, on obtient

[Math. 45] et donc

[Math. 46]

Les matrices Pj ont des supports disjoints si la somme de tous les N j est égale àN :

[Math. 47]

La matrice de sélection P j est telle qu'elle permet de sélectionner N j valeurs appartenant à la fenêtre j d'un vecteur de taille N.

P j est une matrice contenant N j lignes et N colonnes, elle est définie par :

Pourj=l :

[Math. 48] Pj(io, i 1 ) = 0.

Pour j>l

[Math. 49]

Pj(io, i 1 ) = 0.

A la différence de la méthode Vold-Kalman connue où le bruit est supposé stationnaire, ce qui est rarement le cas du bruit aérodynamique, ce mode de réalisation permet de prendre en compte la non-stationnarité du bruit, et ce en analysant le signal à travers les transformées de Fourier S j à court terme mentionnées ci-dessus. L'axe des temps est divisé en J intervalles ou fenêtres. Les paramètres du modèle peuvent changer d'une fenêtre à l'autre. Vu que les paramètres de ces lois sont inconnus (on ne connaît pas le niveau du bruit aérodynamique et on ne connaît pas les instants de couplage entre les composantes ni l'amplitude des composantes interférentes), ils vont être adaptés aussi automatiquement. Des lois a priori sont choisies pour modéliser l'information a priori sur ces paramètres. En général, si on ne dispose pas d'informations a priori, on peut choisir des lois dites non-informatives.

Modèle du bruit l(n) non structuré Le modèle du bruit interfèrent /(n) non structuré peut être réalisé par l'un des modes de réalisation décrits ci-dessous.

Suivant un mode de réalisation de l'invention, le spectre du bruit interfèrent l(n) non structuré est défini par un modèle ayant une distribution de probabilité du type Gaussien.

Suivant un mode de réalisation de l'invention, le spectre du bruit interfèrent l(n) non structuré est défini par un modèle ayant une distribution de probabilité NC(0, 2 σ -1 ) normale complexe de moyenne nulle et de variance 2 σ -1 , par exemple selon l'équation suivante :

[Math. 50] II peut s'avérer que le bruit de fond soit non-stationnaire. Dans ce cas, on peut procéder par segments temporels j comme on l'a fait pour le signal interfèrent e(n). On définit des transformées Sj de Fourier, limitées à un segment temporel de taille Nj pour l ( n ) comme suit : [Math. 51] ont les définitions mentionnées ci-dessus, désigne une distribution normale complexe de moyenne nulle, dont la variance est égale à

Dans chaque fenêtre temporelle j, le bruit interfèrent l(n) non structuré est Gaussien et a une variance On écrit alors

[Math. 52] Modèle du bruit total b (n)

Le modèle du bruit interfèrent total b (n) peut être réalisé par l'un des modes de réalisation décrits ci-dessous. Suivant un mode de réalisation de l'invention, le spectre du bruit interfèrent total b (n) résulte des modèles du bruit interfèrent structuré et du modèle du bruit non structuré. En effet, puisque e(n) et l(n) sont indépendants, la densité de probabilité de la somme b(n)= e(n)+ l(n) est la convolution de la densité de probabilité du bruit interfèrent structuré e(n) avec celle du bruit non structuré l(n) . Suivant un mode de réalisation de l'invention, par exemple, dans le cas d'une loi de Bemoulli-Gaussien pour le bruit structuré, le spectre du bruit interfèrent total b (n) est défini par un modèle ayant un mélange de deux distributions du type Gaussien ayant chacune une pondération associée et de variances différentes.

Suivant un mode de réalisation de l'invention, le spectre du bruit interfèrent total b (n) est défini par un modèle ayant une distribution donnée par le barycentre entre la distribution du type Gaussien pondérée par la pondération 1 — β j et une distribution du type Gaussien pondérée par la pondération β j , par exemple selon l'équation suivante : [Math. 53] ont les définitions mentionnées ci-dessus, désigne une distribution normale complexe de moyenne nulle, dont la variance est égale à

Dans ce qui suit,

[Math. 54] Ainsi, dans ce mode réalisation, le spectre du bruit interfèrent total b ( n ) est défini par un modèle ayant une distribution donnée par le barycentre entre la distribution du type Gaussien pondérée par la pondération 1 — β j et une distribution du type Gaussien pondérée par la pondération β j , par exemple selon l'équation suivante :

[Math. 55] ont les définitions mentionnées ci-dessus, désigne une distribution normale complexe de moyenne nulle, dont la variance est égale à 2n j -1 , avec n j < T j . τ j et n j sont également appelés paramètres de précision du bruit.

La loi d'observation y R , y 1 du signal y(t) de mesure est donnée alors par la Gaussienne N selon le modèle hiérarchique suivant :

[Math. 56] où la matrice C est la matrice diagonale formée par les termes exp(jΦ k (n)) pour n allant de 1 à N et k allant de 1 à K, a est le vecteur des enveloppes, diag(d j ) est une matrice dont la diagonale est formée des coefficients d j et dont les autres coefficients sont nuis, d j est une variable qui suit une distribution donnée par le barycentre entre une distribution de Dirac

[Math. 57] positionnée en l'abscisse X j 1 et pondérée par la pondération 1 — β,· et une distribution de Dirac

[Math. 58] positionnée en l'abscisse n j -1 et pondérée par la pondération β j .

Suivant un mode de réalisation de l'invention, la loi a posteriori de d j est une distribution donnée par le barycentre entre une distribution de Dirac

[Math. 59] positionnée en l'abscisse X j et pondérée par la pondération 1-p j (n) et une distribution de Dirac [Math. 60] positionnée en l'abscisse et pondérée par la pondération p j (n), par exemple selon l'équation suivante :

[Math. 61]

[Math. 62] [Math. 63] Dans ce qui précède,

[Math. 64]

[Math. 65]

[Math. 66]

De plus, diag(d j ) ou diag(d j (n)) est une matrice de covariance diagonale du bruit.

Itération de calcul des enveloppes complexes a k (n) des composantes sinusoïdales

Suivant l'invention, le dispositif 1 de surveillance comporte un module 12 de calcul configuré pour, à chaque itération respective effectuée : - au cours d'une première étape Cl de l'itération respective effectuée mettre à jour les lois (calculer ou donner une valeur aux paramètres (moyenne et variance) de la loi) des enveloppes complexes a k (n) des composantes sinusoïdales, au cours d'une deuxième étape C21 de l'itération respective effectuée (postérieure à la première étape Cl de l'itération respective effectuée) mettre à jour les paramètres du modèle du bruit b (n) et/ou e(n) et/ou l(n)) du signal (y(t)) de mesure à partir des enveloppes complexes a k (n) mises à jour ou sachant les enveloppes complexes a k (n) mises à jour,

- tester au cours d'une troisième étape C4 de l'itération respective effectuée (postérieure à la deuxième étape C21 de l'itération respective effectuée) si le modèle converge de l'itération précédente (t — 1) à l'itération actuelle (t), pour : o dans la négative (NON à la figure 2), refaire une nouvelle itération des étapes Cl, C21, C4 décrite ci-dessus, o et dans l'affirmative (OUI à la figure 2), effectuer lors de la quatrième étape D postérieure à la troisième étape C4, un calcul des enveloppes complexes a k ( n ) à partir des itérations effectuées, pour obtenir une extraction automatique des composantes sinusoïdales a k (n ).

La mise à jour des paramètres du bruit peut être faite d'une manière indépendante sur les j allant de 1 à /. A chaque itération t, étant donné les lois a posteriori des variables à l'itération précédente t — 1 et l'observation y, les lois a posteriori des différents variables sont mises à jour par le module de calcul 12, par exemple selon le digramme montré à la figure 3.

La troisième étape C4 de test permet de vérifier si les lois des différents variables du problème à estimer à une itération donnée sont bien optimales. Suivant un mode de réalisation de l'invention, le test de la troisième étape C4 comprend le calcul automatique d'un ou plusieurs indicateurs par le module 12 de calcul à partir des enveloppes a k (n ) ayant été calculées (par exemple leur valeur absolue).

Suivant un mode de réalisation de l'invention, le test de la troisième étape C4 peut être effectué par le module 12 de calcul qui compare la valeur absolue de la différence entre la moyenne de la loi de l'enveloppe a k (n ) calculée de l'itération actuelle et la moyenne de la loi de l'enveloppe a k (n ) calculée de l'itération précédente à un seuil prescrit à l'avance pour déterminer si cette différence est inférieure à ce seuil (affirmative de convergence ci-dessus (OUI à la figure 2) lorsque cette différence est inférieure au seuil ; négative de convergence ci-dessus (NON à la figure 2) lorsque cette différence n'est pas inférieure au seuil). Par exemple, ce seuil peut être prescrit entre 10 -3 et 10 -6 fois la norme de la moyenne de la loi de l'enveloppe a k (n ) calculée de l'itération précédente ou actuelle. Le test de la troisième étape C4 peut également être effectué par le module 12 de calcul sur l'une et/ou plusieurs des variables décrites du modèle du bruit, comme par exemple une moyenne, une variance ou autres, par exemple de l'enveloppe a k (n) calculée, pour tester si cette variable converge de l'itération précédente à l'itération actuelle. Tant que les statistiques ne sont stabilisées (NON à la figure 2), l'algorithme continue à mettre à jour les variables d'une façon itérative.

Suivant un mode de réalisation de l'invention, lors de la quatrième étape D, le calcul des enveloppes complexes a k (n) peut être effectué par le module 12 de calcul à partir de la moyenne de la loi estimée. Des statistiques supplémentaires peuvent être aussi données comme l'intervalle de confiance.

Suivant un mode de réalisation de l'invention, un nombre maximal d'itérations est prescrit à l'avance dans le module 12 de calcul.

Paramètres du modèle du bruit

Suivant un mode de réalisation de l'invention, ces paramètres peuvent être un ou plusieurs parmi les paramètres paramétrant la loi Gaussienne du bruit non- structuré (moyenne, précision) et ceux paramétrant la loi à queue lourde du bruit structuré (moyenne, précision, paramètre de forme etc) et cela pour une ou plusieurs fenêtres parmi les / fenêtres.

Suivant un mode de réalisation de l'invention, dans le cas d'un modèle de Bemoulli-Gaussien ces paramètres peuvent être un ou plusieurs ou tous parmi les X j pour j allant de 1 à J, et/ou un ou plusieurs ou tous parmi les β j pour j allant de 1 à J, et/ou un ou plusieurs ou tous parmi les n j pour j allant de 1 à J, et/ou un ou plusieurs ou tous parmi les ( j pour j allant de 1 à/. Ce mode de réalisation des paramètres peut être réalisé par l'un ou plusieurs des autres modes de réalisation décrits ci-dessous.

Suivant un mode de réalisation de l'invention, une loi de probabilité a priori est prédéfinie pour ces paramètres. Suivant un mode de réalisation de l'invention, des lois conjuguées sont choisies. Cela facilite les calculs par le module de calcul 12.

Suivant un mode de réalisation de l'invention, dans le cas d'un modèle Bemoulli-Gaussien pour le bruit structuré, une loi de probabilité uniforme sur l'intervalle allant de 0 à 1 est affectée à un ou plusieurs ou tous les β j , c'est-à-dire [Math. 67]

Suivant un mode de réalisation de l'invention, une loi de probabilité du type loi Bêta est affectée à un ou plusieurs ou tous les β j .

Pour rappel, par convention et d'une manière générale, une loi Bêta (notée Beta(ff, H)) ayant un troisième réglage a (réel strictement positif) et un quatrième réglage H (réel strictement positif) est définie par la fonction de densité de probabilité f(x; a, H) suivant l'équation ci-dessous :

[Math. 68]

Suivant un mode de réalisation de l'invention, la loi a posteriori de un ou plusieurs ou tous les β j résultant d'une loi a priori uniforme et sachant les observations y ((n) est une loi Beta(2n 1 + 1, 2 n 2 + 1) selon l'équation suivante :

[Math. 69]

[Math. 70] avec

[Math. 71]

[Math. 72] où la fonction < > utilisée ci-dessus signifie d'une manière générale < P > = 1 si P est vrai, ou sinon < P > = 0.

Cette loi de probabilité Beta(2n 1 + 1, 2 n 2 + 1) est la loi a posteriori qui résulte de la loi a priori (loi uniforme) et de la loi des observations (loi de y). Suivant un mode de réalisation de l'invention, une loi de probabilité de

Gamma est affectée à un ou plusieurs ou tous les t,·, par exemple [Math. 73] ayant [Math. 74] comme premier réglage et [Math. 75] comme deuxième réglage.

Pour rappel, par convention et d'une manière générale, une loi de Gamma (notée Gamma (H, a)) ayant un premier réglage H (réel strictement positif) et un deuxième réglage a (réel strictement positif) est définie par la fonction de densité de probabilité f(x; H, α) suivant l'équation ci-dessous :

[Math. 76] où est la fonction gamma d'Euler selon l'équation ci-dessous

[Math. 77]

Suivant un mode de réalisation de l'invention, une loi de probabilité de Gamma est affectée à un ou plusieurs ou tous les n j , par exemple [Math. 78]

Ayant [Math. 79] comme premier réglage et [Math. 80] comme deuxième réglage.

Les réglages [Math. 81]

[Math. 82]

[Math. 83] [Math. 84] des lois de Gamma sont chacun strictement positifs. Suivant un mode de réalisation de l'invention, [Math. 85]

[Math. 86]

[Math. 87]

[Math. 88] des lois de Gamma sont chacun prescrits à l'avance dans le module de calcul 12. Ainsi,

[Math. 89]

[Math. 90]

[Math. 91]

[Math. 92] des lois de Gamma sont chacun un hyper-paramètre quantifiant les informations a priori sur les premiers paramètres à estimer. Suivant un mode de réalisation de l'invention, chaque couple (a, h) de premier réglage et de deuxième réglage, tel que par exemple

[Math. 93]

, est prescrit à l'avance dans le module de calcul 12 selon une information a priori sur la distribution du premier paramètre θ considéré avec

[Math. 94]

. Ces réglages peuvent être prescrits à l'avance par l'utilisateur. Si l'utilisateur ne prescrit pas ces réglages, ces réglages sont fixés à des valeurs très faibles proches ou égales à zéro (voir l'exemple ci-dessous) pour avoir des lois non-informatives. Cette information peut être une moyenne de la loi considérée pour Q et/ou une variance de la loi considérée pour Q. Suivant un mode de réalisation de l'invention, si on connaît approximativement la valeur de Q notée par θ e > 0 avec une incertitude égale à eθ e où e > 0, alors

[Math. 95] et

[Math. 96] sont prescrits à l'avance. Si on n'a aucune information a priori sur le paramètre Q, alors a = b = 0 sont prescrits à l'avance (la loi est alors non- informative). La présence d'informations a priori sur les paramètres permet le plus souvent d'accélérer la convergence.

Suivant un mode de réalisation de l'invention, la loi a posteriori de un ou plusieurs ou tous les n j résultant de sa loi a priori Gamma de paramètres (

[Math. 97] [Math. 98]

) et le signal mesuré y (n) est une loi [Math. 99] ayant

[Math. 100] comme premier réglage et [Math. 101] comme deuxième réglage selon l'équation suivante :

[Math. 102]

[Math. 103]

Suivant un mode de réalisation de l'invention, la loi a posteriori de un ou plusieurs ou tous les τ j résultant de sa loi a priori Gamma de paramètres (

[Math. 104]

[Math. 105]

) et le signal mesuré y (n) est une loi Gamma [Math. 106] ayant

[Math. 107] comme premier réglage et [Math. 108] comme deuxième réglage selon l'équation suivante :

[Math. 109]

[Math. 110]

Dans ce qui précède,

[Math. 111]

[Math. 112]

[Math. 113] [Math. 114]

Paramètres y k de régularisation

Suivant un mode de réalisation de l'invention, le module 12 de calcul est configuré pour, au cours d'une cinquième étape C22 de l'itération respective effectuée (postérieure à la deuxième étape C21 de l'itération respective effectuée et antérieure à la troisième étape C4) mettre à jour des deuxièmes paramètres y k de régularisation du modèle à partir des enveloppes complexes a k (n) mises à jour ou sachant les enveloppes complexes a k (n) mises à jour. Ces deuxièmes paramètres y k de régularisation sont des contraintes supplémentaires imposées afin de favoriser la régularité des enveloppes complexes a k (n). Ces deuxièmes paramètres y k de régularisation caractérisent la suppression de singularités dans les composantes a k (n). Ce mode de réalisation des deuxièmes paramètres y k de régularisation peut être réalisé par l'un des autres modes de réalisation décrits ci-dessous.

La mise à jour des deuxièmes paramètres peut être faite d'une manière indépendante sur les k allant de 1 à K.

Suivant un mode de réalisation de l'invention, la densité p de probabilité (a priori) des enveloppes complexes a k (n) en présence de ces deuxièmes paramètres y k de régularisation est la suivante, pour k allant de 1 à K:

[Math. 115] où C a est une constante indépendante des y k et des enveloppes complexes a k (n). Suivant un mode de réalisation de l'invention, une loi de probabilité de Gamma est affectée à un ou plusieurs ou tous les deuxièmes paramètres y k de régularisation, par exemple [Math. 116] ayant

[Math. 117] comme premier réglage et [Math. 118] comme deuxième réglage. Suivant un mode de réalisation de l'invention, le couple (

[Math. 119]

[Math. 120]

) est prescrit à l'avance dans le module de calcul 12 selon une information a priori sur la distribution du premier paramètre y k Cette information peut être une moyenne de la loi considérée pour y k et/ou une variance de la loi considérée pour y k . Ces réglages peuvent être prescrits à l'avance par l'utilisateur. Si l'utilisateur ne prescrit pas ces réglages, ces réglages sont fixés à des valeurs très faibles proches ou égales à zéro (voir l'exemple ci-dessous) pour avoir des lois non-informatives. Suivant un mode de réalisation de l'invention, si on connaît approximativement la valeur de y k notée par θ e > 0 avec une incertitude égale à e θ e où e > 0, alors [Math. 121] et

[Math. 122] sont prescrits à l'avance. Si on n'a aucune information a priori sur le paramètre θ, alors [Math. 123] sont prescrits à l'avance (la loi est alors non-informative). La présence d'informations a priori sur les paramètres permet le plus souvent d'accélérer la convergence.

Suivant un mode de réalisation de l'invention, une loi de probabilité a posteriori de y k (résultant de la loi a priori Gamma de paramètres de réglage ( [Math. 124]

[Math. 125]

), et de la loi a priori des enveloppes) est aussi une loi Gamma [Math. 126] ayant

[Math. 127] comme premier réglage et

[Math. 128] comme deuxième réglage, par exemple selon l'équation suivante :

[Math. 129]

Dans ce qui précède, [Math. 131]

Suivant un mode de réalisation de l'invention, une loi a posteriori des composantes sinusoïdales d'intérêt est calculée lors de la première étape Cl par le module de calcul 12 à partir de la loi d'observation du bruit b(n ) et de la loi d'observation du signal y(t) de mesure. Suivant un mode de réalisation de l'invention, la loi a posteriori de la partie réelle a R des composantes sinusoïdales d'intérêt est définie par un modèle ayant une distribution de probabilité de type Gaussien de moyenne (partie réelle de m x ) et de variance selon l'équation suivante :

[Math. 133]

Suivant un mode de réalisation de l'invention, la loi a posteriori de la partie imaginaire a 1 des composantes sinusoïdales d'intérêt est définie lors de la première étape Cl par un modèle ayant une distribution de probabilité de type Gaussien de moyenne (partie imaginaire de m x ) et de variance par exemple selon l'équation suivante :

[Math. 134]

La partie réelle a R des composantes sinusoïdales d'intérêt et la partie imaginaire a 1 des composantes sinusoïdales d'intérêt peuvent être mises à jour d'une manière indépendante lors de la première étape Cl. Dans le cas d'un nombre très grand de composantes à estimer, cette première étape Cl peut être faite d'une façon parallèle grâce à une architecture multi-processeurs du module 12.

Dans ce qui précède,

[Math. 135]

[Math. 136] diag(d j ) est une matrice dont la diagonale est formée des coefficients d j et dont les autres coefficients sont nuis, diag(y k ) est une matrice dont la diagonale est formée des coefficients y k et dont les autres coefficients sont nuis,

A p est la matrice des différences discrètes V p d'ordre p de diag(y k ),

A est la matrice bloc diagonale formée par K matrices A p et est la matrice contenant K blocs de-4 p, la matrice C est la matrice diagonale formée par les termes exp(jΦ k (n)) pour n allant de 1 à N.

Variable auxiliaire v de découplage des enveloppes

Suivant un mode de réalisation de l'invention, on calcule lors de la sixième étape C3 (postérieure à la deuxième étape C21 de l'itération respective effectuée et/ou à la cinquième étape C22 et antérieure à la troisième étape C4) par le module 12 de calcul une variable auxiliaire v de découplage des enveloppes complexes a k (n) des composantes entre elles dans le modèle, à partir des enveloppes complexes a k (n) mises à jour, des premiers paramètres (pouvant être un ou plusieurs ou tous parmi les τ j pour j allant de 1 à /, et/ou un ou plusieurs ou tous parmi les β j pour j allant de 1 à J, et/ou un ou plusieurs ou tous parmi les g j pour j allant de 1 à J, et/ou un ou plusieurs ou tous parmi les pour j allant de 1 à J) et des deuxièmes paramètres y k de régularisation. Un problème connu des algorithmes d'optimisation est l'inversion d'une matrice pleine de grande dimension, le plus souvent mal conditionnée. Ce mode de réalisation permet de résoudre ce problème en découplant les composantes par ajout de la variable auxiliaire v. La variable auxiliaire v prend en compte les informations concernant les corrélations entre les composantes a k (n). Grâce à la variable auxiliaire v, les composantes sont estimées indépendamment. Le couplage intervient implicitement à travers la variable auxiliaire.

Ce mode de réalisation permet de résoudre les problèmes numériques en grande dimension, liés au premier problème précité. Ce mode de réalisation est d'autant plus adapté qu'un grand nombre K de composantes sont souhaitées. En effet, la méthode Vold-Kalman connue souffre du problème d'inversion de matrice de grande dimension surtout quand la matrice en question est très mal conditionnée ce qui est le cas pour des valeurs de r k très grandes. Cette inversion pouvait se faire dans la méthode Vold-Kalman connue grâce à un algorithme itératif de type gradient conjugué. Ce mode de réalisation n'est pas contrarié par ce problème numérique, vu que les différentes composantes sont découplées dans l'algorithme et donc peuvent être estimées d'une manière disjointe. Le couplage entre les différentes composantes est pris en compte implicitement par la variable auxiliaire v ajoutée dans la mise à jour de chaque composante.

Ce mode de réalisation de l'invention sur la variable auxiliaire v de découplage peut être réalisé par l'un des autres modes de réalisation décrits ci- dessous.

Suivant un mode de réalisation de l'invention, la loi de la variable auxiliaire v de découplage est définie par un modèle ayant une distribution de probabilité du type loi normale complexe circulairement symétrique ayant la moyenne m v et la variance On a ainsi pour v une loi conditionnelle

[Math. 137]

La loi normale complexe circulairement symétrique a la moyenne

[Math. 138] et de matrice de covariance

[Math. 139] et est définie par [Math. 140]

Dans ce qui précède,

[Math. 141]

[Math. 142] avec m étant un réel prescrit non nul vérifiant [Math. 143]

Cela permet de pallier aux problèmes d'inversion de matrices de grandes dimensions en particulier pour des grandes valeurs de N et K et permet une augmentation des données permettant le découplage entre les différentes composantes

Calcul des composantes sinusoïdales ar d'intérêt lors de la première étape Cl

Suivant un mode de réalisation de l'invention, la loi a posteriori de la partie réelle a k R des composantes sinusoïdales d'intérêt est définie par un modèle ayant une distribution de probabilité du type loi normale ayant la moyenne m k R et la variance par exemple selon l'équation suivante : [Math. 144]

Suivant un mode de réalisation de l'invention, la loi a posteriori de la partie imaginaire des composantes sinusoïdales d'intérêt est définie par un modèle ayant une distribution de probabilité du type loi normale ayant la moyenne et la variance par exemple selon l'équation suivante :

[Math. 145]

Les définitions ci-dessus de sont des lois a posteriori des composantes.

Dans ce qui précède,

[Math. 146]

[Math. 147]

Ainsi, a k R et a k 1 peuvent être estimées indépendamment. Le couplage est pris en compte seulement à travers la variable auxiliaire v intervenant dans la moyenne m k . Il est intéressant de noter que même si on a besoin de faire l'inversion ici, cete opération est moins coûteuse. D'un côté, la dimension de la matrice est K fois plus petite que la matrice D'un autre côté, si on opte pour des conditions de bords circulants (ou bords à zéro en faisant un padding dans les bords par zéro) dans la matrice de différences A p , la matrice est alors circulante et donc diagonalisable dans le domaine de Fourier. Suivant un mode de réalisation de l'invention,

[Math. 148] où Qk est la matrice diagonale contenant les coefficients de Fourier du filtre associé à la matrice de différence discrète A p qui est égale à la première colonne de A p . La mise à jour de chaque m k peut aussi être faite en passant dans le domaine de Fourier. De plus, on n'a pas besoin de stocker ni son inverse, mais simplement les valeurs diagonales de la matrice

[Math. 149]

Suivant un mode de réalisation de l'invention, le module 12 de calcul est configuré pour effectuer une septième étape E de décision postérieure à la quatrième étape D. Suivant un mode de réalisation de l'invention, la décision lors de septième étape E est prise automatiquement par le module 12 de calcul à partir d'analyses effectuées automatiquement par le module 12 de calcul sur les composantes extraites a k (n ) ayant été obtenues à la quatrième étape D, et/ou à partir d'une caractérisation vibratoire en banc d'essais (par exemple identification des modes propres) effectuée par le module 12 de calcul sur les composantes extraites a k (n ) ayant été obtenues à la quatrième étape D de la machine tournante 100. Suivant un mode de réalisation de l'invention, la septième étape E est mise en œuvre par le module 12 de calcul sur une décision de surveillance et/ou sur le déclenchement d'une alarme AL ou alerte AL de défauts de la machine tournante 100 (par exemple : balourd, engrenage, désalignement, excentricité de masse, roulements usés, arbres tordus, arbres mal alignés, défaut d'alignement parallèle d'un arbre d'un rotor de la machine tournante, arbre cassé, usure prématurée des paliers, joints, arbres et accouplements de la chaîne cinématique, résonance). Suivant un mode de réalisation de l'invention, le dispositif 1 de surveillance est embarqué sur l'aéronef. Suivant un mode de réalisation de l'invention, le module 12 de calcul commande automatiquement la communication de l'alerte, de l'alarme AL ou de la décision qu'il a prise lors de septième étape E sur la sortie physique SP, par exemple sur l'écran d'affichage, et/ou à une station terrestre fixe SP. La mise à jour d'une loi dans le dispositif proposé peut se faire de deux manières.

D'une part, elle peut être approximée par des échantillons. On parle alors des méthodes Monte Carlo Markov Chain (MCMC). Ainsi, mettre à jour une loi revient à tirer un échantillon de cette loi. Notons que toutes les lois sont simples à échantillonner. Pour les paramètres de la loi d'observation et du paramètre de régularisation, cela est évident. Pour la loi a posteriori des enveloppes, la génération de variables aléatoires peut se faire en passant dans le domaine de Fourier. Et vu la forme particulaire de å v , la génération des variables aléatoires auxiliaires peut aussi se faire directement en utilisant le fait que Après un nombre suffisant d'itérations, les variables aléatoires générées par l'algorithme vont suivre les lois a posteriori recherchées ; ainsi les statistiques (moyenne, variance etc) peuvent être approximés en utilisant des estimateurs empiriques avec ces échantillons.

D'autre part, la loi a posteriori peut être approximé par une autre plus simple, par exemple par une approche d'approximation Bayésienne -Variationnelle (ABV). Ainsi les statistiques des lois d'intérêt peuvent être approximées par celles des lois approximantes après un nombre suffisant d'itérations.

On donne ci-dessous un exemple d'application de l'invention sur un signal vibratoire d'un hélicoptère aux figures 4a à 4f. Le but est de surveiller la composante N2 du signal vibratoire y(t) de l'arbre haute pression d'une turbomachine 100, par exemple de l'arbre de la turbine libre dans le moteur de l'hélicoptère.

La figure 4a montre le spectrogramme du signal mesuré y(t) en indiquant pour chaque instant en abscisses et chaque fréquence sur la bande de fréquence [0,700] Hz en ordonnées un pixel, dont le niveau de gris correspond selon l'échelle de valeurs ECH de droite à une amplitude spectrale pour cette fréquence (et de manière analogue pour les spectrogrammes des figures 4c, 4d, 4e et 4f). La figure 4b montre en ordonnées l'amplitude de l'enveloppe a k (n ) ayant été calculée suivant l'invention en fonction du temps en abscisses (le bruit est considéré non-stationnaire et non gaussien et le paramètre de régularisation est estimé automatiquement). On observe à la figure 4b une augmentation significative de l'amplitude à partir de l'instant 5s, ce qui est signe d'un passage par une résonance. L'évaluation de l'amplitude maximale pendant cette résonance est de grande importance pratique pour le contrôle vibratoire de l'aéronef en régime transitoire. La figure 4c montre le spectrogramme du signal résiduel (le signal mesuré y(t) auquel on a soustrait la composante N2). On voit à la figure 4c qu'on a correctement estimé la composante N2 d'intérêt. Les figures 4d, 4e et 4f montrent les spectrogrammes du signal résiduel en utilisant la méthode connue de Vold-Kalman avec différents paramètres de régularisation réglés à la main (ce sont des valeurs couramment utilisées dans la littérature). On peut voir dans la figure 4d une sous-estimation de l'amplitude de la composante N2 d'intérêt au moment du passage par la résonance pour une grande valeur de r k (l'amplitude estimée de la composante N2 d'intérêt est moins importante que celle réelle, l'enveloppe estimée est très lissée d'une manière excessive, ce qui pénalise d'une façon très agressive les grandes variations de l'enveloppe). La figures 4e et 4f montrent qu'on a une surestimation de l'amplitude de la composante N2 d'intérêt, pour une petite valeur de r k (l'amplitude estimée de la composante N2 d'intérêt est plus importante que l'amplitude réelle, l'enveloppe estimée est très affectée par le bruit). En effet, la composante N2 d'intérêt ayant été estimée aux figures 4e et 4f est mal filtrée et contient beaucoup de bruit de fond. De plus, une suresrtimation de l'enveloppe peut résulter d'une matrice mal-conditionnée et donc de problèmes numériques (voir Figures 1 et 4f par exemple).

Suivant un autre mode de réalisation, β j est égal à zéro. Dans ce cas on omet la contribution du bruit e(n ) à bande étroite. L'avantage principal est la réduction du nombre de paramètres à estimer ainsi que le coût de calcul. Certes, cela peut être une bonne approximation si aucune harmonique parasite du signal résiduel n'interfère directement avec les signaux d'intérêt. Cela revient à

Dans ce cas particulier, l'avantage par rapport à la méthode connue de Vold- Kalman est le réglage automatique du paramètre y k de régularisation. En fait, le paramètre de Vold-Kalman r k est lié aux paramètres de l'invention par Un inconvénient de cette approximation pourra se manifester, par exemple, en présence d'interception entre les harmoniques d'intérêt et les harmoniques résiduelles. En effet, des artéfacts d'estimation peuvent avoir lieu aux instants d'interférences. Cependant, si on connaît approximativement les instants d'interférences (par exemple avec une inspection du spectrogramme), on peut supposer que le bruit est Gaussien sauf aux alentours de couplage, c'est-à-dire aux intervalles j ∈ (1, ... ,J} contenant les instants de couplages. Cela revient à imposer β j = 0 dans les autres intervalles ne contenant pas les moments de couplage. Bien entendu, les modes de réalisation, caractéristiques, possibilités et exemples décrits ci-dessus peuvent être combinés l'un avec l'autre ou être sélectionnés indépendamment l'un de l'autre.

Si on utilise un modèle parcimonieux autre que le modèle Bemoulli-Gaussien, on obtient les mêmes formes des lois. Seulement, les paramètres des modèles du bruit vont changer. La loi des enveloppes, des paramètres de régularisation et la variable auxiliaire auront la même forme. Le modèle du bruit choisi va seulement agir sur les valeurs de d j .