Early fault feature extraction of rolling bearing based on ICD and tunable Q-factor wavelet transform

Mechanical Systems and Signal Processing 86 (2017) 204–223

Contents lists available at ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/ymssp

Early fault feature extraction of rolling bearing based on ICD and tunable Q-factor wavelet transform

cross

⁎

Yongbo Lia,b, Xihui Liangb, Minqiang Xua, , Wenhu Huanga a

Department of Astronautical Science and Mechanics, Harbin Institute of Technology (HIT), No. 92 West Dazhi Street, Harbin 150001, People's Republic of China Departments of Mechanical Engineering, University of Alberta, Edmonton, Alberta, Canada T6G 2G8

b

A R T I C L E I N F O

ABSTRACT

Keywords: Signal decomposition Intrinsic characteristic-scale decomposition (ICD) Tunable Q-factor Wavelet transform (TQWT) Fault detection

When a fault occurs on bearings, the measured bearing fault signals contain both high Q-factor oscillation component and low Q-factor periodic impact component. TQWT is the improvement of the traditional single Q-factor wavelet transform, which is very suitable for separating the low Q-factor component from the high Q-factor component. However, the accuracy of its decomposition heavily depended on the selection of Q-factors. There is no reported simple but effective method to select the Q-factors with enough accuracy. This study aims to develop a strategy to diagnostic the early fault of rolling bearings. In this paper, a characteristic frequency ratio (CFR) is used to optimize Q-factors of TQWT (OTQWT). However, directly application of OTQWT is difficult to extract fault signatures at early stage due to the weak fault symptoms and strong noise. A strategy of combination of intrinsic characteristic-scale decomposition (ICD) and TQWT is proposed. ICD owns significant advantages on computation efficiency and alleviation of mode mixing. The effectiveness of the proposed strategy is tested with both simulated and experimental vibration signals. Meanwhile, comparisons are conducted between the proposed method and other methods like: envelope demodulation and EEMD-TQWT. Results show that the proposed method has superior performance in extracting fault features of defective bearings at an early stage.

1. Introduction A rolling bearing is one of the most important and commonly used components in rotating machines. Since rolling bearings often work under high loading and severe environment, local defect often develops gradually into failure on bearing. If there are no early warnings, the bearing failure could lead to unpredicted productivity loss for production and significant maintenance costs. Therefore, bearing fault detection and diagnosis have attracted substantial attention in recent years [1–3]. Vibration signal analysis is most widely used to diagnose bearing faults. When a bearing operates with local defect, it would excite the resonance frequency modulated by characteristic frequencies caused by bearing defects [4]. This would cause the measured vibration signal with the nonlinear and non-stationary characteristics [5]. However, at the early stage, the fault symptoms generated by a localized damage has low amplitude, which are often immersed in heavy background noises and rotation frequencies. Common used indexes, such as: RMS and kurtosis, are effective to describe the degradation process of a bearing. However, the bearing fault type cannot be identified easily using these indexes and the certainty of the fault occurrence requires to be double assured. As such, extraction of the weak fault characteristics of bearings at early stage is an important part of maintenance procedures [6–8]. ⁎

Corressponding author. E-mail address: [email protected] (M. Xu).

http://dx.doi.org/10.1016/j.ymssp.2016.10.013 Received 5 March 2016; Received in revised form 30 August 2016; Accepted 8 October 2016 0888-3270/ © 2016 Elsevier Ltd. All rights reserved.

Mechanical Systems and Signal Processing 86 (2017) 204–223

Y. Li et al.

Many vibration signal processing techniques have been developed to extract the fault features [6–12]. The turntable Q-factor wavelet transform (TQWT) was original proposed by Selesnick, and the most advantage of TQWT is that the Q-factor is easily and continuously adjustable [13,14]. Later, Luo [15] and He [16] provide a deep and detailed analysis of TQWT in terms of the filter bank, decomposition levels and Q-factors. From then on, TQWT theory was rapidly developed. When a fault occurs on bearings, the measured vibration signals contain both the fault transient impact component and sustained oscillatory component. Therefore, the TQWT has the high Q-factor and low Q-factor simultaneously, which is very suitable to separate the fault related components from the sustained oscillatory components. Increasing attentions have been aroused in applying TQWT for fault diagnosis of rolling bearings [15–18]. Wang et al. [7] recently proposed ensemble empirical mode decomposition (EEMD) and TQWT method to extract the weak fault features of rolling bearing in an accelerated bearing life test. In his method, EEMD is applied to decompose the vibration signals into a series of IMFs and then the TQWT is employed to separate the main IMF into high Q-factor component and low Q-factor component to diagnose the early fault. However, there are two main shortcomings in the Wang's method. First, although the TQWT is an interesting attempt to extract the fault characteristics, the selection of the proper Q-factors is still a problem [16]. As the filtered signal processed by Q-factors is far away from its theoretical central frequency and bandwidth is often heavily polluted, it is difficult to represent the factual features clearly. We believe that there are optimal Q-factors, which can reflect the characteristics of a fault better than the experiences. It is reasonable to find out these optimal Q-factors for fault feature extraction. Based on this consideration, we firstly propose a practical criterion for Q-factor determination, called characteristic frequency ratio (CFR), to select the optimal Q-factors that match the fault periodical impulses well. The CFR guided solution for Qfactor selection is called optimized TQWT (OTQWT). Second, EEMD is time consuming in decomposition. This paper will address this problem of EEMD. In this paper, a new weak fault feature extraction method is proposed based on intrinsic character-scale decomposition (ICD) and OTQWT. ICD method is proposed by Li [19] to decompose the vibration signals and the decomposed vibration signal can reduce the noises, which has the merits in suppressing mode mixing phenomenon and improving computation efficiency. To verify the effectiveness of the proposed method based on ICD and OTQWT in the early fault diagnosis, two accelerated life tests of rolling bearing are used to simulate the natural degradation process and pick out early fault stage. These accelerated run-to-failure experiments can not only avoid the inaccuracy of mature and or simulated bearing damages but also save the time in the normal operation. This paper combines the advantages of ICD and OTQWT to extract the fault feature of rolling bearing at an early stage of accelerated life tests. Results show that the proposed method has super performance in the weak feature extraction compared with envelope demodulation and EEMD-TQWT. The rest of this paper is organized as follows. Section 2 describes the ICD method. Section 3 introduces TQWT method, then, an optimization criterion is proposed to select the appropriate Q-factors. Section 4 provides the detailed procedures of the proposed early fault diagnosis approach and its effectiveness is validated via numerical simulations. Section 5 applies the proposed method to identify the fault types of the rolling bearings. Finally, conclusions are drawn in Section 6.

