Time-varying demodulation analysis for rolling bearing fault diagnosis under variable speed conditions

Journal of Sound and Vibration 400 (2017) 71–85

Contents lists available at ScienceDirect

Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi

Time-varying demodulation analysis for rolling bearing fault diagnosis under variable speed conditions Zhipeng Feng a,n, Xiaowang Chen a, Tianyang Wang b a b

School of Mechanical Engineering, University of Science and Technology Beijing, Beijing 100083, China Department of Mechanical Engineering, Tsinghua University, Beijing 100084, China

a r t i c l e i n f o Article history: Received 13 November 2016 Received in revised form 9 March 2017 Accepted 27 March 2017 Handling editor: K. Shin Keywords: Rolling bearing Fault diagnosis Nonstationary Synchrosqueezing transform ConceFT

abstract Rolling bearings often work under variable speed conditions, resulting in nonstationary vibrations. How to effectively extract the time-varying fault frequency from nonstationary vibration signals is a key issue in rolling bearing fault diagnosis. To address this issue, a quality time–frequency analysis of excellent time–frequency readability and robust to noise is necessary. To this end, the concentration of frequency and time (ConceFT) method is exploited. Based on this time–frequency analysis method, and considering the modulation feature of rolling bearing vibrations, we propose joint time-varying amplitude and frequency demodulated spectra to reveal the time-varying fault characteristic frequency. Firstly, the optimal frequency band sensitive to rolling bearing fault is selected by spectral kurtosis. Then, both the amplitude envelope and instantaneous frequency of the sensitive signal component within the selected optimal frequency band are calculated. Next, the ConceFT method is applied to the amplitude envelope and instantaneous frequency to generate the time-varying amplitude and frequency demodulated spectra. Finally, rolling bearing fault can be diagnosed by analysis of the time-varying frequency revealed by the time-varying demodulated spectra. This method is free from complex time-varying sidebands, and is robust to noise interference. It is illustrated by numerical simulated signal analysis, and is further validated via lab experimental rolling bearing vibration signal analyses. The localized defects on both inner and outer race are successfully diagnosed. & 2017 Elsevier Ltd All rights reserved.

1. Introduction Rolling bearings are widely used in many sorts of machinery, and they often work under variable speed conditions in practice. In such circumstances, their vibration signals will exhibit strong nonstationarity. Therefore, how to effectively extract fault features from nonstationary signals is an important topic for rolling bearing fault diagnosis [1–3]. This challenging topic has attracted researchers' attention. Most of their works exploit order tracking to remove the effect of speed time variation. To name a few for example, Borghesani et al. [4] designed a reversed sequence squared envelope spectrum and a reversed sequence envelope spectrogram based on computed order tracking and cepstrum pre-whitening, to remove transient condition effects for rolling bearing diagnostics under nonstationary conditions. Zhao et al. [5] devised an envelope order analysis method by exploiting the generalized demodulation and computed order tracking to detect n

Corresponding author. E-mail address: [email protected] (Z. Feng).

http://dx.doi.org/10.1016/j.jsv.2017.03.037 0022-460X/& 2017 Elsevier Ltd All rights reserved.

72

Z. Feng et al. / Journal of Sound and Vibration 400 (2017) 71–85

rolling bearing fault under variable speeds. Wang et al. [6] presented a fault characteristic frequency order analysis to remove the smearing effect due to speed variation, and thereby extract rolling bearing fault features from nonstationary signals. Mishra et al. [7] denoised signals via wavelet transform and then used envelope order tracking to diagnose rolling bearing fault under variable speeds. Order tracking essentially relies on an equi-angular resampling, but the resampling suffers from interpolation error. More importantly, tachometers or encoders are necessary to provide a reference signal for resampling in angular domain, whereas they are not always available in real applications due to cost concerns and design reasons. Time–frequency analysis can reveal the frequency components and their temporal evolution of nonstationary signals. Some researchers have also tried it to extract the time-varying fault features of rolling bearing vibrations. For example, Li et al. [8] used time–frequency ridges to further enhance the time–frequency analysis via generalized synchrosqueezing transform, and applied it to rolling bearing fault feature extraction in both lower and resonance frequency bands under variable speeds. Shi et al. [9] defined a generalized stepwise demodulation based on iterative generalized demodulation, and combined it with synchrosqueezing transform for time-varying fault feature extraction of rolling bearings under timevarying speed conditions. These publications have enriched the literature on and made contributions to rolling bearing fault diagnosis. Rolling bearing fault signatures are usually contaminated or even overwhelmed by interfering noise. Moreover, under variable speeds, the fault features are often manifested by highly time-varying frequency components. In such circumstances, a higher quality time–frequency analysis method (better time–frequency resolution and readability, and robust to strong noise) is needed. However, reported studies in this direction have not been adequately explored, thus the topic under variable speeds still needs further in-depth investigation. Usually, rolling bearing fault diagnosis relies on effective identification of the fault characteristic frequency from vibration signals. When a localized fault occurs to a rolling bearing, the defective point will strike the mating parts during running, resulting in impulses. Under constant speed, such impulses will be generated almost periodically. The repeating frequency of such periodical impulses links to the defective part of a rolling bearing, thus is called fault characteristic frequency. Given the rolling bearing geometry and running speed, the fault characteristic frequency can be calculated following explicit equations [10,11]. In practice, rolling bearing localized fault is usually detected and located by the presence of fault characteristic frequency in vibration signals. Rolling bearing fault vibration signals have both amplitude modulation (AM) and frequency modulation (FM) features. The fault induced impulses excite the resonance of bearing assembly repetitively. In terms of amplitude, the resonance damps out rapidly before the next resonance arises, resulting in the AM feature. In terms of frequency, in one cycle, the resonance exists in the early portion of the cycle and the instantaneous frequency equals approximately the resonance frequency, while in the later portion, the resonance vanishes due to damping and the instantaneous frequency becomes 0 (as illustrated in Fig. 3(a2)). This means the instantaneous frequency changes repetitively, resulting in FM feature. Because of the AM-FM nature of rolling bearing fault vibration signals, it is somewhat difficult to pinpoint the fault characteristic frequency from the Fourier spectrum of raw signals, even under constant speed conditions. In time domain, the relationship between the AM part and the resonance carrier signal is multiplicative. According to the convolution property of Fourier transform, this results in complex spectral characteristics, i.e. intricate sidebands around the resonance frequency. One has to figure out fault characteristic frequency according to the sideband spacing. This leads to difficulty in rolling bearing fault diagnosis via traditional Fourier spectrum analysis. Amplitude demodulation is one of the most prevalent and effective methods for rolling bearing fault diagnosis [12], because it can avoid the complex sidebands around the resonance frequency and can extract the fault characteristic frequency directly. The amplitude modulating frequency of rolling bearing fault vibration signals links to the fault characteristic frequency. Hence, envelope spectrum is developed to pinpoint the modulating frequency and thereby the fault characteristic frequency. However, envelope spectrum still involves some sidebands around the fault characteristic frequency due to the additional AM effect caused by loading zone passing. The sideband spacing is irrelevant to the rolling bearing fault, such as the shaft rotating frequency in the inner race fault case, and the cage rotating frequency in the rolling element fault case. These irrelevant frequencies may mislead fault diagnosis. Frequency demodulation can also reveal the fault characteristic frequency effectively. The instantaneous frequency of fault induced impulse series repeats a cycle every time when the defective point strikes the mating parts. More importantly, the frequency modulating frequency corresponds to the fault characteristic frequency only, independent of loading zone passing effect. Hence, frequency demodulation may reveal the fault characteristic frequency directly, free from irrelevant frequency interferences. Unfortunately, the FM information is usually discarded and has not been fully exploited. Under variable speed conditions, the rolling bearing fault characteristic frequency is not constant, but time-varying, because it is proportional to the running speed and follows its time-varying profile. Therefore, the key to success in rolling bearing fault diagnosis under variable speed conditions lies in effective extraction of the time-varying fault characteristic frequency. Meanwhile, both the amplitude envelope and the instantaneous frequency of rolling bearing fault vibration signals under variable speed conditions will exhibit time-varying features. To reveal the temporal evolution information of constituent components in nonstationary signals, one effective approach is time–frequency analysis. We can treat the amplitude envelope and the instantaneous frequency as signals and extend time–frequency analysis to them, thus resolving the time-varying modulating frequency and thereby identifying the fault characteristic frequency. This inspires the idea of

