D Abi-Abdallah
email: dima.abiabdallah@utc.fr

A Drochon

V Robin

P Poulet

O Fokapu

REMOVING THE MHD ARTIFACTS FROM THE ECG SIGNAL FOR CARDIAC MRI SYNCHRONIZATION

Blood flow in high static magnetic fields induces elevated voltages that disrupt the ECG signal recorded simultaneously during MRI scans for synchronization purposes. This is known as the magnetohydrodynamic (MHD) effect, it increases the amplitude of the T wave, thus hindering correct R peak detection. In this paper, we present an algorithm for extracting an efficient reference signal from an ECG contaminated by the Nuclear Magnetic Resonance (NMR) environment, that performs a good separation of the R-wave and the MHD artifacts. The proposed signal processing method is based on subband decomposition using the wavelet transform, and has been tested on human and small rodents ECG signals acquired during MRI scans in various magnetic field intensities. The results showed an almost flawless trigger generation in fields up to 4.7 Tesla during the three tested imaging sequences (GE, FSE and IRSE)

Introduction

Magnetic Resonance Imaging (MRI) of a moving organ, such as a beating heart, requires facilities to trigger successive scans and perform synchronization duties. Since an image cannot be wholly obtained during one heart cycle, consecutive acquisitions have to be accurately combined with the heart's phase motion. The easiest way to achieve this purpose consists in simultaneously recording the electrocardiogram (ECG) and detecting its R-waves to trigger the acquisitions. Thus, achieving correct synchronization depends greatly on accurate peak detection in order to guarantee that accurate data collection always starts at the same point of the cardiac cycle.

However the actual technical requirement needed to obtain good image quality (high static magnetic field, intense gradients, short rising time…) constrain the synchronization process, by generating high levels of artifacts that get superimposed on the ECG signal making the R peak detection very difficult. Most of the current synchronization methods fail in a significant number of patients, resulting in poor image quality and low scan efficiency. Artifacts generated by gradient switching, radiofrequency pulses and patient movements can be considerably reduced by the use appropriate acquisition and filtering devices [START_REF] Laudon | Minimizing Interference from Magnetic Resonance Imagers During Electrocardiography[END_REF], [START_REF] Felblinger | Restoration of electrophysiological signals distorted by inductive effects of magnetic field gradients during MR sequences[END_REF]. Nevertheless, previously proposed solutions for NMR artifacts suppression seem to be inefficient for the magnetohydrodynamic (MHD) artifact removal.

The MHD effect is of physiological origins and is due to the blood motion in the magnetic environment. In fact, the blood's charged particles flowing in the static magnetic field get deflected by the Lorentz force and as a result a Hall potential is generated across the vessel walls. The magnitude of these voltages is a function of the blood flow rate and vessel diameter, their amplitude is maximal when the blood flow is perpendicular to the magnetic field. These superimposed voltages, induced by the electrodynamic interactions of the static field with arterial blood flow, especially in the ascending aorta, alter the waveform of the recorded ECG. The largest magnetically induced voltages occur during the blood ejection phase in the aorta which coincides with the ventricle repolarization (T wave). Hence the major MHD-caused change observed on the simultaneously recorded ECG is an increase of the T wave amplitude during in vivo MRI observations. Magnetic field interactions with blow flow have been demonstrated by multiple authors throughout in vitro animal models experiments [START_REF] Keltner | Magnetohydrodynamics of blood flow[END_REF] as well as in vivo studies [START_REF] Gaffey | Alterations in the rat electrocardiogram induced by stationary magnetic fields[END_REF], [START_REF] Tenforde | Cardiovascular Alterations in Macaca Monkeys Exposed to Stationary Magnetic Fields: Experimental Observations and theoretical Analysis[END_REF] where ECG alterations have been observed. For small animals, such as rats, the T wave elevation could be quantitatively measured at field levels above 0.3 T [3], however for larger species such as monkeys, who have larger aortas and greater flow rates, the MHD effect could be noted for much lower field intensities starting 0.1 T [START_REF] Tenforde | Cardiovascular Alterations in Macaca Monkeys Exposed to Stationary Magnetic Fields: Experimental Observations and theoretical Analysis[END_REF]. At strong intensity fields the T wave might reach amplitudes as high as the QRS complex and may be mistaken for an R peak, hence leading to incorrect acquisition triggering.

In this paper, we present an algorithm, based on wavelet sub-band decomposition that allows correct trigger extraction form an ECG contaminated by the NMR environment.

Materials and Methods

Experimental Setup:

The signals were recorded using a MR compatible ECG sensor which technical details are provided in [START_REF] Felblinger | Electrocardiogram acquisition during MR examinations for patient monitoring and sequence triggering[END_REF] and a digital-signal-processing system placed outside the Faraday screen. Signals were collected during MRI scans in three different MRI platforms: 2T and 4.7T magnets for in vivo small animal exploration, and a 1.5T Signa GE MRI scanner.

The Algorithm:

The frequency distribution of the ECG wave forms is exploited here to extract the necessary information for synchronization. The MHD artifact is superimposed on the T wave and is essentially located in the low frequencies. The main idea is to isolate the frequency band that comprises mostly of QRS energy. This is done by decomposing the contaminated ECG into frequency sub-bands using wavelets than constructing a reference signal by keeping only the details containing the maximum QRS energy. The algorithm, illustrated in figure 1, includes three main steps: (i) Frequency sub-band decomposition: The corrupted ECG signal is decomposed onto N levels using filter banks, based on the classical orthogonal discrete wavelet transform (DWT). The number of scales was fixed according to ECG spectral taking into account the frequency components of the contaminating artifacts. For an ECG signal, sampled at 1 kHz, the decomposition would consist of 8 frequency sub-bands resulting in the following filter banks bandwidth: 1.96, 3.91, 7.81, 15.63, 31.25, 62.5, 125, and 250 Hz. For each level, the detail signals are afterwards reconstructed using the inverse discrete wavelet transform (IDWT). As for the choice of the analyzing wavelet, a bibliographical study associated with an analysis of noise samples collected during three imaging sequences (GE, FSE and IRSE) has enabled the selection of 14 wavelets of three families (Daubechies, Symlets and Coiflets) to be tested [START_REF] Fokapu | Extracting a reference signal for cardiac MRI gating: Experimental study for wavelet functions choice[END_REF]. (ii) Reference signal extraction: an ECG model where the QRS features have the most significant frequency components is then constructed. Given that most of the QRS power lies in the 3-20 Hz band [START_REF] Thakor | Estimation of QRS complex power spectra for design of a QRS filter[END_REF], the reference signal is created by summing the sub-bands details that cover that interval. Details 6 and 7 were used when the signal was sampled at 1 kHz, covering the 4-16 Hz band. (iii) Trigger production: A simple double thresholding operation on the reference signal then suffices for QRS pulses extraction, thus producing a synchronization signal that gives the cardiac rhythm and can be used to trigger image acquisitions.

Results and Discussion

The algorithm performances were evaluated on 2 sets of ECG signals collected during MRI observations using three imaging sequences: Gradient Echo (GE), Fast Spin Echo (FSE) and Inversion Recovery with Spin Echo (IRSE). 15 signals were recorded on 5 healthy subjects in 1.5T magnet, and 15 other ECGs were collected on mice and rats in 2T and 4.7T magnets. Human ECGs recorded in a 1.5 T and mice ECGs recorded in a 4.7 T showed a prominent MHD effect, while the T wave elevation on mice signals recorded in a 2T magnet was less perceptible (figure 2) TP is the number of true positives, FN the number of false negatives, and FP the number of false positives. The sensitivity Se is the percentage, among QRS beats, of those that were correctly detected by the algorithm. However, the positive predictive value +P represents the percentage, among the QRS detections of those corresponding to real beats. Wavelet performances appeared to be different according to the noise type. For each sequence, some wavelets were more efficient than others in the removal of NMR noise, nevertheless all wavelets were able to greatly reduce the MHD effect. The global algorithm sensitivity is 99.6% and the positive predictive value is 99.9%. Unlike the method proposed in [START_REF] Fischer | Novel real-time R-wave detection algorithm based on the vectorcardiogram for accurate gated magnetic resonance acquisitions[END_REF] where multiple derivations are necessary for eliminating the MHD artifact, only one ECG derivation was used in our study, thus it is easier to place the electrodes on obese people, and less cumbersome for use on small animals. The algorithm performances will further be tested on signals acquired during more complex imaging sequences.

(n) d2(n) d3(n) d4(n) d5(n) d8(n) a8(n) DWT IDWT + ECG d1(k) d2(k) d3(k) d4(k) d5(k) d6(k) d7(k) d8(k) a8(k) d6(n) d7(n) Sref(n) Cardiac trigger

Monitoring parameters

With its real-time implementation, and association with respiratory blanking, this method will establish an interesting solution for synchronization of cardiac MRI [START_REF] Abi-Abdallah | Real-time wavelet-based algorithm for cardiac and respiratory MRI gating[END_REF].

Conclusions

The perturbation of the ECG signals collected during MR imaging greatly obstructs synchronization tasks. Moreover, the higher the requirements for good spacetime or frequency resolution are, the more important the difficulties they induce. Some NMR generated artifacts may overlap the ECG spectrum, and are hard to filter out, such as the MHD induced voltages. These voltages alter the ECG waveform by elevating the amplitude of the T wave and might lead to R peak misdetections. In this work we presented a method that extracts a signal from a contaminated ECG for image acquisition triggering. First the ECG is decomposed into sub-band details using the wavelet transform. A reference signal, containing the necessary information for synchronization i.e. the QRS complexes, is then constructed. Finally a trigger signal is produced by subjecting the reference signal to a double thresholding operation. The simple model proposed in this study has provided a good trigger signal extraction, perfectly synchronized with cardiac events for both human and small rodents. The low error rates have proven its efficiency for eliminating the MHD effect and performing correct R wave detection in the presence of NMR artifacts.

Figure 1 :

 1 Figure 1: QRS pulse generation with the elaborated extraction algorithm.

Figure 2 :

 2 Figure 2: ECG signals recorded in a magnetic field. (a) Human ECG in a 1.5 T magnet (b) mouse ECG recorded in a 2 T magnet (c) mouse ECG recorded in a 4.7 T magnet A total of 30 reference signals were extracted using the 3 wavelet families (Daubechies, Coiflets, and Symlets).Figure 3 illustrates trigger extraction from a human ECG contaminated by the MHD artifact during a FSE sequence. The details d 1 to d 5 contain artefacts due to RF and gradient switching, d 8 contains the MHD artefact, a 8 is mainly composed of the baseline drift and the respiration signal.

Figure 3

 3 illustrates trigger extraction from a human ECG contaminated by the MHD artifact during a FSE sequence. The details d 1 to d 5 contain artefacts due to RF and gradient switching, d 8 contains the MHD artefact, a 8 is mainly composed of the baseline drift and the respiration signal.

  d1

Figure 3 :

 3 Figure 3: (a) 8 level wavelet based sub-band decomposition of a human ECG subjected to the MHD artifact. (b) The extracted reference signal (sum of details 6 and 7) and the computed triggerAccurate QRS detection implies that a very high percentage of the R waves are detected with a low number of false R-waves being detected. In order to accomplish a quantitative performance evaluation of the present QRS detection algorithm two parameters were defined: The sensitivity (Se) and the positive predictive value (+P) given by:


D Abi-Abdallah
email: dima.abiabdallah@utc.fr

