Rolling bearing fault diagnosis based on variational mode decomposition and weighted multidimensional feature entropy fusion

2.1. VMD

The VMD algorithm is a non-recursive adaptive signal decomposition method, with its decomposition process shown in 0. The signal is decomposed into a number of AM-FM IMFs, and the process can be expressed as:

where $t$ is the current time; the phase ${\phi }_k\left(t\right)$ is a non-decreasing function of $t$; $u_k\left(t\right)$ is the $k$th decomopsed IMF component; $A_k\left(t\right)$ is the instantaneous amplitude and is a non-negative envelope function; $K$ is the number of decomposition.

The process of constructing the variational mode is as follows:

(1) The analytic signal of IMF is obtained by Hilbert Transform, and thus the one-sided spectrum of the signal is obtained.

(2) The exponential term is used to adjust the center frequency of each IMF and transform each IMF frequency to the fundamental frequency band.

(3) The bandwidth of each IMF component is calculated using the Gaussian smoothing demodulated signal. The constructed variational mode is:

where: $f\left(t\right)$ is the input signal; $\delta \left(t\right)$ is the unit impulse function; ${\omega }_k$ is the center frequency; $*$ is the convolution operation.

An extended Lagrangian function given by the following expression is introduced to solve the optimal solution of the variational mode:

where: $\alpha$ is the penalty factor; $\lambda \left(t\right)$ is the Lagrange multiplier.

The Alternate Direction Method of Multipliers (ADMM) is used to update $u_{k,n+1}$, ${\omega }_{k,n+1}$ and ${\lambda }_{n+1}$ to obtain the saddle point of the Lagrangian function, hence the optimal solution of Eq. (3). The Fourier Transform is used to update Eq. (3) from the time domain to the frequency domain:

where: is the fidelity factor; $n$ is the number of iterations; $\omega$ is the frequency; $\wedge$ denotes the Fourier Transform.

Repeat Eq. (4) through Eq. (6) until satisfying the following iteration stopping condition where the update can come to a stop:

where: $\epsilon$ is the discriminatory accuracy which takes 10^-6.

By the time the iteration stops, the frequency-domain characteristics of the signal have been decomposed adaptively, and the modulated signal ${\widehat{u}}_k\left(\omega \right)$ has been transformed into the time-domain IMF component using the inverse Fourier Transform.

Fig. 1Flowchart of VMD algorithm

2.2. WMFE

Based on the RCMSE, three entropy methods, RCMDE, RCMFE and RCMPE, have been established in the fields of biomedical engineering and mechanical fault diagnosis. All these methods have the same coarse-grained process with a scale factor of 3, as shown in 0, Suppose a time series $\{x_i,i=\mathrm{1,2},...N\}$, where $N$ is the total length of the signal. Given a scale factor$k$, $\left(k=\mathrm{1,2},...,\tau \right)$, the $k$th coarse-grained sequence $y_k^{\left(\tau \right)}=\{y_{k,1}^{\tau },y_{k,2}^{\tau },...,y_{k,(i+1)/2}^{\tau }\}$ of $x_i$ is given by:

The steps to calculate the RCMDE are described as follows.

(1) Let $k=$ 1. The DE of the time series $y_1^{\left(\tau \right)}$ is calculated as follows.

a) Map the coarse-grained time series $y_1^{\left(\tau \right)}$ into the class $c$ ($c$ is a positive integer). This process needs to be implemented in two substeps: First, map the coarse-grained series $y_1^{\left(\tau \right)}$ onto the interval (0, 1) by Eq. (9) for the normal cumulative distribution function (NCDF):

where: $\mu$ and $\sigma$ denote the mean and standard deviation of $y_1^{\left(\tau \right)}$, respectively.

Next, map $b_j$ to class $c$ in the form of $s_j^c=R\left(c\cdot b_j+0.5\right)$ by linear transformation, where $R\left(*\right)$ is the rounding operation.

b) Given the embedding dimension $m$ and the time delay $d$ reconstruct the time series $s_i^{m,c}=\left\{s_i^c,s_{i+d}^c,\cdots ,s_{i+(m-1)d}^c\right\}$, where $i=\mathrm{1,2},...,N-\left(m-1\right)d$.

c) Assuming that each time series $s_i^{m,c}$ corresponds to a dispersion pattern, let $s_i^c=v_0$, $s_{i+d}^c=v_1$, …, $s_{i+(m-1)d}^c=v_{m-1}$, then $s_i^{m,c}$ corresponds to a dispersion pattern ${\pi }_{v_0{\nu }_1\cdots v_{m-1}}$, where the largest dispersion pattern does not exceed $c^m$.

