5. Application de la cyclostationnarité aux signaux de roulement

De nombreuses méthodes de diagnostic sur les roulements peuvent être trouvées dans la littérature [#!howard:review!#] :

• dans le domaine temporel, des méthodes statistiques peuvent être employées (valeur efficace, kurtosis, ...).,
• dans le domaine fréquentiel, il est possible d'observer l'évolution du spectre, ...,
• il est également possible d'utiliser l'analyse d'enveloppe qui consiste à employer un filtre passe bande pour isoler une composante issue du roulement, à démoduler cette dernière, puis à observer l'enveloppe du signal démodulé et son spectre,
• des développements récents ont tiré parti du kurtosis spectral [#!antoni:robust!#] [#!vrabie:application!#] ;
• de nombreuses autres méthodes existent : analyse cepstrale, analyse temps-fréquence, ainsi, que les statistiques d'ordres supérieurs.

Nous nous intéresserons dans ce chapitre à l'apport de la cyclostationnarité dans l'analyse des signaux de roulements : dans un contexte de pré-traitement puisque nous ne disposons pas de signaux avec une dégradation progressive.

Les méthodes de re-échantillonnage angulaire présentées au chapitre nous ont permis d'induire la propriété de cyclostationnarité sur des signaux de machines tournantes.

Les réducteurs comportent non seulement des roues mais également des roulements permettant d'assurer la liaison entre les axes et la carcasse métallique. Les roulements étant en contact avec les roues, tout signal accéléromètrique contient un mélange des contributions de ces deux éléments ainsi que du bruit. En général, la contribution des roues plus énergétique masque le signal de roulement. Il est alors nécessaire de séparer ces signaux préalablement au diagnostic.

Nous nous intéresserons dans un premier temps à l'élimination du signal issu de l'engrenage. Dans les réducteurs classiques, il suffit d'utiliser un simple filtre passe bande qui exploite le caractère ``basse fréquence'' des signaux issus des roues [#!mcfadden:vibration!#]. On peut ensuite utiliser des outils tels que l'analyse d'enveloppe pour mettre en évidence les défauts. Cependant, dans le cas des boîtes à vitesse à forte réduction (hélicoptères) où les roues génèrent un signal plus large bande, il est nécessaire de séparer les signaux à l'aide d'une autre méthode. Il est possible d'employer une technique de séparation non supervisée (appelée aussi technique de soustraction de bruit ou Self Adaptive Noise Cancellation par les initiateurs de la méthode) pour améliorer les résultats d'analyse d'enveloppe [#!ho:optimisation!#] (ici le bruit est le signal d'engrènement). Cette technique exploite le caractère périodique et donc prédictible des signaux d'engrenage par opposition aux signaux de roulement. On détermine alors un filtre capable de prédire l'échantillon courant à partir d'une série d'échantillons passés. Ce filtre permet alors d'estimer et donc de supprimer la contribution de l'engrenage.

La contribution des roues est périodique en angle et non en temps. Dès lors, cette technique est efficace uniquement lorsque les fluctuations de vitesse sont faibles, sinon ses performances se dégradent et la partie périodique ne peut plus être extraite. C'est pourquoi, nous avons amélioré cette méthode en la combinant au re-échantillonnage angulaire a posteriori [#!bonnardot:enhanced!#]. La figure montre le principe de cette méthode. Il aurait également été possible d'utiliser la moyenne synchrone afin d'estimer cette composante périodique [#!mcfadden:revised!#]. Néanmoins, la séparation non supervisée a été retenue puisque la moyenne synchrone ne peut séparer qu'une seule période cyclique à la fois et qu'elle nécessite la connaissance précise des fréquences cycliques.

Une fois la composante périodique retranchée, on obtient la contribution des roulements ainsi qu'un bruit additif.

Dans un deuxième temps, nous montrerons comment réduire ce bruit présent dans la contribution. Les signaux de roulement étant cyclostationnaires à l'ordre $2$ (sur une durée limitée), ils présentent une redondance spectrale (c'est-à-dire qu'il existe une corrélation entre certains canaux fréquentiels). Nous exploitons alors cette redondance en essayant de reconstruire le signal de roulement à partir de versions décalées en fréquence. Les décalages correspondant aux fréquences caractéristiques des roulements et à leurs harmoniques.

Figure: Principe de l'extraction du signal de roulement

5.1 Séparation fréquentielle

5.1.1 Exploitation par l'analyse d'enveloppe

L'analyse d'enveloppe [#!mcfadden:model!#,#!ho:optimisation!#] exploite la nature impulsive de la source $d\left(t\right)$ (voir équation ) : son énergie est distribuée sur tout le spectre et pondérée par la réponse impulsionnelle. La figure illustre cette technique :

Figure: Principe de l'analyse d'enveloppe

1. Dans les boîtes de vitesse classiques (à l'exception des très forts rapports de réduction comme dans les hélicoptères), la contribution des roues occupe une bande fréquentielle plus basse que celle des roulements. La contribution ``haute fréquence'' des roulements est donc extraite à l'aide d'un passe-bande $\left[f_1,f_2\right]$. Ce filtre peut être déterminé en examinant l'évolution du spectre entre la situation actuelle et une situation sans défaut. Comme le montre le modèle précédent, la source $d\left(t\right)$ passe par un filtre correspondant à la structure mécanique. Ainsi, le filtrage autour d'une résonance où l'atténuation par la structure mécanique est plus faible, permet d'exploiter au mieux la dynamique de la carte d'acquisition et des capteurs utilisés.

2. Le signal filtré est ensuite démodulé afin de le ramener dans la bande fréquentielle $\left[0,f_2-f_1\right[$. Pour cela, on calcule la transformée de Fourier du signal et on déplace la bande $\left[f_1,f_2\right]$ en $\left[0,f_2-f_1\right]$. Cette bande fréquentielle est alors doublée par zéro-padding avant de calculer la transformée de Fourier inverse, pour rendre le signal analytique. D'autres méthodes de démodulation, l'effet du zéro padding, ainsi que les effets de masquage sont étudiés plus en détail dans [#!ho:optimisation!#].

3. Le module du signal démodulé $dem\left(t\right)$ est ensuite extrait : $r\left(t\right)=\left\vert dem\left(t\right)\right\vert^m$ (avec $m=1$). En utilisant la transformée de Fourier (en bas à droite), il est alors possible d'étudier les propriétés fréquentielles de l'enveloppe.

A l'issue de ce traitement, on obtient une enveloppe ainsi que sa transformée de Fourier qui témoigne de l'état du roulement. Afin de pouvoir prendre une décision, il convient de comparer ce signal à un signal de référence ayant subi les mêmes traitements. Ce dernier doit être enregistré lorsque le roulement fonctionne parfaitement.

Plusieurs améliorations ont été proposées pour cette technique :

• Il est possible d'exploiter les propriétés homomorphiques du logarithme [#!randall:diagnostics!#]. Le logarithme permet de transformer les modulations intervenant dans le modèle (distance source capteur variable) en contributions additives. Alors, la détérioration de la résolution dans le spectre par la modulation est réduite.

• Si les motifs associés au défaut peuvent être localisés précisément dans le signal, il est possible de les isoler à l'aide d'une fenêtre. Dès lors, le bruit présent dans le reste du signal est supprimé.
  

La transformée de Fourier de l'enveloppe élevée au carré ($m=2$) correspond à l'intégrale de la corrélation spectrale selon les fréquences [#!randall:relationship!#]. En effet, on a  :

$$\int_{\mathbb{R}} \speccorr{x}\left(\alpha,f\right) df = \int_{\mathbb{R}} \lim_{W\rightarrow \infty} \frac{1}{W} \int_{\m... ...x}\left(t,\tau\right) \cdot e^{-2\pi j \left(f \tau+\alpha t\right)}dt d\tau df$$ (5.1)

$$= \lim_{W\rightarrow \infty} \frac{1}{W} \int_{\mathbb{R}} \int_{-W... ...ght) \cdot e^{-2\pi j \alpha t}dt \int_{\mathbb{R}} e^{-2\pi j f \tau} d\tau df$$ (5.2)

$$= \lim_{W\rightarrow \infty} \frac{1}{W} \int_{\mathbb{R}} \int_{-W... ...}\left(t,\tau\right) \cdot e^{-2\pi j \alpha t}dt \delta\left(\tau\right) d\tau$$ (5.3)

$$= \lim_{W\rightarrow \infty} \frac{1}{W} \int_{-W/2}^{W/2} \autocorr{x}\left(t,0\right) \cdot e^{-2\pi j \alpha t}dt$$ (5.4)

$$= \lim_{W\rightarrow \infty} \frac{1}{W} \int_{-W/2}^{W/2} \esp{\left\vert x\left(t\right)\right\vert^2} \cdot e^{-2\pi j \alpha t}dt$$ (5.5)

Où $\esp{\left\vert x\left(t\right)\right\vert^2}$ correspond à l'espérance de l'enveloppe. Ce spectre contient donc des pics aux fréquences caractéristiques des roulements (fréquences cycliques dues aux liens entre les harmoniques). Il contient également des pics correspondant aux autres fréquences cycliques (lien entre les bandes latérales par exemple).

5.1.2 Problèmes posés par l'analyse d'enveloppe

Cette technique est difficile à appliquer sur les signaux issus de réducteurs à engrenage d'hélicoptères où la contribution des roues est ``large bande'' (plus haute fréquence d'engrènement $11\;kHz$ - plus faible fréquence de rotation $4.3\;Hz$). C'est pourquoi nous allons présenter une technique de séparation exploitant l'aspect cyclostationnaire des signaux d'engrenage ainsi que la cyclostationnarité floue induite par les signaux de roulement (fréquence non multiple résultant du glissement).

5.2 Séparation engrènement roulement

5.2.1 Contribution des roues

Les vibrations des roues sont causées principalement par l'erreur de transmission. L'erreur de transmission est l'écart entre la position réelle de l'arbre de sortie et celle qu'il occuperait si l'engrenage était parfait. L'erreur de transmission est générée par la variation de raideur périodique des éléments constituant le réducteur et les défauts. Comme elle est pratiquement périodique, les vibrations engendrées (et donc l'accélération) sont périodiques et déterministes. La contribution des roues est donc prédictible à un instant donné à l'aide des valeurs passées.

La figure en bas montre la contribution des roues après séparation des contribution roues/roulement. Il est relativement facile de prédire le contenu du rectangle de gauche à partir de celui de droite. L'autocorrélation de ce signal, à droite, met bien en évidence la nature répétive de ce signal.

Figure: Décomposition du signal

5.2.2 Contribution des roulements

Il se produit un phénomène de glissement entre la cage intérieure et la cage extérieure du roulement à billes : en cas de défaut sur une bille, les impacts ne se reproduisent jamais à la même position. La charge étant liée à l'angle (figure ), chaque impact se produit dans des conditions différentes. Il existe aussi des différences (mineures) entre chaque bille. Tous ces phénomènes induisent des petites fluctuations aléatoires autour de la période de choc moyenne. Dès lors, la contribution produite par les roulements est aléatoire et donc imprévisible à l'aide des valeurs passées (voir figure ). Il est conclu dans [#!antoni:differential!#] que les signaux de roulement sont des signaux aléatoires non stationnaires.

La figure en haut montre la contribution des roulements. Il apparait difficile de prédire l'instant et surtout l'amplitude des ``chocs'' dans le rectangle de droite. L'autocorrélation de ce signal, à droite, souligne la nature aléatoire de ce signal.

5.2.3 Séparation par la méthode SANC (Self adaptive noise cancellation)

Le signal $x\left(t\right)$ d'un système engrenage-roulement peut être décomposé en deux parties : une partie périodique et prévisible $p\left(t\right)$ provenant des roues, et, une autre partie imprévisible $r\left(t\right)$ issue des roulements, nommée signal résiduel.

Le but de la soustraction de bruit non supervisée est de prédire la partie déterministe discrète de $x\left(t\right)$ : $\est{x}\left(n\right)$. Comme seule la contribution périodique est prédictible, $\est{x}\left(t\right)~=~\est{p}\left(t\right)$.

5.2.3.1 Approche temporelle

Soit $\vect{x}$ le vecteur contenant les $N$ valeurs précédentes de $x\left(n\right)$ :  $\vect{x}\left(n\right)~=~\left[x\left(n-N+1\right) \; \cdots \; x\left(n\right)\right]^T$. Le principe est d'estimer $\est{p}\left(t\right)$ à partir d'une version du signal décalée de $\Delta$ échantillon : $x\left(n-\Delta\right)$. Le délai $\Delta$ doit être choisi suffisamment grand pour supprimer toute corrélation entre le résidu et le signal retardé. Néanmoins, comme le signal n'est pas strictement périodique, $\Delta$ doit être suffisamment petit pour que les variations restent négligeables. On aura alors :

$$\hat{x}\left(n\right)=\vect{h}^T \cdot \vect{x}\left(n-\Del... ... la réponse impulsionnelle du filtre de taille $N\times 1$}$$ (5.6)

Comme le signal issu du roulement est imprévisible (on suppose que $\Delta$ est suffisamment grand pour supprimer toute corrélation locale), on peut alors l'estimer par :

$$\est{r}\left(n\right)=x\left(n\right)-\est{p}\left(n\right)$$ (5.7)

Théoriquement, la taille du filtre $N$ doit correspondre au moins au double du nombre de fréquences présentes dans $p\left(t\right)$. En présence de bruit, une valeur plus grande de $N$ permet d'augmenter la qualité des résultats.

La figure décrit la mise en uvre de cette méthode basée sur l'algorithme du gradient stochastique :

$$h^{n+1}=h^n+\mu \cdot \est{r}\left(n\right)\cdot \vect{x}\left(n-\Delta-N\right)$$ (5.8)

Cette approche [#!antoni:unsupervised1!#] nécessite le choix délicat du facteur d'oubli qui doit être suffisament faible pour garantir la convergence [#!widrow:adaptive!#,#!haykin:adaptive!#].

5.2.3.2 Approche fréquentielle

Ces équations peuvent être résolues plus rapidement à l'aide d'un algorithme fonctionnant dans le domaine fréquentiel [#!antoni:unsupervised2!#]. L'équation () peut s'écrire :

$$\interspectre{x\left(n\right)}{x\left(n-\Delta\right)}\left... ...ht)=H\left(f\right)\cdot \spectrev{x}{n-\Delta}\left(f\right)$$ (5.9)

où $\interspectre{x}{y}\left(f\right)$ désigne l'inter-spectre entre $x$ et $y$ (apodisation incluse), et, $H\left(f\right)$ désigne la transformée de Fourier du filtre $h\left(t\right)$.

Comme indiqué sur la figure , le signal est tout d'abord découpé en blocs avec recouvrement. On construit également une autre série décalée de $\Delta$. Des spectres et inter-spectres ``locaux'' sont calculés pour chaque couple de blocs. Les spectres et inter-spectres sont ``estimés'' par moyennage. On peut ensuite déduire le filtre en divisant l'inter spectre par le spectre. Pour cela nous supposons qu'il n'existe pas de composante fréquentielle de puissance nulle car les signaux réels contiennent du bruit.

Il est important de noter que par rapport à l'approche temporelle, l'approche fréquentielle évite le choix délicat du pas d'adaptation $\mu$.

Figure: Mise en oeuvre de l'algorithme SANC

5.2.4 Séparation non supervisée améliorée (E-SANC)

Face au problème posé par les fluctuations de vitesse pour la méthode SANC, l'échantillonnage angulaire semble assez intéressant. Comme les signaux sont rarement acquis suivant l'angle, nous avons proposé dans [#!bonnardot:enhanced!#] d'utiliser le re-échantillonnage a posteriori basé sur la démodulation autour d'une harmonique, avant d'appliquer la méthode $SANC$ comme l'illustre la figure . L'effet des fluctuations de vitesse est alors théoriquement annulé (voir chapitre ) et les résultats sont améliorés (on satisfait l'hypothèse de périodicité de la contribution des roues faite dans la SANC). Il est ensuite possible d'utiliser des méthodes de diagnostic [#!howard:review!#].

L'idée du re-échantillonnage a été également exploitée dans [#!groover:removal!#]. Néanmoins, la technique de re-échantillonage employée, beaucoup plus classique, fait appel à un signal tachymétrique (codeur optique). De plus, la suppression des composantes issues des roues n'est pas aveugle comme avec la méthode $SANC$ mais repose sur la suppression directe des pics dans la Transformée de Fourier.

Figure: Amélioration de la SANC

5.3 Application aux signaux d'hélicoptères (roulements planétaires)

Les signaux issus de roulements planétaires fournissent un des cas les plus difficiles pour la détection et le diagnostic des défauts. En effet, ces derniers passent par un chemin tortueux et changeant en fonction du temps, pour arriver aux points de mesure où ils sont détectés. Dans le cas des boîtes de vitesse d'hélicoptère, cette détection est rendue bien plus ardue car le bruit (généré entre autre par les roues du réducteur) est large bande (relativement au domaine acoustique). Par exemple, un réducteur convertit une fréquence de rotation d'environ $350\;Hz$ issue de la turbine à gaz en une fréquence de rotation de $5\;Hz$ [#!howard:review!#].

5.3.1 Origine des signaux

Nous allons maintenant appliquer ces méthodes sur des signaux accéléromètriques d'hélicoptères Sea Hawk. Ces derniers proviennent d'un banc d'essai du US Naval Air Warfare Center, Aircraft Division, Trenton (New Jersey). Le signal a été acquis à une fréquence d'échantillonnage $f_s=100\;kHz$ pendant une durée de $60\;s$ ($6\cdot 10^6$ échantillons) avec une charge de $400\;ft-lb$ ( $\approx 542\;Nm$).

La figure montre la transmission de l'hélicoptère ainsi que le défaut identifié sur un des rouleaux des roulements planétaires. Malheureusement la position du capteur accéléromètrique nommé ``PortRing'' n'a pas été fournie. Les fréquences caractéristiques du roulement planétaire sont indiquées dans le tableau .

Figure: Transmission de l'hélicoptère et défaut

Table: Fréquences caractéristiques du roulement planétaire (à gauche) et du système complet (à droite)

Désignation	Fréquence	tours/min
Défaut bague interne	$122 \;Hz$	$7322$
Défaut de bague externe	$90.4\;Hz$	$5422$
Défaut de bille	$77.359\;Hz$	$4642$
Vitesse de la cage (p)	$5.0\;Hz$	$300$
Vitesse de rotation	$11.8\;Hz$	$708$
Axe de la couronne (g)	$20.1\;Hz$	$1207$
Engrenement couronne (g) [75]	$1509\;Hz$	$90540$
Porte satellite (k)	$4.3\;Hz$	$258$

Note: le nombre entre crochets correspond au nombre de dents

5.3.2 Examen du signal accéléromètrique

La figure montre les $10\;000$ premiers échantillons ($100\;ms$) du signal. Les $7$ chocs (défaut de bille à $77.359\;Hz$) ne peuvent être détectés visuellement car la contribution du roulement est complètement masquée par le reste du signal (roue, autres roulements, bruit, ...). Une séparation des contributions du signal est donc nécessaire.

Figure: Premiers échantillons

La figure montre la densité spectrale de puissance du signal ainsi qu'un agrandissement. Les marqueurs correspondent à des fréquences de rotation, des fréquences d'engrènement ainsi qu'à leurs harmoniques. Les marqueurs ont été volontairement limités afin de ne pas surcharger la figure. Le but de cette figure n'est pas d'identifier avec précision les éléments associés à chacun des pics mais plutôt de montrer que les fréquences issues des roues occupent une large partie du spectre rendant difficile une analyse d'enveloppe sans débruitage préalable. On pourra noter la présence d'un marqueur à plus de $30\;kHz$.

Le spectre met en évidence un problème de repliement après $35\;kHz$, nous n'utiliserons donc pas ces ``hautes fréquences''.

Figure: Spectres

Nous allons tout d'abord utiliser la $SANC$, puis la comparer à l'$E-SANC$.

5.3.3 SANC

La séparation basée sur la $SANC$ nécessite de choisir deux paramètres : le décalage $\Delta$ et la longueur du filtre. Nous avons utilisé des écarts $\Delta$ variant de $10\;ms$ ($\Delta=1024$ échantillons) à $2.6\;s$ ( $\Delta=262\;144$ échantillons) avec une fenêtre pour le filtrage de $0.66\;s$ ($65\;536$ échantillons). Quelque soit la valeur de $\Delta$ utilisée, les résultats sont identiques. La figure compare le spectre avant et après la SANC (le signal original a été déplacé de $30\;dB$ pour une comparaison plus aisée). Il apparait que les résultats sont mauvais et que la SANC n'apporte rien dans ce cas. Plusieurs essais ont confirmé cet échec.

Figure: Spectre après la SANC

La figure montre l'estimation de la vitesse instantanée de la couronne, le filtre passe bande ($1.5\;kHz$ à $1.52\;kHz$) utilisé pour cette estimation supprime les fluctuations de vitesse supérieure à $20\;Hz$. La fluctuation de vitesse représente $0.9\;\%$ de la vitesse moyenne. Le spectre montre que ces fluctuations sont liées au porte satellite. L'élément mécanique associé aux fréquences de $6\;Hz$ et $11\;Hz$ n'a pas pu être mis en évidence.

Figure: Vitesse instantanée de crownwheel

Figure: Spectre après la ESANC

De telles fluctuations de vitesse peuvent expliquer le dysfonctionnement de la SANC et suggère donc l'utilisation de la E-SANC.

5.3.4 E-SANC

Le re-échantillonnage angulaire a posteriori a été effectué en utilisant la fréquence d'engrènement de l'arbre de la couronne à $1\;509\;Hz$. La SANC a été ensuite appliquée comme préconisé sur la figure .

La figure montre l'effet du débruitage sur le spectre après l'ESANC. Dans la zone basse fréquence, il y a une réduction des pics associés aux roues entre $5\;dB$ et $25\;dB$. Après $18\;kHz$, il y a juste une atténuation du signal, et aucune suppression de pics. La ESANC apporte donc une amélioration pour les basses fréquences. Dès lors, afin de bénéficier de ce débruitage, il est nécessaire de travailler dans la bande améliorée $\left[0;18\;kHz\right]$.

5.3.5 Analyse d'enveloppe

En général, la bande pour l'analyse d'enveloppe est choisie en examinant l'évolution du spectre pour différents niveaux de dégradation. Comme un seul enregistrement est disponible, nous avons essayé plusieurs filtres passe-bande avec une largeur de bande de $800\;Hz$ (cette bande doit être suffisament large pour permettre d'observer l'effet des modulations) et une fréquence centrale variant par pas de $500\;Hz$ (le pas a été choisi empiriquement). La bande de $600\;Hz$ à $1\,400\;Hz$ s'est avérée être la seule intéressante pour l'analyse d'enveloppe.

La figure montre l'apport bénéfique de la $ESANC$ sur l'analyse d'enveloppe, alors que les résultats de la $SANC$ sont similaires aux résultats sur le signal non traité. On peut observer que le pic associé au défaut est nettement plus visible. En effet, les fréquences associées aux roues qui prédominaient ont été fortement réduites. Ainsi, le diagnostic est facilité. Malheureusement, comme un seul signal est disponible, il n'a pas été possible d'observer l'évolution de ce pic en fonction du défaut.

Figure: Apport de l'ESANC sur l'analyse d'enveloppe

La figure montre la transformée de Fourier de l'enveloppe au carré dans la zone. Elle met en évidence à l'aide de traits pointillés les bandes latérales de $4.3\;Hz$ correspondant à la fréquence de rotation du porte satellite. Ces dernières peuvent être expliquées par la variation du chemin de transmission entre le défaut et le capteur à la fréquence du porte satellite.

Il apparaît que la fréquence du défaut ($77.359\;Hz$) est proche de la $18^\text{\tiny ème}$ harmonique de la fréquence de rotation du porte satellite ( $18\times4.3=77.4\;Hz$). Il est possible que cette harmonique aide à faire ressortir le pic de défaut : la fréquence de défaut (avec une fluctuation due au glissement) serait si proche de l'harmonique de la fréquence de rotation du porte satellite qu'elle se ``synchroniserait'' sur cette dernière plus énergétique. Une fois la partie périodique enlevée, seule la contribution fluctuante du défaut subsisterait. Ce lien entre la fréquence de défaut et le $18^\text{\tiny ème}$ harmonique du porte satellite reste tout de même assez gênante pour faire du diagnostic sur ce type d'hélicoptère. Néanmoins, nous avons d'une part montré que l'ESANC améliore nettement l'analyse d'enveloppe par rapport à la SANC quand il y a des fluctuations de vitesse non négligeables, et d'autre part nous avons testé avec succès la méthode de re-échantillonnage a posteriori sur des signaux d'hélicoptère réputés complexes.

Figure: Transformée de Fourier de l'enveloppe au carré du signal après ESANC

La SANC et l'ESANC ont pour but de supprimer la partie périodique du signal, mais, le signal résiduel contiend toujours du bruit en plus du signal provenant des roulements. Dans le paragraphe suivant nous utilisons la cyclostationnarité à l'ordre $2$ pour atténuer ce bruit.

5.4 Débruitage basé sur la cyclostationnarité à l'ordre 2

Le débruitage basé sur l'exploitation des propriétés de cyclostationnarité à l'ordre $2$ permet de réduire le bruit additif toujours présent dans le signal résiduel mais ne permet pas pas de supprimer la composante périodique générée par l'engrenage. Dès lors, il conviendra d'extraire le signal résiduel à l'aide de la ESANC, par exemple, avant d'employer cette méthode. On travaillera alors sur un signal qui sera non cyclostationnaire à l'ordre $1$. Nous utiliserons à présent $x\left(\theta\right)$ pour désigner le signal résiduel.

Tout comme la séparation à l'ordre $1$ qui exploitait les similitudes d'un cycle à l'autre, le débruitage à l'ordre $2$ exploitera les similitudes d'une fréquence cyclique à l'autre. Pour faire le débruitage à l'ordre $2$, nous allons utiliser un algorithme issus du domaine des télécommunications portant le nom de FRESH [#!gardner:cyclic!#].

Nous travaillerons sur le modèle :

$$x\left(t\right)=c\left(t\right)+b\left(t\right)$$ (5.10)

Où $c\left(t\right)$ est le signal issu des roulements. Il est supposé suffisament court ($45$ fois la période moyenne du défaut dans notre exemple) pour justifier le choix du $1^{er}$ modèle. Dans ce cas, il est possible de définir une période moyenne autour de laquelle le défaut fluctue (fluctuation très inférieure à la période). Le terme $b\left(t\right)$ représente un bruit qui contient le reste du signal. Nous considérons comme bruit tout ce qui n'est pas corrélé avec $c\left(t\right)$. Nous rapellons l'expression de $c\left(t\right)$ définie dans l'équation () :

$$c\left(t\right)=r\left(t,\tau\right) \ast \left[ A\left(t\right) \sum_i \dirac{t}{\text{\tiny$T_i$}} \right]$$ (5.11)

Nous ne considérons pas les ordres supérieurs à $2$. En effet, la méthode employée exploite un critère d'énergie (ordre $2$), elle est donc incapable de faire la différence entre deux candidats possédant une erreur de même énergie mais de dissymétries différentes.

5.4.1 Mise en évidence des liens à l'ordre 2

La cyclostationnarité à l'ordre $2$ se traduit par l'existence de liens entre les fréquences. Parmi toutes les représentations existantes, la plus adaptée pour faire apparaître ces liens est la densité de corrélation spectrale. La corrélation spectrale est définie par :

$$\speccorr{x}\left(f,\alpha\right)df d\alpha=\esp{dX\left(f+... ...c{\alpha}{2}\right)\cdot dX^*\left(f-\frac{\alpha}{2}\right)}$$ (5.12)

Pour l'évaluer, on calcule une corrélation entre deux incréments spectraux du signal préalablement décalées en fréquence : $X\left(f+\frac{\alpha}{2}\right)$ et $X\left(f-\frac{\alpha}{2}\right)$. La valeur $\alpha$ correspond à la fréquence cyclique. Pour $\alpha=0$, la corrélation spectrale correspond au spectre du signal $x\left(t\right)$.

La figure montre le module de la transformée de Fourier d'un signal issu du modèle de roulement () simplifié (où le filtre $r\left(t,\tau\right)$ a été remplacé par un filtre $r\left(\tau\right)$). Pour le simuler, on génère au hasard les positions des impulsions $T_i$ conformément à l'équation (). Le caractère aléatoire est introduit à l'aide d'une loi normale centrée pour les jitters $e_i$. On génère ensuite une série aléatoire $A_i$ correspondant à l'amplitude associée à la $i^{\text{\tiny ème}}$ impulsion suivant également une loi normale. Ensuite, on génère la série de pics d'amplitude $A_i$ à la position $T_i$. Puisque les positions ne sont pas entières, chaque pic est obtenu à l'aide d'un décalage circulaire d'un pic à l'origine. Cette position non entière permet de modéliser l'effet d'un échantillonnage des signaux non synchrones de la position des pics. On filtre ensuite le signal par la réponse impulsionnelle de la structure que nous avons modélisée par deux résonances. On ajoute également un bruit blanc.

Figure: Spectre d'un signal de roulement (synthétisé)

Pour la figure , nous avons utilisé une période moyenne $T_{moy}$ de $360$ échantillons avec une fluctuation $e_i$ d'écart type $10$ échantillons. L'amplitude moyenne est $1$ avec un écart type $0.1$. La puissance du bruit est $-7\;dB$. Les résonances sont placées aux fréquences réduites $0.11\;Hz$ et $0.16\;Hz$ et correspondent respectivement à un amortissement de $0.08$ et $0.06$. Le signal comporte environ $300$ cycles.

On retrouve nettement sur le spectre les deux résonances correspondant à la réponse impulsionnelle de la structure. On retrouve également une série de pics espacés de $1/360\;Hz$ et une distribution sur ces pics correspondant à la fluctuation autour de la fréquence moyenne.

Pour le calcul de la corrélation spectrale, cette transformée de Fourier est décalée respectivement de $-\alpha/2$ et de $+\alpha/2$. Lorsque la fréquence cyclique $\alpha$ est multiple de la période moyenne du défaut, les ``pics'' sont face à face et l'on obtient une énergie plus importante que lorsque les pics sont disjoints (un lien existe entre chacun des pics puisque le signal est cyclostationnaire flou). Il apparaît alors des ``lignes'' dans la corrélation spectrale comme l'illustre la figure . Ces lignes sont espacées de la fréquence cyclique $\alpha$ associée à la période moyenne du défaut. On retrouve ainsi le résultat établi pour le $1^{er}$ modèle (jitter) dans [#!randall:relationship!#] : la corrélation spectrale est continue suivant $f$ et discrète suivant $\alpha$.

Figure: Corrélation spectrale d'un signal de roulement (synthétisé)

5.4.2 Principe de la reconstruction

Ce paragraphe montre comment on peut reconstruire le signal à partir de versions décalées en fréquence puis filtrées.

La figure est l'élément clef que nous allons exploiter pour le débruitage : puisqu'il existe un lien entre le signal et ses versions décalées en fréquence, nous allons essayer de l'estimer à partir de plusieurs versions décalées, par un filtre de Wiener Cyclique [#!gardner:cyclostationarity!#].

La figure , montre sur un exemple très simple le principe de ce filtre. Le signal $x\left(t\right)$ de transformée de Fourier $X\left(f\right)$ est construit à partir de deux motifs aléatoires et décorrellés :

• un bruit rose de bande limitée centré sur la fréquence $f_1$ et modulé par un sinus à la fréquence $f_3-f_1$,
• un bruit blanc de bande limitée centré sur $f_2$ et modulé par un sinus à la fréquence $f_3-f_2$.

Un bruit blanc (large bande) est ensuite ajouté. Le choix des motifs (rectangle pour le bruit blanc à bande limitée, et demi-disque pour le bruit rose à bande limitée) sur la figure permet de bien mettre en évidence le lien entre les différentes fréquences.

En général, le but du filtre est de reconstruire le signal sans le bruit blanc (large bande). Nous supposerons que les fréquences cycliques sont connues. En effet, bien qu'en théorie il soit possible d'exploiter tous les décalages, en pratique, le traitement d'un décalage se révèle assez coûteux en temps de calcul. D'après la figure, les fréquences cycliques (ou décalages) à exploiter pour la reconstruction sont $\vect{\alpha/2}=~\left[f_1+f_1,f_1+f_2,f_1+f_3,-f_2,-f_3\right]$.

Figure: Principe du filtre de Wiener Cyclique

Chaque signal décalé est ensuite filtré (un filtre différent par signal décalé - voir colonne $\left\vert G_i\left(f\right)\right\vert^2$ de la figure ). Les signaux résultants (dernière colonne) sont ensuite sommés pour obtenir une estimation du signal débruité $c\left(t\right)$. Dans le cas général, les filtres sont estimés par une procédure proche du filtrage de Wiener classique (minimisation de l'écart quadratique entre le signal estimé et le signal bruité $x\left(t\right)$).

Les filtres ont deux rôles :

• la mise en forme des signaux (dans l'exemple il s'agit d'un simple gain),
• supprimer le bruit où il n'y a pas de signal corrélé.

Après le décalage fréquentiel, le bruit associé à chacun des motifs est différent. Ainsi, la combinaison des signaux décalés-filtrés est destructive pour le bruit (effet de moyennage). Il sera donc intéressant d'utiliser le maximum de fréquences cycliques ainsi que leurs harmoniques (en pratique on sera limité par le temps et la mémoire nécessaire). Dans le cas général, le décalage nul (c'est-à-dire le signal original) ne devra pas être utilisé lors de l'estimation car après filtrage on retrouvera une estimation du signal original.

Il est également important de noter que cette méthode ne supprime pas le bruit mais le réduit (il s'agit donc d'une estimation de $c\left(t\right)$).

La majorité des remarques précédentes sont communes avec le filtre de Wiener classique.

La figure montre le filtre de Wiener cyclique sous la forme d'un diagramme. L'utilisation des cosinus ou sinus permet de travailler sur des filtres réels et ne nécessite pas l'utilisation de fréquences négatives. Puisque $x\left(t\right)$ est réel, ces deux formulations sont équivalentes :

0

$$\Re\left[x\left(n\right) e^{2\pi j\alpha_k n} \ast g\left(t\right)\right] = \Re \left[x\left(n\right) \cos\left(2\pi \alpha_k n\right) \ast g... ...j x\left(n\right) \sin\left(2\pi \alpha_k n\right) \ast g\left(n\right) \right]$$ (5.13)

$$= x\left(n\right) \cos\left(2\pi \alpha_k n\right) \ast h\left(n\right) + x\left(n\right) \sin\left(2\pi \alpha_k n\right) \ast l\left(n\right)$$ (5.14)

Avec

$$g\left(n\right)=h\left(n\right)-j l\left(n\right)$$ (5.15)

Ce diagramme permet de mettre en évidence le filtrage LPTV (Linéaire Variant Périodiquement dans le Temps). Plusieurs représentations équivalentes des filtres LPTV ont été recensées dans [#!ferrara:frequency!#], néanmoins, comme la partie la plus coûteuse en temps de calcul est la détermination des filtres et non le filtrage en lui-même nous ne comparerons pas ces implémentations.

Figure: Implémentation du filtre de Wiener Cyclique

Nous allons présenter deux méthodes pour la recherche des coefficients. Ces deux méthodes sont basées sur la minimisation de l'erreur quadratique entre le signal prédit $\est{c}\left(n\right)$ et le signal $c\left(n\right)$. La première méthode minimise ce critère dans le domaine temporel alors que la deuxième fait la minimisation dans le domaine fréquentiel.

5.4.2.1 Approche temporelle

La figure peut se formaliser par l'équation suivante :

$$\est{c}\left(n\right)=\sum_{k=1}^{K}\left\{h_k\left(n\right... ...lpha_k n\right)\right]\right\} \text{pour $\alpha_k\neq 0$}$$ (5.16)

Où $h_k\left(n\right)$ et $l_k\left(n\right)$ sont les filtres optimaux au sens des moindres carrés associés à la $k^\text{\tiny ième}$ fréquence cyclique $\alpha_k$. Les fréquences cycliques $\alpha_k$ correspondent, par exemple, dans le cas d'un signal de roulement à la fréquence du défaut $1/T_{moy}$ ainsi qu'a ces harmoniques. En général, on utilise toutes les fréquences pour lesquelles la corrélation spectrale a établi un lien.

Afin de minimiser l'erreur quadratique $\varepsilon_x^2\left(n\right)=\left\vert c\left(n\right)-\est{c}\left(n\right)\right\vert^2$, nous allons maintenant chercher à formuler l'équation sous forme matricielle.

Le filtre LPTV (figure ) peut être vu comme un système MISO (entrée multiple sortie unique) dont les entrées sont :

$$\left\{\begin{array}{lcl} a_k\left(n\right) & = & x\left(n... ...\cdot \sin\left(2\pi \alpha_k n\right). \end{array} \right.$$ (5.17)

Afin d'utiliser la formulation vectorielle du filtrage, nous allons définir les vecteurs associés au filtre $h_k\left(n\right)$ et $l_k\left(n\right)$ de longueur $L=\tau_L-\tau_1+1$ :

$$\left\{ \begin{array}{lcl} \vect{h}_k & = & \left[h_k\lef... ...dots \¤ l_k\left(\tau_L\right)\right]^T \end{array} \right.$$ (5.18)

Où les limites du filtre $\tau_1$ et $\tau_L$ sont choisies par l'utilisateur. Cette notation permet notamment de concevoir des filtres non causaux. En conséquence, il est nécessaire de mémoriser les valeurs passées et futures des entrées :

$$\left\{ \begin{array}{lcl} \vect{a}_k\left(n\right) & = &... ...dots \¤ b_k\left(n-\tau_1\right)\right] \end{array} \right.$$ (5.19)

L'estimée $\est{c}\left(n\right)$ s'écrit donc :

$$\est{c}\left(n\right)=\sum_{k=0}^K \left[\vect{a}_k\left(n\... ... \vect{h}_k + \vect{b}_k\left(n\right)\cdot \vect{l}_k\right]$$ (5.20)

Où sous forme vectorielle :

$$\est{c}\left(n\right) = \vect{r} \cdot \vect{s}$$ (5.21)

Avec :

$$\vect{r} = \left[\vect{a}_1\left(n\right) \¤ \cdots \¤ \vect{a}_K\left(n\right) \¤ \vect{b}_1\left(n\right) \¤ \cdots \¤ \vect{b}_K\left(n\right)\right]$$ (5.22)

$$\text{et } \vect{s} = \left[\vect{h}_1^T \¤ \cdots \¤ \vect{h}_K^T \¤ \vect{l}_1^T \¤ \cdots \¤ \vect{l}_K^T \right]^T$$ (5.23)

Cette équation comporte $u=2KL$ inconnues (coefficients du filtre), il faut donc au moins $u$ équations soit $N_{min}=u+\tau_L-\tau_1$ échantillons pour trouver le filtre (un plus grand nombre d'équations permettra de faire une régression et donc d'améliorer les résultats). Pour un signal de $N$ échantillons, ($N>N_{min}$), nous avons :

$$\est{\vect{c}} = \mat{m}\cdot \vect{s}$$ (5.24)

Avec :

$$\est{\vect{c}}=\left[\est{c}\left(N+\tau_1\right) \¤ \cdots \¤ \est{c}\left(1+\tau_L\right)\right]^T$$ (5.25)

$$\mat{m} = \left[ \begin{array}{c} \vect{r}\left(N+\tau_... ...s \\ \vect{r}\left(\tau_L+1\right) \end{array} \right]$$ (5.26)

La minimisation de $\varepsilon_x\left(n\right)^2$ conduit donc au système d'équation :

$$\mat{\autocorr{m}} \cdot \vect{s}=\vect{\crosscorr{\mat{m}}{\vect{c}}}$$ (5.27)

Où $\mat{\autocorr{m}}=\mat{m}^H\cdot \mat{m}$, et $\vect{\crosscorr{\mat{m}}{\vect{c}}}=\mat{m}^H\cdot \vect{c}$.

Cette équation ne peut pas être résolue directement puisque $c\left(n\right)$ est inconnu. Le bruit $b\left(n\right)$ est décorrélé du signal $x\left(n\right)$ et de ces versions décalées. Dès lors, le vecteur de corrélation pourra être estimé en substituant le signal $x\left(n\right)$ à $c\left(n\right)$.

$$\vect{\crosscorr{\mat{m}}{\vect{x}}} = \vect{\crosscorr{\ma... ...rr{\mat{m}}{\vect{b}}} = \vect{\crosscorr{\mat{m}}{\vect{c}}}$$ (5.28)

Le système d'équation à résoudre devient donc :

$$\mat{\autocorr{m}} \cdot \vect{s}=\vect{\crosscorr{\mat{m}}{\vect{x}}}$$ (5.29)

Ce système peut être résolu à l'aide de la pseudo inverse.

L'erreur de prédiction est donnée par :

$$\varepsilon^2 = \autocorr{c}\left(0\right) \left[1-\rho\right]$$ (5.30)

$$\text{avec } \rho = \frac{\vect{\crosscorr{\vect{r}}{c}} \mat{\autocorr{m}}^{-1} \vect{\crosscorr{\mat{m}}{\vect{c}}}}{\autocorr{c}\left(0\right)}$$ (5.31)

5.4.2.2 Approche fréquentielle

Pour éviter d'alourdir les notations, nous avons écrit $X\left(f\right)$ au lieu de $dX\left(f\right)$. Il est également possible de reformuler les équations () dans le domaine fréquentiel :

$$X\left(f\right) = C\left(f\right) + B\left(f\right)$$ (5.32)

Où $X\left(f\right)$, $C\left(f\right)$, $B\left(f\right)$ sont les transformées de Fourier discrètes des signaux $x\left(n\right)$, $c\left(n\right)$, $b\left(n\right)$ de durée finie.

La figure peut alors se formaliser par l'équation suivante :

$$\hat{C}\left(f\right)=\sum_{k=1}^K G_k \left(f\right) \cdot X\left(f-\alpha_k\right)$$ (5.33)

Où $G_k\left(f\right)$ est la réponse fréquentielle du filtre associé au décalage en fréquence $\alpha_k$.

Il est possible de reformuler cette équation sous forme matricielle en créant des vecteurs étendus :

$$\vect{X}_e\left(f\right)=\left[X\left(f-\alpha_1\right) \¤ ... ...\left[G_1\left(f\right) \¤ \cdots \¤ G_k\left(f\right)\right]$$ (5.34)

On aura donc :

$$\hat{C}\left(f\right) = \vect{G}\left(f\right) \cdot \vect{X}_e\left(f\right)$$ (5.35)

$$\varepsilon\left(f\right) = C\left(f\right) - \est{C}\left(f\right) = C\left(f\right) - \vect{G} \left(f\right) \cdot \vect{X}_e \left(f\right)$$ (5.36)

Nous voulons minimiser $\varepsilon^2\left(f\right)=\esp{\varepsilon^* \left(f\right)\cdot \varepsilon\left(f\right)}$. En utilisant le théorème de la projection orthogonale qui établit que l'erreur $\varepsilon\left(f\right)$ doit être orthogonale à l'entrée $\vect{X}_e\left(f\right)$ et après quelques calculs, on obtient :

$$\mat{\spectre{\vect{x}_e}}\left(f\right) \cdot \vect{G}\left(f\right) = \vect{\interspectre{c}{\vect {x}_e}}\left(f\right)$$ (5.37)

Avec :

$$\vect{\interspectre{c}{\vect {x}_e}}\left(f\right) = \esp{C... ...sp{\vect{X}_e\left(f\right) \cdot \vect{X}_e^H\left(f\right)}$$ (5.38)

Où $\vect{\interspectre{c}{\vect {x}_e}}\left(f\right)$ est un vecteur $\left(1 \times K\right)$, et, $\mat{\spectre{\vect{x}_e}}\left(f\right)$ est la matrice inter-spectrale de taille $\left(K \times K\right)$.

La recherche du filtre optimal se fait donc en adaptant les équations de Wiener-Hopf dans le domaine fréquentiel. Comme dans le domaine temporel, $C\left(f\right)$ est inconnu. La décorrélation entre le bruit et le signal recherché sera exploitée afin de substituer $X\left(f\right)$ à $C\left(f\right)$ :

$$\esp{X\left(f\right) \cdot X^*\left(f-\alpha_k\right)} = \esp{\left[C\left(f\right)+N\left(f\right)\right]\cdot X^*\left(f-\alpha_k\right)} \¤ \alpha_k \neq 0$$ (5.39)

$$= E\left\{C\left(f\right)\cdot X^*\left(f-\alpha_k\right)\right\}$$ (5.40)

Cette équation montre l'importance de ne pas inclure la fréquence cyclique nulle (c'est-à-dire le signal original) dans l'algorithme de débruitage. Le filtre optimal est solution de l'équation :

$$\mat{\spectre{\vect{x}_e}}\left(f\right) \cdot \vect{G}\lef... ...nterspectre{x}{\vect {x}_e}}\left(f\right) \¤ \alpha_k \neq 0$$ (5.41)

L'erreur de prédiction associée à une fréquence donnée est :

$$\varepsilon^2\left(f\right) = \spectre{c}\left(f\right)\left[1-\rho\left(f\right)\right]$$ (5.42)

$$\text{avec }\rho\left(f\right) = \frac{\vect{\interspectre{c}{\vect {x}_e}}\left(f\right) \mat{\sp... ... \vect{\interspectre{\vect {x}_e}{c}}\left(f\right)}{\spectre{c}\left(f\right)}$$ (5.43)

A la différence de la méthode temporelle, la méthode fréquentielle effectue la minimisation sur chaque canal fréquentiel indépendamment. Une seule minimisation étant faite dans le domaine temporel, il serait nécessaire de blanchir les signaux afin de s'assurer d'une pondération identique de chacune des fréquences.

5.5 Application de la séparation à l'ordre 2

Nous allons d'abord traiter un exemple à vocation didactique avec des signaux très simples. Ils permettront notamment de construire des filtres intuitifs et de les comparer aux filtres optimaux. Ensuite, nous observerons les résultats sur un modèle de signal de roulement. Enfin, l'algorithme sera testé sur des signaux réels.

Tous les signaux exploités ici sont soit à moyenne synchrone nulle, soit des signaux auxquels on aura au préalable retiré la moyenne synchrone.

5.6 Cas simple

5.6.1 Modèle

Soit le signal :

$$x\left(n\right) = b_1\left(n\right) \cdot \left[\lambda + \sin\left(2\pi f_1 n\right)\right] + b_2\left(n\right)$$ (5.44)

$$= c\left(n\right) + b_2\left(n\right)$$ (5.45)

Où :

$$\left\{ \begin{array}{ll} b_1\left(n\right) & \text{est u... ...blanc et stationnaire}, \\ \lambda=2. \end{array} \right.$$ (5.46)

La différence entre la puissance de $c\left(n\right)$ et $b_2\left(n\right)$ est de $3\;dB$.

Ce signal est généré à partir d'un bruit à bande limité qui est modulé en fréquence et également multiplié par un gain. Cette modulation assure la présence d'un lien à la fréquence cyclique $\alpha=f_1$ entre la contribution à $f_1$ et la contribution non modulée. L'utilisation d'un bruit permet de garantir la non cyclostationnarité à l'ordre 1. Un bruit additif est ensuite ajouté. Les figures g et h montrent respectivement le spectre de chacune des composantes ( $c\left(n\right)$ et $b_2\left(n\right)$) ainsi que le spectre du signal.

L'équation a été volontairement limitée à une seule fréquence cyclique pour réduire le nombre de filtres à étudier (les paragraphes suivant traiteront de signaux plus complexes). Le bruit $b_1\left(n\right)$ a également été filtré pour pouvoir d'une part le discerner dans le spectre avec plus de facilité et d'autre part observer l'effet de l'algorithme dans les bandes fréquentielles ne contenant aucun signal ``utile''.

Avant d'observer les résultats de l'algorithme, nous allons déterminer un filtre empirique qui sera ensuite comparé aux filtres optimaux (au sens des moindres carrés).

5.6.2 Choix du filtre

Après le décalage en fréquence du signal par un sinus et un cosinus à la fréquence $f_1$ (figure ), nous avons :

$$x\left(n\right) \cdot \vect{md_1}\left(n\right) = b_1\left(... ...\right) +b_2\left(n\right) \cdot \vect{md_1}\left(n\right)$$ (5.47)

Avec :

$$\vect{md_m}\left(n\right)=\left[\cos \left(2\pi m f_1 n\right) \¤ \sin\left(2\pi m f_1 n\right)\right]^T$$ (5.48)

Où $md_1\left(n\right)$ représente les fonctions cosinus et sinus utilisées pour moduler le signal et $md_2\left(n\right)$ les harmoniques générées par cette modulation. Comme cette expression est ``simple'' et que $b_1\left(n\right)$ a été choisi de telle sorte qu'il n'y ait pas de recouvrement fréquentiel entre ses versions modulées, il est possible de déterminer les filtres nécessaires à la reconstruction du signal en observant l'équation (). La reconstruction sera basée sur le signal modulé par le sinus. Si l'on introduit un déphasage, il sera indispensable d'exploiter le signal modulé par le sinus et le cosinus.

La figure montre le signal original en trait continu ainsi que le signal décalé en pointillé. Le filtre a trois rôles indiqués par des flèches sur le schéma :

1. Supprimer le bruit situé à la fréquence $2f_1$,
2. Fournir les gains nécessaires pour retrouver la même forme que le signal initial (soit $2\lambda$ en ``basse fréquence'' pour transformer le terme $\frac{b_1\left(n\right)}{2}$ en $\lambda\cdot b_1\left(n\right)$ et $1/\lambda$ en ``haute fréquence'' pour transformer le terme $b_1\left(n\right) \lambda \sin\left(2\pi f_1 n\right)$ en $b_1\left(n\right) \sin\left(2\pi f_1 n\right)$,
3. Supprimer le bruit modulé $b_2\left(n\right)$ pour les fréquences où il est disjoint de $c\left(n\right)$.

Figure: Choix du filtre

Cet exemple améliore uniquement le signal pour les fréquences supérieures à $f_1/4$. Néanmoins, le bruit sera amplifié ailleurs. Ceci est du à la simplicité de l'exemple. En effet, pour donner de bons résultats, cette méthode nécessite l'utilisation de plusieurs fréquences cycliques.

5.6.3 Simulation

La figure montre le résultat de l'algorithme temporel sur des signaux synthétiques pour un filtre de longueur $L=200$. Les courbes bleues correspondent aux estimations, les courbes vertes, rouges et noires correspondent respectivement au signal $x\left(n\right)$ avant traitement, à la composante $c\left(n\right)$, et, du bruit ajouté. La courbe (a) montre 10 cycles du signal $x\left(n\right)$, la courbe (b) montre la variance synchrone. En théorie, la caractérisation de la cyclostationnarité à l'ordre $2$ nécessite l'observation d'une courbe en trois dimensions (corrélation spectrale, corrélation, ...). Nous avons néanmoins choisi d'afficher la variance synchrone pour caractériser la cyclostationnarité à l'ordre $2$ car elle est très proche de l'analyse d'enveloppe (pour des signaux centrés, seul un moyennage est ajouté). Si la puissance du bruit $b_i\left(n\right)$ est $\sigma_i^2$, l'expression de la variance synchrone est :

$$vs\left(n\right)=\sigma_{b_1}^2\left[\lambda^2+2\lambda\sin... ...}-\frac{\cos\left(4\pi f_1 n\right)}{2}\right]+\sigma_{b_2}^2$$ (5.49)

Sur les figures b et d, on peut constater que la variance synchrone est périodique mais qu'il subsiste une légère modulation (due à la lente convergence de son estimateur).

Les figures c et e montrent les signaux après débruitage. On peut remarquer que la variance synchrone du bruit estimé (figure f) n'est pas constante mais qu'il subsiste encore un terme cyclostationnaire (fluctuations périodiques). Le terme cyclostationnaire se trouve introduit artificiellement : si on utilise le filtre $h\left(n\right)$ déterminé au paragraphe (utilisant uniquement la composante sinusoïdale), l'équation () montre que le bruit $b_2\left(n\right)$ est également modulé à la fréquence $f_1$, aussi, si nous notons $\est{b}_2\left(n\right)$ le bruit estimé, nous avons :

$$\est{b}_2\left(n\right)=b_2\left(n\right)-\left[b_2\left(n\... ...\cdot \sin\left(2\pi f_1 n\right)\right] \ast h\left(n\right)$$ (5.50)

Alors, un lien à la fréquence cyclique est créé entre le bruit et sa version décalée qui explique la composante cyclostationnaire présente dans l'estimation du bruit $\est{b}_2\left(n\right)$. Ce défaut n'est pas très gênant si l'on s'intéresse uniquement à $c\left(n\right)$. En effet, $\est{c}\left(n\right)$ contiendra le même type de bruit mais décalé en fréquence et atténué. Néanmoins, il faudra être prudent lors de l'interprétation de $\est{b}_2\left(n\right)$.

La figure g montre le spectre des différentes composantes du signal. La figure h montre le signal estimé par l'algorithme. L'algorithme réduit bien le bruit dans la zone fréquentielle contenant uniquement du bruit (fréquence au delà de $0.04\;Hz$). Ailleurs, le signal reconstruit a un spectre similaire au signal original.

La figure i montre les filtres utilisés pour la reconstruction. Contrairement à notre démarche du paragraphe précédent (), les filtres n'ont pas été déterminés en connaissant le modèle mais en minimisant l'erreur quadratique. Cette différence de critère pourra expliquer une différence entre le filtre attendu et le filtre réel. En observant les filtres, on note une atténuation de $40\;dB$ dans la zone $\left[0,0.01\right]\;Hz$ ainsi qu'une atténuation de $60\;dB$ vers $0.02\;Hz$, soit une différence de $20\;dB$. Cet écart est conforme à notre gabarit de la figure puisque l'atténuation doit être de $20 \log_{10}\left(2\lambda^2\right)$ c'est-à-dire de $18\;dB$ avec $\lambda=2$. L'atténuation de $40\;dB$ résulte de l'implémentation du filtre.

Figure: Débruitage dans le cas simple

5.7 Influence de l'allure du signal

Le but de ce paragraphe est d'examiner l'influence de l'allure du spectre du signal sur la qualité de la reconstruction ou plus simplement du facteur $\lambda$ du paragraphe précédent. Pour évaluer cette influence, on utilise un signal de la forme :

$$x\left(n\right) = c\left(n\right)+b_2\left(n\right)/2$$ (5.51)

$$c\left(n\right) = b_1\left(n\right) \left[0.5+ \sin\left(2\pi f_1 n\right) + \lambda \cdot \sin\left(2\pi f_2 n\right)\right]$$ (5.52)

Ce signal contient un terme modulé supplémentaire par rapport au paragraphe précédent comme l'indique son spectre (figure ). Nous allons évaluer l'erreur de reconstruction $\varepsilon\left(n\right)=\est{c}\left(n\right)-c\left(n\right)$ dans la bande fréquentielle associée au motif central $[0.08;0.12]\;Hz$. Lors de la simulation, cette erreur sera mesurée par moyennage du spectre de l'erreur dans cette bande fréquentielle.

La fréquence centrale sera estimée à partir des deux autres motifs, on aura donc dans la bande fréquentielle $B=[0.08;0.12]\;Hz$ après les décalages, filtrages et combinaisons :

$$\varepsilon_B = \left[\frac{b_1\left(n\right)+b_3\left(n\right)}{2} + \frac{\lamb... ...4\left(n\right)}{2\lambda} -b_1\left(n\right)\right]\cos\left(2\pi f_1 t\right)$$ (5.53)

$$= \left[\frac{b_3\left(n\right)}{2} + \frac{b_4\left(n\right)}{2\lambda} \right] \cos\left(2\pi f_1 t\right)$$ (5.54)

Où $b_3\left(n\right)$ et $b_4\left(n\right)$ représentent les bruits modulés lors du filtrage $LPTV$. Comme ces deux bruits sont modulés à des fréquences différentes, leur intercorrélation est nulle. Leur écart-type $\sigma_b$ est identique. En exploitant ces deux propriétés, on peut écrire :

$$\esp{\varepsilon_B^2}=\sigma_B^2\left(\frac{1}{4}+\frac{1}{4\lambda^2}\right)$$ (5.55)

Dès lors, plus le facteur $\lambda$ sera important plus le bruit associé au motif central reconstruit sera faible. La réduction ne pourra néanmoins excéder $6\;dB$. Pour avoir de meilleurs résultats, il serait nécessaire d'utiliser d'autres fréquences cycliques. Par exemple pour $\lambda=1$ et $K$ fréquences cycliques générées en ajoutant des termes du type $b_1\left(n\right) \sin\left(2\pi f_k n\right)$ au modèle, l'écart type du bruit sera réduit d'un facteur $\sqrt K$.

Bien qu'il paraisse intéressant d'avoir un terme $\lambda$ grand, il n'en est rien. En effet, le motif situé à la fréquence cyclique $0.22\;Hz$ sera reconstruit à partir des deux motifs. A l'aide de calcul similaire on obtient :

$$\esp{\varepsilon_C^2}=\sigma_B^2\left(\frac{\lambda^2}{2}\right)$$ (5.56)

Dans ce cas les résultats sont mauvais pour $\lambda$ grand. Par ailleurs, comme il s'agit de la contribution la plus énergétique l'erreur est nettement visible dans le domaine temporel. Une bonne reconstruction de toutes les composantes nécessite que chacune d'entre elles ait la même puissance (partie ``utile'' de $c\left(n\right)$ blanche). Il ne faut pas perdre de vue qu'un blanchiment préalable du signal $x\left(n\right)$ ne résoud pas ce problème car il intervient également sur le bruit $b_2\left(n\right)$.

La figure montre l'évolution de l'erreur pour la bande fréquentielle associée au motif central et au motif de droite. Cette figure confirme nos calculs (tout particulièrement pour le motif de droite) et les valeurs de $\lambda>1$. On peut remarquer que les courbes théoriques se croisent à $\lambda=1$.

Figure: Influence de l'allure du spectre du signal sur la reconstruction

Nous allons maintenant exploiter le modèle de signaux de roulement.

5.8 Signaux de roulement

Les signaux sont basés sur l'équation simplifiée que nous rappelons ici :

$$c\left(n\right) = r\left(t\right) \ast \left[ A\left(t\right) \sum_i \dirac{t}{\text{\tiny$T_i$}} \right]$$ (5.57)

$$x\left(n\right) = x\left(n\right) + b\left(n\right)$$ (5.58)

Ce modèle est très intéressant pour tester l'algorithme face aux fluctuations de vitesses.

La figure montre le spectre du filtre lors de l'utilisation d'une seule fréquence cyclique. Le spectre du filtre a une forme proche du spectre du signal (g). La figure permet d'expliquer cette ressemblance. Elle montre le spectre signal ainsi que sa version décalée. Le spectre est composé de pics régulièrement espacés (signal proche d'un signal périodique) avec une énergie légèrement distribuée autour des pics à cause des fluctuations de vitesse. Cette figure montre uniquement les premiers points du spectre afin de mettre en évidence les pics. A une échelle plus faible, ces pics seront difficilement discernables et le spectre mettra plutôt en évidence les résonances. Le décalage place ``face à face'' chacun des ces pics, dès lors le module du filtre correspondra à la différence (sur une échelle logarithmique) entre le module (en dB) des deux harmoniques.

Figure: Filtre pour un seul décalage

La figure montre les résultats pour $7$ fréquences (fréquence cyclique et ces harmoniques). Comme le montre la figure g, La puissance du bruit est inférieure de $4\;dB$ à la puissance maximale de $c\left(n\right)$. Un tel bruit rend la détection des pics difficile sur la figure a. La figure c montre l'estimation de $c\left(n\right)$. Le bruit entre les pics ayant été fortement réduit, il est possible de localiser clairement les défauts. Cette méthode apparaît assez intéressante pour le débruitage. La forte atténuation du bruit est confirmée par la figure h : le débruitage fait nettement apparaitre la résonance.

Figure: Signal - 7 fréquences - bruit important

5.8.1 Signaux réels

5.8.1.1 Banc d'essai d'engrenage

Les signaux proviennent d'un banc d'essai d'engrenage de l'UNSW présentant un défaut sur la bague externe. Ce défaut se traduit par une succession de chocs. Un signal accéléromètrique ainsi qu'un signal tachymétrique (codeur optique) ont été acquis pendant une durée correspondant à $45$ chocs. Le banc d'essai comporte un engrenage, il est nécessaire de supprimer la contribution issue des roues. Pour cela le signal a été re-échantillonné angulairement a posteriori à l'aide du signal tachymétrique afin d'estimer puis de retrancher la moyenne synchrone. La figure a montre le signal $x\left(n\right)$ après suppression de la moyenne synchrone sur un intervalle de $3$ chocs. Il est facile de distinguer le dernier choc, les deux premiers sont noyés dans le bruit. La variance synchrone (figure b) met en évidence ces chocs mais elle est assez bruitée. Après avoir utilisé l'algorithme de débruitage, on distingue nettement les trois chocs sur le signal (figure c) ainsi que sur la variance synchrone (figure d). En comparant les spectres avant et après débruitage (figure h), on peut remarquer qu'il y a très peu d'atténuation dans la bande $\left[0.18,0.27\right]\;Hz$. Cela signifie que dans cette bande, la puissance du bruit est faible devant la puissance du signal utile (défaut). Il est apparu intéressant de l'isoler par filtrage passe bande comme lors de l'analyse d'enveloppe (figure ). Cette figure confirme que c'est dans cette bande de fréquences que l'on trouve la contribution issue du défaut de roulement.

Figure: Débruitage de signaux de roulement

Figure: Signaux filtrés passe bande

5.8.1.2 Hélicoptère

Dans ce paragraphe nous montrons l'apport de la méthode de débruitage à l'ordre $2$ sur les signaux d'hélicoptère (les mêmes que pour le chapitre sur l'$ESANC$). Pour des raisons matérielles nous nous sommes limités à 1 millions de points ($0.1\;s$). A titre indicatif le calcul a nécessité plus de $700\;Mo$ de mémoire (en plus du système d'exploitation) et $1h30$ de calculs sur un Athlon XP 2400+ ($1\;Go$ de mémoire). Nous avons utilisé les $9$ harmoniques du défaut pour la reconstruction. La figure montre les résultats du débruitage à l'ordre $1$ seul ($ESANC$) et à l'ordre $1$ et $2$ ($ESANC$+débruitage). Les fréquences associées au défaut apparaissent nettement.

Figure: Transformée de Fourier de l'enveloppe au carré du signal débruité

5.9 Bilan et perspectives

Dans ce chapitre nous avons montré comment utiliser la cyclostationnarité pour le débruitage. D'abord la cyclostationnarité à l'ordre 1 en exploitant la périodicité du signal : on recherche un filtre permettant de prédire le signal à l'aide des valeurs précédentes ($SANC$). Compte tenu du caractère prédictible des signaux périodiques, ce filtre extrait les signaux (poly)cyclostationnaires à l'ordre 1. Malheureusement, les méthodes de type $SANC$, bien que supposant le signal périodique utilisent un signal cyclostationnaire flou. Dès lors, nous avons proposé de re-échantillonner les signaux préalablement ($E-SANC$). Ces méthodes de séparation ont été illustrées à l'aide de signaux d'hélicoptère contenant une composante cyclostationnaire à l'ordre 1 suivant l'angle associée aux roues ainsi qu'une composante cyclostationnaire floue associée aux roulements. Alors que la SANC était incapable de séparer les signaux de roulement à cause des fluctuations de vitesses trop importantes, la E-SANC a permis d'une part de réaliser cette séparation et d'autre part de mettre en évidence la fréquence associée au défaut.

Après séparation à l'ordre 1 de la composante périodique, le reste du signal contient à la fois le signal de roulement et du bruit. Nous avons alors recherché une méthode permettant de réduire ce bruit. Cette méthode, généralisant le filtrage de Wiener, exploite la redondance spectrale engendrée par la cyclostationnarité. Son but est de reconstruire le signal à l'aide d'une somme de contributions obtenues en décalant fréquentiellement puis en filtrant ce même signal. Les exemples ont montré que la réduction de bruit fonctionnait mieux avec un nombre important de décalages (6, 7 ou plus selon les capacités de l'ordinateur) ainsi que pour un signal ``blanc''. Cette technique a permis d'améliorer notablement les signaux et de faire apparaître nettement le défaut de roulement sur le signal d'hélicoptère.

Bien que ces méthodes soient prometteuses, il serait intéressant de disposer d'un ensemble de mesures sur un hélicoptère pour différents niveaux de dégradation afin d'étudier sa sensibilité face à l'évolution du défaut.