Time-Varying Degradation Model for Remaining Useful Life Prediction of Rolling Bearings under Variable Rotational Speed

Abstract

It is difficult to accurately extract the health index of non-stationary signals of rolling bearings under variable rotational speed, which also leads to greater prediction error for bearing degradation models with fixed parameters. For this reason, an angular domain unscented particle filter model with time-varying degradation parameters is proposed to deal with the remaining useful life (RUL) prediction of rolling bearings. Order analysis is first performed to transform the variable-speed signal from time domain to angular domain for extracting the health index in the angular domain representation. To track the bearing degradation state, a real-time finite element model is established to guide the parameters updating the procedure of the performance degradation model. Finally, the bearing degradation state is estimated by the unscented particle filter (UPF), and then the remaining useful life of the bearing is predicted. In this way, the time-varying degradation model is developed by considering both non-stationary feature extraction and dynamic state tracking for rolling bearings. The proposed method is verified by both benchmarks: bearing experimental data, and a bearing accelerated life experiment. Compared with state-of-the-art prognostic methods, the present model can predict the bearing remaining useful life (RUL) more accurately under variable rotational speed.

1. Introduction

As one of the important parts of rotating machinery, rolling bearings need to con-sume a lot of manpower and material resources for maintenance in the case of failure [1,2,3]. Generally, the working conditions of the bearing are complex, and the working environment is harsh. The performance and health of the bearing will inevitably degrade during the production process [4,5,6]. If real-time monitoring of rolling bearing operating status and remaining useful life (RUL) prediction can be performed, early fault warning changes post-maintenance into condition-based maintenance. This can effectively reduce maintenance costs and avoid major safety accidents [7,8,9].

RUL prediction methods include physical model-based and data-driven ones [10,11]. The data-driven method combines historical data and condition monitoring data with machine learning technology, such as artificial neural networks [12], support vector machines (SVMs) [13], and Gaussian process regression [14]. It establishes a prediction model to predict RUL. The advantage of data-driven methods is that they can directly understand the potential degradation trend of bearings from the available sensor data. As such, users do not need to know the exact failure mechanism of bearings. Soualhi et al. [15] used the Hilbert–Huang transform to extract fault features. The obtained characteristic frequency was used to construct prediction index that was fed into SVMs for RUL prediction. Shivani et al. [16] used an overall empirical mode decomposition strategy to calculate time domain features, and constructed a support vector regression model to predict the health status of rolling bearings. Ali et al. [17] developed a prediction method based on an artificial neural network for bearing RUL prediction. However, the above data-driven method mainly establishes the relationship between the degradation state and the system through a large amount of historical data. The prediction accuracy depends not only on the quantity of historical data, but also on the quality of historical data. In the process of RUL prediction, it is assumed that the decline process of the system follows certain rules without any mutation. Therefore, it is not suitable for RUL prediction of rolling bearings under varying working conditions and other uncertain factors.

On the other hand, the model-based method depends on a physical model to describe the overall decline trend of the bearing through the failure mechanism, and then uses statistical estimation techniques, such as the Kalman filter (KF) [18], particle filter (PF) [18], and the unscented particle filter (UPF) [19] to update and estimate the bearing state and predict the bearing RUL. Lei et al. [20] used the maximum likelihood algorithm to initialize the parameters of the Paris–Erdogan model, where the PF model was employed to recursively predict the bearing state. Yan et al. [21] used the multipored autoregressive model to establish the bearing degradation model, and combined this with enhanced PF for bearing RUL prediction. Kan et al. [22] constructed a nonlinear state space model, and used the unscented Kalman filter (UKF) algorithm to update and estimate the state, which realized the bearing performance evaluation and RUL prediction. The main advantage of the model-based method is that the prediction results are more intuitive, and the full combination of expert knowledge and real-time information from the machine can effectively carry out the bearing RUL prediction. However, the above methods mainly aim for the RUL prognostics of bearings under constant conditions, and they are not involved in varying working scenarios, such as variable rotational speed.

