Application of EMD-WVD and particle filter for gearbox fault feature extraction and remaining useful life prediction

2.1. Short overview of EMD

The EMD is a widely-used time-frequency analysis method in the signal processing field, which mainly includes two steps. Firstly, the EMD is used to decompose nonlinear and non-Gaussian time series into a number of IMFs and a residual. Second, each IMF is Hilbert transformed to calculate the time-frequency distribution. The first step is the most important. It makes the original signal linear and stationary time series, as well as produces a series of data sequences with different characteristics, that each sequence is an IMF. The EMD method based on such assumptions: (1) the original signal has at least two extreme points, which are a maximum and a minimum, (2) the characteristic time scale is defined as a time interval between the local extreme value points, (3) if the original signal does not have an extreme point but there is an inflection point, the extreme point can be obtained by differential calculus once or more times, then the decomposed results will be get by integral calculus. Therefore, we can see that EMD is a self-adaptive signal processing method, which especially suits for the nonlinear and non-Gaussian signals.

The EMD algorithm exists in two stop conditions: one is aiming at the IMF, the repeating subtracting process will stop until the IMF satisfied the condition which is called component termination condition, the other one is the EMD decomposition termination condition for the original signal, which is namely decomposition termination condition. The EMD decomposition is terminated until the residual is less than the predetermined value or becomes a monotonic function, the EMD decomposition process is over.

The steps of EMD decomposition are shown as follows.

(1) Identify all the local minima and local maxima of the given signal $x\left(t\right)$.

(2) According to the local extrema of the signal, connect the local maxima and minima to construct the upper and lower envelopes $u\left(t\right)$ and $v\left(t\right)$ by a cubic spline interpolation.

(3) Calculate the mean of signal: $m\left(t\right)=\left(u\left(t\right)+v\left(t\right)\right)/2$.

(4) Obtain the new series: $h\left(t\right)=x\left(t\right)-m\left(t\right)$.

(5) If $h\left(t\right)$ satisfies the two conditions, take it as the first IMF. Otherwise take $h\left(t\right)$ as the original signal and the above processes are repeated until it satisfies the conditions of IMF.

(6) The residual $r\left(t\right)$ is got after separating $c_1\left(t\right)$ from the $x\left(t\right)$, until the final residue becomes a constant value, or $x\left(t\right)$ is replaced by $r\left(t\right)$ and the above process is repeated.

The signal $x\left(t\right)$ can be decomposed into $n$ IMFs $c_i\left(t\right)$ and a residual $r_n\left(t\right)$ via EMD. The original signal $x\left(t\right)$ can be finally expressed as:

In the process of EMD decomposition, components of different time scales can be found in the same one of the IMFs when decompose a complex signal directly, which is called as mode mixing. In simple terms, there is a rapid change time scale in the signal which makes the different mode composition not to be able being separated effectively according to the characteristic time scales. Then the same IMF contains multiple mode composition inside, so the original signal can’t be clearly reflected the intrinsic nature. And the reasons for this phenomenon has the relationship with the EMD algorithm, signal sampling rate and signal change in a certain time scale.

2.2. Short overview of WVD

Time-frequency distributions have been used extensively for nonlinear and non-Gaussian signal analysis. The Wigner-Ville distribution is the best known one for a time-frequency analysis. The Wigner-Ville distribution can be viewed as the signal energy time-frequency joint distribution, and is an important tool for analyzing nonlinear and non-Gaussian signal. The auto-Wigner-Ville distribution of continuous signal $x\left(t\right)$ is defined as:

Accordingly, assumed another continuous signal $y\left(t\right)$, the cross-Wigner-Ville distribution of it is defined as:

where $\omega$ represents frequency, $\tau$ represents time variables. $x^{*}$ and $y^{*}$ are the complex conjugate of continuous signal $x\left(t\right)$ and $y\left(t\right)$.

The disadvantages of Wigner-Ville distribution are the cross-term interference, which is caused by the time and frequency interference caused by nonlinear transformation. When the original signal contains multiple components, there will be no physical meaning of oscillation component at the midpoint of the time-frequency center coordinates between each of the two components in the Wigner-Ville distribution, which provides a false energy distribution and is affected by the physical interpretation of Wigner-Ville distribution. Take the signal $x\left(t\right)=x_1\left(t\right)+x_2\left(t\right)$ as an example, the Wigner-Ville distribution is defined:

According to the definition of WVD, $WV{D}_{x1,x2}\left(t,\omega \right)=WV{D}_{x2,x1}*\left(t,\omega \right)$.

$WV{D}_{x1,x2}\left(t,\omega \right)$ and $WV{D}_{x2,x1}\left(t,\omega \right)$ is the conjugate, then:

This suggests that the WVD of two signals summed is not the sum of their each WVD, in addition to two auto-WVDs, there is a cross-WVD which is named as a cross-term interference. Because the cross-term interference is usually an oscillation, which amplitude can be reached twice as much the auto-WVD, which caused the time-frequency signal characteristics is blur.

2.3. The integration of EMD and WVD for gearbox fault feature extraction

The Wigner-Ville distribution can deal with nonlinear and non-Gaussian signals to get the time–frequency distribution, which is a powerful tool for signal processing. In order to reduce the effects of its cross-term interference, the WVD and EMD are combined together, as well as, this method can also reduce the mode mixing of EMD. Firstly, the original signal is decomposed into IMFs which frequency components are relatively single via EMD, then some IMFs that are selected according to the specific rules are processed by WVD. As we know, the original vibration signal can be decomposed into an IMF component by the EMD, and the IMF component signal frequency is single, so the WVD of the IMFs can restrain the frequency interference. Meanwhile, the Wigner-Ville distribution makes the different mode composition of IMF separated, resulting the mode mixing restrained effectively. Finally, the method can be used to improve the effectiveness of vibration signal feature extraction.

Based on this, a fault feature extraction method based on EMD-WVD is proposed in this paper to realize the nonlinear and non-Gaussian signal analysis, the steps are as follows:

(1) The original signal is decomposed via the EMD to get a set of IMFs which contains the different frequencies.

(2) Selecting the appropriate IMFs, they will be analyzed by WVD to get the signal $C_{c_i}$, $i=$ 1, 2,….

(3) The $C_{c_i}$ will be beaded together to get the Wigner-Ville distribution of signal, which is shown as $C_x=\sum _{i=1}^n{C}_{c_i}$.

During the signal processing, the WVD analysis was conducted for each of the IMFs, so in the final results the restraint for cross-term interference caused by the EMD method will still retain. Such the Wigner-Ville distribution of the superposition original signal can effectively characterize the characteristics in time and frequency domain.

The IMFs which are analyzed by WVD are selected according to kurtosis in this paper. The kurtosis is the degree of a signal waveform peak parameters, which is defined as:

where $\sigma$ represents the standard deviation of signal $x$, $\mu$ represents the mean of signal $x$.

When the gearbox is in the normal state, the distribution of its monitored vibration signal can approximate to the normal one; at this moment, its kurtosis is approximately equal to 3. When the gearbox degenerates to a failure, the monitored vibration signal is no longer than the normal distribution because of an impact component of gearbox, which can lead to a bigger kurtosis. On the other hand, the degradation information of gearbox is often contained in the modulation signal of amplitude which with a more impact component. The vibration caused by the state degradation will cause the gearbox shaking at a different frequency, and each of the IMFs contains different natural vibration frequency caused by degradation. Therefore, the IMF is selected to get the fault feature of gearbox vibration signal according to the kurtosis criterion in this paper.

3.1. Bayesian estimation for particle filter

Particle filter is a hybrid prediction method combined with date-driven and model-driven two kinds of methods. The particle filter employs Monte Carlo simulation of a state dynamic model and Bayesian estimation for estimating the posterior probability density function (PDF) of the state of a degrading component at future times. Therefore, the method is an effective tool to predict the remaining life of nonlinear and non-Gaussian systems.

The adaptive particle filtered Bayesian inference probability index is established for gearbox remaining life prediction in this paper. The state equation and measurement equation of particle filter are expressed as follows:

where $x_k$ is the hidden state to be estimated at time $k$, $y_k$ is the observation at time $k$, $f$ and $h$ represent the known process and observation functions respectively, $\left\{n_{k-1},k\in N\right\}$ and $\left\{{\omega }_k,k\in N\right\}$ are the dimensions of the state and process noise.

When the initial PDF of the state is $p\left(x_{k-1}|z_{0:k-1}\right)$, in order to obtain the posterior probability density function $p\left(x_k|z_{0:k}\right)$ of gearbox RUL, the state prediction is defined as:

The new observations $z_k$ are collected at time $k$. According to the rule of Bayesian, the posterior probability density for system state is updating over time, so posterior probability distribution of the current state $x_k$ is obtained:

where the constant is $p\left(z_k|z_{0:k-1}\right)=\int p\left(x_k|z_{0:k-1}\right)p\left(z_k|x_k\right)d{x}_k$.

Above is the prediction process and update process for Bayesian estimation of particle filter. However, the gearbox vibration signal is nonlinear and non-Gaussian, which makes it difficult to obtain the distribution. In this paper, we utilize Bayesian filtering and Monte Carlo algorithm to obtain the distribution.

Assuming that update monitoring information $z_{0:k}$ is known, the posterior probability density of system state $x_{0:k}$ can be expressed as follows:

where $\delta \left(\cdot \right)$ is the Dirac delta measure, $p\left(x_{0:k}|z_{0:k}\right)$ is the posterior probability function, and $p\left({\xi }_{0:k}|z_{0:k}\right)$ is the priori probability function. If the real posterior probability function $p\left(x_{0:k}|z_{0:k}\right)$ can be obtained, then Eq. (9) can be calculated according to Eq. (10):

where $x_{0:k}^i$, $i=$ 1, 2,..., $N_s$ are independent random samples sampling from $p\left(x_{0:k}|z_{0:k}\right)$.

$p\left(x_{0:k}|z_{0:k}\right)$ are multi-type and non-standard, and it is difficult to get the results. Then importance function sampling method is proposed, which means the probability distribution and $p\left(x_{0:k}|z_{0:k}\right)$ is the same, and the PDF distribution $\pi \left(x_{0:k}|z_{0:k}\right)$ is known. Therefore, Eq. (10) can be transformed as follows:

where $w_k^{*i}=p\left(z_{0:k}|x_{0:k}^i\right)p\left(x_{0:k}^i\right)/p\left(z_{0:k}\right)\pi \left(x_{0:k}^i|z_{0:k}\right)$ is the weight of the system states, $x_{0:k}^i$, $i=$ 1, 2,..., $N_s$, which can be calculated by $\pi \left(x_{0:k}^i|z_{0:k}\right)$ sampling. $p\left(z_{0:k}|x_{0:k}^i\right)$ are likelihood value of observation sequence. In order to obtain the weight, we assume that $p\left(z_{0:k}\right)=\int p\left(z_{0:k}|x_{0:k}\right)p\left(x_{0:k}\right)d{x}_{0:k}\text{,}$ and the probability distribution $p\left(x_{0:k}|z_{0:k}\right)$ can be expressed as follows:

where:

The state distribution at time $t_k$ can be estimated according to the system state distribution $p\left(x_{0:k-1}|z_{o:k-1}\right)$. According to the theory of Bayesian filtering, $p\left(x_{0:k}|z_{o:k}\right)$ can be expressed as follows:

The implementation of the prediction algorithm based on PF is shown in Fig. 1.

Fig. 1The flowchart for the implementation of PF

3.2. Simulated signal analysis compared with KF

In this section, we take the KF as an example and compare the state estimation performance with particle filter, which prove that PF is more effective under the condition of nonlinear and non-Gaussian. A simulation test proposed with a discrete time adaptive noise of high-order non-linear state space model, and the function is expressed as follows:

where $w$ and $v$ obey the standard normal distribution. Simulation repeats 50 times to get 50 real value and observation value of simulation state. The process is repeated twice, and the results are shown in Fig. 2.

Fig. 2KF and PF simulation comparison

As shown in Fig. 2, KF state estimation results is bad, and has obvious deviation with the real data. On the other hand, PF has a higher precision, which indicates that PF has a better prediction performance for the nonlinear and non-Gaussian state space.

4.1. The life cycle test of gearbox

The test object gearbox is shown in Fig. 3. The test bed is mainly composed of gearbox, speed and torque sensor (3056B4, produced by DYTRAN Company), magnetic powder electromotor (YCT180-4A), a load (FZJ-5) and other parts. Among them, the test bed is a steel platform (1500×600×700 mm) with a longitudinal $T$ slot; the output shaft axis of electromagnetic motor is parallel with the gearbox input shaft, where the power transfer is done through a coupling (6305). The gearbox output shaft is connected with the input shaft of the load through a coupling (6309). The vibration signal is acquired by the PXI-1031 case that produced by NI Company and PXI-4472B data acquisition card.

