State Space Formulation of Nonlinear Vibration Responses Collected from a Dynamic Rotor-Bearing System: An Extension of Bearing Diagnostics to Bearing Prognostics

Abstract

Bearings are widely used in various industries to support rotating shafts. Their failures accelerate failures of other adjacent components and may cause unexpected machine breakdowns. In recent years, nonlinear vibration responses collected from a dynamic rotor-bearing system have been widely analyzed for bearing diagnostics. Numerous methods have been proposed to identify different bearing faults. However, these methods are unable to predict the future health conditions of bearings. To extend bearing diagnostics to bearing prognostics, this paper reports the design of a state space formulation of nonlinear vibration responses collected from a dynamic rotor-bearing system in order to intelligently predict bearing remaining useful life (RUL). Firstly, analyses of nonlinear vibration responses were conducted to construct a bearing health indicator (BHI) so as to assess the current bearing health condition. Secondly, a state space model of the BHI was developed to mathematically track the health evolution of the BHI. Thirdly, unscented particle filtering was used to predict bearing RUL. Lastly, a new bearing acceleration life testing setup was designed to collect natural bearing degradation data, which were used to validate the effectiveness of the proposed bearing prognostic method. Results show that the prediction accuracy of the proposed bearing prognostic method is promising and the proposed bearing prognostic method is able to reflect future bearing health conditions.

1. Introduction

In most rotary machines, bearings are widely used in various industries to support the rotating shafts. Their failures accelerate failures of other adjacent components and result in unexpected machine breakdowns. Analyses of nonlinear vibration responses collected from a dynamic rotor-bearing system are widely used to diagnose different bearing faults [1]. A bearing consists of an outer race, an inner race, a cage, and several rollers. When there is a defect on the surface of either the outer race or the inner race, an impulse generated by a roller striking the defect is observed in a nonlinear vibration response collected from the casing of a dynamic rotor-bearing system [2,3,4,5]. Moreover, when a shaft rotates, more and more impulses are observed in the nonlinear vibration response. Experimental results have proved that these impulses are not purely periodic but slightly random, which causes difficulties of bearing diagnostics [6].

In recent years, many advanced methods have been proposed to analyze the nonlinear vibration responses caused by bearing faults [7,8,9]. One of the most attractive methods for bearing diagnostics is spectral kurtosis [10], which aims to calculate kurtosis values of signals obtained by different filters with different frequency supports. Based on the idea of spectral kurtosis, some improvements, such as the improved Kurtogram [11], the enhanced Kurtogram [12], the sparsogram [13], etc., have been designed to detect bearing faults in different cases. However, these methods [10,11,12,13] can only provide rough results for bearing fault diagnosis. To improve their visual inspection performances and achieve optimal filtering [14], different optimization algorithms, such as genetic algorithm [13], differential evolution, etc., and different objective metrics, such as kurtosis [15], entropy [16], smoothness index [17], sparsity measurement [13], etc., have been thoroughly investigated.

Besides bearing diagnostics, bearing performance degradation assessment is an emerging and hot topic in recent years and attracts much attention. The bearing performance degradation assessment focuses on evaluations of the current bearing health condition [18] and its health evolution over time. In other words, the bearing performance degradation aims to assess how far the current bearing health condition deviates from a normal bearing health condition. Some methods, such as wavelet filtering [19,20], hidden Markov model [21], support vector data description [22], Gaussian mixture model [23], etc., have been proposed to track the current bearing health condition. However, the above methods belong to supervised learning. They require sufficient historical data to train the associated models. Therefore, their deficiencies are apparent when historical data are insufficient or unavailable.

Recently, a novel earth mover distance-based bearing performance degradation assessment method [24] was proposed by one of the co-authors to describe bearing deterioration over time. Even though the aforementioned method is effective in evaluating the current health condition of bearings, it has some limitations, illustrated as follows. Firstly, artificial bearing defects were introduced to normal bearings prior to the bearing acceleration life testing. This means that the nonlinear vibration responses collected from the dynamic rotor-bearing system are not natural bearing degradation data, which may not well reflect bearing performance degradation from normal to failure. Secondly, for their proposed health index, it is assumed that one of the bearings used for supporting the shaft is in good condition. In other words, a pair of bearings does not degrade simultaneously. It is natural to argue that such an assumption is not completely true. Thirdly, their proposed method only focuses on evaluations of the current bearing health conditions rather than predictions of the future bearing health conditions and predictions of bearing remaining useful life (RUL).

Particle filtering [25] has been used to solve any nonlinear state space models by using a number of random particles and their associated weights. The key step used in this kind of conventional particle filtering is the design of a proper importance function so as to update the weights of the random particles. Even though conventional particle filtering was used in fault prognosis for gas turbines [26], lithium ion batteries [27] etc., the proper procedure for designing the function to update the weights is still questionable. Unscented transform is a numerical method to find limited deterministic sigma points so as to estimate the mean and variance of a transformed variable. Aiming to target the shortcomings of conventional particle filtering, the use of unscented transform combined with particle filtering or the so-called unscented particle filtering was proposed [28].