Z. Feng et al. / Journal of Sound and Vibration 400 (2017) 71–85

73

time-varying amplitude and frequency demodulated spectra. For this purpose, a quality time–frequency analysis method is necessary. Time–frequency analysis is capable of revealing the constituent frequency components of nonstationary signals and their time variation features, thus providing an effective approach to nonstationary signal analysis. To date, various time–frequency analysis methods have been proposed [13]. Among them, linear and bilinear time–frequency representations are widely used. They have advantages but also suffer from drawbacks. For instance, linear time–frequency representations, such as short time Fourier transform and wavelet transform, are subject to Heisenberg uncertainty principle, and therefore suffer from imperfect time–frequency resolution [13,14]. Bilinear time–frequency representations, e.g., Cohen and affine class distributions based on Wigner-Ville distribution, suffer from outer (cross term) interferences when analyzing multicomponent signals, as well as inner interference when analyzing mono-component signals with nonlinearly time-varying instantaneous frequencies [13–16]. To overcome drawbacks of the aforementioned time–frequency representations, several methods have been proposed to improve the readability of time–frequency representations [13]. One of the most promising methods is the synchrosqueezing transform proposed by Daubechies et al. [17]. It is essentially a special case of time–frequency reassignment to enhance time–frequency concentration. By locally ‘squeezing’ time–frequency distribution along frequency axis, it achieves better time–frequency readability, and in the meantime retains invertibility for signal component reconstruction. Nevertheless, the performance of synchrosqueezing transform degrades for more complex signals and under strong noise interferences. To address this issue, Daubechies et al. [18] further proposed the concentration of frequency and time (ConceFT) method to enhance the time–frequency readability via an extended multitaper synchrosqueezing transform recently. ConceFT inherits the merits from both synchrosqueezing and multitapering, thus having better time–frequency readability and stronger robustness to noise interferences. In this sense, it may give a good insight into the time–frequency structure of nonstationary signals, thus providing a better approach to time-varying amplitude and frequency demodulation analysis for rolling bearing fault diagnosis under variable speed conditions. It is important to select an optimal frequency band sensitive to rolling bearing fault for amplitude envelope and instantaneous frequency estimation. Rolling bearing fault is mainly manifested by impulsive components carried by resonance vibration. Since resonance frequency is independent of running speed, fault impulses distributes in a constant frequency band, no matter how the running speed changes over time. Therefore, effective identification of an optimal frequency band sensitive to impulsive components is a key step in rolling bearing fault diagnosis. Spectral kurtosis is capable of revealing the optimal frequency band containing the impulsive transient components in signals [19–22]. This offers a solution to the problem of optimal frequency band selection. In this paper, we propose joint time-varying amplitude and frequency demodulated spectra based on spectral kurtosis and ConceFT method, to extract the time-varying fault characteristic frequency for rolling bearing fault diagnosis under variable speed conditions. We exploit spectral kurtosis to find the optimal frequency band for fault sensitive component separation, and apply ConceFT to the amplitude envelope and instantaneous frequency of the sensitive component. The resultant time-varying amplitude and frequency demodulated spectra are free from complex time-varying sidebands and can directly reveal the time-varying modulating frequency. Therefore, we can easily pinpoint the time-varying fault characteristic frequency and thereby diagnose rolling bearing fault. Hereafter, this paper is structured as follows. In Section 2, we introduce the principle of spectral kurtosis and ConceFT methods, and propose the time-varying amplitude and frequency demodulated spectra. In Section 3, a numerical simulated rolling bearing vibration signal is analyzed to illustrate the principle and to demonstrate the performance of time-varying amplitude and frequency demodulated spectra. In Section 4, the lab experimental vibration signals of rolling bearings with localized fault on inner and outer race are analyzed to further validate the proposed method. Finally, conclusions are drawn in Section 5.