A Drochon

V Robin

O Fokapu

REAL -TIME WAVELET-BASED ALGORITHM FOR CARDIAC AND RESPIRATORY MRI GATING

A real time algorithm for cardiac and respiratory gating, which only requires an ECG sensor, is proposed here. Three ECG electrodes are placed in such a manner that the modulation of the recorded ECG by the respiratory signal would be maximal; hence, given only one signal we can achieve both cardiac and respiratory MRI gating. First, an off-line learning phase based on wavelet decomposition is run to compute an optimal QRS filter. Afterwards, on one hand the QRS filter is used to accomplish R peak detection, and on the other, a low pass filtering process allows the retrieval of the respiration cycle so that the image acquisition sequences would be triggered by the R peaks only during the expiration phase.

Introduction

The constantly evolving MR techniques have resulted in an increasing number of cardiac MRI scans in order to gain important insights into the functional and metabolic bases of heart disease. However, MRI observations of a moving organ, such as the heart, require synchronization so that the image acquisition might be accurately combined with its phase motion, thus eliminating the movement related blurring effects. Efficient synchronization techniques must associate cardiac gating with respiratory blanking, where the heart beat governed acquisitions are only enabled during atrest periods (between breaths) to overcome the respiratory motion artefacts [START_REF] Wiesmann | High Resolution MRI with cardiac and respiratory gating allows for accurate in vivo atherosclerotic plaque visualization in the murine aortic arch[END_REF], [START_REF] Cassidy | Assessment of motion gated strategies for mouse magnetic resonance at high magnetic fields[END_REF]. Cardiac gating is usually done by detecting the R peaks on the simultaneously recorded electrocardiogram, which are then used to trigger the consecutive image acquisition sequences. However accomplishing an accurate R peak detection, and therefore a correct synchronization, is often obstructed by the high corruption levels of the ECG signal due to electromagnetic effects, especially when high resolution imaging is performed where high static magnetic fields as well as strong and fast switching gradients are needed. In fact, the signal S(t), recorded in the NMR environment, can be decomposed [START_REF] Felblinger | Restoration of electrophysiological signals distorted by inductive effects of magnetic field gradients during MR sequences[END_REF] as follows: (t) represents the signal to be analyzed (ECG). S flow (t) is induced by the magnetohydrodynamic effect that results mainly in an elevation of the T wave, which might reach amplitudes as high as the R wave, thus leading to misgating. S move (t) is due to patient-related sensor motions in the magnetic field such as respiration, heart beat and voluntary motion, which become more pronounced with increasing static magnetic field strength. S rf (t), a high frequency component resulting from the radio frequency pulses, and S mr (t) generated by the temporal variations of the magnetic field gradients, are spurious signals that often resemble the QRS spikes and cause erroneous detection. In the case of small animal imaging, these artefacts are further aggravated by the exceptionally high heart rate (400-600 beats per minute), the high respiratory rate (30-60 breaths per minute) [START_REF] Cassidy | Assessment of motion gated strategies for mouse magnetic resonance at high magnetic fields[END_REF], and the low amplitude of the ECG signal (just a few millivolts). In order to combine cardiac gating with respiratory blanking, the respiration motion must be detected. This is usually done using motion detectors; however this is not always comfortable for some patients and is cumbersome for very small animals. Nevertheless, the respiration motion information is enclosed within the electrocardiogram as a modulating signal [START_REF]Amplitude demodulation of the electrocardiogram signal for respiration monitoring and compensation during MR examinations[END_REF]: In fact, the respiratory motion of the thorax, heart and diaphragm changes the vectorial orientation of the electric heart axis, resulting in an amplitude modulation of the ECG which is itself a vectorial projection of the electrical activity of the heart. In a previous work [START_REF] Fokapu | Extracting a reference signal for cardiac MRI gating: Experimental study for wavelet functions choice[END_REF] we developed an off-line algorithm for cardiac triggering based on wavelet subband decomposition. A reference signal was extracted from the ECG, considering the sub-bands where the QRS energy was maximal, and was then subjected to thresholding to produce a trigger signal. The algorithm has proven its efficiency even for very low SNR signals.

) ( ) ( ) ( ) ( ) ( ) ( t S t S t S t S t S t S mr rf move flow el      S el
In the present work we exploit the previously elaborated algorithm to produce an adapted filter, for real time use, that yields maximum QRS energy using the most appropriate wavelet. Furthermore, by applying an adequate positioning of the ECG electrodes that maximizes the respiratory dependent modulation, we make use of the respiration information contained in the ECG signal in order to generate, in real time, a synchronization signal that takes into account the respiratory movements as well.

The Algorithm

ECG signals are characterized by a cyclic occurrence of patterns (QRS complexes, P and T waves) with different frequency contents. Power spectral analysis of the ECG show that, P and T waves present a significant spectral density up to 10 Hz only, and that most of the QRS power lies in the 3-20 Hz band, beyond which the complex energy decreases gradually [START_REF] Thakor | Estimation of QRS complex power spectra for design of a QRS filter[END_REF]. Thus, a strategical approach for detecting heartbeats would be to analyze different sub-bands of the ECG [START_REF] Afonso | ECG beat detection using filter banks[END_REF] using wavelet based filter banks. Figure 1 shows an illustration of a previously developed algorithm [START_REF] Fokapu | Extracting a reference signal for cardiac MRI gating: Experimental study for wavelet functions choice[END_REF] that consists in decomposing the ECG signal into multiple scales (8 for a signal sampled at f=1Khz), reconstructing the details and keeping only those that contain the maximum of the QRS energy in order to compose a reference signal. The algorithm was tested with simulated as well as small rodents ECGs that had been contaminated by one of three imaging sequences (Gradient Echo, Fast Spin Echo and Inversion Recovery Spin Echo) and some analyzing wavelets revealed having outstanding performances according to the signal and the sequence type.

Even though the wavelet decomposition algorithm gave proof of an almost infallible QRS detection, it could not be applied in real time for image synchronization, because the decomposition into 8 levels followed by the reconstruction of the 6th and 7th details induces a significant delay. Equation (1) gives the time delay necessary for the reconstruction of the n th detail using a wavelet with M coefficients

   1 1 2 1   M f n (1)
For instance if we use the Db6 wavelet, to reconstruct the 7th detail of a signal sampled at 1000Hz we get a 0,635 seconds delay, which is almost equal to one heart cycle. In order to achieve real time synchronization, we elaborated a two phased algorithm that consists firstly of a training stage where an optimal QRS filter is computed based on wavelet decomposition, followed by a real time trigger generation phase. In the learning phase, we record a portion (~5sec) of signal S ori (t) while the imaging sequence gradients are on, and then process it to produce a reference signal S ref (t) based on the previous method using the most appropriate wavelet. Afterwards, with a Recursive Least Square (RLS) algorithm we compute the coefficients of the filter that would yield the same reference signal if applied to S ori (t) (minimizing the mean squared error, see [START_REF] Haykin | Standard recursive leastsquares estimation[END_REF] for details). Thus we obtain a filter F(n) that is adapted to both the ECG signal and the applied imaging sequence. The second phase, illustrated in figure 2, consists of real time acquisition and processing of the ECG signal: on one hand, the recorded signal is filtered by the calculated F(n) then subjected to thresholding to produce the cardiac trigger; on the other hand a low pass filter 1,2Hz is applied to that same recorded signal to produce the respiratory modulation signal, which then provides the respiration trigger by a simple threshold at the mean value . The final synchronization signal is a combination of both triggers (logical AND) where the cardiac trigger is only taken into account during the low respiratory phase.

Results and Discussion

The method has been applied on anaesthetized mice and rats. The experimental setup incorporates a 2 Tesla OXFORD 85/310 horizontal cryomagnet equipped with a 50 mT/m gradient system, a MR compatible ECG sensor which technical details are provided in [START_REF] Felblinger | Electrocardiogram acquisition during MR examinations for patient monitoring and sequence triggering[END_REF], and a digital-signal-processing system. The ECG active electrodes were placed along the heart axis (45°) in order to obtain maximum modulation. The respiration modulated ECG signals were thus detected during MR scanning according to the three sequences mentioned above. The algorithm was able to generate effective online cardiac and respiratory triggering for all 3 sequences. An illustration is given in figure 3. The images acquired with synchronization revealed a very clear enhancement compared to cardiac MRIs taken without cardiac or respiratory triggering (figure 4 and5) Respiration signal prediction to avoid setbacks: Due to the narrow bandwidth of the low pass filter, the respiratory signal is greatly delayed. However, given that the respiration is somewhat periodic, the delay was compensated by a prediction of the upcoming signal along a period which value is continuously updated.

Adaptive improvement: For images which acquisition takes a relatively long time to complete, better results would be obtained when the filter coefficients are recalculated during the synchronization operation and updated every few minutes, in order to adapt the processing to the signal changes.

Conclusions

Cardiac and respiratory synchronization is essential for achieving good quality cardiac MR images. However the perturbation of the ECG signal acquired in the NMR environment, especially during small animal imaging, greatly obstructs gating tasks. In this work, a method for triggering MR imaging sequences on cardiac and respiration rate is presented. An adapted filtering, based on wavelet sub-band decomposition, is applied to the contaminated simultaneously recorded ECG to rid it of the NMR environment artefacts while emphasizing the QRS peaks. This allows the extraction of a cardiac trigger that is completely synchronized with the heart beats. In addition, based on the fact that the amplitude of the ECG varies according to respiration motions, a demodulation of this "cardiac respiratory" signal yields a respiration curve which is used to build the respiratory triggering signal. A final trigger is then obtained by combining both cardiac and respiratory triggers and used for MRI gating. The algorithm has proven great efficiency for small rodents' cardiac imaging. The trigger extraction was almost flawless and the image quality was greatly enhanced.

Figure 1 :

 1 Figure 1: Reference signal extraction: The contaminated ECG is decomposed into 8 scales using the discrete wavelet transform. The details are then reconstructed and a reference signal is formed by summing details 6 and 7, hence covering the 4-16 Hz sub-band.

Figure 2 :

 2 Figure 2: Real time cardiac and respiratory trigger extraction: The signal is subjected to an adapted filtering to enable QRS detection, as well as a low pass filtering to extract the respiration signal. The cardiac and respiratory triggers are then combined to produce a synchronization signal that controls the image sequence triggering

Figure 3 :

 3 Figure 3: Synchronization signal extraction: Cardiac and respiratory trigger from a mouse ECG recorded during the FSE sequence.

Figure 4 :Figure 5 :

 45 Figure 4: Mouse cardiac MRIs acquired with an IRSE sequence: (a) with no synchronization, (b) with cardiac and respiratory gating


MSc Thanh Trung Le

Multi-Branch Hidden Semi-Markov Modeling for RUL Prognosis

Keywords: RUL, HMM, HSMM, MB-HSMM, predictive maintenance, deterioration modeling, diagnosis, prognosis, FCG

Deterioration modeling and remaining useful life (RUL) estimation of equipment are key enabling tasks for the implementation of a predictive maintenance (PM) policy, which plays nowadays an important role for maintaining engineering systems. Hidden Markov Models (HMM) have been used as an efficient tool for modeling the deterioration mechanisms as well as for estimating the RUL of monitored equipment. However, due to some assumptions not always justified in practice, the applications of HMM on real-life problems are still very limited. To tackle this issue and to relax some of these unrealistic assumptions, this paper proposes a multi-branch Hidden semi-Markov modeling (MB-HSMM) framework.

The proposed deterioration model comprises several different branches, each one being itself an HSMM. The proposed model offers thus the capacity to 1) explicitly model the sojourn time in the different states and 2) take into account multiple co-existing and competing deterioration modes, even within a single component. A diagnosis and RUL prognosis methodology based on the MB-HSMM model is also proposed. Thanks to its multiple branches property, the MB-HSMM model makes it possible not only to assess the current health status of the component but also to detect the actual deterioration mechanism. Based on the diagnostic results, the component RUL can then be calculated.

