Renée El

Melhem Bât

Pascal Blaise

M Stéphane

ÉLECTROTECHNIQUE E E A Électronique

M Philippe Delachartre
email: philippe.delachartre@insa-lyon.fr

M Philippe
email: philippe.normand@univ-lyon1.fr

Mme Sylvie

M Hamamache Kheddouci

M Stéphane Benayoun
email: stephane.benayoun@ec-lyon.fr

Stéphanie Cauvin

M Jocelyn Bonjour

M Christian Montes

Keywords: Intelligence Maintenance, Rolling element bearing, Integrated auto-framework, Fault detection, Fault type identification, Fault signal reconstruction, Fault size characterization maintenance du renseignement, Roulement à élément roulant, Auto-framework intégré, Détection de fautes, Identification du type de défaut, Reconstruction du signal de défaut, Caractérisation de

Today, the most fundamental issue of condition monitoring in most industrial plants is fault diagnostics and prognostics. One of the most effective approaches to investigate this issue is condition monitoring based on vibration signal analysis. With the development of industry, multi-threaded maintenance and multi-channel acquisition are becoming more widespread in the current, which put forward higher requirements for maintenance. Based on this observation, it is proposed in this thesis one automated diagnosis framework for the rolling element bearing that integrates the successive steps of fault detection, fault type identification, fault signal reconstruction and fault size characterization. The advantage is that the complete diagnosis process is completed at once, while involving only one key hyperparameter, which improves the degree of automation of current Condition Based Maintenance (CBM) and liberating human participation.

In the presence of incipient fault, vibrations of rolling element bearings show symptomatic signatures in the form of repetitive impulses. This can be seen as a non-stationary signal whose statistical properties switch between two states. The proposed maintenance strategy models such characteristics with an explicit-duration hidden Markov model (EDHMM) and uses the estimated model parameters to perform integrated diagnosis without requiring the user's expertise. The detection of a fault is first achieved by means of a likelihood ratio test built on the EDHMM parameters. One statistical counting approach and posterior probability spectrum are then used for identifying the fault type automatically. In order to obtain the fault signal in some cases, one Bayesian filter based on the EDHMM parameters is constructed. Finally, the fault size is estimated from the duration times returned by EDHMM.

Subsequently, the capability of the integrated auto-diagnosis framework is illustrated on different experimental datasets. The first validation is forced on the vibration data iii Cette thèse est accessible à l'adresse : https://theses.insa-lyon.fr/publication/2022LYSEI045/these.

First of all, I would like to thank my supervisor, Prof. Jérôme Antoni, for his invaluable guidance, patience and support in all aspects of my research. His enthusiasm for research, his creative and critical way of thinking greatly affected me during my Ph.D pursuit. I am also grateful for the excellent example he has personified as a mentor and professor. I would also like to thank my CST jury members, Prof. Didier Remond and Prof.

Quentin Leclere, for their review on my research and encouraging and helpful suggestions, which really improved the quality of my Ph.D research. I would like to thank all my current and former colleagues of the LVA for their companionship and support, making my time in the office memorable and entertaining.

Especially thanks for Dr. Alice Dinsenmeyer, Dr. Fulbert Mbailassem, Dr. Xiaowen Zhu, for their kind support and help that have been of great value in my Ph.D period. A big thank for my bro Dr. Xiaokang Zhang who is not only a good colleague to work together, but also a bosom buddy to share the life. It has been a pleasure knowing you all and creating happy moments with you.

Finally my gratitude goes to my family members, relatives for their constant supports and encouragements throughout the Ph.D period. Especially to my fiancée Junli who provided emotional understanding and unconditional support during these challenging i

Résumé

Aujourd'hui, le problème le plus fondamental de la surveillance d'état dans la plupart des installations industrielles est le diagnostic et le pronostic des défauts. L'une des approches les plus efficaces pour étudier ce problème est la surveillance de l'état basée sur l'analyse des signaux de vibration. Avec le développement de l'industrie, la maintenance multithread et l'acquisition multicanal se généralisent, ce qui met en avant des exigences de maintenance plus élevées. Sur la base de cette observation, il est proposé dans cette thèse un cadre de diagnostic automatisé pour le roulement qui intègre les étapes successives de détection de défaut, d'identification du type de défaut, de reconstruction du signal de défaut et de caractérisation de la taille du défaut. L'avantage est que le processus de diagnostic complet est réalisé en une seule fois, tout en impliquant un seul hyperparamètre clé, ce qui améliore le degré d'automatisation de la maintenance basée sur les conditions (CBM) actuelle et libère la participation humaine.

En présence de défaut naissant, les vibrations des roulements présentent des signatures symptomatiques sous forme d'impulsions répétitives. Cela peut être vu comme un signal non stationnaire dont les propriétés statistiques basculent entre deux états. La stratégie de maintenance proposée modélise ces caractéristiques avec un modèle de Markov caché à durée explicite (EDHMM) et utilise les paramètres estimés du modèle pour effectuer un diagnostic intégré sans nécessiter l'expertise de l'utilisateur. La détection d'un défaut est d'abord réalisée au moyen d'un test de rapport de vraisemblance construit sur les paramètres de l'EDHMM. Une approche de comptage statistique et une probabilité sont ensuite utilisées pour identifier automatiquement le type de défaut. Afin d'obtenir le signal de défaut dans certains cas, un filtre bayésien basé sur les paramètres EDHMM est construit. Enfin, la taille du défaut est estimée à partir des durées renvoyés par EDHMM.

Par la suite, la capacité du cadre d'autodiagnostic intégré est illustrée sur différents v

Résumé étendu en français Introduction Contexte et motivations

En tant que composants essentiels, les roulements sont utilisés dans toutes sortes de systèmes mécaniques rotatifs, par exemple les éoliennes, l'aérospatiale, le transport, etc. Ce sont des composants d'entraînement très sollicités qui, avec le temps, deviennent sensibles à l'usure et à l'écaillage. Les performances normales des machines tournantes dépendent entièrement de l'état de santé des roulements, qui représentent près de 45 à 55% de ces défaillances d'équipement. La défaillance d'un seul roulement peut souvent provoquer l'arrêt de systèmes de fabrication entiers, réduisant ainsi la fiabilité et la disponibilité du système. Cela augmente à son tour les temps d'arrêt de la production causant une perte financière massive à l'organisation et peut même s'avérer dangereux pour la sécurité des travailleurs. Par conséquent, les stratégies de maintenance efficaces pour minimiser l'impact des défaillances des roulements sont très importantes pour l'industrie.

La stratégie de maintenance actuelle est la maintenance prédictive, qui peut être réalisée grâce à une série de tests et d'analyses sans qu'aucun dommage ne se produise. Elle est organisée et mis en oeuvre en fonction de l'état de fonctionnement de l'équipement et est appelée Maintenance Based Maintenance (CBM). Cependant, ce type de maintenance repose encore principalement sur l'expertise humaine, qui sera effectuée à un coût élevé et pour une quantité limitée de données seulement. De plus, la surveillance de l'état du roulement n'est pas simple, ni la détection ou l'identification du type de défaut lorsqu'un défaut s'est produit. En réalité, le diagnostic des roulements est un sujet important et de nombreux sous-problèmes doivent être résolus. Cependant, les techniques actuelles sont presque toutes conçues pour des problèmes spécifiques, plutôt que pour l'ensemble vii du diagnostic, ce qui est très exigeant pour l'utilisateur et nécessite des connaissances en algorithmes et une expérience d'expert. Tels sont les obstacles sur la route de la maintenance.

Objectifs

Sur la base de l'importance et des défis mentionnés ci-dessus de la surveillance de l'état des roulements, l'objectif de cette thèse est de présenter un cadre d'autodiagnostic intégré pour les roulements. L'importance et la contribution de cette recherche sont le développement d'un cadre d'autodiagnostic intégré grâce à la modélisation stochastique du signal pour réaliser à la fois la détection de défaut, la reconstruction du signal de défaut, l'identification du type de défaut et la caractérisation de la taille du défaut. Pour atteindre ces objectifs, un modèle stochastique, le modèle de Markov caché à durée explicite (EDHMM), est sélectionné pour modéliser dans un premier temps la distribution tempsfréquence du signal de vibration. Les paramètres estimés du modèle EDHMM sont utilisés pour effectuer ces tâches sans avoir besoin d'autres informations préalables ni de l'expertise de l'utilisateur.

Les principales contributions du présent travail sont résumées comme suit. Il introduit une méthodologie complète qui:

• réalise la détection, l'identification et la caractérisation des défauts à la fois, dans un cadre intégré qui dépend d'un seul hyperparamètre clé ; cela vient avec une test du rapport de vraisemblance pour la détection d'un défaut un filtre bayésien pour reconstruire le signal de défaut un spectre de probabilité a posteriori qui est une alternative robuste au spectre d'enveloppe standard, exempt de tout pré-traitement, et une autre technique d'identification automatique d'un point de vue statistique sans aucun examen visuel une technique simple basée sur la régression linéaire pour estimer la taille du défaut • est non supervisé, dans le sens où elle s'applique sans besoin de données historiques et sans apprendissage, viii • est automatisée, en ce sens que ces tâches sont effectuées consécutivement sans nécessiter l'intervention de l'utilisateur.

Chapitre 2

Ce chapitre décrit les caractéristiques de base des vibrations des roulements et l'état de l'art sur le diagnostic des défauts des roulements. Premièrement, les caractéristiques des roulements sont décrites dans la section 2.2, y compris la structure de base, les caractéristiques de vibration et les caractéristiques de défaut naissant. La section 2.3 décrit quelles techniques actuelles sont appliquées pour le diagnostic des défauts de roulements de différents points de vue, y compris la détection, l'amélioration du signal de défaut, l'identification et les techniques d'estimation de la taille des défauts. La littérature sur les méthodes de modélisation stochastique est passée en revue, en particulier pour le HMM et l'EDHMM. Enfin, les techniques de diagnostic automatique sont également passées en revue. Sur la base de cette recherche bibliographique, quelques points sont dignes de mention. 1). À travers la littérature disponible, comme indiqué dans la section 2.3, il est clair que diverses techniques avancées ont été proposées pour répondre aux différentes tâches de diagnostic. Cependant, ces techniques sont presque toutes indépendantes et nécessitent que les utilisateurs aient un niveau élevé d'expérience et de connaissances lorsqu'ils traitent différents problèmes de diagnostic. 2). En ce qui concerne les techniques qui tentent d'aller dans le sens d'un diagnostic intégré et automatisé dans la section 2.5, il n'existe encore aucune recherche capable d'intégrer tous les sous-problèmes de diagnostic dans une seule solution. Encore moins de travaux ont essayé d'atteindre ces objectifs de manière automatisée, c'est-à-dire sans réglage manuel des algorithmes par un expert. Donc toutes ces raisons nous poussent dans une direction, l'autodiagnostic intégré, ce qui est exactement ce que poursuit la thèse.

Chapitre 3

Ce chapitre examine le modèle stochastique et comment l'utiliser pour modéliser la distribution temps-fréquence du signal de vibration du roulement, y compris le paramétrage, l'estimation et l'analyse. La section 3.2 introduit d'abord la chaîne de Markov, ouvrant la voie aux HMM et EDHMM suivants respectivement dans les sections 3.3 et 3.4. Le cadre ix proposé dépend principalement de l'EDHMM. Il vise à capturer l'évolution temporelle des coefficients STFT et à calculer la durée dans différents états. L'estimation de ces paramètres nécessite des algorithmes spécifiques, qui sont détaillés dans la section 3.5. Certains paramètres clés de l'EDHMM sont également abordés dans la section 3.5. L'efficacité de ces paramètres et les caractéristiques qu'ils révèlent sont illustrées par un signal synthétique de défaut de roulement. Enfin, l'hyperparamètre de ce modèle et l'initialisation des paramètres sont discutés. Il s'avère qu'un seul hyperparamètre, c'est-à-dire la longueur de la fenêtre N w , doit être défini à l'avance. Cela rend l'approche potentiellement plus robuste et mieux adaptée aux applications pratiques en l'absence de connaissances préalables.

Chapitre 4

Ce chapitre présente le cadre de diagnostic intégré, y compris la détection des défauts, la reconstruction du signal de défaut, l'identification du type de défaut et la quantification de la taille du défaut. Ce chapitre est le travail de base de cette thèse. Les paramètres EDHMM obtenus dans le chapitre précedent sont utilisés dans le cadre de diagnostic proposé. La section 4.2 introduit le test du rapport de vraisemblance (LRT) et prête les paramètres EDHMM pour détecter la présence nécessitant un pré-traitement des données. Si un défaut est détecté, l'étape d'identification renvoie le spectre de probabilité postérieur (PPS), l'équivalent d'un spectre d'enveloppe, qui peut être analysé visuellement ou automatiquement. Afin d'automatiser entièrement, une méthode simple est proposée basée sur un point de vue statistique pour calculer la probabilité postérieure de différents types de défauts. De même, le signal de défaut peut être extrait grâce à une filtre de variable en temps construit par la matrice de covariance de l'observation. Indépendamment de ce dernier, l'étape de caractérisation du défaut utilise le paramètre de Poisson de l'EDHMM pour évaluer la taille du défaut par régression linéaire, offrant ainsi une solution simple à une tâche difficile. Il est à nouveau souligné que l'algorithme complet ne repose que sur l'hyperparamètre critique, N w , l'inverse de la résolution fréquentielle. En particulier, il ne s'appuie ni sur des données historiques, ni sur des informations préalables, ni sur une intervention manuelle.

x

Chapitre 5

Ce chapitre présente la validation du cadre de diagnostic intégré à travers différents jeux de données. La contribution de cette partie est de déterminer les capacités du cadre d'autodiagnostic intégré proposé pour les roulements. Deux scénarios expérimentaux différents sont présentés dans cet article. La première validation expérimentale est relative à des roulements avec différents dommages, fonctionnant à différentes vitesses et sous différentes charges. Les résultats des 36 ensembles de données de roulement associés à trois cas de défaillance différents montrent l'efficacité de ce cadre. Le deuxième type de données expérimentales rapporte le comportement d'un seul roulement endommagé soumis à un long test à vitesse et charge constantes. Lors de la validation, la méthode proposée a également été comparée à l'approche de pointe, le kurtogram. Pour la première validation, quelques commentaires s'imposent 1) Il s'avère que la technique de détection est efficace pour détecter le temps d'occurrence des défauts, comme le montre la figure 5.13(c) dans le deuxième cas associé aux données de dégradation accélérée. 2) La fréquence caractéristique du roulement a été facilement modulée par la fréquence de l'arbre produisant une série de bandes latérales et d'harmoniques, ce qui a une grande influence sur l'identification.

Cependant, les deux techniques d'identification proposées surmontent bien ce problème, vu la Fig. 5.3, la Fig. 5.7 et la Fig. 5. [START_REF] Bulent Ayhan | Integrated Bearing Diagnostics and Fault Severity Estimation Using HMM[END_REF]. 3) Elles sont également capables de quantifier la taille de défaut définé dans la section 5.2, et de fournir des informations sur la propagation des fissures et suivre quantitativement l'évolution de la fissure vue à la section 5.3.

Conclusion

L'objectif principal de cette thèse est l'étude de méthodes de traitement automatique du signal pour différentes tâches de diagnostic dans un cadre intégré. Les avantages remarquables sont qu'elles fonctionnent sans données historiques, sans beaucoup d'hyperparamètres, même sans aucune intervention manuelle dans le processus. La revue

Background and Motivation

As crucial components, rolling bearings serve almost in all kinds of rotating mechanical systems, e.g., wind turbine, aerospace, transportation etc. They are highly-stressed drive components which, over a period of time, become susceptible to wear and spall. The normal performance of rotating machinery is entirely dependent upon the health state of the rolling bearings, which accounts for almost 45-55% of these equipment failures [START_REF] Benali | Importance of the fourth and fifth intrinsic mode functions for bearing fault diagnosis[END_REF][START_REF] Nandi | Condition monitoring and fault diagnosis of electrical motors-A review[END_REF]. The failure of just a single bearing can often cause the stoppage of entire manufacturing systems thereby reducing the reliability and availability of the system. This in turn increases the production downtime causing a massive financial loss and may even prove dangerous to the safety of the workers. Therefore the effective maintenance strategy to minimize the negative impact of bearing failures is very important for the industry.

Fortunately, the maintenance theory is getting more and more advanced to meet the existing large amount of industrial monitoring data. The development of maintenance theory has gone through several important stages in the past decades, from Run-to-Failure maintenance, Scheduled maintenance, Predictive maintenance to Intelligence maintenance.

The Run-to-Failure maintenance is the oldest definition of maintenance, which means that the maintenance only occurs after a sudden or even catastrophic failure. Therefore, such maintenance will cause a lot of time and money costs and even life-threatening. The Based Maintenance (CBM). However, human expertise is an outstanding solution but at a high cost and for a limited quantity of data only, the analysis being time-consuming.

Recently, with the advent of Industry 4.0, the digital factory offers many alternatives to human monitoring. The Intelligence Maintenance is a promising direction, which is dedicated to improving the degree of automation of maintenance and liberating human participation. Figure 1.1 shows the main steps in the intelligence maintenance system.

In the current maintenance strategy, one of the most effective approaches is condition 2/129 Therefore recognizing that kind of transient from non-stationary signals generated from bearing needs helps from the specific signal processing tools.

In addition, as shown in Figure 1.1, in the condition monitoring of bearings, it is not simple to detect or identify the fault type when it occurs. In reality, bearing diagnosis is a big topic, and there are many sub-issues that need to be solved. However, current techniques are almost designed for specific issues, rather than toward the whole diagnosis, which is high demanding for the user and requires different algorithms knowledge and expert experience. These are the obstacles on the road to intelligence maintenance. Based on the above mentioned importance and challenges of bearing condition monitoring, the 3/129 motivation is to investigate further the use of vibration analysis in solving efficiently some key problems in condition monitoring, so as to present an integrated auto-diagnostic framework for the rolling element bearing.

Objectives and scope of the research

According to the identified research problem and the above-mentioned motivations, during the envisage phase of diagnosis solution, it is necessary to consider different factors, i.e., simplicity, practicality, automatically. The following questions can be formulated as input to this thesis.

• How to achieve integrated diagnosis?

Is it possible to integrate all the diagnosis issues (detection, reconstruction, identification, characterization) into one-package framework? This means completing all the diagnosis issues at once without manual tuning of the algorithms.

• How to achieve unsupervised diagnosis?

Is it possible to employ an advanced signal processing technique and stochastic modeling which requires only one measurement and does not need to be trained by a sequence of data? This means without need of historical data and without training.

• How to achieve automated diagnosis?

Would it be possible to apply such a technique automatically? This means these tasks are completed consecutively without requiring the user's intervention, and as few hyperparameter settings as possible.

Based on the concerned questions above, this research focuses on the development • achieves detection, identification, and characterization of faults at once, in an integrated framework that depends on only one key hyperparameter; this comes with a likelihood ratio test for detection of a fault [START_REF] Antoni | A statistical methodology for the design of condition indicators[END_REF] a Bayesian filter for reconstruction of the fault signal [START_REF] Xin | Semi-automated diagnosis of bearing faults based on a hidden Markov model of the vibration signals[END_REF] 4/129 • is automated, in the sense that these tasks are completed consecutively without requiring the user's intervention.

Fault Detection

Fault signal Reconstruction

Fault type Identification

Fault size Characterization

Bearing CBM The estimation of its parameters requires specific algorithms, which are detailed in Section 3.5. In order to have better performance, the model needs some parameters to be set in advance, which is discussed in Section 3. proposes a posterior probability spectrum (PPS) of the states; it is used as a robust alternative to the state-of-the-art envelope spectrum for identifying the fault frequencies.

Organization of the thesis

In order to avoid visual examination, a simple method based on a statistical point of view for identifying automatically the fault type is proposed. Finally, the fault size quantification through a simple idea of linear regression is described in Section 4.6. In this section, the time duration parameter of the EDHMM is used to assess the fault size, thus offering a simple solution to a challenging task.

Introduction

This chapter provides an overview of bearing and the current vibration monitoring approaches. The chapter begins with the introduction of various characteristics of rolling element bearings. Since bearing vibration monitoring consists of predefined steps, several state-of-the-art approaches for each diagnostic module are discussed subsequently, including fault detection, fault signal enhancement, fault identification and fault size characterization.

Markov process is a desired tool that models the nonstationary nature of the vibration signal and uncovers the intrinsic structure in signals. Hence, in this chapter, a review of Hidden Markov model and its variant is also given. Vibration-based monitoring of rolling element bearings is a complex topic, which involves several issues, then the integrated auto-diagnosis techniques are investigated based on current literature. Furthermore, a preliminary conclusion concerning the limitation of current research and the focus of this thesis is drawn in the last subsection.

Rolling element bearing signal characteristics

As a special structure, a rolling element bearing has its own unique characteristics, including vibration characteristics and incipient fault characteristics. These characteristics enable bearing to generate useful vibration information when rotating, which offers the possibility of diagnosing potential bearing faults. This section considers some basic but important knowledge, containing the structure of rolling element bearings, the common bearing fault types and how they occur, also the mechanism of vibration generation in rolling element bearings.

Basic structure of rolling element bearing

Most rolling bearings consist of rings with a raceway (inner ring and outer ring), rolling elements and a cage as shown in Figure 2.1. In reality, the outer ring is often fixed on the bearing seat or mechanical body. The inner ring is connected to the drive shaft and driven by the shaft rotation. Rolling elements geometrically contact with the raceway surfaces of the inner and outer rings at "points", whose role is to convert motion from sliding to rolling. The cage separates the rolling elements at regular intervals, holds them in place within the inner and outer raceways, and allows them to rotate freely and reduces

Vibration characteristics of rolling element bearings

Vibration monitoring has now become a well accepted part of many planned maintenance regimes and relies on the well known characteristic vibration signatures which rolling bearings exhibit as the rolling surfaces degrade. However, in fact, vibrations occur all the time. As a rolling element passes over a location on a race, the surface of the race is subjected to a high effective stress for a short time during element passage. This cycling of stress occurs once per element passage at the location. Simultaneously, a location on the rolling element itself is also subjected to a cycle of subsurface stress twice per revolution of the element. The vibrations are generated when surfaces interact through a combination of rolling and sliding. The bearing natural vibration is originated from elastic factors; so whether the bearing fails or not, it will produce vibration.

When the bearing experiences a localised defect during operation, each interaction will cause an instantaneous impact and excite the system to perform high-frequency free attenuation vibration response according to its natural frequency. Such impact has higher energy than the natural excitation. Due to the bearing geometrical symmetry and dynamic regularity, the impact will show typical signal signature.

For a stationary outer ring and rotating inner ring, from the bearing geometry the 11/129 fundamental frequencies are derived as follows [START_REF] Randall | Rolling element bearing diagnostics-A tutorial[END_REF]:

Bearing characteristic frequencies:

• Ballpass frequency, outer race:

f BP F O = Zfr 2 [1 -d D cosα]
• Ballpass frequency, inner race:

f BP F I = Zfr 2 [1 + d D cosα] • Fundamental train frequency (cage speed): f F T F = fr 2 [1 -d D cosα]
• Ball (roller) spin frequency:

f BSF = Dfr 2d [1 -( d D cosα) 2 ]
where, f r is the shaft speed, Z is the number of rolling elements, and α is the angle of the load from the radial plane. These bearing equations assume that there is no sliding and that the rolling elements roll over the raceway surfaces. However, in practice this is rarely the case and due to a number of factors the rolling elements undergo a combination of rolling and sliding. As a consequence the actual characteristic defect frequencies may differ slightly from those predicted. Therefore, it is necessary to take the nearby frequency value as the fault frequency according to the actual situation.

There are some unique characteristics for the faulty bearing vibrations, such as impulsiveness [START_REF] Antoni | The infogram: Entropic evidence of the signature of repetitive transients[END_REF], cyclostationarity [START_REF] Antoni | The infogram: Entropic evidence of the signature of repetitive transients[END_REF] and modulation [START_REF] Randall | Rolling element bearing diagnostics-A tutorial[END_REF].

• Impulsive characteristic

The impulsive waveform in the faulty bearing vibration signal is due to the response of the structure to excitation forces whenever a faulty point passes through the load zone. Such transient pulse is followed by decaying oscillations with natural frequency corresponding to structural resonances. However, the impulsive waveform caused by the fault is easily submerged in strong background noise at the incipient fault stage, resulting in a low signal-to-noise ratio, which brings great difficulty to the feature extraction of rolling bearings.

• Cyclostationary characteristic

Due to the bearing vibration generation form, i.e. combination of rolling and sliding, the vibration signals from rolling-element bearings often exhibit high levels of 12/129 cyclostationarity, especially in the presence of localised faults [START_REF] Antoni | Cyclostationarity by examples[END_REF]. Traditional signal processing methods based on Fourier transform encounter certain limitations or are powerless in some cases. This asks for the development of specific signal processing tools [START_REF] Antoni | Cyclic spectral analysis of rolling-element bearing signals: Facts and fictions[END_REF].