2. Intrinsic characteristic-scale decomposition method 2.1. A description of ICD method The ICD method is proposed by Li [19], which can decompose the vibration signal into a set of product components (PCs), each of which is the product of an amplitude envelope signal and a purely frequency modulated signal, which can be expressed as follows. n

x (t ) =

∑ PCi (t ) + un (t )

(1)

i =1

where u (t ) is the residual. A flowchart of ICD method is presented in Fig. 1. There are four main steps of ICD. First, identify the local maxima and minima points for a given signal. Second, calculate the reference points Ak to obtain the local mean points and local envelope points. Third, connect all the mean points and envelope points to construct the local mean functions and corresponding the envelope functions. Last, conduct the local mean decomposition using the obtained local mean functions and envelope functions to generate a series of PCs. Detailed explanations of this self-adaptive decomposition method are given in Ref. [19]. After the original signal x(t) is decomposed by ICD method, multiple PCs can be obtained, which are arranged from high frequency to low frequency. However, some of PCs contain mainly fault information, while others are related to the fault independent components and the noise. Since kurtosis index is sensitive to periodical impulse, and the PC with highest kurtosis value implies most impulsive characteristics [20]. Kurtosis index is utilized to select the principal PC component from the ICD decomposition results in this paper.

2.2. Simulated signal analysis To validate the effectiveness of ICD method, a comparison between the ICD method and EEMD method is conducted using a synthetic signal as follows: 205

Mechanical Systems and Signal Processing 86 (2017) 204–223

Y. Li et al.

Fig. 1. A flowchart of ICD method.

⎧ x = x1 (t ) + x2 (t ) + w (t ) ⎪ ⎨ x1 (t ) = (1 + 0.5 sin(8πt ))sin(120πt ) ⎪ x2 (t ) = sin(70πt ) e−t /2 ⎩

(2)

where t = 0: 1/1000: 2 , and the sampling frequency is 1000 Hz. w (t ) is additive white Gaussian noise (AWGN) with standard deviation σ = 0.2 . The simulation signal consists of two amplitude-modulated signals and a noise signal. The time domain waveforms of x (t ), x1 (t ) and x2 (t ) are shown in Fig. 2. The simulation signal x (t ) is decomposed using EEMD [7] and ICD [19] methods, respectively. The decomposition results are shown in Figs. 3 and 4, respectively. We can find that the PC2(t) and IMF2(t) correspond to x1 (t ), while PC3(t) and IMF3(t) correspond to x2 (t ). EEMD and ICD can both separate the noisy signal and provide a clear representation of the two components of the original signal. According to Ref. [7], the amplitude of the white noise added in the EEMD is 0.2σ and the number of the noise added in EEMD is the 1/10 of the number of the data points. To compare the decomposition performance of the two methods, root mean squared error (RMSE) and CPU running time consuming are considered. A laptop with 3.3 GHz i3-Core CPU, 4.0 GB RAM and MATLAB (R2010b) is used to conduct the simulation. The definition of RMSE is given as follows: 206

Mechanical Systems and Signal Processing 86 (2017) 204–223

Y. Li et al.

Fig. 2. The time domain of waveform of the signal x (t ),x1 (t ) and x2 (t ) in Eq. (2).

Fig. 3. EEMD decomposition results of the simulation signal x (t ) .

RMSE =

E (s (t ) − PC (t ))2

(3)

where s(t) and PC(t) are the original defined component and the corresponding decomposed component, respectively. The smaller RMSE value indicates the better decomposition performance. The comparison results using the two assessing indicators are shown in Table 1. It can be found from Table 1 that the capability of the two methods is almost the same even the RMSE values of EEMD are a bit smaller. However, the CPU running time of ICD is much less than EEMD. ICD took less than half second while the EEMD used more than 2 min. ICD improves the computation efficiency greatly, which is suitable in on-line monitoring and diagnosis of bearings. 3. Tunable Q-factor wavelet transform 3.1. Tunable Q-factor wavelet transform It is known that wavelet bases are characterized by their Q-factors, which is defined as follows:

Q=

f0 BW

(4)

where f0 is the center frequency and BW is the bandwidth of a pulse signal. Since the measured fault bearing vibration signal contains both the fault transient impact component (with low Q-factor) and the sustained oscillation component (with high Q-factor). The single Q-factor wavelet transform is not suitable for analyzing bearing vibration signals with both high and low Q-factor components. To overcome the above shortcoming, TQWT is proposed by Selesnick to separate the transient impact component related to fault information from the sustained oscillation components related to the rotation frequency and background noises [13,14]. As shown in Fig. 5, TQWT can be implemented using two-channel filter banks: the low-pass scaling (with the parameter α) and high-pass scaling (with the parameter β). To achieve a wavelet transform with the desired Q-factor and oversampling rate r, it is important to select the appropriate filter bank parameters α and β. The relationship between the two parameters (α, β, the Q-factor 207

Mechanical Systems and Signal Processing 86 (2017) 204–223

Y. Li et al.

Fig. 4. ICD decomposition results of the simulation signal x (t ) . Table 1 The comparison results between the EEMD and ICD methods. Methods

RMSE

EEMD ICD

CPU running time (s)

x1 (t)

x2 (t)

0.1021 0.1032

0.0932 0.0954

150.72 0.39

and the redundancy r) is expressed in Eq. (5).

