Fault identification and severity assessment of rolling element bearings based on EMD and fast kurtogram

2.1. EMD method

The EMD method is based on such a simple assumption that any signal consists of different simple oscillatory modes. The definition and algorithm of EMD has been described in detail by Huang et al. [25]. According, to the EMD algorithm, any original signal $x\left(t\right)$ can be represented by a sum of IMFs:

where $c_i\left(t\right)$ denotes $i$th IMF and $r_n\left(t\right)$ is the final residue after all$\mathrm{}$IMFs are extracted. Each IMF should satisfy the following two conditions [27]:

a) In the whole data set, the number of extrema and the number of zero-crossings must either be equal or differ at most by one.

b) At any point, the mean value of the envelope defined by local maxima and the envelope defined by the local minima is zero.

The IMFs can be achieved by the following steps [25]:

(1) Find out all the local extreme values of the original signal and connect the local maxima and minima respectively to construct the upper and lower envelopes by a cubic spline interpolation.

(2) Specify the mean of the two envelopes as $m_1\left(t\right)$, and the difference between $x\left(t\right)$ and $m_1\left(t\right)$ is designated as $h_1\left(t\right)$, i.e.:

Ideally, $h_1\left(t\right)$ will be the first component of $x\left(t\right)$ if it is an IMF.

(3) If $h_1\left(t\right)$ is not an IMF, take it as the new original signal and repeat the steps (1)-(2); then a new mean $m_{11}\left(t\right)$ is calculated, so does the $h_{11}\left(t\right)$, i.e.:

Judge whether $h_{11}\left(t\right)$ is an IMF, if not, repeat sifting, i.e. up to $k$ times, until $h_{1k}\left(t\right)$ becomes an IMF, i.e.:

Set $h_{1k}\left(t\right)$ as $c_1\left(t\right)$ (the first IMF), i.e.:

(4) Subtract $c_1\left(t\right)$ from$\mathrm{}\mathrm{}x\left(t\right)$, the first residual $r_1\left(t\right)$ is got:

(5) Treat the residual $r_1\left(t\right)$ as the original signal, and repeat the steps (1)-(4) $n\mathrm{}$times and the other IMFs ($c_i$, $i=2,\dots ,n$.) will be obtained .The EMD decomposition will be stopped until the residual $r_n\left(t\right)$ becomes a monotonic function:

Finally, the original signal $x\left(t\right)$ is decomposed into $n$-emprical modes $c_1$,…, $c_n$ including different frequency bands ranging from high to low and a residual $r_n\left(t\right)$ which is the mean trend of $x\left(t\right)$.

2.2. Numerical simulation

To certify the ability of EMD to decompose the original signal into a sum of IMFs which can reveal the internal structure of the signal, a simulated signal containing multiple frequency components and some normal distributed random noise is constructed as Eq. (8):

where $f_1=\mathrm{}$100 Hz, $f_2=\mathrm{}$200 Hz, $f_1=\mathrm{}$300 Hz, and $\phi \left(t\right)$ is the Gaussian white noise.

The IMFs gained by EMD and the corresponding frequency spectrum are illustrated in Fig. 1. According to Fig. 1, it is obvious that the three frequency components belonging to the signal have been extracted successfully by EMD and are clearly shown in the first two IMFs. However, the last six IMFs are irrelevant with the original signal. These pseudo components are generated by over decomposition of EMD and may interfere with further signal processing such as feature extraction and fault diagnosis. Therefore, the IMF selection method is necessary to remove the pseudo components [28].

Fig. 1Demonstration of the IMFs

a) Waveforms in time domain

b) The corresponding frequency spectrums

2.3. IMF selection based on correlation analysis

EMD is based on the local characteristic timescales of a signal and can self-adaptively decompose the complicated signal into some IMFs reflecting the oscillatory mode embedded in the signal. However, in many cases, the high frequency components obtained from EMD are removed directly as noise, which may leave out some important information. In particular, for rolling element bearings, the impulsive signals related to the faults usually lie in the high frequency band. Therefore, selecting appropriate IMFs from EMD results is very necessary for faults diagnosis in different conditions.