•

Incipient fault characteristics of rolling element bearing

The rolling-contact bearing is an element of machinery with a very important role, and it dominates the performance of the machine. However, a rolling element bearing will eventually fail through a fatigue mechanism that can be greatly accelerated by poor installation, overloading, improper lubrication, or contamination. If the abnormality can be detected timely in the early incipient stage, not only the maintenance cost will be reduced, but even catastrophic accidents can be effectively avoided. Therefore, the early detection of bearing faults is extremely important. When the bearing has an initial failure during operation, the vibration signal will have two important characteristics [START_REF] Antoni | A statistical methodology for the design of condition indicators[END_REF]; the first one can be described as an increase in the impulsiveness of the signal, i.e. a deviation from Gaussianity, and/or the second one as a shift of the statistical behaviour of the signal from stationarity to non-stationarity. However, in the incipient fault stage, the level of the vibration will often show unobvious signal features, weak impulsive energy and low signal to noise ratio, etc. These characteristics depend upon many factors, summarised by the following points,

• in the incipient stage, the impact energy of the small fault will be weak, which makes the impact easily submerged in the noise;

• the point at which the vibration is measured may correspond to complicated vibration propagation path, and noise increases during the process of propagation;

13/129

• the operating environment of the bearing and other rotating components, i.e. electric motors, gears, shaft, etc., may modify the bearing features;

• the acquisition and conversion error of the sensor is another non-negligible factor.

Therefore, in order to perform better maintenance, one good measurement environment is very important. In addition, the recovery of the fault signature in a low signal-to-noise ratio and the link to the stochastic and nonstationary nature of the fault signature are the critical problems in the incipient fault diagnosis.

General review of diagnosis of rolling element bearings 2.3.1 Bearing fault detection

During operation, the bearings are subjected to heavy and dynamic loadings generated by machines and transmitted through the components of rolling element bearings. Hence, the healthy condition of bearings is very important in the rotating machinery system. Any failure in the bearings must be detected on time to avoid increase in downtime, production time and catastrophic failure of the machinery. Thus, the detection of these defects is of vital importance for condition monitoring and quality inspection of bearings. As described in Chapter 2.2.2, an instantaneous impact will appear during each interaction. These pulses produce vibration which are the main signatures for proving the fault existence. In this

SK = L t=1 (x t -x) 3 (L -1)σ 3 , (2.1)
where, x t , t = 1, 2, ..., L is the tth sampling point of the raw signal x, L is the number of sampling points, x is the mean value of the sample x defined as

x = 1 L L t=1 x t ,
and σ is the standard deviation of the sample x defined as

σ = 1 L -1 L t=1 (x t -x) 2 .

Kurtosis (KU) [53]

KU = L t=1 (x t -x) 4 (L -1)σ 4 ;
(2.2)

Crest indicator (CI) [START_REF] Dron | Improvement of the sensitivity of the scalar indicators (crest factor, kurtosis) using a de-noising method by spectral subtraction: application to the detection of defects in ball bearings[END_REF]

CI = max|x t | 1 L L t=1 (x t ) 2 ; (2.3) Clearance indicator (CLI) [133] CLI = max|x t | ( 1 L L t=1 |x t |) 2 ; (2.4) Shape indicator (SI) SI = 1 L L t=1 (x t ) 2 1 L L t=1 |x t | ; (2.5)
Impulse indicator (IMI)

IM I = max|x t | 1 L L t=1 |x t | . (2.6) 15/129
The six dimensionless parameters are extracted from the raw vibration signals to reflect the general change of the healthy condition of bearings. Among these dimensionless indicators, the Kurtosis is more sensitive to impact signals and most widely used in the fault detection of rolling bearings. However, the above indicators can be considered as the impulsiveness measures. The cyclic behavior is very an important characteristic in the rotating machine vibration signal. Antoni proposed the indicator that used the negentropy of the square envelope spectrum to indicate the cyclostationary behavior [START_REF] Antoni | The infogram: Entropic evidence of the signature of repetitive transients[END_REF]. In addition, some indicators based statistical point of view are designed to identify the cyclostationary behavior in Ref. [START_REF] Antoni | A statistical methodology for the design of condition indicators[END_REF],

Indicator for testing the Gaussian cyclostationary hypothesis (GCS) against the Gaussian stationary (GS)

I GCS/GS (x) = ln < s 2 (t) > -< ln s 2 (t) >, ( 2.7)
where, s 2 (t) is N-periodic component of the squared envelope of the signal, defined as

s 2 (t) = 1 K K-1 k=0 |x(t + kN )| 2
with K = L/N is the number of N-sample periods of in the L-sample long signal. The operation < s 2 (t) > stands for the time average value of s 2 (t).

Indicator for testing the generalized Gaussian cyclostationary (GGCS) hypothesis against the generalized Gaussian stationary (GGS) is

I GGCS/GGS (x) = 2κ -1 0 ln < s κ 0 (t) > -2κ -1 1 < ln s κ 1 (t) > +C, (2.8)
where, s κ (t), κ = κ 0 , κ 1 , is N-periodic component of the κ-power envelope of the signal and C stands for an additive constant. s κ (t) is defined as

s κ (t) = 1 K K-1 k=0 |x(t + kN )| κ
with κ 0 , κ 1 are the estimates of the shape parameters of generalized Gaussian distribution.

Some researchers have also proposed band pass filter techniques based on time domain in which impulse loading excites structural resonance in the high frequency zone that can be detected by a transducer. The shock pulse method [START_REF] Alfredson | Time domain methods for monitoring the condition of rolling element bearings[END_REF] is based on band pass filtering 16/129

techniques. This method is widely used in industries but cannot effectively detect defects at low speed. Dwyer [START_REF] Van Dyer | Detection of Rolling Element Bearing Damage by Statistical Vibration Analysis[END_REF] initially introduced the concept of spectral kurtosis and Antoni [START_REF] Antoni | The spectral kurtosis: a useful tool for characterising non-stationary signals[END_REF][START_REF] Antoni | The spectral kurtosis: application to the vibratory surveillance and diagnostics of rotating machines[END_REF] implemented it in the field of fault diagnosis of rolling element bearings. Spectral Kurtosis (SK) is a statistical tool that identifies the presence of a series of transients and its position in the frequency domain. Antoni [START_REF] Antoni | Fast computation of the kurtogram for the detection of transient faults[END_REF] in his paper listed various properties of SK for the first time. More recently, Antoni [START_REF] Antoni | The infogram: Entropic evidence of the signature of repetitive transients[END_REF] summarized evidence studies of transients and proposed a new spectral extraction method, namely infogram, which considers both the two characteristics and uses negentropy of the square envelope and its spectrum to depict the two characteristics of transients.

Stochastic model

As described in Chapter 2.2.2, some signal changes will appear in the process of bearing failure; one is the increase in the impulsiveness of the signal, another is a shift of the statistical behaviour of the signal from stationarity to non-stationarity, i.e. the occurrence of repetitive transients. Although the traditional time-domain indicators are effective and sensitive to impulses, they have no ability to differentiate them. The probabilistic model is another common used way to detect bearing faults. In reference [START_REF] Antoni | A statistical methodology for the design of condition indicators[END_REF],

a statistical methodology based on the maximum likelihood ratio is introduced as a general framework to design condition indicators. The proposed indicator is optimal in the sense that it coincides with the statistics of the likelihood ratio test (LRT) which maximizes the probability of true detection given a fixed probability of false alarm (erroneous rejection of the null hypothesis). HMM is also a parametric statistical method that has the capability of pattern classification and is suitable for dynamic time series of signals that are nonstationary. The reference about bearing fault detection based on HMM can be seen in [START_REF] Ocak | HMM-based fault detection and diagnosis scheme for rolling element bearings[END_REF][START_REF] Ogbemudia | Fault detection in roller bearing operating at low speed and varying loads using Bayesian robust new hidden Markov model[END_REF]. In reference [START_REF] Zhang | Anomaly detection: A robust approach to detection of unanticipated faults[END_REF], the researchers introduced a fault detector based on the and integrate it into the defect detection architecture. A parametric adaptation algorithm was then introduced. The features and fault progression patterns are then defined in a particle filter framework where the current feature distribution is compared to its baseline counterpart to detect a gap or divergence between the two.

Bearing fault signal enhancement

Another significant and often unavoidable problem in signal processing is the presence The problem can be simplified as only estimating the spectral amplitude |X(f )|, which is based on an underlying fact that the short-time spectral amplitude rather than phase that is important for vibration signal. The spectral amplitude estimation algorithm is one kind of the classic approaches for de-noising, which can be classified as minimum mean square error (MMSE), spectral subtraction, wiener filter based on their estimation method.

The spectral subtraction estimation approach developed by Boll [START_REF] Boll | Suppression of acoustic noise in speech using spectral subtraction[END_REF] is designed for enhancing signals degraded by uncorrelated additive noise. It is an approach for estimating the magnitude frequency spectrum of the underlying clean signal by subtracting the noise magnitude spectrum from the noisy spectrum. The power spectrum subtraction [START_REF] Soo | Enhancement and bandwidth compression of noisy speech[END_REF] has the same idea, but instead use the power spectrum. These two basic methods are derived in Chapter 2.3.2.1. Some other variations on the basic spectral subtraction approach have been proposed. Notable work in this area is that of McAulay and Malpass [START_REF] Mcaulay | Speech enhancement using a soft-decision noise suppression filter[END_REF], who 18/129 formulated the spectral subtraction approach as a maximum likelihood estimation problem of the variance of each spectral component of the original clean signal. Other popular modifications are those that involve averaging or smoothing of the sample spectrum estimator, controlling the amount of subtracted noise [START_REF] Soo | Enhancement and bandwidth compression of noisy speech[END_REF][START_REF] Berouti | Enhancement of speech corrupted by acoustic noise[END_REF].

Wiener filter [START_REF] Levinson | The Wiener (root mean square) error criterion in filter design and prediction[END_REF] is applied either in the time domain or frequency domain to obtain an estimate of the undergraded signal. In the time domain the Wiener filtering operation can be represented as a convolution, and in the frequency domain as multiplication by the complex filter gain function. At a single frequency, the Wiener filter spectrum estimate [START_REF] Soo | Enhancement and bandwidth compression of noisy speech[END_REF][START_REF] Mcaulay | Speech enhancement using a soft-decision noise suppression filter[END_REF] is therefore a gain function times the corrupted spectrum Y (f ). The Wiener filter is also the optimal MMSE estimator under a Gaussian assumption.

Ephraim et al. [START_REF] Ephraim | Speech enhancement using a minimum-mean square error short-time spectral amplitude estimator[END_REF] proposed one algorithm assuming that the short-time spectrum follows the Gaussian distribution, and utilized MMSE on the short-time spectral amplitude (STSA) of the signal for enhancing the noisy signal. In the literature [START_REF] Ephraim | Speech enhancement using a minimum meansquare error log-spectral amplitude estimator[END_REF], Ephraim took the measure of mean-square error of the spectra into account, and extended the STSA estimator which minimizes the mean-square error of the log-spectra in enhancing noisy signal. These two STSA estimators were derived under the Gaussian assumption. Porter and Boll [START_REF] Porter | Optimal estimators for spectral restoration of noisy speech[END_REF] gave a modification and improvement for the optimal MMSE estimator calculated directly from noisy data in a way that avoids the need for assuming a specific form of the distribution.

Spectral Subtraction

Spectral subtraction [START_REF] Boll | Suppression of acoustic noise in speech using spectral subtraction[END_REF] is a simple and efficient method. The approach is to estimate the magnitude frequency spectrum of the underlying clean spectrum by subtracting the noise magnitude spectrum from the noisy spectrum. As the signal in practice tends to be non-stationary, then the STFT is applied here. The signal model can be rewritten in short time spectrum as

Y (n, k) = X(n, k) + N (n, k),
where n means the window index, and

Y (n, k), X(n, k), N (n, k) represent
∠N (n, k) of N (n, k) is replaced by the phase ∠Y (n, k) of Y (n, k
). These substitutions result in the spectral subtraction estimator X(n, k):

X(n, k) = [|Y (n, k)| -|N (n, k)|]e j∠Y (n,k) or X(n, k) = H 1 • Y (n, k),
where,

H 1 = 1 - |N (n, k)| |Y (n, k)| . (2.9)
In the literature [START_REF] Soo | Enhancement and bandwidth compression of noisy speech[END_REF][START_REF] Mcaulay | Speech enhancement using a soft-decision noise suppression filter[END_REF], one similar idea called power spectrum subtraction is presented. The author estimated the spectral amplitude firstly, and then used the phase of the original signal ∠Y (n, k) to form the estimate X(n, k). The square of the Fourier coefficient modulus can be written as

|Y (n, k)| 2 = |X(n, k)| 2 + |N (n, k)| 2 + X(n, k) † N (n, k) + X(n, k)N † (n, k). (2.10)
From the observed data y(t), |Y (n, k)| 2 can be obtained directly. The terms

|N (n, k)| 2 , X(n, k) † N (n, k) and X(n, k)N † (n, k
) cannot be obtained exactly. In the power spectrum subtraction technique, they are approximated by

E[|N (n, k)| 2 ], E[X(n, k) † N (n, k)] and E[X(n, k)N † (n, k)]. Because the additive noise N (n, k) is statistically independent with the signal of interest X(n, k), then the quantity E[X(n, k) † N (n, k)] and E[X(n, k)N † (n, k)]
are zero. The Eq.(2.10) can be rewritten as

| X(n, k)| 2 = |Y (n, k)| 2 -E[|N (n, k)| 2 ], (2.11)
where, the

E[|N (n, k)| 2 ] can be replaced with |N (n, k)| 2 if noise is stationary. The phase of X(n, k) is from ∠Y (n, k), so that, X(n, k) = | X(n, k)| • e j∠Y (n,k) .
As pointed in reference [START_REF] Soo | Enhancement and bandwidth compression of noisy speech[END_REF], the estimated quantity is not guaranteed to be non-20/129 negative since the right-hand side of Eq.(2.10) or Eq.(2.11) can become negative, and a number of somewhat arbitrary choices have been made. In many studies, the negative values are set to zero.

Wiener filtering

In the previous section, the basis for estimating the short-time spectral magnitude through a process of spectral subtraction was described. In this section, we review the wiener filter technique, in which a frequency weighting for an "optimum" filter is first estimated from the noisy observation. Here, Y (n, k), X(n, k), N (n, k) again denote the short-time spectra associated with the windowed time functions y(t), x(t) and n(t). The estimate X(n, k) of X(n, k) that minimizes the mean-square error is obtained by filtering y(n) with the noncausal Wiener filter as

X(n) = H 2 • Y (n, k),
where, the frequency weighting H 2 can be represented by the power quantity P y , P x and P n as

H 2 = P x P x + P n = E[|X(n, k)| 2 ] E[|X(n, k)| 2 ] + E[|N (n, k)| 2 ]
.

(2.12)

Based on the independence assumption and the relationship Eq.(2.10

), E[|Y (n, k)| 2 ] can be simplified as the sum of E[|X(n, k)| 2 ] and E[|N (n, k)| 2 ],
Then Eq.(2.12) can be also written as

H 2 = 1 - E[|N (n, k)| 2 ] E[|Y (n, k)| 2 ] . ( 2

.13)

The result Eq.(2.13) derived from the MMSE estimator can be also seen as the Wiener filter. Similar with the previous section, the spectral magnitude |N (n, k)| here can also be replaced by its average value |N (n, k)| taken during normal vibration (without fault). In addition, the quantities

E[|Y (n, k)| 2 ], E[|X(n, k)| 2 ], and E[|N (n, k)| 2 ]
can also be understood as the corresponding covariance matrix. Therefore, the covariance matrix will be used to denoise thanks to Eq.(2.13), which will be introduced detailed in the Chapter 4.3.

21/129

Minimum Mean-Square Error

In this subsection, one STSA estimation through minimum mean-square error is described, and the spectral components are assumed as statistically independent Gaussian random variables with zero mean. This approach was first proposed by Ephraim in [START_REF] Porter | Optimal estimators for spectral restoration of noisy speech[END_REF],

and later used in various fields. Based on the formulation of the estimation problem given in the previous section, our task is to estimate the modulus

|X(n, k)| from the observed signal y(t). The MMSE estimator | X(n, k)| of |X(n, k)| is obtained through the Bayesian theory as, | X(n, k)| = E[|X(n, k)||y(t)] = E[|X(n, k)||Y (n, k)] = X(n,k) |X(n, k)|p(X(n, k)|Y (n, k)) = X(n,k) |X(n, k)| p(Y (n, k)|X(n, k))p(X(n, k)) p(Y (n, k)) = X(n,k) |X(n, k)|p(Y (n, k)|X(n, k))p(X(n, k)) X(n,k) p(Y (n, k)|X(n, k))p(X(n, k)) . ( 2.14)
If one substitutes the Gaussian distribution into Eq.(2.14) and performs some manipulations (see literature [START_REF] Porter | Optimal estimators for spectral restoration of noisy speech[END_REF]), the estimate | X(n, k)| can be obtained as

| X(n, k)| = ξ(n, k) 1 + ξ(n, k) |Y (n, k)|,
where ξ(n, k) is the covariance ratio at frequency k of the nth spectral component of the clean signal and the noise ξ

(n, k) = E[|X(n,k)| 2 ] E[|N (n,k)| 2 ]
. Similarly, ∠Y (n, k) is used as the phase of X(n, k), and substituting ξ(n, k) into the above formula as

X(n, k) = E[|X(n, k)| 2 ] E[|X(n, k)| 2 ] + E[|N (n, k)| 2 ] • Y (n, k). (2.15)
It is found that the MMSE estimator of the magnitude of the kth signal spectral component is in fact same as the Wiener estimator. For this reason, Eq.(2.15) is referred to as a Wiener amplitude estimator.

22/129

Bearing fault identification

The aim of the identification step is to locate the localized fault through different Cyclostationary characteristic widely exists in the rotating or reciprocating machinery signal, especially the bearing vibration signal, which exhibits some hidden periodicity in its energy flow. McCormick and Nandi [START_REF] Mccormick | Cyclostationarity in rotating machine vibrations[END_REF] appears to be the first to apply the secondorder cyclostationary statistics in the monitoring of rotating machinery, and a number of cyclostationary tools such as the spectral correlation density, the Wigner-Ville spectrum and the degree of cyclostationarity were compared against traditional stationary tools to diagnose bearings. The literature [START_REF] Jerome Antoni | Cyclostationary modelling of rotating machine vibration signals[END_REF][START_REF] Bonnardot | Extraction of secondorder cyclostationary sources-application to vibration analysis[END_REF] of sidebands spaced at the speed of the gear on which the local fault is located [START_REF] Randall | A history of cepstrum analysis and its application to mechanical problems[END_REF].

The majority of such sidebands are only visible on a spectrum with a log amplitude scale, and so the cepstrum is an ideal way to extract the frequency information hidden in the sidebands. So the cepstrum can detect and quantify families of periodically spaced spectral components. Borghesani et al. [START_REF] Borghesani | Application of cepstrum pre-whitening for the diagnosis of bearing faults under variable speed conditions[END_REF] demonstrated the effectiveness of a cepstral based technique for prewhitening of the bearing signals thereby providing a relatively less complicated tool for identifying the bearing faults than the spectral kurtosis based methods.

Park et al. [START_REF] Park | Early fault detection in automotive ball bearings using the minimum variance cepstrum[END_REF] used a modified form of cepstrum analysis known as the minimum variance cepstrum obtained by liftering a logarithmic power spectrum by applying the minimum variance algorithm, for the early detection of faults in the ball bearings. It was observed that the minimum variance cepstrum could easily predict the fault period regardless of system frequency response or choice of optimal resonance bands subjugated mostly by the wavelet analysis methods. There are also some articles [START_REF] Jiang | Sparse dictionary design based on edited cepstrum and its application in rolling bearing fault diagnosis[END_REF][START_REF] Liu | Asymmetric penalty sparse model based cepstrum analysis for bearing fault detections[END_REF] that provide some more insight into the cepstrum based analysis and the application on rolling element bearings fault diagnosis.

Matrix decomposition technique is usually acted as one preprocessing tool to extract the signal feature or reduce the dimension in the vibration signal analysis. As reference [START_REF] Li | Feature extraction for rolling element bearing fault diagnosis utilizing generalized S transform and two-dimensional non-negative matrix factorization[END_REF] stated, the multidimensional data reduction methods such as independent component analysis, principal component analysis, and singular value decomposition fail to fully understand the multivariate character of bearing signals data. The authors [START_REF] Li | A new feature extraction and selection scheme for hybrid fault diagnosis of gearbox[END_REF][88] adopted an innovative method using the generalized S-transform for time frequency representation of the signal and a two-dimensional non-negative matrix factorization (NMF)

to reduce the dimension of the time frequency matrix and extract the desired features from a bearing signal. Such combination is a very common way in vibration signal analysis, as found in other similar methods such as NMF+SVM [START_REF] Yang | Application of the multikernel non-negative matrix factorization on the mechanical fault diagnosis[END_REF], NMF+kNN [START_REF] Yang | Discriminative non-negative matrix factorization (DNMF) and its application to the fault diagnosis of diesel engine[END_REF], NMF based 25/129 feature selection for classification [START_REF] Liang | Feature selection for machine fault diagnosis using clustering of non-negation matrix factorization[END_REF], etc. The above research only took advantage of dimensionality reduction decomposition. But matrix decomposition technique is more than just this. The analysis could take a further step for recovering the time domain denoised signal. In reference [START_REF] Wang | A Novel Signal Separation Method Based on Improved Sparse Non-Negative Matrix Factorization[END_REF], the improved sparse NMF was proposed to separate and extract the compound fault features of rolling bearings, and the time-frequency information was transformed into the time domain by using an inverse STFT. However, the phase recovery process was not introduced. Wodecki et al. [START_REF] Wodecki | Impulsive source separation using combination of Nonnegative Matrix Factorization of bi-frequency map, spatial denoising and Monte Carlo simulation[END_REF] used the cyclic spectral coherence as the two dimensional representation map. Then NMF was used to analyze such map in two ways: first, it helped to initially separate cyclic components by producing a set of filters for input vibration data, and second, to identify proper damage-related frequency components in envelope spectrum.

Signal decomposition mainly include empirical mode decomposition (EMD) [START_REF] Huang | The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF],

local mean decomposition (LMD) [START_REF] Smith | The local mean decomposition and its application to EEG perception data[END_REF], variational mode decomposition (VMD) [START_REF] Dragomiretskiy | Variational mode decomposition[END_REF]. In recent years, feature extraction based on signal decomposition technology has been widely used in fault diagnosis. Guo et al. [START_REF] Guo | Faulty bearing signal recovery from large noise using a hybrid method based on spectral kurtosis and ensemble empirical mode decomposition[END_REF] notified the EMD may be difficult to recover impulses from large noise. Further, the conventional EMD was reported to introduce mode mixing effect and the distortion of the faulty impulses. Based on these two observations, the authors therefore develop a hybrid signal processing method that combines spectral kurtosis (SK) with ensemble empirical mode decomposition (EEMD). Because the end effect and mode mixing problems also plague LMD, the authors [START_REF] Liu | Improved local mean decomposition for modulation information mining and its application to machinery fault diagnosis[END_REF][START_REF] Liu | Time-frequency representation based on robust local mean decomposition for multicomponent AM-FM signal analysis[END_REF] proposed a soft screening stop criterion that enables LMD to automatically find an optimal number of iterations for each screening process. In the proposed method, an objective function that considers two characteristics (the root mean square and the excess kurtosis) of the target signal was defined to automatically determine the optimal number of iterations.

The recently proposed VMD is to transform the signal decomposition problem into a constrained optimization problem [START_REF] Dragomiretskiy | Variational mode decomposition[END_REF]. So it overcomes the problems of end effect and mode component aliasing in EMD and LMD methods and has a more solid mathematical theory foundation.

Blind deconvolution is the technique to find an inverse filter, which can recover the original impact signal by reacting on the transmission path of the signal. At the same time, it can also enhance the impact component in the signal. Artificial intelligence is another powerful pattern identification tool, which has attracted great attention from many researchers and shows promise in bearing fault identification or classification applications. The step of fault identification amounts to mapping the information obtained in the feature space to bearing faults in the fault space.

Numerous artificial intelligence tools or techniques have been used, including mathematical optimization, as well as classification-and probability-based methods. Specifically, classifiers and statistical learning methods have been widely used in fault diagnosis of rotating machinery, that includes, k-NN algorithms [START_REF] Wang | K-nearest neighbors based methods for identification of different gear crack levels under different motor speeds and loads: Revisited[END_REF], Bayesian classifier [START_REF] Baraldi | Comparing the treatment of uncertainty in Bayesian networks and fuzzy expert systems used for a human reliability analysis application[END_REF], SVM [START_REF] Zhang | Intelligent fault diagnosis of roller bearings with multivariable ensemble-based incremental support vector machine[END_REF] and artificial neural network (ANN). Most recently, deep learning approaches have also began to be applied in the field of fault diagnosis [START_REF] Lei | Applications of machine learning to machine fault diagnosis: A review and roadmap[END_REF]. The effectiveness of these approaches largely depends upon the quality of features extracted from the bearing signals.

