Novel complete ensemble EMD with adaptive noise-based hybrid filtering for rolling bearing fault diagnosis

1. Introduction

Bearing is a critical component in rotary machine. The reliability directly affects machine health status. Thus, an effective condition monitoring (CM) and fault diagnosis (FD) method is required to detect and recognize potential abnormalities which reduce risk of unexpected severity and other disastrous events [1-4].

Vibration analysis technology is a key topic in CM and FD research domain. This technology may extract fault feature of rolling bearing effectively. Conventional methods, such as fast Fourier transform (FFT), are used to address stationary and linear signals [5] widely. However, FFT is unsuitable to extract fault feature on account of characteristic of non-linear and non-stationary for vibration signal. Wavelet transform is a typical application in bearing fault diagnosis [6-8]. This technique has attracted considerable attention of researchers in the past two decades. However, wavelet transform is difficult to select wavelet base function. Recently, empirical mode decomposition (EMD) has been proposed to process complex non-linear and non-stationary signal [9]. EMD is an adaptive decomposition method that decomposes complicated signals into a set of multi-component with intrinsic mode function (IMF). Every IMF indicates a mono-component amplitude- and frequency-modulated signal. EMD has been applied to bearing [10], gearbox [11], and rubbing [12] fault diagnoses. However, EMD have several limitations, such as mode mixing and end effect. To solve the problem, Huang [13] proposed an ensemble EMD (EEMD) to alleviate the limitations of EMD. This method reflected the frequency of every IMF accurately and effectively. Moreover, EEMD has attracted significant attention in fault diagnoses domain. For example, Wang [14] proposed a method based on multi-scale principal component analysis to extract fault frequency of slewing bearing. Chen [15] improved EEMD by adding white noises and adaptively determining an ensemble number to extract fault features. However, EEMD has caused spurious mode and residual noise because of interaction of noise and signal. Subsequently, Torres [16] introduced a complete EEMD with adaptive noise (CEEMDAN) to overcome this drawback. CEEMDAN has been applied to fault diagnoses [17, 18] successfully.

For vibration signal being immersed by heavy noise, it is important to develop a de-noised method to extract the early detection of weak fault. However, filtering method based on EMD or its improved versions (EEMD and CEEMDAN) do not filter out noise in vibration signal fully. The reason is that IMF do not effectively reflect frequency component of fault because of noise interference. Recently, the hybrid filtering method based on singular-value decomposition (SVD) and EMD or its improved versions have been proposed to extract weak fault feature. Cheng [19] used EMD to decompose impulse signal, and used SVD and support vector machine to extract fault feature and classify fault pattern of gearbox and bearing. Yang [20] adopted EEMD with sample entropy and SVD to diagnose fault by acoustic emission signal. It is found that SVD is an effective method to filter out noise and extract fault feature. However, important fault information losses because of spurious mode and heavy noise interference in IMF.

To solve this problem, this study introduces a multi-layer filtering (hybrid filtering) method to extract weak fault feature of bearing. In the first layer, vibration signal is decomposed to obtain IMFs by CEEMDAN. Then, a de-trended fluctuation analysis (DFA) is utilized to identify effective IMF and spurious mode adaptively. In second layer, a non-local mean (NLM) de-noised method is used to filter out noise in each IMF. In last layer, SVD is used to remove residual noise in vibration signal. The hybrid filtering method is applied to simulated and measured vibration signal from rolling bearing. The result indicates that this method extracts weak fault effectively and operate well for bearing fault diagnose.

The remainder of this paper is organized as follows: Section 2 introduces the CEEMDAN algorithm. Section 3 describes the hybrid filtering algorithm. Section 4 verifies the effectiveness of the proposed algorithm for fault extraction. Section 5 presents the conclusion.

2. CEEMDAN algorithm

Traditional EEMD have high computational cost and residue noise because of interaction of signal and added noise. To solve this problem, CEEMDAN is introduced in reference [16]. Principle of CEEMDAN is as follows:

(1) Generate the noisy signal $x^i\left(t\right)$, $x^i\left(t\right)=x\left(t\right)+{\omega }_0{\xi }^i\left(t\right)$ ($i=$1,…, $N$), where ${\xi }^i\left(t\right)$ is white noise with unit variance, and ${\omega }_0$ is coefficient of added white noise.

(2) Obtain the first IMF ($c_1^i\left(t\right)$) of each noisy signal $x^i\left(t\right)$ by EMD;

(3) Obtain the first decomposed IMF ($c_1\left(t\right)$) by taking average of each $c_1^i\left(t\right)$:

(4) Compute the first residue by Eq. (2);

