A fault feature extraction method for single-channel signal of rotary machinery based on VMD and KICA

2.1. VMD

VMD is a new signal decomposition method, which is the process of solving the variational problem based on classical Wiener filter, Hilbert transform and frequency mixing. In this method, the frequency center and bandwidth of each decomposition component are determined by searching the optimal solution of the variational model, and the signal can be adaptively decomposed into a component with sparsity. It is assumed that each mode is a limited bandwidth with a central frequency, the central frequency and bandwidth are constantly updated during the decomposition process, and VMD is finding the mode function $f_k\left(t\right)$, $k\in \{\mathrm{1,2},\cdots ,K\}$ with the smallest sum of $K$ estimation bandwidths, and the sum of the mode functions is the input signal $S$. The bandwidth of each mode function is determined through the following steps:

(1) In order to obtain the analytic signal of the mode function, Hilbert transform is performed for each mode function $f_k\left(t\right)$, that is:

where $t$ is time, $\delta \left(t\right)$ is impact function, $\left\{f_k\right(t\left)\right\}=\left\{f_1\right(t),f_2(t),\cdots ,f_K(t\left)\right\}$ is $K$ components.

(2) Modulate the spectrum of each mode to the corresponding baseband, that is:

where $\left\{{\omega }_k\right\}=\{{\omega }_1,\cdots ,{\omega }_K\}$ is the center frequency of $f_k\left(t\right)$.

(3) Estimate the bandwidth of each mode component. The corresponding constraint variational model expression is as follows:

In order to obtain the above constraint variational problems, the penalty factor $\alpha$ of the quadratic term and the Lagrange multiplication operator $\lambda \left(t\right)$ are introduced. $\alpha$ is a large enough positive number, which can guarantee the reconstruction precision of the signal under the presence of Gaussian noise. The Lagrange operator makes the constraint conditions strict, and the extended Lagrange expression is:

The above variational problem is obtained by using alternate direction method of multipliers (ADMM), and the saddle points of the extended Lagrange expression can be obtained by alternately updating $f_k^{n+1}\left(t\right)$, ${\omega }_k^{n+1}$ and ${\lambda }^{n+1}\left(t\right)$. $f_k^{n+1}\left(t\right)$ is the mode function at the $n+1$th cycle, and ${\omega }_k^{n+1}$ is the center of the power spectrum of the current modal function. ${\lambda }^{n+1}\left(t\right)$ is the multiplication operator of the $n+1$th cycle. The solution process of $f_k^{n+1}\left(t\right)$ can be expressed as:

where $i\in \left\{\mathrm{1,2},\cdots ,K\right\}\wedge i\ne k$. By using the Fourier isometric transformation, Eq. (5) can be transformed to the frequency domain and replaced $\omega$ with $\omega -{\omega }_k$, and the result is transformed into the form of negative frequency interval integral. The solution of the optimization problem is:

According to the same process, get the center frequency update mode:

where ${\widehat{f}}_k^{n+1}\left(\omega \right)$ is equivalent to the Wiener filtering of $\widehat{S}\left(\omega \right)-\sum _{i\ne k}{\widehat{f}}_i\left(\omega \right)$. The Fourier transform to $\left\{{\widehat{f}}_k\left(\omega \right)\right\}$, the actual part is $\left\{f_k\left(t\right)\right\}$.

The specific implementation process of the VMD method is as follows:

(1) Initialize mode function $\left\{{\widehat{f}}_k^1\right\}$, center frequency $\left\{{\omega }_k^1\right\}$, Lagrange multiplication operator, initial cycle times $n=$ 0.

(2) Execute the cycle $n=n+1$.

(3) Update $f_k^{n+1}$, ${\omega }_k^{n+1}$.

(4) Update ${\widehat{\lambda }}^{n+1}\left(\omega \right)$. ${\widehat{\lambda }}^{n+1}\left(\omega \right)={\widehat{\lambda }}^n\left(\omega \right)+\tau \left(\widehat{S}\right(\omega )-\sum _k{\widehat{f}}_k^{n+1}\left(\omega \right))\text{,}$ where $\tau$ is a time constant, usually taken as 0.

