Bearing incipient fault diagnosis based upon maximal spectral kurtosis TQWT and group sparsity total variation denoising approach

2.1. Maximal spectral kurtosis tunable $\mathit{Q}$-factor wavelet transform

The resonance characteristics of oscillatory signal can be depicted by quality factor $Q$, $Q=f_c/B_w$, where $f_c$ represents center frequency of signal and $B_w$ is bandwidth. The tunable $Q$-factor wavelet transform (TQWT) was proposed by Selesnick as a flexible discrete wavelet transform for oscillatory signal processing [16]. The TQWT method includes three easily changeable parameters, i.e. quality factor $Q$, total over-sampling rate $r$ and number of levels of decomposition $j$.

Fig. 1 illustrates the wavelet time domain waveform and frequency response curves with different $Q$-factors (e.g., $j=$ 2, $Q=$ 1, 2, 3, 4, 5, 6). Fig. 2 illustrates the wavelet waveform and frequency response curves with fixed $Q$-factors and different $j$ scales (e.g., $Q=$ 2.5, $j=$ 1, 2, 3, 4, 5, 6). From Fig. 1 and Fig. 2, it can be concluded that the quality factor $Q$ controls the oscillatory behavior and waveform shape of wavelet waveform, and the decomposition level $j$ controls the expansion extent and bandpass location of wavelet waveform.

Fig. 1Wavelet waveform and frequency response curves with different Q-factors (e.g., j= 2, Q= 1, 2, 3, 4, 5, 6)

a) Wavelet waveform

b) Frequency response curves

Fig. 2Wavelet waveform and frequency response curves with same Q-factors and different j scales (e.g., Q= 2.5, j= 1, 2, 3, 4, 5, 6)

a) Wavelet waveform

b) Frequency response curves

For each level of TQWT decomposition, the signal $s\left(n\right)$ with sampling frequency $f_s$ can be decomposed into low-pass and high-pass sub-band signals with sampling frequencies $\alpha fs$ and $\beta fs$, respectively, the parameters $\alpha$ and $\beta$ are scaling values. Moreover, the low-pass filter $F_0\left(\omega \right)$ and low-pass scaling $LPS\alpha$ are used to obtain low-pass sub-band. Similarly, the high-pass sub-band can be generated by high-pass filter $F_1\left(\omega \right)$ and high-pass scaling $HPS\alpha$ accordingly. The TQWT approach uses following given low-pass and high pass filters:

where $\theta \left(\omega \right)$ is the frequency response of the Daubechies filter which has two vanishing moments. The $\theta \left(\omega \right)$ can be defined as follows:

The $Q$-factor $Q$ and parameters $r$ can be formulated via filter bank parameters $\alpha$ and $\beta$ as follows [16]:

where $f_c$ represents center frequency of signal and $B_w$ is bandwidth of sub-band $j$. For more details about decomposition and reconstructed operations can be found from Selesnick [16].

For maximal spectral kurtosis TQWT method, it applies the principle of maximum kurtosis, the kurtosis of each scale wavelet coefficient can be calculated once the original vibration signals are decomposed by TWQT method, then the calculated Kurtosis value can be used for selecting the optimal frequency band, finally, the corresponding reconstructed signal is obtained by inverse TQWT method.

2.2. Group sparsity total variation denoising

Total variation denoising (TVD) assumes that the noisy signal $y$ is of the form [20-22]:

where $x$ is unknown signal and $w$ white Gaussian noise. The unknown signal $x$ can be obtained by solving the following optimization problem:

The regularization parameter $\lambda >$ 0 controls how much smoothing is performed, and matrix $D$ is first-order difference matrix, i.e.:

Note that $v=Dx$ and penalty function $\varphi \left(v\right)$ can be defined by:

Here, the parameter $k$ represents the group size.

If $k=$ 1, the penalty function $\varphi \left(v\right)={‖v‖}_1$ and Eq. (6) is the standard 1-D TVD problem. If $k>$ 1, the penalty function $\varphi \left(v\right)$ is a convex measure of group sparsity. In this paper, we discuss group-sparse total variation denoising (GS-TVD) method.