α=1−

β 2 , β= γ Q+1

(5)

To make a perfect reconstruction, the low-pass filter H0(w) and the high-pass filter H1(w) are constructed as follows:

⎧ w ≤ (1 − β ) π 1, ⎪ w + (β − 1) π H0 (w ) = ⎨ θ ( α + β − 1 ), (1 − β ) π ≤ w < απ ⎪ απ ≤ w < π 0, ⎩

(6)

⎧ w ≤ (1 − β ) π 0, ⎪ απ − ω H1 (w ) = ⎨ θ ( α + β − 1 ), (1 − β ) π ≤ w < απ ⎪ απ ≤ w < π 1, ⎩

(7)

where θ (⋅) is a function, which can be expressed as θ (ν ) = 0.5(1 + cosν ) 2 − cos ν , υ ≤ π . It can be found from Eqs. (5)–(7) that TQWT is easily parameterized by its Q-factor and its oversampling rate (redundancy r). Different (Q, r) are associated with different wavelet basis and different frequency responses [7,14]. The number of the decomposition levels (or stages) is denoted by L. When implement a L-stage TQWT decomposition, L+1 sub-bands are generated: the high-pass filter output signal of each filter bank {Vhj (n ), j = 1, ⋯l}, and the low-pass filter output signal of the final filter bank Vll (l ). A illustration of the L-stage TQWT decomposition is shown in Fig. 6. Meanwhile, the maximum number of allowable decompositions Lmax can be calculated using the (Q, r) [13].

Fig. 5. Block diagram of analysis and synthesis filter banks for the application of TQWT.

208

Mechanical Systems and Signal Processing 86 (2017) 204–223

Y. Li et al.

Fig. 6. Application the TQWT with L-stage decomposition of TQWT.

⎢ ⎥ log(N /4(Q + 1)) L max = ⎢ ⎥ ⎣ log(Q + 1/(Q + 1 − 2)/ r ) ⎦

(8)

where N represents the data length and ⌊⋅⌋ represents round-down integer value. Because the observed signal x is the combination of an oscillatory signal x1 and a non-stationary signal x2, which can be expressed as x= x1+ x2. Considering that x1 consists mostly of sustained oscillations and x2 consists mostly of non-oscillatory transients, we aim to decompose the signal x into two main components using high Q-factor and low Q-factor. In this case, these decompositions cannot be accomplished by applying frequency-based filtering [21]. With the help of morphological component analysis (MCA) [22], two components with different oscillatory behaviors can be separated successfully. Denote TQWT1 and TQWT2 as the TQWT with two different Q-factors (high Q-factor and low Q-factor). Then, the sought decomposition can be achieved by solving the constrained optimization problem [22].

⎧ arg min λ1 w1 + λ2 w2 ⎪ w1 w 2 ⎨ ⎪ x = TQWT−1(w ) + TQWT−1(w ) ⎩ 1 2 1 2

(9)

For greater flexibility, sub-band-dependent regularization is given as follows: j2 +1 ⎧ arg min ∑ j1 +1 λ w 1j 1 + ∑ j =1 λ2j w2j ⎪ j =1 1j ⎨ w1 w2 ⎪ x = TQWT1−1(w1) + TQWT−1 ⎩ 2 (w 2 )

1

(10)

where Wij in Eq. (10) denotes sub-band j of TQWTi (i=1,2). After w1 and w2 are obtained, we set

x1 = TQWT1−1(w1),

x2 = TQWT−1 2 (w 2 )

(11)

In the practical situation, the background noise is inevitably occurred with the x1 (low Q-factor transient impact component) and x2 (the high Q-factor sustained oscillation component). Therefore, the practically fault signal can be expressed as: (12)

y = x1 + x2 + Noise Then the problem can be converted into solving the Eq. (13): j1 +1

arg min y − φ1 w1 − φ2 w2 + w1, w 2

∑ λ1j j =1

j2 +1

w1j

1

+

∑

λ2j w2j

1

j =1

(13)

where φ1 and φ2 correspond to high Q-factor and low Q-factor, respectively. The regularization parameters λ1 and λ2 are selected according to the power of the noise. It should be noted that the above Eqs. (10)–(13) can be referred to [14,23]. 3.2. The parameters of TQWT and existing problems In TQWT method, three parameters determine its decomposition performance: the Q-factors, the redundancy r, and the number of decomposition levels L. Q-factor reflects the oscillations of the wavelet basis, which plays a significant role in its decomposition performance. The value of the Q-factors can be chosen in a wide range and the high Q factor should be greater than the low Q factor [7]. According to the Ref. [14], the redundancy r must be strictly greater than 1. However, for a given Q-factor, redundancy r can affect the overlap between adjacent frequency responses. For a larger r, a larger L value is required to cover the same frequency range due to the increased overlap [17]. In practical application, it is suggested to set r ≥ 3. Considering the computation efficiency, we set r = 3 in this paper. The number of decomposition levels L determines the number of filter banks. In this paper, we calculate the decomposition levels L according to Eq. (8). Above all, the accuracy of TQWT decomposition results heavily depend on the selection of Q-factors [7,18]. The main obstacle lies in the selection of proper Q-factors. However, the reported Q-factor selection methods are based on the experience of the operators and prior knowledge of the data. The application of CFR in Q-factor selection is presented in the following section. 3.3. The characteristic frequency ratio (CFR) based solution for Q-factors selection We aim to design a process to quantify the accuracy of the TQWT decomposition results and find the optimal one. Therefore, a selection criterion called CFR is established in this paper. The CFR is defined based on the Hilbert demodulation analysis with a 209