d) The probability of each dispersion pattern is calculated as follows:

where the term $number\left({\pi }_{v_0{v}_1\cdots v_{m-1}}\right)$ indicates the number of $s_i^{m,c}$ corresponding to the dispersion pattern.

e) Therefore, DE can be calculated as:

(2) Repeat step (1) for $k=\mathrm{2,3},...,\tau$ to obtain DE at multi-scale. Thus, the RCMDE expression is given by:

where $\stackrel{-}{p}\left({\pi }_{v_0{v}_1\cdots v_{m-1}}\right)=\frac{1}{\tau }\sum _k^{\tau }{p}_k^{\tau }$.

The steps to calculate the RCMFE are described as follows.

(1) Let $k=$ 1. The FE calculation substeps for the time series $y_1^{\left(\tau \right)}$ are as follows.

a) Given the embedding dimension $m$ and time delay $d$, reconstruct the coarse-grained time series $y_1^{\left(\tau \right)}$ into the following time series $x_i^m=\left\{x_i^{},x_{i+d}^{},\cdots ,x_{i+(m-1)d}^{}\right\}-x_0\left(i\right)$, where $i=\mathrm{1,2},\dots ,N-\left(m-1\right)d$, $x_0\left(i\right)=\frac{1}{m}\sum _{t=0}^{m-1}x(i+td)$.

b) The maximum absolute distance between $x_i^m$ and $x_j^m$ is calculated using $d_{ij}^m=d\left[x_i^m,x_j^m\right]=\mathrm{m}\mathrm{a}{\mathrm{x}}_{v\in (0,m-1)}\left|\left[\left(x(i+v)-x_0\left(i\right)\right)-\left(x(j+v)-x_0\left(j\right)\right)\right]\right|$ where, $i,j=\mathrm{1,2},\dots ,(N-k+1)-m\mathrm{}$, $i\ne j$.

c) Defining the fuzzy weight as $n$ and the similarity tolerance as $r$, the similarity $R_{ij}^m$ between $x_i^m$ and $x_j^m$ is calculated by the fuzzy relationship function $\mu \left(d_{ij},n,r\right)=\exp \left(-{\left(d_{\frac{ij}{r}}\right)}^n\right)$.

d) Where by:

e) FE is obtained by the natural logarithm ratio of the continuous function in substep d):

(2) Repeat step (1) for $k=\mathrm{2,3},...,\tau$ at multi-scale to obtain FE. Thus, the RCMFE expression is given by:

where ${\stackrel{-}{\varphi }}^{m+1}(n,r)=\frac{1}{\tau }\sum _k^{\tau }{\varphi }^{m+1}(n,r)$, ${\stackrel{-}{\varphi }}^m(n,r)=\frac{1}{\tau }\sum _k^{\tau }{\varphi }^m(n,r)$.

Fig. 2Coarse-grained process of time series

The steps to calculate the RCMPE are described as follows.

(1) Let $k=$ 1. The PE of the time series $y_1^{\left(\tau \right)}$ is calculated as follows.

a) Given the embedding dimension $m$ and the time delay $d$, reconstruct the coarse-grained time series $y_1^{\left(\tau \right)}$ into the time series $x_i^m=\left\{x_i^{},x_{i+d}^{},\cdots ,x_{i+(m-1)d}^{}\right\}$, which is arranged in ascending order according to the numerical value of the elements. There are $m!$ possible permutations for the patterns, and for each permutation of pattern –$\pi$, the relative frequency is obtained as follows:

b) PE is calculated as follows:

Fig. 3Flowchart of WMFE algorithm

(2) Repeat step (1) for $k=\mathrm{2,3},...,\tau$ to obtain PE at multi-scale. Thus, the RCMPE expression is given by:

where $\stackrel{-}{p}\left(\pi \right)=\frac{1}{\tau }\sum _k^{\tau }{p}_k^{\tau }$.

The steps to calculate the WMFE are described as follows, with the flowchart of the algorithm shown in 0.

(1) The formula for calculating the weights using the standard deviation method is given by:

where, $i$ is RCMDE, RCMFE or RCMPE; ${\sigma }_i$ is the variance of the $i$th entropy.

(2) The variance formula for calculating entropy is given by:

where, $E_i$ is the entropy at the $i$th scale; ${\overline{E}}_i$ is the mean value of entropies.

Therefore, the WMFE expression is given by:

2.3. Algorithm parameter setting

2.3.2. Parameter setting of RCMDE, RCMFE and RCMPE

From Eqn. (12), four parameters are needed to calculate RCMDE: number of classes $c$, embedding dimension $m$, time delay $d$, and scale factor $\tau$. It is generally recommended that the value of $c$ be set to an integer between 4 and 8. For the embedding dimension $m$, although a larger value of m can reconstruct the dynamic process in more detail, an overly large value of $m$ will require many data points, while the length of most data in the real world is always finite; on the other hand, too small a value of $m$ may result in insufficient information. In this paper, $m$ is kept 2. To prevent information loss, the time delay $d$ is usually set as 1. Considering that too large a value of$\tau$increases the computational quantity, while too small a value of$\tau$falls short of extracting enough information for preventing the feature vector from growing too large and for ensuring a balance between performance and computational efficiency [34], in this paper the scale factor $\tau$ is kept 20, i.e., the RCMDE are extracted by a scale factor of 20 for data samples. In addition, the embedding dimension $m$ of RCMFE is the same as that of RCMDE, which is kept 2. And $n$ determines the width and gradient of the fuzzy affiliation function boundary. Too narrow a fuzzy affiliation function will make the final estimation of statistical properties inaccurate and sensitive to noise, while too wide a fuzzy affiliation function may lead to a loss of detailed information, so the value of $n$ is kept 2. The similarity tolerance $r$ determines the similarity of matching. The larger the value is, the more data information will be lost, while the smaller it is, the more sensitive it is to noise. And $r$ is generally recommended to be 0.1-0.25 times the SD (Standard Deviation) of the original data, and thus kept 0.15 SD in this paper [35]. Time delay $d$ and scale factor $\tau$ are the same for RCMDE. In calculating RCMPE, too large a value of $m$ makes it difficult to identify the dynamic changes in the time series, while too small a value of $m$ may impede the RCMPE from working due to the much smaller number of different states (symbols) [36]. In this paper, $m$ is kept 3, and the parameters of time delay $d$ and scale factor $\tau$ are the same as above.

Fig. 4Central frequency curves

2.4. Rolling bearing fault diagnosis based on VMD and WMFE Fusion

Based on the good robustness of VMD, end effect suppression, fewer model pseudo-components, and the advantages of WMFE being capable to extract multidimensional signal feature vectors adaptively, the fault diagnosis process is designed as shown in 0 The rolling bearing fault diagnosis based on VMD and WMFE fusion includes three major parts: VMD of the original signal, extraction of WMFE characteristic information, and CNN neural network classification. The procedure starts with decomposing the rolling bearing vibration signal under variable working conditions by VMD; next, RCMDE, RCMFE and RCMPE feature vectors are extracted from IMFs and weights are obtained using SDM, which are then assembled into WMFE data grid in matrix form and input into CNN. After the convolution layer of the multi-scale convolution kernel and the maximum pooling and dimension reduction of the multi-scale features in high dimension, the final input fully connected layer for classification by the Softmax activation function.

Fig. 5Flowchart of rolling bearing fault diagnosis method

3.1. Construction of a data set for rolling bearing simulation signals

Based on the structural characteristics of rolling bearings, a component of the bearing, when damaged in the bearing operation process where other components collide with each other, excites a shock signal in the form of high-frequency damped oscillation and an inherent vibration at high frequency. The vibration signal simulation model for rolling bearing fault diagnosis is expressed as follows [37]:

where: ${\tau }_k$ is the time fluctuation of the $k$th shock interval with respect to the shock period $T$; $s\left(t\right)$ is an exponentially decaying sinusoidal signal vibrating at the intrinsic frequency $f_n$; $C$ is the system damping ratio; $A_k$is the amplitude modulated signal;${\phi }_k$ and $c_k$ are random numbers; $\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{d}n\left(t\right)$ is a random signal with mean 0; $f_m$ is the modulation frequency, which depends on the fault type. In case of outer race fault, $f_m$ is 0; in case of inner race fault, $f_m$ is the axis rotation frequency; in case of ball fault, $f_m$ is the cage rotation frequency.

