Bearing fault diagnosis based on improved VMD and DCNN

2.1. Variation mode decomposition theory

VMD can divide the frequency domain signal adaptively and separate each component effectively. Assuming that VMD decomposes the original input signal $f$ into $K$ mode functions, steps for constructing the variational problem are as follows:

Step 1: Perform Hilbert transformation on sub-signals (modes) $u_k$to calculate the unilateral frequency spectrum of analytic signal as shown in Eq. (1):

where $\delta \left(t\right)$ is the impulse function, $\left\{u_k\left(t\right)\right\}=\left\{u_1\left(t\right),u_2\left(t\right),\dots ,u_k\left(t\right)\right\}$ is mode ensemble, and $K$ is the mode number.

Step 2: $e^{-jwt}$ is multiplied by the analytic signal of each $u_k\left(t\right)$, and modulating the frequency spectrum of each modal component to the corresponding fundamental frequency band, modal function is as follows:

Step 3: By calculating the squared $L^2$ norm of each demodulation signal gradient after translation and estimating the bandwidth of each mode function, the mathematical description of the final constrained variational problem is shown as follows:

where ${\partial }_t$ is the partial derivative of the function, $\left\{{\omega }_k\right\}=\left\{{\omega }_1\right.,{{\omega }_2,{\omega }_3,\dots ,\omega }_k\}$ is center frequencies of $u_k$.

In order to solve the constrained variational problem, the Lagrangian multiplier $\lambda \left(t\right)$ and the quadratic penalty term $\alpha$ are introduced to render the problem into an unconstrained variational problem. The extended Lagrangian expression is shown as follows:

Aiming at the problem Eq. (4), VMD uses the alternate direction method of multipliers (ADMM) to get the optimal solution of the constrained variational model.

At first, the decomposition number of the modal $K$ needs to be set in advance, and ${\widehat{u}}_k^1$, ${\omega }_k^1$, ${\widehat{\lambda }}^1$ should be initialized. Then, according to Eqs. (5) and (6), the mode function $u_k$and their central frequency ${\omega }_k$ are updated alternately:

The Lagrangian multiplier is also updated according to the following Eq. (7):

Repeat the above updating iteration until the convergence condition Eqs. (8) is satisfied:

2.2. Improved variational mode decomposition

In VMD, the decomposition number $K$ needs to be preset artificially. When $K$ is too small, decomposition will cause the loss of some crucial information. When $K$ is too large, the calculation amount of VMD is increased, and frequency aliasing is easy to occur. Therefore, selecting suitable number of decomposition modes is vital in VMD.

This paper proposes a method to determine the number of $K$ based on the change of instantaneous frequency mean of modal component. Instantaneous Frequency (IF) is the derivative of the phase of complex analytic signal, which is an important physical quantity to describe the characteristics of nonstationary signals. IF is based on the Hilbert transform and has instantaneous validity. Assuming that the original signal is decomposed by VMD to obtain $i$ modal components, each of mode has $N$ points, the expression of mean of IF is established as follows:

where $M_{ij}$ denotes the $j$th point of the $i$th mode component, $Q$ is the number of instantaneous frequencies obtained by instfreq function after Hilbert transformation of the $i$th mode component.

For illustrating and verifying the method, the outer ring fault model is used to simulate the periodic impact signal generated by bearing fault, and white noise is added into it. The expression of simulation signal is constructed as follows:

In the above equation, the carrier frequency is $f_n=$3000 Hz, the displacement of constant $x_0=$3, the Damping coefficient is $\xi =$0.05, the period of impact failure is $T=$0.01 s, the sampling frequency is $f_s=$12 kHz, the sampling points is $N=$4096, the sampling time is $t$, and the white noise is $n\left(t\right)$. The simulated time domain waveform is shown in Fig. 1.

Fig. 1Simulated vibration signal of rolling bearing with outer race fault

On the premise that the quadratic penalty term is constant, when $K$ is 1 to 9, the IF mean curves of the simulation signal decomposed by VMD is shown in Fig. 2.

When $K\ge \text{4}$, the curves all appear an obvious downward bending. It can be explained that when $K$ is too large, frequency aliasing will occur in the modal components and intermittent phenomenon will appear. Especially in the high frequency, it causes the mean of IF to drop suddenly, which is the root of the curve bending phenomenon. When $K$ is less than 4, the curve is approximately a straight line, indicating that there is no mode mixing phenomenon in the decomposed components. However, when $K$ is too small, VMD algorithm will filter out some important information in the original signal. Therefore, the number of decomposed modes is determined to be the critical $K=$3 before the curve is bending down.

Fig. 2Instantaneous frequency mean curves of modal components under different K values

