# Bearing Remaining Useful Life Prediction Based on a Nonlinear Wiener Process Model

Authors: Juan Wen; Hongli Gao; Jiangquan Zhang
Journal: Shock and Vibration
 (2018)

---

## Abstract

Prognostic is an essential part of condition-based maintenance, which can be employed to enhance the reliability and availability and reduce the maintenance cost of mechanical systems. This paper develops an improved remaining useful life (RUL) prediction method for bearings based on a nonlinear Wiener process model. First, the service life of bearings is divided into two stages in terms of the working condition. Then a new prognostic model is constructed to reflect the relationship between time and bearing health status. Besides, a variety of factors that cause uncertainties toward the degradation path are considered and appropriately managed to obtain reliable RUL prediction results. The particle filtering is utilized to estimate the degradation state, qualify the uncertainties, and predict the RUL. The experimental studies show that the proposed method has a better performance in RUL prediction and uncertainty management than the exponential model and the linear model.


---

## Body

## 1. Introduction

The condition-based maintenance (CBM) is a widely used maintenance policy which schedules the maintenance for components or systems according to the condition monitoring data [1]. Generally, the CBM focuses on two aspects of work: diagnostics and prognostics [2–6]. The condition monitoring data, such as vibration, current, acoustic emission, and temperature, can be utilized to implement the CBM. Diagnostics attempt to detect and identify faults in a specified system, while prognostics are detecting changes in system states and predicting the remaining time before the occurrence of failure [1, 7]. Therefore, reasonable maintenance policies often schedule maintenance plans according to prognostics results to avoid unnecessary maintenances or replacements, which can significantly reduce the overall life-cycle costs and increase the reliability and availability of mechanical systems. As a result, more and more importance is attached to the prediction of remaining useful life (RUL), and it triggers a certain amount of research in this area.Roughly speaking, the prognostics approaches for mechanical equipment can be classified as data-driven methods and model-based methods [8]. For data-driven methods, no prior knowledge about systems is needed, and the relationship between the RUL and the historical failure data is constructed with machine learning techniques. Artificial neural network (ANN) is one of the most commonly used tools. Tian et al. [9] proposed an ANN based prognostics method using both failure data and suspension data. Besides, support vector machine (SVM) is also a widely selected technique for prognostics. Tran et al. [10] forecasted the RUL through the SVM and time-series methods. However, data-driven methods do not use the useful information among the degradation process, and they are restricted to the quantity and quality of the collected data.Differently, model-based methods deal with the problem of prognostics by building mathematical models which can describe the deterioration process. The model-based approaches can be generally divided into two categories, physics model-based methods and statistical model-based methods [11]. Physics model-based methods construct the models based on the failure mechanisms of equipment [12]. The Paris law is one of the most commonly used physics model in describing the degradation process of machinery, and many investigations employed it and its variants for RUL prediction [6, 7, 13]. Besides, other physics models were also used in prognostics, such as the Forman law and the Norton law [14, 15]. The physics model-based techniques can forecast the RUL accurately when there is enough knowledge about the physics of damage. However, it is difficult to completely understand the failure mechanisms of some complex systems and new equipment.For statistical model-based methods, they often build stochastic models based on empirical knowledge to achieve machinery prognostics [16]. A lot of work has been done on how to estimate the RUL with stochastic models. Lawless and Crowder [17] constructed a gamma process model to fit crack growth data. Wang and Xu [18] used the inverse Gaussian process model to process the degradation data and estimate the parameters with the expectation maximization algorithm. Besides, the Wiener process models were also widely used in RUL prediction. There are three different Wiener process models, i.e., the linear model, the exponential-like model, and the nonlinear model. Investigations of the linear model for RUL prediction can be found in [19–21]. However, most mechanical systems often experience nonlinear deterioration process, which limits the application of the linear model. The exponential-like model is a special nonlinear model, which can be transformed into the linear model. Gebraeel et al. [22] put forward an exponential degradation model for bearing RUL prediction. This model was extensively used and improved by many researchers [23–26]. However, it is restricted to the exponential degradation process and cannot describe other nonlinear degradation process. As a result, the nonlinear model is used to make up for this deficiency. Si [27] presented an adaptive prognostics method based on the nonlinear model and applied it to battery RUL prediction. However, these models are usually not suitable for practical mechanical systems whose deterioration process is more complicated than the battery degradation.In this study, we attempt to predict the RUL of bearings with a Wiener process based nonlinear model. Generally, there are three stages during the entire life of bearings, i.e., the normal stage, the degradation stage, and the failure stage [28]. In the failure stage, the bearing is out of order and should be replaced to guarantee the safety of the whole mechanical system. Consequently, we focus on the former two stages. In the normal stage, the health status of bearings is relatively stable, and the condition monitoring data rarely change during this period. In contrast, the condition monitoring data vary a lot from time to time in the second stage. To predict the RUL accurately, it is necessary to manage different stages differently. For degradation modeling with switches in different stages, Si et al. [29, 30] and Zhang et al. [31] focused on this topic and provided applications for the prognostics of gyroscopes, the inertial navigation system, and li-ion batteries. Based on the previous work, in this paper, we attempt to properly separate the bearings’ two working stages and model the degradation process to avoid unnecessary computational cost and obtain accurate RUL prediction results.In addition, no model can perfectly describe the stochastic degradation process of a specific system due to various uncertainties resulting from the materials, load, environment, and sensors [32, 33]. Uncertainties due to different factors can be reflected by the model parameters, and, by updating the parameters with the real-time observations, reliable prognostics results can be obtained. In general, it is more appropriate to include all the parameters related to various uncertainties as multiple states in state space modeling and update them simultaneously.Many investigations were conducted to reduce the uncertainties of the deterioration process in structure reliability analysis and health monitoring [34, 35]. An et al. [36] studied the influence of noise and bias in updating the parameters of Paris law based model. Sankararaman et al. [37] qualified the uncertainties with first-order reliability methods. Sun et al. [7] estimated the parameters related to the Paris law, the process noise, and the measurement noise.Regarding the work of uncertainty management in the mechanical area, Zhao et al. [6] used the Paris law to describe the development of the gear crack and updated the two correlated parameters in the model. But it was a simulation based framework and the uncertainties resulting from process noise and measurement error were neglected. Leiet al. [13] used a Paris law based model to predict the RUL of bearings, but only one parameter was updated simultaneously with the health state. As for previous studies about prognostics based on the nonlinear Wiener process, the drift coefficient was the only state in the state space model [27, 38, 39]. Consequently, another objective of this work is to consider all the parameters related to the nonlinear model during the state estimation process. Specifically, the particle filtering (PF) algorithm is utilized to incorporate the model and measurements for parameters estimation as well as RUL prediction. Unlike Kalman filtering techniques, the PF is not restricted to the Gaussian assumption, and it can also deal with nonlinear scenarios.This study proposes a new Wiener process based RUL prediction approach for bearings, and the major contributions are summarized as follows.(1)
The service life of bearings is divided into two stages, i.e., the normal stage and the degradation stage, and the two stages are managed differently during the prognostics. Based on the health stage division, an improved Wiener process based nonlinear model is developed to describe the degradation path of the degradation stage, and various uncertainties are included in this model.(2)
Besides the drift coefficient, other parameters in the nonlinear model are also considered as hidden states in state space modeling. With the PF, all the parameters are updated simultaneously.The remainder of this paper is organized as follows. The Bayesian inference and the PF algorithm are briefly introduced in Section2. The proposed method for bearing prognostics is illustrated in Section 3. In Section 4, a numerical example and an application in bearing prognostics are provided to manifest the effectiveness of the presented approach. Conclusions are drawn in Section 5.


## 2. Theoretical Background of the PF

### 2.1. Bayesian Inference