The most important part of the experiment is a secondary transmission gearbox, as shown in Fig. 4.

Fig. 3Image of test bed. (1 – load, 2, 3 – speed and torque sensor, 4 – electromotor, 5 – test bed, 6 – gearbox system)

Fig. 4Gearbox structure sketch map

The maximum capacity for lubricant is 2$L$, the rated speed is 1450 rpm, and rated transmission power is 0.75 kW. At particular parts of the gearbox there are four acceleration sensors numbered 1-4, the test load for gearbox is 2-2.5 times as much as the rated load. The signal is collected every five minutes with the sampling frequency for 20 kHz. The main parameters of the gearbox are shown in Table 1.

The life cycle test for gearbox is carried out; the vibration signal is acquired by four sensors. According to the monitoring vibration information, a failure of severe wear of the gear tooth surface occurs when the gearbox is working during about 400 hours, and the result is shown in Fig. 5. The incomplete data of worn-teeth is collected, which means that the wear condition is monitored only at some run time, the wear data obtained is shown in Table 2.

Fig. 5Comparison of gear after the experiment

a) Gear before experiment

b) Gear after experiment

4.2. Fault feature extraction for gearbox vibration signal

First of all, a vibration signal of gearbox is analyzed, as the original vibration signal of channel 1 at a different monitoring time is shown in Fig. 6.

Fig. 6Original vibration signal of channel 1 at different monitoring time

a)$t=$ 1 h

b)$t=$ 100 h

c)$t=$ 200 h

d)$t=$ 300 h

e)$t=$ 400 h

f)$t=$ 442 h

However, the health status of gearbox can’t be directly identified by the time domain analysis. Thus, fault feature extraction of vibration signal based on EMD-WVD is proposed in this paper.

The vibration signal of the gearbox shown in Fig. 6 at a different monitoring time is decomposed by the EMD, and the IMFs and a residual are obtained, as shown in Fig. 7. $c_1$-$c_7$ are the IMFs decomposed from the vibration signal of gearbox and $c_8$ is the residual.

The kurtosis of the IMFs $c_1$-$c_8$ obtained by EMD is shown in Table 3. $c_1$ and $c_2$, which have bigger kurtosis, are selected to be analyzed by the WVD. All the time-frequency distributions of different states are obtained, as shown in Fig. 8.

Fig. 7Gearbox vibration signal EMD results at t= 200 h

Fig. 8EMD-WVD results of channel 1 vibration signal at different condition time

a)$t=$ 1 h

b)$t=$ 100 h

c)$t=$ 200 h

d)$t=$ 300 h

e)$t=$ 400 h

f)$t=$ 442 h

It can be seen in Fig. 8 that the amplitude is increasing within 0.17 kHz to 0.37 kHz over time, so this frequency band is taken as a specific frequency band of gearbox, and the special frequency band energy collected from vibration signal is also taken as the health status value of gearbox.

In order to confirm the effect of feature extraction for vibration signal of gearbox, taking the vibration signal of the gearbox under the normal state without EMD as an example, the original signal is analyzed via WVD directly. For comparison, the result of WVD is given according to the contour, as shown in Fig. 9.

As compared with Fig. 9(b) and Fig. 9(a), it can be seen that the energy concentration of Fig. 9(b) whose distribution is got from the signal obtained by EMD is more obvious than Fig. 9(a) whose distribution is got from the original signal in the time-frequency domain. This is mainly because the original vibration signal is decomposed into IMFs with a single frequency by EMD, so the signal WVD obtained by EMD can restrain frequency interference and mode mixing effectively.

Fig. 9EMD-WVD comparison results of gearbox vibration signal at t= 200 h

a) WVD of the original signal

b) WVD of the signal obtained by EMD

Analyzing the result of Fig. 8, we can see that the amplitude of frequency band (0.71 kHz-0.37 kHz) is increasing with the gearbox working longer, the specific frequency band energy of the vibration signals is extracted based on energy analysis using wavelet packet as the special frequency band of gearbox, and the special frequency band energy of vibration signal $\left\{E_i,i=\text{1,}\text{2,...}\text{,}\text{442}\right\}$ is taking as the health status characteristic of gearbox, which is shown in Fig. 10.