Bearing fault size characterization

And last but not least, bearing fault size is also an important information during the degradation process. A knowledge of when a bearing will fail -that is, its remaining useful life (RUL) -can serve as supplement to maintenance decision-making such as determining in advance the time an equipment needs to be taken out-of-service and that and is shown to be accurate for estimating a range of defect geometries over a range of shaft speeds and applied loads. The change in the static stiffness of the bearing and the applied load effect were also researched in [START_REF] Larizza | The importance of bearing stiffness and load when estimating the size of a defect in a rolling element bearing[END_REF]. However, as stated in reference [START_REF] Larizza | Defect size estimation in rolling element bearings with angled leading and trailing edges[END_REF], almost all previous researches are all based on the tested bearings with defects that have sharp 28/129 reasons why it has been so popular in the last decades. First the models are very rich in mathematical structure and hence can form the theoretical basis for use in a wide range of applications. Second the models, when applied properly, work very well in practice for several important applications [START_REF] Lawrence R Rabiner | A tutorial on hidden Markov models and selected applications in speech recognition[END_REF]. Although the number of reported researches based on HMM's is too large to discuss in detail here, a survey on the development of modeling the time-frequency features of nonstationary signal through HMM is given here.

HMM is a statistical model of the time series, and its basic theory was founded in the early 1960s. The HMM was firstly implemented in speech processing in the 1970s, and prosperous in the 1980s. Taking into account the unique characteristics of speech signals, some mainstream features, such as, perceptual linear prediction [START_REF] Juang | On the hidden Markov model and dynamic time warping for speech recognition-A unified view[END_REF], mel-frequency cepstral coefficients [START_REF] Chen | Cepstral domain talker stress compensation for robust speech recognition[END_REF], Vector Quantization [START_REF] Lawrence R Rabiner | On the application of vector quantization and hidden Markov models to speaker-independent, isolated word recognition[END_REF] etc., have been widely used to be modeled through HMMs. In addition to these features, short-time Fourier spectrum is also a good tool. This is because the short-time Fourier transformation and its inverse establish a oneto-one mapping between the time domain and the frequency domain, and FFT algorithms can be implemented efficiently [START_REF] Yang | Human action learning via hidden Markov model[END_REF]. Also, the short time Fourier transformation preserves information from the original signal, and ensures that important features are not lost as a 29/129 result of the transformation. Furthermore, the nonstationary signal can be considered as stationary over a short time interval. In reference [START_REF] Andrew | Comparison of time-frequency signal analysis techniques with application to speech recognition[END_REF], the authors had a comparison on these time-frequency features with their application to speech recognition. Juang

and Rabiner [START_REF] Hwang | Hidden Markov models for speech recognition[END_REF] firstly used the short-time Fourier spectral vectors to characterize the sequence of sounds, and applied HMM in speech recognition and observed an accuracy rate higher than 95% in speaker-independent tasks. Flores [START_REF] Flores | Continuous speech recognition in noise using spectral subtraction and HMM adaptation[END_REF] describes a scheme for robust speech recognition at poor signal to noise ratios. It consists of the spectral subtraction scheme followed by hidden markov model compensation framework. The estimated clean speech signal is then modeled through HMM for improving the recognition. Deng et al. [START_REF] Deng | Speech recognition using hidden Markov models with polynomial regression functions as nonstationary states[END_REF] modeled the time-varying spectral information as a deterministic trend which is essentially a low order regression (state-dependent) of the spectral coefficients with parameters estimated as an integral component of the HMM training process. One kind of continuous density HMM (CDHMM) [START_REF] Lawrence R Rabiner | A tutorial on hidden Markov models and selected applications in speech recognition[END_REF][START_REF] Tokuda | Hidden Markov models based on multi-space probability distribution for pitch pattern modeling[END_REF] is very common in the speech recognition.

In the CDHMM model, it is assumed that the feature vector distribution encountered at a given state can be modeled by a Gaussian distribution whose underlying mean represents the "true value" of the feature vector at that state. A more general approach is to associate with each state a weighted sum, or mixture of Gaussian distributions.

A complete statement on the speech processing through observation vectorization and CDHMM is described in the paper [START_REF] Knill | Hidden Markov models in speech and language processing[END_REF]. Each observation can be modeled as the single component Gaussian probability density or a more generalized Gaussian mixture model [START_REF] Knill | Hidden Markov models in speech and language processing[END_REF]. The combination of time-frequency features and HMM approaches have also been successfully for animal vocalizations analyses [START_REF] Somervuo | Parametric representations of bird sounds for automatic species recognition[END_REF][START_REF] Vieira | Call recognition and individual identification of fish vocalizations based on automatic speech recognition: An example with the Lusitanian toadfish[END_REF].

This analysis method is also applied in other fields, for example, robot application.

In order to judge the best trajectory of the robot, Tso [START_REF] Shiu | Demonstrated trajectory selection by hidden Markov model[END_REF] proposed to apply STFT as the feature extraction tool for HMM technique to measure human behaviour consistency. A similar method is used in human behavior recognition and robot control [START_REF] Yang | Human action learning via hidden Markov model[END_REF]. Furthermore, in the paper [START_REF] Wu | Fall detection based on sequential modeling of radar signal time-frequency features[END_REF], the HMM was utilized to characterize the time-frequency features of radar signals for the purpose of fall detection, and the feasibility of early warning was investigated. The transient signal detection through HMM is another interesting subject.

In paper [START_REF] Kundu | Transient sonar signal classification using hidden Markov models and neural nets[END_REF], the sonar signal is divided into a sequence of time-frequency frames, and each frame is represented by a feature in HMM for recognizing sonar transient signals.

Chen [START_REF] Chen | Detection of hidden Markov model transient signals[END_REF] modeled the transient-present observations as an HMM and the transientabsent as independent random noise, and combined it with a sequential probability ratio The first kind of approach is concerned about modeling the total life process and describe the degradation process for prediction [START_REF] Kwan | A novel approach to fault diagnostics and prognostics[END_REF][START_REF] Zhang | An integrated approach to bearing fault diagnostics and prog-128/129[END_REF][START_REF] Bulent Ayhan | Integrated Bearing Diagnostics and Fault Severity Estimation Using HMM[END_REF][START_REF] Tobon-Mejia | Hidden Markov models for failure diagnostic and prognostic[END_REF][START_REF] Sloukia | Bearings prognostic using mixture of gaussians hidden markov model and support vector machine[END_REF][START_REF] Sedira | Rotating machinery Diagnostic Using Hidden Markov Models (HMMs)[END_REF]. In the other kind of approach, several HMM models for every fault mode are trained, and then the trained models are exploited to do fault mode classification [START_REF] Ashar Aditiya | Fault diagnosis system of rotating machines using Hidden Markov Model (HMM)[END_REF][START_REF] Boutros | Detection and diagnosis of bearing and cutting tool faults using hidden Markov models[END_REF][START_REF] Xiao | A multichannel fusion approach based on coupled hidden Markov models for rolling element bearing fault diagnosis[END_REF][START_REF] Geramifard | An HMM-based semi-nonparametric approach for fault diagnostics in rotary electric motors[END_REF][START_REF] Chen | Fault diagnostics of helicopter gearboxes based on multi-sensor mixtured hidden Markov models[END_REF]. In [START_REF] Xin | An exploring study of hidden Markov model in rolling element bearing diagnostis[END_REF],

Xin proposed a method based on short-time Fourier spectrum and HMM for separating

Explicit duration hidden Markov model

As a rich mathematical structure and good models for real world signal, HMMs have became increasingly popular in recent decades. It is an important class of models that are successful in many application areas. However, the standard HMM allows self-transition process, then the time the Markov chain spends in a state is a random variable that is statistically described by a geometric probability mass distribution [START_REF] Howard | Dynamic probabilistic systems: Markov models[END_REF]. A method that does not require the estimation of the state duration time for prognosis was proposed in [START_REF] He | Probabilistic model based algorithms for prognostics[END_REF]. Dong et al. [START_REF] Dong | A novel approach to equipment health management based on autoregressive hidden semi-Markov model (AR-HSMM)[END_REF] presented a novel statistical learning approach based on auto-regressive EDHMM that took the correlation between observations into account for fault classification and RUL prediction. In the next few years, Dong et al. proposed a series of improved methods, for example the non-stationary EDHMM [START_REF] Dong | Equipment PHM using non-stationary segmental hidden semi-Markov model[END_REF], aging factors integrated into EDHMM [START_REF] Peng | A prognosis method using age-dependent hidden semi-Markov model for equipment health prediction[END_REF], joint with sequential Monte Carlo method [START_REF] Liu | A novel method for online health prognosis of equipment based on hidden semi-Markov model using sequential Monte Carlo methods[END_REF][START_REF] Liu | Online health management for complex nonlinear systems based on hidden semi-markov model using sequential monte carlo methods[END_REF], etc., for equipment health prediction. Some researchers consider that the system is subjected to degradation due to different operating conditions, leading to different rates in deterioration evolution. To take into account such problem, multi-branch or multi-sensor EDHMM [START_REF] Dong | Hidden semi-Markov model-based methodology for multi-sensor equipment health diagnosis and prognosis[END_REF][START_REF] Liu | A novel method using adaptive hidden semi-Markov model for multi-sensor monitoring equipment health prognosis[END_REF][START_REF] Trung | Multi-branch Hidden semi-Markov modeling for RUL prognosis[END_REF][START_REF] Trung | Prognosis based on Multi-branch Hidden semi-Markov Models: A case study[END_REF][START_REF] Chen | Fault prognosis of complex mechanical systems based on multi-sensor mixtured hidden semi-Markov models[END_REF] were proposed. In addition, another strategy is the combination of EDHMM with various assistant methods, such as the sequential Monte Carlo method for decreasing the computational and space complexity [START_REF] Liu | A novel method for online health prognosis of equipment based on hidden semi-Markov model using sequential Monte Carlo methods[END_REF][START_REF] Liu | Online health management for complex nonlinear systems based on hidden semi-markov model using sequential monte carlo methods[END_REF], AIC for evaluating the correct model configuration [START_REF] Cartella | Hidden semi-Markov models for predictive maintenance[END_REF], particle filter for predicting the online state [START_REF] Xiao | Online machine health prognostics based on modified duration-dependent hidden semi-Markov model and 127/129 high-order particle filtering[END_REF], etc. These application results show that better performance has been achieved by the combination strategy. In the articles of Dong [START_REF] Dong | Equipment health diagnosis and prognosis using hidden semi-Markov models[END_REF][START_REF] Dong | A novel approach to equipment health management based on autoregressive hidden semi-Markov model (AR-HSMM)[END_REF], the wavelet decomposition coefficients were utilized for the observations of EDHMM. But the goal is to develop trained EDHMMs to recognize N different states of a component for a given failure mode. Therefore, for the diagnosis of the fault, it is necessary that a separate EDHMM is trained for all possible fault types in addition to the EDHMM for normal conditions. After investigation, it is found that there is no article that really does fault diagnosis through the combination of EDHMM and time-frequency observations.

Integrated automatic diagnosis

The previous sections have focused mainly on the diagnosis techniques from different viewpoints. As mentioned before, the objective of this thesis is to present an integrated auto-diagnostic framework that includes fault detection, fault signal reconstruction, fault type identification, and fault size characterization. In this section, therefore, the automated diagnostic framework through the models introduced before is discussed and reviewed.

Bearing diagnosis is not simply to identify the fault type when a fault occurred. In the reality, bearing diagnosis is a big topic, and there are many sub-issues to be solved.

34/129

Integrated automatic diagnosis

For example, the health status of the collected vibration signal is often unknown. In the first step, therefore, the signal needs to be checked whether it is collected from a faulty bearing or healthy one. If the bearing is determined to be faulty, the fault identification is the next subtask. In some cases, the fault signal extraction is also an important task.

These are not enough, and the engineers also want to know the fault degree and whether the bearing can still be used. These problems all belong to the category of rolling element bearing CBM. Therefore, integrated diagnosis can be understood in general as "a field of research and application which aims at making use of present vibration signal in order to detect its degradation, diagnose faults, assess and pro-actively manage its failures". Based on this viewpoint, related references are reviewed.

The automation of fault diagnosis was primarily handicapped in the past due to the lack of appropriate techniques to represent the knowledge-based, symbolic reasoning of an expert diagnostician [START_REF] Boutros | Fault detection and diagnosis in machining processes and rotating machinery using fuzzy approach and hidden Markov model[END_REF]. The authors of [START_REF] Boutros | Fault detection and diagnosis in machining processes and rotating machinery using fuzzy approach and hidden Markov model[END_REF] presented a two-stage fault detection and diagnosis approach. The first stage is a generic and fast technique for detection using an index fusion and fuzzy logic approach; the second stage describes the development of a fault diagnosis method through the HMM 'codebook' training. In reference [START_REF] Ocak | A new bearing fault detection and diagnosis scheme based on hidden Markov modeling of vibration signals[END_REF][START_REF] Ocak | HMM-based fault detection and diagnosis scheme for rolling element bearings[END_REF],

Ocak proposed one method that can realize automatically fault detection and diagnosis by training HMM. Faults can be detected online by monitoring the probabilities of the pretrained HMM for the normal case given the features extracted from the vibration signals, and the HMM for which the probability is maximum determines the condition of the bearing. The bearing prognosis issue was integrated into Ocak's PhD thesis in [START_REF] Ocak | Fault detection, diagnosis and prognosis of rolling element bearings: frequency domain methods and hidden markov modeling[END_REF]. In reference [START_REF] Kwan | A novel approach to fault diagnostics and prognostics[END_REF], a fault diagnostics and prognostics algorithm based on HMM was proposed.

The algorithm trained different HMM models, and then achieved the fault diagnostics and prognostics in a unified framework. Similar research can be found in many references [START_REF] Zhang | An integrated approach to bearing fault diagnostics and prog-128/129[END_REF][START_REF] Bulent Ayhan | Integrated Bearing Diagnostics and Fault Severity Estimation Using HMM[END_REF][START_REF] Tobon-Mejia | Hidden Markov models for failure diagnostic and prognostic[END_REF][START_REF] Tobon-Mejia | A mixture of gaussians hidden markov model for failure diagnostic and prognostic[END_REF]. In reference [START_REF] Bulent Ayhan | Integrated Bearing Diagnostics and Fault Severity Estimation Using HMM[END_REF], the authors proposed an integrated framework based on HMM to incorporate fault diagnostics and fault severity estimation. EDHMM can also be used to achieve the integrated auto-diagnosis framework. In reference [START_REF] Dong | A segmental hidden semi-Markov model (HSMM)-based diagnostics and prognostics framework and methodology[END_REF], the proposed segmental EDHMM-based framework combines diagnostics and prognostics in an integrated manner. Reference [START_REF] Kass | Self-running bearing diagnosis based on scalar indicator using fast order frequency spectral coherence[END_REF] is another attempt, where the authors proposed an indicator based on the order-frequency spectral coherence to achieve self-running diagnosis, including detection, identification and classification of typical rolling element bearing faults.

Zuo [START_REF] Moghaddass | An integrated framework for online diagnostic and prognostic health monitoring using a multistate deterioration process[END_REF] proposed an integrated framework, and showed how to model the degradation process of a condition monitored device through the EDHMM and then presented how 35/129 to employ real-time condition monitoring data for online diagnostics and prognostics.

Reference [START_REF] Martin | Automated machine health monitoring at an expert level[END_REF] offers a recent discussion on integrated and automated machine health monitoring. After investigation, the limitation of current studies is the lack of diagnosis by modeling the nonstationary vibration signal. In other words, in most works in the literature, they trained the different HMMs in advance, and then use them for detecting or identifying or classifying. The second limitation is that most of them did not integrate the sub-issues mentioned above into one diagnostic framework.

Discussion

Based on the above literature survey, there are some points of worthy concern.

Through the available literature discussed in Chapter 2.3, it is clear that various advanced is capable of more automation to address different sub-issues.

z 2 z 3 z N-1 z N Hidden states z 1 Y N-1 Y N Y 3 Y 1 Y 2 ......

Conclusion

In this chapter, the basic bearing knowledge and bearing condition monitoring techniques have been reviewed firstly. Afterward, the advantages, disadvantages, and performance of the described different types of techniques have been discussed in detail.

37/129

While it is impossible to investigate all of the existing techniques, this chapter still aims to shed some light on the strengths and drawbacks of different techniques in each diagnostic issue. Based on the investigation provided in this chapter, it is evident that effective techniques already exist in different issues, but that there is still a gap in the integrated

39/129

Introduction

As we know, in most situations bearing fault characteristics cannot be measured directly as the bearing vibration signature is modified by the machine structure and this situation is further complicated by other vibration coming from other equipment on the machine. Such interferences often make the interpretation of signals difficult. In probability theory, stochastic processes are the main mathematical tool for describing uncertainty.

Compared with deterministic models, the stochastic models have gifted flexibility and adaptability to cope with the complexity of non-stationary phenomena. In this chapter, we now formally introduce the stochastic model, hidden Markov model and its variant.

We start by reviewing the basic definitions and concepts pertaining to Markov chains, and then the hidden Markov model and its special extension, explicit duration hidden Markov model. Based on the vibration signal, a detailed estimation and analysis of the model parameters is given subsequently. Some issues with respect to the input model parameters are discussed at last.

Markov chains

In reality, sequential data are ubiquitous, for example, the rainfall measurement or temperature on successive days at a particular location, the vibration signature at successive time frames, etc. To express such data in a probabilistic model, one of the simplest ways is offered by the Markov model. It is a stochastic model describing a sequence of possible events in which the probability of each event depends only on the observation attained in the previous event [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF]. It can also be called first-order Markov chain, which is depicted as a graphical model in Figure 3.1. Here we consider a sequential data Y 1:N ∈ M representing the specific stochastic process. The set M is the sample space, and the element in M is called an observation of the process. Then the Markov process can be expressed mathematically as

p(Y n+1 = j|Y n = i, Y n-1 = i n-1 , ..., Y 1 = i 1 ) = p(Y n+1 = j|Y n = i) = a ij , (3.1)
where, i, j, i 1 , ..., i n-1 ∈ M . The observations are assumed to be independent and identically distributed. It is noted that a ij is a fixed probability independent of time, representing the probability that the Markov process will make a transition to observation j given the 40/129

3.2. Markov chains current observation i. Since probabilities are nonnegative and since the process must make a transition into a specific observation, clearly one has an important constraint as

a ij ≥ 0, j a ij = 1, i, j ∈ M .
Let A denote the transition probability matrix of first-order Markov chain, so that

A =          a 11 a 12 • • • a 1M a 21 a 22 • • • a 2M . . . . . . . . . . . . a M 1 a M 2 • • • a M M          . (3.2) Y 1 Y 2 Y 3 Y 4
p(Y 1 , Y 2 , ..., Y N ) = N n=1 p(Y n |Y n-1 , ..., Y 1 ) = p(Y 1 ) N n=1 p(Y n |Y n-1 ), (3.3)
where p(Y 1 ) is the initial probability, and one constrain for the initial probability is

Y 1 ∈M p(Y 1 ) = 1.
Thus if we use such a model to predict the next observation in a sequence, the distribution of predictions will depend only on the value of the previous observation and will be independent of all earlier observations. Although the first-order Markov chain is very simple, it provides an important modeling idea and theoretic basis for afterword probabilistic model.

Here an example related to sequential data is used to illustrate the first-order Markov 41/129 chain. Suppose that the chance of rain tomorrow (denoted as observation Y n+1 ) depends on only whether or not it is raining today (denoted as Y n ), but not on past weather conditions. Suppose also that if it rains today, then it will rain tomorrow with probability p 1 ; and if it does not rain today, then it will rain tomorrow with probability p 2 . If we say that the process is in state 0 when it rains and state 1 when it does not rain, then the above-mentioned process is a two-state Markov chain whose transition probability matrix is given by,

A =    p 1 1 -p 1 p 2 1 -p 2    , (3.4)
and the transition probabilities are

p 1 = p(Y n+1 = 0|Y n = 0), p 2 = p(Y n+1 = 1|Y n = 0).

Hidden Markov model

However, for some cases, the observations are complex and infinite, and the corre- It is a constant matrix that provides a Markovian connection network between the hidden states, and models the variations in the duration of the signals associated with each state.

• Initial probability vector

This is different from a transition probability as it is the first state, and so it has a marginal distribution represented by a vector of probabilities π with elements

π i = p(z 1 = i), i ∈ M .

• Emission probability

This is a conditional probability from hidden state to the observed variable, p(Y n |z n ).

If the observation is discrete, this probability will be a matrix that models the transition from the state to the observation. If the observation is continuous, this probability further depends on the specific parameter φ according to a probability distribution. Because Y n is observed, the distribution p(Y n |z n ) consists of a vector of N numbers corresponding to the N possible states.

These HMM parameters can be described as the parameter set θ = {π, A, φ}, where, φ represents the emission probability matrix or the parameter in the observation distribution. For Gaussian observations, the distribution parameter φ include the covariance and the mean value. Then, the joint distribution for this model is given by,

p(Y 1:N , Z 1:N |θ) = p(z 1 ) N n=2 p(z n |z n-1 ) N n=1 p(x n |z n ). (3.5)
A useful way of interpreting and using HMMs is to consider each state of an HMM as a model of a segment of a stochastic process. One example about the vibration signal The estimation task of prime interest for HMMs is to estimate the model parameters θ and the hidden sequence {z n } from the observed sequence {Y n } in some optimal way.

We will discuss how to estimate the parameters in the later section of this chapter, and the convergence of the algorithm in the appendix.

Vibration signal modeling

HMM is a versatile stochastic model for nonstationary phenomena that involves transitions between hidden states. It is therefore well suited for describing the time-frequency spectrum coefficients of a faulty signal, whose probability distribution switches between one state corresponding to the occurrence of impulses and another one corresponding to dead zones [START_REF] Xin | An exploring study of hidden Markov model in rolling element bearing diagnostis[END_REF]. However, in the standard HMM the transition probability matrix is assumed constant, which implies that the time spent in a state is a random variable that 44/129 follows a geometric probability distribution [START_REF] Howard | Dynamic probabilistic systems: Markov models[END_REF]. For bearing signals, this is an unrealistic assumption. To fix this drawback, it is proposed to introduce an explicit duration hidden

Markov model [START_REF] Freguson | Variable duration models for speech[END_REF] (also known as "segmental HMM [START_REF] Russell | A segmental HMM for speech pattern modelling[END_REF]" and "variable-duration HMM [START_REF] Tang | Learning latent temporal structure for complex event detection[END_REF]", as part of the larger theory of hidden semi-Markov models (HSMM) [START_REF] Yu | Hidden semi-Markov models[END_REF]). In this section we introduce one useful variant, EDHMM. We will give a detailed introduction about EDHMM, and how to use it to model vibration signal.

Signal model and STFT decomposition

In

y(t) = n(t),
whereas in the faulty state

y(t) = x(t) + n(t)
is a superposition of transients and noise. The faulty signal will be modeled hereafter as a nonstationary signal that switches between two states: an "active" state where an impulse occurs and an "inactive" state where noise only is observable.

As discussed above, the continuous observation need to be segmented before modeling.

How to segment the vibration signal is an essential component of modeling, often having an important impact on the final result. The main concern in the segment action is the information extraction completeness. In the present work, the STFT is used as a versatile transform for decomposing nonstationary vibration signals with the ability of evidencing possible transitions between two different states. This is because the STFT decomposition and its inverse transformation establish a one-to-one mapping between the time domain and the time-frequency domain, and its oscillating basis functions are well-fitted to represent vibration signals. Also, a nonstationary signal can be considered as quasi-stationary over a short time interval. As it preserves information from the original signal, the STFT ensures that important features are not lost as a result of the transformation. In order to better use the STFT coefficients Y (n, k), three assumptions concerning Y (n, k) will be needed:

• i) Y (n, k) and Y (n , k ) are independent for any k = k and n = n ,

• ii) the expected value of the STFT coefficients is zero,

• iii) the real {Y (n, k)} and imaginary {Y (n, k)} parts are independent and identically distributed.

The assumptions are made to simplify the model, but without losing reasonableness.

The windowed signal will tend to be stationary, which means the frequency bin k and k under one window are uncorrelated. Then the reasonableness of assumption i is only depend on the window shift R. If the shift is large enough, n and n will be independent.

Assumption ii and iii are also due to the stationary characteristic of windowed signal.

Under mild conditions, the vector of STFT coefficients

Y n = [Y (n, 1), ..., Y (n, N k )]
converges in distribution to the complex-valued Gaussian distribution according to the central limit theorem [START_REF] Brillinger | Fourier analysis of stationary processes[END_REF], as expressed in Eq.(3.7). More reasons justifying the Gaussian assumption of the STFT coefficients can be found in reference [START_REF] Pearlman | Source coding of the discrete Fourier transform[END_REF]. As will be seen shortly, this property is advantageous in following modeling.