In general, a dynamic system can be represented by a process model (1) and an observation model (2) as follows:(1)zk=fzk-1,ωk-1,(2)sk=hzk,τk,wherezk andsk represent the actual states and the corresponding measurements at current timetk, respectively;f is the state transition function, reflecting the relationship between the states attk and the states at tk-1, whereash is the observation function, indicating the relationship betweensk andzk; ωk-1 is the i.i.d. (independent and identically distributed) process noise, and τk is the i.i.d. measurement noise.With the application of Bayesian inference, the posterior probability density function (PDF)p(z0:k∣s1:k) of the states z0:k={z0,…,zk} fromt0 totk can be estimated with a series of observations s1:k=s0,…,sk. Specifically,p(z0:k∣s1:k) can be calculated recursively by the following equations:(3)pzk∣s1:k-1=∫pzk∣zk-1pzk-1∣s1:k-1dzk-1,(4)pzk∣s1:k=psk∣zkpzk∣s1:k-1psk∣s1:k-1,wherep(zk∣s1:k-1) is the prior PDF defined by the process model,p(sk∣zk) is the likelihood determined by the observation model, andp(sk∣s1:k-1) is the evidence, which can be described as follows:(5)psk∣s1:k-1=∫psk∣zkpzk∣s1:k-1dzk.


### 2.2. The PF Algorithm

Generally, it is impractical to solve (3) and (4) analytically. Hence, the PF is utilized to approximate the posterior PDF by a group of particles with corresponding weights [40, 41]. Specifically, the estimated posterior PDF can be denoted as follows:(6)pzk∣s1:k≈∑i=1Nwkiδzk-zki,where zkii=1N represents the particles, wkii=1N are the associated weights, N is the number of particles, and δ(·) stands for the Dirac function.Particleszkii=1N are sampled from an importance PDF q(zk∣s1:k), and their weights can be calculated by the following: (7)wki=wk-1ipsk∣zkipzki∣zk-1iqzki∣zk-1i,s1:k.Usually, the importance PDF is defined asq(zk∣s1:k)=p(zk∣zk-1) for an easy implementation. Therefore, (7) becomes the following: (8)wki=wk-1ipsk∣zki.The procedures of standard PF are concluded as follows [42].(  1) Initialization. Setk= 0, and sample particles z0ii=1N from the initial distribution p(z0). Also, initialize the weight of each particle as w0i =1/N.(  2) Importance Sampling. Set k=k+1, calculate the prior PDF p(zk∣zk-1) with (1), and draw zkii=1N from p(zk∣zk-1).(  3) Weights Updating. Update the weights with the newly obtained measurement with (7), and normalize the weights by(9)wki=wki∑i=1Nwki.(  4) Resampling. Remove particles with small weights and copy those with large weights. The implementation of the resampling is shown in detail as follows [43]:(a)
Setd=1:N, and generate a random data ud from the uniform distribution U(0,1).(b)
Setj=1:N, and calculate the cumulative distribution function (CDF) of weights as ∑i=1jwki. If ∑i=1jwki≥ud, duplicate zkj as a new particle z^kd with a weight of 1/N, and go back to Step a; otherwise, set j=j+1, and return to Step b.(  5) State Estimation. Estimate the current states with the resampled particles and the associated weights:(10)z^k=1N∑l=1Nz^kd.Then turn to Step 2 and repeat Steps 2-5 for the next inspection time.


## 2.1. Bayesian Inference

In general, a dynamic system can be represented by a process model (1) and an observation model (2) as follows:(1)zk=fzk-1,ωk-1,(2)sk=hzk,τk,wherezk andsk represent the actual states and the corresponding measurements at current timetk, respectively;f is the state transition function, reflecting the relationship between the states attk and the states at tk-1, whereash is the observation function, indicating the relationship betweensk andzk; ωk-1 is the i.i.d. (independent and identically distributed) process noise, and τk is the i.i.d. measurement noise.With the application of Bayesian inference, the posterior probability density function (PDF)p(z0:k∣s1:k) of the states z0:k={z0,…,zk} fromt0 totk can be estimated with a series of observations s1:k=s0,…,sk. Specifically,p(z0:k∣s1:k) can be calculated recursively by the following equations:(3)pzk∣s1:k-1=∫pzk∣zk-1pzk-1∣s1:k-1dzk-1,(4)pzk∣s1:k=psk∣zkpzk∣s1:k-1psk∣s1:k-1,wherep(zk∣s1:k-1) is the prior PDF defined by the process model,p(sk∣zk) is the likelihood determined by the observation model, andp(sk∣s1:k-1) is the evidence, which can be described as follows:(5)psk∣s1:k-1=∫psk∣zkpzk∣s1:k-1dzk.


## 2.2. The PF Algorithm

Generally, it is impractical to solve (3) and (4) analytically. Hence, the PF is utilized to approximate the posterior PDF by a group of particles with corresponding weights [40, 41]. Specifically, the estimated posterior PDF can be denoted as follows:(6)pzk∣s1:k≈∑i=1Nwkiδzk-zki,where zkii=1N represents the particles, wkii=1N are the associated weights, N is the number of particles, and δ(·) stands for the Dirac function.Particleszkii=1N are sampled from an importance PDF q(zk∣s1:k), and their weights can be calculated by the following: (7)wki=wk-1ipsk∣zkipzki∣zk-1iqzki∣zk-1i,s1:k.Usually, the importance PDF is defined asq(zk∣s1:k)=p(zk∣zk-1) for an easy implementation. Therefore, (7) becomes the following: (8)wki=wk-1ipsk∣zki.The procedures of standard PF are concluded as follows [42].(  1) Initialization. Setk= 0, and sample particles z0ii=1N from the initial distribution p(z0). Also, initialize the weight of each particle as w0i =1/N.(  2) Importance Sampling. Set k=k+1, calculate the prior PDF p(zk∣zk-1) with (1), and draw zkii=1N from p(zk∣zk-1).(  3) Weights Updating. Update the weights with the newly obtained measurement with (7), and normalize the weights by(9)wki=wki∑i=1Nwki.(  4) Resampling. Remove particles with small weights and copy those with large weights. The implementation of the resampling is shown in detail as follows [43]:(a)
Setd=1:N, and generate a random data ud from the uniform distribution U(0,1).(b)
Setj=1:N, and calculate the cumulative distribution function (CDF) of weights as ∑i=1jwki. If ∑i=1jwki≥ud, duplicate zkj as a new particle z^kd with a weight of 1/N, and go back to Step a; otherwise, set j=j+1, and return to Step b.(  5) State Estimation. Estimate the current states with the resampled particles and the associated weights:(10)z^k=1N∑l=1Nz^kd.Then turn to Step 2 and repeat Steps 2-5 for the next inspection time.


## 3. The Proposed Method

The framework of the proposed method is illustrated in Figure1, which consists of three parts: health stage division, Bayesian inference, and RUL prediction. At the beginning, data collected from the normal stage is used to determine the deterioration threshold (DT). Then newly obtained measurements are used to detect the degradation. Once the bearing begins to deteriorate, the PF is employed to update model parameters and estimate the bearing’s health status. For a specific bearing, its health state and model parameters can be initialized by failure history data or experience. At each inspection time in the degradation stage, the updated health state and parameters can be input into the degradation model for RUL prediction in the next part. The update and prediction process will repeat until the estimated health state exceeds a prespecified failure threshold (FT).Figure 1
Flowchart of the proposed method.


### 3.1. Health Stage Division

To conduct the proposed method, an important task is to divide the health stages precisely. The statistical feature Root Mean Square (RMS) is used as the health indicator (HI) of bearings in this study due to its sensitivity to the global damage of bearings [44]. To detect the degradation, the HIs obtained from the healthy state are used to determine the DT. Specifically, the mean μn and the standard deviation σn of these HIs are calculated, respectively, and μn + 3σn is defined as the DT in this study. Then the health state of bearings is monitored by comparing the DT with the newly obtained HI.However, there are some abnormal values due to random noises, which could mislead the recognition of health states. To solve this problem, some measures are taken to find the start point of degradation, as shown in Figure2. The process of the degradation detection mechanism is described as follows:(1)
At inspection timetk, calculate the RMS of the new measured data.(2)
Compare the RMS attk (HIk) with the DT (μn + 3σn). If HIk > μn + 3σn, definetk astm; otherwise, setk =k + 1, return to Step 1.(3)
Compute the RMS attm+1 andtm+2, respectively. If these two values of RMS satisfy HIm+q-μn>3σnq=1,2, considertm as the beginning point of the degradation stage and conduct RUL prediction; otherwise, setk =k + 1, and return to Step 1 and repeat Steps 1-3.Figure 2
The degradation detection mechanism.Practically, the RMS of a bearing may exceed the DT because of random noises rather than the degradation. However, these abnormal scenarios are difficult to arise three times consecutively in the normal stage. Therefore, we define the bearing begins to degrade when three consecutive values of RMS are bigger than the DT.


