The enhancement of fault detection for rolling bearing via optimized VMD and TQWT based sparse code shrinkage

2.1. Theoretical basis of VMD

As a non-recursive signal decomposition method, VMD has its solid theoretical foundation, which is the generalization of classical Wiener filter in multiple adaptive frequency bands. VMD decomposes the input signal into a series of sub-modes $u_k$ with specific sparse characteristics. In frequency domain, the bandwidth of each sub-mode $u_k$ is compact near the frequency center, and the bandwidth is estimated by the $L^2$ norm of the gradient. The resulting constrained variational problem is the following [19]:

where $\left\{u_k\right\}=\left\{u_1,u_2,\dots u_k\right\}$ and $\left\{{\omega }_k\right\}=\left\{{\omega }_1,{\omega }_2,\dots {\omega }_k\right\}$ are the set of sub-modes and their center frequencies, respectively. Equally, $\delta \left(t\right)$ is the impulse function and $k$ is the decomposition level.

In VMD algorithm, the quadratic penalty factor $\alpha$ and Lagrange multiplication operator $\lambda$ are introduced to solve the constrained variational problem. The quadratic penalty factor can ensure the reconstruction accuracy of the signal, and the Lagrange multiplication operator can ensure the accurate execution of the constraints. The enhanced Lagrangian operator $L$ can be described as [19]:

The original minimization problem Eq. (1) is transformed into problem Eq. (2), which is solved by alternate direction method of multipliers (ADMM). Firstly, the number of decomposition modes is determined, and the sub-modes ${\widehat{u}}_k^1$, the corresponding center frequency ${\omega }_k^1$ and Lagrange operator ${\lambda }^1$ are initialized. Then, the sub-mode ${\widehat{u}}_k$ and the center frequency${\omega }_k$are updated by Eq. (3) and Eq. (4), respectively:

The Lagrange operator is updated by Eq. (5):

where $\tau$ denotes the tolerance parameter of noise.

While the following convergence condition is satisfied, the iteration is terminated:

where $\epsilon$ is the convergence error. The noise tolerance $\tau$ and convergence error $\epsilon$ have little influence on the decomposition result, so the default values are usually used.

2.2. Parameter optimization method

As a parametric signal decomposition method, VMD involves two key parameters, which have a crucial impact on the decomposition results. A key parameter is the decomposition level $k$. If decomposition level is too large or too small, it is not conducive to the signal decomposition, thus affecting the recognition accuracy of information features. Another key parameter is penalty factor $\alpha$. The penalty factor determines the bandwidth of each mode component. To obtain the impulse characteristics of bearing fault as much as possible, we synthesize multiple dimensionless indexes to form a composite characteristic index, so as to optimize and screen two key parameters. In other words, the largest average composite index of mode components corresponds to two key parameters of optimization. Three dimensionless indexes are selected, which are kurtosis $K$, peak factor $P$ and margin factor $M$. These three dimensionless indexes measure the impulse characteristics from different aspects, and can well obtain the optimal mode decomposition of the VMD. To maintain the balance between indicators and increase the comparability, we use the relative ratio to construct the composite index $C$:

where $K_u$, $P_u$ and $M_u$ are the corresponding indexes of VMD decomposition components, and $K$, $P$ and $M$ are the corresponding indexes of original signal respectively.

After the optimization parameters are determined, the sensitive component of VMD is then judged according to the kurtosis maximization criterion. It is easy to see that the extracted sensitive component contains more bearing fault information, which is conducive to the subsequent impulse feature enhancement.

3.1. Tunable Q-factor wavelet transform

According to the above analysis, the optimized VMD can adaptively extract the mode components with more impulse information, but the extracted impulse will still contain a lot of noise interference. It is well known that impulse signal has non-Gaussian statistical characteristics, while noise in vibration signal is generally considered as Gaussian distribution. Therefore, a new TQWT based sparse code shrinkage algorithm is proposed for sparse representation of the mode component to further highlight the periodic impulses.

The difference between TQWT and dyadic wavelet transform is that TQWT realizes $Q$-factor adjustment and redundant operation by iterating filter banks through low pass scale factor $\alpha$ and high pass scale factor $\beta$, shown in Fig. 1. This design concept can not only extract the oscillation mode, but also obtain the detailed characteristics, which is very suitable for matching the attenuation oscillation mode of rolling bearing fault signal.

Fig. 1Analysis and synthesis filter bank for the tunable-Q wavelet transform

a) Analysis filter bank

b) Synthesis filter bank

TQWT directly designates quality factor $Q$ and redundancy factor $r$ to design wavelet, which further increases the flexibility of quality factor selection and makes wavelet acquisition more convenient. After $Q$ and $r$ are selected, $\alpha$ and $\beta$ can be obtained by [20]:

The expressions of low pass filter and high pass filter are as follows:

where $\theta \left(\omega \right)=0.5\left(1+\mathrm{c}\mathrm{o}\mathrm{s}\omega \right)\sqrt{2-\mathrm{c}\mathrm{o}\mathrm{s}\omega }$, $\left|\omega \right|\le \pi$.

Furthermore, the decomposition level $J$ is limited by the length of the signal $N$. Their relation can be expressed as:

When $0<\alpha \le 1$, the input signal is set as $X\left(\omega \right)$, and the low-scale filtering characteristic of the corresponding output signal $Y\left(\omega \right)$ is expressed as follows:

If $\alpha \ge 1$:

For the high scale, when $0<\beta \le 1$, the filtering characteristic is given by:

If $\beta \ge 1$:

According to the frequency domain characteristics of high-low-pass filter and its subsequent high-low-pass scale transformation, combining with the iterative operation of the basic filter banks under the multi-scale decomposition, the equivalent frequency response function for the $j$th stage is given by [20]:

3.2. Sparse code shrinkage

The SCS algorithm [27] proposed by Hyvarinen uses the statistical characteristics of non-Gaussian components to get the threshold shrinkage function, and denoises the measured signal by removing the noise with the same or similar frequency as the target signal. The sparse representation of bearing vibration signal can be realized by SCS algorithm, because the fault-generated impulse is strongly non-Gaussian in statistical characteristic. Actually, the local damage vibration signal of rolling bearing is characterized by Non-Gaussian property, repetitive transient impulse waveform and typical sparseness.

Assume $x_f$ is the original fault signal and $v$ is Gaussian noise of zero mean and variance ${\sigma }^2$, then the observed signal $x$ is given by:

To represent a sparse distribution, Hyvarinen proposes the following probability density function of a sparse signal [27]:

where $d$ is the standard deviation of fault impulse signal $x_f$ and $\alpha$ is a parameter controlling the sparseness of the probability density function. In this article, we set $\alpha$ = 0.1 to adapt bearing fault feature extraction [5].

For the sparse signal with this distribution, the following thresholding rule can be obtained according to the maximum likelihood principle [27]:

where $a=\sqrt{\alpha \left(\alpha +1\right)/2}$, $\sigma$ is the standard deviation of the noise, and $x_f$ is set to zero in the case that the square root in above equation is imaginary. The SCS algorithm is used as TQWT threshold and the noise standard deviation is estimated for each scale [28-29]:

where $mad$ denotes the median absolute deviation of wavelet coefficients.

Assume $x$ is zero mean, for $x$ and $v$ are uncorrelated, the standard deviation $d$ of $x_f$ can be estimated by [28-29]:

It is worth noting that the standard deviation $d$ needs to be calculated separately on each wavelet scale.

3.3. Adaptive sparse representation

As a new explicit wavelet construction theory in frequency domain, TQWT has the advantages of matching specific oscillation behavior of signal components and fast implementation by FFT algorithm. SCS is a useful tool for signal sparse representation. Naturally, the TQWT along with SCS can efficiently implement multi-scale sparse decomposition of dynamic signals.

For TQWT, as mentioned before, three parameters need to be determined: the $Q$-factor, the decomposition level $J$ and redundancy $r$. The specified value of the redundancy must satisfy that $r>1$. In this work, the redundancy $r$ is set to 3 by considering the translation invariance and complex calculations. The decomposition level $J$ only affects the frequency domain decomposition performance in low frequency region. Too many decomposition levels will lead to high computational cost, and may lead to excessive decomposition or redundant decomposition of fault characteristic frequency band information. The selection of $J$ should meet the requirements of Eq. (11). $Q$ is an important parameter, which directly affects the extent of wavelet sustained oscillation. We know that the bearing fault signal has local attenuation oscillation behavior. Therefore, the selection of $Q$ has an important influence on matching the bearing fault signal.

The signal processed by TQWT based sparse code shrink has obvious sparsity, and the periodic impulse characteristics are significantly enhanced. Since cross-correlation coefficient directly reflects the similarity of two signals, it is utilized to adaptively determine $Q$ and $J$. Cross-correlation coefficient $C$ is formulated as:

where $x\left(n\right)$ is the TQWT reconstruction signal in optimization process, and $u\left(n\right)$ is the sensitive component of optimized VMD.