The performance of the proposed model and prognosis method is evaluated through a numerical study. A Fatigue Crack Growth (FCG) model based on the Paris-Erdogan law is used to simulate deterioration data of a bearing under different operation conditions. The results show that the proposed MB-HSMM gives a very promising performance in deterioration mode detection as well as in the RUL estimation, especially in the case where these deterioration modes exhibit very different dynamics.

MOTIVATION AND PROBLEM STATEMENT

Predictive Maintenance (PM) plays nowadays an important role in the maintenance of engineering systems since it has the ability to recommend in advance proper maintenance actions by based on the relevant information collected from condition monitoring [START_REF] Jardine | A review on machinery diagnostics and prognostics implementing condition based maintenance[END_REF]. Within a PM program, diagnosis and prognosis are two key enabling steps: diagnosis deals with fault detection, isolation and identification as well as health assessment of the monitored system; prognosis aims at predicting the evolution of the fault in the future and estimate the remaining time before a failure occurs. An accurate estimation of the RUL could provide ample time for maintenance engineers to schedule a repair or to acquire replacement components before they actually fail [START_REF] Zhang | An integrated approach to bearing fault diagnostics and prognostics[END_REF]. Due to the stochastic nature of the deterioration mechanisms, stochastic processes are well-suited for modeling the deterioration processes and have shown promising results within the PM framework [START_REF] Si | Remaining useful life estimation-a review on the statistical data driven approaches[END_REF][START_REF] Baruah | HMMs for diagnostics and prognostics in machining processes[END_REF][START_REF] Dong | Equipment health diagnosis and prognosis using hidden semi-Markov models[END_REF][START_REF] Geramifard | A physically segmented hidden Markov model approach for continuous tool condition monitoring: Diagnostics and prognostics[END_REF][START_REF] Dong | A segmental hidden semi-Markov model (HSMM)-based diagnostics and prognostics framework and methodology[END_REF][START_REF] Huynh | Adaptive condition-based maintenance decision framework for deteriorating systems operating under variable environment and uncertain condition monitoring[END_REF].