Y n ∼ CN (µ, C, M). (3.7)
The covariance of this distribution is

C = E[(Y n -µ)(Y n -µ) † ],
where † means the complex conjugate of the transpose, and the pseudo-covariance is

M = E[(Y n -µ)(Y n -µ) T ].
According to assumption i, the covariance of two random variables with different frequency is zero; according to assumption ii, the mean value µ is zero. Therefore, the covariance C can be simplified as a diagonal matrix with zeros in the off-diagonal elements as given by 46/129 Eq.(3.8):

C = E[(Y n -µ)(Y n -µ) † ] = E[Y n Y † n ] =          C 11 0 • • • 0 0 C 22 • • • 0 . . . . . . . . . . . . 0 0 • • • C N k N k          , (3.8)
where

C kk = | (Y (n, k))| 2 + | (Y (n, k))| 2 , k ∈ {1, 2, ..., N k } is
p(Y n ) = exp{-Y † n C -1 Y n } π N k |C| . ( 3.9)
It is worth noting that the STFT is used here to pre-process the data and not as a visual diagnostic tool -i.e. as with time-frequency based approaches; therefore, the uncertainty principle to which it is subjected -which limits the time-frequency resolutionis not really an issue. It is noted here that other transforms are also possible. The columns of the coefficient matrix Y n will be used as the observation sequence of the following modeling.

Explicit duration hidden Markov model

EDHMM [START_REF] Freguson | Variable duration models for speech[END_REF][START_REF] Yu | Hidden semi-Markov models[END_REF] is an extension of hidden Markov models that includes the case where the duration time in a state is described by a specific distribution, which is the origin of the term "explicit duration". It is used in this section to describe the temporal relationship among the STFT coefficients.

p(d i = k) = λ k i k! e -λ i , k = 0, 1, ..., (3.10)
where, the Poisson parameter λ i represents the average window number in state i, which

48/129

The transition probability matrix A in EDHMM is an 2 × 2 dynamic matrix that depends on τ n (i), rather than a constant matrix as in the standard HMM. The diagonal elements, a ii (τ n ), represent the recurrent state transition probabilities. These two elements are defined as the probability of remaining in the current state i at the next time step, given the time, τ n (i), already spent in the current state i. Since the duration time follows the Poisson distribution, the recurrent probability a ii (τ n ) will decrease with τ n . Then the recurrent probability a ii (τ n ) can be expressed through the Poisson cumulative distribution function (CDF) as,

a ii (τ n (i)) = p(z n+1 = i|z n = i, τ n (i)) = p(z n+1 = i|z n = i, ..., z n-τn(i)+1 = i, z n-τn(i) = i) = p(z n+1 = i, ..., z n-τn(i)+2 = i|z n-τn(i)+1 = i, z n-τn(i) = i) p(z n = i, ..., z n-τn(i)+2 = i|z n-τn(i)+1 = i, z n-τn(i) = i) = p(d i > τ n (i)) p(d i > τ n (i) -1) = 1 -p(d i ≤ τ n (i)) 1 -p(d i ≤ τ n (i) -1) = 1 -F (τ n (i)|λ i ) 1 -F (τ n (i) -1|λ i ) , (3.11)
where F (τ n (i)|λ i ) represents the CDF of Poisson distribution with parameter λ i . The numerator, 1 -F (τ n (i)|λ i ), represents the probability that the explicit duration Markov process, at datum n, has stayed in state i for at least τ n (i) samples. [START_REF] Antoni | The spectral kurtosis: a useful tool for characterising non-stationary signals[END_REF], it is found that the recurrent transition probability a ii (τ n (i)) will gradually decrease as τ n (i) increase. Eventually, the Markov chain will transit to state j from state i, and then go to the next Poisson process.

According to the recurrent probability defined above and the constraint j a ij = 1, the dynamic transition probability matrix A(τ n (i)) can be written as where, the two off-diagonal elements, a ij (τ n )(i = j), is defined as the transition probabilities between two different states, given the previous state i has lasted τ n (i). For the two-states case, the relationship between the recurrent probability and non-recurrent probability is fixed,

A(τ n (i)) = [a ij (τ n (i))] = [p(z n+1 = j|z n = i, τ n (i))] =    a 11 (τ n (1)) 1 -a 11 (τ n (1)) 1 -a 22 (τ n (2)) a 22 (τ n (2))    . ( 3
a ij (τ n ) = 1 -a ii (τ n (i)), under the constraint j a ij = 1.
In this definition, τ n (i) can be adaptively optimized through the duration distribution parameters in each datum n. Therefore, the dynamic matrix A(τ n (i)) will also vary regularly according to the Poisson parameter, which does not need to be estimated in the following re-estimation section. Then the state sequence generated by this dynamic matrix will have In EDHMM, transitions between states start from an initial probability (n = 1) denoted π i , i ∈ {1, 2}, defined as π i = p(z 1 = i|Y 1:N ) with the constraint i π i = 1. The probabilities of different states are then updated according to transitions and observations.

The modeling stage is completed up to now, and the EDHMM is completely described by the set of parameters θ = {π 1 , λ, C} (note that π 2 = 1 -π 1 ), where λ = [λ 1 , λ 2 ] and 51/129

C = [C 1 , C 2 ]
. The next section addresses the estimation of these parameters, and explains what characteristics of the signal these parameters reveal.

Parameter estimation and analysis

Parameter estimation

The goal is to find the optimal model parameter set θ that maximizes the likelihood p(Y 1:N |θ) from the observations Y 1:N . The forward-backward recursion formulas are first introduced to calculate the likelihood. The forward and backward variables at datum n are defined as

α n (i) = p(Y 1:n , z n = i|θ), (3.13) β n (i) = p(Y n+1:N |z n = i, θ). (3.14)
The quantity α n (i) represents the joint probability of observing all of the given data up to time n and the value of state i, where β n (i) represents the conditional probability of all future data from time n + 1 up to N given the value of state i. Therefore, α and β are the forms of 2×N matrices for the two states model.

In order to calculate the forward and backward variables, the recursion relations are derived using the conditional independence properties together with the sum and product rules. The forward quantity α n (i) is expressed in terms of α n-1 (i) as,

α n (j) = p(Y 1 , Y 2 , ..., Y n |z n = j)p(z n = j) = p(Y n |z n = j)p(Y 1 , Y 2 , ..., Y n-1 |z n = j)p(z n = j) = p(Y n |z n = j) z n-1 =i p(Y 1 , Y 2 , ..., Y n-1 , z n-1 = i, z n = j) = p(Y n |z n = j) z n-1 =i p(Y 1 , Y 2 , ..., Y n-1 , z n |z n-1 = i)p(z n-1 = i) = p(Y n |z n = j) z n-1 =i p(Y 1 , Y 2 , ..., Y n-1 |z n-1 = i)a ij (τ n-1 )p(z n-1 = i) = p(Y n |z n = j) z n-1 =i p(Y 1 , Y 2 , ..., Y n-1 , z n-1 = i)a ij (τ n-1 ) = p(Y n |z n = j) i α n-1 (i)a ij (τ n-1 ) . (3.15)

52/129

With the same manipulation, the recursion relation for quantity β n (i) can be formulated as,

β n (i) = j β n+1 (j)a ij (τ n )p(Y n+1 |z n+1 = j). (3.16)
It is seen these two recursion formulas contains the dynamic duration time τ n (i). Therefore, to calculate the formulas above, the update formula for τ n (i) is required [START_REF] Azimi | Offline and online identification of hidden semi-Markov models[END_REF], as follows,

τ n (i) = p(z n = i|Y 1:n , θ)τ n-1 (i) + 1 = α n (i) i α n (i) τ n-1 (i) + 1. (3.17)
From now on, several other posterior probabilities are defined for facilitating the subsequent re-estimation. First, the posterior probability of the current state z n being equal to i at datum n, given the observations Y 1:N and parameter set θ, denoted as γ n (i), is expressed as

γ n (i) = p(z n = i|Y 1:N , θ) = p(Y 1:N |z n = i, θ)p(z n = i|θ) p(Y 1:N |θ) = p(Y 1:n , z n = i|θ)p(Y n+1:N |z n = i, θ) p(Y 1:N |θ) = α n (i)β n (i) p(Y 1:N |θ) . ( 3

.18)

The sequence γ 1:N (i) will be used for identifying the fault type in Chapter 4.4. According to this equation, there are some useful remarks worth noting. On one hand, if setting n = N and replacing α N (i), we can obtain one relationship as

γ N (i) = p(z N = i|Y 1:N , θ) = p(Y 1:N , z N = i|θ)β N (i) p(Y 1:N |θ) = p(z n = i|Y 1:N , θ)β N (i), (3.19)
from which it can be seen that β

N (i) = γ N (i)
p(zn=i|Y 1:N ,θ) = 1 for all the states i. This inference will be useful in the derivation of likelihood function. As we know, the likelihood function p(Y 1:N |θ) is the indicator for the convergence criterion of the algorithm, so it drives the 53/129 decision as whether the iterative algorithm is to be stopped or not. On the other hand, if we sum both side of Eq. (3.18), and use the fact that the left-side is a normalized distribution, the likelihood function p(Y 1:N |θ) can be written as

p(Y 1:N |θ) = i α n (j)β n (i).
(3.20)

Then the likelihood function can be estimated for any convenient choice of n. In order to simplify the algorithm, n = N is selected for removing β n (i) according to the inference form Eq.(3.19), i.e, β N (i) = 1. The likelihood function can be simplified as

p(Y 1:N |θ) = i α N (i). (3.21)
Next, another useful posterior probability, denoted as ξ n (i, j), is the probability of a transition from state i to state j given the observations Y 1:N , which can be expressed in terms of the forward-backward variables and model parameters as,

ξ n (i, j) = p(z n = i, z n+1 = j|Y 1:N , θ) = p(Y 1:N |z n = i, z n+1 = j, θ)p(z n = i, z n+1 = j|θ) p(Y 1:N |θ) = α n (i)β n+1 (j)a ij (τ n (i))p(Y n+1 |z n+1 = j, θ) p(Y 1:N |θ) . ( 3

.22)

According to the above formulas, the parameters θ = {π 1 , λ, C} can be estimated as follows. The initial state probability

π i estimate is πi = γ 1 (i)/ i γ 1 (i), (3.23)
where the numerator γ 1 (i) can be understood as the probability of occurrence of z 1 = i, and the denominator, i γ 1 (i), as the possibility of z 1 over all states. The estimate of the In the Poisson distribution, the maximum likelihood estimate of the parameter from the samples is the sample mean. Therefore, the parameter λ in the duration distribution is estimated as λi =

covariance C i is (See details in appendix A), Ĉi = N n=1 γ n (i)Y n Y † n N n=1 γ n (i) . ( 3
N -1 n=1 ξ n (i, j)τ n (i) (i = j), (3.26)
where the weight ξ n (i, j) represents the probability of transition from state i to state j, 3) Check for the stop criterion of either the logarithmic likelihood or the parameters. If the stop criterion is not satisfied, the let θ old ← θ new and return to step 2).

and
The pseudo-code is summarized in Algorithm 1. It is noted that it involves two parameters: max_iter for the maximum number of iterations and Threshold for the expected relative tolerance between the logarithmic likelihood probabilities of two consecutive times.

These two parameters can be easily set by default. Also, the STFT parameters, N w , R, and the EDHMM parameter θ = {π 1 , λ, C} will be discussed in detail in Chapter 3.6.

Parameter analysis

The parameters of EDHMM be estimated in the previous section, will be used to

(i) = αn(i) i αn(i) τ n-1 (i) + 1 γ n (i) = αn(i)βn(i) p(Y 1:N |θ) ξ n (i, j) = αn(i)β n+1 (j)a ij (τn(i))p(Y n+1 |z n+1 =j,θ) p(Y 1:N |θ) → Update the parameter θ = {π 1 , λ, C} πi = γ 1 (i)/ i γ 1 (i) Ĉi = N n=1 γn(i)YnY † n N n=1 γn(i) λi = N -1 n=1 ξ n (i, j)τ n (i) (i = j) → Stop criteria if log p k (Y 1:N )-log p k-1 (Y 1:N ) | log p k (Y 1:N )|+| log p k-1 (Y 1:N )| <Threshold then Stop else k ← k + 1; θ old ← θ new Endif

Covariance Ĉ

The positive diagonal matrix Ĉi is the covariance of the observations in state i, as described in Eq. (3.8). The diagonal elements,

[ Ĉi ] kk = | (Y (n, k))| 2 + | (Y (n, k))| 2 ,
k ∈ {1, 2, ..., N k }, are the estimates of the variance at frequency k. Therefore, this parameter is a good indicator for displaying the energy distribution in the frequency domain [START_REF] Xin | Semi-automated diagnosis of bearing faults based on a hidden Markov model of the vibration signals[END_REF]. It can also be used as a criterion for deciding whether the signal is healthy is because the posterior probability is not sensitive to signal amplitude, so it will avoid the influence of signal modulation and lead to a clear spectrum. This advantage will be demonstrated in the validation chapter 5. In addition, it is found that there are two transient pulses wrongly estimated in Figure 3.9(b) as the roller was in the non-load zone so that the amplitude was small and covered by the noise. But the state sequence incompleteness will affect less the diagnosis result, identification and characterization, which will be discussed in Chapter 4.6. From this figure, some remarks with respect to the parameter λ are worthy of attention. 1) Both λ 1 and λ 2 increase when the fault size increases, which is counterintuitive. As the time period (intervals between the impacts) is fixed, λ 1 should theoretically decrease as the fault increases. This may be due to the transient pulses not being detected, and being identified wrongly as noise states.

Poission parameter

2) It is found that, under the same fault size, λ 2 of the outer race is greater than the λ 2 of the inner race. It is also found that the rotating speed has an impact on duration time, but not much. From the figure, we can see that the features (λ 1 , λ 2 ) are good at separating different samples (different fault type, different fault size). In this sense, they could be used as the features for classification. 60/129

Input parameter settings

In order to be fully automated, a methodology should rely on a minimum of hyperparameters. The hyperparameters introduced so far relate to the STFT and the EDHMM.

On the one hand, the STFT is characterized by the window type, its length N w , and its shift R. It is known to be relatively insensitive to the window type, provided the window is a smooth function. Given a window type, the value of the shift should be small enough to prevent any loss of information. A typical choice is a Hann window with R = N w /4. Hence, the window length N w is the only critical parameter. It directly controls the frequency resolution ∆f = F s /N w . As discussed in the reference [START_REF] Xin | Semi-automated diagnosis of bearing faults based on a hidden Markov model of the vibration signals[END_REF], the window length should be shorter than the fault period (interval between two impacts).

Hence, this condition can possibly be met given a prior knowledge of the potential faults and their associated characteristic frequencies. A lower bound is surely more difficult to find; it is data-dependent in general as it corresponds to the coarsest frequency resolution required to analyze the signal.

On the other hand, the EDHMM involves initial values for the unknown paramaters π 1 , C 1 , and λ i , i = 1, 2 (remember that π 2 and C 2 are easily deduced from π 1 and C 1 according to formula (3.25)). There are various possible initialization strategies. An ad hoc strategy is discussed in [START_REF] Xin | Semi-automated diagnosis of bearing faults based on a hidden Markov model of the vibration signals[END_REF] for a two-state HMM of the STFT coefficients, which also applies here. Another popular strategy is to first cluster the data with the K-means algorithm and then get estimates from the obtained clusters. A last possibility is to resort to random initialization. Intensive numerical experiments have shown that the proposed EDHMM is robust to initialization, whatever the strategy used, because it involves only two states. It is therefore concluded that initialization of the EDHMM will not significantly affect the finial results. In conclusion, N w remains the sole critical hyperparameter of the proposed methodology.

Conclusion

This Subsequently, the identification of the fault type is achieved through one statistical counting method, and the posterior probability spectrum is also introduced. Section 4.5 introduces the characterization of the fault size using the duration time in different states. One discussion based on this integrated framework is given in Section 4.6.

Fault detection

Early detection is capable for examination of the fault occurrence and prediction of the fault evolution of the failing component before the fault progresses to a state that endangers the system's operational integrity. Timely and correct decision making is important for the real-time monitoring of crucial system. Therefore, the first step is to check whether the recorded vibration signal belongs to a healthy bearing or if it indicates the presence of a fault, in which case subsequent fault identification is needed.

As seen in the investigation of Chapter 2. indicators, the likelihood ratio naturally follows a chi-squared distribution, which provides the mathematical basis for finding the alarm threshold. The detailed process is formulated as follows. First, the two alternative hypotheses, H 0 and H 1 , are characterized by the following likelihood functions:

                            
Null hypothesis H 0 (healthy signal with state 1 only):

p(Y 1:N |H 0 ) = N n=1 CN (Y n ; 0, C H 0 )
Alternative hypothesis H 1 (faulty signal with intermittent states 1 and 2):

p(Y 1:N |H 1 ) = N n=1 [γ n (1) • CN (Y n ; 0, C H 1 1 ) + γ n (2) • CN (Y n ; 0, C H 1 2 )].
Since the theoretical covariance matrices are unknown, they have to be estimated under the two alternative hypotheses. On the one hand, under H 1 , the covariance matrices

C H 1 i , i = 1, 2
, are estimated by Eq.(3.23) as explained in Chapter 3.5. On the other hand, under H 0 , the estimate of

C H 0 is ĈH 0 = N n=1 Y n Y H n N = ĈH 1 1 + ĈH 1 2 , (4.1)
where the last equality results from Eq. (3.24). As explained in Chapter 3.5.2, the covariance associated to one state is a good indicator for displaying the energy distribution in the frequency domain. Therefore, the covariances C H 1 i and C H 0 will be similar for the null hypothesis, while very different under the alternative hypothesis. Based on this feature, the statistical indicator, generalized likelihood ratio (GLR), can be constructed by plugging in these estimates into the probabilities p(Y 1:

N |H i ), i = 0, 1 as, Λ = ln p(Y 1:N |H 1 ) p(Y 1:N |H 0 ) = ln N n=1 [γ n (1) • CN (Y n ; 0, ĈH 1 1 ) + γ n (2) • CN (Y n ; 0, ĈH 1 2 )] N n=1 CN (Y n ; 0, ĈH 0 ) . ( 4.2)
The principle of the test is to accept the hypothesis H 1 if the GLR Λ is greater than a given alarm threshold. According to Wilk's theorem, twice the GLR asymptotically follows a chi-squared distribution with number of degrees of freedom, denoted v, equal to the difference between the number of unknown parameters under hypotheses H 1 and H 0 ,

65/129 that is 2Λ ∼ χ 2 v . (4.3)
In the null hypothesis H 0 , there are only N k unknown elements in the covariance C H 0 , whereas for hypothesis H 1 , 2N k for the two covariance matrices and N for the posterior probability γ n (noting that γ n (1) + γ n (2) = 1). So the difference in degrees of freedom is

v = N k + N . Therefore, the null hypothesis H 0 is rejected if Λ > χ 2
v,1-α /2 at the risk α, where χ 2 v,1-α is the quantile of the chi-squared distribution with probability 1 -α.

Fault signal reconstruction

Spectral amplitude estimation forms the basis of many signal restoration systems, such as for vibration signal denoising. For restoration of time-domain signals, an estimate of the instantaneous magnitude spectrum is combined with the phase of the noisy signal and then transformed via an inverse discrete Fourier transform to the time domain. The 66/129 related methods are covered in Chapter 2.3.2. Bayesian spectral amplitude estimation methods offer substantial performance improvements on spectral subtraction by utilising the probability density functions of the signal and noise process. The hidden states associated with signal spectra provide useful information for the evolution of signal spectra and the correlation between those spectra. Therefore, HMM can provide a mathematically tractable model basis for signal enhancement. In this section, the EDHMM based denoising method is formulated in the following.

The basic goal of feature enhancement is to extract the clean signal x(t) from the measurement y(t) with additive background noise n(t). Firstly, according to the Bayesian theory and EDHMM theory, the posterior probability distribution of the nth observation X n can be written as

p(X n |Y n , θ) ∝ p(Y n |X n , θ)p(X n |θ),
where, X n represents the Fourier coefficients of the fault signal at nth window. For the specific frequency k, the fault signal spectrum posterior probability is

p(X(n, k)|Y (n, k), θ) ∝ exp{ -|Y (n,k)-γnX(n,k)| 2 Ĉ1 (k) } exp{ -|X(n,k)| 2 Ĉ2 (k)-Ĉ1 (k) } π 2 Ĉ1 (k)( Ĉ2 (k) -Ĉ1 (k)) ∝ exp{-[γ 2 n ( Ĉ2 (k)-Ĉ1 (k))+ Ĉ1 (k)]X(n,k) 2 -2γn( Ĉ2 (k)-Ĉ1 (k))Y (n,k)X(n,k)+( Ĉ2 (k)-Ĉ1 (k))X(n,k) 2 Ĉ1 (k)( Ĉ2 (k)-Ĉ1 (k)) } π 2 Ĉ1 (k)( Ĉ2 (k) -Ĉ1 (k)) , (4.4)
where, Ĉ1 and Ĉ2 are the covariances of the observation in different states, and they can be estimated through Eq.(3.23). Then, the estimation of covariance of the signal of interest X n is obtained as Ĉ2 -Ĉ1 . From another way, the STFT coefficients of the fault signal x(t) converge in distribution to the complex-valued Gaussian distribution according to the central limit theorem [START_REF] Brillinger | Fourier analysis of stationary processes[END_REF]. According to the certain assumptions, the STFT coefficients X(n) at a given datum n follow the complex Gaussian distribution. For the specific frequency point k, the posterior probability distribution can be written as,

p(X(n, k)|Y (n, k), θ) = exp{-|X(n,k)-µx(k)| 2 Cx(k) } πC x (k) . ( 4.5)
Combining the two formulas (4.4) and (4.5), the corresponding mean µ x (k) and 67/129 covariance matrix C x (k) of the posterior probability distribution can be given as

             C x (k) = Ĉ1 (k)( Ĉ2 (k) -Ĉ1 (k)) γ 2 n ( Ĉ2 (k) -Ĉ1 (k)) + Ĉ1 (k) µ x (k) = C x (k) Ĉ1 (k) γ n Y (n, k).
If one substitutes the covariance C x (k) into the mean µ x (k), the expectation of the fault signal spectrum at datum n can be obtained as,

E{X n |Y n , θ} = µ x = γ n ( Ĉ2 -Ĉ1 ) γ 2 n ( Ĉ2 -Ĉ1 ) + Ĉ1 Y n . (4.6)
Finally, the time-domain signal of interest, x(t), is obtained from Eq.(4.6) by using the inverse STFT. Compared with the Wiener filter as Eq.(2.10), this can be seen as a time-varying filter, as which takes account of the posterior probability γ n .

The technique of fault signal reconstruction is assessed based on its ability to recover the fault signal when the signal is disturbed by noise. Two different SNRs were chosen, 0 dB and -5 dB. For higher SNRs, the bearing fault signal impacts are more prominent and for lower SNRs the bearing fault impacts are easily masked by background noise.The noise immunity of the method will be discussion in the last section.

Fault type identification

Once a fault is detected, the type and severity of the fault will become the primary issue. Spectral analysis of vibration signal is widely used in bearing fault analysis. Most advanced techniques all rely on demodulation and spectrum analysis for finding the bearing

Posterior probability spectrum

The components in the vibration signal of a faulty bearing mainly consist of damped and repetitive impulses on the top of stationary background noise. This has been modeled in Chapter 3.4.2 as statistical properties switch between two states (the active and inactive states). The transition frequency between these two states relates to the bearing fault frequency, which is itself expressed in terms of the geometrical parameters of the bearing and is specific to each type of fault. Therefore, the idea is to identify the type of fault by performing the spectral analysis of the posterior probability of the active state,

S(k) = |F{γ n (2)}| 2 = | N n=1 γ n (2)e -j2πk n N | 2 , (4.7)
where F{•} represents the discrete Fourier operator. The advantage of the PPS is that it is only related to the state information of the vibration signal, which is more conducive to spectrum analysis than the original signal. In addition, the PPS is not sensitive to the amplitude of the transient pulses, so it can avoid the influence of modulation, thereby reducing the sidebands in the spectrum. This advantage had been illustrated by an exemplary signal with inner race fault shown in Fig. 4.4. It is shown in the subplot (a)

that the characteristic frequency f BP F I and its harmonic are very clear. Compared to the squared envelop spectrum in subplot (b), most sidebands or amplitude modulation have been removed in the PPS. This advantages will be further explained in Chapter 3.5.2.