Due to the complexity of long-term high-speed working conditions and the special characteristics of their structure, rolling bearings are prone to compound faults, and multiple fault characteristics are superimposed on and interfere with each other, increasing the difficulty of feature extraction from compound faults. $x_1\left(t\right)$ (outer race), $x_2\left(t\right)$ (inner race), and $x_3\left(t\right)$ (ball) are the vibration simulation signals for single faults of rolling bearings, whereas $x_4\left(t\right)$ (outer race & inner race), $x_5\left(t\right)$ (outer race & ball), and $x_6\left(t\right)$ (inner race & ball) are the vibration simulation signals for compound faults of rolling bearings:

where: $C_o$ = 600, $C_i$ = 800, $C_b$ = 600, $f_{no}$ = 3000 Hz, $f_{ni}$ = 5000 Hz, $f_{nb}$ = 2000 Hz, $f_r$ = 25 Hz, $f_c$ = 5 Hz, $T_o$ = 1/50, $T_i$ = 1/90, $T_b$=1/40, sampling frequency $f_s$ = 1.6 kHz.

Figure 6 shows the temporal waveforms and FFT spectra of the 6 simulation signals to which the Gaussian white noise has been added with a signal-to-noise ratio (SNR) of –4 dB. The SNR equation is given as follows. According to in 0, it is impossible to distinguish fault types by directly observing the six temporal signals due to noise interference. Besides, from the FFT spectra in Fig. 6, the characteristic information of all other frequency bands but the resonance frequency band is quite similar, which means it is necessary to adopt an effective technique to process the information [38]:

where: $P_s$ is the input signal power, and $p_n$ is the noise power.

Gaussian white noise signals with distinct SNRs are added to $x_1\left(t\right)$, $x_2\left(t\right)$, $x_3\left(t\right)$, $x_4\left(t\right)$, $x_5\left(t\right)$, and $x_6\left(t\right)$, forming the simulation vibration signals for rolling bearing faults with SNRs ranging within [–4, 8]. The number of training sets and test samples at each $S/N$ ratio of these signals are 40 and 10, respectively, with sample length $N=$ 5000, summing up to 1680 training set samples and 420 testing set samples. The data sets are described in Table 2.

Fig. 6Temporal waveform and FFT spectrum of 6 simulation signals (SNR = –4 dB)

3.2. Comparative analysis of simulation

First, each rolling bearing simulation signal is decomposed by VMD into 4 IMF components. These IMFs with distinct center frequencies contain all the characteristic information of the simulation signal, as shown in Fig. 7. Then the RCMDE, RCMFE, and RCMPE of the IMFs are calculated via SDM to obtain the WMFE of the IMFs with entropy weights, as shown in Fig. 8, which correspond to the curves of different entropies of the 4 IMFs. It can be seen that the separability of the RCMDE, RCMFE, and RCMPE values of each IMF after weighting has been enhanced, which verifies that the WMFE method delivers higher accuracy than do the existing entropy methods in estimating the complexity of the signal at each scale. Compared with RCMDE, RCMFE, and RCMPE, as shown in Fig. 9, the recognition rate of the proposed WMFE method converges fast, up to 99.52 % at the third time and 100 % at the fifth time, while the recognition rates of RCMDE and RCMPE finally reach 99.29 % and 99.76 %, respectively; the effect of RCMPE is slightly worse, and its recognition rate finally reaches 99.29 %. Therefore, WMFE outperforms the other entropy methods in recognition rate, convergence speed, and noise immunity according to the simulation data.

Fig. 7The VMD of rolling ball fault simulation signal (ball fault –4 dB)

Fig. 8Curves of different entropies of 4 IMFs: a), b) IMF1, c), d) IMF2, e), f) IMF3, and g), h) IMF4

a)

b)

c)

d)

e)

f)

g)

h)

Fig. 9Accuracy of CNN training process by different methods

4.1. Description of dataset

As shown in Fig. 10, the experimental platform consists of four components: a motor, a torque transducer, a power tester, and a controller. The faulty bearing under test is the drive end bearing of model SKF6205 motor. The electric spark method was used to machine single-point grooves with damage diameters of 0.007 inches, 0.014 inches, and 0.021 inches on the surfaces of the inner race, ball and outer race, respectively, to simulate the wear process of the rolling bearing in real operation. At a sampling frequency of 12 kHz in this experiment, the acceleration data sets were collected at the speed of 1772 r/min, 1750 r/min, and 1730 r/min, corresponding to the load of 1 hp, 2 hp, and 3 hp, respectively; the collected data were divided into 10 types according to the location and degree of different damages. As shown in 0, 5000 sampling points are set for each segment, and 1,000 samples are selected randomly from the 10 types and divided into the training samples and test samples at a ratio of 4:1. The data set for the experimental variable working condition is described in Table 3.