Mechanical Systems and Signal Processing 86 (2017) 204–223

Y. Li et al.

definite physical meaning. The Hilbert transform (HT) technique has been proven to be an effective tool for the detection of localized defects in bearings [7,23]. It extracts bearing characteristic defect frequencies along with the modulation. Meanwhile, the localized defects emerged on rolling bearings can cause the transient impacts. In TQWT method, the low Q-factor component has fewer sign changes and consists of fewer oscillatory cycles, which is more suitable to extract the fault-induced impulses. As such, the HT method is used to demodulate the low Q-factor component of TQWT to calculate the instantaneous amplitude (IA) signal of x(t) through the following way:

y (t ) =

1 π

∞

∫−∞ tx−(τ τ) dτ

(14)

Then, perform the inverse transform to obtain the x(t) as

x (t ) =

1 π

∞

∫−∞ ty−(τ )τ dτ

(15)

Coupling the x(t) and y(t), we can have an analytic signal z(t), as

z (t ) = x (t ) + iy (t ) = a (t ) eiφ (t )

(16)

Calculate the IA a(t) of x(t) as follows:

a (t ) =

x 2 (t ) + y 2 (t ) ,

⎛ y (t ) ⎞ φ (t ) = arctan ⎜ ⎟ ⎝ x (t ) ⎠

(17)

where φ (t ) is the instantaneous phase of x (t ). Based on HT, we operate the FFT on the IA to obtain the corresponding envelope spectrum. The defect frequencies of bearings with different fault types can be calculated using the calculation equations listed in Table 2. Since the fault type is unknown to us when an incipient fault occurred on the bearings, the proposed criterion should involve all the fault frequencies mentioned above. Based on the Hilbert envelope spectrum analysis, the CFR can be defined as follows: M

CFR =

∑ (Akfo

N

+ Akfi + Akfe + Akfc )/ ∑ A f j

k =1

(18)

j =1

where Akfo , Akfi , Akfe and Akfc represent the amplitude of the kth harmonic of the outer race fault frequency, inner race fault frequency, element fault frequency and cage fault frequency, respectively. M is the number of harmonics of the fault frequency and N is the number of frequency components in the frequency spectrum. CFR index is a novel measure of the periodic impulses, which is sensitive to the fault features excited by a localized defect on bearings. Therefore, CFR can be taken as a quantitative index to estimate the decomposition performance and select the optimal Q-factors. A larger CFR value implies a higher capability of fault feature extraction. In this work, the M is set to be 3 and N is set to be 500. Combined with CFR index, the optimized TQWT (OTQWT) is proposed to select most suitable wavelet basis that matches the oscillatory behavior of fault signatures well, and the optimization process is implemented to maximize the CFR value. The calculation steps of CFR for a given signal are as follows: (1) (2) (3) (4)

Decompose the original vibration signal into high Q-factor component and low Q-factor component; Perform HT method to the low Q-factor component to obtain the IA; Conduct FFT in Step (2) and obtain the amplitudes of the fault characteristic frequencies in Table 2; Compute the CFR value based on the obtained fault characteristic frequencies.

Table 2 Calculation equations of bearing characteristic frequencies. Bearing characteristic frequency

Equations

Ball pass frequency of the outer race (BPFO)

1 ⎛ fo = Z ⎜1 − 2 ⎝

⎞ d cosα⎟ fs D ⎠

Ball pass frequency of the inner race (BPFI)

1 ⎛ fi = Z ⎜1 + 2 ⎝

⎞ d cosα⎟ fs D ⎠

⎡

Ball pass frequency of the rolling element (BPFR)

Cage defect frequency (CDP)

⎤ ⎛ d ⎞2 − ⎜ ⎟ cos 2α ⎥ fs D⎠ ⎝ ⎥⎦ ⎣

fe =

D⎢ 1 2d ⎢

fc =

1⎡ ⎢1 2⎣

Note that fs is the rotating frequency, D is the pitch diameter, d is the ball diameter, Z is the number of the balls, and α is the contact angle.

210

−

⎤ d cosα ⎥ fs D ⎦

Mechanical Systems and Signal Processing 86 (2017) 204–223

Y. Li et al.

To find out the optimal high and low Q-factors in TQWT simultaneously, a discrete 2D exhaustive searching method is utilized. According to the studies in [13,18], the parameter searching ranges are set as follows: high Q-factor (4,12) and low Q-factor (1,3). The optimization step size is selected to be 0.5 according to Ref. [11]. Additionally, the decomposition level L is determined by the Eq. (8) and parameter r is set to 3. 4. The proposed fault diagnosis method 4.1. The fault feature extraction based on ICD and OTQWT Motivated by the advantages of ICD and OTQWT, the ICD-OTQWT technology for early fault diagnosis of rolling bearings is proposed. The work steps are summarized as follows: (1) Apply ICD to decompose the measured vibration signal into a set of PC components; (2) Compute the kurtosis value of each PC obtained in step (1) and select the main component according to the maximum kurtosis value; (3) Determine the range of high and low Q-factors of TQWT: the range for the low Q-factor is (1,3) and that for the high Q-factor is (4,12), with the increase of 0.5; (4) Conduct the TQWT operation; (5) Calculate the CFR of each TQWT at each pair of Q-factors in step (4). The larger value of CFR implies a better extraction performance; (6) Pick out the optimal Q-factors according to the largest CFR value and only this pair of Q-factors will be used to generate the high and low Q-factor components; (7) Implement the envelope analysis to the low Q-factor component for the early fault diagnosis. It is noteworthy that no parameter needs to be defined beforehand for the proposed method. The proposed approach is easy for application as the optimal Q-factors are easy to get. In addition, one advantage of this method is that it is unnecessary to know the fault types in advance to select the optimal Q-factors. This advantage can make the method be effective in the industrial applications. A flow chart of the proposed algorithm is presented in Fig. 7. 4.2. Algorithm validation with synthetic signal To validate the effectiveness of the proposed algorithm, a bearing fault model in Ref. [6] is employed to simulate the rolling bearing with localized defect. The bearing model consists of the effects of bearing geometry, shaft speed, the load, decaying exponential and so on [24,25]. The mathematical model of a rolling bearing can be used to generate the vibration signal with local defect as follows:

⎧ x (t ) = ∑M A s (t − iT − τ ) + w (t ) i i =1 i ⎪ ⎨ Ai = A0 cos(2πQt + ϕA) + CA ⎪ s (t ) = e−Bt sin(2πfn t + ϕw ) ⎩

(19)

where Ai is the amplitude modulation signal with a period time 1/ Q , A0 is the amplitude of the signal, CA is a constant with restriction

Fig. 7. Framework of the proposed algorithm.

211

Mechanical Systems and Signal Processing 86 (2017) 204–223

Y. Li et al.

Fig. 8. The waveforms of original vibration signal in Eq. (19).

Fig. 9. The obtained results of the simulated fault signal using the proposed method: (a) high Q-factor component; (b) low Q-factor component; (c) residual component.

CA > A0 , s (t ) is the discrete oscillating impulse signal with an interval time T between two adjacent impacts, B is the damping coefficient, fn is the natural frequency of the system, τi is the time lag derived from the random slip of rolling elements and w (t ) is the white Gaussian noise. We use the mathematical model of a rolling bearing as given in Eq. (19) to simulate the signals of a rolling bearing with weak fault. The parameters are given as follows: the natural frequency of the system fn is 3000 Hz, the rotating frequency fr is 33 Hz, the fault frequency of outer race fo is 124 Hz, and the data length L is 4096 with the sampling frequency fs =12 kHz. The maximal slip ratio of its period T=0.01. Meanwhile,w (t ) is the WGN with SNR=−12 dB. The above setting is used to simulate the early bearing fault under strong noisy environment. The time domain waveforms of the simulated signal and the decomposition results of the proposed method are displayed in Fig. 8 and Fig. 9, respectively. Then, the envelope spectrum of the low Q-factor component is shown in Fig. 10. It can be clearly seen from Fig. 10 that the rotation frequency fr and the fault characteristic frequency of outer race fo and its harmonics can be extracted. This validates that the proposed method can extract the early fault characteristics of outer race fault under strong background noise effectively. 5. Early fault detection of rolling bearings At present, most researches on bearing fault diagnosis focus on studying bearing with mature faults or simulated damages [8– 12]. However, these test bearings could not exhibit the degradation propagation at early stage of bearing defect. To assess the diagnostic ability of the proposed method, we select two accelerated experiments (Experiment 1 and Experiment 2) under constant speed and load conditions. Experiment 1 and Experiment 2 are both accelerated life tests, which run from their normal conditions to the fault initiation and then to the final failure. The two experiments have different working conditions (i.e. rotation speed and radial load) and different types of rolling bearings. The detailed introduction of the two experiments is given later in Section 5.1 and

Fig. 10. The envelope spectrum of the low Q-factor component derived from the proposed method.

212

Mechanical Systems and Signal Processing 86 (2017) 204–223

Y. Li et al.

Fig. 11. (a) Bearing test experimental system; (b) the diagram of the system.

Section 5.2, respectively. For comparison, the EEMD and TQWT in Ref. [7] (EEMD-TQWT), the ICD and TQWT (ICD-TQWT) and the ICD and OTQWT (ICD-OTQWT) are all applied to analyze the same vibration signal and extract the early fault features. We use characteristic frequency intensity coefficient (CFIC) [26] and CPU running time to evaluate the performance of the three methods. 5.1. Experiment 1 The first experiment data was contributed by Intelligent Maintenance System, University of Cincinnati [27], and the bearing test experimental system is shown in Fig. 11. As shown in Fig. 11(a), the experiment performed bearing accelerated life tests under constant load conditions on a specially designed test rig. The bearings have 16 rollers in each row, a pitch diameter of 2.815 in., roller diameter of 0.331 in., and a tapered contact angle of 15.17◦. On each bearing, two PCB 353B33 High Sensitivity Quartz ICP accelerometers were installed for data acquisition (one vertical and one horizontal) as shown in Fig. 11(b). Meanwhile, the vibration signals were measured at an interval of every 10 min. Besides, the sampling frequency is 20000 Hz and the data length is 20480 points. Also the shaft rotating speed of the motor was 2000 rpm. A radial load of 2721.6 kg is applied on bearings. Four bearings named Bearing 1–4 were installed on one shaft as shown in Fig. 11. At the end of the accelerated experiment, an outer race fault was discovered in test bearing 1. More detailed information about this experiment is shown in literature [27]. Since kurtosis index is more sensitive than RMS to the impulses when a localized fault occurs on rolling bearings [7,28]. The kurtosis index is employed to reflect the degradation processes over its whole life. The expression of kurtosis is given as follows: N

kurtosis =

(1/ N ) ∑n =1 (x (n ) − x )4 N

[(1/ N ) ∑n =1 (x (n ) − x )2 ]2

(20)

where x(n) is a raw experimental vibration signal and x is the mean value of all amplitudes. The obtained result using kurtosis index is shown in Fig. 12. Here, we use group as the x-axis unit as did in [6]. One group corresponds to 10 min in Fig. 12. According to the calculation equations in Table 2, the fault characteristic frequencies of the tested rolling bearings can be computed and listed in Table 3. Seen from Fig. 12, it is clearly observed that the kurtosis index increases abruptly at 643rd group. According to [6,7], the 643rd group can be seen as the severe fault occurring time. To achieve the requirement of condition monitoring and device maintenance, the weak fault features of rolling bearing should be extracted [6]. Before we perform the proposed method, we should select the vibration signal prior to its completed failure as the early fault initiation of rolling bearing. Since the kurtosis value at 534th group is almost the same as the kurtosis values of its normal period 1533 groups), the 534th group can be considered as the early fault stage of rolling bearings. The waveform of the vibration signal and its corresponding envelope demodulation spectrum are shown in Fig. 13(a) and (b), respectively.