(5) Given discriminant accuracy $\epsilon >$ 0, repeat the above steps until the iteration stop condition is satisfied:

2.2. KICA

The mathematical description of independent component analysis (ICA) is: $\mathbf{X}=\mathbf{A}\mathbf{S}+\mathbf{N}$, where $\mathbf{S}=[s_1,s_2,\cdots ,s_N]$ is a vector of $N$ source signals, $\mathbf{X}=[x_1,x_2,\cdots ,x_M]$ is a vector of M mixed signals, The mixing matrix $\mathbf{A}$ is a $M$×$N$ dimensional matrix, $\mathbf{N}$ is noise vector. The meaning of ICA is that when the mixing matrix and source signals are unknown, we only determine the remove-mixing matrix $\mathbf{W}$ based on the observed data, and the output $\widehat{\mathbf{S}}=\mathbf{W}\mathbf{X}$ is the estimation of source signals.

KICA is not a simple nucleation of ICA, but a new ICA method. The idea of kernel technology is to use nonlinear mapping $\varphi :R^m\to R$, and map the nonlinear variable $y^i\in R^m\mathrm{}(i=\mathrm{1,2},\cdots ,N)$ of original input space into a kernel feature space $R$ to linearize it, and then analyse the mapped data in this feature space. Thus, the linear blind source separation in the space $R$ is equivalent to nonlinear blind source separation in the original space. One of the important characteristics of this technology is that kernel function can be used instead of inner product between two vectors to realize nonlinear transformation without specific form being considered. The kernel function of radial basis function (RBF) $K(x_i,y_i)=\mathrm{e}\mathrm{x}\mathrm{p}(-||x_i-y_i||/{\delta }^2)$ is chosen here.

KICA is characterized by the use of nonlinear functions as a contrast function in reproducing kernel Hilbert space (RKHS), the signal is mapped from low-dimensional space to the high-dimensional space, and this method is used to find out the minimum value of contrast function in this space. This function has a certain correlation with mutual information and has better mathematical properties. Moreover, this function space is suitable for various source signals. Therefore, compared with traditional ICA, KICA has better flexibility and robustness.

The KICA’s contrast function is constructed by calculating the correlation of a set of random variables directly. Let $F$ be a real vector function space, for simple calculations, $s_1$ and $s_2$ are two unary random variables of space $F$. Define the correlation coefficient ${\tau }_F$ of $s_1$ and $s_2$ as the maximum correlation coefficient between random variables $f_1\left(s_1\right)$ and $f_2\left(s_2\right)$:

${\tau }_F$ is also called a contrast function between random variables. Obviously, if $s_1$ and $s_2$ are independent, so, ${\tau }_F=$ 0. If space $F$ is large enough and ${\tau }_F=$ 0, $s_1$ and $s_2$ are also independent of each other. The procedure of KICA is as follows:

Input: Data vector $x_1$, $x_2$, …, $x_n$ and kernel function $k(x,s)$.

Step 1: Whiten data vector $x_1$, $x_2$, …, $x_n$.

Step 2: Use Cholesky decomposition to find the Gram matrix (${\mathbf{G}}_1$, ${\mathbf{G}}_2$, …, ${\mathbf{G}}_m)$ of original independent data, where $s_i=\mathbf{W}{x}_i$ ($\mathbf{W}$ is the remove-mixing matrix).

Step 3: Define ${\lambda }_H$ be the maximum eigenvalue of Eq. (10):

Step 4: Minimize $M_{{\lambda }_H}=-0.5\mathrm{l}\mathrm{o}\mathrm{g}{\lambda }_H$ for $\mathbf{W}$.

Output: $\mathbf{W}$.

This algorithm keeps running repeatedly between Step 2 and Step 4 until the convergence condition is satisfied, so that the remove-mixing matrix $\mathbf{W}$ can be obtained. According to $s=\mathbf{W}x$, for a set of observed data $x_1$, $x_2$, …, $x_n$, the original independent source signals can be estimated effectively through the remove-mixing matrix $\mathbf{W}$.

