Remaining useful life prediction of rolling bearings by the particle filter method based on degradation rate tracking

2.1. Optimal Bayesian filtering

Supposing that the dynamic system can be described by a mathematical model, the state equation and measurement equation are expressed as follows:

where $f(\bullet )$: ${\mathbf{R}}^{n_x}\times {\mathbf{R}}^{n_n}\to {\mathbf{R}}^{n_x}$ is the nonlinear function of the system state ${\mathbf{x}}_{k-1}$; $\{{\mathbf{n}}_{k-1},k\in \mathbf{N}\}$ is stationary noise; $n_x$ and $n_n$ are the dimensions of the state and process noise, $h(\bullet )$: ${\mathbf{R}}^{n_x}\times {\mathbf{R}}^{n_v}\to {\mathbf{R}}^{n_y}$ is the nonlinear function of the system state ${\mathbf{x}}_k$; $\{{\mathbf{\omega }}_k,k\in \mathbf{N}\}$ is stationary noise; $n_y$ and $n_v$ are respectively the measurement dimension and observation noise dimension.

If the initial PDF of the state is known as $p\left({\mathbf{x}}_0|{\mathbf{y}}_0\right)=p\left({\mathbf{x}}_0\right)$, then the state prediction equation is:

The state updating equation is:

where:

2.2. Sampling

With Monte Carlo method, the integral calculation can be converted into a summation operation of limited samples. If $q\left({\mathbf{x}}_{0:k}|{\mathbf{y}}_{1:k}\right)$ is the first-order Markov process with easy sampling, its probability density distribution is known and the same to $p\left({\mathbf{x}}_{0:k}|{\mathbf{y}}_{1:k}\right)$. Supposing that we can extract $N$ samples $\{{\mathbf{x}}_{0:k}^i,\mathrm{}i=\text{1, 2}\text{,…,}N\}$${(\mathbf{x}}_{0:k}=\left\{\left({\mathbf{x}}_0,\dots ,{\mathbf{x}}_k\right\}\right)$ from the probability density distribution $q\left({\mathbf{x}}_{0:k}|{\mathbf{y}}_{1:k}\right)$, then the posterior probability density distribution of the state is obtained approximately by Eq. (5):

where $w_k\left({\mathbf{x}}_{0:k}^i\right)$ is the weight of ${\mathbf{x}}_{0:k}^i$ ($i=\text{1, 2,…,}N\text{)}$; $\delta (\bullet )$ is the Dirac function; ${\tilde{w}}_k\left({\mathbf{x}}_{0:k}^i\right)$ is the weight of normalized samples:

If $q\left({\mathbf{x}}_k|{\mathbf{x}}_{0:k-1},{\mathbf{y}}_{1:k}\right)=q\left({\mathbf{x}}_k|{\mathbf{x}}_{k-1},{\mathbf{y}}_k\right)$, then $w_k$ just depends on ${\mathbf{x}}_{k-1}$ and ${\mathbf{y}}_k$:

In the standard filtering algorithm, take:

Then, the weight $w_k^i$, normalized weight ${\tilde{w}}_k^i$, and posterior probability density $p\left({\mathbf{x}}_k|{\mathbf{y}}_{1:k}\right)$ can be respectively expressed as follows:

2.3. Resampling

One of the main difficulties in the development of PF is the degeneracy problem. As the algorithm evolves in time, the weight variances increase and the importance weight may be concentrated on a small number of particles. Therefore, the predicted value becomes progressively more skewed from the real value. In order to solve the degeneracy problem and adapt the dynamic process of the system, a resampling algorithm is performed to eliminate particles with small weights and duplicate those with large ones. The number of effective particles is defined as:

Whenever $N_{eff}<N_{th}$, a resampling algorithm is performed. Smaller $N_{eff}$ means more serious degeneracy. $N_{th}$ is a fixed threshold and is generally set to be $\text{2}N/\text{3}$. $N$ is the number of particles. The stratified resampling algorithm proposed by Carpenter et al. can improve the accuracy of PF algorithm and cut down computational time [17]. If considering residual information, the stratified resampling algorithm will be more reasonable. Therefore, in this paper, stratified resampling algorithm considering residual information is used for PF.

2.4. The implementation steps of the prediction algorithm based on particle filter

In the prediction strategy of system state with PF, the prediction model is determined firstly, and then model parameters are filtered based on known observations. Through iterative updating, model parameters are approximated to the real value as possible. Finally, the estimated values of model parameters are substituted into the prediction model to extrapolate and predict the trend of system degradation state. The steps are described as follows:

Step 1. Extracting degradation features. According to the raw measurement data, the proper degradation features are chosen and extracted for predictive modeling.

Step 2. Obtaining the priori model knowledge. With the extracted degradation features, based on historical data or experiences, the prediction model as well as the initial distribution of the model parameters is determined.

Step 3. Initialization of particle samples. The particle population is randomly generated with the priori probability $p\left({\mathbf{x}}_0\right)$, and the weight of every particle is $\text{1}/N$.

Step 4. Recursive prediction. The state at time $k+\text{1}$ is predicted with the state equation.

Step 5. Updating the particles’ state. Given the current observation, the likelihood function of particles and the observation value are calculated. Accordingly, the particle weights are updated and normalized to obtain the least square estimate of state $\mathbf{x}$ at time $k$ as ${\tilde{\mathbf{x}}}_k\approx \sum _{i=1}^N{w}_k^i{\mathbf{x}}_k^i$.

Step 6. Resampling. Resampling is performed according to the particle weights to obtain a new set of particles.

Step 7. At the time $k+\text{1}$, if there is a new observation, go back to Step 2, otherwise go to Step 8.

Step 8. The statistical results of the prediction. All particles are extrapolated for prediction with the prediction model and estimations of model parameters until the specified time or a preset threshold is reached. The statistics of the predicted value ${\mathbf{x}}_i$ or the time when the specified threshold of $N$ particles is reached were performed to obtain PDF of $\mathbf{x}$ at specified time or the PDF of time to exceed a threshold.

3.1. Shortcomings of the basic particle filter prediction method

Although the basic PF prediction method has been applied in different domains, in the process of its implementation, there are still some problems.

1) The “quality” of degradation features directly affects the accuracy of the prediction results. The “quality” refers to the correlation degree of degradation feature to the real degradation trend of the target system. In most cases, features directly measured or extracted would be affected by a variety of factors, such as environmental noise. These features can only reflect the approximate degradation trend of the system. If these features are used for prediction in particle filtering, observation noise would disturb the iterative updating of the model, and the model parameters are prone to diverge.