(5) Apply EMD to decompose $r_1\left(t\right)+w_1{E}_1\left({\epsilon }^i\left(t\right)\right)$, and extract the first IMF to obtain the decomposed second IMF ($c_2\left(t\right)$), where $E_m\left(\cdot \right)$ signifies that EMD is used to decompose $x^i\left(t\right)$ to get the $m$th IMF:

(6) Compute the $m$th residue ($m=$ 2,…, $K$), and extract the first decomposed IMF to obtain $\left(m+1\right)$th decomposed IMF by Eq. (4):

(7) Repeat above steps until the residue $R\left(t\right)$ contain less than two extrema:

The signal $x\left(t\right)$ can be written as:

3. Hybrid filtering method

3.1. Selection of spurious IMFs

When CEEMDAN is used to decompose vibration signal, the obtained IMF contains periodic impulsive component (IPIC) and non-IPIC (NIPIC). The IMF with impulsive component can reflect physical meaning of vibration signal, and the residue is a non-effective IMF, which is called “spurious mode” or “spurious IMF”. It is necessary to develop a method to distinguish IPIC and NIPIC effectively. It is found that there is obvious fluctuation difference between IPIC and NIPIC. Recently, DFA [21-23] has been regarded as analysis tool of scale and is suitable to analyze fluctuation of non-linear and non-stationary time series. Thus, DFA is proposed to evaluate the fluctuation of IPIC and NIPIC. The following is the principle of DFA:

(1) Generate the integrated time series, where $\stackrel{-}{x}$ indicates the mean of time series $x\left(i\right)$:

7

$y\left(k\right)=\sum _{i=1}^k\left[x\left(i\right)-\stackrel{-}{x}\right],k=\mathrm{1,2},\dots ,N.$

(2) Divide $y\left(k\right)$ into $n$ length sections;

(3) Use least square method to determine the local trend $y_n\left(k\right)$;

(4) Get fluctuation function $F\left(n\right)$ by $y_n\left(k\right)$:

8

$F\left(n\right)={\left(\frac{1}{N}\sum _{k=1}^N{\left[y\left(k\right)-y_n\left(k\right)\right]}^2\right)}^{1/2}.$

(5) Obtain different $F\left(n\right)$ using different length segments;

(6) Calculate slope between $\mathrm{l}\mathrm{o}\mathrm{g}\left(F\left(n\right)\right)$ and $\mathrm{l}\mathrm{o}\mathrm{g}k$, where the slope is called Hurst exponent $\alpha$:

9

$F\left(n\right)\propto n^{\alpha }.$

Reference [24] indicated that the larger value of Hurst exponent ($\alpha$) is, the smoother a time series is. Numerous experiments found that NIPIC is smoother than IPIC. It means that Hurst exponent ($\alpha$) decreases with the increase of IMF index. When the Hurst exponent ($\alpha$) of the $\left(k+\text{1}\right)$th IMF become large suddenly, and the IMF that corresponds to the index is a spurious mode. Then, effective vibration signals with impulsive components are adaptively constructed by subtracting first $k$th IMFs. The identified steps are summarized as follows:

(1) Decompose signal $x\left(t\right)$ into a set of IMFs by CEEMDAN:

10

$x\left(t\right)=\sum _{i=1}^{k}IM{F}_i+\sum _{i=k+1}^{N}IM{F}_i+R_N\left(t\right).$

(2) Calculate Hurst exponent (${\alpha }_i$) of every IMF by DFA ($i=$ 1,…, $N$);

(3) Identify index of the relevant $k$th from ${\alpha }_i$:

11

$k=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}\left({\alpha }_i\right).$

(4) Extract all IMFs up to index $k$th.

To verify the effectiveness of proposed method, the following simulated vibration signal [25] is used to identify IPIC and NIPIC:

12

$y\left(t\right)=x\left(t\right)+n\left(t\right),$

where, $y\left(t\right)$ is noisy signal, $x\left(t\right)$ indicates the real signal and $n\left(t\right)$ is the noise. Eq. (13) and (14) is expression of $x\left(t\right)$:

13

$x\left(t\right)=A{e}^{-at}\mathrm{s}\mathrm{i}\mathrm{n}\left(2\pi f_{c1}kT\right)+B\mathrm{s}\mathrm{i}\mathrm{n}\left(2\pi f_{r}kT\right)\mathrm{s}\mathrm{i}\mathrm{n}\left(2\pi f_{c2}kT\right),$

14

$t=\mathrm{m}\mathrm{o}\mathrm{d}\left(kT,1/f_m\right),$