Because of interpolation error, boundary effects and over decomposition of EMD, the pseudo components, which have nothing to do with the original signal, may appear in EMD. This may interfere the diagnosis results. So a method should be taken to find out these pseudo components and remove them. Wang et al. [28] introduced the correlation-based method which will be used in this paper.

Ideally, a signal $s\left(t\right)$ consists of a series of IMFs $c_i\left(t\right)$, which are all embedded in the original signal in theory, i.e.:

In practice, considering the fact that the pseudo components may be produced by EMD, the form of Eq. (9) will be changed as follows:

where $\overline{c_i}\left(t\right)$ and $x_j\left(t\right)$ denotes the true and pseudo IMF respectively.

The cross-correlation coefficient between $s\left(t\right)$ and $c_i\left(t\right)$ can be calculated by Eq. (11):

Because of the orthogonality of IMFs, Eq. (12) can be established:

Obviously, the pseudo components have little correlation with the original signal, that is, the correlation coefficient between $s\left(t\right)$ and $x_k\left(t\right)$ should be close to zero as Eq. (13) shows:

In conclusion, the correlation coefficient between the original signal and the true IMFs is approximate to the auto-correlation coefficient of the true IMFs themselves, while the correlation coefficient between the original signal and the pseudo IMFs is very small, close to zero. Therefore, it is very easy to distinguish the true IMFs and pseudo components.

2.4. IMF Energy percentage

The energy of the vibration signal will change with the conditions in which the bearings is operating. From Section 2.1, it can be seen that EMD can decompose the vibration signal into a series of IMFs which has the frequency band from high to low. So each ${IMF}_i$ ($i=$1,…, $n$) component corresponds to an energy $E_i$ ($i=$1,…, $n$) which forms an energy distribution in the frequency domain, and then the corresponding IMF energy percentage is designated as:

To demonstrate that energy distribution of IMF components will vary from the conditions in which the bearings is operating, four cases including normal, inner raceway fault (IRF), outer raceway fault (ORF), ball fault (BF) are considered. The vibration acceleration signals corresponding to each case is decomposed firstly by EMD. Next, the IMF energy percentages are calculated. Fig. 2 shows the original vibration signals and its first two IMFs of four cases respectively. The corresponding energy percentage is shown in Table 1.

From Table 1, it is seen that the IMFs energy distribution of four cases is different distinctly. For IMF1, the energy percentage increases from Normal to ORF, while, for IMF2, it is just the reverse. In particular, the energy percentage will change suddenly from the normal to the fault, such as the energy percentage of IMF1 being just 0.2304 in the normal, while it becomes 0.9574 when the bearings has ORF. So it can be concluded that energy percentages of IMFs based on EMD can basically reflect the work condition of the bearings.

Fig. 2The original vibration signal and its first two IMFs

a) Normal

b) IRF

c) ORF

d) BF

2.5. Fast Kurtogram

Kurtogram appeared firstly in [20], which derived from the so called SK. SK is very sensitive to the transient components in signal, and can also indicate accurately their types. According to Antoni [21], considering any stationary random process $Y\left(t\right)$, its Wold-Cramer representation can be defined as follows:

where $H\left(t,f\right)$ stands for the complex envelope of $Y\left(t\right)$ at time $t$ and frequency $f$, and $\mathrm{}dX\left(f\right)$ refers to the spectral increment associated with $X\left(t\right)$.($\mathrm{}X\left(t\right)$, it is a process which has a flat spectrum everywhere).

The fourth-order spectrum cumulant of $Y\left(t\right)$ is expressed as follows:

where $S\left(f\right)$ denotes the instantaneous spectrum moment, i.e.:

Thus, the SK can then be obtained by the normalized fourth-order cumulant:

In practical applications, the vibration signals$\mathrm{}x\left(t\right)$ measured from rotating machinery usually contain two parts: the fault-related signal $f\left(t\right)$ and the additive noise $n\left(t\right)$, i.e.:

Therefore, the SK of $x\left(t\right)$ can be given as follows:

where $\rho \left(f\right)=S_n\left(f\right)/S_f\left(f\right)$ is the noise-to-signal ratio, $S_n\left(f\right)$ and $S_f\left(f\right)$ is the power spectral density of $n\left(t\right)$ and $f\left(t\right)$ respectively.

According to Eq. (20), it can be concluded that when the noise is very weak, i.e. $\rho \left(f\right)$ is quite small, $K_x\left(f\right)\approx K_f\left(f\right)$; whilst the noise is very strong, i.e. $\rho \left(f\right)$ is very large, and $K_x\left(f\right)$ is close to zero. Therefore, the whole frequency domain can be detected when the SK is applied locally to different frequency bands, which indicates that the SK can not only detect effectively fault-related signals, but also locate their types in the frequency domain.

Fast Kurtogram, as the name suggests, represents a kurtosis vs frequency diagram produced by a fast algorithm. This fast algorithm [23] bases on the assertion that “each type of transient is associated with an optimal dyad (frequency/frequency resolution), which maximizes its kurtosis, and hence its detection”. To learn about dyadic wavelet decomposition algorithm, firstly, the original signal is decomposed into a series of frequency-band signals, the kurtosis value on each frequency band is calculated in the following; Secondly, a kurtosis vs frequency diagram is generated, from which where there is the biggest kurtosis value , and there is the optimal center frequency and bandwidth ;Thirdly, the signal is filtered based on the optimal center frequency and bandwidth ,and demodulated by square enveloping; Finally the characteristic frequency is gained by spectrum analysis and the faults can be diagnosed.

6.1. Bearings fault diagnosis

The experimental rigs and data records are described in detail in Section 4. In this paper, the inner raceway fault with three severities (0.1778 mm, 0.3556 mm, 0.5332 mm) at four running speeds (1797, 1772, 1750, 1730 rpm) are considered. The sampling frequency is 12 kHz and the data length is 65536.The vibration signals with four cases (case_0: Normal, case_1: 0.1778 mm, case_2: 0.3556 mm, case_3: 0.5332 mm) for IRF at 1730 rpm are decomposed into 16 IMFs by EMD firstly, and the cross-correlation coefficient (R) between ${IMF}_i$, $i=\mathrm{}$1,…, 16 and the original signal are calculated subsequently. Table 4 shows the result of the first 8 IMFs. It is known that, the correlation degree is low at $\left|R\right|<\mathrm{}$0.3 [30], which can be neglected. It can be seen from Table 4 that the first three IMFs for case_1 (0.2703 ≈ 0.3), the first two IMFs for case_2, and the first IMF for case_3, satisfy the condition of $\left|R\right|\ge \mathrm{}$0.3. So the representative IMFs for the three cases are the first three, the first two, and the first component(s) respectively. The reconstruction signals (RS) are generated by the addition of the representative IMFs for each case. The results applying Fast Kurtogram method to RS for each case are shown as Fig. 6, and the corresponding envelope analysis are shown in Fig. 7, Fig. 8, and Fig. 9 respectively. The maximum kurtosis value $K_{max}$, the optimum bandwidth $B_w$ and the optimum center frequency $f_c$ of the four cases are shown in Table 5. In Fig. 7(c), Fig. 8(c), Fig. 9(c), the characteristic frequencies (154.9 Hz, 155.9 Hz, 155.6 Hz) and their multiplications are marked, and this indicates clearly the inner raceway fault (BPFI=156.1 Hz at 1730 rpm in Table 3).