2) It is not easy to determine the prediction model. In the basic PF prediction method, the selection of the prediction model is an essential and important process and the accuracy of the prediction results depends heavily on the prediction model. Even subtle change of the prediction model may lead to the significant differences in predicted results and it is very disadvantageous to the robustness of the prediction method [15].

3) With more parameters, the prediction model makes particle filtering of multi-dimensional state more difficult. Usually, the selection of the prediction model depends on the fitting error of the model for historical data. Simple models are often not able to accurately indicate the degradation trend of historical data. However, an accurate prediction model is relatively complex with more model parameters, which increase the implementation difficulty of the method and the computation load of the algorithm.

4) Prediction results are much sensitive to the initial distribution of model parameters. Because the parameters of usual prediction model have no actual observation in the updating phase, the initial distribution set of these parameters directly determine whether the parameters can be converged, and control the convergence rate. Although a suitable initial distribution of model parameters was proposed according to the Dempster-Shafer theory in ref. [13], in fact this approach does not always give the most proper results.

3.2. Improvements

To overcome the above disadvantages, we proposed a DRT-PF method, including three improvements based on the basic PF prediction framework as follows.

1) Moving average smoothing and monotonizing for the degradation features sequentially. Because we only care about the degradation trend of the system reflected by degradation features while excluding other details contained in the features, it is necessary to remove the baneful influence of noise or other external factors as possible. Through the above pretreatment, we can obtain the features which indicate the degradation trend of the system more intuitively. Fig. 1 illustrates the comparison between raw degradation features and the pretreated features.