2.3. Correlation coefficient

The formula for the correlation coefficient ${\sigma }_{x{y}_i}$ of component signal to the source signal is as follows:

where $\mathbf{X}$ is the mixed source signal. ${\mathbf{Y}}_i$ is the $i$th component signal. $D\left(\mathbf{X}\right)$ and $D\left({\mathbf{Y}}_i\right)$ are the variance of signal $\mathbf{X}$ and ${\mathbf{Y}}_i$. $\mathrm{c}\mathrm{o}\mathrm{v}(\mathbf{X},{\mathbf{Y}}_i)$ is the covariance between the signal $\mathbf{X}$ and ${\mathbf{Y}}_i$.

2.4. Evaluation index $\mathit{\rho }$

The fault signal extracted from the real signal generally contains some noise. In the power spectrum of the extracted fault signal, the power of the fault frequency is usually the highest, so other frequencies can be considered as noise. If the power of noise is very close to or even higher than the power of fault frequency, the fault type cannot be accurately determined. Therefore, in order to compare the extraction effect of each fault diagnosis method, this paper proposes an evaluation index $\rho$:

where $P_f$ is the power of fault frequency. $P_s$ is the power of the noise of the highest amplitude. Therefore, the bigger the value of $\rho$, the better the effect of the extraction.

4.1. Faulty rolling bearing experiment

In order to verify the validity of the proposed method for faulty rolling bearing vibration signal, the test data from the bearing database of Case Western Reserve University are selected for analysis. The rolling bearing used in the test is SKF6205. The motor speed is 1800 rpm, that is, the rotation frequency $f_r=$ 30 Hz. The sampling frequency $f_s=$ 12000 Hz and the number of sampling points is 12000. The inner and outer rings of the rolling bearing are respectively machined with tiny pitting pits to simulate the fault of inner ring and outer ring. The fault frequency of inner ring can be expressed as:

The fault frequency of outer ring can be expressed as:

where $N$ is the number of roller, $d$ is the diameter of roller, $D$ is the pitch diameter of bearing, and $\alpha$ is bearing contact angle. According to Eq. (13) and Eq. (14), $f_{ic}=$ 162 Hz and $f_{oc}=$ 107 Hz.

The time-domain waveform and power spectrum of single-channel vibration signal of inner ring fault and outer ring fault mixed are shown in Fig. 2.

The fault features of inner ring and outer ring cannot be found from time-domain waveform and power spectrum in Fig. 2. Therefore, so as to verify the high efficiency of the proposed method, the effect of feature extraction of wavelets, EEMD+FastICA and the proposed method are compared as follows.

Fig. 2The time-domain waveform and power spectrum of the observation signal

A. Wavelets.

In this method, the observation signal is decomposed with six-layer wavelet decomposition, and then the power spectrum of each component can be obtained. The result is shown in Fig. 3. $f_{ic}$ and $f_{oc}$ are not extracted, so the evaluation index $\rho$ can be considered as $-\mathrm{\infty }$. Therefore, wavelets cannot extract the feature of fault signals of rolling bearing.

B. EEMD + FastICA.

In this method, the observation signal is decomposed by EEMD first, and each intrinsic mode function (IMF) can be obtained. Then, calculate the correlation coefficient of each IMF to the observation signal. IMFs with high correlation coefficient were selected to reconstruct the new observation signal. Finally, the fault signals can be extracted by FastICA. The power spectrum of the extracted fault signals is shown in Fig. 4.

In Fig. 4(b), the frequency of 107.7 Hz which is very close to $f_{oc}$ is extracted and its double frequency is also extracted. The frequency of 538.3 Hz is the noise of the highest amplitude. The frequency of 161.9 Hz which is very close to $f_{ic}$ can be seen in Fig. 4(a), but the frequency of 107.7 Hz is also obvious, and it can be considered as noise now. Therefore, mode mixing occurs at this time. According to Eq. (12), ${\rho }_{oc}=$ 4.687, ${\rho }_{ic}=$ 0.229.