Fig. 10Experimental apparatus for bearing: a) physical diagram, b) structural sketch

a)

b)

Fig. 11Time-domain waveforms of CWRU bearing data for 10 different states

4.2. Experimental comparative analysis

First, considering that VMD features adaptive signal decomposition, the rolling bearing vibration signal is decomposed by VMD into 4 IMF components. These IMF components with distinct center frequencies contain all the characteristic information of the rolling bearing vibration signal, as shown in Fig. 12.

Fig. 12The VMD of rolling ball fault signal(ball)

Then the RCMDE, RCMFE, and RCMPE of the IMFs are calculated via SDM to obtain the WMFE of the IMFs with entropy weights, as shown in Fig. 13, which correspond to the curves of different entropies of the 4 IMFs. It can be seen that the separability of the RCMDE, RCMFE, and RCMPE values of each IMF after weighting has been enhanced, which verifies that the WMFE method delivers higher accuracy than do the existing entropy methods in estimating the complexity of the signal at each scale. In this paper, the high-dimensional data are represented by the low-dimensional distribution of the T-SNE method. As shown in Fig. 14, the T-SNE visualizations of the WMFE, RCMDE, RCMFE, and RCMPE methods are presented for the 200 validation sets classified with the signal features extracted from them each. Evidently, the T-SNE visualization of WMFE for extracting features from rolling bearing signals has very distinctive classification features; the same categories are clustered at the same location, and the inter class distance is relatively much great; the 10 categories are independently separated from each other. In contrast, the classification results of the RCMDE, RCMFE and RCMPE methods are cluttered in the two-dimensional space and are much less effective than the WMFE method. The features extracted by the WMFE method are definitely far more sensitive than by other entropy methods, and thus the performance of WMFE is verified. This performance upgrade is owing to WMFE taking advantage of the complementary nature of the differences between different feature entropies and combining the advantages of the RCMDE, RCMFE and RCMPE methods with more effectiveness, to build a more comprehensive representation for the signal information of the high-dimensional feature set and for the characteristics of fault type information.

To verify the effectiveness of the method under real working conditions, the proposed method is compared with the RCMDE, RCMFE, and RCMPE methods, with the results shown in Fig. 15. The recognition rate of the WMFE method has converged to 98 % at the 3rd iterations, and stabilized to 100 % at the 4th iteration, while the maximum recognition rates of RCMFE and RCMDE are both 100 %. Nevertheless, they show unstable trends and fluctuations during the iterative process. Although RCMPE can also reach 100 % recognition rate eventually, its convergence speed is far inferior to that of the WMFE method, so WMFE significantly outperforms the other three methods. The combined VMD and WMFE methods are compared with WMFE, to find that the results of fault diagnosis by the former method are relatively satisfactory, probably because VMD decomposes and separates the vibration signal with distinct amplitude and frequency characteristics, thereby essentially enhancing the original signal. The processing results via CNN are visualized in Fig. 16. The VMD and WMFE combined feature extraction method is capable to completely separate the 10 states with the optimal intra class distance and inter class spacing compared with the previous 3 methods.

Fig. 13Curves of different entropies of 4 IMFs: a), b) IMF1, c), d) IMF2, e), f) IMF3, and g), h) IMF4

a)

b)

c)

d)

e)

f)

g)

h)

Fig. 14Visualization feature distribution map using T-SNE

Fig. 15Accuracy of CNN training process by different methods

Fig. 16Visualization feature distribution map using T-SNE

Considering that the performance of the WMFE method may differ by the size of data samples, three sets of 300, 500, and 1000 faulty samples were created for comparison experiments and divided into training sets and testing sets at a ratio of 4:1. According to the results shown in Fig. 17, the WMFE method delivers the highest accuracy among the 300, 500, and 1000 samples. In particular, for small sample sets, WMFE can extract rolling bearing signal features with higher accuracy for taking advantage of the complementary nature of the differences between different entropies to build a more comprehensive representation for the signal information of the high-dimensional feature set. Therefore, the performance and robustness of the method proposed in this paper are superior to the counterparts of RCMDE, RCMFE and RCMPE in extracting signal features.

Fig. 17Accuracy of different methods in different sample sets