As doubly embedded stochastic processes, Hidden Markov Models (HMM) have been successfully applied in many domains, such as speech recognition, genes and deoxyribonucleic acid (DNA) analysis, etc. [START_REF] Rabiner | A tutorial on hidden Markov models and selected applications in speech recognition[END_REF] thanks to their strong mathematical basis. In the PM context, HMM classify equipment conditions among several meaningful discrete states, such as "normal", "minor defect", "maintenance required" and "failure" and therefore can give easy-tointerpret results for maintenance implementation [START_REF] Si | Remaining useful life estimation-a review on the statistical data driven approaches[END_REF]. These models can also characterize the stochastic relationship between the features extracted from condition monitoring data and the actual health states of the equipment. For these reasons, HMMs in the recent year are being more and more investigated to be used as an efficient tool for deterioration modeling as well as for RUL estimation [START_REF] Wang | CHMM for tool condition monitoring and remaining useful life prediction[END_REF]. For example, Baruah et al. employed HMMs for carrying out both the diagnosis and the prognosis for metal cutting tools [START_REF] Baruah | HMMs for diagnostics and prognostics in machining processes[END_REF]. In [START_REF] Zhang | An integrated approach to bearing fault diagnostics and prognostics[END_REF], HMMs are used to build a health/degradation index representing the current system health status. This index is then fed into an adaptive prognostic scheme in order to estimate the RUL of a bearing. However, due to its inherent Markov property, the sojourn time in a state of the HMM model is geometrical or exponential distributed, which is not often true in practice. The semi-Markov based models could overcome this problem by allowing arbitrary probability distributions such as Gaussian or Gamma for the state sojourn time, [START_REF] Rabiner | A tutorial on hidden Markov models and selected applications in speech recognition[END_REF]. From this point, several extensions have been adopted in the literature, such as Hidden semi-Markov models [START_REF] Dong | Equipment health diagnosis and prognosis using hidden semi-Markov models[END_REF], segmented HMMs [START_REF] Geramifard | A physically segmented hidden Markov model approach for continuous tool condition monitoring: Diagnostics and prognostics[END_REF], segmented Hidden Semi-Markov Models [START_REF] Dong | A segmental hidden semi-Markov model (HSMM)-based diagnostics and prognostics framework and methodology[END_REF], etc.

To the best of our knowledge, almost all of the Markovbased modeling frameworks for diagnosis and prognosis in the literature deal with the mono-mode case, meaning that only one degradation mechanism is allowed to take place at one time. For example, the authors in [START_REF] Dong | Equipment health diagnosis and prognosis using hidden semi-Markov models[END_REF] use different HSMMs to represent different health states of a pump but these states belong to only one deterioration mode. Or in [START_REF] Geramifard | A physically segmented hidden Markov model approach for continuous tool condition monitoring: Diagnostics and prognostics[END_REF], only one degradation mechanism of a milling machine is considered and modeled by one segmented HMM. In real cases, however, several degradation modes can co-exist in competition even within a single component. For instance, different internal elements of bearings can be subjected to degradation due to different operating conditions, leading to different rates in deterioration evolution of the bearings. To take into account such phenomena, we propose a semi-Markov based structure called multi-branch Hidden semi-Markov Model (MB-HSMM). A simulation study is carried out in order to evaluate the effectiveness of the proposed model with respect to a standard HSMM one.