In order to automatically judge the bending down of the curve, second derivative is calculated and adopted as a numerical criterion. If signs of second derivative at two adjacent points are opposite and their absolute difference is greater than 0.05, an obvious bending down phenomenon can be found in the curve. It should be noted that 0.05 is a threshold obtained by many trials, which can avoid subtle fluctuations in the curve. To explain this criterion, the instantaneous frequency mean curve of $K=$4 in Fig. 2 is taken as an example and shown in Fig. 3.

Fig. 3The instantaneous frequency mean curves of K= 4

In Fig. 3, A, B, C and D are the points of the instantaneous frequency mean of the modal components $u_1$, $u_2$, $u_3$ and $u_4$. The value of second derivative $f^{\text{'}\text{'}}\left(C\right)$ and the second derivative $f^{\text{'}\text{'}}\left(D\right)$ can be calculated as 0.023 and –0.032. $f^{\text{'}\text{'}}\left(C\right)$ and $f^{\text{'}\text{'}}\left(D\right)$ are opposite signs and their absolute difference is greater than 0.05. Therefore, according to the change of curvature, $K=$4 is the inflection point where the curve has been bending down. So, a suitable decomposition number $K$ can be determined to be 3.

4.1. Processing of original vibration data

Since DCNN fault diagnosis model requires more samples for training, sliding window shown in Fig. 5 is adopted to intercept the original vibration data to construct samples.

Fig. 5Sample construction with sliding window

4.2. Diagnostic process

The flow of rolling element bearing fault diagnosis method based on IVMD and DCNN is shown in Fig. 6, and specific steps are shown as follows:

Steps 1: Process the original vibration data of rolling element bearings by sliding window, and the vibration data samples in different states are obtained.

Steps 2: The number of modes $K$ is determined according to the IF mean curve. Then, use VMD to decompose vibration data into $K$ modal components.

Steps 3: Stack all the decomposed components into a multi-channel sample in order. Then the whole dataset is divided into training set and test set.

Steps 4: Train the designed DCNN model with training dataset. Finally, diagnosis results can be obtained by inputting the test data into the trained DCNN model.

Fig. 6Flow of bearing fault diagnosis method based on IVMD and DCNN

5.1. Description of experimental data

In order to verify the effectiveness and feasibility of the proposed method, the rolling element bearing experimental data provided by Bearing Data Center of Case Western Reserve University (CWRU) [19] is used in this study. The basic layout of the test rig is shown in Fig. 7, it consists of a 2 HP motor, a torque sensor, dynamometer and control electronics. The vibration data with inner race, rolling ball and outer race fault were collected at the drive end bearing. The sampling frequency of the experiment data is 12 kHz, and the rotating speed of shaft is 1727 r/min, namely the rotating frequency $f_r$ is 28.78 Hz.

Bearing state is divided into four types: normal, inner ring fault, outer ring fault and rolling element fault. The defect sizes are 0.007 in, 0.014 in, and 0.021 in respectively, and ten different rolling element bearing states are obtained. The sliding window is adopted to segment the original vibration data. The sliding length is set as 1200, and the moving step length is set as 50. Eventually, there are 1500 samples in each state, and each sample contains 1200 data points.

Fig. 7Test rig of CWRU

5.2. Data preprocessing based on IVMD algorithm

Firstly, IVMD is used to process rolling bearing signals in ten states. For brevity, Fig. 8 shows the decomposition results of four health conditions with fault defects of 0.14 in, including normal condition, inner race fault, outer race fault and rolling element fault conditions. It can be seen that when $K=$4, all the instantaneous frequency curves show obvious downward bending. Similar phenomenon can also be observed in other states. Then, according to the change of curvature, the downward bending point of the bearing signal in ten states is judged to be at $K=$4.

Fig. 8Instantaneous frequency mean curves of modal components under different K

a) Normal

b) IR014

c) OR014

d) B014

In order to further explain the cause of the curve bending down, the signal of inner ring fault is decomposed by the VMD under different number of modes $K$, and the spectrums of original signal and each mode component are shown in the Fig. 9.

Fig. 9The original signal and spectrum of various modal components at different K values

a) The spectrum of original signal

b) Spectrums of mode components when $K=$2

c) The spectrums of mode components when $K=$3

d) The spectrums of mode components when $K=$4

It can be seen that when $K=$2, the frequency band 1000 Hz-1500 Hz are mostly discarded, so it has a lack of information. When $K=$4, the central frequencies of the modal components $u_3$ and $u_4$ are close to each other, it may cause the frequency aliasing. Therefore, the modal decomposition number is determined as $K=$3. And according to the experimental debugging, the penalty parameter $\alpha$ is set to 1000.