Fig. 1Comparison between raw degradation features and the pretreated features

2) The following general prediction model is used to replace the previous specific model, which needs to be identified, to simplify the prediction process:

where ${\mathbf{a}}_k$ is the degradation rate of the feature and defined as the ratio of the state at time $k+\text{1}$ to the state at the previous time $k$. In general prediction model, the content of every particle is simplified and a particle includes only the state $\mathbf{x}$ and the degradation rate $\mathbf{a}$, thus reducing computational complexity considerably. Moreover, we can compute real-time observations of $\mathbf{a}$ in the iterative updating phase. It means that the resampling can be used in the estimation of $\mathbf{a}$, significantly reduce the influence of the initial distribution set of model parameters on prediction results.

3) The priori information of the degradation rate $\mathbf{a}$ should be utilized completely. The degradation rate $\mathbf{a}$ reflects the system degradation rate from the normal state to failure, and the estimation of $\mathbf{a}$ determines the prediction of RUL for the system. Therefore, the accurate estimation for the evolution trend of degradation rate in the prediction phase is very important. Because the degradation behavior of equipment is a one-way variation over time, the range of the degradation rate after smoothing and monotonizing pretreatment is easily determined. According to the signification of degradation feature, a boundary $\mathbf{a}\ge \text{1}$ or $\mathbf{a}\le \text{1}$ can be determined firstly. Another boundary can be obtained from the statistics of the extreme value of $\mathbf{a}$ in historical data population. The distribution map of the degradation rate with different features is computed based on historical data statistics to guide the updating of particle state in the prediction phase.

The flowchart of the DRT-PF prediction method is shown in Fig. 2.

Fig. 2Framework of the DRT-PF prediction algorithm

4.1. Extraction of degradation features

From the experimental data, we tried to extract different features in the time and frequency domains. However, neither the time domain features, such as root mean square (RMS), crest factor, and kurtosis of vibration signals, nor the frequency domain features, such as the magnitude at bearing defect frequencies, exhibited a consistent trend of bearing degradation. Therefore, according to the rating standard of vibration quality ISO 2372-1974, we selected the vibration intensity which combines the advantages of both time and frequency domain features as the degradation feature for RUL prediction [19]. The vibration intensity is defined as the RMS value of vibration velocity in the frequency range of 10-1000 Hz. It is a comprehensive description of equipment vibration and reflects the quantity of total vibrational energy including the energy of each harmonic. The vibration velocitycan be expressed as:

where $T$ is time length of the signal, $\mathbf{v}\left(t\right)$ is the vibration velocity of the object system. If the signal is discreet, Eq. (14) can be written as:

According to the properties of the Fourier transform, vibration intensity can be calculated in the frequency domain. For the vibration signal $\mathbf{x}\left(n\right)$ of $N$ points, the sampling frequency is $f_s$, the vibration intensity can be calculated with DFT as:

The single-sided amplitude spectrum of vibration signal can be expressed as:

The corresponding frequency is:

The minimum integer greater than $N\times f_a/f_s$ is denoted as $k_a$. The maximum integer less than $N\times f_b/f_s$ is denoted as $k_b$. If $\mathbf{x}\left(n\right)$ is a vibration acceleration signal, then its vibration intensity in the frequency range of $f_b{f}_b$ can be calculated as:

Fig. 4a) Raw vibration signals and b) their vibration intensities

a)

b)

The vibration intensity of some specific locations can be used as a criterion to judge whether equipment is operating normally. The vibration intensity is in a low level and stable state when the system is operating normally. Once fatigue damage or other faults occur at a certain component, the device enters a depletion phase and the fault will be reflected as the increase of vibrational energy. When vibration intensity increases to a certain extent, it will affect the operational safety of the system. Consequently, we can calculate the bearing vibration intensity with raw vibration data to obtain the degradation feature for RUL prediction.