According to our recent literature review in fault prognosis, the use of unscented particle filtering in bearing prognosis is seldom reported. This paper reports how unscented particle filtering could be tailor-made to predict the remaining useful life (RUL) of bearing. The results of the feasibility study have been stated in this paper. With the success of the feasibility study, several new contributions have been made to the research field of bearing prognosis. Firstly, a new setup for testing the bearing acceleration life was designed and the bearing’s degradation data was naturally collected from real running bearings. Secondly, based on the bearing degradation data, its nonlinear vibration responses were analyzed and then a bearing health indicator (BHI) was designed to indicate the severity of a damaged bearing. Thirdly, the state space formulation of the nonlinear vibration responses collected from a dynamic rotor-bearing system was designed to intelligently estimate the bearing’s RUL. That is, a state space model of the BHI was constructed to mathematically track the health evolutions of the BHI. With the help of the above methods and the verification results generated from the experiments, the proposed BHI was able to predict the RUL of real operating bearings with good accuracy.

The rest of this paper is organized as follows. In Section 2, the basic theories of the conventional particle filtering and unscented transform are introduced. In Section 3, the design of state space formulation of nonlinear vibration responses collected from a dynamic rotor-bearing system is reported. The procedures to intelligently predict the RUL of bearings are also discussed in this section. In Section 4, the experimental results are presented and verification is provided to demonstrate the effectiveness of the proposed bearing prognostic method. Conclusions are drawn in Section 5.

2. The Basic Theories of Unscented Particle Filtering

2.1. Particle Filtering

Particle filtering [29] aims to use $N_s$ random particles ${\left\{x_k^i\right\}}_{i=1}^{N_s}$ and their associated weights ${\left\{{\omega }_k^i\right\}}_{i=1}^{N_s}$ to implement a recursive Bayesian filter and to characterize the posterior probability density function $p(x_k|z_{1:k})$. Here, $x_k$ is the system state and $z_{1:k}$ are measurements. Therefore, particle filtering provides a solution for solving a linear or nonlinear state space model. Ideally, it is expected to draw random particles from $p(x_k|z_{1:k})$, which can be approximated by:

$$p(x_k|z_{1:k})\approx {\displaystyle \sum _{i=1}^{N_s}{\omega }_k^i\delta (x_k-x_k^i)},$$

(1)

where $\delta (\cdot )$ is the Dirac delta function. In many cases, it is impossible to directly draw random particles from the true posterior density function $p(x_k|z_{1:k})$. Nevertheless, it is possible to draw random particles from another proposal distribution $q(x_k|z_{1:k})$ called an importance function. Suppose that the true posterior density function $p(x_k|z_{1:k})$ is proportional to an analytical function $\pi (x_k|z_{1:k})$. The weight of the $p(x_k|z_{1:k})$ can be represented by [29,30]:

$${\omega }_k^i\propto \frac{p(x_k^i|z_{1:k})}{q(x_k^i|z_{1:k})}\propto \frac{\pi (x_k^i|z_{1:k})}{q(x_k^i|z_{1:k})}.$$

(2)

If the importance function is factorized as [29,30]

$$q(x_k^i,x_{k-1}^i|z_{1:k})=q(x_k^i|x_{k-1}^i,z_{1:k})q(x_{k-1}^i|z_{1:k-1})$$

(3)

and the $p(x_k^i|z_{1:k})$ is represented by [29,30],

$$p(x_k^i|z_{1:k})\propto p(z_k|x_k^i)p(x_k^i|x_{k-1}^i)p(x_{k-1}^i|z_{1:k-1}),$$

(4)

the weight ${\omega }_k^i$ can be iteratively updated by:

$${\omega }_k^i\propto \frac{p(x_k^i|z_{1:k})}{q(x_k^i|z_{1:k})}\propto \frac{p(z_k|x_k^i)p(x_k^i|x_{k-1}^i)p(x_{k-1}^i|z_{1:k-1})}{q(x_k^i|x_{k-1}^i,z_{1:k})q(x_{k-1}^i|z_{1:k-1})}\propto {\omega }_{k-1}^i\frac{p(z_k|x_k^i)p(x_k^i|x_{k-1}^i)}{q(x_k^i|x_{k-1}^i,z_{1:k})}.$$

(5)

If $q(x_k^i|x_{k-1}^i,z_{1:k})=q(x_k^i|x_{k-1}^i,z_k)$ is satisfied, calculation of the weight ${\omega }_k^i$ can be simplified as:

$${\omega }_k^i={\omega }_{k-1}^i\frac{p(z_k|x_k^i)p(x_k^i|x_{k-1}^i)}{q(x_k^i|x_{k-1}^i,z_k)}/\left({\displaystyle \sum _{i=1}^{N_s}{\omega }_{k-1}^i\frac{p(z_k|x_k^i)p(x_k^i|x_{k-1}^i)}{q(x_k^i|x_{k-1}^i,z_k)}}\right).$$

(6)

In terms of Equation (6), the importance function depends on the previous system state $x_{k-1}^i$ and the current measurement $z_k$. However, in many applications of particle filtering, to simplify Equation (6), the importance function is usually chosen as a prior state transition function, namely $q(x_k^i|x_{k-1}^i,z_k)=p(x_k^i|x_{k-1}^i)$. Then, Equation (6) is simplified as:

$${\omega }_k^i={\omega }_{k-1}^{i}p(z_k|x_k^i)/\left({\displaystyle \sum _{i=1}^{N_s}{\omega }_{k-1}^{i}p(z_k|x_k^i)}\right).$$

(7)

The major disadvantage of Equation (7) is that the current measurement is not considered in the importance function to update the weight ${\omega }_k^i$. Therefore, it is necessary to incorporate the current measurement to the importance function.

2.2. Unscented Transform

Unscented transform [25] aims to find limited deterministic sigma points to estimate the mean and variance of a transformed variable $z=f(x)$. Here, $x$ is a multivariate Gaussian distribution with n-dimensional mean vector $m$ and $n\times n$ covariance matrix $P$, and $f(\cdot )$ is a linear or nonlinear function. The major difference between the particle filtering and the unscented transform is that an amount of random particles are used in particle filtering, whereas 2n + 1 deterministic sigma points are used in unscented transform. Simply speaking, the unscented transform can sufficiently provide a rough solution for solving any linear or nonlinear state space models. The working principle of the unscented transform is introduced as follows [25]:

Step 1. Determine 2n + 1 sigma points:

$$\begin{array}{l}{\chi }^0=m,\\ {\chi }^i=m+\sqrt{n+\lambda }{\left[\sqrt{P}\right]}_i,\\ {\chi }^{i+n}=m-\sqrt{n+\lambda }{\left[\sqrt{P}\right]}_i,i=1,\dots ,n,\end{array}$$

(8)

where ${[\cdot ]}_i$ takes the i-th column of the matrix. $\lambda =3{\alpha }^2-n$ is the parameter to control the spread of the sigma points around the mean. Moreover, $\sqrt{P}{\sqrt{P}}^T=P$.

Step 2. Propagate the determined sigma points via a non-linear function:

$$Y^i=f({\chi }^i),i=0,\dots ,2n,$$

(9)

where $Y^i$ is called the transformed sigma points.

Step 3. Estimate the mean and covariance of the transformed variable by:

$$\mathrm{E}[g(x)]=u_{\mathrm{U}}={\displaystyle \sum _{i=0}^{2n}{W}_i^m{Y}^i},$$

(10)

$$\mathrm{Cov}[g(x)]=S_{\mathrm{U}}={\displaystyle \sum _{i=0}^{2n}{W}_i^c(Y^i-u_{\mathrm{U}}){(Y^i-u_{\mathrm{U}})}^{\mathrm{T}}},$$

(11)

where $W_i^m$ and $W_i^c$ are respectively defined as:

$$\begin{array}{l}{W}_0^m=\frac{\lambda }{\lambda +n},\\ W_0^c=\frac{\lambda }{\lambda +n}+(1-{\alpha }^2+\beta ),\\ W_i^m=\frac{1}{2\left(\lambda +n\right)},i=1,\dots ,2n,\\ W_i^c=\frac{1}{2\left(\lambda +n\right)},i=1,\dots ,2n,\end{array}$$

(12)

where $\beta$ controls the prior information of $x$. According to [25], $\alpha$ and $\beta$ are recommended to be the default values 1 and 0, respectively.

3. State Space Formulation of Nonlinear Vibration Responses for the Bearing Prognostics of a Dynamic Rotor-Bearing System

In this section, state space formulation of nonlinear vibration responses collected from a dynamic rotor-bearing system is built to predict bearing RUL. A flowchart of the proposed bearing prognostic method is given in Figure 1. The procedure of the proposed bearing prognostic method is summarized as follows. Firstly, analyses of the nonlinear vibration responses should be conducted to extract a bearing degradation trend. This means that a BHI should be found to describe the current health condition of the bearing. Secondly, a state space model of the BHI is constructed to mathematically track the health evolutions of the BHI. Finally, bearing RUL is predicted using unscented particle filtering. The details of the above steps are introduced in the following sections.

3.1. Bearing Performance Degradation Assessment

Prior to bearing performance degradation assessment, a new bearing acceleration life testing is introduced. The aim of the acceleration life testing is to collect natural bearing degradation data. The experiments of bearing prognosis were designed by the first author. The experiment setup, shown in Figure 2a, is a rotary machine testing platform called “Machinery Fault Simulators (MFS)” that is manufactured by SpectraQuest Inc. (Richmond, VA, USA). The MFS is designed for conducting experimental work to reveal the signatures of common machinery faults without compromising production schedule or profits. This bench-top system has a larger baseplate, a powerful driver as a motor, various gearboxes, bearings, belt drives, reciprocating mechanisms, induction motors, pumps, and compressors. Then, various faults can be introduced either individually or jointly in a totally controlled environment. The types of faults can be rotor balancing, alignment, resonance, bearing defects, crack shafts, fan and mechanical rub, pump defects, etc.

Four accelerometers were attached to the dynamic rotor-bearing system. The sensitivities of the accelerometers from 1 to 4 are 0.093 V/g, 0.095 V/g, 0.093 V/g, and 0.092 V/g, respectively. Accelerometers 1 and 2 were used to measure nonlinear vibration responses from the vertical and horizontal directions of the right bearing. Accelerometers 3 and 4 were used to measure nonlinear vibration responses from the vertical and horizontal directions of the left bearing. Specifications of the right and left bearings and associated bearing fault characteristic frequencies are listed in