Definition: A $k$-point group can be denoted by the vector $v$:

Thus, the penalty function $\varphi \left(v\right)$ can be modified and written by:

To calculate the optimization of $F\left(x\right)$, the optimization of penalty function $\varphi \left(v\right)$ is calculated firstly. Not that for all $v$ and $u\ne$ 0 with equality when $u=v$, we have:

The Eq. (9) is used for each group, thus the optimization of $\varphi \left(v\right)$ is presented as follows:

with $g(v,u)\ge \varphi \left(v\right)$, $g(u,u)=\varphi \left(u\right)$. Note that $g(v,u)$ is quadratic equation on $v$. The Eq. (10) can be written as:

where $c$ does not depend on the value of $v$, and $\mathrm{\Lambda }\left(u\right)$ is a diagonal matrix, i.e., $\mathrm{\Lambda }\left(u\right)=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\left|v_k\right|\right)$. After manipulation, we have:

Therefore, by using Eq. (11), the optimization of $F\left(x\right)$ can be written as:

i.e., $G(x,u)\ge F\left(x\right)$, $G(u,u)=F\left(u\right)$.

In order to minimize the function $F\left(x\right)$, the majorization-minimization (MM) algorithm [23] defines an iterative algorithm via:

where $i$ is iteration exponent. The iteration is initialized with $x^{\left(0\right)}$. Thus the Eq. (14) is written as:

Solving the above equation, we have:

Note that the problem with update Eq. (16) is that as the iterations progress, some values of $\mathrm{\Lambda }\left(D{x}^{\left(k\right)}\right)$ will generally go to zero and becomes sparse, and therefore some entries of ${\mathrm{\Lambda }}^{-1}\left(D{x}^{\left(k\right)}\right)$ in Eq. (16) will go to infinity, leading the result is not accurate. We use the matrix inverse lemma from Ref. [24] to address this issue. Generally, the matrix inverse lemma has several forms, a common form is:

Thus, the $[I+\lambda D^T\mathrm{\Lambda }(D{x}^{\left(i\right)})D{]}^{-1}$ in Eq. (16) can be written as:

Finally, using Eq. (18), the update Eq. (16) can be written as:

Therefore, the group-sparse total variation denoising (GS-TVD) problem can be solved by the iterative algorithm in Eq. (19).

In order to investigate and compare the denoising performance of GS-TVD and general TVD method, Fig. 3 illustrates group-sparse TV denoising on a blocks signal. The test block signal and its corrupted signal are provided in Fig. 3(a). Fig. 3(b) and (c) respectively illustrate the denoising performance of GS-TVD and general TVD method. Moreover, denoising errors are shown in Fig. 3(d). By comparison, the simulation results indicate that the significant improvement in denosing is obtained with GS-TVD approach compared to the TVD method. Here, the group size was set to $k=$ 2, the regularlization parameter is 2, and the number of iteration is 20.

Fig. 3Comparison of simulation results between GS-TVD compared to the TVD method

a) The test signal and noisy test signal

b) The test signal and denoising signal using TVD method

c) The test signal and denoising signal using GS-TVD method

d) Denoising error using TVD and GS-TVD method

4.1. Case 1 – experimental setup

The bearing fault vibration signals were generated by the NSFI/UCR Center for Intelligent Maintenance Systems (IMS) [25, 26]. The experimental setup is shown in Fig. 5. The rotation speed is 2000 rpm and the sampling rate is set at 20 kHz. Vibration acceleration signals were collected every 10 mins by a NI 6062E-DAQ-Card. The accelerated life test experiment was carried out successively for 8 days until the magnetic plug exceeds a certain level and causes an electrical switch to close. Meanwhile, the severe wear failure of outer race was detected in bearing 1. The geometric parameters and ball pass frequency outer (BPFO) race of the tested bearing are presented in Table 1.