Fig. 12. The kurtosis values of the test bearing in Experiment 1 over its whole life.

213

Mechanical Systems and Signal Processing 86 (2017) 204–223

Y. Li et al.

Table 3 The theoretical characteristic frequencies values in two experiments. Experiments

Rotation frequency fr (Hz)

BPFO fo (Hz)

BPFI fi (Hz)

BPFR fe (Hz)

BPFC fc (Hz)

1 2

33.3 50

236.4 203.9

296.9 296.1

139.9 262.0

14.8 20.4

Fig. 13. (a) the waveforms of vibration signal of the bearing with early fault (b) the corresponding envelope spectrum of the original signal.

We cannot find the outer race fault frequency fo = 236.4Hz in Fig. 13(b), which implies the ineffectiveness of application envelop demodulation on the original early weak fault vibration signal directly. This can be attributed to that the weak period impulse excited by the outer race fault is totally merged in the background noise. Then the proposed method is applied to extract the fault features. First, apply ICD to decompose the original vibration signal. The obtained decomposition results are displayed in Fig. 14. Second, calculate the kurtosis values of PCs and select the PC2(t) that contains the most fault information according to its maximum kurtosis value. Third, conduct the OTQWT (the optimal high and low Q-factors are selected as 1 and 6.5, respectively) on PC2(t) and the obtained results are shown in Fig. 15. Last, apply the envelope demodulation on the low Q-factor component, and the obtained Hilbert spectrum of the low Q-factor component is illustrated in

Fig. 14. ICD decomposition results of the vibration signal with early fault.

214

Mechanical Systems and Signal Processing 86 (2017) 204–223

Y. Li et al.

Fig. 15. The decomposition results of PC2 using TQWT (a) high Q-factor component; (b) low Q-factor component; (c) residual component.

Fig. 16. The envelope spectrum of the low Q-factor component in Fig. 15.

Fig. 16. The outer race fault frequency fo = 234Hz and its side frequencies fo + fr , fo −fr , fo + 2fr , fo −2fr can all be extracted perfectly. These frequency components clearly indicate outer inner race fault on the bearings. To demonstrate the advantages of the proposed method, the EEMD-TQWT [7] and ICD-TQWT are also employed to analyze the original weak fault being signal. The decomposition results of EEMD are shown in Fig. 17 (same as did in [7], the amplitude and number of the white noise added in EEMD are 0.2σ and 1/10 of the number of the data points, respectively), and then TQWT is employed to decompose IMF3 as shown in Fig. 18. The envelope spectrum of the low Q-factor component in Fig. 18 is presented in Fig. 19. Meanwhile, the decomposition result of ICD-TQWT and its corresponding envelope spectrum of its low Q-factor component are shown in Fig. 20 and Fig. 21, respectively. We can find that the outer race fault fo = 234Hz can both be detected in Fig. 19 and Fig. 21. To evaluate the performance of the three methods, a characteristic frequency intensity coefficient (CFIC) [26] is employed in this paper. A larger value of CFIC indicates a better fault feature extraction performance. The definition of CFIC is expressed as follows: N

CFIC =

1 ∑i =1 Aifc

N

2 ∑ j =1 A fj

(21)

where Aifc and A f j are the amplitude of the ith harmonic of the fault characteristic frequency and the frequencies analyzed, respectively. N1 and N2 are the number of harmonics of the fault characteristic frequency and all frequency components, respectively. In this paper, we set N1=1 and N2=350. Here, we take five characteristic frequencies to calculate CFIC as: fo , fo + fr , fo − fr , fo + 2fr , fo − 2fr . The CFIC values and the CPU running time of the three methods are shown in Table 4. From Table 4, we can find that the proposed method has the biggest CFIC value, which implies the best performance in the early fault extraction among the three methods. Compared with ICD-TQWT, higher CFIC value validates the advantage of the optimization process in the proposed OTQWT method. From Table 4, it is obvious that ICD-OTQWT takes much less time than EEMD-TQWT method. ICD-TQWT method has the highest computation efficiency. Therefore, the proposed method can not only improve the ability of the weak feature extraction but also takes much less time in comparison with EEMD-TQWT method.

215

Mechanical Systems and Signal Processing 86 (2017) 204–223

Y. Li et al.

Fig. 17. EEMD decomposition results of the vibration signal with early fault.

Fig. 18. The decomposition results of IMF3 using TQWT (a) high Q-factor component; (b) low Q-factor component; (c) residual component.

Fig. 19. The envelope spectrum of the low Q-factor component in Fig. 18.

5.2. Experiment 2 The second accelerated bearing life tester (ABLT-1) is provided by the Hangzhou Bearing Test & Research Center (HBRC) in China. The designed tester has four bearings on one shaft driven by an AC motor and coupled by rub belt, the sketch of the bearing 216

Mechanical Systems and Signal Processing 86 (2017) 204–223

Y. Li et al.

Fig. 20. The decomposition results of PC2 using TQWT (a) high Q-factor component; (b) low Q-factor component; (c) residual component.

Fig. 21. The envelope spectrum of the low Q-factor component in Fig. 20.