The Kalman filter (KF) [18] is an optimal recursive data Processing algorithm, which can optimally estimate the future state of dynamic system. However, KF is more effective for a linear system, and is not suitable for dealing with serious nonlinear signals. KF assumes a linear system dynamic model with Gaussian noise in measurement, which is not always realistic in practical application. The extended Kalman filter (EKF) [23] is an extension of the Kalman filter to nonlinear system dynamics. EKF linearizes the nonlinear equation locally through the Taylor series expansion formula, so as to adapt to the nonlinear system. However, EKF must solve the Jacobi matrix of nonlinear function, which is complex and error prone. For systems with strong nonlinearity, the stability is relatively poor, which is easy to lead to the decline of filtering results, and can only meet the requirements of some weak nonlinear environments. The unscented Kalman filter (UKF) [22] is a combination of the unscented transform and the standard Kalman filter system. It uses the linear regression of sigma points to linearize the nonlinear function of variables, avoid the operation of derivatives, and make it more suitable for real strong nonlinear systems. However, in numerical calculations, UKF encounters some problems, such as rounding error, and negative definite covariance matrix, which leads to filter divergence. The derivative-free nonlinear Kalman filter [24] is a QR decomposition of the state variance matrix in the unscented Kalman filtering process, which can ensure the positive definiteness of the covariance matrix and overcome the divergence problem caused by the accumulation of calculation errors in the filtering processing No matter how the KF is improved, it still estimates the system under the linear and Gaussian assumptions. When the linearity of the system is very poor, the effect of the Kalman filter will be greatly reduced. Different from the Kalman filter method, particle filter (PF) [18] uses a group of particles with specific weights to approximate the state estimation. This method is not limited to Gaussian hypothesis and can approximate any system state distribution. It is especially suitable for nonlinear and non-Gaussian application scenarios, and has been widely used in residual life prediction. The standard PF algorithm uses the prior distribution of the system as the importance sampling distribution. After a certain number of iterative steps, the weight of many particles becomes minimal, resulting in particle degradation and further affecting the accuracy of residual life prediction. The unscented particle filter (UPF) [19] combines the advantages of the particle filter and the unscented Kalman filter, uses the UKF algorithm to obtain an appropriate importance sampling probability distribution, and makes full use of the latest observations, which can effectively solve the problem of particle degradation.

In practical engineering applications, the environmental conditions of the bearing may change in real time, and its own operating conditions may also change, i.e., variable working conditions. Affected by environmental conditions, operating conditions and other uncertain factors, a variety of failure modes may occur during bearing operation, resulting in non-stationary vibration signals. When using the time-frequency analysis method to analyze non-stationary signals, there are some problems, such as low time-frequency resolution and false component interference, which make it difficult to extract and identify fault features. Under the condition of time-varying speed, the vibration impact interval and impact amplitude of the bearing are time-varying, resulting in the distortion of the vibration signal envelope spectrum obtained by the conventional envelope demodulation method, which makes it difficult to identify the fault characteristics. Order analysis [25], generalized demodulation (GD) [26], time frequency representation and phase-space dissimilarity measure [27] are common methods for vibration signal analysis under variable working conditions, which can accurately extract the degradation characteristics of the bearing life cycle, and are suitable for predicting the remaining useful life of rolling bearings under variable working conditions. However, the existing research mostly focuses on RUL prediction under single working condition, ignoring the consideration of environmental conditions and operating conditions to a certain extent. Therefore, it lacks generalization in practical engineering application, which affects the effectiveness of prediction.

To sum up, it is difficult to accurately extract features containing bearing degradation information from non-stationary vibration signals in a time domain. In addition, the establishment of the degradation model usually requires clear prior knowledge and extensive empirical data. The degradation model is not updated in real time according to the actual degradation state of the bearing, which has a large error with the actual degradation trend of the bearing and affects the RUL prediction accuracy.

To solve the above problems, this paper proposes a RUL prediction method dealing with rolling bearings under variable rotational speed, based on angular domain transformation and time-varying degradation model parameters. Firstly, the angle domain transform is used to process the original time domain signal and extract the real-time degradation features of bearings. Secondly, a finite element method (FEM)-based time-varying parameter degradation model consistent with the bearing performance degradation state is established. The unscented particle filter algorithm (UPF) is used to estimate the bearing degradation state for predicting its RUL. This method not only alleviates the limitations of RUL prediction methods based on data-driven and physical models, but also predicts the future degradation trend of bearings more accurately.