Fig. 3The power spectrum of the decomposed signals after wavelets: a) the sixth layer of low-frequency component, b) the sixth layer of high-frequency component, c) the fifth layer of high-frequency component, d) the fourth layer of high-frequency component

a)

b)

c)

d)

Fig. 4The power spectrum of the extracted signals after EEMD + FastICA

a)

b)

C. The proposed method.

In this method, the observation signal is decomposed by VMD ($K=$ 4, $\alpha =$ 1000) first and each component signal can be obtained. Then, calculate the correlation coefficient of each component signal to the observation signal and their values are in Table 1. Component signals with high correlation coefficient (Threshold > 0.5) were selected to reconstruct the new observation signal. Finally, the fault signals can be extracted by KICA. The power spectrum of the extracted fault signals is shown in Fig. 5.

The frequency of 161.9 Hz which is very close to $f_{ic}$ can be seen in Fig. 5(a). The frequency of 60.06 Hz is the noise of the highest amplitude. In Fig. 5(b), the frequency of 107.7 Hz which is very close to $f_{oc}$ is extracted and its double frequency is also extracted. The frequency of 30.03 Hz is the noise of the highest amplitude. Furthermore, the fault features are very obvious. According to Eq. (12), ${\rho }_{oc}=$ 14.597, ${\rho }_{ic}=$ 1.046.

Fig. 5The power spectrum of the extracted signals after VMD + KICA

The evaluation index of the three methods are in Table 2. The $\rho$ of the proposed method is much bigger than other two methods. Therefore, the proposed method has a better effect of extraction. The result proves that the proposed method is efficient to extract the features of single-channel vibration signal of rolling bearing with the fault of inner ring and outer ring mixed.

4.2. Comprehensive experiment with faulty rolling bearing and faulty gear

The experimental device is shown in Fig. 6. The whole device is driven by a 550 W (220 V-50 Hz) AC motor and drives the shaft system with couplings. The shaft section between two bearing seats is equipped with a belt wheel, and the belt drives the active gear shaft of gear box. One acceleration transducer is installed vertically on the shell near the rolling bearing. There are two rolling bearings on the shaft system, and one of them is machined with tiny pitting pits on the outer ring to simulate the fault of outer ring. The sampling frequency $f_s=$ 8192 Hz and the number of sampling points is 8192. The motor speed is 850 rpm, that is, the rotation frequency $f_r=$ 14.2 Hz.

According to Eq. (14), $f_{oc}=$ 70 Hz. Destroy one tooth of the driving gear to simulate the fault of tooth breaking. The number of gear teeth $Z=$ 20. The number of the broken teeth $Z_b=$ 1. According to $f_c=Z\times f_r$, the gear mesh frequency$f_c=$ 283.3 Hz. When the gear is affected by broken tooth fault, the rotation frequency $f_r$ and its frequency multiplication will be the main features in frequency-domain, and the frequency $f_z=Z_b\times f_c$ will also exist. The time-domain waveform and power spectrum of single-channel mixed signal are shown in Fig. 7.

Fig. 6The experimental device

Fig. 7The time-domain waveform and power spectrum of the mixed signal

a)

b)

Fig. 8The power spectrum of the decomposed signals after wavelets: a) the seventh layer of low-frequency component, b) the seventh layer of high-frequency component, c) the sixth layer of high-frequency component, d) the fifth layer of high-frequency component

a)

b)

c)

d)

The fault features of outer ring and gear breaking cannot be found from time-domain waveform and power spectrum in Fig. 7. Therefore, the effect of feature extraction of wavelets, EEMD + FastICA and the proposed method are compared as follows:

A. Wavelets.

In this method, the mixed signal is decomposed with seven-layer wavelet decomposition, and then the power spectrum of each component can be obtained. The result is shown in Fig. 8. The fault features are not extracted, so the evaluation index $\rho$ can be considered as $-\mathrm{\infty }$.

B. EEMD + FastICA.

In this method, the mixed signal is decomposed by EEMD first, and each intrinsic mode function (IMF) can be obtained. Then, calculate the correlation coefficient of each IMF to the mixed signal. IMFs with high correlation coefficient were selected to reconstruct the new observation signal. Finally, the fault signals can be extracted by FastICA. The power spectrum of the extracted fault signals is shown in Fig. 9.