### 3.2. Development of the Bearing Degradation Model

There are four major variability sources related to the uncertainties of machinery deterioration: the temporal variability, the measurement variability, the unit-to-unit variability, and the nonlinear variability [39]. To manage the uncertainties in a degradation process appropriately, it is necessary to consider all these factors. Accordingly, the degradation process of bearings can be illustrated by a nonlinear form as follows:(11)St=a+λtb+σBt,whereS(t) is the health state of a bearing at timet,a represents the initial health state of the degradation stage;λtb denotes the degradation process, λ reflects the degradation speed, andb represents the nonlinear characteristic;B(t) is the standard Brownian motion (BM), and σ is the diffusion coefficient.For convenience, (11) is transformed into(12)Xt=St-a=λtb+σBt.Accordingly, we can obtain the following difference equation of the degradation process:(13)xk=xk-1+λk-1tkbk-1-tk-1bk-1+ηk,where ηk=σ(B(tk)-B(tk-1)) follows a normal distribution N(0,σ2Δtk), and Δtk=tk-tk-1.In practice, vibration signals are the most commonly selected condition monitoring data, for the advantage that it can well indicate the degradation process of machinery. However, there exists a difference between the actual health state and the HI extracted from the vibration signal due to the noises from the sensors, acquisition systems, and data processing. Thus the measurement variability should be considered, and the observed health state can be described as follows [39]:(14)Zt=St+ν,whereZ(t) represents the HI extracted from the vibration signal at timet, and ν is the measurement noise, which follows a normal distribution N(0,γ2). Also, (14) can be derived as follows:(15)Yt=Zt-a=Xt+ν.In this paper, we useΘ = [λ,b, σ, γ] to denote the unknown model parameters.


### 3.3. RUL Prediction

The RULLk at time tk can be defined as follows:(16)Lk=inf⁡lk:Xtk+lk≥κ∣Y1:k,Θ,where inf⁡· stands for the inferior limit of a variable,lk is the time betweentk and the failure time, and κ is the predefined FT;Y1:k represents the available observations from the initial point of degradation stage to timetk, and the estimated health state at timetk+lk is denoted as follows:(17)Xtk+lk=xk+λtk+lkb-tkb+σBtk+lk-Btk.According to the dependent increment properties of the standard BM,B(tk+lk)-B(tk) is also a standard BM. Thus the predicted PDF of RUL at timetk can be illustrated as follows [27]:(18)fLk∣xk,Θlk∣xk,Θ≅κk-λφlk-lkμlkσ2πlk3×exp⁡-κk-λφlk22σ2lk,lk≥0,where κk=κ-xk,  φ(lk)=(lk+tk)b-tkb, and μ(lk)=b(lk+tk)b-1.


### 3.4. Implementation with the PF Algorithm

In this work, the PF algorithm is utilized to integrate the vibration signals with the nonlinear degradation model described by (12) and (15). Due to the existing of uncertainties for the specified system under a certain working condition, the model parameters are not deterministic. Condition monitoring data like vibration signals contain specific information corresponding to the concerned system. By updating the parameters with on-line collected data, the uncertainties can be reduced. Therefore, the RUL can be forecasted more accurately with the updated parameters.The model parameterb is regarded as a constant for a simplified case in reported studies [27, 38, 39]. However, two parameters λ and b are statistically dependent random variables for the degradation process. For a more realistic demonstration, in this paper, λ andb are both updated with the available data. In practical, there is no certain information about the BM and the measurement noise for a specific system. Hence, parameters related to the BM and the measurement noise are also updated with λ and b, and we can get the following state space model:(19)xk=xk-1+λk-1tkbk-1-tk-1bk-1+ηkηk~N0,σk2Δtkλk=λk-1bk=bk-1σk=σk-1γk=γk-1yk=xk+νkνk~N0,γk2.Letzk=[xk,λk,bk,σk,γk]T denote the states of the system and execute the procedures of the standard PF explained in Section 2.2; we can get the estimated states z^k and the corresponding particles z^kii=1N at each inspection time.Based on the state estimation, we can achieve ap-step ahead state prediction by propagating the current distribution with the state space model (19) recursively. Specifically, it can be implemented with the following formula [40]: (20)p~zk+p∣y1:k=∫p~zk∣y1:k∏j=k+1k+ppzj∣zj-1dzk:k+p-1.However, it is both formidable and time-consuming to solve these integrals. By applying the PF algorithm, we can use the particlesz^kii=1N as the initial condition p~(zk∣y1:k) for prediction. If the particles z^kii=1N can well represent the posterior PDF at current time, the predicted PDF at time tk+β (β=1, …,p) can be approximated by (21).(21)pzk+β∣z1:k=∫pzk+β∣zk+β-1pzk+β-1∣z1:kdzk+β-1≈∫∑i=1Nwk+β-1iδzk+β-1-zk+β-1ipzk+β∣zk+β-1dzk+β-1=∑i=1Nwk+β-1ipzk+β∣zk+β-1i=p^zk+β∣y1:kTo calculate (21), the particles at each time instant are regarded as the prior state of the latter prediction at the next time instant. Specifically, the particles are recursively propagated by the state transition function:(22)z^k+βi=fz^k+β-1i.Since this approach can reduce computational time, it is suitable for on-line forecasting cases. Besides, we can derive the failure PDF for future time instants withp-step ahead prediction. The number of particles isN, which means that there areN predicted state transition paths. Consequently, the number of the HIs at time tk+β isN. If a HI exceeds the FT, the time tk+β is considered as the failure time for this prediction path, which means that we can calculate the failure PDF at any future time instant with the following:(23)p^ftk+β=∑i=1NPr⁡x^k+βi≥κ.Particularly, what we concern in this paper is to obtain the RUL distributionp(RUL∣y1:k). As described above,N degradation paths can be produced withN particles, and each path has a failure time. Hence, the number of RULs is also N. Let L^ki denote the RUL obtained by particle i, and for each degradation path, calculate its future states recursively by (22). When xk+βi>κ, let L^ki=β. Then we can get the PDF of RUL by L^kii=1N at the given time tk.The detailed procedures of the parameters update and RUL prediction process with the PF are illustrated in Table1.Table 1
Parameters update and RUL prediction process.
Step 1: Monitor the bearing’s health state using the method described in Section3.1. When the degradation starts, define the beginning time as t0. Step 2: At timet0, define the initial distribution of z0. Then sample particles z0ii=1Nfrom the initial distribution and initialize particle weights with w0i =1/N. Step 3: At any inspection timetk (k≥1), sample particles zkii=1N from the importance PDF defined by the state space model (19). Step 4: Update particle weightswkii=1N with the newly obtained HI using equation (8). Step 5: Resamplez^kdi=1N from zkii=1Naccording to wkii=1N. Step 6: Update the HI and the model parameters at timetk with equation (10). Step 7: Calculate the future states recursively by equation (22) for each particle, producing N degradation paths. For each degradation path, its failure time F^ki is defined as the time when the predicted HI exceeds the FT, and L^ki=F^ki-tk is the RUL of this degradation path. Consequently, the failure time PDF and the RUL PDF can be obtained from all the values of failure time F^kii=1N and RUL L^kii=1N, respectively. Step 8: For the next inspection interval, letk= k+ 1, then turn to Step 3 and repeat Step 3-7 until the estimated HI xk>κ.


## 3.1. Health Stage Division

To conduct the proposed method, an important task is to divide the health stages precisely. The statistical feature Root Mean Square (RMS) is used as the health indicator (HI) of bearings in this study due to its sensitivity to the global damage of bearings [44]. To detect the degradation, the HIs obtained from the healthy state are used to determine the DT. Specifically, the mean μn and the standard deviation σn of these HIs are calculated, respectively, and μn + 3σn is defined as the DT in this study. Then the health state of bearings is monitored by comparing the DT with the newly obtained HI.However, there are some abnormal values due to random noises, which could mislead the recognition of health states. To solve this problem, some measures are taken to find the start point of degradation, as shown in Figure2. The process of the degradation detection mechanism is described as follows:(1)
At inspection timetk, calculate the RMS of the new measured data.(2)
Compare the RMS attk (HIk) with the DT (μn + 3σn). If HIk > μn + 3σn, definetk astm; otherwise, setk =k + 1, return to Step 1.(3)
Compute the RMS attm+1 andtm+2, respectively. If these two values of RMS satisfy HIm+q-μn>3σnq=1,2, considertm as the beginning point of the degradation stage and conduct RUL prediction; otherwise, setk =k + 1, and return to Step 1 and repeat Steps 1-3.Figure 2
The degradation detection mechanism.Practically, the RMS of a bearing may exceed the DT because of random noises rather than the degradation. However, these abnormal scenarios are difficult to arise three times consecutively in the normal stage. Therefore, we define the bearing begins to degrade when three consecutive values of RMS are bigger than the DT.