There are three main contributions of this paper. (i) Aiming at the difficulty in extracting features of non-stationary signals under variable operating conditions, this paper proposes an angular domain unscented particle filter (UPF) RUL prediction method. The angle domain transform is used to process the signals, which can accurately extract the features containing important information, and combined with the UPF algorithm to update and estimate the bearing state, thereby improving the accuracy of the remaining useful life prediction of the bearing. (ii) In the process of establishing the bearing degradation model based on the Paris–Erdogan theorem, the parameters of the model are not selected according to empirical formulas, but by establishing a finite element model consistent with the bearing degradation model to update model parameters in real time, to make it more in line with the real degradation trend of the bearing. (iii) The effectiveness of the proposed method is verified by two bearing experiments. Compared with the state-of-the-art methods, the developed method obtained a superior performance.

The rest of this paper is arranged as follows. Section 2 introduces the basic theory of angular domain transformation, the unscented particle filter, and the time-varying parameter degradation model. Effectiveness of the proposed method is verified by a benchmark dataset experiment and an accelerated life prognosis experiment as presented in Section 3, respectively. Conclusions are obtained in Section 4.

2. Materials and Methods

2.1. Angular Domain Transform

Signals collected by the accelerometer and the tachometer are employed as the input of the angular domain transform. Numerical integration is then performed to convert the rotational speed pulse signal into a turning angle signal, which is then converted by equal-angle resampling [28]. A non-stationary signal in the time domain can be converted into a stationary one in the angular domain. Then the characteristic indexes (kurtosis, root mean square, etc.) are extracted from the angle domain signal. The time interval of equal angle sampling is determined by the rotation angle of the reference axis. Therefore, the time domain non-stationary signal is converted into the angular domain stationary signal by using the sampling method of computed order tracking and equal angle incremental resampling [29]. In a short time, the rotation axis accelerates uniformly at an angle, and the relationship between the rotation angle and time is given by [30].

$$\theta \left(t\right)=a_0+a_{1}t+a_2{t}^2$$

(1)

where $a_0$, $a_1$ and $a_2$ are polynomial coefficients.

$$\left[\begin{array}{c}\theta \left(t_1\right)\\ \theta \left(t_2\right)\\ \theta \left(t_3\right)\end{array}\right]=\left[\begin{array}{ccc}1& t_1& t_1{}^2\\ 1& t_2& t_2{}^2\\ 1& t_3& t_3{}^2\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\end{array}\right]$$

(2)

Supposing the angular interval between two adjacent key phase pulses be a fixed value $2\pi$, according to Equation (2), the undetermined coefficients can be obtained. By substituting Equation (2) into Equation (1), the time corresponding to the rotation angle change can be obtained as

$$t\left(\theta \right)=\frac{1}{2a_2}\left[\sqrt{4a_2\left(\theta -a_0\right)+a_1{}^2}-a_1\right]$$

(3)

The rotation angle is discretized according to the equal angle sampling interval. The above equation is therefore recast as

$$T_n=\frac{1}{2a_2}\left[\sqrt{4a_2\left(n\Delta \theta -a_0\right)+a_1{}^2}-a_1\right]$$

(4)

where the time of equiangular sampling point is $T_n$ and $\Delta \theta$ is the sampling interval.

According to the equal angular sampling time, the Lagrange interpolation operation is performed on the signal $y(t)$, and then the resampled angular domain signal is obtained.

2.2. Unscented Particle Filter

The degradation law of dynamic system state can be described by a set of dynamic models including a state equation and observation equation expressed by [19].

$$x_t=f_t(x_{t-1},w_t)$$

(5)

$$y_t=h_t(x_t,v_t)$$

(6)

where $f_t$ is the state transfer function from the previous time to the current time, $x_t$ and $w_t$ are the state value of the system and the system noise at the current time, $h_t$ is the functional relationship between the current state and the observed value, $y_t$ is the observed value of the current system state, and $v_t$ is the observation noise at the current time.

RUL prediction is to estimate the state value of the next moment from the current state value and the observed value of the rolling bearing. According to Bayesian theory, all information of the system state estimation is included in the posterior probability distribution, given by:

$$p\left(x_t|y_{1:t-1}\right)={\displaystyle \int p\left(x_t|x_{t-1}\right)}p\left(x_{t-1}|y_{1:t-1}\right)d{x}_{t-1}$$

(7)

$$p\left(x_t|y_{1:t}\right)=\frac{p\left(y_t|x_t\right)p\left(x_t|y_{1:t-1}\right)}{p\left(y_t|y_{1:t-1}\right)}$$

(8)