The frequency of 70 Hz which is very close to $f_{oc}$ can be seen in Fig. 9(a). The frequency of 20 Hz is the noise of the highest amplitude. In Fig. 9(b), the frequency of 14 Hz which is very close to $f_r$ and its double frequency are extracted, and the frequency $f_c$ is also extracted. The frequency of 72 Hz is the noise of the highest amplitude. According to Eq. (12), ${\rho }_{oc}=$ 0.146, ${\rho }_g=$ –0.220.

Fig. 9The power spectrum of the extracted signals after EEMD + FastICA

C. The proposed method.

In this method, the mixed signal is decomposed by VMD ($K=$ 6, $\alpha =$ 1000) first and each component signal can be obtained. Then, calculate the correlation coefficient of each component signal to the mixed signal and their values are in Table 3. Component signals with high correlation coefficient (Threshold > 0.5) were selected to reconstruct the new observation signal. Finally, the fault signals can be extracted by KICA. The power spectrum of the extracted fault signals is shown in Fig. 10.

Fig. 10The power spectrum of the extracted signals after VMD + KICA

The frequency of 70 Hz which is very close to $f_{oc}$ can be seen in Fig. 10(a). The frequency of 31 Hz is the noise of the highest amplitude. In Fig. 10(b), the frequency of 14 Hz which is very close to $f_r$ and its double frequency are extracted, and the frequency $f_c$ is also extracted. The frequency of 72 Hz is the noise of the highest amplitude. Furthermore, the fault features are obvious. According to Eq. (12), ${\rho }_{oc}=$ 0.517, ${\rho }_g=$ 0.197.

The evaluation index of the three methods are in Table 4. The $\rho$ of the proposed method is much bigger than other two methods. Therefore, the proposed method has a better effect of extraction. The results demonstrate that the proposed method is efficient to extract the fault features of single-channel vibration signal of complex machinery.

4.3. Faulty rotary shaft experiment

A faulty rotary shaft experiment is applied to verify the effectiveness of the proposed method. The schematic diagram of test system is Fig. 11 and the test rig and signal acquisition system is shown in Fig. 12.

Fig. 11The schematic diagram of test system

Fig. 12The test rig and signal acquisition system of faulty rotary shaft

In this experiment, the artificial misalignment fault and rub-impact fault are created. Only one eddy current displacement sensor is used to collect the mixed vibration signal of faulty shaft. The sampling frequency is 1000 Hz and the number of sampling points is 5000. The motor speed is 2000 rpm, that is, the rotation frequency $f_r=$ 33.3 Hz. According to experience, the feature of misalignment fault is that a higher amplitude of $\text{2}\times f_r$appears, even exceeding the amplitude of $f_r$ in frequency-domain. The feature of rub-impact fault is emerging 1/$n$ of $f_r$, where $n$ is equal to 2, 3, or 4. The time-domain waveform of the collected signal is shown in Fig. 13.

Fig. 13The time-domain waveform of the collected signal

In the proposed method, the collected signal is decomposed by VMD ($K=$ 4, $\alpha =$ 1000) first and each component signal can be obtained. Then, calculate the correlation coefficient of each component signal to the collected signal and their values are in Table 5. Component signals with high correlation coefficient (Threshold > 0.5) were selected to reconstruct the new observation signal. Finally, the fault signals can be extracted by KICA. The power spectrum of the extracted fault signals is shown in Fig. 14.

Fig. 14The power spectrum of the extracted signals after VMD + KICA

Both the frequency of 33.45 Hz and the frequency of 16.48 Hz which is equal to half of $f_r$ can be found in Fig. 14(a). It can be considered as the feature of rub-impact fault. In Fig. 14(b), the amplitude of 66.89 Hz is larger than 33.45 Hz, which can be considered as the feature of misalignment fault. Furthermore, the fault features are significantly obvious. It demonstrates that the proposed method has a good performance of the fault feature extraction of rotary shaft with the misalignment fault and rub-impact fault mixed.