Table 1. Specifications of bearings used in our experiment.

Specifications of Bearings	Parameters
Bearing model	MB ER-12K
Number of rolling elements	8
Rolling element diameter	7.9375 mm
Pitch diameter	33.4772 mm
Contact angle	0°
Fundamental train frequency (FTF)	11.3 Hz
Ball pass frequency outer (BPFO)	91.4 Hz
Ball pass frequency inner (BPFI)	148.5 Hz
Ball spin frequency (BSF)	59.8 Hz

. During the experiment, the running speed of the motor or the shaft rotation frequency was set to 30 ± 0.5 Hz. The sampling frequency was set to 36,864 Hz, or collecting 36,864 vibration samples per second. The length of each measurement collected from each of the four accelerometers was set to 1.0 s. Measurements were collected every 30 min until the bearings generated loud sound and failed. A total of 145 vibration measurements was collected from each accelerometer or a total of 5,345,280 vibration samples was collected from each accelerometer. For the four accelerometers, the number of vibration samples collected was 21,381,120 samples. Since the correlation between vibration and rotor position can be affected by disturbances in the spindle speed [31,32], in order to obtain a steady operating condition, the motor rotating speed or the shaft rotation frequency was set to 30 ± 0.5 Hz.

To accelerate bearing life from a normal health condition to failure, a hydraulic jack with a continuously applied force was exerted to a location of the outer race of the right bearing through an adapter. When the hydraulic jack was working or applied to the surface of the bearing’s outer race, a force of 128 kg (20 kg/cm²) was applied to the bottom of that bearing’s outer race. That is, the hydraulic jack acted as an external force continuously exerted on the bearing so as to generate a wear status to the bearing. The schematic working principle of the hydraulic jack with the adapter is shown in Figure 2b–d. Moreover, the right bearing was running without proper lubricant. The experiment stopped until the rotor-bearing system had heavily vibrated and created a loud noise. There were 145 timely measurements collected until the shutdown of the rotor-bearing system. In this paper, the vibration data collected from accelerometer 1 were used as this bearing suffered from extra external force generated by the hydraulic jack and continuous wear in that location of the bearing was expected.

It is known that impacts generated by bearing defects excite resonant frequencies of structures and adjacent components, and they induce modulation phenomenon. For bearing fault data produced by a single artificial defect, it is more likely to use some advanced signal processing methods to extract bearing fault features, such as bearing fault characteristic frequencies, as shown in Table 1. For natural bearing degradation data, bearing fault types are unknown and multiple bearing faults may occur simutaneously. In our preliminary analyses, some advanced signal processing methods, such as the fast Kurtogram and the improved Kurtogram, fail to provide any information related to bearing faults. Considering these points, it is preferable to use a high-pass filter to retain some resonant frequency bands and to remove the interruptions from low-frequency vibration components. A vibration measurement collected from the vertical direction of the casing of the right bearing at measurement number 50 and its corresponding frequency spectrum are plotted in Figure 3a,b, respectively. Note that a high-pass filter with a cutoff frequency of 6000 Hz to 18,000 Hz was used to retain the high-frequency components generated by the tested bearing. That is, the high-frequency range from 6000 Hz to 18,000 Hz, which contains high vibration energy, was retained and analyzed.

Since the statistical root mean square (RMS) can characterize the raw degradation data in the view of vibration energy, the RMS is used as an important fault indicator. If faults have occurred in a bearing, the generated vibration energy is expected to be much higher than that from a normal running bearing. The measured RMS was preprocessed by a high-pass filter and its result is shown in Figure 4a, in which it is shown that the RMS fluctuates so much that it is not easy to use for bearing prognostics. Consequently, it is necessary to construct a BHI that has a direct connection with the RMS. In terms of the central tendency of the RMS, a BHI is defined as follows:

$$\mathrm{BHI}(j)=\frac{{\displaystyle \sum _{k=1}^j\mathrm{RMS}(k)}}{j},j=1,2,\dots ,145,$$

(13)

where RMS(k) is the i-th RMS after the high-pass filter was conducted on a vibration signal stored at the i-th measurement number. The BHI is plotted in Figure 4b, in which it can be seen that the BHI increases continuously and sharply until it reached measurement number 100. After that, the increase rate of the BHI was slower than before. This is because, when the tested bearing was suffering from the continuous high force exerted by the hydraulic jack, wear was created [31,32].

Initially only small spall or cracks would be formed, but the damage to the surface became severe over time. Continuously applying such force to the same location would create wear to the bearing at an outer race, balls and then an inner race and a cage. That is, all parts at this point of the bearing would be damaged in time. Moreover, once the balls were damaged, they would cause other parts of the bearing to be damaged when the damaged balls rotated to the other part of the bearing. At a later stage, the outer race, the inner race, the balls, and the cage were all worn at different degrees. The conditions would deteriorate over time with a continuous increase in vibration energy, as shown in Figure 4b. At a certain point, because of the continuous rubbing of different surfaces of the bearing, some of the damaged surfaces became smoother than before, causing the rate of increasing in vibration energy to be slower or the slope of degradation trend becomes flat. In Figure 4b, it is seen that the proposed BHI can monitor the current health conditions of the right bearing.