where $p\left(x_{1:t}|y_{1:t}\right)$ is the posterior probability and $p\left(x_t|y_{1:t-1}\right)$ is the prior probability of the state, which describes the prior knowledge of the state before the observation value is obtained.

The traditional PF algorithm is based on Bayesian theory. It estimates the probability density function (PDF) of the system state through the observation data, and uses the sequential importance sampling (SIS) algorithm to replace PDF with a priori probability approximation. However, the PF algorithm easily falls into the problem of particle degradation, which leads to the divergence of estimation results. The unscented particle filter algorithm (UPF) uses UKF to construct the importance sampling distribution of the particle filter. It integrates the latest observation information into the importance distribution function. It strengthens the tracking ability of particles, improves the accuracy of target tracking, and effectively solves the problem of particle degradation.

The first step is the initialization. Setting t = 0, sample N particles from the initial distribution to generate the original particle set $x_0^i=p\left(x_0\right)$. The weight coefficient corresponding to each particle is

$$w_0^i=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$N$}\right.,i=1,2,N$$

(9)

For each particle, UKF algorithm is used to estimate the state of each particle. The particle $x^i$ is obtained according to the Gaussian distribution $N(\stackrel{-}{x_t^i},\stackrel{\wedge }{p_t^i})$ as

$$\stackrel{-}{x_t^i}={\displaystyle \sum _{i=1}^N{w}_t^i}{x}_t^i$$

(10)

$$\stackrel{\wedge }{p_t^i}={\displaystyle \sum _{i=1}^N{w}_t^i}\left(x_k^i-\stackrel{-}{x_k^i}\right){\left(x_k^i-\stackrel{-}{x_k^i}\right)}^T$$

(11)

The weight of each particle is updated according to the measured value at the current time. The weight is normalized as

$$w_t^i=w_{t-1}^i\frac{p\left({\mathrm{y}}_t|{\mathrm{x}}_t^i\right)p\left({\mathrm{x}}_t^i|{\mathrm{x}}_{t-1}^i\right)}{q\left({\mathrm{x}}_t^i|{\mathrm{x}}_{0:t-1}^i,{\mathrm{y}}_{1:t}\right)}$$

(12)

$$w_t^i=\frac{w_{t-1}^i}{{\displaystyle \sum _{i=1}^N{w}_t^i}},i=1,2,N.$$

(13)

$$N_{\mathrm{eff}}=1/{\displaystyle \sum _{i=1}^N{(w_t^i)}^2},i=1,2,N.$$

(14)

where $N_{\mathrm{eff}}\le N_{thr}$ ($N_{thr}$ represents the threshold). This means that the particles have been severely degraded and need to be resampled according to the importance weight, and the particles are mapped into N particles of equal weight. In this way, the updated particle and its weight can be used to estimate the state at the current time as

$$x_t=\frac{1}{N}{\displaystyle \sum _{i=1}^N{W}_t^i}{x}_t^i,i=1,2,N.$$

(15)

2.3. Time-Varying Parameters of the Degradation Model

When the initial fault appears in a bearing, according to the Paris–Erdogan theorem, the relationship between crack propagation and the number of stress cycles can be derived as [31]

$$\frac{dx}{dL}=C{\left(\Delta K\right)}^m$$

(16)

where $x$ is the crack length, L is the number of stress cycles (fatigue life), C and m are constants related to materials, and $\Delta K$ is the amplitude of stress intensity factor, which changes with the propagation of bearing cracks.

To obtain the relationship between the stress intensity factor and the crack length in the actual degradation process of bearing, a finite element method (FEM) is used to analyze the crack of rolling bearing, as shown in Figure 1.

Firstly, an experimental bearing FEM is established, in which the crack length and depth information are set to simulate the fatigue crack propagation. Constraints and loads are imposed according to experimental conditions. Because the shape function of ordinary element cannot represent the singularity of stress and strain at the crack front, a 12-node singular element is used for the tetrahedral mesh. A singular element is a kind of quadratic element whose nodes are located at the quarter, with its shape function given by

$$u=2\left(\sqrt{v}/\sqrt{l}\right)\left[\left(1-\eta \right)u_2+0.5u_3\right]$$

(17)

$$\frac{\partial u}{\partial v}=\frac{1}{\sqrt{v}}\ast \frac{1}{l}\left[-2\eta u_2+\left(0.5+\eta \right)u_3\right]$$