The remainder of the paper is organized as follows. Section 2 is devoted to review some basic theories of HMM and HSMM as well as the extension to MB-HSMM model. Section 3 describes a framework based on the proposed model for diagnosis and prognosis. To illustrate the methodology, a comparative study is introduced in the section 4. Conclusion follows in Section 5.

THEORETICAL BACKGROUND

Basic theory of Hidden Markov Model

Hidden Markov Model consists of an underlying Markov chain which is not directly observable but can be revealed through observation [START_REF] Rabiner | A tutorial on hidden Markov models and selected applications in speech recognition[END_REF]. Within the deterioration modeling framework, the hidden states represent the health states of equipment, while the observations can be either measurable raw data such as vibration signals or high-level features extracted from condition monitoring data. The relation between the hidden states and the observations is often represented by probabilistic models.

Mathematically, the complete specification of an HMM consists of the following elements [START_REF] Rabiner | A tutorial on hidden Markov models and selected applications in speech recognition[END_REF]  N , the number of states in the model  Q , the number of distinct observation symbols per state  A finite set of hidden states, i.e.,

  1 2 , ,..., N S S S S   A state transition probability distribution, i.e.,   ij A a  , where   1 ij t j t i a P q S q S     |  An observation symbol probability distribution, i.e.,     j B b k  where     j k t j b k P Y q S   | , 1 k Q    An initial state probability distribution, i.e.,   i p p  where   i t i p P q S   , 1 i N   A complete HMM is often represented by a compact notation   , , p A B   for convenience.
There exist several types of HMMs, such as ergodic or left-right (Bakis) HMMs [START_REF] Rabiner | A tutorial on hidden Markov models and selected applications in speech recognition[END_REF]. Since the deterioration processes are often irreversible in practice, meaning that equipment cannot come back from a worse health state to a better one as time progresses, the left-right topology is chosen for the deterioration modeling purpose in this paper.

In the framework of machine condition monitoring, measures are usually continuous signals, such as the vibration amplitudes, motor current, etc. The observation sequences should therefore be modeled by continuous probability density functions (pdf) to obtain a better performance in the diagnosis and the prognosis [START_REF] Wang | CHMM for tool condition monitoring and remaining useful life prediction[END_REF]. For this reason, the Gaussian distribution will be adopted for the observations in this paper.

Hidden semi-Markov Model (HSMM)

In standard HMM models, due to the Markovian assumption the sojourn time in a state is either geometrical distributed in the discrete time case or exponential distributed in the continuous time case. As argued in [START_REF] Russell | Explicit modeling of state occupancy in hidden Markov models for automatic speech recognition[END_REF], this can be a source of inaccurate duration modeling with the HMMs since most real-life systems do not exhibit this property. To overcome this problem, the underlying Markov process can be replaced with a semi-Markov one resulting in a more general model called Hidden semi-Markov Model (HSMM).

Unlike the standard HMM, the state sojourn time in HSMMs can follow any arbitrary distribution, such as Gamma or Gaussian, etc. Once entering a state, the process stays in the state for a period of time whose length is determined by the arbitrary distribution. At the end of this period, the process randomly moves to another state according to the state transition probability matrix A . Clearly, HMMs are a special class of HSMMs.

As in the previous section, the underlying semi-Markov process in the HSMM is assumed to be strictly left-right and the observations are Gaussian distributed. Furthermore, the sojourn time in a state is also modeled by continuous Gaussian distributions. Similarly to the standard HMM, an HSMM can be characterized by the following parameters: the initial state distribution p , the transition model A , the state duration distribution denoted by D , and the observation model B . Thus, an HSMM can be written as

 

, , , p A D B   .

Multi-branch Hidden semi-Markov model (MB-HSMM)

As already stated, standard HSMMs are only suitable for modeling one deterioration mode at a time. To deal with the co-existence of multiple competing deterioration modes, the concept of multi-branch model has recently been introduced [START_REF] Le | Hidden Markov Models for diagnostics and prognostics of systems under multiple deterioration modes[END_REF]. However, due to the Markovian assumption in the model, limitation in the sojourn time modeling still exists. In this section, we extend the model in [START_REF] Le | Hidden Markov Models for diagnostics and prognostics of systems under multiple deterioration modes[END_REF] to a so called multibranch HSMM for modeling multiple deterioration mechanisms of equipment. An example of a left-right multibranch HSMM is illustrated in Figure 1.

The proposed model consists of several branches, each one representing one deterioration mode. Some important assumptions have been made in constructing the model. Firstly, it is assumed that once the degradation has started following a given path, the equipment follows this deterioration mode until it reaches the end of life, i.e. there is no allowed switch between branches. The second assumption is that, the normal condition can last very long with no abnormal sign in the observations. The observations emitted in this state are not hence of interest in the framework of deterioration modeling. Furthermore, once the equipment has failed, it is stopped and no more measures can be obtained. From these points, the initial 0 S as well as the failed f S states are considered as dummy states with no emission of observation in our model. Therefore, each branch can be seen as a single HSMM model and is assigned an a priori probability, Figure 1.

  / k k k P p K K    , 1 1 M k k p    (2)
where k  denotes the constituent HSMM for branch k , k K is the number of training data sequences corresponding to the mode k , K is the number of the total training data sequences and M is the total number of the branches.

MB-HSMM FOR DIAGNOSIS AND PROGNOSIS

In this section, a two-phase framework on the MB-HSMM model is proposed for the diagnosis and the prognosis (c.f. Figure 2). It consists of two phases: off-line and on-line.

Off-line phase