## 3.2. Development of the Bearing Degradation Model

There are four major variability sources related to the uncertainties of machinery deterioration: the temporal variability, the measurement variability, the unit-to-unit variability, and the nonlinear variability [39]. To manage the uncertainties in a degradation process appropriately, it is necessary to consider all these factors. Accordingly, the degradation process of bearings can be illustrated by a nonlinear form as follows:(11)St=a+λtb+σBt,whereS(t) is the health state of a bearing at timet,a represents the initial health state of the degradation stage;λtb denotes the degradation process, λ reflects the degradation speed, andb represents the nonlinear characteristic;B(t) is the standard Brownian motion (BM), and σ is the diffusion coefficient.For convenience, (11) is transformed into(12)Xt=St-a=λtb+σBt.Accordingly, we can obtain the following difference equation of the degradation process:(13)xk=xk-1+λk-1tkbk-1-tk-1bk-1+ηk,where ηk=σ(B(tk)-B(tk-1)) follows a normal distribution N(0,σ2Δtk), and Δtk=tk-tk-1.In practice, vibration signals are the most commonly selected condition monitoring data, for the advantage that it can well indicate the degradation process of machinery. However, there exists a difference between the actual health state and the HI extracted from the vibration signal due to the noises from the sensors, acquisition systems, and data processing. Thus the measurement variability should be considered, and the observed health state can be described as follows [39]:(14)Zt=St+ν,whereZ(t) represents the HI extracted from the vibration signal at timet, and ν is the measurement noise, which follows a normal distribution N(0,γ2). Also, (14) can be derived as follows:(15)Yt=Zt-a=Xt+ν.In this paper, we useΘ = [λ,b, σ, γ] to denote the unknown model parameters.


## 3.3. RUL Prediction

The RULLk at time tk can be defined as follows:(16)Lk=inf⁡lk:Xtk+lk≥κ∣Y1:k,Θ,where inf⁡· stands for the inferior limit of a variable,lk is the time betweentk and the failure time, and κ is the predefined FT;Y1:k represents the available observations from the initial point of degradation stage to timetk, and the estimated health state at timetk+lk is denoted as follows:(17)Xtk+lk=xk+λtk+lkb-tkb+σBtk+lk-Btk.According to the dependent increment properties of the standard BM,B(tk+lk)-B(tk) is also a standard BM. Thus the predicted PDF of RUL at timetk can be illustrated as follows [27]:(18)fLk∣xk,Θlk∣xk,Θ≅κk-λφlk-lkμlkσ2πlk3×exp⁡-κk-λφlk22σ2lk,lk≥0,where κk=κ-xk,  φ(lk)=(lk+tk)b-tkb, and μ(lk)=b(lk+tk)b-1.


## 3.4. Implementation with the PF Algorithm

In this work, the PF algorithm is utilized to integrate the vibration signals with the nonlinear degradation model described by (12) and (15). Due to the existing of uncertainties for the specified system under a certain working condition, the model parameters are not deterministic. Condition monitoring data like vibration signals contain specific information corresponding to the concerned system. By updating the parameters with on-line collected data, the uncertainties can be reduced. Therefore, the RUL can be forecasted more accurately with the updated parameters.The model parameterb is regarded as a constant for a simplified case in reported studies [27, 38, 39]. However, two parameters λ and b are statistically dependent random variables for the degradation process. For a more realistic demonstration, in this paper, λ andb are both updated with the available data. In practical, there is no certain information about the BM and the measurement noise for a specific system. Hence, parameters related to the BM and the measurement noise are also updated with λ and b, and we can get the following state space model:(19)xk=xk-1+λk-1tkbk-1-tk-1bk-1+ηkηk~N0,σk2Δtkλk=λk-1bk=bk-1σk=σk-1γk=γk-1yk=xk+νkνk~N0,γk2.Letzk=[xk,λk,bk,σk,γk]T denote the states of the system and execute the procedures of the standard PF explained in Section 2.2; we can get the estimated states z^k and the corresponding particles z^kii=1N at each inspection time.Based on the state estimation, we can achieve ap-step ahead state prediction by propagating the current distribution with the state space model (19) recursively. Specifically, it can be implemented with the following formula [40]: (20)p~zk+p∣y1:k=∫p~zk∣y1:k∏j=k+1k+ppzj∣zj-1dzk:k+p-1.However, it is both formidable and time-consuming to solve these integrals. By applying the PF algorithm, we can use the particlesz^kii=1N as the initial condition p~(zk∣y1:k) for prediction. If the particles z^kii=1N can well represent the posterior PDF at current time, the predicted PDF at time tk+β (β=1, …,p) can be approximated by (21).(21)pzk+β∣z1:k=∫pzk+β∣zk+β-1pzk+β-1∣z1:kdzk+β-1≈∫∑i=1Nwk+β-1iδzk+β-1-zk+β-1ipzk+β∣zk+β-1dzk+β-1=∑i=1Nwk+β-1ipzk+β∣zk+β-1i=p^zk+β∣y1:kTo calculate (21), the particles at each time instant are regarded as the prior state of the latter prediction at the next time instant. Specifically, the particles are recursively propagated by the state transition function:(22)z^k+βi=fz^k+β-1i.Since this approach can reduce computational time, it is suitable for on-line forecasting cases. Besides, we can derive the failure PDF for future time instants withp-step ahead prediction. The number of particles isN, which means that there areN predicted state transition paths. Consequently, the number of the HIs at time tk+β isN. If a HI exceeds the FT, the time tk+β is considered as the failure time for this prediction path, which means that we can calculate the failure PDF at any future time instant with the following:(23)p^ftk+β=∑i=1NPr⁡x^k+βi≥κ.Particularly, what we concern in this paper is to obtain the RUL distributionp(RUL∣y1:k). As described above,N degradation paths can be produced withN particles, and each path has a failure time. Hence, the number of RULs is also N. Let L^ki denote the RUL obtained by particle i, and for each degradation path, calculate its future states recursively by (22). When xk+βi>κ, let L^ki=β. Then we can get the PDF of RUL by L^kii=1N at the given time tk.The detailed procedures of the parameters update and RUL prediction process with the PF are illustrated in Table1.Table 1
Parameters update and RUL prediction process.
Step 1: Monitor the bearing’s health state using the method described in Section3.1. When the degradation starts, define the beginning time as t0. Step 2: At timet0, define the initial distribution of z0. Then sample particles z0ii=1Nfrom the initial distribution and initialize particle weights with w0i =1/N. Step 3: At any inspection timetk (k≥1), sample particles zkii=1N from the importance PDF defined by the state space model (19). Step 4: Update particle weightswkii=1N with the newly obtained HI using equation (8). Step 5: Resamplez^kdi=1N from zkii=1Naccording to wkii=1N. Step 6: Update the HI and the model parameters at timetk with equation (10). Step 7: Calculate the future states recursively by equation (22) for each particle, producing N degradation paths. For each degradation path, its failure time F^ki is defined as the time when the predicted HI exceeds the FT, and L^ki=F^ki-tk is the RUL of this degradation path. Consequently, the failure time PDF and the RUL PDF can be obtained from all the values of failure time F^kii=1N and RUL L^kii=1N, respectively. Step 8: For the next inspection interval, letk= k+ 1, then turn to Step 3 and repeat Step 3-7 until the estimated HI xk>κ.


## 4. Experimental Studies

### 4.1. Simulation