(18)

where $u$ and $v$ are the node coordinate, $u_2$ and $u_3$ are the corresponding node coordinate, $\eta$ is a constant, and l is the unit side length.

The finite element analysis is used to simulate the crack propagation, and a set of stress intensity factors varying with the crack length are obtained by J-integral analysis to simulate the fatigue crack propagation. The mathematical relationship of $\Delta K$ with the square root of the crack length is obtained by numerical interpolation fitting. The time-varying physical parameters obtained from the Paris-Erdogan theorem and FEM are transferred to the degradation model of the bearing. The stress intensity factor is taken as the time-varying physical parameter in the bearing degradation model.

2.4. The Present RUL Prediction Method

As shown in Figure 2, the addressed time-varying degradation model for bearing RUL prediction can be divided into the following eight steps.

Step 1: The original non-stationary time domain signal is transformed into the angular domain stationary signal by equal angle resampling.

Step 2: The bearing life data are transformed into the angular domain under variable speed conditions towards a complete angular dataset. The stable angular domain signal of the bearing is extracted at each time, and the degradation characteristic index (kurtosis, root mean square, etc.) of the rolling bearing can be obtained.

Step 3: Kurtosis is used to monitor the occurrence of bearing faults. The starting point of bearing RUL prediction is determined. The specific process is as follows. Firstly, the mean value $\mu$ and standard deviation ${\sigma }_k$ of kurtosis under the normal working condition of the bearing are calculated. Then, the kurtosis interval $\left[\mu -3{\sigma }_k,\mu +3{\sigma }_k\right]$ is defined as the normal working interval of the bearing. When the kurtosis of the bearing at a certain time exceeds the defined normal interval, it is judged that the bearing enters the degradation stage. This time is regarded as the starting point for the bearing RUL prediction is then started.

Step 4: The envelope analysis of the vibration signal at the starting point of the prediction of the RUL of the bearing is carried out to obtain the fault characteristic. According to the results of envelope analysis, the fault frequency is obtained, and the bearing fault location is determined by calculation. Finally, the initial crack at the bearing fault location is set in the finite element model.

Step 5: The relative root mean square value (RRMS) is employed as the input of the degradation model to predict the bearing RUL. The specific definition of RRMS is expressed by

$$X_{RRMS}\left(t\right)=\frac{X_{RMS}\left(t\right)-X_{RMS}\left(T\right)}{X_{RMS}\left(T\right)}$$

(19)

where $X_{RRMS}\left(t\right)$ is the RRMS value of the signal at time t, $X_{RMS}\left(t\right)$ is the RMS value of the signal at time t, and $X_{RMS}\left(T\right)$ is the RMS value of the signal at RUL starting time T.

Step 6: The real-time finite element model of the bearing is used to obtain the numerical relationship between physical parameters $\Delta K$ and the crack length $x$.

$$\Delta K=\Delta \sigma \sqrt{x}$$

(20)

Through integrating L above both sides of Equation (16), the function can be obtained by

$$\begin{array}{l}\frac{x_t-x_{t-1}}{L_t-L_{t-1}}=\alpha x_{t-1}{}^{\beta }\\ x_t=x_{t-1}+\alpha x_{t-1}{}^{\beta }\Delta L\end{array}$$

(21)

where $\alpha =C{\left(\Delta \sigma \right)}^m,\beta =3m/2$, C and m are constants related to the material, $\Delta \sigma$ is a constant, $x$ is the crack length, and L is stress cycles.

The relationship between stress cycles L and the running time t of the bearing is linear, and the coefficient is simplified to 1. According to references [32,33,34], under ideal conditions, the characteristic index of vibration signal has a linear relationship with the crack length. In order to simplify, the coefficient is set to 1. Combining Equations (5), (6) and (21) leads to the state transition equation and observation equation in the bearing degradation process as

$$\begin{array}{l}{x}_t=x_{t-1}+\alpha x_{k-1}^{\beta }\Delta t+w_t\\ y_t=x_t+v_t\end{array}$$

(22)

where $w_t,v_t$ are the system random noise, $\Delta t$ is the time interval, $x_t$ is state value, and $y_t$ is the measured value.

Step 7: The degradation index (RRMS) is brought into the degradation model by using the steps of the UPF algorithm, and the model parameters are updated by the measured values to estimate the degradation index value at the current time. The degradation model is used to transfer the existing distribution and predict the degradation index at the next time.