In respect of spectrum analysis, it is similar to the envelope spectrum, a cutting-edge tool widely used in vibration-based condition monitoring [START_REF] Randall | Rolling element bearing diagnostics-A tutorial[END_REF]. However, contrary to the latter, the PPS does not rely on manual pre-processing for the selection of informative bands in the signal, as typically achieved by means of finely tuned bandpass filters (see e.g. [START_REF] Antoni | A Critical Overview of the "Filterbank-Feature-Decision" Methodology in Machine Condition Monitoring[END_REF], [START_REF] Hebda-Sobkowicz | Selection of the Informative Frequency Band in a Bearing Fault Diagnosis in the Presence of Non-Gaussian Noise-Comparison of Recently Developed Methods[END_REF], [START_REF] Niehaus | NIC Methodology: A probabilistic methodology for improved informative frequency band identification 122/129 by utilizing the available healthy historical data under time-varying operating conditions[END_REF]). The advantage of the PPS will be demonstrated through validation in the next chapter.

The spectral analysis of the state sequence is a valuable tool for identifying the bearing fault type, yet it still requires a visual inspection. If this task is to be automated as well, one can still resume the spectrum by scalar indicators (see e.g. [START_REF] Borghesani | The relationship between kurtosis-and envelope-based indexes for the diagnostic of rolling element bearings[END_REF]) or use automatic spectral analysis methods as introduced in [START_REF] Firla | Automatic characteristic frequency association and all-sideband demodulation for the detection of a bearing fault[END_REF]. Another strategy is proposed hereafter in order to automatically identify the fault type among a given set of candidates.

Statistical analysis

Spectrum analysis is one common and effective method for identifying the bearing fault type. But the disadvantage is that it is not automatic. Finally, it still depends on the visual examination of the end-user. Although some methods can realize automatic scanning, they are difficult to implement. In order to completely escape the manual intervention, this section introduces one simple method based on a statistical point of view for identifying automatically the fault type. Since the EDHMM returns the time τ (i) already spent in the current state i at each instant n, the average period between two successive impulses can also be estimated, thus giving an indication of the characteristic fault frequency. More specifically, let t s be the elapsed time between the s-th and the (s + 1)-th impulses (see Figure 4.5) and let f c be the expected characteristic frequency of a fault type c. Then, accounting for the possible presence of uncertainties in the actual value of f c (unpredictable dependence on load, presence of slippage, etc.) and in the measurement of t s (estimation errors, effect of noise), one can count the number of times t s falls in a narrow frequency interval

[f c -σ, f c + σ]
of width 2σ centered on f c :

N c = card({t s : 1/t s ∈ [f c -σ, f c + σ], s ∈ {1, 2, ..., S}}), ( 4.8)
where card(A) means the number of elements in the set A. If a number C of exclusive fault types are considered (e.g. inner-race fault, outer-race fault, ball-fault), the posterior 72/129 probability of observing fault c, c = 1, ..., C given t 1:S is then estimated by

p(f c |t 1:S ) = N c C n=1 N n . ( 4.9)
Since the equations to calculate the characteristic frequencies are based on bearing geometry and speed alone, variations due to loading and slipping are not considered.

σ is the frequency error, which extends the characteristic frequency to a characteristic interval; it is typically about 5% of the fault characteristic frequency. Therefore, automatic identification of the type can be achieved by selecting the highest empirical posterior probability, i.e. finding the maximum probability is to identify the fault type. In addition, the samples 1/t 1:S can be also represented in the form of a histogram, and then the fault identification issue is equivalent to finding the highest bin.

Fault size characterization

The fault type can be determined through the techniques defined in the previous section, but the fault severity is still unknown. In the characterization stage, the size of the bearing fault is of concern. A proportional relationship between the transient pulse duration and fault size seems a reasonable assumption. In reference [START_REF] Sawalhi | Vibration response of spalled rolling element bearings: Observations, simulations and signal processing techniques to track the spall size[END_REF]

≈ λ 2 R F s , (4.10)
with F s the sampling frequency of the signal. The value of ∆t reflects the duration of the passage of the fault in the contact zone of the bearing. In practice, it is likely to be slightly over-estimated since it includes the coda, say ∆e (in seconds), of the impulse response of the structure. However, if needed, the latter can possibly be measured by other means and subtracted from ∆t. For instance, a reasonable estimate of the coda can be obtained by measuring the response of the structure after being impacted by a hammer. Another simple method -used hereafter -is to deduce the value of ∆e from measurement taken at different rotation speeds. Besides, the precision of ∆t is on the order of the window shift, i.e. R/F s . This is surely less accurate than the precision that could be achieved with other methods such as blind deconvolution [START_REF] Sawalhi | The enhancement of fault detection and diagnosis in rolling element bearings using minimum entropy deconvolution combined with spectral kurtosis[END_REF][START_REF] Randall | Signal processing tools for tracking the size of a spall in a rolling element bearing[END_REF][START_REF] Buzzoni | Blind deconvolution based on cyclostationarity maximization and its application to fault identification[END_REF][START_REF] Buzzoni | Development of an indicator for the assessment of damage level in rolling element bearings based on blind deconvolution methods[END_REF]], yet the latter usually comes with a high demand on the user's expertise. Sacrificing a degree of precision is the price to pay for having an automated method.

For a faulty bearing with stationary outer race, the fault size l (in meter) is related to the transient duration ∆t as follows,

Outer race fault:

l ≈ πf F T F (∆t -∆e) • D (4.11)
Inner race fault:

l ≈ πf F T F (∆t -∆e) • D(1 + d/D) (4.12)
Rolling element fault:

l ≈ πf F T F (∆t -∆e) • D(1 + d/(2D)) ≈ πf BSF (∆t -∆e) • d (4.13)
where D and d stand for the pitch and element diameters of the bearing, respectively, and f F T F for the fundamental train frequency. It is worth noting that a fault on a rolling 74/129 element actually produces impacts of slightly different durations on the outer and inner races. Formula (4.11) applies to the former case and formula (4.12) to the latter case.

Formula (4.13) is an average, when both races are impacted. Noting that when several rotation speeds are available, say f j r , j = 1, 2, ..., the above formulas can be cast in the general form ∆t j = l • c(f j r ) + ∆e, where c(f j r ) is a function of f j r . Therefore, the fault size l and the coda error ∆e can be estimated as the slope and the intercept, respectively, of the fitted lines between the duration times ∆t j and c(f j r ), as shown in Figure 4.6. One advantage of this technique is that regression will decrease the error produced in the estimation and transformation of ∆t, and this advantage will become more obvious as the fitted samples increase. Although this technique still cannot accurately quantify the size, it will be useful for condition monitoring as demonstrated hereafter in the next chapter.

0 c( f r 1 ) c( f r 2 ) c( f r 3 ) c( f r 4 )
Rotating speed

t 1 t 3 t 2 t 4
Duration time e Linear regression and returns elements of information on 1) detection, 2) reconstruction, 3) identification and 4) size characterization of a potential fault. The detection step returns a binary output, which can be possibly substantiated by the likelihood ratio Λ (or a p-value). The higher the value of Λ with respect to a threshold, the lower the risk of false alarm (see Eq.( 4.3)). If a fault is detected, the identification step returns the PPS, the equivalent of an envelope spectrum, which can be analyzed visually or automatically. In order to be fully automatic, a simple method that calculates the posterior probability of different fault types for identifying the fault type is also introduced. Eventually, the fault signal can be extracted through a time-varying filter constructed from the covariance matrix of the observations. Independently of the latter, the fault characterization step returns an estimate of the fault size.

It is emphasized again that the complete algorithm relies only on one critical hyperparameter, N w , the inverse of the frequency resolution. In particular, it does neither rely on historical data, nor on finely tuned pre-filters, nor on trial-and-error manipulation of time-frequency distributions. In addition, these four tasks are independent of each other, in the sense that any one of the tasks can be accomplished effectively without the aid of other. Therefore, it can be understood as a modular integrated framework. However, one high risk in this diagnosis framework is the strong dependency on the stochastic model.

Only if the model is well established, the integrated diagnosis will perform well.

In practical applications, the vibration signal is always distorted by all kinds of noise.

For example, when a roller goes to the non-load zone, the magnitude of the transient pulse will be low and easily covered by noise, which may be estimated wrongly as noise state. As shown in Figure 3.9(b), it is seen that two transient pulses with low energy that have been estimated wrongly. But the wrong estimation of a state or the state sequence incompleteness affect little the final result. The reason is because the statistical analysis is based on statistics comparison, it is rather unaffected by the hidden state incompleteness.

The example in Figure 3.10 is used here to explain. In this synthetic signal, there are some transient pulses that are totally covered by noise. However, for the statistical analysis based identification, the probability of inner race fault is the highest in the posterior 76/129

Fault detection

Compute the GLR and compare to the threshold Accept H 0

Yes Stop

No

Fault signal Reconstruction

Construct the time-vary filter

Fault type Identification

Compute the fault type probability or posterior probability

Fault size Characterization

Estimate the fault size l Hz, which includes the characteristic frequency f BP F I . The second highest bin is [76 86] Hz, which is due to the doubling period because of the wrong estimation of the state.

With respect to the characterization, the missing of some transient pulses will also not affect the estimation of λ 2 . Because the estimated λ2 (the expected number of Fourier 77/129 windows over transient state) is only related to τ n (2) according to Eq. 3.26. Therefore, hidden state incompleteness will not affect the estimation of λ 2 . But overcompleteness (the noise state is estimated wrongly as transient state) will affect the estimation. Such influence will be decreased after the linear regression. In conclusion, wrong estimation of a state or the state sequence incompleteness affect little the finial result. x(t). The composite signal is shown in Figure 4.9 (b), in which, the fault signal x(t) is totally covered by the noise. However, in Figure 4.10, the PPS of the composite signal displays the obvious information of signal x(t). This example proves the proposed method is immune to noise that locates in different resonance frequency band. In another case, when they have the same resonance frequency, the effectiveness of the method depends on the magnitude of the noise. If the signal of interest x(t) is totally covered by the noise in the same resonance frequency band, for example, the strong noise from the same mechanical system, the proposed method will fail.

Variable condition data

In order to illustrate the benefit of the integrated auto-diagnostic framework, experimental data from the Case Western Reserve University Bearing Data Center [START_REF]Case Western Reserve University Bearing Data Center Website[END_REF] are examined. In Cases 1 and 2, 24 sets of data associated to the drive end with three bearing fault sizes are selected: 12 sets for the inner race fault and 12 sets for the outer race fault. There are four different motor speeds and motor loads in each fault size. In Case 3, the rolling element fault data are from Technology Research Group of the Polytechnic University of Madrid [START_REF] César | Bearing Database[END_REF]. The dataset comprises vibration signals from a double-row spherical roller bearing (FAG22205E1KC3) with four different fault depths and three speeds of rotation.

Case 1: Inner race fault data

In this case, 12 sets of inner race fault data acquired from drive end are used, and the sampling frequency is 48 kHz. Table 5.1 gives a summary of the data number and the

p 1 = p(f BP F I |t 1:S ), p 2 = p(f BP F O |t 1:S ), p 3 = p(f BSF |t 1:S ), p 4 = p(f F T F |t 1:S ).
From the results, it is seen that the inner race fault probability, p 1 , is maximum for almost all the records, proving the effectiveness of the statistical analysis for automated identification method. In the "Characterization" column, the fault size quantification results are given according to the algorithm of Chapter 3.5. Note that the results for the record #174 are meaningless because this dataset is likely to have been corrupted due to some unknown reasons, e.g. mechanical looseness as suggested in the reference [START_REF] Smith | Rolling element bearing diagnostics using the Case Western Reserve University data: A benchmark study[END_REF]. The assessment of the fault size is illustrated in 5.4. The corrected duration time of the fault is eventually obtained by subtracting the coda error ∆e, i.e. ∆t -∆e, with ∆t calculated from Eq.(4.9). Figure 5.1 shows that the estimate of the fault size varies slightly with motor speed, but is close to the actual size within an acceptable error. In this figure, the estimated fault sizes are in the 90% confidence intervals of the actual values obtained by assuming a uniform error with standard deviation equal to R/F s . In order to assess its performance for bearing fault identification, the PPS S(k) 85/129 Envelope spectrum of bandpass signal In order to illustrate the sensitivity of the proposed method with respect to the window length, the effect of taking different values of N w is investigated when estimating the empirical posterior probability p 1 -see Figure 5.6. It is seen that only the case with N w = 128 has probability greater than 0.5, yet the probability remains significant for at least three octaves of N w . This demonstrates that the temporal structure of the signal can 87/129 be revealed in a reasonably large range of values of the hyperparameter N w . In addition, it is found that the window shift R almost has no effect on the result.

84/129

f BPFI f BPFI 2 f BPFI 3 f BPFI 3 f BPFI 2 f BPFI

Case 2: Outer race fault data

In this case, the outer race fault datasets from the Case Western Reserve University are used to validate the integrated diagnostic method. Table 5.5 lists the basic information of the used 12 data sets. The sampling frequency is 48 kHz.

The same STFT parameters are used, so the quantile in the fault detection is also same with the first case, i.e. χ 2 v,1-α = 3253.1. For record #201, the presence of a fault is correctly detected, but the probability of outer race fault is not the maximum. Records #202 -#204 are detected wrongly as healthy signal. It is noted here that these signals were already recognized as difficult to diagnose in the reference [START_REF] Smith | Rolling element bearing diagnostics using the Case Western Reserve University data: A benchmark study[END_REF]. A visual inspection of these signals actually confirms that they strongly resemble healthy signals. They have been analyzed here to fairly assess the limits of the proposed method. It is emphasized that the latter has not been designed to improve 89/129 the detection capability of state-of-the-art methods (e.g. as given in [START_REF] Smith | Rolling element bearing diagnostics using the Case Western Reserve University data: A benchmark study[END_REF]), but rather to cover the four consecutive tasks of diagnosis in an automated way. The fact that it does so nearly as well as finely hand-tuned methods is remarkable. There are four different fault depths and three speeds of rotation. The sampling frequency is 40 kHz. Each combination was recorded three times during 30 s to ensure that small differences due to uncontrollable variables were distributed evenly across all records. Table 5.7 reports the rolling element fault data sets and their corresponding basic information, and Table 5.8 lists the potential fault frequencies. More details can be found in reference [START_REF] Ricardo | Low-Cost, High-Frequency, Data Acquisition System for Condition Monitoring of Rotating Machinery through Vibration Analysis-Case Study[END_REF]. This case aims at analyzing all signals recorded by the accelerator at 6 o'clock on the bearing casing. In this case, it is less efficient than the PPS. The reason is that the signal contains one very strong impulse, much higher than the other ones (an observation quite typical of ball faults), which jeopardizes its spectral analysis (note that this could be fixed on an expert level by using the robust envelope analysis technique described in Ref. [START_REF] Borghesani | CS2 analysis in presence of non-Gaussian background noise -Effect on traditional estimators and resilience of log-envelope indicators[END_REF]). By comparison, the PPS is robust to outliers.

Similarly, the posterior probabilities p(f |t 1:S ) of the fault type displayed in Table 5.9. It is seen that p(f BSF |t 1:S ) is the highest, thus pointing out that the fault comes from the rolling element. One advantage of the latter is that it performed automatically, without visual examination.

The zoomed output of the EDHMM is shown in Envelope spectrum of bandpass signal in Ref. [START_REF] Ricardo | Low-Cost, High-Frequency, Data Acquisition System for Condition Monitoring of Rotating Machinery through Vibration Analysis-Case Study[END_REF], the fault covers an area of 11.05 mm 2 and extends over the full length of the roller, that is 6.5 mm, the fault width is roughly 1.70 mm. It is found that the estimated value is not too far from the actual value. Now, the capability of the proposed method to reconstruct the fault signal is validated on the record #12.

2 f BSF f BSF 3 f BSF f BSF 3 f BSF

Accelerated degradation data

The previous section has focused on the failure resulting from an artificial damage.

In this section, the life time dataset of rolling element bearings is used to gauge the performance of the proposed method in detecting a growing potential failure of a new 94/129 of failure: 10%). According to [START_REF]Rolling Bearings-Dynamic Load Ratings and Rating Life[END_REF], the theoretical expected fatigue life can be calculated using the following equation,

L h = 10 6 60n C P , (5.1)
where, L h is the expected fatigue life (h), n is the rotating speed (rpm), C and P are the dynamic load rating (kN) and the applied dynamic load (kN) respectively. The exponent of the life equation is 3 for a ball bearing. The actual life time tends to be longer than the rating life. For this example, the expected life time can be calculated as 14.6 h, while the the actual situation is 42. From the time the fault is found in the vibration signal to the end of the degradation test, the fault size increases by about a factor 1.6.

Compared to the traditional time domain indicators (Kurtosis, Skewness, RMS, crest, etc.), the advantage of the GLR Λ is that it follows the chi-squared distribution, which provides a theoretical quantile that can be used to make a decision. This is extremely important in engineering applications. In addition, in the proposed framework, the stochastic model also provides the mean transient pulse duration time, λ 2 , which is able to track the development of faults in a proportional manner. This enables to characterize quantitatively the fault size, not just a qualitative analysis like the traditional time domain indicators.

After the fault occurrence point is detected, the subsequent signal can be saved and analyzed in details. Figure 5.15 displays the signal at 39.6 h, which shows bearing failure occurrence. Since the fault is incipient, the impulses are still weak, and strongly masked by background noise. This signal has a kurtosis level of 3.16, which complies with a Gaussian characteristic. This shows that this signal is challenging and difficult to diagnose. The identification approach of the proposed integrated framework is applied to identify the 99/129 In order to confirm this statistical analysis, the posterior probability spectrum is shown in Figure 5.17 At last, the capability of the proposed method to reconstruct the fault signal is validated on the signal at the fault occurrence point. In the accelerated degradation data, the failure is natural, not artificial like the first experimental validation. On this set of data, it is important to detect an incipient fault in 103/129 advance before it continues to develop. After validation, it is found that the technique is effective for detecting the fault occurrence time, and is also capable of providing information on its propagation and tracking quantitatively the evolution of the crack.

Conclusion

104/129

Conclusions and perspectives Conclusion

The main objective of this thesis is the investigation of automatic signal processing

Perspectives

So far, it is time to close this thesis. The previous chapters described how to model the signal, how to estimate the model parameters, how to use the parameters for integrated diagnosis and the performance in the experimental validation. Although everything seems perfect, there are several aspects that have not been considered in the thesis. In this part, some of the future work directions related to this research are suggested.

• Simplicity and reliability

One of the most crucial aspects that should be kept in mind is that the ultimate goal of scientific research is to find practical application. In the practical CBM, the reliability of the technique rather than its sophistication is the primary pursuit. It is often required that part or even all of the proposed technique should be suitable for the slave computer or microcomputer, which means the method should be performed in a simple and automated way. Therefore, computational optimization of this integrated diagnosis framework is the next step to be done.

• Ability for variable speed conditions

In addition, in most industrial environments, the speed is often an uncontrollable variable, e.g. as with wind turbines. However, the current framework is only valid at constant speed. Angle-time analysis [START_REF] Abboud | Angle time cyclostationarity for the analysis of rolling element bearing vibrations[END_REF][12] could potentially be jointed with the Maximization with respect to π is achieved using appropriate Lagrange multipliers. First, setting the derivatives of the equation (A.3) with respect to π i to zero, together with the constrain i π i = 1, we can obtain γ 1 (i) π i + ρ = 0. Joint with the the constrain i π i = 1, the formula can be written as i γ 1 (i) = -ρ. Then the parameter π is given by,

π i = γ 1 (i) i γ 1 (i) . (A.4)
We take the derivative with respect to the parameter C i as follow, Spécialité : Génie Mécanique RESUME : Aujourd'hui, le problème le plus fondamental de la surveillance d'état dans la plupart des installations industrielles est le diagnostic et le pronostic des défauts. L'une des approches les plus efficaces pour étudier ce problème est la surveillance de l'état basée sur l'analyse des signaux de vibration. Avec le développement de l'industrie, la maintenance multithread et l'acquisition multicanal se généralisent, ce qui met en avant des exigences de maintenance plus élevées. Sur la base de cette observation, il est proposé dans cette thèse un cadre de diagnostic automatisé pour le roulement qui intègre les étapes successives de détection de défaut, d'identification du type de défaut, de reconstruction du signal de défaut et de caractérisation de la taille du défaut. L'avantage est que le processus de diagnostic complet est réalisé en une seule fois, tout en impliquant un seul hyperparamètre clé, ce qui améliore le degré d'automatisation de la maintenance basée sur les conditions (CBM) actuelle et libère la participation humaine.

∂Q(θ, θ old ) = - N n=1 γ n (i) ∂[ln |C i | + Y † n C -1 i Y n ] ∂C i = - N n=1 γ n (i) tr(C -1 i ∂C i ) -Y † n C -1 i ∂C i C -1 i Y n = - N n=1 γ n (i) tr(C -1 i ∂C i ) -tr(Y † n C -1 i ∂C i C -1 i Y ) = - N n=1 γ n (i) tr(C -1 i ∂C i ) -tr(C -1 i Y Y † C -1 i ∂C i ) = - N n=1 γ n (i) tr([C -1 i -C -1 k Y Y † C -1 i ]∂C i ) =⇒ ∂Q(θ, θ old ) ∂C i = - N n=1 γ n (i)[C -1 i -C -1 i Y Y † C -1 i ] = 0 =⇒ C i = N n=1 γ n (i)Y Y †

I

  would like to initially express my deepest gratitude to China Scholarship Council, which provides me sufficient supports to complete my Ph.D project (No.CSC201806070149) in Laboratory of vibration and acoustics (LVA) of INSA-Lyon from 2018 to 2022. The wonderful experience in France will last forever in my heart. During the Ph.D research, I believe I have gradually grasped how to become a good researcher and gained invaluable understandings about myself. This growth in my research and mindset cannot happen without the help of the professors, colleagues and friends.