2. Time-varying demodulated spectra Rolling bearing fault information is mainly contained in the amplitude envelope and the instantaneous frequency of vibration signals. Hence, joint amplitude and frequency demodulation analysis is an effective approach to resolve the fault characteristic frequency for fault diagnosis. Under nonstationary conditions, the amplitude envelope and instantaneous frequency also exhibit nonstationarity, therefore we propose time-varying amplitude and frequency demodulated spectra to identify the time-varying fault characteristic frequency. To render effective amplitude and frequency demodulation of rolling bearing vibration, we select an optimal frequency band sensitive to fault with spectral kurtosis method. To better reveal the time evolution feature of amplitude envelope and instantaneous frequency, we conduct time–frequency analysis via ConceFT. 2.1. Spectral kurtosis Kurtosis is a commonly used index for detecting impulses in a signal. The more impulsive the signal is, the larger the kurtosis becomes. However, it is a statistical index in time domain only, thus cannot reveal the distribution of impulse in frequency domain. In order to reveal the distribution of impulsive signal components over frequency domain, spectral

74

Z. Feng et al. / Journal of Sound and Vibration 400 (2017) 71–85

kurtosis (SK) is proposed [19–22]. It is defined as a 4th order normalized cumulant in time–frequency domain

SK(f ) =

X 4 (t , f ) X2(t , f )

2

−2 (1)

where 〈⋅〉 stands for an average operator in time domain, X (t , f ) is the time–frequency envelope of band-pass filtered component of signal x(t ) around frequency f , the constant 2 stems from the fact that X (t , f ) is complex. If the impulses concentrate more in a certain frequency band, the kurtosis of the filtered signal will be higher accordingly. Therefore, SK can be used as an index to find the most impulsive component of the signal in frequency domain. For highly nonstationary signals such as transients, their SK values strongly depend both on frequency and on frequency resolution. There may exist an optimal combination of a frequency f and a frequency resolution Δf which maximizes the SK. This leads to kurtogram, a representation of SK on the frequency–bandwidth/frequency resolution ( f , Δf ) plane [21]. Various algorithms have been proposed to calculate the kurtogram, such as short time Fourier transform, wavelet transform and the fast kurtogram [19–22]. In this paper, the fast kurtogram algorithm [21] is used because of its high computational efficiency. The algorithm is implemented by signal decomposition via an arborescent filter bank constructed based on a quasi-analytic band-pass filter. In different decomposition levels, the signal component is represented based on “binary-1/3 tree”. The signal is decomposed into an n level binary tree and an n − 1 level 1/3 tree separately. For each level k , the decomposition at level k of the 1/3 tree is inserted between levels k and k + 1 of the binary tree. The algorithm is detailed as follows [21]: Firstly, construct the following quasi-analytic low-pass and high-pass filters respectively based on a low-pass filter h(n) with a normalized cut-off frequency fc = 1/8 + ε , ε ≥ 0

h0(n) = h(n)exp( jnπ /4)

(2a)

h1(n) = h(n)exp( j 3nπ /4)

(2b)

where h0(n) and h1(n) represent the low-pass and high-pass filters with normalized frequency pass band of [0,1/4] and [1/4,1/ 2], respectively. Then, iterate the above filter bank construction process in a pyramidal manner to realize the binary-tree like decomposition. In this way, the signal is decomposed into multiple levels. At each level k , the signal is separated into 2k frequency bands. At next level k + 1, the signal decomposition can be iteratively derived as

Ck2+i 1(n) = h0(m)⁎Cki(2n)

(3a)

Ck2+i +11(n) = h1(m)⁎Cki(2n)

(3b)

where the asterisk ⁎ means filtration convolution, the coefficient Cki(n) actually represents the complex envelope of signal

x(n) in the pass band with a center frequency f cSK = (i + 2−1)2−k − 1 and a bandwidth Δfk = 2−k − 1, and Cki(n) is the correi sponding envelope signal. Next, extend the dyadic filter bank to a 1/3-binary tree structure, in order for a higher frequency resolution. It is constructed using three additional quasi-analytic band-pass filters g0(n), g1(n) and g2(n) with normalized frequency pass-bands of [0, 1/6], [1/6, 1/3] and [1/3, 1/2] respectively. As such, the frequency–frequency resolution plane can be sampled more finely via the 1/3-binary tree filter bank. Decompose further the dyadic binary signal decomposition coefficient Cki(n) via the 1/3-binary tree filter bank, obtaining three sub-sequence Cki+.6j(n) ( j = 0, 1, 2) corresponding to the low-, medium-, and highfrequency parts of the frequency interval [i⋅2−k − 1, (i + 1)⋅2−k − 1] respectively. The subscript k.6 means that there are

3 × 2k ≈ 2k + 0.6 sub-sequences Ckl.6(n) inserted between levels k and k + 1 of the binary tree. Add the 1/3-binary tree filter bank to the existing levels to obtain a higher frequency resolution decomposition. Finally, calculate SK f cSK and a , Δfk of each filtered signal envelope in the pass band with a center frequency f cSK i i bandwidth Δfk

(

SK

(

f cSK , i

)

Δfk =

)

Cki(n) Cki(n)

4

2 2

−2 (4)

Rolling element bearing fault induced vibration usually features impulses. The more impulsive the fault induced vibration, the larger the SK. On the frequency–frequency resolution plane, the locations with higher SK values indicate which frequency band the impulses most likely concentrate on. Therefore, according to the center frequency and bandwidth pertaining to the larger SK, an optimal filter can be designed, and the impulsive component sensitive to bearing fault can be filtered out from the raw signal.

Z. Feng et al. / Journal of Sound and Vibration 400 (2017) 71–85

75

2.2. Concentration of frequency and time Synchrosqueezing transform can effectively improve the time–frequency readability, but it suffers from noise and its quality degrades as noise level increases. To improve the robustness to noise, Daubechies et al. [18] extended multitaper to synchrosqueezing transform, and further proposed the ConceF T method. 2.2.1. Synchrosqueezing transform Recently, Daubechies et al. [17] proposed the synchrosqueezing transform to enhance time–frequency analysis. It eliminates the blur in scale/frequency domain by squeezing the time–frequency representation in a proper scale/frequency range, thus achieving better time–frequency readability. Here, we introduce briefly the synchrosqueezing transform based on the wavelet transform. The wavelet transform of an arbitrary signal x(t ) can be written as

1 a

WT(a, b) =

+∞

∫−∞