Step 8: The predicted value of the degradation index RRMS and its time to reach the failure threshold are further mapped to the RUL of the rolling bearing.

3. Results

3.1. Experimental Study of a Benchmark Bearing

3.1.1. Benchmark Dataset

The bearing dataset [35] is a benchmark provided by the Xi’an Jiao tong University (XJTU) and the Sum young Technology Co., Ltd. (SY) (Sum young, Huzhou, China), it is defined as XJTU-SY. As shown in Figure 3, the bearing testbed is composed of an alternating current (AC) induction motor, a motor speed controller, a support shaft, two support bearings (heavy duty roller bearings), a hydraulic loading system and so on. This testbed was designed to conduct the accelerated degradation tests of rolling element bearings under different operating conditions (i.e., different radial force and rotational speed). The radial force was generated by the hydraulic loading system and applied to the housing of tested bearings. The rotational speed was set by the speed controller of the AC induction motor.

The type of tested bearings is LDK UER204, and the detailed parameters are given in

Table 1. Parameters of the tested bearings.

Parameter	Value	Parameter	Value
Outer race diameter	39.80 mm	Inner race diameter	29.30 mm
Bearing mean diameter	34.55 mm	Ball diameter	7.92
Number of balls	8	Contact angle	0 rad
Load rating (static)	6.65 kN	Load rating (dynamic)	12.82 kN

. A total of three different operating conditions were set in the accelerated degradation experiments. Under each operating condition, five bearings were tested. The operating conditions include: 2100 rpm (35 Hz) and 12 kN; 2250 rpm (37.5 Hz) and 11 kN; 2400 rpm (40 Hz) and 10 kN. The dataset contains complete run-to-failure data of 15 rolling element bearings acquired by conducting many accelerated degradation experiments. The sampling frequency is 25.6 kHz, a total of 32,768 data points (i.e., 1.28 s) are recorded for each sampling, and the sampling period is equal to 1 min. Five bearings in condition 1 are selected for verification.

Table 2. XJTU-SY bearing dataset.

Operating Condition	Bearing Dataset	Number of Files	Bearing Lifespan	Fault Element
Condition1 (37.5 Hz/ 11 kN)	Bearing1_1	123	2 h 3 min	Outer race
	Bearing1_2	161	2 h 41 min	Outer race
	Bearing1_3	158	2 h 38 min	Outer race
	Bearing1_4	122	2 h 2 min	Cage
	Bearing1_5	52	52 min	Inner and outer race

lists the detailed information of each test bearing, including the number of CSV files, bearing life and fault elements.

3.1.2. Results Analysis

As shown in Figure 4, health indicators including root mean square and kurtosis of five bearings are used to monitor the degradation starting time. The relative root mean square (RRMS) of the moment is taken as the failure threshold of the bearing. Bearing 1_1 begins to degrade from the 78th minute, bearing 1_2 begins to degrade from the 126th minute, bearing 1_3 begins to degrade from the 110th minute, bearing 1_4 begins to degrade from the 80th minute, and bearing 1_5 begins to degrade from the 39th minute. They are taken as the starting points of RUL prediction, respectively. The degradation characteristic index RRMS is introduced into the time-varying parameter degradation model. The prediction results are shown in Figure 5, which are further mapped to the rolling bearing RUL through RRMS.

3.1.3. Discussion and Comparison

To show the superiority of the developed bearing degradation model in this paper, it is compared with the fixed degradation model established by the traditional empirical formula. The degradation characteristic index RRMS is used as the input of the model to calculate the state prediction values of five bearings, respectively. In order to quantitatively evaluate the performance of the proposed method, the prediction error is designed and calculated, which can be evaluated by the root-mean-square error (RMSE) using the following equation:

$$RMSE=\sqrt{\frac{1}{M}{\displaystyle \sum _{i=1}^M{\left(x_p(i)-x_r(i)\right)}^2}}$$

(23)

where $x_p$ is the predicted data point, $x_r$ is the real data point, and M is the number of predicted data points.