4 . 5 1 . 1

 4511 de la littérature sur cette question montre que les techniques existantes se concentrent au contraire sur des tâches diagnostique spécifiques. De plus ces techniques sont presque toutes indépendantes et nécessitent que les utilisateurs aient un niveau élevé d'expérience et de connaissances. Certaines méthodes tentent de répondre à une telle aspiration, mais n'atteignent pas la dose automatisée et intégrée comme dans cette thèse. La méthode proposée surmonte ces lacunes comme démontré sur les données expérimentales du chapitre xi Contents 2.3.1.2 Stochastic model . . . . . . . . . . . . . . . . . . . . . . . 2.3.2 Bearing fault signal enhancement . . . . . . . . . . . . . . . . . . . 2.3.2.1 Spectral Subtraction . . . . . . . . . . . . . . . . . . . . . 2.3.2.2 Wiener filtering . . . . . . . . . . . . . . . . . . . . . . . . 2.3.2.3 Minimum Mean-Square Error . . . . . . . . . . . . . . . . 2.3.3 Bearing fault identification . . . . . . . . . . . . . . . . . . . . . . . 2.3.4 Bearing fault size characterization . . . . . . . . . . . . . . . . . . . 2.4 Hidden Markov model and its variants . . . . . . . . . . . . . . . . . . . . 2.4.1 Hidden Markov model . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.2 Explicit duration hidden Markov model . . . . . . . . . . . . . . . . 2.5 Integrated automatic diagnosis . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Markov stochastic process 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Markov chains . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Hidden Markov model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Vibration signal modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Signal model and STFT decomposition . . . . . . . . . . . . . . . . 3.4.2 Explicit duration hidden Markov model . . . . . . . . . . . . . . . . 3.5 Parameter estimation and analysis . . . . . . . . . . . . . . . . . . . . . . 3.5.1 Parameter estimation . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5.2 Parameter analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Input parameter settings . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Integrated auto-diagnostic framework 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Fault detection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Fault signal reconstruction . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Fault type identification . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.1 Posterior probability spectrum . . . . . . . . . . . . . . . . . . . . . 4.4.2 Statistical analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv Contents Fault size characterization . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Experimental validation 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Variable condition data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.1 Case 1: Inner race fault data . . . . . . . . . . . . . . . . . . . . . . 5.2.2 Case 2: Outer race fault data . . . . . . . . . . . . . . . . . . . . . 5.2.3 Case 3: Rolling element fault data . . . . . . . . . . . . . . . . . . . 5.3 Accelerated degradation data . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Intelligence maintenance framework. . . . . . . . . . . . . . . . . . . . . . 1.2 Components of the integrated diagnosis framework. . . . . . . . . . . . . . 2.1 The diagrammatic sketch of a rolling bearing. . . . . . . . . . . . . . . . . 2.2 Vibration signal analysis techniques. . . . . . . . . . . . . . . . . . . . . . 2.3 Illustration of HMM-based signal analysis. . . . . . . . . . . . . . . . . . . 3.1 A first-order Markov chain of observation {Y n }, in which the distribution p(Y n |Y n-1 ) of a particular observation Y n is conditioned on the value of the previous observation Y n-1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 The graphical structure of a hidden Markov model, in which each observation Y n is condition on the state of the corresponding hidden state and the state sequence is a first order Markov chain. . . . . . . . . . . . . . . . . . . . . 3.3 Example of a binary state Markov chain in a noisy signal. The continuous signal is segmented by short windows and conditioned on the corresponding hidden state. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Explicit duration hidden Markov model, in which one state emits d i values of observations, and each observation follows a complex-valued Gaussian distribution, i, j ∈ {1, 2}. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Illustration of Poisson distribution. The upper subplot is the PMF with the horizontal axis as the duration time τ n (i), i.e. the number of remainings in current state i. λ i is the expected duration time in state i; The lower subplot is the CDF showing discontinuous at the integers of τ n (i) because a variable that is Poisson distributed takes on only integer values. . . . . . . xvii 3.6 Illustration of the difference between a dynamic and constant transition matrix. (a) The EDHMM state sequence and (b) its Fourier spectrum; (c) the corresponding HMM state sequence and (d) its Fourier spectrum. . . . 51 3.7 The synthetic signal with SNR=0 dB and sampling frequency f s = 4.8kHz over one second. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.8 The spectrogram of the synthetic signal with resonance frequency f 0 = 5kHz and fault frequency f BP F I = 161Hz. . . . . . . . . . . . . . . . . . . . . . 58 3.9 Zoom of simulated signal in interval [0 0.1] s. (a) Raw vibration signal; (b) posterior probability γ n (2); (c) duration time sequence τ n (2) of the active state; (d) duration time sequence τ n (1) of the inactive state. . . . . . . . . 59 3.10 Illustration of the discrepancy of the parameters λ in different data sets. . 60 4.1 Illustration of fault detection based on chi-squared distribution with parameter v = 10 and α = 0.05. . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 4.2 Illustration of the capability of fault signal reconstruction technique. (a) Simulated signal with 1800 RPM; (b) Corrupted signal with SNR= 0 dB Gaussian noise; (c) Reconstructed fault signal. . . . . . . . . . . . . . . . . 69 4.3 Illustration of the capability of fault signal reconstruction technique. (a) Simulated signal with 1800 RPM; (b) Corrupted signal with SNR= -5 dB Gaussian noise; (c) Reconstructed fault signal. . . . . . . . . . . . . . . . . 70 4.4 Illustration of the capability of PPS. (a) Posterior probability spectrum S(k); (b) Envelope spectrum of original vibration signal. . . . . . . . . . . 72 4.5 Illustration of the impulse cycles. The first sequence is the state posterior probability; second is the vibration signal; The last two subplots are the duration time sequence for different states. . . . . . . . . . . . . . . . . . . 73 4.6 The illustration of fault size characterization through linear regression, with the horizontal axis as the function c(f j r ) and vertical axis as the duration time ∆t j , in which the estimated fault size l and the coda error ∆e correspond to the slope and the intercept of the fitting line. . . . . . . 75 4.7 The flowchart of the integrated diagnosis framework. . . . . . . . . . . . . 77 4.8 Histogram of the transient pulse cycle frequency, showing the red dash line relating to the fault characteristic frequency located in the highest bin. . . 78 4.9 Illustration of the synthetic signal. (a) The noise n(t) and the signal of interest x(t); (b) the composite signal. . . . . . . . . . . . . . . . . . . . . 79 xviii 4.10 The posterior probability spectrum of the composite signal. . . . . . . . . . 5.1 Quantification of the fault size together with 90% confidence intervals based on the assumption of a uniform error with standard deviation equal to R/F s . 5.2 Spectrogram of the record #175 and the corresponding covariance of two states, C 1 and C 2 , in the right subplot. . . . . . . . . . . . . . . . . . . . . 5.3 Spectral analysis of record #175. (a) Posterior probability spectrum S(k); (b) Envelope spectrum of bandpass signal with maximum kurtosis in band [5; 6] kHz selected from kurtogram. . . . . . . . . . . . . . . . . . . . . . . 5.4 Zoom of #175 in interval [0.5, 0.6] s. (a) Raw vibration signal; (b) the corresponding spectrogram; (c) posterior probability γ n (2); (d) duration time sequence τ n (2) of the active state; (e) duration time sequence τ n (1) of the inactive state. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Fault signal reconstruction. (a) Raw vibration signal #175; (b) Reconstructed fault signal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Empirical posterior probability p 1 for different values of the hyperparameter N w . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Quantification of the fault size together with 90% confidence intervals based on the assumption of a uniform error with standard deviation equal to R/F s . 5.8 Spectral analysis of dataset # 12. (a) Posterior probability spectrum S(k); (b) Envelope spectrum of bandpass signal with maximum kurtosis in band [8.8; 10] kHz selected from kurtogram. . . . . . . . . . . . . . . . . . . . . 5.9 Zoom of #12 in interval [0.5 0.6] s. (a) Vibration signal; (b) posterior probability γ n (2); (c) duration time sequence τ n (2) of active state; (d) duration time sequence τ n (1) of inactive state. . . . . . . . . . . . . . . . . 5.10 Fault signal reconstruction. (a) Raw vibration signal #12; (b) Reconstructed fault signal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.11 Signal acquisition settings. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.12 Full life vibration signal in the horizontal direction of the bearing 3_1. . . 5.13 Bearing 3_1 vibration signal indicators; RMS: Root mean square; KU: Kurtosis; CI: Crest indicator; SI: Shape Indicator. . . . . . . . . . . . . . . xix 5.14 Illustration of fault detection and fault development tracking. (a) accelerated degradation vibration signal; (b) interval of fault occurrence; (c) generalized likelihood ratio Λ for detecting fault occurrence; (d) Poisson distribution parameter λ 2 for tracking the fault development. . . . . . . . . . . . . . . . 100 5.15 Signal with the incipient fault @ 39.6 h. . . . . . . . . . . . . . . . . . . . 101 5.16 Histogram of the transient pulse cycle frequency, showing the red dash line relating to the fault characteristic frequency located in the highest bin. . . 101 5.17 Spectral analysis of the signal at 39.6 h. (a) Posterior probability spectrum S(k); (b) Envelope spectrum of bandpass signal with maximum kurtosis in band [3.2; 4.3] kHz selected from kurtogram. . . . . . . . . . . . . . . . . . 102 5.18 Fault signal reconstruction. (a) Raw vibration signal at 39.6 h; (b) Reconstructed fault signal. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 xx Nomenclature Acronyms and abbreviations

  Scheduled maintenance refers to schedule the maintenance plan in advance through the prior knowledge about the average service life of the crucial component, and intervene in regularly to the system based on the plan. So it implies taking the risk of replacing the healthy operational components. With the development of science and technology, it makes the Predictive Maintenance possible. This kind of maintenance can be achieved through a series of tests and analyses without the occurrence of damage. It is arranged and implemented based on the operating status of the equipment, and is called Condition

Figure 1 . 1 :

 11 Figure 1.1: Intelligence maintenance framework.

  of bearing intelligence CBM and proposed one integrated auto-diagnosis framework including fault detection, fault signal reconstruction, fault type identification and fault size characterization, which are as shown in Figure 1.2. The main contributions of the present work are resumed as follows. It introduces a complete methodology that:

Figure 1 . 2 :

 12 Figure 1.2: Components of the integrated diagnosis framework.

Chapter 3

 3 This thesis is divided into six chapters. The contributions and contents of each chapter are explained separately as follows:Chapter 1 explains the background and motivations behind this research, as well as the research questions. It also describes the main objectives and contributions of this research and outlines an overview of the dissertation. Chapter 2 outlines the state of the art on the bearing fault diagnosis. Firstly, the characteristics of rolling element bearing are described in Section 2.2, including basic structure, vibration characteristics and incipient fault characteristics. Section 2.3 describes what current techniques are applied in bearing fault diagnostics from different viewpoints, including detection, fault signal enhancement , identification and the fault size estimation techniques. Literature on stochastic modeling method are reviewed, especially the hidden Markov model and explicit duration hidden Markov model. In the final, the automatic diagnosis techniques are also reviewed. The chapter closes with one discussion based on the shortcomings of existing techniques, thereby leading to the research content of this dissertation. considers the stochastic models, and how to use them to model the timefrequency distribution of the bearing vibration signal, including the parameter setting, estimation and analysis. Section 3.2 introduces firstly the Markov chain, paving the way for the following HMM and EDHMM in Section 3.3 and Section 3.4, respectively. The proposed framework mainly depends on the EDHMM. It aims at capturing the time evolution of the STFT coefficients and computing the duration time in different states.

6 . 4 129 1. 3 .

 641293 Chapter presents the integrated diagnosis framework including fault detection, fault signal reconstruction, fault type identification and fault size quantification. This chapter 6/Organization of the thesis is the core work of this thesis. The EDHMM parameters obtained in the previous chapter will be used in this proposed diagnosis framework. Section 4.2 introduces the likelihood ratio test (LRT) and provides the EDHMM parameters for detecting the presence of fault without requiring any data pre-preprocessing. A time-varying filter based on Bayesian theory is proposed to extract the fault signal in full-band in Section 4.3. Section 4.4

Chapter 5

 5 presents the validation of the integrated diagnosis framework through different bearing datasets. Two different experimental scenarios are reported in this paper. In Section 5.1, the vibration signal relative to bearings with different damages, running at different speeds and under different loads are used. The second scenario (Section 5.2) addresses the accelerated degradation data. In the validation, the proposed method is also compared with the state-of-the-art kurtogram. In the last section, results of the different experimental datasets are discussed. The last part completes the thesis with general conclusions about the proposed integrated auto-diagnosis framework. Ideas for future research challenges are discussed together with suggestions for improvements of the proposed techniques. The Appendix gives some mathematical deduction used in this thesis, including the detailed process of the stochastic model parameter estimation in Appendix A, and the algorithm convergence in Appendix B. 7/129

Figure 2 . 1 :

 21 Figure 2.1: The diagrammatic sketch of a rolling bearing.

  estimation of features or condition indicators from sensor data. The key attribute of the features is that they characterize a large number of targeted faults. Their selection and extraction are critical processes that affect unequivocally the success or failure of the detector. Another critical component is the fault progression model describing the degrading state of the system. In reality, however, this model is often unknown and severely hinders the application of model-based detection method. With the absence of data, it is impossible to build a data-driven model. The authors [168] therefore developed a model based on the Paris fatigue law and modified it to adapt to the fault mode of interest 17/129

  of background noise due to adverse recording environments, as well as convolutional noise due to sensor and signal propagation variability. Various signal processing techniques involving time, frequency and statistical methods have been used to detect incipient faults. These techniques require a high signal-to-noise-ratio (SNR), where the faulty component vibrations are higher than the background noise. Standard spectral features are highly sensitive to noise, which can decrease the following modeling effectiveness and give misleading results. It is important, therefore, to incorporate signal enhancement and/or robust spectral extraction techniques that enhance the feature while suppressing background noise. The basic goal of feature enhancement is to extract the clean signal x(t) from the measurement y(t) with additive background noise n(t). There are various ways to denoise in the pretreatment, and this part mainly concentrates here on the class of denoising methods that capitalize on time-frequency feature. The signal model can be written as y(t) = x(t) + n(t), and the corresponding Fourier transform as Y (f ) = X(f ) + N (f ).

  techniques using the data collected from the bearing. During the past few decades, a significant body of research has been carried out for addressing directly or indirectly fault identification. The state-of-the-art signal processing techniques used recently are listed in Figure 2.2, mainly including 1) signal decomposition: empirical mode decomposition and local mean decomposition, blind deconvolution, matrix decomposition techniques;2) signal modeling: cyclostationary, hidden Markov model; 3) signal transformation: STFT, Wigner-Ville decomposition, wavelet transformation, cepstral analysis; 4) artificial intelligence, etc. In the final analysis, these advanced processing techniques all need to return back to demodulation and spectrum analysis for finding the bearing fault frequency so as to identify the fault type. These methods are reviewed as follow.

Figure 2 . 2 :

 22 Figure 2.2: Vibration signal analysis techniques.

  explores the cyclostationary characteristics of the defective bearing signals. The vibration signals of rolling bearings consist of random and periodic components. The autocorrelation function of these signals exhibits time varying, periodic and cyclostationary temperament. Sawalhi and Randall[START_REF] Sawalhi | Simulating gear and bearing interactions in the presence of faults: Part I. The combined gear bearing dynamic model and the simulation of localised bearing faults[END_REF][START_REF] Sawalhi | Simulating gear and bearing interactions in the presence of faults: Part II: Simulation of the vibrations produced by extended bearing faults[END_REF] studied the extended inner race and outer race fault in rolling element bearings in the presence of gear interaction. The cyclostationary properties such as spectral correlation function of the system were utilized to separate out the bearing faults from the gear faults. Antoni[START_REF] Antoni | Cyclostationarity by examples[END_REF] formulated systematically and comprehensively the theory of cyclostationarity, including the definition, cyclostationary tools and the application. Chen[START_REF] Dong | Noise resistant time frequency analysis and application in fault diagnosis of rolling element bearings[END_REF] used the cyclic spectral density in fault diagnosis of rolling element bearings; this method is described as the Fourier transform of the cyclic autocorrelation function obtained by modifying the time varying autocorrelation function of a Wigner-Ville distribution based on the assumption that the bearing signals possess a cyclostationary character. Abboud et al. [2] extended cyclostationary theory to nonstationary operating conditions through the interaction between time and angle, which is materialized by the angle-periodicity of the correlation 24/129 measure of two versions of the signals shifted by a constant time-lag. Cepstrum is one of the oldest techniques used in the field of signal processing and has been improved by several authors to apply it to the rolling element bearings fault investigation. Cepstrum and logarithmic power spectrum are a pair of Fourier transforms, put differently, cepstrum is the inverse Fourier transform of the power spectrum in logarithmic coordinates. Local faults in gears give an impulsive modulation of the gearmesh signals (both amplitude and frequency modulation) resulting in large numbers

27 / 129 can

 27129 alternatively allow for sufficient lead time for maintenance planning as well. The current research are almost based on the double impact detection, and the time between them is used to estimated fault size.The earliest research of bearing fault characterization and quantification using vibration analysis was by Epps[START_REF] Epps | An investigation into vibrations excited by discrete faults in rolling element bearings[END_REF]. He mainly investigated the relationship between the waveform of the vibration source and the localized defect. As he stated in the reference, the vibration characteristics associated with ball passage over a raceway fault could consist of a step response for the ball entry into the fault and an impulse response for the ball impact on the trailing edge of the fault. The entry into the fault appears like a step response, with mainly low frequency content, while the impact on exit excites a much broader band impulse response. Based on this observation, Sawalhi and Randall[START_REF] Sawalhi | Vibration response of spalled rolling element bearings: Observations, simulations and signal processing techniques to track the spall size[END_REF] had a further investigation with both simulated and seeded faults. In order to enable a clear separation of the two events, and produce an averaged estimate of the size of the fault, some advanced signal processing techniques, such as pre-whitening filtering, wavelet and cepstrum analyses, are used to enhance the entry event while keeping the impulse response for fault size estimation. In reference[START_REF] Wang | Size estimation for naturally occurring bearing faults using synchronous averaging of vibration signals[END_REF], Wang et al. presented a technique based on synchronous signal averaging with the bearing fault characteristic frequency obtained from raw vibration signatures. The averaged signal represents the vibration characteristics within one period of impact produced by the bearing fault. When the fault size is smaller than the diameter of the balls, the features associated with the ball's entry into and exit from the fault may be extracted and the fault size could be derived. Subsequently, in reference[START_REF] Sawalhi | Vibration signal processing for spall size estimation in rolling element bearings using autoregressive inverse filtration combined with bearing signal synchronous averaging[END_REF], the squared envelope of the bearing synchronous averaging is jointed with autoregressive inverse filtration, this combination gives a superior enhancement to the step response and balances it with the impulse response. Ahmadi et al.[START_REF] Alireza Moazen Ahmadi | The path of rolling elements in defective bearings: observations, analysis and methods to estimate spall size[END_REF] conducted some experiments, and observed the vibration characteristics under different applied load and rotating speed. The experimental results show strong dependency of the duration time of the low frequency entry and exit transient events on the applied load. It is evident that changes in the speed have almost no effect on duration time of the low frequency entry and exit transient events. Therefore, a new defect size estimation method is proposed

30 / 129 test

 30129 to achieve the transient signal detection. The frequency line detection is a similar problem studied in[START_REF] Streit | Frequency line tracking using hidden Markov models[END_REF][START_REF] Quinn | The estimation and HMM tracking of weak narrowband signals[END_REF][START_REF] Barrett | Frequency tracking using hidden Markov models with amplitude and phase information[END_REF]. Streit and Barrett[START_REF] Streit | Frequency line tracking using hidden Markov models[END_REF] and Quinn et al.[START_REF] Quinn | The estimation and HMM tracking of weak narrowband signals[END_REF] used an HMM-based algorithm to track the peak of the STFT. In the literature[START_REF] Barrett | Frequency tracking using hidden Markov models with amplitude and phase information[END_REF], the author combined the FFTs feature (amplitude and phase values) and HMM tracker as frequency line detectors. The recurrence of periodicities in the airflow trace data is also covered in the literature. The authors in[START_REF] Hadj-Amar | A Spectral Hidden Markov Model for Nonstationary Oscillatory Processes[END_REF] proposed a spectrum-based HMM where the discrete latent state sequence reflects the time-varying changes as well as recurrence of periodic regimes as defined by their spectral properties. HMM was used for fault diagnosis based on vibration signal since the 1990s. There are more than 150 paper focus on fault diagnosis and prognosis from 2000. Looking at the existing literature, the HMM-based methods for fault diagnosis can be roughly classified in two ways.

  the fault signal from the raw vibration signal. They modeled the signal pattern with different distribution in time-frequency domain as the different hidden states, and utilize the advantage of invertible mathematically of STFT to extract the time domain fault signal. After investigation, this is a novel idea for modeling the time-frequency futures of the vibration signal for diagnosis. Then they proposed a complete framework integrating detection, recognition, extraction and diagnosis in literature[START_REF] Xin | Semi-automated diagnosis of bearing faults based on a hidden Markov model of the vibration signals[END_REF].

  techniques have been proposed to address the different diagnostic tasks. However, these techniques are almost all independent, which requires users to have a high level of experience and knowledge when dealing with different diagnostic problems. Concerning the techniques that demonstrate attempts towards integrated and automated diagnosis in Chapter 2.5, there is still no one research that is able to integrate all the diagnostic sub-issues in one solution. Even fewer have tried to achieve these objectives in an automated way, i.e. without manual tuning of the algorithms by an expert. According to the ISO 13374-1:2003 at Condition monitoring and diagnostics of machines, the diagnosis process can be divided into data acquisition, data manipulation, state detection, fault identification, health assessment, prognostic assessment, and advisory generation. So all of these reasons push us in one direction, integrated auto-diagnosis technique. As we known, the ultimate goal of research is to apply in practice, which requires the method performed in a simple and automated way. This is exactly what the thesis is pursuing.Following the goal mentioned above, this thesis is dedicated to propose one stochastic model based diagnostic framework, achieving automation and integration at the same time. As mentioned in the bearing vibration characteristics, the intrinsic uncertainties and nonstationarity which underlie the incipient fault vibration signal and its temporal sequential nature, made EDHMM-based approaches a perfect option for our objectives.In addition, as an advanced model, EDHMM comes with several valuable parameters useful for clarifying the transition law and finding temporal structures in the vibration signals. However, as investigated in Chapter 2.4, the HMMs-based approaches that are 36/129 implemented for the diagnostic task in the literature are basically used as a supervised way which needs massive data to train different models in advance, and then exploiting the trained models to process the testing signal. The computational complexity and the assumption of data available significantly impede the diagnosis applications in a practical scenario. Unlike the supervised way, in this thesis, the stochastic model is used as signal analysis tool to extract the fault information for the following integrated auto-diagnosis.The modeling and analysis process is illustrated in Figure2.3. The repetitive series of transients passes a filter-bank analysis implemented by the STFT to generate observations. For such rotating vibration signals, the Fourier transform have proved good performances, but other transforms are also possible. And then under the theory of stochastic model, the useful information hidden in the vibration signal can be extracted by estimating the model parameters. In this example, there exist two different time-frequency distributions which indicate the transient state and the noisy state, respectively. In the meantime, the distribution parameters in the observation with different states can be estimated, and the duration time in different states can also be obtained. These important information will encourage us to achieve fault detection, fault signal reconstruction, fault identification and fault size characterization at once without historical data and manual intervention. Since this integrated framework does not impose too many prior knowledge and parameters, it

Figure 2 . 3 :

 23 Figure 2.3: Illustration of HMM-based signal analysis.

3 . 1 3 . 2 3 . 3 3 . 4 3 . 5 3 . 6 3 . 7

 31323334353637 diagnosis framework, and one in automatic way. It is found the integrated and automated diagnosis based on stochastic model of vibration signal is still a new and challenging topic. Hence, some comments have been discussed in the end, and the objectives of this thesis have been emphasized once again. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Markov chains . . . . . . . . . . . . . . . . . . . . . . . . . . . . Hidden Markov model . . . . . . . . . . . . . . . . . . . . . . . Vibration signal modeling . . . . . . . . . . . . . . . . . . . . . 3.4.1 Signal model and STFT decomposition . . . . . . . . . . . . . 3.4.2 Explicit duration hidden Markov model . . . . . . . . . . . . . Parameter estimation and analysis . . . . . . . . . . . . . . . . 3.5.1 Parameter estimation . . . . . . . . . . . . . . . . . . . . . . . 3.5.2 Parameter analysis . . . . . . . . . . . . . . . . . . . . . . . . . Input parameter settings . . . . . . . . . . . . . . . . . . . . . . Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3 . 1 :

 31 Figure 3.1: A first-order Markov chain of observation {Y n }, in which the distribution p(Y n |Y n-1 ) of a particular observation Y n is conditioned on the value of the previous observation Y n-1 .

••

  sponding states are finite but invisible. In such cases, each state randomly generates 1 out of every k observations visible to us. We can model this by introducing additional hidden states to permit a general case. Let the hidden state sequence {z n , n = 1, 2, 3, ...} be a Markov chain with the transition probabilities a ij and initial probabilities π i = p(z 1 = i), i ∈ M defined in the previous section. Further, suppose that when the Markov chain enters state i at time n, then the observation Y n is emitted with probability p(Y n |z n = i).Therefore, a model of the previous defined type in which the sequence {Y n , n = 1, 2, 3, ...} is observed, while the sequence of underlying Markov chain states {z n , n = 1, 2, 3, ...} are unobserved, is called a hidden Markov model (HMM). It is an extension of a Markov chain which is able to capture the sequential relations among the hidden variables. The sequential data is represented using a Markov chain of hidden states, with each observation conditioned on the corresponding hidden state, as shown in Figure3.2.According to the above introduction, a HMM can be governed by the following parameters: Number of states N This is usually set to the total number of distinct, or elementary, stochastic events in a signal process. Transition

Figure 3 . 2 :

 32 Figure 3.2: The graphical structure of a hidden Markov model, in which each observation Y n is condition on the state of the corresponding hidden state and the state sequence is a first order Markov chain.

43 / 129 isFigure 3 . 3 :

 4312933 Figure 3.3: Example of a binary state Markov chain in a noisy signal. The continuous signal is segmented by short windows and conditioned on the corresponding hidden state.

  The STFT coefficients Y (n, k) of a discrete-time signal y(t), t = 0, 1, 2, ..., L -1 over 45/129 a time interval of length N w is given byY (n, k) = Nw-1 m=0 y(nR + m)w(m)e -j2πm k Nw ,(3.6) where n = 1, 2, ..., N , N = (L -N w )/R + 1 and k = 1, 2, ..., N k , N k = N w /2 , represent the time and frequency indices, respectively; w(m) denotes a positive and smooth N w -long data window, and R represents the window shift.

  the sum of squares of real and imaginary parts, and represents the energy of the signal at the datum n with frequency k. Similiarly, the pseudo-covarianceM = E[(Y n -µ)(Y n -µ)T ] can be deduced as a zero matrix. Then, the Gaussian random variables Y n are only determined by their covariance matrix C, i.e. Y n ∼ CN (0, C). It is also called a circularly-symmetric Gaussian distribution, whose distribution, p(Y n ), is given by,

  will be a key parameter for the following diagnosis. One advantage of Poisson distribution is that it has only one parameter, and all the distribution information can be obtained through this parameter. In addition, from the characteristic of Poisson distribution, the duration time in different transient pulses is very similar, which makes it well suited for Poisson distribution. The simplified model is shown in Figure3.4. This extends the standard HMM to the domain of modeling the temporal structure of hidden state sequence.From now on, one will need to define the dynamic duration time τ n (i), the time already spent in the current state i at datum n. If the entry and exit in state i occur at time instants n 1 and n 2 , respectively, then τ n 1 -1 (i) = 0, τ n 1 +m (i) = m for m = 0, 1, ..., n 2 -n 1 , and τ n 2 (i) = n 2 -n 1 + 1 is one sample of the discrete variable d i .

Figure 3 . 4 :

 34 Figure 3.4: Explicit duration hidden Markov model, in which one state emits d i values of observations, and each observation follows a complex-valued Gaussian distribution, i, j ∈ {1, 2}.