here, $a$ is the modulated factor, and $T$ is the sample period. $k$ is the sample points ($k=$ 1000), $f_m$ is fault frequency of bearing, $f_r$ is the rotation frequency, $f_{c1}$ and $f_{c2}$ are the carrier frequencies from system resonance. Table 1 are the parameters of simulated vibration signal. Where the sample rate is 1.6 kHz, and the input signal-to-noise ratio $(SN{R}_{in}$) is −10 dB.

Table 1Parameters of the simulated vibration signal

$a$	$T$ (s)	$f_m$ (Hz)	$f_r$ (Hz)	$f_{c1}$ (Hz)	$f_{c2}$ (Hz)
800	1/16000	160	30	2000	7000

Fig. 1(a) and (b) are the simulated noisy signal and the corresponding fast Fourier transform (FFT), respectively. It is found that FFT cannot identify characteristic frequency of simulated fault signal, and extract only the carrier frequencies ($f_{c1}$ and $f_{c2}$).

To extract the fault frequency, CEEMDAN is first used to decompose noisy vibration signals. Fig. 2 are the decomposed IMFs. It is found that there are nine IMFs and residue. The first five IMFs contain PICs obviously, and periodic information of fifth IMF is shown red box in the fifth IMF. However, it is obvious that the sixth IMF mainly contains NPIC. To identify the PIC, DFA is used to calculate Hurst exponent ($\alpha$) of each IMF to distinguish PIC and NPIC. Fig. 3 demonstrates the DFA of IMF1, and Fig. 4 exhibits the Hurst exponent ($\alpha$) for all IMFs. It is found that Hurst exponent ($\alpha$) decreases with increase of IMF index. The Hurst exponent ($\alpha$) of IMF6 rapidly become large when Hurst exponent ($\alpha$) of IMF5 is smallest. It indicates that the spurious modes are from IMF6 to IMF10 (expected residue). The identified result is in accordance with the analysis result in Fig. 2. It means that the proposed DFA method is suitable to identify IPIC and NIPIC. Fig. 5 presents the combination of signals with first five IMFs. It is found that the signal still contains heavy noise.

Fig. 1Noisy signal and corresponding FFT

a) Noisy signal

b) Noisy signal FFT

Fig. 2CEEMDAN decomposed IMFs for simulation vibration signal

a)

b)

Fig. 3DFA of IMF1

Fig. 4Hurst exponent of each IMF

Fig. 5Summed signal of the first five IMFs

3.3. Singular value decomposition de-noising

3.3.1. Singular value decomposition principle

The final layer filtering utilizes SVD to remove extra noise. SVD demonstrates an excellent performance in noise removal and is extensively used in fault extraction of rolling bearing [30-32]. It is assumed that noisy series $x\left(k\right)$ is expressed as:

19

$x\left(k\right)=\left[x\left(1\right),x\left(2\right),\dots x\left(N-1\right),x\left(N\right)\right],\mathrm{}\mathrm{}\mathrm{}k=1,\dots ,N.$

Fig. 6Second filtering layer

Signal $x\left(k\right)$ is consisted of real signal and noise. It is expressed as follows:

20

$x\left(k\right)=s\left(k\right)+n\left(k\right),$

where, $s\left(k\right)$ is real signal, and $n\left(k\right)$ is noise. The following is Hankel matrix $H\left(H\in R^{p\times q}\right)$ of signal $x\left(k\right)$:

21

$H={\left[\begin{array}{cccc}x\left(1\right)& x\left(2\right)& \cdots & x\left(q\right)\\ x\left(2\right)& x\left(3\right)& \cdots & x\left(q+1\right)\\ ⋮& ⋮& \cdots & ⋮\\ x\left(p\right)& x\left(p+1\right)& \cdots & x\left(N\right)\end{array}\right]}_{p\times q},$

where, $p+q-1=N$, $p\ge q$, and $H\left(i,j\right)=x\left(i+j-1\right)$. Eq. (21) is furtherly expressed as:

22

$H={\left[\begin{array}{cccc}s\left(1\right)& s\left(2\right)& \cdots & s\left(q\right)\\ s\left(2\right)& s\left(3\right)& \cdots & s\left(q+1\right)\\ ⋮& ⋮& \cdots & ⋮\\ s\left(p\right)& s\left(p+1\right)& \cdots & s\left(N\right)\end{array}\right]}_{p\times q}+{\left[\begin{array}{cccc}n\left(1\right)& n\left(2\right)& \cdots & n\left(q\right)\\ n\left(2\right)& n\left(3\right)& \cdots & n\left(q+1\right)\\ ⋮& ⋮& \cdots & ⋮\\ n\left(p\right)& n\left(p+1\right)& \cdots & n\left(N\right)\end{array}\right]}_{p\times q}.$