accelerated life test system is shown in Fig. 22(a). In this experiment, 6211 single row deep groove ball bearings are used in the full life test. One group of vibration data with 20480 points are collected every 5 min by data acquisition card. The sampling frequency is 12000 Hz and the load used on the bearings is 1445.2 kg. The analyzed bearing is confirmed with inner race fault after the accelerated life test. Fig. 22(b) shows the rolling bearing with early inner race fault. Similar to Experiment 1, the kurtosis index of the test bearing 2 over its whole life is calculated as shown in Fig. 23. As seen in Fig. 23, it is clearly observed that the kurtosis index increases abruptly at 2438th group, which implies the complete fault occurs. Then we select the 2420th group (earlier than the 2438th group) to test the methods as the kurtosis value of the 2420th group is almost the same as that of normal condition of rolling bearing. The original time-waveform of the analyzed signal and its corresponding envelope spectrum are shown in Fig. 24(a) and Fig. 24(b), respectively. From Fig. 24(b), we cannot find the fault frequency of the inner race fault fi = 296.1Hz , which demonstrates that the envelope demodulation cannot extract the early fault from the original signal. Subsequently, we apply the proposed approach to extract the weak feature. First, apply ICD to decompose the vibration signal, and the obtained results are shown in Fig. 25. Second, calculate the kurtosis values of PCs and select the PC2(t) with the highest kurtosis for further analysis. Third, conduct the OTQWT (the optimal high and low Q-factors are 1 and 4.5, respectively) on the PC2(t) and the obtained results are shown in Fig. 26. Lastly, envelope demodulation technique is employed to analyze the low Q-factor component as shown in Fig. 27. We can see from the Fig. 27 that the inner race fault frequency fi = 296.1Hz together with fi + fr , fi − fr and the first and second harmonic of modulation frequency fr and 2fr are all observed. For comparison purpose, the EEMD-TQWT and ICD-TQWT are also employed to analyze the original weak fault being signal.

Table 4 The CFIC and computation time of the three methods in Experiment 1. Methods

CFIC values

CPU running time

EEMD-TQWT ICD-TQWT ICD-OTQWT

0.0205 0.0212 0.0431

971.57 4.564 244.04

217

Mechanical Systems and Signal Processing 86 (2017) 204–223

Y. Li et al.

Fig. 22. (a) A sketch of the bearing accelerated life test system (b) a rolling bearing with inner race fault.

Fig. 23. The kurtosis values of the test bearing in Experiment 2 over its whole life.

Fig. 24. (a) the waveforms of vibration signal of the bearing with early fault (b) the envelope spectrum of the original signal.

The decomposition results of EEMD are shown in Fig. 28, and then TQWT is employed to decompose IMF2. The results are shown in Fig. 29. The envelope spectrum of the low Q-factor component in Fig. 29 is presented in Fig. 30. Meanwhile, the decomposition result of ICD-TQWT and its corresponding envelope spectrum of its low Q-factor component are shown in Fig. 31 and Fig. 32, respectively. The inner race fault fi = 296.1Hz can be identified both in Fig. 30 (EEMD-TQWT method) and Fig. 32 (ICD-TQWT method). However, the fi is not particularly evident. Therefore, the accuracy of using EEMD-TQWT and ICD-TQWT to extract the 218

Mechanical Systems and Signal Processing 86 (2017) 204–223

Y. Li et al.

Fig. 25. ICD decomposition results of the vibration signal with early fault.

Fig. 26. The decomposition results of PC2 using TQWT (a) high Q-factor component; (b) low Q-factor component; (c) residual component.

Fig. 27. The envelope spectrum of the low Q-factor component in Fig. 26.

early fault feature is actually modest. The CFIC and the CPU running time of the three methods are computed and the comparison results are presented in Table 5. The CFIC in Table 5 also expresses a best performance of the proposed method. Compared with EEMD-TQWT, the proposed method saves much computation time. 219

Mechanical Systems and Signal Processing 86 (2017) 204–223

Y. Li et al.

Fig. 28. EEMD decomposition results of the vibration signal with early fault.

Fig. 29. The decomposition results of IMF2 using TQWT (a) high Q-factor component; (b) low Q-factor component; (c) residual component.

220

Mechanical Systems and Signal Processing 86 (2017) 204–223

Y. Li et al.

Fig. 30. The envelope spectrum of the low Q-factor component in Fig. 29.

Fig. 31. The decomposition results of PC2 using TQWT (a) high Q-factor component; (b) low Q-factor component; (c) residual component.

Fig. 32. The envelope spectrum of the low Q-factor component in Fig. 31.

Table 5 The CFIC and computation time of the three methods in Experiment 2. Methods

CFIC values

CPU running time

EEMD-TQWT ICD-TQWT ICD-OTQWT

0.0113 0.0120 0.0212

1293.80 5.58 253.68

221

Mechanical Systems and Signal Processing 86 (2017) 204–223

Y. Li et al.