The higher the degree of similarity and matching between the two signals is, the greater the $C$ value is. In general, $C$ can not only ensure that the detected impulses are not lost, but also ensure that the local attenuation oscillation waveform matches well. Therefore, the maximized $C$ is selected as the optimization objective, and the defect occurrence can be effectively judged through the analysis of the time domain waveform and envelope spectrum under these parameters.

5.1. Experimental verification

The test rig (ABLT-1A) is shown in Fig. 6. ABLT-1A can test four bearings at a time. Four bearings are installed on one shaft. The four test bearings are type 6309. The rotation speed is kept constant at 3000 r/min. The vibration signal is measured by the probe sensor directly contacting the bearing outer race, and the sampling frequency is 48K. The structural parameters of 6309 are as follows: pitch diameter – 72.5 mm, ball diameter – 17.462 mm, ball number – 8. Thus, the fault frequencies of outer race and inner race are $f_o=$151.829 Hz and $f_i=$248.171 Hz, respectively.

Due to the large noise in test process, the vibration signals are disordered, and the impulse characteristic signals are almost completely submerged, so the useful fault information cannot be obtained. The proposed optimized VMD ($k$ = 3 and $\alpha$ = 200) approach is applied to decompose outer race fault signal from the noisy observation, as shown in Fig. 7. According to kurtosis criterion, $U_2$ is automatically identified as sensitive component.

Fig. 6ABLT-1A test rig

Fig. 7Outer race fault signal and decomposition components of the optimized VMD

Fig. 8Processing results by TQWT based sparse code shrinkage: a) time domain waveform, b) envelope spectrum, c) envelope spectrum of original signal

Furthermore, the processing results by using TQWT ($Q$ = 1.8 and $J$ = 9) based sparse code shrinkage is illustrated in Fig. 8(a) and (b). It can be observed that the useful fault features of outer race can be detected clearly. As expected, the repeated impulse period 1/$fo$ are clearly revealed and the fault frequency $f_o$ and its harmonic components 2$f_o$, 3$f_o$, 4$f_o$, 5$f_o$ are also clearly displayed in envelope spectrum. Although the envelope spectrum of original signal shows significant magnitude at fault frequency $f_o$, many harmonic components are lost, which is disadvantageous to fault identification, as displayed in Fig. 8(c).

Fig. 9Inner race fault signal and decomposition components of the optimized VMD

Fig. 10Processing results by TQWT based sparse code shrinkage: a) time domain waveform, b) envelope spectrum, c) envelope spectrum of original signal

The proposed approach is adopted to analyze inner race fault signal. Inner race fault signal and decomposition components of the optimized VMD ($k$ = 4 and $\alpha$ = 200) are illustrated in Fig. 9. Similarly, U3 is automatically identified as sensitive component. The processing results by TQWT ($Q$ = 2.6 and $J$ = 6) based sparse code shrinkage are displayed in Fig. 10. It can be seen that the result is sparse and has periodic structure. It should be noted that when there is a fault on the inner race surface, the fault frequency $f_i$ will be modulated by the rotating frequency $f_r$ because the damage area is rotating with the rotating shaft. This phenomenon is well reflected in Fig. 10(a) and 10(b), which may improve the accuracy of decisions in fault detection. Nevertheless, it is not easy to get useful fault information through envelope spectrum directly as shown in Fig. 10(c).

From the above experimental signal analysis results, it can be seen that the proposed approach is an effective tool for sparse enhancement and weak impulse extraction.

5.2. Comparative studies

As one of the powerful tools for detecting cyclostationarity signals, spectral correlation can effectively identify the periodic impulses caused by local damage of rolling bearings [7-9]. Fig. 11 illustrates the processing results of the spectral correlation. Although the fault frequency of outer race can be identified in the enhanced envelope spectrum, the effect is far from satisfactory because the harmonic components of fault frequency are not well reflected. In fact, due to the influence of the noisy working environment, the vibration transmission path and the variation of workload, the periodic impulses in vibration signals are overwhelmed by heavy background noise, which makes it difficult to effectively extract the impulse features.

Fig. 11Processing results of vibration signal by Fast-SC

Fig. 12Processing results of vibration signal by MED

In addition, the MED method [11] is also used to process the vibration signal, and the processing results obtained are shown in Fig. 12. Although there is a rich literature by using the MED, its use is still limited for machine condition monitoring and fault diagnosis because the MED is prone to detect dominant impulses, which are usually unrelated with the information of interest. It can be observed that Fig. 12 can only clearly display the fault frequency of outer race; and it is difficult to distinguish the harmonic components of fault frequency from background noise similar to spectral correlation.

Through the above comparative analysis, it is found that compared with the spectral correlation and MED, the presented method can increase the accuracy of periodic impulse features extraction and is beneficial to the fault diagnosis of rolling bearing.