3.2. State Space Formulation of Bearing Performance Degradation

To predict the future health conditions of the right bearing, its state space model should be built. Firstly, it is necessary to find a mathematical model to properly describe the BHI. In Figure 4b, it seems that measurement number 116 is the peak of the trend created by BHI. Such a degradation trend is the typical bearing wear trend.

Because the maximum of the BHI is 0.1537 at measurement number 116, the failure threshold of the BHI could be set to 0.1537. Moreover, the effective bearing degradation length is from 1 to 116. By using goodness of fit, the R-squared values of a linear polynomial function, a quadratic polynomial function, and an exponential function are 0.9866, 0.9925, and 0.9663, respectively. The closer the value of the R-squared is to 1, the better the performance of the function. Consequently, the quadratic polynomial function is used as the measurement function of the bearing state space model. Then, the bearing state space model used in this paper is constructed as follows:

$$a_k=a_{k-1}+v_1$$

(14)

$$b_k=b_{k-1}+v_2$$

(15)

$$c_k=c_{k-1}+v_3$$

(16)

$${\mathrm{BHI}}_k=f(k)+v_4=a_k\times k^2+b_k\times k+c_k+v_4,$$

(17)

where $a_k$, $b_k$, and $c_k$ are the parameters of the BHI. $v_1$, $v_2$, $v_3$, and $v_4$ are the additive Gaussian noises with zero means and different standard deviations ${\sigma }_1,{\sigma }_2,{\sigma }_3$, and ${\sigma }_4$, respectively. Initializations of Equations (14) to (17) are clarified at the beginning of Section 4. At last, ${\sigma }_1^2,{\sigma }_2^2,{\sigma }_3^2$ form a state noise covariation matrix $\mathit{Q}$, in which the diagonal elements are ${\sigma }_1^2,{\sigma }_2^2,{\sigma }_3^2$, respectively.

3.3. Posterior State Parameter Estimation of the Bearing State Space Model Using Unscented Particle Filtering

According to the bearing state space model and the fundamental theories of particle filtering and unscented transform, when a new BHI is available, the three state parameters $a_k$, $b_k$, and $c_k$ can be iteratively updated by the following steps.

Step 1. Draw initial random particles ${\left\{a_0^i\right\}}_{i=1}^{N_s}$, ${\left\{b_0^i\right\}}_{i=1}^{N_s}$, and ${\left\{c_0^i\right\}}_{i=1}^{N_s}$ from the initial distributions of $q(a_0)=N(a_0,{\sigma }_1^2)$, $q(b_0)=N(b_0,{\sigma }_2^2)$, and $q(c_0)=N(c_0,{\sigma }_3^2)$, respectively. If $N_s$ is taken to be large then we will be estimating the true likelihood quite precisely, but its computing will be very expensive. On the other hand, a small value of $N_s$ will result in cheap evaluations, but low computation burden. Here, $N_s$ is equal to 1000, which is manually selected and is deemed to be sufficient for our research’s purpose. Set all initial weights ${\left\{{\omega }_0^i\right\}}_{i=1}^{N_s}$ to a value of $1/N_s$. Suppose that the initial distribution for the use of the unscented transform is $p(x_0^i)=p(a_0^i,b_0^i,c_0^i)$ ≈ $\mathrm{N}(a_0^i,b_0^i,c_0^i|m_0,P_0)$. Initializations of the $m_0$ and the $P_0$ are clarified in Section 4. According to the theory of the unscented transform, only seven sigma points ($2\times 3+1$) are required to propagate via the measurement function:

$$\begin{array}{l}{\chi }_0^0=m_0,\\ {\chi }_0^j=m_0+\sqrt{3}{\left[\sqrt{P_0+\mathit{Q}}\right]}_j,\\ {\chi }_0^{j+2}=m_0-\sqrt{3}{\left[\sqrt{P_0+\mathit{Q}}\right]}_j,j=1,2,3\end{array}$$

(18)

$$Y_0^j=f({\chi }_0^j),j=0,\dots ,6.$$

(19)

The predicted mean $u_{\mathrm{U}}^0$, the predicted covariance $S_{\mathrm{U}}^0$, and the cross-covariance $C_{\mathrm{U}}^0$ of the state and the measurement are calculated as follows:

$$\mathrm{E}[g(x_0)]=u_{\mathrm{U}}^0={\displaystyle \sum _{j=0}^6{W}_j^m{Y}_0^j}$$

(20)

$$\mathrm{Cov}[g(x_0)]=S_{\mathrm{U}}^0={\displaystyle \sum _{j=0}^6{W}_j^c(Y_0^j-u_{\mathrm{U}}^0){(Y_0^j-u_{\mathrm{U}}^0)}^{\mathrm{T}}}+v_4^2$$

(21)