In this section, a numerical example is utilized to verify the performance of the presented approach. Measurements of the normal stage and the degradation stage are both simulated. A nonlinear degradation path is generated according to the running process of bearings with (10), (12), and (13), which incorporate four uncertainty sources. The degradation path includes the measurement time and the corresponding HI extracted from sensor data. Specifically, the interval between two adjacent observations is 10 s, and the simulated HI is the RMS extracted from the vibration signal. The details of model parameters used in the simulation are shown in Table 2. In the simulation, xkk=1100 represent the hidden and actual health states of a bearing, while ykk=1100 are the corresponding measurements, as illustrated in Figure 3. The FT is set to be 55 g, and therefore the simulated process fails at 1000 s. The beginning point of the degradation stage is 510 s, as marked by the red line in Figure 3.Table 2
Parameters of the simulated degradation path.
Model parameters a λ b σ 2 γ 2 Value 5.00 0.02 2.00 0.25 4.00Figure 3
Simulated degradation path.In this case, the model parametersΘ=[λ,b,σ,γ] are updated together with the HIs {x51:51+h}(h=1, …, 49). In practice, it is difficult to conduct enough life tests to obtain a reliable prior PDF for a specific system. Hence, the prior distributions of the HI and the model parameters are given empirically, as shown in Table 3. The number of particles is 1000.Table 3
Prior distribution of the states.
Variable Distribution x 0 x 0 ~ U(0, 0.001) λ 0 λ 0 ~ U(0, 0.05) b 0 b 0 ~ U(1.5, 2.5) σ 0 σ 0 ~ U(0.1, 1) γ 0 γ 0 ~ U(1, 3)The estimations of the HI are shown in Figure3. It can be observed that the proposed method can track the actual degradation trend accurately. At each inspection time, we can calculate the mean of all the model parameters. Figure 4 illustrates the updating history of parameters during the estimation process. It can be seen that all the parameters fluctuate to some degree. The parameters λ,b,γ tend to approximate their true values after 820 s, while σ still varies slightly at this stage. The final estimated results of λ and b are 0.0172 and 2.055, respectively, which are a little deviated from the real values. The final estimated value of σ is 0.4406, and there is still a small gap between it and the actual value 0.5. The final estimation of γ is 2.0208, which is very close to the predefined value 2.Figure 4
Evolution of the parameters.At each inspection time, after we get the posterior distributions of the HI and all the related parameters, the RUL can be predicted by Step 7 in Table1. The mean values of the predicted RULs are shown in Figure 5. It can be observed that the predicted RULs deviate a lot from the real values at the beginning because of the multiple uncertainties. As more measurements become available, the RUL can be estimated more accurately due to the reduction of uncertainties. Specifically, accurate results can be obtained after 750 s.Figure 5
RUL prediction results.As described above, even the final estimated values of parametersλ,b, σ still differ a little from their real values, the predicted RUL converges to the real RUL as time goes. The reason for this phenomenon is that different combinations of parameters λ,b, and σ in the nonlinear degradation model can result in the same value of HI. Although there is a small difference between the estimations of parameters and their true values, the RUL can still be predicted accurately. A similar conclusion was also drawn in previous work [6]. In a word, this simulation indicates that the proposed method is effective in RUL prediction and uncertainty reduction even with random prior parameters.


### 4.2. Bearing RUL Prediction

In this section, real condition monitoring data collected from accelerated life tests for bearings are used to verify the proposed method’s effectiveness in real application. The datasets were used for the IEEE PHM 2012 Data Challenge [45]. The RMS is extracted from the vibration signal as the HI to indicate the health status of bearings, and then the presented approach is used to forecast the RUL of bearings.


#### 4.2.1. Experimental Setup and Data

The accelerated tests were conducted on an experimental platform called PRONOSTIA, which was designed and realized at FEMTO-ST Institute. The configuration of the platform is shown in Figure6. The experimental system can offer run-to-failure condition monitoring data for ball bearings. With the accelerated test techniques, bearings can fail within only a few hours, which saves a large amount of time and makes it possible to obtain enough run-to-failure data for research.Figure 6
Overview of the experimental platform.During the test, the radial force is set to be higher than the bearing’s maximum dynamic load. As a result, the time for a bearing to fail is reduced. The loading force and the rotating speed are constant during the overall life cycle. Two vibration sensors are installed on the vertical axis and the horizontal axis, respectively. The sample frequency of the vibration is 25.6 kHz, and the length of data is 2560. During the operational life, vibration signals are collected every 10 s.The force exerted horizontally is the main reason for the failure of bearings, so the horizontal signals capture more relevant information about the degradation status of bearings. Consequently, the RMS extracted from vibration signals in this direction are chosen as the HI of bearings in this study. Even in the same operation condition, the life durations of tested bearings vary a lot from each other. Accordingly, the uncertainty for each bearing should be quantified to obtain accurate RUL prediction results.


#### 4.2.2. RUL Estimation of Bearings

Figure7 shows the vibration signals evolution of a typical degradation process during its life time. The values of RMS during the bearing’s service life are illustrated in Figure 8. It can be observed that the bearing works normally for a long time. However, it is not necessary to predict the RUL when the bearing is healthy, and what we concern is the performance of the bearing during its degradation stage. Therefore, in the normal stage, we do not conduct the updating process, while in the degradation stage we incorporate the model developed in Section 3.2 and the vibration signals to track the health status and perform RUL prediction. Besides, prognostics based on the health stage division can improve the estimation accuracy and reduce the computational cost [46, 47]. Hence, the operational process of the bearing is divided into two stages according to the method described in Section 3.1. According to the steps of the degradation detection mechanism, the bearing begins to deteriorate att = 15850 s, as marked by the red line in Figure 8.Figure 7
Typical bearing vibration signals.Figure 8
RMS results of the bearing.When the bearing begins to deteriorate, the PF is utilized to estimate the health status and the parameters. The number of particles is 1000. Then the RUL can be forecasted with the updated results. The initial parameters are given randomly like what we did in Section4.1. The exponential Wiener process model proposed by Gebraeel [22] and the linear Wiener process model are utilized for comparison and demonstrate the performance of the presented framework. The number of particles for the Gebraeel’s model and the linear model are also 1000. The prior parameters for these two models are optimal which are obtained by fitting the RMS results. Specifically, following Steps 3-6 in Table 1, once a new measurement is obtained, we can update the health state and the parameters. The RMS estimations through the three models are shown in Figure 9. From the figure, we can see that the RMS tracking results of the exponential model and the linear model almost coincide with the real RMS, while the results of our method deviate a little from the actual RMS. Besides, the estimated RMS of the proposed approach has fewer fluctuations. The reason for this may be that our method is affected less by the random errors, which can indicate the intrinsic degradation evolution.Figure 9
RMS estimated results.At each inspection interval, we can obtain the RUL distribution by Step 7 in Table1 with the newly updated HI and parameters. Figure 10 illustrates the mean values of predicted RULs in the degradation stage with the three models. As more measurements become available, the prediction of our model converges to the real RUL after 15980 s. However, the predicted RUL by Gebraeel’s model still deviates to some extent until 16210 s. Among the three models, the linear model converges slowest, which do not get close to the actual RUL until the end. To evaluate the prediction accuracy of the three models, the root mean square error (RMSE) is calculated for each method. The RMSE of the Gebraeel’s model and the linear model are 217.79 and 115.05, respectively, which are both larger than 78.11 of our model. Consequently, our model can achieve a higher accuracy than the comparison models.Figure 10
RUL prediction results.To further verify the superiority of the presented method quantitatively, we use a commonly used metric, named Cumulative Relative Accuracy (CRA) [48] to evaluate the performance of the three models. The CRA is expressed as follows: (24)CRAζ=1Tζ∑k=1k=TζWLk∗RAk,where Tζ is the number of time indices, Lk∗ is the real RUL attk, andW(Lk∗) is a weight factor attk, which is a function of Lk∗. Here we defineW(Lk∗)=1 for all the time indexes. The RAk is the Relative Accuracy (RA) of RUL prediction attk, which is defined as follows:(25)RAk=1-Lk∗-LkLk∗,whereLk is the predicted RUL attk.We calculate the CRA values for each model in different intervals, and the results of CRA are shown in Table4. It is seen that the CRA of our model are much larger than that of the Gebraeel’s model for all the time intervals listed in Table 4. The linear model performs better than our model in the first two intervals, which indicates the severity of the degradation increases linearly at the beginning of the deterioration stage. However, our model can provide more accurate results in the latter three intervals, and our model can achieve much higher overall accuracy for the entire degradation stage than the linear model.Table 4
CRA of three models for different intervals.
T ζ Gebraeel’s model Linear model Our model 10 0.3638 0.9281 0.6365 20 0.2739 0.8183 0.7876 30 0.3141 0.6701 0.8443 40 0.4189 0.5064 0.8678 52 0.5034 0.4807 0.8688To sum up, among the three methods, our model converges fastest and provides the most accurate prediction results during the whole prognostics process.In order to demonstrate the effectiveness of the proposed model in uncertainty reduction, we calculate the PDFs of failure time for three approaches at different inspection time according to Step 7 in Table1, and the results are illustrated in Figures 11, 12, and 13. The distributions are all wide at the beginning. According to Table 1, the failure time PDF at a specific time instant is determined by the estimated health state and parameters. At the beginning of the updating process, there exist lots of uncertainties in the updated health state and parameters, which may cause an inaccurate failure time estimation. As the update times increase, the PDFs get closer to the actual failure time and become narrower. Since the distribution of failure time is nonnormal, we use Mean Absolute Deviation (MAD) [48] to measure the variance of distribution, which is shown in Table 5. It is observed that all the variances of our model at the listed updating indices are much smaller than that of the Gebraeel’s model. Besides, the variances of the linear model for most time indexes are bigger than that of our model except the 10th update. Therefore, our model performs best in the uncertainty management.Table 5
MAD of three methods at different time indexes.
T ζ Gebraeel’s model (s) The linear model (s) Our model (s) 10 4847.5296 13.1931 110.9541 20 1007.8243 16.0942 4.6667 30 217.8650 14.4032 5.2765 40 33.8555 12.2644 4.6832 50 27.2808 12.3839 3.7097Figure 11
Updated failure time distribution of our model.Figure 12
Updated failure time distribution of Gebraeel’s model.Figure 13
Updated failure time distribution of the linear model.To further demonstrate the effectiveness of the proposed approach, we calculate the RUL of other bearing datasets with the proposed method. Due to the lack of space, we do not show the pictures of their results and just give the above-mentioned bearing for example. Table6 displays the CRA values during the whole degradation stage for the tested bearings. It is observed that our model performs better than other two models for almost all the tested bearings, which further proves its effectiveness in real application.Table 6
CRA for tested bearings.
No. Gebraeel’s model Linear model Our model Bearing 1 0.6230 0.5890 0.6960 Bearing 2 0.6096 0.7661 0.7530 Bearing 3 0.5411 0.7771 0.8429 Bearing 4 0.6961 0.6402 0.7808 Bearing 5 0.7487 0.7579 0.8324 Bearing 6 0.6876 0.7573 0.7647 Bearing 7 0.3101 0.6159 0.7603The superior properties of the presented method can be explained as follows:(1)
The operational process of the bearing is divided into two stages, which are processed separately. Modeling and prognostics with the information of degradation stage alone can avoid the inference caused by the measurements from the normal stage.(2)
Compared with our model, the exponential model and the linear model estimated the actual health state less accurately because of the random errors. Consequently, there is a relatively large gap between the actual RUL and their estimated results. However, the random errors have less influence on the health state estimation of the presented approach. Thus, the corresponding predicted RUL accords better with the real RUL.(3)
In our model, both BMσ(B(tk)-B(tk-1)) and the measurement noise ν are considered, and the corresponding parameters are updated by Bayesian inference. As a result, it can indicate the temporal variability and the measurement variability in real time. Moreover, the exponential model can be transformed into a linear model. In our model, the nonlinear variability and the unit-to-unit variability are both updated with the condition monitoring data. In a word, the proposed method provides the best uncertainty management among the three approaches. Therefore, our model can achieve higher prognostics accuracy than the two traditional methods.