The standard deviation of the state estimation error (RMSE) of the five bearings is calculated as shown in Figure 6. Taking the result of bearing1_1 as an example, the RMSE of the degradation model based on time-varying parameters is 0.7123, while that of the fixed degradation model based on empirical formula is 0.8323. This shows that the performance of the degradation model based on time-varying parameters is better. The state prediction results of the two models are mapped to the RUL of the bearing as shown in Figure 7. The RUL prediction result of bearing1_1 is also taken as an example to illustrate. The RUL prediction error of the degradation model based on time-varying parameters is 10 min, while that of the fixed degradation model based on empirical formula is 16 min. Therefore, the prediction results of the proposed model are more convergent than those of the traditional empirical model. The accuracy is higher and is closer to the actual RUL.

To illustrate the advantages of the proposed time-varying model compared with the data-driven method and the fixed-model based method, the same dataset was compared with four advanced methods published in the literature, including the mixed prediction method based on RVM and the random degradation model [35], the method based on RVM [36], the method based on PF [37], and the method based on EKF [38]. To quantitatively evaluate the prediction performance of these five forecasting methods, a widely-used evaluation index, namely cumulative relative accuracy (CRA), is used in this paper. CRA is able to comprehensively evaluate the accuracy of the prediction method by summarizing the relative prediction accuracy of all forecast times. It is calculated as

$$CRA\left(T_K\right)=1-\frac{\left|ActRUL\left(T_k\right)-RUL\left(T_k\right)\right|}{ActRUL\left(T_k\right)}$$

(24)

where $ActRUL\left(T_k\right)$ is the actual life in the test process, and $RUL\left(T_k\right)$ is the predicted life. The closer the CRA value is to 1, the more accurate the RUL estimation result of the prediction method is.

Table 3. CRA comparisons between five methods.

Dataset	RVM	PF	EKF	PHPA	This Model
Bearing1_1	0.5741	0.6107	0.6209	0.9047	0.9186
Bearing1_2	0.1815	0.7256	0.3500	0.8546	0.8992
Bearing1_3	0.6245	0.4850	0.8010	0.8482	0.8663
Bearing1_4	0.3722	0.2305	0.6839	0.7240	0.7976
Bearing1_5	0.6122	0.4311	0.5042	0.7878	0.8293

shows the prediction results of the five methods. Taking Bearing1_1 as an example, CRA obtained by the proposed method is 0.9186, while that obtained by the method based on RVM is 0.5741, that obtained by the method based on PF is 0.6107, that obtained by the method based on EKF is 0.6209, and that obtained by the method based on PHPA is 0.9186. Compared with other four state-of-the-art methods, the proposed prediction model in this paper has the best performance and the highest prediction accuracy.

3.2. Run-to-Failure Bearing Experiments

3.2.1. Experimental Set-Up

To illustrate the effectiveness of the method proposed in this paper, the vibration data of a run-to-failure experiment of a rolling element bearing were collected by a test bench, as shown in Figure 8.

The test bearing is a single-row deep groove ball bearing (ERK16) whose main geometric parameters are shown in

Table 4. Main geometric parameters of rolling bearing ERK16.

Type	Pitch Diameter	Ball Diameter	Contact Angle	Rated Dynamic Load	Numbers of Rollers
ERK16	D (mm)	d (mm)	(rad)	Cr (kN)	Z
	39.04	7.94	0	14	9

. According to different working conditions of the bearing, experiments were carried out under both constant speed and variable speed. This paper mainly discussed the RUL prediction of the rolling bearing under variable rotational speed. Hence, the test conditions are as follows: radial load 8500 N rotational speed fluctuation range [1450, 1550] r/min, sampling frequency 12,800 Hz, sampling interval 5 min, and sampling duration 4 s.

3.2.2. Experimental Set-Up

The variable speed experimental data are transformed in the angular domain. The spectrum of the original non-stationary signal is shown in Figure 9a, and the stationary signal after angular domain transformation is shown in Figure 9b.

As shown in Figure 10a, the characteristic indexes (kurtosis, root mean square) are extracted from the angle domain signal. The degradation start time of the bearing is 1150 min, and the failure time is 1450 min. The RRMS of each moment was calculated, and the degradation characteristics of rolling bearings were predicted by the method presented in this paper. It was further mapped to the bearing RUL, as shown in Figure 10b.