⎛ t − b⎞ x(t )ψ ⁎⎜ ⎟dt ⎝ a ⎠

(5)

where ψ (t ) is the reference wavelet, a is the scale, b is the time offset, and * stands for the complex conjugate. Subjected to the Heisenberg uncertainty principle, the signal energy distribution spreads over a vicinity around the true instantaneous frequency on time–scale/frequency plane, causing the ‘blur’ phenomenon. The smearing distribute only along scale dimension a , while its oscillatory behavior in time dimension b points to the true frequency, regardless of the scale value a . Thus the instantaneous frequency at any time–scale location (a, b) can be calculated as

ω(a, b) =

−j ∂WT(a, b) ∂b WT(a, b)

(6)

where ∂(⋅) stands for partial derivative and WT(a, b) ≠ 0. To suppress the blur along the scale, the spread energy is reallocated to the instantaneous frequency, by summing the contribution within successive bins [ωl − (1/2)Δω, ωl + (1/2)Δω]. Then the synchrosqueezing transform is achieved as

ST(ωl , b) =

1 ∑ WT(ak, b)ak−3/2(Δa)k Δω a k

(7)

where ak is the kth discrete scale point satisfying ω(ak , b) − ωl ≤ Δω/2, ωl is the lth center frequency, Δω = ωl − ωl − 1 and (Δa)k = ak − ak − 1. The synchrosqueezing transform condenses the blur of time–frequency distribution along the scale/frequency dimension into a concentrated region, thus considerably improving the time–frequency readability. However, its performance degrades in the presence of noise, because the spurious noise distribution obscures the concentration of true signal components. This can be reduced by multitapering. 2.2.2. Multitaper synchrosqueezing transform In the time–frequency representation based on synchrosqueezed continuous wavelet transform, the dominant component always has similar time–frequency structure, even under different reference wavelets [21]. However, the noise distribution varies from one reference wavelet to another, because wavelet transform is essentially a convolution with a reference wavelet, so that the wavelet transforms of independent and identically distributed noise based on different orthonormal reference wavelets are independent. This property leads to the idea of a multitaper synchrosqueezing transform algorithm. For a signal x(t ), given I orthonormal reference wavelets ψi(t ), i = 1, ⋯ , I , we can calculate the instantaneous frequency ω(ψi)(a, b) and the synchrosqueezing transform ST (ψi)(ω, b) under each wavelet ψi(t ), according to Eqs. (6) and (7). Then the multitaper synchrosqueezing transform is defined as the average of the individual ST (ψi)(ω, b)

MST(ω, b) =

1 I

I

∑ ST (ψi)(ω, b) i=1

(8)

Averaging over a large number of orthonormal reference wavelets cancels out the time–frequency distribution artifacts caused by the noise. It seems that the larger the number of orthonormal reference wavelets, the better the noise artifacts suppressed. However, the smeared area of the signal information on the time–frequency plane also increases with the number of orthonormal reference wavelets. A balance needs to be set, and this limits how many different orthonormal reference wavelets can be used. 2.2.3. Concentration of frequency and time To overcome the limitations of multitaper synchrosqueezing transform, Daubechies et al. [18] recently generalized the mutlitaper method and proposed a ConceFT method, by exploiting the nonlinearity of synchrosqueezing transform. I For each linear combination of orthonormal reference wavelets ψ (t ) = ∑i = 1 riψi(t ) , where the weight ri is real, the

76

Z. Feng et al. / Journal of Sound and Vibration 400 (2017) 71–85

synchrosqueezing transform can be applied, but the result is not a linear combination of the synchrosqueezing under each I orthonormal wavelet ψi(t ), i.e. ST (ψ )(ω, b) ≠ ∑i = 1 ST (ψi)(ω, b), because of the nonlinearity of synchrosqueezing operation. Therefore, for different weight vector r = (r1, ⋯ , rI ) , the time–frequency distribution artifacts due to noise appear at sufficiently different locations; hence, averaging over many choices of weight vector r successfully suppresses noise artifacts. This inspires the idea of ConceFT algorithm. Its details based on wavelet transform are as follows: (1) Select I orthonormal reference wavelets ψi(t ), i = 1, ⋯ , I , with good concentration on the time–frequency plane. (2) Pick N random vectors of unit norm rn , n = 1, ⋯ , N . I (3) For each n = 1, ⋯ , N , define ψn(t ) = ∑i = 1 rniψi(t ) , and calculate the corresponding wavelet transform WT (ψn)(a, b)= I (ψi ) ∑i = 1 rni WT (a, b) . (4) Calculate the instantaneous frequency ω(ψn)(a, b) and the corresponding synchrosqueezing transform ST (ψn)(ωl , b). (5) Averaging over the random vectors rn , yields the ConceFT representation

CFT(ω, t ) =

1 N

N

∑ ST (ψn)(ω, t )

(9)

n= 1

In practice, to reduce the smeared area on the time–frequency plane, the number of orthonormal reference wavelets can be set as small as 2, while to average out the noise artifacts, the number of weight vectors N can be chosen as large as one wishes. By this ConceFT method, the spurious distribution due to noise is effectively suppressed, and meanwhile good time– frequency readability is retained. 2.3. Analysis procedure Rolling bearing faults are manifested by impulse series in vibration signals, and such impulse series are carried by the resonance. However, they are usually of low amplitudes and buried in high level background noise in practice. Hence, it is necessary to find an optimal frequency band so as to separate the impulse series sensitive to bearing fault via filtering. Even after an optimal filtration, the separated component is still contaminated by the noise within the filter pass band. Moreover, inevitable estimation error also exists in the amplitude envelope and instantaneous frequency, acting as interfering noise. Meanwhile, under nonstationary conditions, rolling bearing vibration signals have strong nonstationarity. Particularly under variable speed conditions, the fault characteristic frequency follows the speed variation profile. Therefore, an appropriate time–frequency analysis of good time–frequency readability and robust to noise is needed to resolve the time-varying fault characteristic frequency. In order to effectively identify the time-varying fault characteristic frequency from noisy signals, we propose joint timevarying amplitude and frequency demodulated spectra based on spectral kurtosis and ConceFT methods, by exploiting respectively their capability of identifying fault sensitive frequency band, as well as good time–frequency readability and robustness to interfering noises. For a real rolling bearing vibration signal x(t ), the time-varying amplitude and frequency demodulated spectra can be generated following the procedure listed below: (1) Calculate the fast kurtogram of the raw vibration signal x(t ). The peaks with high SK value in the kurtogram most likely indicate the concentration frequency band of bearing fault induced vibration impulses. Identify the center frequency and bandwidth of the highest (or one of the first a few highest) SK peak on the frequency–frequency resolution plane for further analysis. (2) Design a band pass filter based on the center frequency and bandwidth pertaining to the selected SK peak in the kurtogram. Separate the sensitive component via the designed filter. (3.1) Calculate the amplitude envelope of band-pass filtered signal x^(t ) via the Hilbert transform