Fig. 10Specific band energy of gearbox vibration signal

As we can see in Fig. 10, the green curve represents the real degradation process of gearbox, which is the wear of the gear tooth. The red curve and blue curve represent the fault feature extracted based on the EMD-WVD and fault feature extracted based on the WVD respectively. They all correspond with the degradation trend of gearbox, which demonstrates the accuracy of the methods. In addition, the fault feature based on the proposed method has a better stability compared with the fault feature based on the WVD, especially during the early failure and serious failure occurred. In the proposed method, some IMFs with larger kurtosis were selected, and the WVD was applied to these selected IMFs to attain the time-frequency distribution, so the noise, mode mixing and the cross-term interference were restrained. However, the WVD method can’t avoid the influence of these factors. Therefore, the proposed method has a better accuracy and effectiveness than the traditional method.

4.3. Gearbox remaining life prediction based on particle filter

The training samples for gearbox remaining useful life prediction based on particle filter are specific frequency band energy $E$ shown in Fig. 10 and gear teeth wear $wl$ shown in Table 2. To describe the degradation process of gearbox, the state space model based on Wiener process (Wiener-SSM) is proposed, and gearbox remaining useful life is obtained based on the degradation process model. The state space model for the gear teeth wear is as follows:

where $x_1$ is the gearbox worn state, $x_2$ is the unknown parameters related to the state which can’t be observable, $y\left(t\right)$is extracted specific frequency band energy, noise $w~N(0,Q)$ and $e~N(0,R)$is Gaussian distributed. Eq. (15) is used to describe the gearbox worn over time, parameters updating during the prediction process and the relationship between the gearbox worn and the fault feature, respectively. Take gearbox teeth wear at 442 hours as the failure threshold, namely $W{L}_f=\text{0.3}\text{mm}$. Assuming that the initial model parameters are as follows: $A=\text{1.00}$, $H=\text{1.00}$, $Q=\text{0}$, $R=\text{0}$, $x_0=\text{8000}$.

Taking the parameters estimation at 300 hours as an example, initial value $x_0=\text{8000}$ and historical data are both iterative based on PF to get the best estimation for the gearbox state $(x_{300}^s,p_{300}^s)$. Afterwards, all parameters of the state space model are calculated at 300 hours, as shown in Fig. 11.

Fig. 11The updated parameters at monitoring time 300 h

Model parameters are convergence to the optimum value quickly, which is the response of the gearbox SSM degradation trend, and will be used in the gearbox state feature prediction.

The parameters have randomness during the solving process, 1000 times repeated calculation is carried through to obtain the mean and standard deviation of parameters in this section, as shown in Table 4.

The probability density can be obtained based on the model parameters:

where $p\left(x_k\right|x_{k-1}^s)$ is the posterior probability under the conditions of $x_{k-1}^s$ at time $t_k$, and $w_{k-1}^s$ is the corresponding weight of particle $x_{k-1}^s$. According to the Eq. (13), the estimated parameters are shown in Table 5.

The probability density function of RUL and mean remaining life can be obtained based on the model parameters at different time, as shown in Fig. 12. The mean remaining life of gearbox is closer to the real remaining life with the increase of monitoring data. The variances of the probability density function are smaller over time, with more condition monitoring information is used, which means that the prediction results are more and more accurate.

Fig. 12Remaining useful life probability density function of gearbox

4.4. Performance evaluation

In order to evaluate the performance of remaining useful prediction based on PF, we take RMSE, MAE and VAE as indicators of absolute error, and take MARE and VRE as indicators of relative error. These indicators are shown as follows:

where $x_t$ is the real state, ${\widehat{x}}_t$ is the predicted state. The indicator is smaller means higher prediction precision.

The predicted results based on KF and PF are compared, as shown in Fig. 13, and the calculation results according to Eqs. (17-21) is shown in Table 6. The main conclusions are as follows. The RUL prediction method based on PF has a better precision and accuracy than based on KF. The prediction method based on PF considers degradation feature extracted from vibration signal, gearbox teeth wear and nonlinear and non-Gaussian system, which make it more suitable for gearbox remaining life prediction. However, the KP model is simple, and is easier to obtain the prediction results.

Fig. 13Comparison of two RUL prediction methods