Figure 3 .

 3 5 shows one example of Poisson distribution with the parameter λ = 5. From the CDF in the lower subplot and Eq.(3.

Figure 3 . 5 :

 35 Figure 3.5: Illustration of Poisson distribution. The upper subplot is the PMF with the horizontal axis as the duration time τ n (i), i.e. the number of remainings in current state i. λ i is the expected duration time in state i; The lower subplot is the CDF showing discontinuous at the integers of τ n (i) because a variable that is Poisson distributed takes on only integer values.

Figure 3 . 6 :

 36 Figure 3.6: Illustration of the difference between a dynamic and constant transition matrix. (a) The EDHMM state sequence and (b) its Fourier spectrum; (c) the corresponding HMM state sequence and (d) its Fourier spectrum.

  .[START_REF] Berouti | Enhancement of speech corrupted by acoustic noise[END_REF] 54/129 In the case of two-state EDHMM, one has the useful relationships π1 + π2 = 1 and Ĉ1 + Ĉ2 =

1 ) 2 )

 12 the condition (i = j) means the weighted average of τ n (i) is calculated only at the transition time. Eqs.(3.23) to (3.26) complete one updating process for estimating the parameters θ = {π 1 , λ, C}, which is then iterated. The iteration stop criterion depends on the increase rate of the likelihood function p(Y 1:N |θ), shown in Eq.(3.21). The proof of this iterative convergence is demonstrated in Appendix B. The full estimation procedure is summarized as follows, Assume an initial parameter set θ old . Compute the forward-backward variables α n (i) and β n (i) through Eqs.(3.15) and (3.16). Next, calculate the forward-backward variables and the related probability γ n (i) and ξ n (i, j) through Eqs.(3.18) and (3.22). Finally, etimate the new parameter set θ new through Eqs.(3.23) to (3.26).

Figure 3 . 7 :

 37 Figure 3.7: The synthetic signal with SNR=0 dB and sampling frequency f s = 4.8kHz over one second.

(Figure 3 . 8 :

 38 Figure 3.8: The spectrogram of the synthetic signal with resonance frequency f 0 = 5kHz and fault frequency f BP F I = 161Hz.

Figure 3 .Figure 3 . 9 :

 339 Figure 3.9: Zoom of simulated signal in interval [0 0.1] s. (a) Raw vibration signal; (b) posterior probability γ n (2); (c) duration time sequence τ n (2) of the active state; (d) duration time sequence τ n (1) of the inactive state.

Figure 3 . 10 :

 310 Figure 3.10: Illustration of the discrepancy of the parameters λ in different data sets.

2 4 . 3 4 . 4 4 . 5 4 . 7 129 4. 1 Introduction

 2434445471291 chapter introduced and reviewed firstly the basic theory of Markov model, including Markov chain and hidden Markov model. It then mainly introduced the used stochastic model, explicit duration hidden Markov model, for modeling the vibration 61/129 signal, from vibration signal segmentation, stochastic modeling, parameter estimation to parameter analysis.First, the STFT was formulated as a decomposition tool and it was explained how to apply it to the time-frequency transform of vibration signals, and how the observed time-frequency coefficients are modeled. Then, the iterative estimation algorithm was used to find the optimal parameter θ to maximize the likelihood p(Y 1:N |θ). Afterward, some valuable parameters provided by EDHMM were introduced, which are useful for the vibration signal diagnosis, such as the posterior probability sequence γ 1:N (i), covariance Ĉ, duration time sequence τ 1:N , Poisson parameter λ, etc. The effectiveness of these parameters and what characteristics they reveal were illustrated through a synthetic bearing fault signal. For example, the elements in the covariance Ĉ measure the energy distribution in the frequency domain, which will be a good parameter for fault detection; and the probability sequence γ 1:N (i) removes the influence of other high amplitude components like resonances and amplitude modulation; it is for some occasions able to increase the effect of spectrum analysis; the parameter λ in the Poisson distribution estimates the expected value of the duration time of different states, which is able to assess the fault severity. Finally, the hyperparameter of this model and the parameters initialization were discussed. It is found that only one hyperparameter, i.e. window length N w , needs to be set in advance. This makes it potentially suited for practical applications in the absence of prior knowledge. Fault detection . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fault signal reconstruction . . . . . . . . . . . . . . . . . . . . . Fault type identification . . . . . . . . . . . . . . . . . . . . . . 4.4.1 Posterior probability spectrum . . . . . . . . . . . . . . . . . . 4.4.2 Statistical analysis . . . . . . . . . . . . . . . . . . . . . . . . . Fault size characterization . . . . . . . . . . . . . . . . . . . . . 4.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63/In modern applications, the vibrations are often acquired in continuous flow on multiple channels and machines simultaneously, which brings the difficulty to process and analyze them. So it is not practical or possible anymore to have an expert dedicated to only looking at analysis results and condition indicators continuously. Such practical problems demand for the development of integrated and automatic diagnosis techniques. In this Chapter, to achieve this goal, an automated diagnosis framework is proposed, including detection, identification, fault size characterization and fault signal reconstruction. The estimated parameters of the EDHMM model are used to complete these tasks without need of other prior information and user's expertise. This Chapter is organized as follows. Section 4.2 describes the statistical hypothetical test algorithm called likelihood ratio test for detecting the vibration signal. In Section 4.3, a time-varying filter is designed to recover the fault signal from corrupted signal.

3 . 1 ,

 31 the current techniques, either cannot accurately find the failure occurrence point as there is no corresponding alarm threshold, or take up a lot of resources for data training. In this section, one technique based on the likelihood ratio test (LRT) is introduced. According to Neyman-Pearson, the LRT is the most powerful test, which means it has the largest possible test power (probability of detection) given the significance level α (probability of false alarm). Unlike other 64/129

Figure 4 .

 4 [START_REF]Rolling Bearings-Dynamic Load Ratings and Rating Life[END_REF] shows one example of Chi-squared distribution with the freedom degree v = 5 and alarm α = 0.05. Therefore, any confidence interval can be obtained based on distribution with known parameters, which is able to provide a mathematical basis for fault detection.

Figure 4 . 1 :

 41 Figure 4.1: Illustration of fault detection based on chi-squared distribution with parameter v = 10 and α = 0.05.

Figure 4 .

 4 2(a) depicts a simulated bearing inner race fault signal at rotating speed of 30 Hz, and the simulated fault characteristic frequency f BP F I = 161 Hz. The simulated signal was further corrupted by Gaussian noise with SNR= 0 dB which is shown in Figure 4.2(b).

Figure 4 .

 4 2(c) shows that this technique has the capability of removing Gaussian noise with a 0 dB SNR, recovering a clean signal as close as possible to the simulated signal (bearing fault signal). The same simulated signal is corrupted by the Gaussian noise with SNR=-5 dB, shown in Figure 4.3(b).

Figure 4 .

 4 3(c) shows that this technique has the capability of removing Gaussian with a -5 dB SNR and recovering a clean signal as close as possible to the source signal. The bottom subplot shows a consistent impulsive signal of a damaged bearing with an average time interval of 6.23 ms very close to the characteristic defect.The kurtosis of the corrupted signal was found to be 4.19, which is very close to a Gaussian distribution; while the kurtosis of the reconstructed fault signal was 20.73 and it can be seen that the recovered signal was improved.

Figure 4 . 2 :

 42 Figure 4.2: Illustration of the capability of fault signal reconstruction technique. (a) Simulated signal with 1800 RPM; (b) Corrupted signal with SNR= 0 dB Gaussian noise; (c) Reconstructed fault signal.

  fault type. However, apart from evaluation of vibration signal spectrum to identify specific bearing component frequencies, vibration signals can also be analysed by other means to identify the fault type. In this section, two different approaches are proposed based on the EDHMM parameters, i.e, posterior probability spectrum and statistical analysis-based identification.

Figure 4 . 3 :

 43 Figure 4.3: Illustration of the capability of fault signal reconstruction technique. (a) Simulated signal with 1800 RPM; (b) Corrupted signal with SNR= -5 dB Gaussian noise; (c) Reconstructed fault signal.

γ n ( 2 )

 2 = p(z n = 2|Y 1:N , θ), with respect to time datum n. This defines the posterior 70/129 probability spectrum (PPS) S(k) as

Figure 4 . 4 :

 44 Figure 4.4: Illustration of the capability of PPS. (a) Posterior probability spectrum S(k); (b) Envelope spectrum of original vibration signal.

Figure 4 . 5 :

 45 Figure 4.5: Illustration of the impulse cycles. The first sequence is the state posterior probability; second is the vibration signal; The last two subplots are the duration time sequence for different states.

,

  Randall made a very thorough study of the vibration signatures due to the entry and exit of the rolling elements from the fault region and the duration time between them was used to estimate the fault size. Similarly, in EDHMM the duration time of the transient pulse state can 73/129 be used as an indicator for the fault size. The parameter λ 2 of the Poisson distribution (shown in the Eq.(3.25)) represents the average number of windows that cover the transient pulse. The corresponding average duration of the transient (in second) is ∆t

Figure 4 . 6 : 129 4. 6 DiscussionFigure 4 . 7

 46129647 Figure 4.6: The illustration of fault size characterization through linear regression, withthe horizontal axis as the function c(f j r ) and vertical axis as the duration time ∆t j , in which the estimated fault size l and the coda error ∆e correspond to the slope and the intercept of the fitting line.

Figure 4 . 7 :

 47 Figure 4.7: The flowchart of the integrated diagnosis framework.

Table 4 . 1 :

 41 The posterior probabilities of fault type given t 1:S Fault type BPFI BPFO BSF FTF Posterior probability p(f BP F I |t 1:S ) p(f BP F O |t 1:S ) p(f BSF |t 1:S ) p(f F T F |t 1:S )

Figure 4 . 8 :

 48 Figure 4.8: Histogram of the transient pulse cycle frequency, showing the red dash line relating to the fault characteristic frequency located in the highest bin.

Figure 4 . 129 4

 4129 Figure 4.9 (a). It is found the amplitude of the noise is higher than the signal of interest

Figure 4 . 9 :

 49 Figure 4.9: Illustration of the synthetic signal. (a) The noise n(t) and the signal of interest x(t); (b) the composite signal.

Figure 4 . 10 : 5 . 1 5 . 2 5 . 3 5 . 4 129 5. 1 Introduction

 410515253541291 Figure 4.10: The posterior probability spectrum of the composite signal.

  corresponding information. According to the frequency sampling and rotation frequency, the window length and window shift are set to N w = 128 samples and R = 16, which gives a frequency resolution of 375 Hz. Therefore, there are N = 48000-Nw R + 1 = 2993 windows 82/129 covering one second signal. The number of Fourier transform (NFFT) is set as 256, then N k = N F F T /2 = 128. The quantile in Eq.(4.3) can be calculated as χ 2 v,1-α = 3253.1, with risk α = 0.05 and v = N + N k = 3121 degrees of freedom. The bearing details and its characteristic frequencies are listed in Table 5.2.

Figure 5 . 1 .

 51 As explained in the previous chapter, the parameter λ 2 returns the average duration of the transient pulse (in window length) augmented by the coda of the impulse response of the structure. As explained in Chapter 4.5, the latter can be measured by the vertical-intercept of the fitted lines between the duration time ∆t and the rotation frequency f r . The results of coda ∆e are reported in Table

Figure 5 . 1 :Figure 5 . 2 .

 5152 Figure 5.1: Quantification of the fault size together with 90% confidence intervals based on the assumption of a uniform error with standard deviation equal to R/F s .

10 4Figure 5 . 2 :

 1052 Figure 5.2: Spectrogram of the record #175 and the corresponding covariance of two states, C 1 and C 2 , in the right subplot.

Figure 5 . 3 :

 53 Figure 5.3: Spectral analysis of record #175. (a) Posterior probability spectrum S(k); (b) Envelope spectrum of bandpass signal with maximum kurtosis in band [5; 6] kHz selected from kurtogram.

2 2 Figure 5 . 4 :

 2254 Figure 5.4: Zoom of #175 in interval [0.5, 0.6] s. (a) Raw vibration signal; (b) the corresponding spectrogram; (c) posterior probability γ n (2); (d) duration time sequence τ n (2) of the active state; (e) duration time sequence τ n (1) of the inactive state.

Figure 5 . 5 :

 55 Figure 5.5: Fault signal reconstruction. (a) Raw vibration signal #175; (b) Reconstructed fault signal.

16 Figure 5 . 6 :

 1656 Figure 5.6: Empirical posterior probability p 1 for different values of the hyperparameter N w .

Figure 5 . 7 :

 57 Figure 5.7: Quantification of the fault size together with 90% confidence intervals based on the assumption of a uniform error with standard deviation equal to R/F s .

Figure 5 . 9 .

 59 According to the duration time sequence τ n (1) and τ n (2) shown in Figure 5.9(c)-(d), the corresponding mean values λ 1 and λ 2 are 20.50 and 6.76 windows. The corresponding duration time are 16.40 ms and 5.40 ms, respectively. The corrected impulse duration time is calculated by subtracting the coda error ∆e from ∆t, equal to 3.79 ms. Following the lines of Chapter 4.5, the fault size l estimate returned by EDHMM is about 1.76 mm. According to the indication given 93

Figure 5 . 8 :

 58 Figure 5.8: Spectral analysis of dataset # 12. (a) Posterior probability spectrum S(k); (b) Envelope spectrum of bandpass signal with maximum kurtosis in band [8.8; 10] kHz selected from kurtogram.

Figure 5 .

 5 Figure 5.10(a) has a kurtosis level of 6.98, which indicates a bearing damage. The kurtosis of the recovered signal was 23.02, which is even a stronger indication of a damaged bearing. In the signal of subplot (a), it is not possible to identify all the spikes of a rolling element fault and some transient pulses with low energy are masked by background noise. After signal reconstruction, subplot (b) shows a typical sequence of spikes with a time interval of about 22.2 ms, very close to the characteristic defect frequency.

Figure 5 . 10 :

 510 Figure 5.10: Fault signal reconstruction. (a) Raw vibration signal #12; (b) Reconstructed fault signal.

Figure 5 . 129 RecordingFigure 5 . 11 :

 5129511 Figure 5.11: Signal acquisition settings.

Figure 5 . 12 :

 512 Figure 5.12: Full life vibration signal in the horizontal direction of the bearing 3_1.

4 SIFigure 5 . 13 :

 4513 Figure 5.13, it is seen that there is no clear indication of failure in the first 40 h, where the bearing operation is stable as seen in Figure5.12. Data after about 40 h shows increasing

Figure 5 .

 5 14 (b) shows the vibration signal from 39 h, which is the interval of interest where the fault appear. The GLR Λ is calculated and shown by the blue line in Figure 5.14 (c). Each data point in the blue line represents the Λ for 1.28 seconds duration of the signal. The 98/129 time data is taken at 1 minute interval. Thus for 3.3 hours a total of 198 Λ points are calculated. The red dashed line represents the quantile level χ 2v,1-α /2. The higher the value of Λ with respect to the threshold χ 2 v,1-α /2, the higher the probability of detection, and the lower the risk of false alarm. It is seen in Figure5.14 (c) that the two lines cross around 39.6 h, which reflects the starting point of the failure. Compared with the time domain indicators plotted in Figure5.13, GLR Λ has the ability to provide a threshold for comparison and deciding the health status.

Figure 5 .

 5 14 (d) displays the evolution of the Poisson parameter λ 2 . This parameter controls the distribution of the transient pulse duration time, which can be estimated as the average number of windows that cover the transient pulse. A proportional relationship between the transient pulse duration and fault size seems a reasonable assumption. So λ 2 can be used as a quantity to characterize the fault size. In Figure 5.14 (d), each asterisk represents λ 2 for 1.28 seconds duration of the signal every two minutes. Then there are 81 obtained values for tracking quantitatively, rather than qualitatively like the other time domain indicators, the development of the failure. From the figure, it is seen that the λ 2 sequence gradually increases with the development of the failure, from the first value at about 5.0 to the end value at about 7.9.

2 Figure 5 . 14 :

 2514 Figure 5.14: Illustration of fault detection and fault development tracking. (a) accelerated degradation vibration signal; (b) interval of fault occurrence; (c) generalized likelihood ratio Λ for detecting fault occurrence; (d) Poisson distribution parameter λ 2 for tracking the fault development.

Figure 5 .Figure 5 . 15 :Figure 5 . 16 :

 5515516 Figure 5.16 shows the histogram of the frequency data of the incipient fault signal, f = 1/t s , s ∈ {1, 2, ..., S}. It is seen that most of the impulsive cycles fall within the frequency bin [112 135] Hz. This bin contains the outer race characteristic defect frequency f BP F O , indicating that the transient pulses in the signal were generated to a certain extent by the outer race. The other posterior probabilities of different fault types are listed in Table 5.11. It is seen the probability of the outer race p(f BP F O |t 1:S ) is the highest, confirming the outer race fault. Such statistical analysis can be realized automatically without the visual examination of a spectrum.

Figure 5 . 17 :

 517 Figure 5.17: Spectral analysis of the signal at 39.6 h. (a) Posterior probability spectrum S(k); (b) Envelope spectrum of bandpass signal with maximum kurtosis in band [3.2; 4.3] kHz selected from kurtogram.

Figure 5 . 129 5Figure 5 . 18 :

 5129518 Figure 5.18: Fault signal reconstruction. (a) Raw vibration signal at 39.6 h; (b) Reconstructed fault signal.

  This chapter has illustrated the application of the proposed integrated auto-diagnostic technique for detection, identification, quantification of a fault and fault signal reconstruction. Two kinds of experimental data have been considered in this chapter. The results on 36 sets of bearing data associated to three different fault cases has shown the efficiency of this framework. Almost all faults have been detected through the stochastic indicator GLR. Although, the bearing characteristic frequency was modulated by the shaft frequency thus producing a series of sidebands and harmonics, which has a big influence on the identification, the proposed identification technique overcame this problem well. For the signal reconstruction technique, the transients interval became clearer after fault signal reconstruction. In addition, from Figure5.1 and 5.7, the fault size has been assessed close to the real size, and within the 90% confidence intervals, thus validating the technique.

  methods for different diagnostic tasks in one integrated framework. The remarkable advantages are that it performs without the historic data, without few hyperparameters, and without any manual intervention in the process. The literature review shows that the most existing techniques are focused on one specific diagnostic task. But these techniques are almost independent, which requires users to have a high level of experience and knowledge when dealing with different diagnostic problems. Some methods attempt to fulfill such aspiration, yet not reaching the automated and integrated level as this thesis does.The thesis introduced an automated diagnosis framework rooted on a stochastic model of the vibration signals produced by rolling element bearing faults. The proposed method can achieve the fault detection, fault type identification, fault signal reconstruction and fault size characterization at once. The distinguishing feature of this method is that it is performed automatically, which means that no human intervention or manual tuning of the algorithms is required for solving these subtasks. The whole process depends mainly on only one key hyperparameter, i.e. the window length N w in the STFT. Firstly, a technique based on hypothesis test is introduced to detect the healthy status of the vibration signal. The test statistics is the GLR constructed through the EDHMM parameter, i.e. the covariance matrix C in different states, and it asymptotically follows a chi-squared distribution according to Wilk's theorem. Unlike the supervised approach, this technique does not require prior knowledge or training data. Secondly, the fault signal is extracted from the noisy vibration signal based on the 105/129 Bayesian theorem and EDHMM parameters. After modeling through EDHMM, the STFT coefficients related to the failure can be singled out and estimated as the active state, and then the inverse STFT is used to reconstruct the fault signal. Compared with the Wiener filter, this technique can be seen as a time-varying filter, which takes account of the posterior probability γ n . Subsequently, two different techniques are proposed based on the EDHMM parameters to identify the fault type. The PPS is the spectrum of the posterior probability, a robust alternative to the standard envelope spectrum, free of any pre-processing. Unlike the signal spectrum, the PPS is not sensitive to the amplitude of the transient pulses, so it can avoid the influence of modulation, thereby reducing the sidebands in the spectrum. However, in order to completely escape visual inspection, a simple statistical method is introduced for identifying automatically the fault type. It is based on the duration time of each transient period, to calculate the probabilities of different fault types. Then the fault identification is achieved by a probability comparison. At last, the fault size quantification is addressed. This technique is based on the hypothesis that a proportional relationship exists between the transient pulse duration and fault size. So we can get the linear relationship between the duration time and fault size based on the geometry structure of the rolling element bearing. Therefore, after some mathematical manipulation, the estimated fault size l and the coda error ∆e can be estimated as the slope and the intercept, respectively, of the fitted lines with respect to machine speed. The Poisson parameter λ 2 , i.e. the expected duration time of transient pulse, can also track quantitatively the fault size evolution. The above techniques for different diagnostic tasks are all based on the EDHMM parameters. A nonstationary bearing signal cannot be entirely modeled by a deterministic model, while EDHMM, as one stochastic model, has gifted flexibility and adaptability to cope with the complexity of non-stationary phenomena. Finally, it should be mentioned that the proposed integrated auto-diagnostic framework strongly depends on this stochastic model, which poses a risk to the diagnosis. This will be considered in future research. In the last chapter, the capability of the proposed automated diagnosis framework has been validated through two different experimental scenarios. The results of the experimental validation not only meets the expectation of diagnosis, but also reflect its unique advantages, i.e. the automation and integration. 106/129

•

  EDHMM in future work to achieve the integrated diagnosis for variable rotating speed signals. • Diversity of applicable object This research attempted to establish an integrated diagnosis framework for rolling element bearings. So another future work is the application of integrated diagnosis in condition monitoring applications for other rotational machine components such as gears, pumps, blades etc. Comprehensiveness of the framework One last recommendation for future scientific work is related to the embedding of the RUL forecast. The reliable RUL estimation is a challenging task and also very useful in real CBM. In the current framework, the Poisson parameter λ 2 (i.e. the expected duration time of transient pulse) can quantitatively track the evolution of the fault size. Since this consideration could lead to increased prognostic performance, further 107/129 research effort in this field is suggested.108/129L(θ) in each iteration. As maximizing L(θ) with respect to the parameters θ, then we remove the items that are not related to the parameters θ. And using the Q(θ, θ old ) as the new objective function, then it can be simplified asQ(θ, θ old ) = Z p(Z|Y , θ old ) ln p(Y , Z|θ). (A.2)We substitute the joint probability ln p(X, Z|θ) into Q(θ, θ old ) and make use the definitions of γ and ξ, Q(θ, θ old ) can be written as,Q(θ, θ old ) = i γ 1 (i) ln π i + N n=2 i j ξ n (i,j) ln a ij + N n=1 i γ n (i) ln p(Y n |θ k ). (A.3)

  'UNIVERSITE DE LYON OPEREE AU SEIN DE L'INSA LYON NOM : JIN DATE de SOUTENANCE : 25/05/2022 Prénoms : Yaqiang TITRE : Integrated auto-diagnosis based on stochastic model for rolling element bearings NATURE : Doctorat Numéro d'ordre : 2022LYSEI045 Ecole doctorale : Mécanique, énergétique, Génie Civil, Acoustique (MEGA)

Chapter 2 State of the art Contents 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Rolling element bearing signal characteristics . . . . . . . . .

	2.2.1 Basic structure of rolling element bearing . . . . . . . . . . . .
	2.2.2 Vibration characteristics of rolling element bearings . . . . . .
	2.2.3 Incipient fault characteristics of rolling element bearing . . . .

2.3 General review of diagnosis of rolling element bearings . . .

	2.3.1 Bearing fault detection . . . . . . . . . . . . . . . . . . . . . . .
	2.3.2 Bearing fault signal enhancement . . . . . . . . . . . . . . . . .
	2.3.3 Bearing fault identification . . . . . . . . . . . . . . . . . . . .
	2.3.4 Bearing fault size characterization . . . . . . . . . . . . . . . .

2.4 Hidden Markov model and its variants . . . . . . . . . . . . . 2

  .4.1 Hidden Markov model . . . . . . . . . . . . . . . . . . . . . . . 2.4.2 Explicit duration hidden Markov model . . . . . . . . . . . . .

2.5 Integrated automatic diagnosis . . . . . . . . . . . . . . . . . . 2.6 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	9/129

Modulation characteristic

	Rolling bearing vibration signals have typical multi-modulation characteristics. Mc-
	fadden, et al. [106] proposed firstly the vibration model of a single point defect in
	a rolling element bearing. It can be seen from the model that the vibration signal
	caused by the localised fault can be regarded as a combination of a low-frequency
	part (faulty impact) and a high-frequency part (natural resonance frequency). In
	diagnosis, the actual resonance and other modulation are not of interest, but only
	the fault impulse demodulated signature.

2.3.1.1 Time domain indicators

	usually used to monitor the health condition of bearings. The time-domain dimensionless
	indicators are not only independent on loads and speeds of rotating machinery, but they
	are also able to effectively indicate early faults occurring in the rotating machinery [84].
	There are six indicators usually used for the fault diagnosis of rolling element bearing,
	shown as follows:
	Skewness (SK) [114]
	subsection, two different kinds of detection techniques are reviewed, i.e. the time-domain
	indicators and the stochastic indicators.
	Because the time domain indicators have the advantage of intuition, are easy to
	understand, convenient to calculate, and high efficiency, they have been widely used in the
	field of bearing detection. The time domain indicators can be divided into two categories,
	i.e. dimensional indicators and dimensionless indicators. The dimensional indicators mainly
	include Root Mean Square (RMS), Standard Deviation, Effective value etc. However,
	these indicators are usually closely related to the load and speed of rotating machines, and
	sensitive to the operating conditions. In practical applications, dimensionless indicators are
	14/129

  MED), as a linear analysis method, was first proposed by Wiggins and applied to the fault diagnosis of bearings by Sawalhi[START_REF] Sawalhi | The enhancement of fault detection and diagnosis in rolling element bearings using minimum entropy deconvolution combined with spectral kurtosis[END_REF]. Cheng et al.[START_REF] Cheng | Application of an improved minimum entropy deconvolution method for railway rolling element bearing fault diagnosis[END_REF] proposed an improved MED method, which sloves the filter coefficients by the standard particle swarm optimization algorithm, assisted by a generalized spherical coordinate transformation. The proposed method delivered better deconvolution performance than the classical MED method. McDonald et al.[START_REF] Mcdonald | Maximum correlated Kurtosis deconvolution and application on gear tooth chip fault detection[END_REF] proposed the Maximum Correlation Kurtosis deconvolution.Based on correlation kurtosis as an evaluation index, the deconvolution was realized through an iterative process. In this way, continuous impulse sequences that are submerged by noise in the signal were highlighted. Recently, a novel blind deconvolution method rooted