a(t ) =

{

1 2

}

2 x^ (t ) + H2[x^(t )]

(10)

where H (⋅) stands for Hilbert transform. (3.2) Envelope normalize the band-pass filtered signal x^(t )

x¯ (t ) =

x^(t ) a(t )

(11)

to make it a pure FM one, thus removing the effect of AM on instantaneous frequency estimation and fulfilling the requirement of the Bedrosian and Nuttall theorem [23,24]. Estimate the instantaneous frequency of envelope normalized signal x¯ (t ) via Hilbert transform

f (t ) =

⎧ H[x¯ (t )] ⎫ 1 d ⎬ arctan⎨ ⎩ x¯ (t ) ⎭ 2π d t

(12)

(4) Apply ConceFT method to the amplitude envelope a(t ) and the instantaneous frequency f (t ) respectively, to generate the time-varying amplitude and frequency demodulated spectra.

Z. Feng et al. / Journal of Sound and Vibration 400 (2017) 71–85

77

Since rolling bearing fault information is mainly contained by the amplitude envelope of vibration signals, the time variation features of fault characteristic frequencies can be extracted from the time-varying amplitude demodulated spectrum. More importantly, as the amplitude envelope alone does not involve multiplication with the carrier signal (i.e. resonance vibration), the derived time-varying amplitude demodulated spectrum is free from complex sidebands around the resonance frequency, thus enabling fault feature extraction in a more explicit way than the time–frequency representation of raw signals. However, it is still subject to additional AM effect due to loading zone passing. As a result, sidebands may appear around the fault characteristic frequency with a spacing equal to the shaft rotating frequency (in the inner race fault case) or the cage rotating frequency (in the rolling element fault case). Rolling bearing fault signature is also carried by the instantaneous frequency of vibration signals, and the fault characteristic frequency can be revealed from the time-varying frequency demodulated spectrum too. More importantly, as the loading zone passing does not have any effect on the FM part, the derived time-varying frequency demodulated spectrum is free from irrelevant frequencies, and is dominated by the fault characteristic frequency and its harmonics only, thus enabling fault feature extraction in an even better way than the time–frequency representations of both raw signals and their amplitude envelope.

3. Simulated signal analysis In this section, to illustrate the performance of proposed method, a numerical simulated rolling bearing vibration signal with a localized outer race fault under variable speed is generated and analyzed. Considering the fast-decaying fault excited bearing resonance, rolling element random slippage and speed time variation, rolling bearing vibration signals under variable speed conditions can be modeled as [6] M

∑

x(t ) =

A m exp{ − βωn[t − υm(t )]} sin {ωn[t − υm(t )] + ϕ}u[t − υm(t )] + n(t )

m =−M

(13)

m

υm(t ) =

∑

Ti + τi

i =−M

(14)

where m denotes the impulse index, Am the amplitude of mth impulse, β the damping characteristic factor, ωn the bearing resonance frequency, Ti the instantaneous repeating period of the ith impulse at occurrence time ti and equal to the reciprocal of the fault characteristic frequency at time ti , τi represents the random slippage effect of rolling elements, a uniformly distributed random variable with a zero mean and a standard deviation of 0.01Ti~0.02Ti , u(⋅) a unit step function, and n(t ) is a white Gaussian noise at a signal-to-noise ratio of 0 dB to mimic background interferences. Assume the bearing outer race is fixed, the inner race rotates with a shaft at a time-varying frequency f (t ) = 20 + 10 sin(πt ), the outer race fault characteristic frequency fouter (t ) = 3.7f (t ) (so that T (t ) = 1/fouter (t )), the damping characteristic factor β = 0.15, the bearing resonance frequency ωn ¼2048 Hz, and the amplitude of each impulse Am = f [υm(t )] /20 for simulating the time variation of impulse magnitude. Suppose P fault impulses are generated per shaft revolution (here P = 3.7), the beginning time of the vibration signal is t0 , then the occurrence time of the ith impulse can be recursively calculated as

ti = τi +

1 f (ti − 1)P

(15)

Under such settings, a signal is generated and is sampled at a frequency 8000 Hz. Fig. 1(a) and (b) shows its waveform and Fourier spectrum. Fig. 1(d) shows the Morlet wavelet transform scalogram of the signal. The major signal energy distributes in the vicinity of the bearing resonance frequency 2048 Hz. Meanwhile, some time-varying sidebands appear in a region of 1500–2500 Hz. However, it is difficult to identify the sideband spacing and thereby the time-varying fault characteristic frequency, because of obstacles due to the time variation of sidebands and the strong background noise interferences. To better resolve the fault characteristic frequency, the signal is analyzed via the proposed time-varying amplitude demodulation method. Firstly, spectral kurtosis is used to find the optimal frequency band for demodulation analysis. Fig. 2 (a) shows the kurtogram, where higher SK value corresponds to a frequency band 1500–2500 Hz, the same as marked in Fig. 1(d). Then, the sensitive component within the optimal frequency band is separated via a band pass filter. Fig. 2 (b) shows its amplitude envelope. Applying the ConceFT method to the amplitude envelope, yields the time-varying amplitude demodulated spectrum, as displayed by Fig. 2(c1) and (c2). It clearly reveals the outer race fault characteristic frequency and its harmonics up to the third order, despite of their time-variation and the strong noise interference. This finding confirms the outer race fault. For comparison study, the time-varying amplitude demodulated spectrum is also generated via other time–frequency analysis methods, i.e. Morlet wavelet transform, synchrosqueezing transform, and multitaper synchrosqueezing transform. Fig. 2(d)–(f) present the results respectively. Although they also roughly reveal the time-varying fault feature, they suffer from strong noise interferences. For example, in the Morlet wavelet transform scalogram, Fig. 2(d), artifacts due to noise distribute almost everywhere on the time–frequency plane. Although the synchrosqueezing transform improves the time–