$$\mathrm{Cov}[x_0,g(x_0)]=C_{\mathrm{U}}^0={\displaystyle \sum _{j=0}^6{W}_j^c({\chi }_0^j-m_0){(Y_0^j-u_{\mathrm{U}}^0)}^{\mathrm{T}}}.$$

(22)

When the first BHI is available, the $p(a_1^i,b_1^i,c_1^i|{\mathrm{BHI}}_1)$ is approximated by the multivariate Gaussian distribution with the following parameters:

$$m_1=m_0+C_{\mathrm{U}}^0{\left(S_{\mathrm{U}}^0\right)}^{-1}\left({\mathrm{BHI}}_1-u_{\mathrm{U}}^0\right)$$

(23)

$$P_1=P_0-C_{\mathrm{U}}^0{\left(S_{\mathrm{U}}^0\right)}^{-1}{S}_{\mathrm{U}}^0{\left(C_{\mathrm{U}}^0{\left(S_{\mathrm{U}}^0\right)}^{-1}\right)}^{\mathrm{T}},$$

(24)

where $C_{\mathrm{U}}^0{\left(S_{\mathrm{U}}^0\right)}^{-1}$ is the filter gain and can be determined by Equations (21) and (22).

The importance function $q(x_k^i|x_{k-1}^i,z_k)$ used in Equation (6) is replaced with the joint posterior probability density function $p(a_1,b_1,c_1|{\mathrm{BHI}}_1)$. Then, the weight updating formula is rewritten as:

$${\omega }_1^i=\frac{{\omega }_0^i\frac{p(z_1|a_1^i,b_1^i,c_1^i)p(a_1^i|a_0^i)p(b_1^i|b_0^i)p(c_1^i|c_0^i)}{p(a_1^i,b_1^i,c_1^i|{\mathrm{BHI}}_1)}}{{\displaystyle \sum _{i=1}^{N_s}{\omega }_0^i\frac{p(z_1|a_1^i,b_1^i,c_1^i)p(a_1^i|a_0^i)p(b_1^i|b_0^i)p(c_1^i|c_0^i)}{p(a_1^i,b_1^i,c_1^i|{\mathrm{BHI}}_1)}}}.$$

(25)

Compared with the weight updating shown in Equation (7), Equation (25) considers the first BHI when the importance function is used. Then, the posterior probability density functions of the three parameters are respectively represented by:

$$p(a_1|{\mathrm{BHI}}_1)\approx {\displaystyle \sum _{i=1}^{N_s}{\omega }_1^i\delta (a_1-a_0^i)}$$

(26)

$$p(b_1|{\mathrm{BHI}}_1)\approx {\displaystyle \sum _{i=1}^{N_s}{\omega }_1^i\delta (b_1-b_0^i)}$$

(27)

$$p(c_1|{\mathrm{BHI}}_1)\approx {\displaystyle \sum _{i=1}^{N_s}{\omega }_1^i\delta (c_1-c_0^i)}.$$

(28)

To avoid the degeneracy problem of the original particle filtering, which is that most of the weights become negligible after a few iterations, the systematic resampling is required to redrawn random particles from the posterior probability density functions calculated by the original particle filtering according the size of ${\omega }_1^i$ if the ${\left({\displaystyle \sum _{i=1}^{N_s}{\left({\omega }_1^i\right)}^2}\right)}^{-1}$ is smaller than half of the N_s.

Step 2. Draw $N_s$ random particles ${\left\{a_{k-1}^i\right\}}_{i=1}^{N_s}$, ${\left\{b_{k-1}^i\right\}}_{i=1}^{N_s}$, and ${\left\{c_{k-1}^i\right\}}_{i=1}^{N_s}$ from Equations (26)–(28), respectively. Suppose that the posterior probability density function $p(x_{k-1}^i|{\mathrm{BHI}}_{1:k-1})=p(a_{k-1}^i,b_{k-1}^i,c_{k-1}^i|{\mathrm{BHI}}_{1:k-1})$ ≈ $\mathrm{N}(a_{k-1}^i,b_{k-1}^i,c_{k-1}^i|m_{k-1},P_{k-1})$ at iteration k − 1 is estimated by using the unscented transform. Only seven sigma points are required to propagate via Equation (17):

$$\begin{array}{l}{\chi }_{k-1}^0=m_{k-1},\\ {\chi }_{k-1}^j=m_{k-1}+\sqrt{3}{\left[\sqrt{P_{k-1}+\mathit{Q}}\right]}_j,\\ {\chi }_{k-1}^{j+2}=m_{k-1}-\sqrt{3}{\left[\sqrt{P_{k-1}+\mathit{Q}}\right]}_j,j=1,2,3\end{array}$$

(29)

$$Y_{k-1}^j=f({\chi }_{k-1}^j),j=0,\dots ,6.$$

(30)

The predicted mean $u_{\mathrm{U}}^{k-1}$, the predicted covariance $S_{\mathrm{U}}^{k-1}$, and the cross-covariance $C_{\mathrm{U}}^{k-1}$ of the state and the measurement are calculated as follows:

$$\mathrm{E}[g(x_{k-1})]=u_{\mathrm{U}}^{k-1}={\displaystyle \sum _{j=0}^6{W}_j^m{Y}_{k-1}^j}$$