Fig. 6 shows the Kurtosis curve over the whole life-cycle of bearing 1 and shows that there is a long time in stable or normal operation for the bearing in whole life-cycle and the period of fault occurrence and severity is relatively short. As shown in Fig. 6, there is an obvious transient feature for the incipient fault at the stage of point 647.

Fig. 5Experimental setup for bearing accelerated life test [22, 23]

a)

b)

Fig. 6The Kurtosis curve of whole life-cycle of the rolling bearing 1

However, due to the interference of environment and background noises, the engineers are not sure whether the fault is happened before point 647 or not. Hence, to verify the effectiveness of the proposed method for bearing incipient fault diagnosis, the experimental data at point 535 was chosen which has no obviously wave phenomenon during the whole life-cycle.

4.2. Results and discussion for Case 1

The original vibration signal (2048 sampling points is selected, namely, approximately 0.1 s), the amplitude spectrum and its corresponding envelope spectrum of outer-race at point 535 are displayed in Fig. 7(a), (b) and (c), respectively. From Fig. 7(a), the cyclic periodicity characteristic of time domain waveform relevant to bearing fault cannot be discovered. Moreover, as shown in Fig. 7(c), although the spectrum peak at 230 Hz which consists with the outer-race fault frequency can be detected without denosing, however, the spectrum peak masked by heavy background noise and features are not be evident enough to detect fault.

The proposed method is then employed to analyze the outer-race vibration signal. First, for parameter $Q$ setting, it can be seen from the amplitude spectrum (see Fig. 7(b)) that the spectrum of the shock component is wider, and if the value of $Q$ is selected higher, due to that the noise component will increase and the underlying wavelets will have more oscillations with a narrower frequency responses, hence, we set the scope of resonant quality factor is $Q\in \left[1,2\right]$ with the interval $∆Q=$ 0.1. The maximal spectral kurtosis TQWT method is applied on the original vibration signal, then the kurtosis values under different scales-$j$ are obtained. The Kurtosis surf plot under different $Q$-factor and decompose level is displayed in Fig. 8(a). It can be seen from Fig. 8(a), the maximum index of kurtosis 4.709 is emerged in $Q$-factor $Q=$ 1.55 and level $j=$ 3.5, as marked by the black rectangle, which means that this area contains much more transient impulse components information than other $Q$-factor and levels areas. Applying the inverse TQWT method to this area and the corresponding reconstructed fault signal is shown in Fig. 8(b). Obviously, the reconstructed signal illustrated in Fig. 8(b) shows that the background noise is suppressed relatively. However, the low-frequency noise components still exist, mainly because of the range of frequencies available band. By computation, the minimum frequency is:

and the maximum frequency is$f_H={\alpha }^{j-1}{f}_s/2\approx$ 4688 Hz, thus the impact component with frequency 4252 Hz (see Fig. 7(b)) and high-frequency noise were contained within the scope [1011 Hz, 4688 Hz]. This conclusion can be verified by the corresponding envelope spectrum of the reconstructed signal as illustrated in Fig. 8(c). From Fig. 8(c), although the outer race fault frequency ($f_{BPFO}$) can be identified, several unrelated frequencies around $f_{BPFO}$ still remain in frequency domain and also the $f_{BPFO}$ information are not very obvious.

Fig. 7Original experimental signal, amplitude spectrum and envelope spectrum

a) Original experimental signal

b) Amplitude spectrum of the original bearing vibration signal

c) Envelope spectrum of the original bearing vibration signal

Furthermore, the group-sparse TVD method is employed to process the reconstructed bearing vibration signal, here, the group size is set to $k=$ 3, the regularization parameter $\lambda =$ 2 and the iteration time is 300. The filtered signal with group-sparse TVD method and its envelope spectrum are illustrated in Fig. 9(a) and Fig. 9(b), respectively. It can be observed that the noise components in reconstructed vibration signal have been removed evidently as shown in Fig. 9(a), and the cyclic periodicity characteristic is prominent in the filtered vibration signal. From Fig. 9(b), an obvious spectrum line 230.1 Hz, which is close to the $f_{BPFO}$ frequency according to the Table 1, and its harmonics are exposed in the envelope spectrum, meanwhile, the amplitudes of fault frequencies illustrated in Fig. 9(b) are evidently enhanced compared with that in Fig. 8(c), thus making the feature extraction becomes easier and accurate. Moreover, besides the outer race fault frequency, it is observed that the spectrum peaks at 70 Hz and 100 Hz are also identified in Fig. 9(b), which consists with the double shaft rotating frequency (70 Hz ≈ 2$f_r=$ 33.33×2) and its triple frequencies (100 Hz ≈ 3$f_r=$ 33.33×3). The above results illustrate that the fault characteristic frequencies can be detected by the proposed method even if the useful failure information is extremely weak.