78

Z. Feng et al. / Journal of Sound and Vibration 400 (2017) 71–85

2.5

0.07

2

0.06 0.05

1

Amplitude

Amplitude

1.5

0.5 0 −0.5

0.04 0.03 0.02

−1 0.01

−1.5 −2 0

0.2

0.4

0.6

0.8

0 0

1

500

1000

1500

2000

2500

3000

3500

4000

Frequency [Hz]

Time [s] 35 4 3.5 3

Frequency [kHz]

Speed [Hz]

30

25

20

2.5 2 1.5 1 0.5

15 0

0.2

0.4

0.6

Time [s]

0.8

1

0

0

200

400

600

800

1000

Time [s]

Fig. 1. Simulated signal: (a) waveform, (b) Fourier spectrum, (c) speed, and (d) Morlet wavelet transform scalogram.

frequency concentration of the true components, the noise artifacts still exist and even intersect with the instantaneous frequency skeletons of the true components, as shown in Fig. 2(e). As can be seen from Fig. 2(f), the multitaper synchrosqueezing transform reduces the noise effect to some extent, but small artifacts due to noise still exist, causing the energy dilution from the true time–frequency ridges, especially at the boundaries. Moreover, their time–frequency concentration is not as fine as the ConceFT. These findings demonstrate the effectiveness of the proposed method, and the advantages of ConceFT over other conventional methods, i.e. fine time–frequency concentration and robustness to noise. To further reveal the time-varying fault feature, the signal is also analyzed via the proposed time-varying frequency demodulation method. After filtering based on the optimal frequency band found by the kurtogram, the separated sensitive component is envelope normalized, and then its instantaneous frequency is calculated, as shown in Fig. 3(a). We treat the instantaneous frequency as a signal, and analyze it with ConceFT, yielding the time-varying frequency demodulated spectrum, as illustrated in Fig. 3(b1) and (b2). Once again, the outer race fault characteristic frequency and its twice are clearly displayed, thanks to the good time–frequency readability and robustness to noise of ConceFT. This validates our proposition to extract fault feature via frequency demodulation. For comparison, the instantaneous frequency is also analyzed with Morlet wavelet transform, synchrosqueezing wavelet transform, and multitaper synchrosqueezing wavelet transform. Fig. 3(c)–(e) show the results. The time-varying fault characteristic frequency is roughly identified, but in a bad time–frequency readability due to the imperfect time–frequency resolution and the noise artifacts. This again demonstrates the advantages of ConceFT over traditional time–frequency analysis methods.

4. Lab experimental evaluation In this section, the proposed method is further validated via lab experimental dataset analysis. 4.1. Experimental settings The tests on rolling bearings are carried out on a SpectraQuest machinery fault simulator at the University of Ottawa lab. The test rig consists of a shaft supported by two rolling bearings, as shown in Fig. 4. Three groups of tests, i.e., baseline, inner race fault and outer race fault cases, are designed. In the baseline case, ER16K bearings are used and both bearings are healthy. In the inner race fault case, ER16K bearings are utilized while an inner race localized fault is introduced to the left one. In the outer race fault case, ER10K bearings are used, and an outer race localized

Z. Feng et al. / Journal of Sound and Vibration 400 (2017) 71–85

79

1.5 0 1 1.6

1

Amplitude

2 2.6 3

0.5

3.6 4 4.6 5 0

500

1000

1500

2000

2500

3000

3500

0 0

4000

0.2

0.4

0.6

0.8

1

Time [s]

400

true IF trajectories

350

Frequency [Hz]

300

3f outer

250 200 2f outer

150 100

f outer

50 0 0

0.2

0.4

0.6

0.8

1

400

400

300

300

Frequency [Hz]

Frequency [Hz]

Time [s]

200 100 0 0

0.2

0.4

0.6

0.8

200 100 0 0

1

0.2

0.4

0.6

0.8

1

Time [s]

Time [s]

Frequency [Hz]

400 300 200 100 0 0

0.2

0.4

0.6

0.8

1

Time [s] Fig. 2. Time-varying amplitude demodulated spectrum: (a) kurtogram, (b) amplitude envelope, (c1) ConceFT 2D, (c2) ConceFT 3D, (d) Morlet wavelet transform scalogram, (e) synchrosqueezing wavelet transform, and (f) multitaper synchrosqueezing wavelet transform.

fault is introduced to the left one. Table 1 lists the configuration parameters of both bearing types. Table 2 lists the relationship between the fault characteristic frequency and the shaft rotating frequency. Given the shaft rotating frequency frotation (t ) at any time instant t , the fault characteristic frequencies of both types can be calculated accordingly. An accelerometer is mounted on top of the left bearing housing to collect vibration signals. In the baseline and inner race fault cases, the shaft accelerates from 1000 rpm to 1500 rpm within 6 s, and the vibration signals are collected at 12,000 Hz.

80

Z. Feng et al. / Journal of Sound and Vibration 400 (2017) 71–85

f

n

2000 1500 1000 500 0 0

f

2500

Frequency [Hz]

Frequency [Hz]

2500

n

2000 1500 1000 500

0.2

0.4

0.6

0.8

0 0.8

1

0.81

0.82

400

0.83

0.84

Time [s]

Time [s]

true IF trajectories

350

Frequency [Hz]

300 250 200 150 100 50 0 0

0.2

0.4

0.6

0.8

1

Time [s] 400 350

400

Frequency [Hz]

Frequency [Hz]

300

300 200

250 200 150 100

100

50

0 0

0.2

0.4

0.6

0.8

0 0

1

0.2

0.4

0.6

0.8

1

Time [s]

Time [s] 400 350

Frequency [Hz]

300 250 200 150 100 50 0 0

0.2

0.4

0.6

0.8

1

Time [s] Fig. 3. Time-varying frequency demodulated spectrum: (a1) instantaneous frequency overall, (a2) instantaneous frequency close-up, (b1) ConceF T 2D, (b2) ConceFT 3D, (c) Morlet wavelet transform scalogram, (d) synchrosqueezing wavelet transform, and (e) multitaper synchrosqueezing wavelet transform.