(31)

$$\mathrm{Cov}[g(x_{k-1})]=S_{\mathrm{U}}^{k-1}={\displaystyle \sum _{j=0}^6{W}_j^c(Y_{k-1}^j-u_{\mathrm{U}}^{k-1}){(Y_{k-1}^j-u_{\mathrm{U}}^{k-1})}^{\mathrm{T}}}+v_4^2$$

(32)

$$\mathrm{Cov}[x_{k-1},g(x_{k-1})]=C_{\mathrm{U}}^{k-1}={\displaystyle \sum _{i=0}^6{W}_i^c({\chi }_{k-1}^i-m_{k-1}){(Y_{k-1}^i-u_{\mathrm{U}}^{k-1})}^{\mathrm{T}}}.$$

(33)

When a new BHI is available, the $p(a_k^i,b_k^i,c_k^i|{\mathrm{BHI}}_{1:k})$ is approximated by a multivariate Gaussian distribution with the following parameters:

$$m_k=m_{k-1}+C_{\mathrm{U}}^{k-1}{\left(S_{\mathrm{U}}^{k-1}\right)}^{-1}\left(BH{I}_k-u_{\mathrm{U}}^{k-1}\right)$$

(34)

$$P_k=P_{k-1}-C_{\mathrm{U}}^{k-1}{\left(S_{\mathrm{U}}^{k-1}\right)}^{-1}{S}_{\mathrm{U}}^{k-1}{\left(C_{\mathrm{U}}^{k-1}{\left(S_{\mathrm{U}}^{k-1}\right)}^{-1}\right)}^{\mathrm{T}}.$$

(35)

The weight updating formula is rewritten as:

$${\omega }_k^i=\frac{{\omega }_{k-1}^i\frac{p(z_k|a_k^i,b_k^i,c_k^i)p(a_k^i|a_{k-1}^i)p(b_k^i|b_{k-1}^i)p(c_k^i|c_{k-1}^i)}{p(a_k^i,b_k^i,c_k^i|{\mathrm{BHI}}_{1:k})}}{{\displaystyle \sum _{i=1}^{N_s}{\omega }_{k-1}^i\frac{p(z_k|a_k^i,b_k^i,c_k^i)p(a_k^i|a_{k-1}^i)p(b_k^i|b_{k-1}^i)p(c_k^i|c_{k-1}^i)}{p(a_k^i,b_k^i,c_k^i|{\mathrm{BHI}}_{1:k})}}}.$$

(36)

If ${\left({\displaystyle \sum _{i=1}^{N_s}{\left({\omega }_k^i\right)}^2}\right)}^{-1}$ is smaller than half of N_s, the systematic resampling is conducted to redrawn random particles from the posterior probability density function according to the size of ${\omega }_k^i$.

Step 3. Increase k to k + 1 and repeat Step 2 until $k>M$ is satisfied. Here, M is the length of the available BHI. Then, the posterior probability density functions of the three parameters at iteration M can be respectively represented by:

$$p(a_M|{\mathrm{BHI}}_{1:M})\approx {\displaystyle \sum _{i=1}^{N_s}{\omega }_M^i\delta (a_M-a_M^i)}$$

(37)

$$p(b_M|{\mathrm{BHI}}_{1:M})\approx {\displaystyle \sum _{i=1}^{N_s}{\omega }_M^i\delta (b_M-b_M^i)}$$

(38)

$$p(c_M|{\mathrm{BHI}}_{1:M})\approx {\displaystyle \sum _{i=1}^{N_s}{\omega }_M^i\delta (c_M-c_M^i)}.$$

(39)

3.4. Bearing Remaining Useful Life Prediction

Once the posterior probability density functions of the three parameters are estimated, extrapolations of the measurement function to a specified failure threshold, such as 0.1537 used in this paper, are used to calculate bearing RUL at measurement number M:

$$\mathrm{RUL}(M)={\displaystyle \sum _{i=1}^{N_s}{\omega }_M^i\delta (\mathrm{RUL}-\inf (kk\in \mathrm{int}:a_M^i\times k{k}^2+b_M^i\times kk+c_M^i\ge 0.1537)+kk)},i=1,2,\dots ,N_s$$

(40)

where the “int” is the set of all integers and the “inf” takes the greatest lower bound of a set and kk is a future possible measurement number. The 50th, 5th, and 95th percentiles of Equation (40) can be regarded as the estimated RUL and its lower and upper bounds.

4. A Case Study of Bearing Prognostics

In this section, nonlinear vibration responses collected from the bearing acceleration life testing are used to illustrate how the proposed bearing prognostic method works and to validate its effectiveness in predicting bearing RUL. If more historical BHI data are available, the initial estimates of $a$, $b$, and $c$ were calculated using the nonlinear least squares regression on historical BHI data. If historical BHI data are not available, the nonlinear least squares regression on part of the available BHI can be conducted. In this paper, based on the BHI from 1 to 20, the initial estimates of $a$, $b$, and $c$ were calculated as −1.81 × 10⁻⁶, 9.11 × 10⁻⁴, and 7.635 × 10⁻², respectively. The 95% confidence bounds of the initial estimates of $a$, $b$, and $c$ were calculated as (−2.189 × 10⁻⁶, −1.431 × 10⁻⁶), (8.654 × 10⁻⁴, 9.569 × 10⁻⁴) and (7.519 × 10⁻², 7.751 × 10⁻²), respectively. Therefore, according to six sigma principle, the standard deviations of the three parameters were artificially set to 1.2633 × 10⁻⁷, 1.5250 × 10⁻⁵, and 3.8333 × 10⁻⁴, respectively. Because the scale of the BHI is small, the standard deviation of the measurement noise was set to 0.005, which results in a small residual in Equation (17).