After the decomposition of each sample in the experimental dataset, three components containing the fault information are obtained. The original signal of inner ring fault and three modal components obtained by IVMD method are given in Fig. 10.

Then, normalize all the decomposition data, and label them with the one-hot structure. The whole dataset is divided into training dataset and test dataset, which are shown in Table 1.

Fig. 10VMD results of bearing vibration signal with inner race fault

5.3. Design of DCNN model

The structural parameters of the designed DCNN model after repeated experiments and debugging are shown in Table 2. CS represents the size of convolution kernel. CN represents the depth of convolution kernel. I denotes the number of input graphs. Strides represents the step length of convolution kernel movement. S represents the width of pooling layer. FS₁ and FS₂ represent the size of the first and second full connection layers, respectively.

RMSProp algorithm is adopted to optimize the DCNN model, and the global learning rate $\lambda$ is set as 0.001. To prevent overfitting, dropout strategy and L2 regularization algorithms are adopted, with coefficients of 0.5 and 0.0012. In this way, the generalization performance of DCNN model is highly improved.

5.4. Analysis of experimental results

A total of 7500 samples in the training set are input into the designed DCNN model. Fig. 11 displays the change of the loss function values of the training set. When the model training reaches 2100 times, the loss function no longer changes significantly, which means that the model training is completed at this time. After multiple iterations, the training accuracy of the model is 100 % and the average diagnostic accuracy in the test set is 99.5 %.

In order to present the diagnostic results of the model in the test set in more detail, a confusion matrix is introduced to analyze the experimental results. In Fig. 12, B007, B014 and B014 represent the rolling element fault with the fault size of 0.07 in, 0.14 in and 0.21 in, respectively.

Except that the diagnostic accuracy of rolling element fault with the fault size of 0.07 in is 96.7 %, identification accuracies of other 9 states is higher than 98.8 %, which indicates that this method can effectively realize the intelligent bearing fault diagnosis and obtain a high diagnostic accuracy.

Meanwhile, in order to compare the advantages of the proposed method in fault diagnosis, the above experimental data is also decomposed by EEMD method, where the intensity of white noise and the number of cycles are respectively set to Nstd = 0.3 and NE = 100. Then, the first three modal components with large kurtosis constitute a data sample set. In the light of Table 2, the data sets processed by EEMD are divided in the same way. Using the training dataset to train the DCNN model. Finally, the diagnosis accuracy of test dataset is 95.07 %.

From the above comparison, it can be seen that the proposed method based on IVMD and DCNN has a higher diagnostic accuracy, indicating that the IVMD algorithm can enhance the fault characteristics more effectively.

Fig. 11Loss function curve for training set

Fig. 12Diagnostic result of test dataset

5.5. t-SNE visual analysis

In order to more intuitively illustrate the adaptive feature learning ability of the DCNN model designed in this paper, t-SNE algorithm [20] is adopted to visualize the features learned from the first layer of the full connection layer of the DCNN model. The test samples are input into the trained DCNN model, and the feature distribution is shown in Fig. 13.

Each color in Fig. 13 represents a fault type, and the features in each state learned by the model are highly separated. The visualization results further show that the DCNN model can adaptively learn different state characteristics from the bearing vibration data after IVMD processing, and thus has a good fault classification capability.

Fig. 15Feature visualization

5.6. Comparison with other deep learning methods

In order to validate the effectiveness and superiority of the proposed method, the proposed method based on DCNN is further compared with other deep learning-based fault diagnosis methods, including deep belief network (DBN), sparse autoencoder (SAE) and one-dimensional CNN (1-D CNN) [21]. In every method, the same dataset is adopted. The structure of the DBN is 27-100-50-10. Its learning rate is 0.1, and the number of iterations is 400. SAE contains three hidden layers, each of which has a structure of 100-60-10. The learning rate is 0.2, and the number of iterations is 200. Structural parameters of 1-D CNN are just similar to those in literature [21]. In all these comparative models, learned signal characteristics are classified by softmax classifier. In order to eliminate the influences of randomness, each algorithm is repeated 20 times. Fault diagnosis results of these methods are summarized in Table 3.

Compared with DBN and SAE models, the proposed method using DCNN model has the highest diagnostic accuracy and lowest standard deviation. Meanwhile, compared with the proposed DCNN-based method, the 1-D CNN model applied in literature [21] only has two convolution layers. Although its recognition rate reaches 99.15 %, it’s still lower than the proposed method. Moreover, the standard deviation of diagnosis results is much higher than the proposed DCNN model. These results proves that the DCNN model with deeper structure can improve the accuracy and stability of diagnosis results.