In the outer race fault case, the shaft accelerates from 1000 rpm to 2200 rpm within 6 s, and the vibration signal is collected at 24,000 Hz. The speed is measured with a tachometer (1 pulse per cycle) at different time resolution for each case, i.e. 0.150 s, 0.013 s, and 0.023 s in the baseline, inner race fault and outer race fault case respectively.

Z. Feng et al. / Journal of Sound and Vibration 400 (2017) 71–85

81

Fig. 4. Rolling bearing test rig.

Table 1 Parameters of bearings. Bearing type

Number of rolling elements

Pitch diameter (mm)

Rolling element diameter (mm)

ER16K ER10K

9 8

38.52 33.50

7.94 7.94

Table 2 Fault characteristic frequencies of bearings. Bearing type

Inner race

Outer race

Rolling element

ER16K ER10K

5.43frotation(t) 4.95 frotation(t)

3.57 frotation(t) 3.05 frotation(t)

2.32 frotation(t) 1.99 frotation(t)

Note: frotation stands for the shaft rotating frequency.

4.2. Baseline signal analysis In the baseline case, both bearings are healthy. Fig. 5 shows the analysis results. No significant sidebands are found in the time–frequency distribution obtained from Morlet wavelet transform, Fig. 5(c). Based on the higher SK value in the kurtogram, Fig. 5(d), an optimal frequency band of 4500–6000 Hz is selected for further analysis. Fig. 5(e) and (f) show the time-varying amplitude and frequency demodulated spectrum respectively, which are generated following the procedure in Section 2.3. Ripples and speckles distribute evenly over the entire time–frequency plane, but no continuous time–frequency ridges correspond to the fault characteristic frequency of any bearing component. These results imply that both bearings are healthy. 4.3. Detection of outer race fault Fig. 6 shows the signal analysis results of the outer race fault case. In the time–frequency distribution, Fig. 6(c), the signal energy mainly distributes within the marked region 0–2500 Hz. Time-varying sidebands with a spacing pertaining to the outer race fault characteristic frequency is expected within the region, but they are not found due to the insufficient time– frequency resolution of wavelet transform. According to the region associated to the higher SK value in the kurtogram, Fig. 6 (d), 1500–2200 Hz is selected as the optimal frequency band. Then the sensitive component within this frequency band is separated via a band pass filter. By applying the ConceFT to the amplitude envelope of the sensitive component, the corresponding time-varying amplitude demodulated spectrum is generated, as shown by Fig. 6(e). On the time–frequency plane, four significant time-varying frequency components appear, following the shaft speed changing profile. Further analysis confirms that they correspond to the first four harmonics of outer race fault characteristic frequency kfouter (k = 1, ⋯ , 4 ). Applying the ConceFT to the instantaneous frequency of the sensitive component, yields the time-varying

82

Z. Feng et al. / Journal of Sound and Vibration 400 (2017) 71–85

1600

0.15

1500

0.1

Speed [rpm ]

2

Amplitude [m/s ]

1400 0.05

0

1300 1200 1100

−0.05

1000 −0.1 0

1

2

3

4

5

900 0

6

Time [s]

1

2

3

4

5

6

Time [s]

6 5 4 3 2 1

0

1

2

3

4

5

6

400

400

350

350

300

300

Frequency [Hz]

Frequency [Hz]

0

250 200 150

250 200 150

100

100

50

50

0 0

1

2

3

Time [s]

4

5

6

0 0

1

2

3

4

5

6

Time [s]

Fig. 5. Baseline signal: (a) waveform, (b) speed, (c) Morlet wavelet transform scalogram, (d) kurtogram, (e) time-varying amplitude demodulated spectrum, and (f) time-varying frequency demodulated spectrum.

frequency demodulated spectrum, as displayed by Fig. 6(f). The dominant time-varying frequency components also follow the time variation profile of the shaft speed, corresponding to the rotating frequency frotation and the outer race fault characteristic frequency fouter respectively. These findings indicate the existence of outer race fault, consistent with the actual experimental settings. 4.4. Detection of inner race fault Fig. 7 presents the inner race fault signal analysis results. In the time–frequency distribution of the raw signal, Fig. 7(c), some impulsive components appear in the frequency band 3000–6000 Hz, as marked by the dashed rectangle. However, it is difficult to extract the fault feature, because of the limited time–frequency resolution of wavelet transform and particularly the irregularity and time-variation of the fault characteristics. Kurtogram is calculated to assist selecting the optimal frequency band for demodulation analysis, as illustrated by Fig. 7(d). According to the highest SK value, 4500–6000 Hz is adopted as the sensitive band, and the associated sensitive component is separated via band pass filter. Fig. 7(e) shows the time-varying amplitude demodulated spectrum. The amplitude envelope is dominated by several prominent time-varying frequency components, which are the shaft rotating frequency frotation , as well as the inner race fault characteristic frequency

Z. Feng et al. / Journal of Sound and Vibration 400 (2017) 71–85

83

Fig. 6. Outer race fault signal: (a) waveform, (b) speed, (c) Morlet wavelet transform scalogram, (d) kurtogram, (e) time-varying amplitude demodulated spectrum, and (f) time-varying frequency demodulated spectrum.

finner , and its harmonics 2finner and 3finner . Besides, prominent time-varying sidebands are also present around the inner race fault characteristic frequency and its harmonics with a spacing equal to the shaft rotating frequency frotation , such as

finner + frotation and 2finner + frotation . The inner race rotates with the shaft, and passes the loading zone at the shaft rotating frequency. This leads to an additional amplitude modulation effect on the inner race fault impulses, and thereby resulting in the sidebands in the amplitude demodulated spectrum. Fig. 7(f) shows the time-varying frequency demodulated spectrum. The shaft rotating frequency and its twice exist, but the inner race fault characteristic frequency also appear with a prominent magnitude. These features imply the inner race fault, again in agreement with the real experimental settings.

