J. Cent. South Univ. (2016) 23: 1682-1689
DOI: 10.1007/s11771-016-3222-x
Rolling element bearing instantaneous rotational frequency estimation based on EMD soft-thresholding denoising and instantaneous fault characteristic frequency
ZHAO De-zun(赵德尊), LI Jian-yong(李建勇), CHENG Wei-dong(程卫东),
WANG Tian-yang(王天杨), WEN Wei-gang(温伟刚)
School of Mechanical Electronic and Control Engineering, Beijing Jiaotong University, Beijing 100044, China
Central South University Press and Springer-Verlag Berlin Heidelberg 2016
Abstract: The accurate estimation of the rolling element bearing instantaneous rotational frequency (IRF) is the key capability of the order tracking method based on time-frequency analysis. The rolling element bearing IRF can be accurately estimated according to the instantaneous fault characteristic frequency(IFCF). However, in an environment with a low signal-to-noise ratio (SNR), e.g., an incipient fault or function at a low speed, the signal contains strong background noise that seriously affects the effectiveness of the aforementioned method. An algorithm of signal preprocessing based on empirical mode decomposition (EMD) and wavelet shrinkage was proposed in this work. Compared with EMD denoising by the cross-correlation coefficient and kurtosis(CCK) criterion, the method of EMD soft-thresholding (ST) denoising can ensure the integrity of the signal, improve the SNR, and highlight fault features. The effectiveness of the algorithm for rolling element bearing IRF estimation by EMD ST denoising and the IFCF was validated by both simulated and experimental bearing vibration signals at a low SNR.
Key words: rolling element bearing; low signal-to-noise ratio; empirical mode decomposition soft-thresholding denoising; instantaneous fault characteristic frequency; instantaneous rotational frequency
1 Introduction
There is prevalent equipment that functions under non-stationary operating conditions(varying speed, varying load, etc.) in petrochemical engineering, mining industries, and other fields [1]. To improve the reliability and the safety of the system and reduce the downtime of the equipment, there has been extensive interest in research on bearing diagnosis under non-stationary operating conditions [2]. Order tracking is one of the most effective algorithms that can be used to remove the effects of speed fluctuation. The essence of order tracking is resampling the original vibration signal at a constant angle increment which converts the non-stationary signal (in time domain) into the stationary one (in angular domain) [3]. Traditional order analysis techniques can be divided into two types: hardware- based order tracking and computed order tracking. The hardware-based order tracking has greater hardware requirements. Compared with hardware-based order tracking, computed order tracking allows a simpler rotational speed measuring device, however it cannot work without devices such as tachometers or encoders. In order to overcome the limits of traditional order analysis techniques with respect to the equipment installation and costs, estimating the instantaneous rotational frequency(IRF) directly from vibration signals without tachometers or encoders under the framework of order analysis techniques has become a new focus for many researchers [4-5].
IRF estimation algorithms can be classified into two types: the first is based on phase demodulation [6-7], and the second is based on time-frequency representation. The second algorithm is more suitable for operating conditions with large speed variation [8]. The principle of this algorithm is to extract the time-frequency component, which is related to the IRF, and estimate the IRF. BOASHASH [9-10] discussed the concept of the instantaneous frequency(IF) in detail and estimated it from a time-frequency representation. GUO et al [11] estimated a motor’s IRF using the IF, which is extracted from the time-frequency representation of the motor’s vibration signal using the peak searching algorithm. ZHAO et al [12] combined the short-time Fourier transform (STFT) and the Viterbi fitting algorithm to extract the IF that is related to the horizontal spiral centrifuge’s IRF. It can be concluded that the accuracy of the extracted IF that is related to the IRF, which will directly decides the applicability of the algorithm for IRF estimation based on a time-frequency representation. However, the algorithm for IRF estimation based on the time-frequency spectrum has no effect on rolling element bearing vibration signals under variable rotational speeds. Therefore, the concept of the instantaneous fault characteristic frequency (IFCF) was proposed by WANG et al [13] and applied to estimate the rolling element bearing IRF. The algorithm for rolling element bearing estimation based on the IFCF firstly extracts the IFCFs from a time-frequency representation of the envelope signal by the peak searching, and then calculates the IRF by dividing the IFCF trend by the fault characteristic coefficient (FCC, Cfc).
However, in many cases, the rolling element bearing has an incipient fault or function at a low rotational speed. Under these conditions, the fault impulses are covered by the background noise, and the SNR of the signal is very low. The accuracy of the IRF, which is extracted using the algorithm for the bearing IRF estimation based on the IFCF, decreases because of the noise. Hence, the aforementioned algorithm has limited applications.
Therefore, this work aims to correct the algorithm for the bearing IRF estimation based on the IFCF using a denoising method, which combines empirical mode decomposition (EMD) and wavelet soft-thresholding (ST) denoising to preprocess raw signals. In the method of EMD ST denoising, the raw signal is firstly decomposed by EMD, which generates a collection of intrinsic mode functions(IMFs). Then the illusive IMFs are removed. The IMFs are selected which are contaminated by noise, and the noise is removed via the soft-thresholding denoising method. Finally the denoising IMFs are added to the rest of the IMFs to reconstruct the signal. The improved algorithm, namely the algorithm for rolling element bearing IRF estimation based on EMD ST denoising and the IFCF, can reduce the noise interference effectively, compensate for the inadequacy of the original algorithm, and expand its applications.
2 Algorithm for IFCF based rolling element bearing IRF estimation and existing problems
When a rolling element bearing has a local fault, it generates an impulse that recurs at a constant interval and repeats at the fault characteristic frequency (FCF) under a constant rotational speed. Relatively obvious amplitudes also recurs in the envelope spectrum of the vibration signal. The FCCs of outer race, inter race and ball are expressed by Eqs. (1)-(3), respectively:
(1)
(2)
(3)
where Cfco, Cfci and Cfcb are the FCCs of the outer race, inner race and ball, respectively; fo, fi and fb are FCFs of the outer race, inner race and ball, respectively; fr is the rotational speed; and n, d, D and α represent the number of rolling elements, diameter of the rolling element, pitch diameter and contact angle, respectively.
The above formulas indicate that when a rolling element bearing has a local fault, its FCF has a fixed relationship (i.e., the FCC) with the rotational speed. When the bearing type is determined, the FCCs of the outer race, inner race and ball are also determined. If the rotational speed is changing, the FCF will change with the rotational speed by fixed values, i.e., the FCC. Each moment of the FCF is called an IFCF. Because the IFCF has relatively obvious amplitudes in the envelope spectrum of the vibration signal, the IFCFs are easily extracted by the peak searching algorithm, and then the IRF of the rolling element bearing is estimated.
The algorithm for rolling element bearing IRF estimation based on the IFCF can be mainly divided into three steps. The first is raw vibration signal filtering by a fast kurtogram. The second involves obtaining the time- frequency representation of the envelope signal by a Hilbert transform and STFT, extracting the IFCFs by an amplitude-sum based peak searching algorithm, and then determining the IFCF trend. The third step is to calculate the IRF by dividing the IFCF trend by the FCC.
When the fault impulse is very mild or the rotational speed is low, the impulse is easily contaminated by the background noise. In low SNR conditions, the algorithm for rolling element bearing IRF estimation based on the IFCF may not work well. After white Gaussian noise is added to the simulated signal, the SNRs are -6 and -9.5 dB, respectively. The IRFs estimated using the algorithm for rolling element bearing IRF estimation based on the IFCF are shown in Fig. 1, where the squares represent the estimated IRF with SNR of -6 dB, the dotted line indicates the estimated IRF with SNR of -9.5 dB, and the solid line is the theoretical IRF. The IRFs for both SNRs deviate from the theoretical value. The IRF accuracy with SNR of -9.5 dB is worse than that with SNR of -6 dB. As previously discussed, the algorithm for rolling element bearing IRF estimation based on the IFCF generates errors at a low SNR, and as the SNR decreases, the errors in the IRF increase. Hence, to improve the accuracy of the IRF estimation, and expand its applications, it is necessary to remove the noise of the raw vibration signal.
Fig. 1 Calculated IRF with SNR of -6 and -9.5 dB, along with theoretical IRF
3 EMD ST denoising method
A variety of denoising methods such as weighted denoising, the median denoising, the Butterworth filter and the Chebyshev filter have been proposed. However, these methods are mainly used for stationary signals, their effectiveness for the non-stationary signals is not obvious. With the emergence of EMD and the wavelet transform, the EMD denoising method and the wavelet shrinkage method have yielded significant effects.
EMD is a self-adaptive analysis method proposed by HUANG et al [14]. It is based on the local features of signals, and involves dividing the signals into many narrowband components, called IMFs. The IMFs are distributed according to their frequencies from high to low, and the frequency of each IMF changes with the raw signal. Unlike traditional denoising methods, EMD is applicable for both nonlinear and non-stationary signals [14].
An EMD denoising method based on the cross- correlation coefficient and kurtosis(CCK) criterion was proposed [15-16]. The cross-correlation coefficient represents the correlation between each of IMFs and the raw signal. A greater cross-correlation coefficient indicates a greater relevance between the IMF and the raw signal. Kurtosis, a dimension parameter, is described as “peakedness”. It is very sensitive to the impulse component. The kurtosis value of a healthy rolling element bearing is approximately three. When a bearing has a fault, its kurtosis value increases noticeably. The EMD CCK criterion denoising method involves selecting the IMFs whose cross-correlation coefficient is high and whose kurtosis value is greater than three, according to the CCK criterion, and then reconstructing the signal. The rest of the IMFs, which are considered as noise and false components, are removed.
The wavelet shrinkage method was proposed [17-18], and it is mainly divided into three steps, as follows. 1) Choose the appropriate wavelet basis function and the optimal decomposition level J, and then decompose the signal; 2) Select a threshold for the wavelet coefficients from 1 to J and quantify them; 3) Reconstruct signal using the low frequency coefficient of the Jth level and the quantized wavelet coefficients from 1 to J.
The noise is typically a high frequency signal, and the rolling element bearing vibration signal is also concentrated in front of the IMFs. Some IMFs leaded by the noise also have impulse components in low SNR conditions. Hence, removing the IMFs whose cross- correlation coefficient is high and whose kurtosis value is greater than three directly according to the CCK criterion may cause the serious distortion of the reconstructed signal. The optimal decomposition level is not easy to determine in the process of wavelet shrinkage, and the denoising effect is not ideal if the wavelet basis function is selected inappropriately.
In the present work, the EMD CCK criterion denoising is combined with wavelet shrinkage, and the cross-correlation coefficient threshold value is set as θ=0.1. The IMFs are removed if their cross-correlation coefficient ρxy is less than 0.1, and the IMFs with ρxy≥0.1 and kurtosis value less than three are subjected to soft-thresholding denoising. Lastly, the signal is reconstructed by adding the denoising IMFs to the rest of the IMFs.
The soft-thresholding function of the wavelet shrinkage [17] is given by Eq. (4).
(4)
where sgn() is the sign function; λj represents the threshold value of the jth level, and 
σj
is the standard deviation at the jth level, 
(1≤j≤J); N is the number of wavelet
coefficients at the jth level, and J is the maximum level of the wavelet decomposition.
According to Eq. (4), the soft-thresholding function 
in this work is calculated as follows:
(5)
where λj is the threshold value of imfj, and can be described as
(6)
where N is the length of the IMF component, and J is the number of IMFs.
The advantages of EMD ST denoising can include the following: 1) the adaptability, multi-resolution and completeness of the EMD, 2) the ability to ensure that the IMFs containing impulse components are not removed and ensure the integrity of the signal, and 3) the ability to avoid selecting the wavelet optimal decomposition level, and the wavelet basis function in the process of the wavelet decomposition.
4 Technique process of proposed algorithm
The steps of the algorithm for rolling element bearing IRF estimation by EMD ST denoising and the IFCF for a low SNR and non-stationary signal can be given as follows. A detailed flowchart is presented in Fig. 2.
1) Decompose raw signal by EMD;
2) Calculate cross-correlation coefficient between IMFs and raw signal, and calculate IMFs’ kurtosis values;
3) Select IMFs that must be denoised and removed according to CCK criterion, respectively;
4) Process the selected IMFs by the EMD soft-thresholding denoising method, then reconstruct signal by the rest of the IMFs and denoising IMFs;
5) Signal filtering via fast spectral kurtosis analysis;
6) Obtain the filtered envelope signal by Hilbert transform, then obtain the time-frequency representation (TFR) of the envelope signal by taking a STFT;
7) Extract the IFCFs using an amplitude-sum based peak search algorithm;
8) Determine the IFCF trend using a least squares fitting algorithm;
9) Calculate the IRF by dividing IFCF trend by the FCC.
Fig. 2 Flowchart of algorithm for rolling element bearing IRF estimation by EMD ST denoising and IFCF
5 Simulation analyses and experimental tests
5.1 Simulation analyses
In order to test the proposed algorithm, a fault bearing signal under time-varying speed is simulated, and the simulation formula [13] is as follows:
(7)
where Am is the amplitude of the mth fault impulse, β is the structural attenuation coefficient, wr represents resonance frequency exited by bearing a local fault, u(t) is a unit step function, and tm is the occurrence time of the mth impulse, calculated as follows:
(8)
where u=0.01-0.02 is the error caused by the random slippage of the rolling elements; f(t) is the time-varying rotational frequency, t0=0, and n is the number of the fault impulse of each rotation.
According to these two formulas, the simulation formula for the fault rolling element under time-varying speed in this work is described as follows:
(9)
where Am is proportional to the time; the instantaneous rotation frequency is f(t)=10(t)+30; Noi is noise the resonance frequency is wr=3500; the attenuation coefficient is β=8; the sampling rate is Fs=5000 Hz; the number of sampling points is 20000; Cfc=2; the noise is white Gaussian noise, and the SNR=-8 dB; the simulated signal in time domain is shown in Fig. 3.
The SNR of the simulated time-varying speed fault bearing signal is low, and the fault impulses are entirely covered by the background noise. The IFCF trend, the IRF calculated by the IFCF trend and Cfc=2 with the uncorrected algorithm, and the theoretical IRFs are presented in Fig. 4 respectively. The calculated IRF reflects the changing trend of the theoretical IRF. However, because of the interference of the noise, it has large errors. The largest deviation is 14.77 Hz, and the average deviation is 9.36 Hz. The errors affect the applicability of the algorithm for rolling element bearing IRF estimation based on the IFCF.
Fig. 3 Simulated signal in time domain
Fig. 4 IFCF trend, calculated IRF without denoising and theoretical IRF
The IRF is estimated using the algorithm for rolling element bearing IRF estimation based on EMD ST denoising and the IFCF. The simulated signal is decomposed by EMD, and 13 IMFs are obtained. The cross-correlation coefficient between the IMFs and the simulated signal and the IMFs’ kurtosis value are calculated. Table 1 presents the calculation results for IMF1-IMF7. According to Table 1 and the cross- correlation coefficient threshold θ=0.1, the cross- correlation coefficients of IMF7-IMF13 are less than 0.1, thus these IMFs are removed. For IMF1-IMF6, the kurtosis values of IMF1, IMF3, IMF5 and IMF6 are less than 3, and then applying soft-thresholding denoising to them. The kurtosis values of 
and 
are shown in Table 1, and they increase noticeably compared with the corresponding kurtosis values, highlighting the fault characteristics of the simulated signal. The signal is reconstructed by adding the denoising IMFs into IMF1 and IMF3. Signal filtering is performed via a fast kurtogram. We can apply theHilbert transform to process filter signal and obtain a time-frequency representation of the envelope signal by STFT. The IFCF trend is determined by fitting the IFCFs, which are extracted using amplitude-sum based peak search algorithm from the time-frequency representation of the envelope signal. We can determine the IRF by dividing the IFCF trend by the FCC. The fitted IFCF trend, the calculated IRF and the theoretical IRF are shown in Fig. 5(a). The calculated IRF is consistent with the theoretical IRF, that is f(t)=10(t)+30. The biggest deviation is 2 Hz, and the average deviation is 1.03 Hz.
Table 1 Cross-correlation coefficients and kurtosis values of IMF1-IMF7, kurtosis values of IMF′
Figure 5(b) displays the fitted IFCF trend, the IRF calculated using EMD CCK criterion denoising, and the theoretical IRF. The calculated IRF contradicts the theoretical IRF. The IRFs calculated using the proposed algorithm, no modificatory algorithm and EMD CCK criterion denoising, respectively, are plotted in Fig. 5(c). It is concluded that EMD CCK criterion denoising can cause severe loss, and EMD ST denoising is better than EMD CCK criterion denoising. On the other hand, compared with the uncorrected algorithm the algorithm for rolling element bearing IRF estimation based on EMD ST denoising and the IFCF can reduce the interference of the noise, yielding an estimated IRF that is closer to the theoretical value.
5.2 Experimental verification
The effectiveness of the presented algorithm is further tested using a faulty bearing vibration signal, which is measured on a test rig. The experimental setup of the bearing is shown in Fig. 6. A motor is controlled by a speed regulator. An acceleration sensor is mounted on the top surface of the faulty bearing. A tachometer is mounted in order to measure the shaft rotational speed. The sampling device is a YE6231 acquisition card and its corresponding software. A spark cut groove is used to simulate the outer fault of the rolling element bearing, and the parameters of the bearing in the experiment are listed in Table 2. According to Table 2, the FCC can be calculated as
Fig. 5 (a) IFCF trend and calculated IRF obtained using proposed algorithm, along with theoretical IRF; (b) IFCF trend and calculated IRF obtained using EMD CCK criterion denoising, along with theoretical IRF; (c) IRFs calculated using proposed algorithm, uncorrected algorithm and EMD CCK criterion denoising respectively, along with theoretical IRF
The raw vibration signal in time domain is shown in Fig. 7. The sampling rate is 24000 Hz, and the sampling time is 2 s.
Figure 8 shows the estimation results for the bearing IRF obtained using different algorithms. The IFCF trend, the IRF calculated using the proposed algorithm and the IRF measured by the tachometer are presented inFig. 8(a). The calculated IRF accurately reflects the changing trend of the measured IRF, the biggest deviation is 6.3 Hz, and the average deviation is 2.6 Hz. The IFCF trend, the IRF calculated using EMD CCK criterion denoising, and the measured IRF are shown in Fig. 8(b). Here the calculated IRF generates distortion and has no reference value. The IFCF trend, the IRF calculated using the uncorrected algorithm and themeasured IRF are shown in Fig. 8(c). Compared with the measured IRF, the calculated IRF has obvious errors, the biggest deviation is 46.1 Hz, and the average deviation is 28.05 Hz. The IRFs calculated using three types of algorithms and the measured IRF are shown in Fig. 8(d). Here, it is observed that the EMD ST denoising is better than the EMD CCK criterion denoising, as it avoids distortion of the reconstructed signal. The proposed algorithm can estimate the IRF accurately, and has extensive applications compared with the uncorrected algorithm.
Fig. 6 Experimental setup for rolling element bearing
Table 2 Parameters of bearing in experiment
Fig. 7 Raw vibration signal in time domain
Fig. 8 (a) IFCF trend, IRF calculated using corrected algorithm and measured IRF; (b)IFCF trend, IRF calculated using EMD CCK criterion denoising and measured IRF; (c) IFCF trend, IRF calculated using uncorrected algorithm, and measured IRF; (d) IRFs calculated using corrected algorithm, uncorrected algorithm and EMD CCK criterion denoising respectively, along with measured IRF
6 Conclusions
In order to solve the problem that the algorithm for rolling element bearing IRF estimation based on the IFCF generates errors in a low SNR conditions, a corrected algorithm for rolling element bearing IRF estimation based on EMD ST denoising and the IFCF is proposed in this work. Analyses of both simulated and experimental bearing vibration signals show that the proposed algorithm has good effectiveness.
1) The EMD ST denoising can ensure the integrity of the signal, and highlight fault features in contrast to the EMD CCK criterion denoising.
2) The algorithm for rolling element bearing IRF estimation based on EMD ST denoising and the IFCF can effectively reduce the interference of noise, improving the accuracy of the estimated IRF, and expanding the applications of the method.
References
[1] AMIRAT Y, BENBOUZIDA M E H, AL-AHMAR E, BENSAKER B,TURRI S. A brief status on condition Monitoring and fault diagnosis in wind energy conversion systems [J]. Renewable and Sustainable Energy Reviews, 2009, 13(9): 2629-2636.
[2] CHENG W D, ROBERT X G, WANG J J, WANG T Y，WEN W G, LI J Y. Envelope deformation in computed order tracking and error in order analysis [J]. Mechanical System and Signal Processing, 2014, 48(1/2): 92-102.
[3] FYFE K R, MUNCK E D S. Analysis of computed order tracking [J]. Mechanical System And Signal Processing, 1997, 11(2): 187-205.
[4] COATS M D, RANDALL R B. Order-tracking with and without a tacho signal for gear fault diagnostics [C]// Proceedings of Acoustics, Fremantle, Australia: Australian Acoustical Society, 2012: 447-454.
[5] BONNARDOT F, BADAOUI M E, RANDALL R B, DANIERE J, GUILLET F. Use of the acceleration signal of a gearbox in order to perform angular resampling (with limited speed fluctuation) [J]. Mechanical Systems and Signal Processing, 2005, 19(4): 766-785.
[6] COMBET F, GELMAN L. An automated methodology for performing time synchronous averaging of a gearbox signal without speed sensor [J]. Mechanical Systems and Signal Processing, 2007, 21(6): 2590-2606.
[7] COMBET F, ZIMROZ R. A new method for the estimation of the instantaneous speed relative fluctuation in a vibration signal based on the short time scale transform [J]. Mechanical Systems and Signal, 2009, 23(4): 1382-1397.
[8] URBANEK J, BARSZCZ T, ANTONI J. A two-step procedure for estimation of instantaneous rotational speed with large fluctuations [J]. Mechanical Systems and Signal Processing, 2013, 38(1): 96-102.
[9] BOASHASH B. Interpreting and estimating the instantaneous frequency of a signal-Part 1: Fundamentals [J]. Proceedings of the IEEE, 1992, 80(4): 520-538.
[10] BOASHASH B. Interpreting and estimating the instantaneous frequency of a signal-Part 2: Algorithms and applications [J]. Proceedings of the IEEE, 1992, 80(4): 540-568.
[11] GUO Yu, QIN Shu-ren, TANG Bao-ping, JI Yue-bo. Order tracking of rotating machinery based on instantaneous frequency estimation [J]. Chinese Journal of Mechanical Engineering, 2003, 39(3): 32-36. (in Chinese)
[12] ZHAO Xiao-ping, ZHAO Xiu-li, HOU Rong-tao, HENG Yu-hui, YANG Wen-ping. A new method for instantaneous frequency estimation of run-up or run- down vibration signal for rotating machinery [J]. Chinese Journal of Mechanical Engineering, 2011, 47(7): 103-108. (in Chinese)
[13] WANG T Y, LIANG M, LI J Y, CHENG W D. Rolling element bearing fault diagnosis via fault characteristic order (FCO) analysis [J]. Mechanical System and Signal Processing, 2014, 45(1): 139-153.
[14] HUANG N E, ZENG S, LONG S R, WU M C, SHIH H H, ZHENG Q, YEN N C, TUNG C C, LIU H H. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis [J]. Proc R Soc Lond A, 1998, 454(1971): 903-995.
[15] SU Wen-sheng, WANG Feng-tao, ZHANG Zhi-xin, GUO Zheng- gang, LI Hong-kun. Application of EMD denoising and spectral kurtosis in early fault diagnosis of rolling element bearing [J]. Journal of Vibration and Shock, 2010, 29(3): 18-21. (in Chinese)
[16] TANG Gui-ji, WANG Xiao-long. Fault diagnosis for roller bearings based on EEMD de-noising and 1.5-dimension energy spectrum [J]. Journal of Vibration and Shock, 2014, 33(1): 6-10. (in Chinese)
[17] DAVID L D, IAIN M J. Threshold selection for wavelet shrinkage of noisy data [J]. Engineering in Medicine and Biology Society, IEEE, 1994, 1: A24-A25.
[18] DONOHO D L, JOHNSTONE I. Wavelet shrinkage asymptotic [J]. Journal of Royal Statistical Society, 1995, 57(2): 301-369.
(Edited by FANG Jing-hua)
Foundation item: Project(51275030) supported by the National Natural Science Foundation of China; Project(2016JBM051) supported by the Fundamental Research Funds for the Central Universities, China
Received date: 2015-04-14; Accepted date: 2015-07-28
Corresponding author: CHENG Wei-dong, Associate Professor, PhD; Tel: +86-10-51687004; E-mail: wdcheng@bjtu.edu.cn