The Eq. (22) is simplified as:

23

$H=H_s+H_n,$

here, $H_s\left(i,j\right)=s\left(i+j-1\right)$ and $H_n\left(i,j\right)=n\left(i+j-1\right)$ are Hankel matrices of the real signal and noisy signal, respectively. Thus, the main task is to found $H_s$ to obtain de-noised signal $\widehat{s}\left(k\right)$. Here, $H$ is first decomposed into the following expression:

24

$H=UD{V}^T=\sum _{i=1}^r{\sigma }_i{u}_i{v}_i^T=\sum _{i=1}^r{\sigma }_i{A}_i,$

here, $U\left(U\in R^{p\times q}\right)$ and $V\left(V\in R^{p\times q}\right)$ are two orthogonal matrices, $D$ is a diagonal matrix and defined as $D=\left[\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\sigma }_1,{\sigma }_2,\cdots ,{\sigma }_q\right),O\right]$, $O$ is a zero matrix, and ${\sigma }_1>{\sigma }_2,\cdots ,{\sigma }_i,\cdots ,>{\sigma }_q>0$. Here, ${\sigma }_i$ is assumed to be last singular value of Hankel matrix $H_s$. Other values are set to zero if the previous $i$ singular values are preserved. Then, reverse re-construction can obtain the de-noised signal. Thus, parameter $i$ must be determined. The next part uses a relative energy difference spectrum of a singular value to identify important parameter $i$.

3.3.2. Relative energy difference spectrum

In principle, real signal $s\left(k\right)$ and noise $n\left(k\right)$ are non-correlated. The energy of the real signal is concentrated. However, the noise is relatively decentralized. It indicates that the singular values of real signal components are rounded to certain numerical value, and singular value of noise components are around other numerical value. It means that the energy of first $i$ singular value (singular value of real signal component) is different from energy of the others ($r-i$) singular value (noise component). If the energy difference of adjacent singular value is done, it is found that the energy difference of singular value are almost same as energy for real signal component or noise component, respectively. However, the energy difference of singular value is biggest in demarcation point of real signal component and noise component. It proves that the relative energy difference spectrum of SVD changes at the demarcation point of singular value of the real and noise signal components abruptly. This paper introduces the following relative energy difference spectrum to identify singular value ${\sigma }_i$:

25

$\rho \left(i\right)=\frac{{\sigma }_i^2-{\sigma }_{i+1}^2}{\mathrm{m}\mathrm{a}\mathrm{x}\left({\sigma }_i^2\right)-\mathrm{m}\mathrm{i}\mathrm{n}\left({\sigma }_i^2\right)}\times 100\mathrm{\%},\mathrm{}\mathrm{}\mathrm{}\mathrm{}i=1,\dots ,r,$

when the singular value of real signal components are extracted by identifying demarcation point by proposed method, the de-noised signal can been obtained by Eq. (24). It means that the signal can been re-constructed by finding the maximum ${\sigma }_i$. Here, the size of Hankel matrix is determined on the basis of Reference [33].

Now, the proposed de-noising method is used to filter out residual noise of vibration signal in Fig. 6. Fig. 7 displays singular value of Hankel matrix. Figs. 8(a) and (b) depict the local large of singular value (first 50 singular values are displayed) and the corresponding relative energy difference spectrum, respectively. It is found that the 10th singular value is the maximum in all different spectra. Thus, the first 10 singular values are selected to re-construct the de-noised signal. Fig. 9 are comparison result of second-layer and third-layer filtered signal. It is found that periodical information of the third layer filtered signal is more obvious than that of the second layer.

Fig. 7Singular value

Fig. 8Singular value and corresponding difference spectrum

a) Singular value

b) Difference spectrum of singular value

Fig. 9Comparing of filtered signal with second the third layer filtering method

Analysis technique of envelope spectrum is used to extract fault feature. Fig. 10 presents the envelope spectrum of the filtered signal. It is found that fault characteristics $\text{1}X$ and $\text{2}X$ (159.7 and 319.4 Hz, respectively) and rotational frequency $\text{2}X$ (63.87 Hz) are extracted in the envelope spectrum. It indicates that the proposed hybrid filtering method is effective in extracting fault frequency of bearing.

Fig. 10Envelope spectrum of hybrid filtered signal

Finally, the three-layer de-noising method is presented as follows:

(1) Decompose vibrated signal $x\left(t\right)$ into a set of IMFs by CEEMDAN;

(2) Select the real and spurious IMF component by DFA adaptively;

(3) Filter out noise in each real IMF by NLM;