Fig. 8Kurtosis surf plot under different Q-factor and decompose level, and the reconstructed vibration signal from the maximum peak

a) Kurtosis surf plot under different $Q$-factor and decompose level

b) The reconstructed vibration signal

c) The envelope spectrum of the reconstructed vibration signal

Fig. 9The filtered bearing signal and its envelope spectrum processed by proposed method

a) The filtered signal with GS-TVD method

b) The envelope spectrum of the filtered signal

To verify the superiority of the proposed approach, the Spectral Kurtogram (SK) combined with EEMD method is employed to analyze the original vibration signal. The SK is a powerful analysis tool to detect transients feature from their noisy vibration signals and has been widely studied and applied in the rotating machine diagnosis, see refs [27] and [28]. First, we apply the EEMD method [29, 30] to decompose the original vibration signal, thus the 12 IMF components are generated and displayed in Fig. 10, which can be helping to distinguish the periodic impulse from the mixed noisy signal. Then, we calculate the Kurtosis values of all the IMFs model orderly, see Table 2, and 2nd IMF model that contains the most fault information according to its maximum kurtosis value.

The Kurtogram of 2nd IMF model is displayed in Fig. 11(a), from which the optimal demodulation frequency band namely [0 Hz-2560 Hz] can be detected, as marked by the black rectangle. Thus, a band-pass filter is designed to extract the potential features from the 2nd IMF model. Last, the envelope spectrum is applied to the filtered signal, the corresponding envelope spectrum is illustrated in Fig. 11(b). As can be seen, the fault characteristic frequencies $f_i=$ 236.4 Hz cannot be detected in the envelope spectrum and it is also hard to distinguish the fault location from the incipient vibration signal. In conclusion, the comparison above re-confirms the efficiency of the proposed approach.

Fig. 10Sub-band components of original signal decomposed by using EEMD

a) IMF1-IMF6

b) The Hilbert envelop spectrum of IMF1-IMF6

c) IMF7-IMF12

d) The Hilbert envelop spectrum of IMF7-IMF12

4.3. Case 2 – experimental setup

The bearing vibration signals were generated by the Curtin University in Australia [31]. The experimental setup (Machinery Fault Simulator, MFS) is shown in Fig. 12. The rotation speed is 29 Hz, sampling frequency is 51200 Hz and length of record is 10 seconds. The bearing information and the fault frequencies of inner race and outer race are summarized in Table 3.

In this work, the inner race fault signal was chosen as an example for the subsequent analysis, and the theoretical inner race fault frequency is 103.59 Hz based on the equation $f_i=\frac{n}{2}\left[1+\frac{BD}{PD}\mathrm{c}\mathrm{o}\mathrm{s}\beta \right]\frac{RS}{60}$, where $BD$ and $PD$ are ball diameter and pitch diameter of the rolling bearing, respectively. Parameter $\beta$ represents contact angle ($\beta =$ 0), $n$ denotes the number of elements, $RS$ represents shaft rotating speed ($RS=$ 29×60 $=$ 1740 rpm).

Fig. 11Kurtogram of second IMF model component and its frequency spectrum

a) Kurtogram of second IMF model component

b) The envelope spectrum of band-pass filtered signal

Fig. 12Experimental setup

4.4. Results and discussion for case 2