Among all the data, the vibration intensity of the data of three bearings (Bearing 1-1, 1-3 and 1-4, respectively referred to as #1, #3, and #4) exhibits a gradual degradation trend, as shown in Fig. 4. These three sets of data are chosen to verify the improved PF prediction method. With the three sets of data, a vibration intensity value of 0.1 mm/s was regarded as the initial point of the prediction and the failure threshold is set to be 1.1 mm/s. Then, through the pretreatment of smoothing and monotonizing, the final degradation features for PF prediction algorithm are obtained, as shown in Fig. 5.

Fig. 5Degradation features after pretreatment

4.2. Historical priori information

4.2.2. The distribution map of the degradation rate

Another work before prediction is the calculation of the distribution map of the degradation rate with different features. The distribution information will be used to guide the estimation of the degradation rate in the prediction phase. In the basic PF prediction method, degradation trajectory of all particles are extrapolated directly according to the updateded model parameters. During this process, a partcile would not be adjusted byany other ways, and the corresponding estimation errors of model parameters would accompany its evolution to the end. Whereas, DFT-PF needs specific degradation rate of the feature for extrapolating, that can control the accidental divergence of predictive error. It is not hard to obtain the distribution map of the degradation rate with different features. Firstly, the degradation rates of different historical data are assigned to the corresponding axis of the degradation feature, and then to calculate the degradation trend curve with the interpolation algorithm. Finally, the normal distribution of the degradation rate is fitted at each feature point and the distribution map of the degradation rate with different features is obtained and shown in Fig. 6. Different colors denote the dfferent values, and the bigger value indicate the more concentrated distribution.

Fig. 6The distribution map of the degradation rate with different features

In Fig. 6, the distribution map can be roughly divided into three stages: (1) $\mathbf{x}<\text{0.4}$, the distribution of $\mathbf{a}$ is relatively decentralized at the begining, but as the increase of $\mathbf{x}$, the distribution converges to a narrower range gradually. That is because the manifestations of different bearings differ widely at the early stage of faults, as severity deepening of faults, however, they would become stable. (2) $\mathbf{x}=\text{0.4}\text{0.9}$, the distribution of $\mathbf{a}$ always focuses on a narrow band. It means that the faults of bearings are in a stably degradation condition at this stage. (3) $\mathbf{x}>\text{0.9}$, distribution of $\mathbf{a}$ becomes wider again, and expands like a funnel. It shows that, as long as the faults deteriorate to certain extent, the differences between different bearings or different faults become more and more obvious. Some of them would sped the deterioration until complete failure, while some others may retain previous level of degradation rate. Because of these complex situations during the degradation process of bearings, it is necessary to track the degradation rate for enhancing prediction performance.

According to Fig. 6, in the prediction phase, if degradation feature $\mathbf{x}$ reaches 0.8, the corresponding degradation rate $\mathbf{a}$ focuses on 1.010-1.015; and if $\mathbf{x}$ reaches 1.0, the range of $\mathbf{a}$ enlarges to about 1.01-1.03. We generate $N$ degradation rate values randomly based on the specific distribution obtained before, and substitute these values to prediction model for prediction of every particles at next step.

In addition, the guidance information for degradation rate prediction also includes the value range. According to the statistical results above, the value range of the degradation rate in the whole degradation process, except for a few outliers, is about [1, 1.04]. Through limiting the degradation rate with the algorithm, the particle divergence can be effectively controlled during the calculation.

4.3. Prognostic results and discussion

Some statistical priori information has been obtained with historical data and is helpful in predicting the degradation rate data of an individual sample. Because the details of every set of data are ignored in the obtained statistical priori information, three sets of data above also can be used to validate the method.

Fig. 7(a) depicts the prediction results of 500 observation points for #3 data. Because the relatively fewer data points are for model updating, the PDF range of RUL prediction is broad and the width of 90 % confidence interval is 3415 s. The error between the prediction median and the actual RUL is 770 s. The prediction results of 1000 observation points for #3 data are illustrated in Fig. 7(b). With the increase in the data points for updating, 90 % confidence interval range becomes smaller and the error between the prediction median and the actual RUL decreases to 550 s. In Fig. 7, the thick solid line denoting the prediction median has been traced to the nonlinear trend of the degradation feature. However, some great oscillations occur between 10000 s and 12000 s, causing a step jump of the degradation feature after pretreatment. Therefore, these details cannot be tracked and the prediction error becomes larger.

Fig. 7Prediction results of a) 500 and b) 1000 observation points for #3 data