(4) Combine the filtered IMFs to obtain effective signal $x^{\mathrm{\text{'}}}\left(t\right)$;

(5) Use SVD to filter out noise to obtain signal $x^{\mathrm{\text{'}}\mathrm{\text{'}}}\left(t\right)$;

(6) Calculate the envelope spectrum of signal $x^{\mathrm{\text{'}}\mathrm{\text{'}}}\left(t\right)$;

(7) Extract the weak feature of vibration signal $x^{\mathrm{\text{'}}\mathrm{\text{'}}}\left(t\right)$;

Fig. 11 is flow chart of the hybrid filtering method.

Fig. 11Flow chart of hybrid filtering method

4. Results and discussions

In this section, the fault data from the website of Case Western Reserve University [34] is used to evaluate the performance of proposed method. Fig. 12 exhibits the experiment rig which consists of a 2 HP motor, a torque transducer, a dynamometer, and control electronics. The test bearing is 6205-2RS SKF deep groove ball bearing. Table 2 is the bearing parameters.

Fig. 12Bearing experiment rig: (1) – motor, (2) – torque transducer, (3) – dynamometer

The theoretical characteristic frequency of inner race $f_{in}$ can be calculated using the following expression:

where, $f_r$ is the rotation frequency.

In this experiment, single point fault is introduced to the test bearing using electro-discharge machining with fault diameters of 0.018, 0.036, 0.054, and 0.072 mm. Fault data from fault diameter of 0.018 mm is used to identify fault recognition as early as possible. The sample rate is 12 kHz. The selected data length is 1000. Figs. 13(a) and (b) display the time domain signal and corresponding FFT, respectively. It is found that FFT of vibration signal do not extract frequency feature of vibration signal, and the vibration signal contain heavy noise. Here, the proposed method is first used to remove noise by flow chart in Fig. 11, then analysis technology of envelope spectrum of vibration signal is employed to extract fault feature. Moreover, typical filtering methods, namely, EMD-SVD [25], SVD, and NLM, are also used to extract fault feature for the same vibration signal. The filtered result and the corresponding envelope spectrum are demonstrated in Figs. 14(b), (c), and (d). It is found that more noise is observed using EMD-SVD, SVD, and NLM methods than proposed method (CEEMAN-DFA-NLM-SVD) for a time domain waveform. Moreover, for envelope spectrum, it also has similar result and contains a noisy frequency component (high-frequency component or interference frequency as shown elliptical annotation in red box in Fig. 14) among EMD-SVD, SVD, and NLM methods (envelope spectrum is only displayed in the first 2 kHz). However, the envelope spectrum of the proposed method does not reflect a high frequency and only reflect fault frequency of bearing (low-frequency information in Fig. 14). It indicates that the proposed method more effectively remove noise than other de-noising methods. The envelope spectrum of the CEEMAN-DFA-NLM-SVD presents the maximum frequency (155.7 Hz) and multiplication frequency (311.4 Hz). The identified frequency is nearly consistent with the theoretical frequency (151.44 Hz) by Eq. (26). For the EMD-SVD, SVD method, although the fault frequency and corresponding multiplication frequency can be reflected in Fig. 14, the multiplication frequency is difficult to identify because of interference frequency. For NLM method, it is not found fault frequency and multiplication frequency in the envelope spectrum, and it also produce the interference frequency in red box in Fig. 14. In addition, the amplitude of envelope spectrum of proposed method is bigger than the other methods. It makes the fault frequency be more identified. Based on above analysis, the proposed method demonstrates an outstanding performance than the other filtering methods.

Fig. 13Noisy signal and corresponding FFT

a) Noisy signal

b) Noisy signal FFT

Fig. 14Comparing of filtered results and feature extraction

5. Conclusions

This study proposes a novel hybrid filtering method based on the CEEMDAN-DFA-NLM-SVD to extract weak fault features under heavy background noise. On the basis of the abovementioned research results, the conclusions are as follows:

1) The first layer filtering introduces criterion of the smallest Hurst exponent to identify the periodic impulse among the IMFs by CEEMDAN adaptively.

2) The second filtering layer uses NLM to filter out noise in each IMF that contains a PIC.

3) In the final layer filtering, the relative energy difference spectrum, which identifies the singular value of Hankel matrix of real signal, is applied to adaptively remove noise for the second-layer filtered signal.

Finally, the proposed method is used to effectively filter out noise and extract weak fault features of the simulation and measured signal. Furthermore, the method is also compared with typical filtering methods, respectively. The result indicates that the proposed method demonstrates the optimum performance among the filtering methods for extracting fault feature of bearing.