This phase aims to learn an MB-HSMM model from the historical data. The idea is to train each branch individually and then combine them to obtain a final model. To do this, we firstly divide the training data set into M different groups in which each one corresponds to a deterioration mode. In this study, we suppose that the number of modes M is known beforehand thanks to the expert knowledge. Clustering techniques such as "kmeans" or "k-NN", etc. can hence be employed for the classification. Each data group is then used to train a constituent HSMM by using the Baum-Welch algorithm [START_REF] Rabiner | A tutorial on hidden Markov models and selected applications in speech recognition[END_REF] in which the forward-backward (FB) procedure is the key process for parameters learning and likelihood estimation. We adopt the FB procedure proposed by Yu and Kobayashi in [START_REF] Yu | Practical implementation of an efficient forward-backward algorithm for an explicitduration hidden Markov model[END_REF], which stands out to be the most efficient one.

After having trained the constituent branches, the a priori probabilities can be estimated by Equation ( 2) and a full MB-HSMM model is obtained.

On-line phase

In this phase, the trained model is used to assess the actual health state of the monitored equipment as well as to estimate its RUL. Due to the co-existence of multiple deterioration mechanisms, one must firstly identify the underlying one that is occurring. Indeed, wrong mode detection can lead to a wrong state assessment and inaccurate RUL estimation as well. In order to minimize the misclassification rate, the branch having the maximum a posterior probability given the data can be considered to correspond to the actual deterioration mechanism [START_REF] Bishop | Pattern recognition and machine learning[END_REF], that is:

  ˆarg max k k k P   O | (3)
Using Bayes' theorem, these probabilities can be given by:

          1 k k k M k k k P P P P P         O O O | | | (4)
where     k P  is calculated by Equation (2). After having identified the deterioration mode, the next diagnosis stage is to recognize the actual health state of the equipment. This can be carried out thanks to the Viterbi algorithm [START_REF] Rabiner | A tutorial on hidden Markov models and selected applications in speech recognition[END_REF]. Given the detected mode, the most probable single state sequence (path) is the one that gives the maximum joint probability of the path and the observations:

  ˆˆâ rg max , Q k k k Q P Q    O | (5)
where k Q represents a possible path under the mode k . The actual health state of the equipment is considered as the final one of the path * Q .

RUL calculation

Given the deterioration mode and the current health state estimated from the diagnostic stage, the RUL can be estimated in the same way as in a single HSMM case. However in the MB-HSMM framework, this estimation depends strongly on the mode detection results. Indeed, detecting the wrong mode may lead to a large bias in the RUL estimation. For example, if a crack is propagating with a high rate, the wrong detection of a low deterioration mode will lead to a quite high RUL value with respect to the actual one. This phenomenon often occurs at the beginning of the defect propagation since observations are insufficient to show an obvious trend of the deterioration. In this study, the Bayesian Model Averaging (BMA) technique [START_REF] Hoeting | Bayesian model averaging: a tutorial[END_REF] is adopted to tackle the problem. The RUL distribution is calculated as an average of the posterior distributions under constituent modes, weighted by their posterior mode probabilities:

      1 RUL RUL , M k k k P P P      O O O | | | (6)
where  

k P  O | is given from Equation (4).
The procedure can now turn to the estimation of the RUL for each individual constituent HSMM model, which is straightforward. Assuming that equipment goes through degraded states [START_REF] Jardine | A review on machinery diagnostics and prognostics implementing condition based maintenance[END_REF]

RUL N t t i i j j i D D      (7)
Intuitively, the residual time is a conditional random variable |

t i i i i i D D D D D   
where i D is the time that the equipment has spent in the state i S which can be calculated thanks to the Viterbi algorithm from the diagnostic stage. Given that i D is Gaussian distributed, e.g.

 

,

i i i D     , it
can be deduced that the residual time | 0

t t i i D D  follows a
Normal distribution with mean i i D   and standard deviation i  but truncated to the left of zero.

Besides that, due also to the Gaussian assumption, the sum element in Equation ( 7) is Normal distributed. Denote

1 N j j i Z D     , we have   , z z Z     where 1 N z j j i       and 2 1 N z j j i    
  . The RUL becomes the summation of a truncated Normal distribution t i D and a Normal distributed variable Z . The cumulative distribution function of the RUL is therefore given by [START_REF] Lipow | Query 2: The sum of values from a normal and a truncated normal distribution[END_REF]:

      1 1 x a RUL x u a F x u du a                                         (8)
where

i i z D a      , i z     and   2 /2 1 2 e      
is the probability density function of the standard normal distribution and     is its cumulative distribution function. After evaluating all the individual branches, the RUL estimation under the MB-HSMM can be accomplished by Equation [START_REF] Baruah | HMMs for diagnostics and prognostics in machining processes[END_REF].

NUMERICAL EXAMPLES

To evaluate the proposed MB-HSMM-based diagnosis and prognosis performance, the data representing several different deterioration paths have been generated by the Fatigue Crack Growth (FCG) model. The FCG model is chosen since it has been widely used in the literature to describe the propagation process of a crack within a bearing [START_REF] Huynh | Adaptive condition-based maintenance decision framework for deteriorating systems operating under variable environment and uncertain condition monitoring[END_REF]. The bearing is supposed to operate under different operation conditions so that the crack can propagate at different rates, which represent different deterioration modes.

Simulated fatigue crack growth data

In this section, we use the discretized and randomized FCG model in [START_REF] Huynh | Adaptive condition-based maintenance decision framework for deteriorating systems operating under variable environment and uncertain condition monitoring[END_REF]:

  1 1 t i e i i i n w t t b t x x e C e x t        (9)
where i t

x represents crack depth at time i t , C and n are constants depending on the material property, the exponential element t i w e helps to take into account the stochastic aspect of the crack propagation, b  is the base stress level of component, the values

    determine the level of extra stress linking with the state e of the environment.