## 4.1. Simulation

In this section, a numerical example is utilized to verify the performance of the presented approach. Measurements of the normal stage and the degradation stage are both simulated. A nonlinear degradation path is generated according to the running process of bearings with (10), (12), and (13), which incorporate four uncertainty sources. The degradation path includes the measurement time and the corresponding HI extracted from sensor data. Specifically, the interval between two adjacent observations is 10 s, and the simulated HI is the RMS extracted from the vibration signal. The details of model parameters used in the simulation are shown in Table 2. In the simulation, xkk=1100 represent the hidden and actual health states of a bearing, while ykk=1100 are the corresponding measurements, as illustrated in Figure 3. The FT is set to be 55 g, and therefore the simulated process fails at 1000 s. The beginning point of the degradation stage is 510 s, as marked by the red line in Figure 3.Table 2
Parameters of the simulated degradation path.
Model parameters a λ b σ 2 γ 2 Value 5.00 0.02 2.00 0.25 4.00Figure 3
Simulated degradation path.In this case, the model parametersΘ=[λ,b,σ,γ] are updated together with the HIs {x51:51+h}(h=1, …, 49). In practice, it is difficult to conduct enough life tests to obtain a reliable prior PDF for a specific system. Hence, the prior distributions of the HI and the model parameters are given empirically, as shown in Table 3. The number of particles is 1000.Table 3
Prior distribution of the states.
Variable Distribution x 0 x 0 ~ U(0, 0.001) λ 0 λ 0 ~ U(0, 0.05) b 0 b 0 ~ U(1.5, 2.5) σ 0 σ 0 ~ U(0.1, 1) γ 0 γ 0 ~ U(1, 3)The estimations of the HI are shown in Figure3. It can be observed that the proposed method can track the actual degradation trend accurately. At each inspection time, we can calculate the mean of all the model parameters. Figure 4 illustrates the updating history of parameters during the estimation process. It can be seen that all the parameters fluctuate to some degree. The parameters λ,b,γ tend to approximate their true values after 820 s, while σ still varies slightly at this stage. The final estimated results of λ and b are 0.0172 and 2.055, respectively, which are a little deviated from the real values. The final estimated value of σ is 0.4406, and there is still a small gap between it and the actual value 0.5. The final estimation of γ is 2.0208, which is very close to the predefined value 2.Figure 4
Evolution of the parameters.At each inspection time, after we get the posterior distributions of the HI and all the related parameters, the RUL can be predicted by Step 7 in Table1. The mean values of the predicted RULs are shown in Figure 5. It can be observed that the predicted RULs deviate a lot from the real values at the beginning because of the multiple uncertainties. As more measurements become available, the RUL can be estimated more accurately due to the reduction of uncertainties. Specifically, accurate results can be obtained after 750 s.Figure 5
RUL prediction results.As described above, even the final estimated values of parametersλ,b, σ still differ a little from their real values, the predicted RUL converges to the real RUL as time goes. The reason for this phenomenon is that different combinations of parameters λ,b, and σ in the nonlinear degradation model can result in the same value of HI. Although there is a small difference between the estimations of parameters and their true values, the RUL can still be predicted accurately. A similar conclusion was also drawn in previous work [6]. In a word, this simulation indicates that the proposed method is effective in RUL prediction and uncertainty reduction even with random prior parameters.


## 4.2. Bearing RUL Prediction

In this section, real condition monitoring data collected from accelerated life tests for bearings are used to verify the proposed method’s effectiveness in real application. The datasets were used for the IEEE PHM 2012 Data Challenge [45]. The RMS is extracted from the vibration signal as the HI to indicate the health status of bearings, and then the presented approach is used to forecast the RUL of bearings.


### 4.2.1. Experimental Setup and Data

The accelerated tests were conducted on an experimental platform called PRONOSTIA, which was designed and realized at FEMTO-ST Institute. The configuration of the platform is shown in Figure6. The experimental system can offer run-to-failure condition monitoring data for ball bearings. With the accelerated test techniques, bearings can fail within only a few hours, which saves a large amount of time and makes it possible to obtain enough run-to-failure data for research.Figure 6
Overview of the experimental platform.During the test, the radial force is set to be higher than the bearing’s maximum dynamic load. As a result, the time for a bearing to fail is reduced. The loading force and the rotating speed are constant during the overall life cycle. Two vibration sensors are installed on the vertical axis and the horizontal axis, respectively. The sample frequency of the vibration is 25.6 kHz, and the length of data is 2560. During the operational life, vibration signals are collected every 10 s.The force exerted horizontally is the main reason for the failure of bearings, so the horizontal signals capture more relevant information about the degradation status of bearings. Consequently, the RMS extracted from vibration signals in this direction are chosen as the HI of bearings in this study. Even in the same operation condition, the life durations of tested bearings vary a lot from each other. Accordingly, the uncertainty for each bearing should be quantified to obtain accurate RUL prediction results.


### 4.2.2. RUL Estimation of Bearings