6. Conclusions Early fault detection of bearings is essential to reduce failure related costs by avoiding unexpected failures. This paper presents a novel fault diagnosis approach based on ICD and OTQWT. The experimental vibration signals collected from two bearing accelerated life test rigs are employed to evaluate the effectiveness of the proposed method. Test results demonstrate that the proposed method can achieve good performance in the extraction of weak fault characteristics. Another big advantage of the proposed method is increasing the computation efficiency greatly over the EEMD-TQWT method. The main contributions of this paper are summarized as follows: (1) The combination of ICD and TQWT saves computation cost. (2) A criterion called CFR is used to optimize Q-factors of TQWT, which improves the fault detection ability of ICD-TQWT. (3) Simulated and experimental signals demonstrate that the proposed method has a better fault detection ability and also the advantage of saving computation cost comparing with the existing method EEMD-TQWT. Generally, more validations are required for a method to be used in real application. In this preliminary study, the proposed method was tested and demonstrated to be effective in bearing fault detection using two sets of data collected from two accelerated life bearing test platforms under constant speed and load conditions. Further test under different speed and load combination conditions will be considered in our future work. Acknowledgements The research is supported by National Natural Science Foundation of China (NO. 11172078) and Important National Basic Research Program of China (973 Program-2012CB720003), and the authors are grateful to all the reviewers and the editor for their valuable comments. References [1] S. Janjarasjitt, H. Ocak, K.A. Loparo, Bearing condition diagnosis and prognosis using applied nonlinear dynamical analysis of machine vibration signal, J. Sound Vib. 317 (1–2) (2008) 112–126. [2] C.-S. Park, Y.-C. Choi, Y.-H. Kim, Early fault detection in automotive ball bearings using the minimum variance cepstrum, Mech. Syst. Signal Process. 38 (2) (2013) 534–548. [3] J. Yu, Local and nonlocal preserving projection for bearing defect classification and performance assessment, IEEE Trans. Ind. Electron. 59 (5) (2012) 2363–2376. [4] Y. Lei, J. Lin, Z. He, Y. Zi, Application of an improved kurtogram method for fault diagnosis of rolling element bearings, Mech. Syst. Signal Process. 25 (5) (2011) 1738–1749. [5] 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. [6] Y. Ming, J. Chen, G. Dong, Weak fault feature extraction of rolling bearing based on cyclic Wiener filter and envelope spectrum, Mech. Syst. Signal Process. 25 (5) (2011) 1773–1785. [7] H. Wang, J. Chen, G. Dong, Feature extraction of rolling bearing's early weak fault based on EEMD and tunable Q-factor wavelet transform, Mech. Syst. Signal Process. 48 (1–2) (2014) 103–119. [8] A. Soleimani, S.E. Khadem, Early fault detection of rotating machinery through chaotic vibration feature extraction of experimental data sets, Chaos Solitons Fractals 78 (2015) 61–75. [9] G.F. Bin, J.J. Gao, X.J. Li, B.S. Dhillon, Early fault diagnosis of rotating machinery based on wavelet packets—Empirical mode decomposition feature extraction and neural network, Mech. Syst. Signal Process. 27 (2012) 696–711. [10] S. Lu, Q. He, F. Kong, Stochastic resonance with Woods–Saxon potential for rolling element bearing fault diagnosis, Mech. Syst. Signal Process. 45 (2) (2014) 488–503. [11] W. He, Y. Zi, B. Chen, F. Wu, Z. He, Automatic fault feature extraction of mechanical anomaly on induction motor bearing using ensemble super-wavelet transform, Mech. Syst. Signal Process. 54–55 (2015) 457–480. [12] A. Tabrizi, L. Garibaldi, A. Fasana, S. Marchesiello, Early damage detection of roller bearings using wavelet packet decomposition, ensemble empirical mode decomposition and support vector machine, Meccanica 50 (3) (2014) 865–874. [13] I.W. Selesnick, Wavelet transform with tunable Q-factor, IEEE Trans. Signal Process. 59 (8) (2011) 3560–3575. [14] I.W. Selesnick, Sparse signal representations using the tunable Q-factor wavelet transform, 8138, 2011, p. 81381U–81381U–15 [15] J. Luo, D. Yu, M. Liang, A kurtosis-guided adaptive demodulation technique for bearing fault detection based on tunable-Q wavelet transform, Meas. Sci. Technol. 24 (5) (2013) 55009. [16] W. He, Y. Zi, B. Chen, S. Wang, Z. He, Tunable Q-factor wavelet transform denoising with neighboring coefficients and its application to rotating machinery fault diagnosis, Sci. China Technol. Sci. 56 (8) (2013) 1956–1965. [17] G. Cai, X. Chen, Z. He, Sparsity-enabled signal decomposition using tunable Q-factor wavelet transform for fault feature extraction of gearbox, Mech. Syst. Signal Process. 41 (1–2) (2013) 34–53. [18] X. Chen, G. Cai, H. Cao, W. Xin, Condition assessment for automatic tool changer based on sparsity-enabled signal decomposition method, Mechatronics 31 (2015) 50–59. [19] Y. Li, M. Xu, Y. Wei, W. Huang, Rotating machine fault diagnosis based on intrinsic characteristic-scale decomposition, Mech. Mach. Theory 94 (2015) 9–27. [20] Y. Li, M. Xu, R. Wang, W. Huang, A fault diagnosis scheme for rolling bearing based on local mean decomposition and improved multiscale fuzzy entropy, J. Sound Vib. 360 (2016) 277–299. [21] I.W. Selesnick, Resonance-based signal decomposition: a new sparsity-enabled signal analysis method, Signal Process. 91 (12) (2011) 2793–2809. [22] J.L. Starck, M. Elad, D.L. Donoho, Image decomposition via the combination of sparse representations and a variational approach, IEEE Trans. Image Process 14 (10) (2005) 1570–1582. [23] N.E. Huang, Z. Shen, S.R. Long, M.C. Wu, H.H. Shih, Q. Zheng, N.-C. Yen, C.C. Tung, H.H. Liu, The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis, Proc. R. Soc. Lond. A: Math., Phys. Eng. Sci. 454 (1998) 903–995. [24] R.B. Randall, J. Antoni, S. Chobsaard, The relationship between spectral correlation and envelope analysis in the diagnostics of bearing faults and other cyclostationary machine signals, Mech. Syst. Signal Process. 15 (5) (2001) 945–962.

222

Mechanical Systems and Signal Processing 86 (2017) 204–223

Y. Li et al.

[25] J. Antoni, F. Bonnardot, A. Raad, M. El Badaoui, Cyclostationary modelling of rotating machine vibration signal, Mech. Syst. Signal Process. 18 (6) (2004) 1285–1314. [26] B. Li, P. Zhang, Z. Wang, S. Mi, Y. Zhang, Gear fault detection using multi-scale morphological filters, Measurement 44 (10) (2011) 2078–2089. [27] J. Lee, H. Qiu, G. Yu, C. Lin, Bearing data set, IMS Univ. Cincinnati NASA Ames Progn. Data Repos. Rexnord Tech. Serv., 2007. [28] Z. Shen, Z. He, X. Chen, C. Sun, Z. Liu, A monotonic degradation assessment index of rolling bearings using fuzzy support vector data description and running time, Sensors 12 (8) (2012) 10109–10135.

223