Since it is often difficult to measure directly and accurately the crack depth in practice, an observation model has been adopted to represent the relationship between the actual values and the measure ones:

Diagnosis and prognosis results

In this section, the bearing is supposed to operate under two different environments and therefore cracks propagate possibly at two different rates. A two-branch HSMM model (cf. figure 1) is implemented to represent these two modes.

The following parameters of the FCG model are chosen:

0.005 C  , 1.3 n  , 1 b   , 1.7 w   , 1 t   , 2 10    ,   0 0.75 T e   , 100 L  , 1 2 0.5 p p  
where L represents the critical level beyond which the bearing is supposed to fail and the superscript T denotes the transpose of a vector. For training purpose, 100 data paths are generated (Figure 3). Based on the times to failure information, this data set can be separated in two different "clusters" thanks to the k-means algorithm. The a priori probabilities for the branches are then calculated (c.f. Equation ( 2)) and we obtain: To determine the number of states for each branch, the Bayesian Information Criterion (BIC) can be adopted as detailed in [START_REF] Le | Hidden Markov Models for diagnostics and prognostics of systems under multiple deterioration modes[END_REF]. The candidate model having the minimum value of BIC will be chosen. In this study, the value 8 N  is obtained for every constituent model.

After having determined the number of states, each model is trained thanks to the Baum-Welch algorithm. An MB-HSMM is obtained and the offline phase is terminated.

We now move to the on-line phase. One data curve serving as the test data is also generated by the FCG model. In this study, the diagnosis and the prognosis procedures are repeated after each 50h as time progresses. The mean values of the estimated RUL associated with the 95% confidence interval are represented in Figure 4. It can be seen that at early time instances, the lack of information due to the limited observations results in the bias and large variance in the RUL estimation. Even though, the real RUL values always lie within the 95% confidence interval. As the time passes, the width of the confidence interval significantly decreases and the estimated RUL values converge to the ones. This demonstrates the accuracy of the proposed RUL estimation method.

Multi-branch vs average HSMM

An interesting question naturally rises in evaluating the advantage of the proposed approach: between the MB-HSMM and an "average" HSMM, which one does give the better RUL prediction results given the same data set. To answer such question, a comparative study is conducted in this section with 100 test deterioration paths. The prognostic performance is evaluated by the root mean squared error (RMSE) metric, where the errors are the differences between the mean estimated RULs and the actual values.

For notation convenience, we denote the average HSMM by AVG-HSMM. The training data set in the previous section is not grouped and is used to train only one AVG-HSMM. The trained model is then used in the on-line phase to estimate the RUL. The test data is the same as in the previous section. for the AVG-HSMM and the MB-HSMM models respectively. It can be concluded that the proposed model give a better performance in RUL estimation than the average model in this case. Another question raises naturally: does the above conclusion still hold in case that the two propagation rates are similar, or in case that they are quite different? To answer this question, it is necessary to define a "distance" between two deterioration modes. In this study, the mode "distance" is defined by the difference between the propagation rates in the FCG model. From Equation ( 9), it can be seen that the propagation rate of a crack is directly proportional to e  : the greater the value of e  is, the more quickly the crack propagates. Therefore, by fixing 1  and varying 2  , the mode distance can be changed. Figure 5 shows that at the beginning when the two deterioration modes are similar, the multi-branch model gives a less accurate result than the one obtained with the average model. In this case, the simple average model is well-adapted to the data. Conversely, the multi-branch model may appear overly complex, leading to over fitting and poor predictive performances. However, once the mode distance increases, the MB-HSMM becomes more and more relevant and its RMSE decreases while it increases with AVG-HSMM. Hence, it can be concluded that the larger the distance between deterioration modes the better the performance in RUL estimation given by the proposed MB-HSMM approach, as compared to an "average" HSMM.

CONCLUSION

The present paper proposes a novel model called multibranch Hidden semi-Markov Model (MB-HSMM) to deal with the deterioration modeling problem in case that there are several competing modes existing in the deterioration mechanisms of an equipment. A diagnostic and prognostic framework based on the proposed model is also given. The actual health state recognition and the RUL estimation can be implemented in the same way as in the traditional HSMM approach. However, by using different HSMM models for different modes, together with the a priori probabilities of each model, the MB-HSMM shows a very promising result in dealing with the multi-mode deterioration mechanism.

Future research will be focused on the extension of the MB-HSMM to allow the switching between the states of the different branches. The proposed model and the extension should also be validated on the data of real systems.

Figure 1 -

 1 Figure 1 -A left-right multi-branch HSMM model.From a training data set, these a priori probabilities can be estimated by:

Figure 2 -

 2 Figure 2 -MB-HSMM framework for diagnosis and prognosis

  of the model k  given the data O and can be calculated through the forward-backward algorithm[START_REF] Yu | Practical implementation of an efficient forward-backward algorithm for an explicitduration hidden Markov model[END_REF].

Figure 3 -

 3 Figure 3 -Two-mode training data

Figure 4 -

 4 Figure 4 -RUL estimation at different time instances.

Figure 5 -

 5 Figure 5 -RMSE at different mode distances.

  and is in the state i S at the current time t . Since the model is strictly left-right, the RUL can be calculated as the summation of two terms: the residual time of staying in the current state and the duration for the equipment to stay in the future health states before entering the failure one. Denote

	2 , ,..., N S S S before reaching the failure state
	i D and t D the random variables representing respectively i S for the total and the residual sojourn time in the state i 1,2, i N   , we have:
	1

f S

ACKNOWLEDGEMENTS

This work is done in the framework of the SUPREME ("SUstainable PREdictive Maintenance for manufacturing Equipment") project, funded by the European Commission in FP7 Program, under the Factories of the Future PPP.