Figure7 shows the vibration signals evolution of a typical degradation process during its life time. The values of RMS during the bearing’s service life are illustrated in Figure 8. It can be observed that the bearing works normally for a long time. However, it is not necessary to predict the RUL when the bearing is healthy, and what we concern is the performance of the bearing during its degradation stage. Therefore, in the normal stage, we do not conduct the updating process, while in the degradation stage we incorporate the model developed in Section 3.2 and the vibration signals to track the health status and perform RUL prediction. Besides, prognostics based on the health stage division can improve the estimation accuracy and reduce the computational cost [46, 47]. Hence, the operational process of the bearing is divided into two stages according to the method described in Section 3.1. According to the steps of the degradation detection mechanism, the bearing begins to deteriorate att = 15850 s, as marked by the red line in Figure 8.Figure 7
Typical bearing vibration signals.Figure 8
RMS results of the bearing.When the bearing begins to deteriorate, the PF is utilized to estimate the health status and the parameters. The number of particles is 1000. Then the RUL can be forecasted with the updated results. The initial parameters are given randomly like what we did in Section4.1. The exponential Wiener process model proposed by Gebraeel [22] and the linear Wiener process model are utilized for comparison and demonstrate the performance of the presented framework. The number of particles for the Gebraeel’s model and the linear model are also 1000. The prior parameters for these two models are optimal which are obtained by fitting the RMS results. Specifically, following Steps 3-6 in Table 1, once a new measurement is obtained, we can update the health state and the parameters. The RMS estimations through the three models are shown in Figure 9. From the figure, we can see that the RMS tracking results of the exponential model and the linear model almost coincide with the real RMS, while the results of our method deviate a little from the actual RMS. Besides, the estimated RMS of the proposed approach has fewer fluctuations. The reason for this may be that our method is affected less by the random errors, which can indicate the intrinsic degradation evolution.Figure 9
RMS estimated results.At each inspection interval, we can obtain the RUL distribution by Step 7 in Table1 with the newly updated HI and parameters. Figure 10 illustrates the mean values of predicted RULs in the degradation stage with the three models. As more measurements become available, the prediction of our model converges to the real RUL after 15980 s. However, the predicted RUL by Gebraeel’s model still deviates to some extent until 16210 s. Among the three models, the linear model converges slowest, which do not get close to the actual RUL until the end. To evaluate the prediction accuracy of the three models, the root mean square error (RMSE) is calculated for each method. The RMSE of the Gebraeel’s model and the linear model are 217.79 and 115.05, respectively, which are both larger than 78.11 of our model. Consequently, our model can achieve a higher accuracy than the comparison models.Figure 10
RUL prediction results.To further verify the superiority of the presented method quantitatively, we use a commonly used metric, named Cumulative Relative Accuracy (CRA) [48] to evaluate the performance of the three models. The CRA is expressed as follows: (24)CRAζ=1Tζ∑k=1k=TζWLk∗RAk,where Tζ is the number of time indices, Lk∗ is the real RUL attk, andW(Lk∗) is a weight factor attk, which is a function of Lk∗. Here we defineW(Lk∗)=1 for all the time indexes. The RAk is the Relative Accuracy (RA) of RUL prediction attk, which is defined as follows:(25)RAk=1-Lk∗-LkLk∗,whereLk is the predicted RUL attk.We calculate the CRA values for each model in different intervals, and the results of CRA are shown in Table4. It is seen that the CRA of our model are much larger than that of the Gebraeel’s model for all the time intervals listed in Table 4. The linear model performs better than our model in the first two intervals, which indicates the severity of the degradation increases linearly at the beginning of the deterioration stage. However, our model can provide more accurate results in the latter three intervals, and our model can achieve much higher overall accuracy for the entire degradation stage than the linear model.Table 4
CRA of three models for different intervals.
T ζ Gebraeel’s model Linear model Our model 10 0.3638 0.9281 0.6365 20 0.2739 0.8183 0.7876 30 0.3141 0.6701 0.8443 40 0.4189 0.5064 0.8678 52 0.5034 0.4807 0.8688To sum up, among the three methods, our model converges fastest and provides the most accurate prediction results during the whole prognostics process.In order to demonstrate the effectiveness of the proposed model in uncertainty reduction, we calculate the PDFs of failure time for three approaches at different inspection time according to Step 7 in Table1, and the results are illustrated in Figures 11, 12, and 13. The distributions are all wide at the beginning. According to Table 1, the failure time PDF at a specific time instant is determined by the estimated health state and parameters. At the beginning of the updating process, there exist lots of uncertainties in the updated health state and parameters, which may cause an inaccurate failure time estimation. As the update times increase, the PDFs get closer to the actual failure time and become narrower. Since the distribution of failure time is nonnormal, we use Mean Absolute Deviation (MAD) [48] to measure the variance of distribution, which is shown in Table 5. It is observed that all the variances of our model at the listed updating indices are much smaller than that of the Gebraeel’s model. Besides, the variances of the linear model for most time indexes are bigger than that of our model except the 10th update. Therefore, our model performs best in the uncertainty management.Table 5
MAD of three methods at different time indexes.
T ζ Gebraeel’s model (s) The linear model (s) Our model (s) 10 4847.5296 13.1931 110.9541 20 1007.8243 16.0942 4.6667 30 217.8650 14.4032 5.2765 40 33.8555 12.2644 4.6832 50 27.2808 12.3839 3.7097Figure 11
Updated failure time distribution of our model.Figure 12
Updated failure time distribution of Gebraeel’s model.Figure 13
Updated failure time distribution of the linear model.To further demonstrate the effectiveness of the proposed approach, we calculate the RUL of other bearing datasets with the proposed method. Due to the lack of space, we do not show the pictures of their results and just give the above-mentioned bearing for example. Table6 displays the CRA values during the whole degradation stage for the tested bearings. It is observed that our model performs better than other two models for almost all the tested bearings, which further proves its effectiveness in real application.Table 6
CRA for tested bearings.
No. Gebraeel’s model Linear model Our model Bearing 1 0.6230 0.5890 0.6960 Bearing 2 0.6096 0.7661 0.7530 Bearing 3 0.5411 0.7771 0.8429 Bearing 4 0.6961 0.6402 0.7808 Bearing 5 0.7487 0.7579 0.8324 Bearing 6 0.6876 0.7573 0.7647 Bearing 7 0.3101 0.6159 0.7603The superior properties of the presented method can be explained as follows:(1)
The operational process of the bearing is divided into two stages, which are processed separately. Modeling and prognostics with the information of degradation stage alone can avoid the inference caused by the measurements from the normal stage.(2)
Compared with our model, the exponential model and the linear model estimated the actual health state less accurately because of the random errors. Consequently, there is a relatively large gap between the actual RUL and their estimated results. However, the random errors have less influence on the health state estimation of the presented approach. Thus, the corresponding predicted RUL accords better with the real RUL.(3)
In our model, both BMσ(B(tk)-B(tk-1)) and the measurement noise ν are considered, and the corresponding parameters are updated by Bayesian inference. As a result, it can indicate the temporal variability and the measurement variability in real time. Moreover, the exponential model can be transformed into a linear model. In our model, the nonlinear variability and the unit-to-unit variability are both updated with the condition monitoring data. In a word, the proposed method provides the best uncertainty management among the three approaches. Therefore, our model can achieve higher prognostics accuracy than the two traditional methods.


## 4.2.1. Experimental Setup and Data

The accelerated tests were conducted on an experimental platform called PRONOSTIA, which was designed and realized at FEMTO-ST Institute. The configuration of the platform is shown in Figure6. The experimental system can offer run-to-failure condition monitoring data for ball bearings. With the accelerated test techniques, bearings can fail within only a few hours, which saves a large amount of time and makes it possible to obtain enough run-to-failure data for research.Figure 6
Overview of the experimental platform.During the test, the radial force is set to be higher than the bearing’s maximum dynamic load. As a result, the time for a bearing to fail is reduced. The loading force and the rotating speed are constant during the overall life cycle. Two vibration sensors are installed on the vertical axis and the horizontal axis, respectively. The sample frequency of the vibration is 25.6 kHz, and the length of data is 2560. During the operational life, vibration signals are collected every 10 s.The force exerted horizontally is the main reason for the failure of bearings, so the horizontal signals capture more relevant information about the degradation status of bearings. Consequently, the RMS extracted from vibration signals in this direction are chosen as the HI of bearings in this study. Even in the same operation condition, the life durations of tested bearings vary a lot from each other. Accordingly, the uncertainty for each bearing should be quantified to obtain accurate RUL prediction results.