a)

b)

Fig. 8 illustrates the statistics of prediction results of various observation points. Fig. 8(a) shows the statistical prediction results for #1 data, which start from 1400 observation points. The RUL is predicted for every 25 points. The dotted line denotes the actual RUL curve and the asterisk line indicates 90 % confidence interval of RUL prediction at the corresponding time. The asterisk positions indicate the interval boundary and the prediction median. Like #1 data, in #3 data, the RUL is predicted for every 25 points from 800 observation points. The statistical results of the two sets of data are similar. Moreover, the errors between prediction results and the actual RUL gradually decrease with the increase of the known observations, while the prediction accuracy gradually increases. With the increase of observations, the width of the confidence interval becomes narrow, indicating that the prediction uncertainty gradually decreases. Moreover, the prediction results of these two sets of data tend to be lower than actual RUL. Because the oscillation at the initial stage of the raw feature is violent and its effect cannot be eliminated completely through pretreatment, prediction results are significantly lower than the actual value. However, the prediction result is very close to the actual RUL.

For #4 data, the prediction starts from 60 observation points and the RUL is predicted for every 10 points. The statistical results are shown in Fig. 8(c). Unlike #1 and #3 obviously, the results between 100 to 180 observation points greatly exceed the actual RUL. After that, the results converge to the actual RUL, and the corresponding confidence interval shrink to a very narrow range quickly. After that, the prediction and their variance are smooth and stable until the bearing failure. The reason caused this phenomenon is the particularity of #4 data, not the method itself. As shown in Fig. 9, the evolution trend ofthe degradation feature points of #4 data which less than 100 or more than 180, have better consistency with the diminishing trend of RUL, and compare with #1 and #3 data, they have better convergence and less fluctuation. But there is an 800 s “platform” phase of feature in the range of 100 to 180 points, which would produce a misleading effect on the calculation and cause prediction errors. When the inspection data number is in the range of the “platform”, such as 150 points, it is more reasonable to estimate that the degradation rate after would be close to the latest one. It means that current prediction results can not capture the turning point from the “platform” to ascending stage, so the predicted RUL is much more than actual RUL definitely. Obviously, this problem cannot be avoided in most of the prediction methods [14]. As long as the inspection data number is out of the “platform” period, the misleading effect disappeared. Moreover, because the prediction features after 2000 s retain a better level of convergence and consistency, so that all particles can rapidly converge to the actual RUL, the results become very satisfactory.

Fig. 8Prediction statistical result of a) #1, b) #3 and c) #4 data

a)

b)

c)

Compare with some traditional prediction approaches, particle filtering has more intuitive physical meaning, and could provide the crucial probability information of prediction results for decision maker. But for basic PF method, besides the shortcomings mentioned in Section 3.1, it lacks extra error control in the prediction phase. So it is usually powerless for the long-term prediction or other situations which prediction feature is dissatisfactory. During the experiment verification, we tried to process the same data with basic PF method. However much time is taken for determination of prediction model and initial distribution of model parameters, but we always can not obtain acceptable results, sometimes the particles are even no convergence. The proposed DRT-PF, through the dedicated pretreatment for degradation features, effectively improve their convergent behaviour and trend consistency, and effectively control the prediction error due to the strategy of degradation rate tracking. Though these measures can not fully remove the influences of the “platform”, they are important to optimize the degradation features extracted from vibration signal of bearings. From the case of #4 data, we can further determine the significance of the quality of degradation features to a prediction method. The good consistency between features and degradation trend of the component could ensure the accuracy of the results, and the convrgence of feature decide the convrgence of resultsin great degree. Therefore, in order to improve the ability of RUL prediction for mechanical components, besides improvement of method itself, trying to find a kind of more satisfactory prediction feature is equally important.

Fig. 9Analysis of predictive degradation trend for #4 data at different observation points