5. Conclusions In order to effective extract time-varying fault features from nonstationary signals of rolling bearings under variable speed conditions, joint time-varying amplitude and frequency demodulated spectra are proposed. They are implemented via spectral kurtosis and ConceFT methods. The optimal frequency band sensitive to rolling bearing fault impulses can be found according to higher spectral kurtosis value in the kurtogram. Rolling bearing fault is mainly contained in the amplitude envelope and the instantaneous frequency, and the fault characteristic frequency links to the modulating frequency

84

Z. Feng et al. / Journal of Sound and Vibration 400 (2017) 71–85

Fig. 7. Inner race fault signal: (a) waveform, (b) speed, (c) Morlet wavelet transform scalogram, (d) kurtogram, (e) time-varying amplitude demodulated spectrum, and (f) time-varying frequency demodulated spectrum.

of both the AM and the FM parts. To finely extract the highly time-varying fault characteristic frequency, ConceFT is applied to the amplitude envelope and the instantaneous frequency of separated sensitive component, thus obtaining time-varying amplitude and frequency demodulated spectra. They have a higher quality time–frequency representation, owing to the merits of ConceFT, i.e. better time–frequency readability and robustness to interfering noise, and more importantly, it is free from complex time-varying sidebands inherent with the time–frequency representation of raw signals. This enables effective fault feature extraction of rolling bearings under variable speed conditions. The proposed method has been demonstrated with numerical simulated signal, and further validated with lab experimental bearing vibration signals. Localized faults on both inner race and outer race have been successfully diagnosed.

Acknowledgements This work is supported by the National Natural Science Foundation of China (51475038), and the Natural Sciences and Engineering Research Council of Canada (RGPIN 121433-2011).

Z. Feng et al. / Journal of Sound and Vibration 400 (2017) 71–85

85

References [1] R.B. Randall, J. Antoni, Rolling element bearing diagnostics-a tutorial, Mech. Syst. Signal Process. 25 (2) (2011) 485–520. [2] A. Jardine, D. Lin, D. Banjevic, A review on machinery diagnostics and prognostics implementing condition-based maintenance, Mech. Syst. Signal Process. 20 (2006) 1483–1510. [3] J. Liu, Shannon wavelet spectrum analysis on truncated vibration signal for machine incipient fault detection, Meas. Sci. Technol. 23 (2012) 1–11. [4] P. Borghesani, R. Ricci, S. Chatterton, P. Pennacchi, A new procedure for using envelope analysis for rolling element bearing diagnostics in variable operating conditions, Mech. Syst. Signal Process. 38 (1) (2013) 23–35. [5] M. Zhao, J. Lin, X. Xu, Y. Lei, Tacholess envelope order analysis and its application to fault detection of rolling element bearing with varying speeds, Sensors 13 (2013) 10856–10875. [6] T. Wang, M. Liang, J. Li, W. Cheng, Rolling element bearing fault diagnosis via fault characteristic order (FCO) analysis, Mech. Syst. Signal Process. 45 (2014) 139–153. [7] C. Mishra, A.K. Samantaray, G. Chakraborty, Rolling element bearing defect diagnosis under variable speed operation through angle synchronous averaging of wavelet de-noised estimate, Mech. Syst. Signal Process. 72–73 (2016) 206–222. [8] J. Shi, M. Liang, D. Necsulescu, Y. Guan, Generalized stepwise demodulation transform and synchrosqueezing for time–frequency analysis and bearing fault diagnosis, J. Sound Vib. 368 (2016) 202–222. [9] C. Li, V. Sanchez, G. Zurita, M.C. Lozada, D. Cabrera, Rolling element bearing defect detection using the generalized synchrosqueezing transform guided by time–frequency ridge enhancement, ISA Trans. 60 (2016) 274–284. [10] P.D. McFadden, J.D. Smith, Model for the vibration produced by a single point defect in a rolling element bearing, J. Sound Vib. 96 (1) (1984) 69–82. [11] P.D. McFadden, J.D. Smith, The vibration produced by multiple point defect in a rolling element bearing, J. Sound Vib. 98 (2) (1985) 263–273. [12] P.D. McFadden, J.D. Smith, Vibration monitoring of rolling element bearings by the high frequency resonance technique-a review, Tribol. Int. 17 (1) (1984) 3–10. [13] Z. Feng, M. Liang, F. Chu, Recent advances in time–frequency analysis methods for machinery fault diagnosis: a review with application examples, Mech. Syst. Signal Process. 38 (1) (2013) 165–205. [14] F. Hlawatsch, G.F. Boudreaux-Bartels, Linear and quadratic time–frequency signal representations, IEEE Signal Process. Mag. 9 (2) (1992) 21–67. [15] L. Cohen, Time–frequency distributions: a review, Proc. IEEE 77 (1989) 941–981. [16] O. Rioul, P. Flandrin, Time–scale energy distributions: a general class extending wavelet transforms, IEEE Trans. Signal Process. 40 (7) (1992) 1746–1757. [17] I. Daubechies, J.F. Lu, H.T. Wu, Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool, Appl. Comput. Harmon. Anal. 30 (2011) 243–261. [18] I. Daubechies, Y. Wamg, H.T. Wu, ConceFT: concentration of frequency and time via a multitapered synchrosqueezed transform, Philos. Trans. R. Soc. A-Math. Phys. Eng. Sci. 374 (2016) 20150193. [19] J. Antoni, The spectral kurtosis: a useful tool for characterizing non-stationary signals, Mech. Syst. Signal Process. 20 (2006) 282–307. [20] J. Antoni, R.B. Randall, The spectral kurtosis: application to the vibratory surveillance and diagnostics of rotating machines, Mech. Syst. Signal Process. 20 (2006) 308–331. [21] J. Antoni, Fast computation of the kurtogram for the detection of transient faults, Mech. Syst. Signal Process. 20 (2007) 108–124. [22] Y. Wang, J. Xiang, R. Markert, M. Liang, Spectral kurtosis for fault detection, diagnosis and prognostics of rotating machines: a review with applications, Mech. Syst. Signal Process. 66–67 (2016) 679–698. [23] A.H. Nuttall, E. Bedrosian, On the quadrature approximation on the Hilbert transform of modulated signals, Proc. IEEE 54 (1966) 1458–1459. [24] N.E. Huang, Z. Wu, S.R. Long, et al., On instantaneous frequency, Adv. Adapt. Data Anal. 1 (2) (2009) 177–229.