In this experiment, the inner race vibration fault signal (5120 points), i.e. 0.1 second data was selected. In order to get a real early failure, a Gaussian white noise with amplitude $A=$ 17 (a high amplitude noise was added for highlighting the effect of proposed method) was added to the original signals. The time domain waveform, amplitude spectrum and its corresponding envelope spectrum are shown in Fig. 13(a), (b) and (c), respectively. However, it can be found that the periodical pulses are submerged in heavy noise in the time domain signal, and the inner fault $f_{BPFI}$ cannot be identified from the envelope spectrum and fault type cannot be determined yet.

The proposed approach is applied to the experimental signal, the algorithm parameters of maximal spectral kurtosis TQWT method are the same as the previous section. The Kurtosis surf plot under different $Q$-factor and decompose level, and the reconstructed bearing vibration signal are shown in Fig. 14. From Fig. 14(a), the maximum index of Kurtosis 12.63 is emerged in $Q$-factor $Q=$ 1.55 and level $j=$ 7.5, as marked by the black rectangle. The fault- induced impulses with cyclic periodicity characteristic can be identified from the reconstructed vibration signal.

Fig. 13Experimental signal of inner race failure, amplitude spectrum and its envelope spectrum

a) Experimental signal of inner race failure

b) Amplitude spectrum

c) Envelope spectrum

Fig. 14Kurtosis surf plot under different Q-factor and decompose level, and the reconstructed vibration signal from the maximum peak

a) Kurtosis surf plot under different $Q$-factor and decompose level

b) The reconstructed vibration signal

Fig. 15The reconstructed bearing signal and its envelope spectrum processed by proposed method

a) The filtered signal with GS-TVD method and its wavelet time-frequency diagram

b) The envelope spectrum of the filtered signal

Similarly, the reconstructed bearing vibration signal is also processed by the group-sparse TVD method, and the group size is set to $k=$ 2, the regularlization parameter is $\lambda =$ 2 and the iteration time is 300. The filtered signal with GS-TVD method and its wavelet time-frequency diagram, and its envelope spectrum are displayed in Fig. 15.

Obviously, it can be observed that the GS-TVD algorithm not only extracts the transient impulse components but also the noise components in reconstructed vibration signal have been removed evidently. As shown in Fig. 15(a), the phenomenon of frequency modulation can be clearly observed by the wavelet time-frequency diagram, and the impact interval of 9.7 ms could be also identified via the wavelet time-frequency diagram and the time-domain waveform, which coincides with the inner race fault frequency. From Fig. 15(b), as can be seen, the characteristic frequency 100.2 Hz, which is agreeable with the theoretical fault frequency of inner race 103.59 Hz, and its harmonic frequencies are clearly detected, hence, the proposed approach has a very powerful effect on the incipient fault diagnosis.

Fig. 16Sub-band components of original signal decomposed by using EEMD

a) IMF1-IMF7

b) The Hilbert envelop spectrums of IMF1-IMF7

c) IMF8-IMF13

d) The Hilbert envelop spectrums of IMF8-IMF13

For a comparison, the EEMD result, i.e., 13 IMF components related to the inner race fault are shown in Fig. 16. Meanwhile, the Kurtogram is also calculated to assist selecting the optimal frequency band for envelope demodulation analysis, as illustrated by Fig. 17(a). Based on highest SK value in Kurtogram, see Table 4, an optimal frequency band namely [0-1066.66 Hz] is selected for the 4th IMF model, as marked by the black rectangle, then the sensitive component within potential fault frequency is separated via a band pass filter. Unfortunately, the inner race fault characteristic frequencies $f_{BPFI}=$ 103.59 Hz cannot be detected in the envelope spectrum, see Fig. 17(b).

Obviously, although the double inner fault frequency (210 Hz) can be identified, the peak amplitude is lower than that of the proposed approach. Therefore, the comparison analysis demonstrates that the proposed maximal spectral kurtosis TQWT and group sparsity TVD approach can not only provide better de-nosing performance but also has superiority in enhancing performance for incipient fault diagnosis.

Fig. 17Kurtogram and its frequency spectrum

a) Kurtogram of 4th IMF model component

b) The envelope spectrum of band-pass filtered signal