It is can be seen from Table 4 that the cross-correlation coefficient between the first IMF (IMF1), and the original signal, for all the three cases, is the highest in the representative IMFs selected above. This indicates that the IMF1 is the most fault-related component in all IMFs decomposed by EMD. So the energy percentage of IMF1 is selected as the indicator of the fault severity assessment. To make the results convincing, in this experiment, three cases (case_1: 0.1778 mm, case_2: 0.3556 mm, case_3: 0.5332 mm) for IRF, and case_0 for the normal at four running speeds (1797, 1772, 1750, 1730 rpm), 16 cases in total (as Table 6), are considered. The length of data for each case is 10240 which is divided into 10 groups with the length of 1024. For each group of data, it is decomposed into 10 IMFs by EMD, and the energy percentage of each IMF is also calculated. Similarly, the energy percentages of each IMF for the other 9 groups can be gained as well. The final result is the average of that of 10 groups of data for each case. Table 6 shows the energy percentage of IMF1 for all the cases. The corresponding bar chart and the trend curve are shown in Fig. 10. As is shown in Fig. 10, it is obvious that the energy percentage of IMF1 increases with the fault severity, meaning the fault severity becoming more serious when the energy percentage of IMF1 increases. The change is relatively smooth from a light fault to a heavy fault, but it changes sharply from normal to the fault, such as, at 1750 rpm, the percentage is just 0.1585 in case_0, while it becomes 0.7949 which is more than five times of the former in case_1. Furthermore, the results are similar for all four running speeds, which may indicate to a certain extent that the result is not affected by the running speed of the bearings.

Fig. 6The Kutogram of RS at 1730 rpm

a) Case_0

b) Case_1

c) Case_2

d) Case_3

According to the application to the inner raceway fault, it demonstrates the effectiveness of the proposed method in the rolling element bearings inner raceway fault type identification and severity assessment. This method shows the great potential in the application of condition monitoring of bearings. It provides us a possibility to judge whether there is a defect on bearing or not and predict the trend of the faults by monitoring the energy percentage of IMF1. For instance, we can try to calculate the threshold of the energy percentage of IMF1 by a certain amount of experimental data of the normal bearings at a certain condition. If the monitoring result exceeds the threshold, the monitored bearing may have a fault. Further, if the monitoring result becomes more and more big, it indicates that the bearing gets worse and we need to replace it.

Fig. 7IRF in case_1 at 1730 rpm

a) RS

b) The envelope

c) The envelope spectrum

6.2. Comparison with Hilbert envelope demodulation

In order to verify the superiority of the proposed method in fault identification, when analyzing the vibration signal shown above, another signal processing technique, envelope demodulation based on Hilbert transform has been applied too, but it did not extract the fault characteristics effectively. For instance, after applying Hilbert envelope demodulation method to the vibration signal at case_1 (0.3556 mm) for IRF at 1730 rpm, the result is shown in Fig. 11. As shown in Fig. 11, although the IRF characteristic frequency (155.9 Hz) and its multiplication (311.8 Hz) can be extracted, there exists some conspicuous and unfathomed frequency components, such as 69.5 Hz and 242.2 Hz marked in Fig. 11, which may easily interfere with diagnostic results and even possibly lead to the wrong ones. By contrast, the characteristic frequencies shown in Fig. 7(c), Fig. 8(c), Fig. 9(c), are very prominent, and the other frequency components are relatively weak which will not affect the diagnostic results.

Fig. 8IRF in case_2 at 1730 rpm

a) RS

b) The envelope

c) The envelope spectrum

Fig. 9IRF in case_3 at 1730 rpm

a) RS

b) The envelope

c) The envelope spectrum

Fig. 10The percentage of IMF1 in different conditions

a) The bar char

b) The trend curve

Besides, another superiority of the proposed method is that it can be not only used to identify the fault patterns but also to assess the fault severities which is also essential for condition monitoring of bearings. To indicate the different severities of fault, a representative indicator is needed to show the variation tendency. However, the envelope demodulation based on Hilbert transform just gives the characteristic frequency which is constant for the same fault with different severities. So Hilbert envelope demodulation fails to assess the fault severities. In contrast, the energy percentage of IMF1 is a good indicator for fault severity assessment.

Fig. 11The Hilbert envelope demodulation spectrum