At measurement number 30, the predictions of the RUL are plotted in Figure 5. In Figure 5a, it is observed that the predictive values obtained by the unscented particle filtering are well matched with the true BHI. This means that the predictions obtained by the unscented particle filtering are suitable for tracking the BHI. The extrapolations of the measurement function are highlighted by the green dashed lines, which indicate the future health evolutions of the BHI. The probability density function of the RUL is plotted in Figure 5b, where the 5th, 50th, and 95th percentiles of the RUL are 68, 78, and 86 measurement numbers, respectively. The error between the estimated RUL 78 and the actual RUL 86 is eight measurement numbers.

At measurement number 50, the predictions obtained by using the proposed bearing prognostic method are presented in Figure 6. In Figure 6a, it is observed that the predictive BHI obtained by the unscented particle filtering is very close to the true BHI. In other words, unscented particle filtering results in high prediction accuracies for the BHI. In Figure 6b, the probability density function of the RUL is plotted. It is found that the 5th, 50th, and 95th percentiles of the RUL are 56, 64, and 95 measurement numbers, respectively. The error between the actual RUL 66 and the predicted RUL 64 is two measurement numbers.

To reduce the length of this paper, more bearing RUL predictions at document numbers from 20 to 110 with an increment of 10 are tabulated in

Table 2. The predictions of the RUL by the proposed bearing prognostic method (Unit: measurement number).

Prediction at File Number	5th Percentile of Predicted RUL	50th Percentile of Predicted RUL	95th Percentile of Predicted RUL	Actual RUL	Error between Actual RUL and 50th Percentile of Predicted RUL
20	67	73	83	96	23
30	68	78	95.5	86	8
40	58	69	90	76	7
50	56	64	95	66	2
60	38	47	71	56	9
70	25	28	33	46	8
80	29	30	32	36	6
90	19	24	37	26	2
100	10	14	20	16	2
110	2	3	4	6	3

, where the proposed bearing prognostic method results in the high prediction accuracies of the RUL.

In the standard particle filtering, the importance function is chosen as a prior state transition function and the calculation of the weights is largely simplified. In other words, the recent BHI is not considered in the importance function. To highlight the advantage of the unscented particle filtering, some comparisons are made by replacing the unscented particle filtering with the standard particle filtering in our proposed prognostic method. The results obtained by the standard particle filtering are tabulated in

Table 3. The predictions of the RUL by replacing the unscented particle filtering with the standard particle filtering in our proposed prognostic method (Unit: measurement number).

Prediction at File Number	5th Percentile of Predicted RUL	50th Percentile of Predicted RUL	95th Percentile of Predicted RUL	Actual RUL	Error between Actual RUL and 50th Percentile of Predicted RUL
20	66	73	85	96	23
30	65	76	100	86	10
40	51	61	82	76	15
50	52	63	89	66	3
60	37	46	69	56	10
70	31	37	58	46	9
80	19	25	35	36	11
90	11	15	22	26	11
100	5	8	16	16	8
110	2	3	7	6	3

. From the prediction errors shown in Table 2 and Table 3, it is found that the prediction errors obtained by the standard particle filtering are larger than the prediction errors obtained by the unscented particle filtering, especially after file number 20. It should be noted that the unscented transform is not considered for comparisons because the unscented transform is just a method that can be used to make fast and approximate parameter distributions using a few deterministic sigma points, while the particle filtering uses a number of random particles to infer parameter distributions and is able to provide higher prediction accuracy.

5. Conclusions

The state space formulation of the nonlinear vibration responses collected from the dynamic rotor-bearing system was reported in this paper for the extensions from bearing diagnostics to prognostics. Prior to prediction of the RUL of bearing, the bearing performance degradation assessment was conducted by constructing the BHI. This indicator was then used to monitor the current health conditions of the bearing. Then, to mathematically describe the future health evolutions of BHI, the state space model of the BHI was constructed. Unscented particle filtering was introduced to calculate the posterior probability density functions of the parameters of the state space model. The extrapolations of the measurement function to the specified failure threshold were used to estimate the RUL of the inspected bearing. To validate the effectiveness of the proposed methods, a new setup for monitoring the bearing operational life was designed. Vibration data related to the deterioration of the bearing were collected naturally from real running bearings. The predicted results of the RUL demonstrated that the bearing prognostic method reported here is effective at inferring future bearing health conditions and the RUL of bearing. Moreover, comparisons with the conventional particle filtering were made to demonstrate that the method of unscented particle filtering can achieve better prediction accuracy than conventional particle filtering methods.