To explain the influence of different fault locations on the stress intensity factor of time-varying model parameters, a plane open crack with length of 1 mm and depth of 1 mm was set in the inner race and outer race of the bearing in FEM. According to the experimental conditions, the inner race of the bearing was fixed, and a radial load of 8500 N was applied to the bearing. The relationship between the stress intensity factor and the square root of the crack length is shown in Figure 11. It can be seen that the crack location has little effect on the stress intensity factor. $\Delta K=\Delta {\sigma }_1\sqrt{x}=5.675\times {10}^5\times \sqrt{x}$ for inner race, and $\Delta K=\Delta {\sigma }_2\sqrt{x}=5.753\times {10}^5\times \sqrt{x}$ for outer race were set in the experiments.

To illustrate the superiority of the bearing degradation model based on time-varying parameters proposed in this paper, RRMS extracted directly from the original time domain was introduced into the bearing degradation model based on the traditional empirical formula. RMSE of state estimation error is obtained as shown in Figure 12. RMSE obtained by directly extracting features from time domain is 0.8327, while the error obtained by extracting features after angular domain transformation becomes 0.2596. Therefore, it can be seen that angular domain transformation has better effect in Processing non-stationary time domain signals with variable rotational speed. RMSE of the degradation model based on the traditional empirical formula is 0.4708, while RMSE of the degradation model based on the time-varying parameters is reduced to 0.2596. Therefore, the degradation model based on the time-varying parameters proposed in this paper is closer to the actual degradation state of the bearing with high prediction accuracy.

We further verify the effectiveness of the developed method in bearing RUL prediction. The predicted degradation characteristic RRMS is mapped to the RUL of the bearing. The time when the bearing reaches the failure threshold is calculated, as shown in

Table 5. Comparison of RUL prediction results of four methods.

	Angle Domain Fixed Model	Developed Model	Time Domain Fixed Model	Time Domain Time-Varying Model
Actual RUL	1450 min	1450 min	1450 min	1450 min
Predicted RUL	1425 min	1462 min	1414 min	1423 min
Error	25 min	12 min	36 min	27 min
CRA	0.9824	0.9917	0.9751	0.9813

. The prediction error of RUL based on time domain analysis and the traditional fixed model is 36 min, while the prediction error of RUL is reduced to 27 min after angular domain transformation. The prediction error of RUL based on angular domain transformation and the empirical fixed degradation model is 25 min, while the prediction error of RUL based on angular domain transformation and time-varying parameter degradation model is 12 min. Therefore, the prediction accuracy of the bearing degradation model based on time-varying parameters proposed in this paper is higher than that of the fixed degradation model based on the traditional empirical formula. Figure 13a shows the life prediction results of two different models under time domain analysis. It can be seen that the prediction results of the degradation model based on time-varying parameters are closer to the real life at each step, and the prediction results are more convergent. Figure 13b shows the life prediction results of two different models after angular domain transformation. It can also be seen that the prediction results of the degradation model based on time-varying parameters are closer to the real life at each step. The prediction results are more convergent than the degradation model based on the empirical formula. These experimental results further verified that the degradation model based on angular domain transformation and time-varying parameters is closer to the real degradation process of the bearing.

4. Conclusions

A time-varying degradation model has been reported in this paper for predicting RUL of bearings under variable rotational speed. For non-stationary time-domain signals, it remains challenging to extract the degradation features accurately. Bearing degradation models established by theoretical knowledge have great error with the real degradation processing For this reason, angular domain transformation and time-varying parameters were integrated as the time-varying degradation model. The non-stationary temporal signal of bearing lifecycle under variable speed was transformed as stationary angular domain representation. Real-time FEM of the bearing was established to track the bearing performance degradation state. The parameters of the bearing performance degradation model were updated according to the time-varying physical commitment calculated by FEM. The degradation characteristic index is brought into the degradation model by using the steps of the UPF algorithm, and the model parameters are updated by the measured values to estimate the degradation characteristic value at the current time. The degradation model is used to transfer the existing distribution and predict the degradation characteristic at the next time. Finally, the remaining useful life of the bearing is predicted by calculating the time when the predicted degradation eigenvalue reaches the failure threshold. For the XJTU-SY experiment, compared with state-of-the-art prognostic method, the accuracy of the present model is improved by 5.3%. For the BPS experiment, compared with the traditional fixed model, the prediction accuracy of the present model is improved by 1.7%. Therefore, the present model can predict the bearing remaining useful life (RUL) more accurately under variable rotational speed.

However, the proposed method still has some limitations. The developed method is only suitable for one-dimensional cracks. Regarding multi-dimensional cracks, which is a crucial and important issue for RUL prognosis, this will be studied in the future.