## 4.2.2. RUL Estimation of Bearings

Figure7 shows the vibration signals evolution of a typical degradation process during its life time. The values of RMS during the bearing’s service life are illustrated in Figure 8. It can be observed that the bearing works normally for a long time. However, it is not necessary to predict the RUL when the bearing is healthy, and what we concern is the performance of the bearing during its degradation stage. Therefore, in the normal stage, we do not conduct the updating process, while in the degradation stage we incorporate the model developed in Section 3.2 and the vibration signals to track the health status and perform RUL prediction. Besides, prognostics based on the health stage division can improve the estimation accuracy and reduce the computational cost [46, 47]. Hence, the operational process of the bearing is divided into two stages according to the method described in Section 3.1. According to the steps of the degradation detection mechanism, the bearing begins to deteriorate att = 15850 s, as marked by the red line in Figure 8.Figure 7
Typical bearing vibration signals.Figure 8
RMS results of the bearing.When the bearing begins to deteriorate, the PF is utilized to estimate the health status and the parameters. The number of particles is 1000. Then the RUL can be forecasted with the updated results. The initial parameters are given randomly like what we did in Section4.1. The exponential Wiener process model proposed by Gebraeel [22] and the linear Wiener process model are utilized for comparison and demonstrate the performance of the presented framework. The number of particles for the Gebraeel’s model and the linear model are also 1000. The prior parameters for these two models are optimal which are obtained by fitting the RMS results. Specifically, following Steps 3-6 in Table 1, once a new measurement is obtained, we can update the health state and the parameters. The RMS estimations through the three models are shown in Figure 9. From the figure, we can see that the RMS tracking results of the exponential model and the linear model almost coincide with the real RMS, while the results of our method deviate a little from the actual RMS. Besides, the estimated RMS of the proposed approach has fewer fluctuations. The reason for this may be that our method is affected less by the random errors, which can indicate the intrinsic degradation evolution.Figure 9
RMS estimated results.At each inspection interval, we can obtain the RUL distribution by Step 7 in Table1 with the newly updated HI and parameters. Figure 10 illustrates the mean values of predicted RULs in the degradation stage with the three models. As more measurements become available, the prediction of our model converges to the real RUL after 15980 s. However, the predicted RUL by Gebraeel’s model still deviates to some extent until 16210 s. Among the three models, the linear model converges slowest, which do not get close to the actual RUL until the end. To evaluate the prediction accuracy of the three models, the root mean square error (RMSE) is calculated for each method. The RMSE of the Gebraeel’s model and the linear model are 217.79 and 115.05, respectively, which are both larger than 78.11 of our model. Consequently, our model can achieve a higher accuracy than the comparison models.Figure 10
RUL prediction results.To further verify the superiority of the presented method quantitatively, we use a commonly used metric, named Cumulative Relative Accuracy (CRA) [48] to evaluate the performance of the three models. The CRA is expressed as follows: (24)CRAζ=1Tζ∑k=1k=TζWLk∗RAk,where Tζ is the number of time indices, Lk∗ is the real RUL attk, andW(Lk∗) is a weight factor attk, which is a function of Lk∗. Here we defineW(Lk∗)=1 for all the time indexes. The RAk is the Relative Accuracy (RA) of RUL prediction attk, which is defined as follows:(25)RAk=1-Lk∗-LkLk∗,whereLk is the predicted RUL attk.We calculate the CRA values for each model in different intervals, and the results of CRA are shown in Table4. It is seen that the CRA of our model are much larger than that of the Gebraeel’s model for all the time intervals listed in Table 4. The linear model performs better than our model in the first two intervals, which indicates the severity of the degradation increases linearly at the beginning of the deterioration stage. However, our model can provide more accurate results in the latter three intervals, and our model can achieve much higher overall accuracy for the entire degradation stage than the linear model.Table 4
CRA of three models for different intervals.
T ζ Gebraeel’s model Linear model Our model 10 0.3638 0.9281 0.6365 20 0.2739 0.8183 0.7876 30 0.3141 0.6701 0.8443 40 0.4189 0.5064 0.8678 52 0.5034 0.4807 0.8688To sum up, among the three methods, our model converges fastest and provides the most accurate prediction results during the whole prognostics process.In order to demonstrate the effectiveness of the proposed model in uncertainty reduction, we calculate the PDFs of failure time for three approaches at different inspection time according to Step 7 in Table1, and the results are illustrated in Figures 11, 12, and 13. The distributions are all wide at the beginning. According to Table 1, the failure time PDF at a specific time instant is determined by the estimated health state and parameters. At the beginning of the updating process, there exist lots of uncertainties in the updated health state and parameters, which may cause an inaccurate failure time estimation. As the update times increase, the PDFs get closer to the actual failure time and become narrower. Since the distribution of failure time is nonnormal, we use Mean Absolute Deviation (MAD) [48] to measure the variance of distribution, which is shown in Table 5. It is observed that all the variances of our model at the listed updating indices are much smaller than that of the Gebraeel’s model. Besides, the variances of the linear model for most time indexes are bigger than that of our model except the 10th update. Therefore, our model performs best in the uncertainty management.Table 5
MAD of three methods at different time indexes.
T ζ Gebraeel’s model (s) The linear model (s) Our model (s) 10 4847.5296 13.1931 110.9541 20 1007.8243 16.0942 4.6667 30 217.8650 14.4032 5.2765 40 33.8555 12.2644 4.6832 50 27.2808 12.3839 3.7097Figure 11
Updated failure time distribution of our model.Figure 12
Updated failure time distribution of Gebraeel’s model.Figure 13
Updated failure time distribution of the linear model.To further demonstrate the effectiveness of the proposed approach, we calculate the RUL of other bearing datasets with the proposed method. Due to the lack of space, we do not show the pictures of their results and just give the above-mentioned bearing for example. Table6 displays the CRA values during the whole degradation stage for the tested bearings. It is observed that our model performs better than other two models for almost all the tested bearings, which further proves its effectiveness in real application.Table 6
CRA for tested bearings.
No. Gebraeel’s model Linear model Our model Bearing 1 0.6230 0.5890 0.6960 Bearing 2 0.6096 0.7661 0.7530 Bearing 3 0.5411 0.7771 0.8429 Bearing 4 0.6961 0.6402 0.7808 Bearing 5 0.7487 0.7579 0.8324 Bearing 6 0.6876 0.7573 0.7647 Bearing 7 0.3101 0.6159 0.7603The superior properties of the presented method can be explained as follows:(1)
The operational process of the bearing is divided into two stages, which are processed separately. Modeling and prognostics with the information of degradation stage alone can avoid the inference caused by the measurements from the normal stage.(2)
Compared with our model, the exponential model and the linear model estimated the actual health state less accurately because of the random errors. Consequently, there is a relatively large gap between the actual RUL and their estimated results. However, the random errors have less influence on the health state estimation of the presented approach. Thus, the corresponding predicted RUL accords better with the real RUL.(3)
In our model, both BMσ(B(tk)-B(tk-1)) and the measurement noise ν are considered, and the corresponding parameters are updated by Bayesian inference. As a result, it can indicate the temporal variability and the measurement variability in real time. Moreover, the exponential model can be transformed into a linear model. In our model, the nonlinear variability and the unit-to-unit variability are both updated with the condition monitoring data. In a word, the proposed method provides the best uncertainty management among the three approaches. Therefore, our model can achieve higher prognostics accuracy than the two traditional methods.


## 5. Conclusions

This paper presents a new RUL prediction method for bearings based on a Wiener process model. First, the condition monitoring data collected in the normal stage is used to define the DT. Then we use the DT to detect whether the bearing begins to deteriorate. When the degradation occurs, the PF is used to integrate the degradation model with measurements for parameters updating and RUL prediction. To manifest the validity of the presented framework, the Gebraeel’s model and the linear model are used for comparison. The results show that our model performs best in RUL prediction and uncertainty management. Additionally, our model can deal with the cases whose initial parameters deviate a lot from their real values, which makes our model applicable in the situation where there is no enough failure history data for parameter initialization.



---

*Source: 4068431-2018-06-26.xml*