Thus, blind deconvolution is very suitable for the monitoring and diagnosis of rolling bearing impact fault. There are different 26/129 blind deconvolution methods based on different optimization criteria. Minimum entropy deconvolution (on the maximization of the cyclostationarity of the excitation -as typically encountered with machine faults -was proposed [34]. This method is based on the generalized Rayleigh quotient and solved by means of an iterative eigenvalue decomposition algorithm. The comparison results reveal superior other blind deconvolution methods existing in the literature.

4 Hidden Markov model and its variants 2.4.1 Hidden Markov model

  90• rectangular edges, which is inconsistent with the facts. Larizza et al. developed a numerical model for a rolling element bearing, and a new defect size estimation method was developed to improve the accuracy in estimating the size of a defect that has sloping leading and trailing edges.In this section, one type of stochastic signal model will be focused, namely the hidden

2.

Markov model (HMM). It is a good tool to model the transient pattern in nonstationary signals. The basic mathematical theory was published in a series of classic papers by Baum and his colleagues

[START_REF] Leonard | Statistical inference for probabilistic functions of finite state Markov chains[END_REF][START_REF] Leonard | An inequality with applications to statistical estimation for probabilistic functions of Markov processes and to a model for ecology[END_REF]

in the late 1960s and early 1970s. Later, it became an important research direction in signal processing; especially in the field of speech recognition. After 2000s, HMMs have become more and more popular in various fields, such as, medical signal modeling, facial expression recognition, gene prediction, gesture recognition, musical composition and vibration signal analysis. There are two strong

Table 2 .

 2 is the probability of taking the self-transition at state i for d times. Because of this property, the Markov chain has no idea of how long it has stayed in a state and when to leave the state. In many cases this property seems to be a strong limitation that results in inaccurate signal modeling[START_REF] Goutsias | Optimal simultaneous detection and estimation of filtered discrete semi-Markov chains[END_REF]. Researchers have proposed a number of techniques to address these limitations. In 1980, Ferguson[START_REF] Freguson | Variable duration models for speech[END_REF] modeled the self-transition process with an explicit duration probability distribution for addressing this problem. 1, following the diagnosis objects, the objective, the type of data and methodology.The table is supposed to provide a more direct view of reported studies to readers who are concerned with the use of EDHMM in prognostics and health management (PHM).

	Dong et al. [48] were the first researchers to apply EDHMM to PHM in 2006. They
	validated the effectiveness of the EDHMM for the fault classification of UH-60A Blackhawk

Suppose the self-transition (from state i to state i) probability is denoted as a ii . The probability of state occupancy decreases exponentially with time, and the probability of d consecutive 31/129 observations in state i can be written as p(d|i) = a d ii (1 -a ii ). (2.16) p(d|i) The explicit duration hidden Markov model can be viewed as a special instance of the standard HMM. What makes waiting time hidden Markov model different from the standard HMM is that there exists a duration time d in each state. The duration distribution p(d|j) can be either a discrete distribution or a continuous density, a nonparametric or a parametric density. When the underlying process enters state j for duration d with the probability p(d|j), there are d observations produced. Therefore, the self-transitions are prohibited in the explicit duration hidden Markov model. This makes it suitable for use in a wide range of applications. The capacity and complexity of EDHMM are analyzed in [72], which points that standard hidden Markov model coupled with a moderate increase in overall topological complexity and state distribution parameter tying, are already well suited to handling nonexponential duration distributions. Since it was proposed, EDHMM has attracted the attention of scholars in various fields. In this section, only the application in fault diagnosis and prognosis is given. The application to fault diagnosis and prognosis appeared gradually after 2000, in particular in the recent ten years. The related literature is synthesized in one table, seen in main transmission planetary carriers and prognosis of a hydraulic pump health monitoring application. He et al. pointed out the drawback of EDHMMs in [66]: the computational 32/129

Table 2 . 1 :

 21 Summary of the use of EDHMM for fault diagnosis for inference and parameter estimation. But duration modeling by parametric probability distributions has been used to alleviate the computational burden.

	Object	Reference		Objective	Type of observation	Methodology
	planetary box	Dong et al. [48]	fault classification	time domain vibration sig-	original EDHMM
						nal
		Fan et al. [57]		diagnosis and prog-	time domain vibration sig-	original EDHMM
					nosis	nal
	hydraulic pump	Dong et al.[48]	prognosis	time domain vibration sig-	original EDHMM
						nal
		Dong et al. [46]	fault classification	time domain vibration sig-	segmental EDHMM
					and prognosis	nal
		Dong et al. [47]	fault classification	wavelet decomposition co-	segmental EDHMM
					and prognosis	efficients
		Dong et al. [45]	fault classification	wavelet decomposition co-	auto-regressive EDHMM
					and prognosis	efficients
		Dong et al. [49]	RUL prediction	wavelet decomposition co-	non-stationary segmental
						efficients	EDHMM
		Dong et al. [97,	RUL prediction	wavelet decomposition co-	joint EDHMM with sequen-
		95]				efficients	tial Monte Carlo method
		Dong et al. [96]	RUL prediction	wavelet decomposition co-	adaptive EDHMM
						efficients
	rotor	He et al. [66]		prognosis	frequency domain vibra-	original EDHMM
						tion signal
	shaft	Bechhoefer et al.	prognosis	frequency domain vibra-	original EDHMM
		[22]				tion signal
	gearbox	Teng et al. [144] RUL prediction	time-frequency domain of	original EDHMM
						vibration signal
		Li et al. [90]		fault detection and	time domain vibration sig-	multivariate	Bayesian
					RUL prediction	nal	control scheme based on
							EDHMM
		Li et al. [91]		fault detection and	time domain vibration sig-	Optimal Bayesian control
					RUL prediction	nal	scheme based on EDHMM
	hydraulic cylin-	Su et al.[142]		RUL prediction	cylinder wear data	original EDHMM
	der
		Huang et al.	fault diagnosis	wavelet decomposition co-	original EDHMM
		[70]				efficients
		Xiao et al. [159]	RUL prediction	dynamic pressure signal	EDHMM combined with
							high-order particle filter
	bearing	Wu et al. [158]	wear prediction	certain data(unspecified in	modified EDHMM
						the paper)
		Chen et al. [41]	prognosis	time domain vibration sig-	mixtured	non-gaussian
						nal	based EDHMM
		Wang	et	al.	prognosis	time domain vibration sig-	Duration-dependent
		[153]				nal	EDHMM
		Cartella et al.	RUL prediction	time domain vibration sig-	EDHMM combined with
		[36]				nal	AIC
		Le et al. [82, 83] RUL prediction	simulation vibration signal	Multi-Branch EDHMM
	draught fan	Lin et al. [94]		prognosis	time domain vibration sig-	original EDHMM
						nal
	milling tool	Liu et al. [98]		RUL prediction	cutting force signal	modified EDHMM
		Duan et al. [52]	RUL prediction	cutting force signal	EDHMM combined with
							vector	autoregressive
							degradation modeling
	transformer	Hao et al. [65]	RUL prediction	vibration signal	time-varying EDHMM
	fused deposition	Wu et al. [156]	Condition monitor-	AE signal	original EDHMM
	modeling				ing
	unknown object	Khaleghei et al.	RUL prediction	simulation data	EDHMM
		[76]
					33/129

  order to motivate the EDHMM, let first consider the signal model. Let y(t) denotes the measured vibration signal, x(t) the series of impulses due to impacts on a

fault and n(t) the background noise, which, by definition, embodies all other sources of vibration that are not of interest for diagnosis. Clearly, in the healthy state,

  Therefore, the hidden states z 1:N take two values, z n ∈ {1, 2}. In each state i, the whole duration is denoted as d i , which means there are d i successive observations of the signal belonging to the state i. Here, the Poisson distribution with parameter λ i ,d i ∼ P (λ i ), is used to model the discrete variable d i , although other distribution choices are equally possible,

	1. "inactive" state, hereafter labeled #1, where noise only is present,
	2. "active" state, labeled #2, where an impulse occurs.
	47/129

Each observation Y n belongs to one unknown state, denoted as z n . Then the corresponding state sequence of the observations Y 1:N is denoted by z 1:N {z 1 , z 2 , ..., z N }. In the present application, the faulty bearing vibration signal comprises two states:

  λ λ i is the parameter of a Poisson distribution that models the number of STFT windows covering the state i, d i . It can be estimated as the expected value of the duration time in different states. The parameter λ 2 represents the average number of windows that The fault size quantification based on λ 2 will be addressed in Chapter 4.5.In order to illustrate the discrepancy of this parameter, 16 different sets of inner race fault data from the CWRU bearing data center are used here, including 4 different rotating speeds and 2 fault sizes. For statistical purposes, it is required to divide each dataset into several short signals through sliding segmentation with a overlap, to obtain more dataset. Then we can obtain a point (λ 1 , λ 2 ) for each segment of the signal, as shown in Figure3.10.

											BPFI; 1796rpm; 0.18mm
											BPFI; 1772rpm; 0.18mm
	2 (duration time of transient pulse)	4 6 8 10 12 2									BPFI; 1748rpm; 0.18mm BPFI; 1721rpm; 0.18mm BPFI; 1796rpm; 0.54mm BPFI; 1772rpm; 0.54mm BPFI; 1748rpm; 0.54mm BPFI; 1721rpm; 0.54mm BPFO; 1796rpm; 0.18mm BPFO; 1772rpm; 0.18mm BPFO; 1748rpm; 0.18mm BPFO; 1721rpm; 0.18mm BPFO; 1796rpm; 0.54mm BPFO; 1772rpm; 0.54mm BPFO; 1748rpm; 0.54mm
											BPFO; 1721rpm; 0.54mm
		0
		0	5	10	15	20	25	30	35	40	45
						1

cover a transient pulse. It will increase as the fault size increases, so it will be a good indicator to assess the fault size. This parameter is proposed first time to used in fault 59/129 size estimation.

Table 5 . 1 :

 51 List of the inner race fault datasets and their corresponding basic information /2, whereas the test statistic Λ of the faulty data are greater than this threshold except for #174. The "Identification" column lists the posterior probabilities of four fault types given the impulse cycle t 1:S , where,

	Dataset Accelerometer location	Fault Size Motor load Motor Speed [mm] [hp] [rpm]
	#097				0	1796
	#098 #099		0		1 2	1772 1750
	#100				3	1730
	#109				0	1796
	#110 #111		0.18	1 2	1772 1750
	#112 #174	Drive end			3 0	1730 1796
	#175 #176		0.36	1 2	1772 1750
	#177				3	1730
	#213				0	1796
	#214 #215		0.54	1 2	1772 1750
	#217				3	1730
	Table 5.2: The basic information bearing details and fault frequencies
		Bearing type	SKF 6205-2RS JEM
		Position		Drive end
			BPFO	3.585
		Characteristic frequencies BPFI	5.415
		(×rotation frequency f r )	BSF	2.357
				FTF	0.3983
	Table 5.3 summarize the inner race fault diagnosis results, including fault detection,
	fault identification, and fault size quantification. In order to demonstrate the effectiveness
	of the detection algorithm, four healthy data records (#097 -#100) are used as a basis
	of comparison. In the "Detection" column, it is found that the test statistic Λ of the
	healthy data are all less than the quantile χ 2 v,1-α 83/129

Table 5 . 3 :

 53 Inner race fault diagnosis results

	Dataset	Λ	Detection Result	p 1	Identification p 2 p 3	p 4	Characterization λ 2 l [mm]
	#097	503.4 Accept H 0		--	--	--
	#098	491.7 Accept H 0		--	--	--
	#099	727.1 Accept H 0		--	--	--
	#100	529.6 Accept H 0		--	--	--
	#109	83606 Reject H 0	1	0	0	0 5.23	0.14
	#110	107329 Reject H 0	1	0	0	0 6.11	0.19
	#111	85730 Reject H 0 0.98 0.02	0	0 6.18	0.20
	#112	98226 Reject H 0 0.96 0.04	0	0 5.99	0.15
	#174	NAN	NAN		--	--	--
	#175	84173 Reject H 0 0.70 0.30	0	0 6.26	0.27
	#176	85226 Reject H 0 0.57 0.09 0.34 0 6.55	0.37
	#177	73326 Reject H 0 0.55 0.06 0.39 0 6.59	0.39
	#213	71309 Reject H 0	1	0	0	0 7.29	0.68
	#214	49428 Reject H 0 0.97 0.03	0	0 7.10	0.57
	#215	51570 Reject H 0	1	0	0	0 7.28	0.66
	#217	26051 Reject H 0 0.93	0	0.07 0 7.11	0.56

Table 5 . 4 :

 54 The coda error ∆e[ms]

	Fault size (mm) 0.18 0.36 0.54
	∆e (ms)	1.92 1.96 2.03

Table 5 .

 5 6 displays the diagnosis results in the case of an outer race fault. It shows that faults of sizes 0.18 mm (#135 -#138) and 0.54 mm (#238 -#241) are detected correctly. They are also correctly identified as outer-race faults with the maximum probability p(f BP F O |t 1:S ). In addition, the last column of the table shows that the estimated fault sizes are reasonable. However, the

	88/129

Table 5 . 5 :

 55 List of the outer race fault data datasets and their corresponding basic information

	Dataset Accelerometer location	Fault Size Motor load Motor Speed [mm] [hp] [rpm]
	#135			0
	#136 #137		0.18	1 2
	#138			3
	#201			0
	#202 #203	Drive end	0.36	1 2
	#204			3
	#238			0
	#239 #240		0.54	1 2
	#241			3

Table 5 . 6 :

 56 Outer racer fault diagnosis results

	Dataset	Λ	Detection Result	p 1	Identification p 2 p 3	p 4	Characterization λ 2 l [mm]
	#135	321543 Reject H 0	0	1	0	0 8.65	0.15
	#136	284131 Reject H 0	0	1	0	0 9.24	0.21
	#137	268814 Reject H 0	0	1	0	0 8.99	0.20
	#138	265902 Reject H 0	0	1	0	0 9.40	0.23
	#201	2797.8 Reject H 0 0.76 0.24	0	0 --	--
	#202	385.2 Accept H 0		--			--	--
	#203	275.2 Accept H 0		--			--	--
	#204	153.8 Accept H 0		--			--	--
	#238	210946 Reject H 0 0.03 0.97	0	0 9.26	0.45
	#239	221400 Reject H 0	0	1	0	0 10.10	0.62
	#240	200886 Reject H 0	0	1	0	0 9.85	0.61
	#241	206363 Reject H 0	0	0.82 0.18 0 10.37	0.68
	5.2.

3 Case 3: Rolling element fault data

  In this subsection, the proposed automated diagnosis framework is illustrated on the dataset from Technology Research Group of the Polytechnic University of Madrid[START_REF] César | Bearing Database[END_REF].

	22205E1KC3), with different fault depths, recorded by accelerometers in three directions.
	The dataset comprises vibration signals from a double-row spherical roller bearing (FAG
	90/129

Table 5 . 7 :

 57 List of the rolling element fault datasets and their corresponding basic information

	Dataset	Accelerometer location	Fault Size Motor Speed [mm] [rpm]
	#10			200
	#11		F1(0.006)	350
	#12			500
	#19			200
	#20		F2(0.014)	350
	#21 #31	6 o'clock of the bearing caeing		500 200
	#32		F3(0.019)	350
	#33			500
	#40			200
	#41		F4(0.027)	350
	#42			500

Table 5 . 8 :

 58 Bearing details and fault frequencies N k = 1130. Similarly, the degrees of freedom for 500 rpm can also be obtained as p 2 = 1747. Accordingly, the statistical threshold for fault detection can be calculated as χ 2 1130,0.95 /2 = 605.2 and χ 2 1747,1-α /2 = 923.2. The results for fault detection, fault identification and fault size quantification are reported in Table5.9. These faults are correctly detected and identified in all cases.Following similar lines as in the first case, the coda ∆e is assessed from the relationship between parameter λ 2 and rotation frequency f r . The estimated fault size based on the corrected duration time is displayed in the "Characterization" column. Figure5.7 illustrates the fault size quantification through the model parameter λ 2 . From this figure, it is found that the estimated value is slightly influenced by the motor rotation speed, thus illustrating the overall uncertainty of the method. The errors caused by fluctuations are located in the 90% confidence intervals based on the assumption of a uniform error with standard deviation equal to R/F s . From the figure, it is found the error of the fault size quantification

	Bearing type	FAG 22205E1KC3
		BPFO	6.1852
	Characteristic frequencies BPFI	8.8148
	(× rotation frequency f r )	BSF	5.4030
		FTF	0.4123
	The duration of the records is one second (40000 samples). According to the sampling
	frequency and rotating frequency, the window length N w and NFFT are set to 512 and
	1024 respectively. Due to the large difference between the rotation speeds in this case,
	different window shifts R are used, 64 for 200 rpm and 350 rpm and 32 for 500 rpm.

Therefore, we can obtain the number of windows in STFT N = 40000-Nw R + 1 = 618 and the N k = N F F T /2 = 512. So the degrees of freedom for 200 rpm and 350 rpm data is 91/129 p 1 = N + probability spectrum S(k) of the posterior probability sequence is displayed in Figure 5.8(a). It is seen that the fault frequency f BSF and its harmonics dominate the spectrum,

  This datasets contain complete run-to-failure data of 15 rolling element bearings, and one dataset named Bearing 3_1 with the end of outer race fault is selected to demonstrate the effectiveness of the proposed integrated diagnosis framework. The type of tested bearings is LDK UER204, and the detailed parameters are given in Table5.10. The sampling frequency is 25.6 kHz. As

	5.3. Accelerated degradation data
		1.95		2	1	2
	(a)	1.9
		1.85
			0.5	0.52	0.54	0.56	0.58	0.6
			2
	(b) 0.8	1		0.9		1
			0		Time [s]
			0.5	0.52	0.54	0.56	0.58	0.6
		15
	(c)	5 10
			0
			0.5	0.52	0.54	0.56	0.58	0.6
		30
	(d)	10 20
			0
			0.5	0.52	0.54	0.56	0.58	0.6
	shown in Figure 5.11, a total of 32768 sampling points (i.e. 1.28 s) are recorded for each
	sampling, and the sampling period is equal to 1 min.
	The total operating life is 42 h 18 min. The rotating speed is constant at 2400 rpm,

Figure 5.9: Zoom of #12 in interval [0.5 0.6] s. (a) Vibration signal; (b) posterior probability γ n (2); (c) duration time sequence τ n (2) of active state; (d) duration time sequence τ n (1) of inactive state. bearing, which is eventually run to failure, and also in tracking fault size evolution. Results from naturally occurring fault of new bearings are more difficult to analysis than artificial defects.

The used datasets are provided by the Institute of Design Science and Basic Component at Xi'an Jiaotong University (XJTU)

[START_REF] Wang | A hybrid prognostics approach for estimating remaining useful life of rolling element bearings[END_REF]

.

and the radial load is 10 kN. Because the load is applied in the horizontal direction, the accelerometer placed in this direction is able to capture more degradation information of 95/129

Table 5 .

 5 10: Bearing details and fault frequencies

	Bearing type		LDK UER204
	Outside diameter	39.80mm Inside diameter 29.30 mm
	Pitch diameter	34.55mm	Ball diameter	7.92 mm
	Number of balls	8	Contact angle	0 •
	Dynamic load rating 12.82 kN Static load rating 6.65 kN
			BPFO	3.0831
	Characteristic frequencies	BPFI	4.9169
	(× rotation frequency f r )	BSF	2.0666
			FTF	0.3854
	the tested bearings. Therefore, the horizontal vibration signal is selected to demonstrate
	the effectiveness of the proposed integrated diagnosis framework.

Table 5 .

 5 11: The posterior probabilities of fault type given t 1:S BP F I |t 1:S ) p(f BP F O |t 1:S ) p(f BSF |t 1:S ) p(f F T F |t 1:S )

	Fault type	BPFI	BPFO	BSF	FTF
	Posterior probability p(f 0.18	0.66	0.14	0.02

  En présence de défaut naissant, les vibrations des roulements présentent des signatures symptomatiques sous forme d'impulsions répétitives. Cela peut être vu comme un signal non stationnaire dont les propriétés statistiques basculent entre deux états. La stratégie de maintenance proposée modélise ces caractéristiques avec un modèle de Markov caché à durée explicite (EDHMM) et utilise les paramètres estimés du modèle pour effectuer un diagnostic intégré sans nécessiter l'expertise de l'utilisateur. La détection d'un défaut est d'abord réalisée au moyen d'un test de rapport de vraisemblance construit sur les paramètres de l'EDHMM. Une approche de comptage statistique et une probabilité sont ensuite utilisées pour identifier automatiquement le type de défaut. Afin d'obtenir le signal de défaut dans certains cas, un filtre bayésien basé sur les paramètres EDHMM est construit. Enfin, la taille du défaut est estimée à partir des durées renvoyés par EDHMM.Par la suite, la capacité du cadre d'autodiagnostic intégré est illustrée sur différents ensembles de données expérimentales. La première validation est réalisée sur les données de vibration pour des conditions spécifiques. Les résultats indiquent une diagnostic robuste et précise du roulement à éléments roulants. De plus, le résultat sur des données de dégradation accélérée montre également l'efficacité de la méthode, en particulier la capacité à détecter l'occurrence d'une défaillance et de suivre quantitativement son développement. Cette technique a un potentiel d'utilisation en CBM. MOTS-CLÉS : maintenance du renseignement; Roulement à élément roulant; Auto-framework intégré; Détection de fautes; Identification du type de défaut; Reconstruction du signal de défaut; Caractérisation de la taille des défauts. Laboratoire (s) de recherche : Laboratoire de vibration et acoustique (LVA)

Directeur de thèse: Jérôme Antoni Président de jury : Composition du jury : GRYLLIAS, Konstantinos (rapporteur) WYLOMANSKA, Agnieszka (rapporteur) EL BADAOUI, Mohamed (examinateur) ANTONI, Jérôme (directeur de thèse)

Cette thèse est accessible à l'adresse : https://theses.insa-lyon.fr/publication/2022LYSEI045/these.pdf © [Y. Jin], [2022], INSA Lyon, tous droits réservés

Acknowledgments

Appendix A

Parameter estimation

In this section, the detail mathematical process about the Eq.(3.23) and (3.25) is demonstrated. The logarithmic likelihood function of the observation is obtained by marginalizing over the hidden variables as,

where, Z = z 1:N . The goal is to find the optimal parameter θ that maximizes the likelihood function L(θ). But in this likelihood function, it is found that the summation over the hidden states appears inside the logarithm. As the hidden variables Z are also unknown, which brings difficulties to direct maximization. We therefore turn to the iterative algorithm to find an efficient framework for maximizing indirectly the likelihood function. Some manipulation on the logarithmic likelihood function L(θ) are performed as follow,

The updated θ new can be found by maximizing the logarithmic likelihood function

109/129

Appendix B

Proof of the convergence

In this section, the detailed proof about the convergence of the iterative algorithm in Chapter 3.4 is given. Intuitively, as long as it can be proved that the new estimated parameter θ new makes the logarithmic likelihood function L(θ new ) greater than the previous iteration, then the iterative algorithm converges. The proof of the convergence is equivalent to prove L(θ new ) > L(θ old ). If proven, the L(θ) will be maximized after a sufficient number of iterations. Some manipulations about the logarithmic likelihood function L(θ) are made as follow,

where,