Fault diagnosis of rolling bearing based on hierarchical discrete entropy and semi-supervised local Fisher discriminant analysis

2.1. Principle of variational mode decomposition method

VMD is a new adaptive time-frequency decomposition method. It integrates the concepts of Wiener filtering and Hilbert transform, outlier demodulation, frequency mixing and signal analysis. Through a series of iterations, the original signal can be decomposed into multiple IMF with different center frequencies and bandwidths. The constraint is that the sum of the decomposed modal components is equal to the input signal and the sum of the bandwidths of the decomposed modal components should be minimized, so the constrained model is shown in Eq. (1):

where, ${\partial }_t$ is the partial derivative of $t$, $\delta \left(t\right)$ is the impulse function, ${\mu }_k$ is the $k$th modal function, ${\omega }_k$ is the center frequency of each mode, and $f$ is the original signal.

On this basis, the quadratic penalty factor $\alpha$ and the Lagrange multiplication operator $\lambda \left(t\right)$ are introduced to transform the above variational model into an unconstrained variational model. The transformed Lagrange expression is shown in Eq. (2):

where, the $\alpha$ parameters are used to ensure the accuracy of the reconstructed signal, which $\lambda \left(t\right)$ can make the constraints more stringent.

The Lagrange expression can be solved using the alternating power multiplier algorithm, and the saddle point of the Lagrange expression can be solved by alternately updating the ${\mu }_k^{n+1}$, ${\omega }_k^{n+1}$ and ${\lambda }^{n+1}$. Among them, it can be expressed by Eq. (3):

The specific process of the VMD algorithm is as follows:

Step 1: Initialize ${\mu }_k^1$, ${\omega }_k^1$, ${\lambda }_k^1$ and $n$;

Step 2: Update ${\mu }_k$:

Step 3: Update ${\omega }_k$:

Step 4: Update $\lambda$:

Step 5: If:

($\epsilon >0$, this paper $\epsilon =$1×10^-7), then stop the iteration, otherwise return to the second step.

2.2. Principle of hierarchical discrete entropy feature extraction method

Based on the advantages of hierarchical entropy and the definition of discrete entropy, the calculation process of hierarchical discrete entropy (HDE) is as follows [17]:

1) Given a time series $\left\{u\left(i\right),i=\mathrm{1,2},...N\right\}$, the hierarchical operators $Q_0$ and $Q_1$ are defined as:

Among them, $N=2^n$ is a positive integer, and the length of operators $Q_0$ and $Q_1$ is $2^{n-1}$.

According to operators $Q_0$ and $Q_1$, the original sequence can be reconstructed as:

When $j=0$ or $j=1$, the matrix operator $Q_j$ can be defined as:

2) Construct a dimension vector $\left[{\gamma }_1,{\gamma }_2,...,{\gamma }_n\right]\in \left\{\mathrm{0,1}\right\}$ that corresponds to a positive integer $e$.

3) For vector $\left[{\gamma }_1,{\gamma }_2,...,{\gamma }_n\right]$, the node components $u\left(i\right)$ decomposed at each level are defined as:

where $u_{k,0}$ and $u_{k,1}$ are the low-frequency and high-frequency parts of the original time series under the decomposition layers of layer $k$.

4) The hierarchical sequence $u_{k,e}$ is mapped to $\left[y_1,y_2,...,y_{2^k}\right]$ through a normal cumulative distribution function $y_i=\frac{1}{\sigma \sqrt{2\pi }}{\int }_{-\infty }^{x_i}{e}^{\frac{-{\left(t-u\right)}^2}{2{\sigma }^2}}dt$, where $y_i\in \left(\mathrm{0,1}\right)$.

5) $y_i$ is assigned to integer $1~c$ through function $z_i^c=round\left(c\mathrm{*}{y}_i+0.5\right)$. $c$ is the number of categories. $z_i^c$ is reconstructed as $z_j^{m,c}$ by embedding dimension $m$ and delay parameter $d$:

6) Construct all possible discrete models ${\pi }_{v_1...v_{c^m}}$ by embedding dimension $m$ and number of categories $c$. $z_j^{m,c}$ matches the discrete model ${\pi }_{v_1...v_{c^m}}$ one by one, and Eq. (14) is used to calculate the frequency of each discrete model ${\pi }_{v_1...v_{c^m}}$ in the reconstruction sequence $z_j^{m,c}$:

7) A single discrete entropy can be expressed as:

8) Finally, $HDE$ can be expressed as:

2.3. Principle of semi-supervised local fisher discriminant analysis feature dimensionality reduction method

Assuming a given sample set contains $D$-dimensional features, $C$ categories, denoted as $X=\left\{x_i\in R^D\left(i=\mathrm{1,2},...,n^{\mathrm{\text{'}}},...,n\right)\right\}$, with labeled sample $x_i\left(i=\mathrm{1,2},\dots ,n^{\mathrm{\text{'}}}\right)$ and category labels denoted as $l_i\in \left\{\mathrm{1,2},...,C\right\}\left(i=\mathrm{1,2},\dots ,n^{\mathrm{\text{'}}}\right)$, the global divergence matrix of PCA is defined as [15]:

The weight $W_{i,j}^{\left(t\right)}=1/n$, then the optimization objective function of PCA is:

The local inter class divergence matrix $S^{\left(lb\right)}$ and local intra class divergence matrix $S^{lw}$ of LFDA can be defined as:

Among them, the weight matrices $W^{\left(lb\right)}$ and $W^{\left(lw\right)}$ are defined as:

Among them, $n_{l_i}^{\mathrm{\text{'}}}$ represents the number of Class $l_i\in \left\{\mathrm{1,2},...,C\right\}\left(i=\mathrm{1,2},...,n^{\mathrm{\text{'}}}\right)$ samples; The $\left(i,j\right)$th element $A_{i,j}\in \left[\mathrm{0,1}\right]$ of the similarity matrix $A$ is used to describe the similarity between two samples $x_i$ and $x_j$, defined as $A_{i,j}=\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{{‖x_i-x_j‖}^2}{{\sigma }_i{\sigma }_j}\right)$, where ${\sigma }_i$ is the local scale of sample point $x_i$, defined as ${\sigma }_i=‖x_i-x_i^{\left(k\right)}‖$, and $x_i^{\left(k\right)}$ is the $k$th nearest neighbor of $x_i$.

The inter class divergence matrix and intra class divergence matrix of SELF can be defined by Eqs. (17) to (22) as:

Among them, the weight coefficients $\beta \in \left[\mathrm{0,1}\right]$ and $I_d$ are the standard matrices. When $\beta =0$, SLEF is equivalent to LFDA, and when $\beta =1$, SELF is equivalent to PCA. Find the optimal projection transformation matrix $T$, that is, solve the problem of maximizing the objective function as follows:

The solution of the transformation matrix in the above equation is equivalent to the problem of obtaining the generalized eigenvectors of $S^{\left(b\right)}\alpha =\lambda S^{\left(w\right)}\alpha$. The transformation matrix $T$ is composed of the generalized eigenvectors $\left({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_d\right)$ corresponding to the first $d$ maximum generalized eigenvalues.

4.1. Case Western Reserve University (CWRU) dataset

The CWRU bearing dataset was collected from the experimental platform in Fig. 1 [18]. The collected bearing vibration signal was conducted at a motor speed of 1797 r/min and a sampling frequency of 12000 Hz. The fault sample at the fan end was selected for analysis. Firstly, establish fault category labels based on different fault sizes and locations, as shown in Table 1. Taking the bearing signal with a fault size of 0.1778 mm as an example, the waveform is shown in Fig. 2.

Fig. 1Case Western Reserve University bearing test bench

Fig. 2Fault waveform diagram

Divide ten different types of bearing fault signals into 50 sets of test set samples and 50 sets of training set samples. Due to the complexity of the vibration signal, the sample is preprocessed using VMD method, and several intrinsic mode components are obtained after decomposition. Taking the rolling element fault signal as an example, different pre-decomposition modes $k$ are selected for VMD, in which the penalty factor a takes the default value 2000. The main difference of different modes lies in the difference of the central frequency. According to the VMD algorithm, the central frequency of each IMF component obtained from the VMD of the vibration signal will be distributed from low frequency to high frequency. If you want to obtain the optimal preset scale $k$, the center frequency of the last order IMF component should be the maximum for the first time. Table 2 shows the central frequencies of each IMF component under different modes.

As can be seen from Table 2, the minimum value of the central frequency is taken from the initial preset scale $k=2$, and the maximum value of the central frequency tends to be stable when $k=4$ and $k=7$. In order to prevent the mode number$k$of VMD from being too large and over-decomposing, the$k$value selected in this paper is 4.

Taking the rolling element fault signal with a fault size of 0.1778mm as an example, the decomposed time-domain and frequency spectrum are shown in Fig. 3. From the time-domain graph, it can be seen that the VMD algorithm can decompose the original signal into multiple intrinsic mode components with different amplitudes through its built-in filtering algorithm. Observing the spectrum graph, it can be seen that the spectra of each component have obvious peaks, and there is no situation where multiple peaks appear simultaneously, avoiding mode mixing and boundary effects. From this, it can be seen that the VMD algorithm can effectively preprocess signals and improve the accuracy of subsequent work.

Fig. 3Time domain and frequency spectrum after VMD

a) Time domain diagram after VMD

b) Spectrum after VMD

Use the HDE method to extract features from the decomposed multiple intrinsic mode components. The parameters in HDE affect the size of the extracted entropy and the overall effect. In order to determine the sensitivity of HDE to the original signal under different parameters, the average and standard deviation at different levels are calculated, and the coefficient of variation is introduced to determine the degree of node dispersion. The coefficient of variation is equal to the standard deviation/mean. As defined, the smaller the coefficient of variation, the higher the stability of the signal [18]. Fig. 4 shows the size of HDE entropy values under different parameters, while Tables 3 to 6 provide the coefficient of variation of HDE entropy values under different parameters.

Due to the characteristics of the HDE algorithm, the final extracted sample length is $2^n$. If the length parameter $n$ is set too small, the extracted entropy value is too small to be meaningful for discussion. If it is set too large, it will lead to excessive extraction, resulting in an invalid entropy value. Therefore, when discussing the sample length parameter $n$, setting the range to [4, 7] will result in extracted entropy sample lengths of 16, 32, 64, and 128. From Fig. 4(a), it can be seen that when the sample length is 16 and 32, the entropy values are distributed within the [1.2, 1.4] range, while larger sample lengths only cause the entropy value to fluctuate and increase. From Table 3, it can be seen that the coefficient of variation is the smallest when the sample length parameter is 5, which also represents that the signal is most stable when the extracted sample length is 32. From Fig. 4(b), it can be seen that different lengths of original samples also affect the size and stability of the final extracted entropy. When the length of the original sample is $N\ge \text{512}$, HDE tends to stabilize. Combined with Table 4, it can be seen that when the length of the original sample is 1024, the coefficient of variation is the smallest and the entropy value is also the most stable.

The different embedding dimensions $m$ and category $c$ also affect the performance of HDE. From Fig. 4(c), it can be seen that as $m$ increases, the fluctuation of entropy value increases. From Table 5, it can also be seen that when $m=2$, the coefficient of variation is the smallest and the entropy value is the most stable. Therefore, the optimal embedding dimension $m=2$ is chosen; Similarly, by observing Fig. 4(d) and Table 6, the optimal number of categories $c=2$ can be obtained.

Fig. 4Entropy values under different conditions

a) Different entropy extraction sample lengths$\mathrm{n}$

b) Different original sample lengths $N$

c) Different embedding dimensions$\mathrm{m}$

d) Different categories$\mathrm{c}$

In order to improve the accuracy of extraction, first extract the entropy information sample from the signal that has not been decomposed by VMD, and then use this sample as a benchmark to perform secondary extraction on the decomposed intrinsic mode components. The HDE parameters selected in this article are: original sample length $N=1024$, embedding dimension $m=2$, number of categories $c=2$, and sample length parameter $n=5$. Due to secondary extraction, the length of each extracted sample increased from 32 to 64. The feature sample set extracted under this parameter is 100×64. In order to compare the effectiveness of the HDE method, the hierarchical range entropy (HRE) and hierarchical sample entropy (HSE) under the same parameters were also used for feature sample extraction. HDE is similar to HRE and HSE, which integrates range entropy and sample entropy on the basis of hierarchical entropy. The difference lies in the different calculation methods for signal sequences, resulting in different sets of extracted feature data. Taking the rolling element fault signal with a fault size of 0.1778 mm as an example, Fig. 5 shows a comparison of three entropy values.

The entropy value can reflect the complexity of the bearing vibration signal, and the larger the entropy value, the higher the complexity of the vibration signal it represents. From the figure, it can be seen that after taking the mean of the extracted 100×64 feature samples, although the samples extracted by the HRE method have stability, their entropy values are between [0, 0.5], which cannot reflect the characteristics of the original fault signal well; The entropy value extracted by HSE method has increased compared to HRE, but its entropy value exhibits significant volatility, which is not conducive to the final fault classification; The entropy value extracted by the HDE method has high stationarity, and the entropy value is also larger than the other two extraction methods, indicating that this method can better demonstrate the complexity of signal features, which is beneficial for subsequent signal processing. To further demonstrate the effectiveness of the HDE method, Pearson correlation coefficient analysis was performed on the feature samples extracted by the three entropy extraction methods under four fault states with a fault size of 0.1778 mm. The results are shown in Fig. 6.

Fig. 5Comparison chart of three entropy values

The Pearson correlation coefficient measures the correlation between the two, with a value range of [–1, 1]. After taking the absolute value of the correlation coefficient, the larger the value, the more relevant the two are. Comparing Fig. 6(a) and Fig. 6(b), it can be seen that the features extracted by HDE and HRE have a high correlation with both inner-outer faults, reaching above 0.6. However, for other types of faults, such as normal-ball faults, inner-ball faults, HDE Pearson correlation coefficients are much smaller than those of the HRE method, with only 0.02998 and 0.01139, while the HRE method only shows a small correlation when comparing ball-outer faults, and the correlation between other faults is not much different from that of the HDE method; Comparing three figures, it can be seen that the correlation between HSE and the extracted samples of inner-outer faults reaches 0.8227, which is much higher than the other two types of entropy values, indicating that the entropy values extracted by HSE for these two types of faults are highly correlated, which is not conducive to the final fault classification. The samples extracted by HSE showed good non correlation between normal-outer faults, as well as between outer-ball faults, but the non correlation between other faults was far inferior to the other two entropy values. In summary, the entropy extracted by HDE has better non correlation between each fault type compared to the other two methods, which can make the final fault diagnosis more accurate.

After feature extraction of the fault signal, the semi-supervised local fisher discriminant analysis (SELF) method was used to reduce its dimensionality. The weight coefficient $\beta =0.5$ of the SELF method, the reduced vector length $d=32$, and the vector length obtained after SELF processing were reduced from 100×64 to 100×32. For comparative analysis, the entropy data extracted by HRE was also subjected to dimensionality reduction using LLE, KPCA combined with LDA, and PCA combined with LDA. Taking the feature samples under four fault states with a fault size of 0.1778 mm as an example, the t-SNE scatter plots of the reduced samples are shown in Fig. 7.

Although the four types of fault samples in Fig. 7(a) have a certain distribution pattern after dimension reduction using the LLE method, their various fault sample points are relatively scattered, and there is no clear boundary between different types of faults. The samples in Fig. 7(b), after dimensionality reduction using the KPCA-LDA method, although the sample points of various types of faults have a certain degree of aggregation compared to the LLE method, there is a mixture of different lengths in the four types of fault samples, and there is a significant separation phenomenon between the outer ring fault and rolling element fault samples, which greatly affects the accuracy of the final fault classification. From Fig. 7(c), it can be seen that the data after dimensionality reduction using PCA-LDA method has a clear boundary for four different types of faults, and there is no mixing of feature samples for each type; The four fault categories in Fig. 7(d) are similar to those in Fig. 7(c), with clear boundaries. However, the clustering degree of the data samples in Fig. 7(d) is closer than that in Fig. 7(c), indicating that the data samples after SELF dimensionality reduction can effectively distinguish different types of fault sample types, and the numerical values of different types of samples are divided into a certain range, resulting in higher discrimination between different types of samples and better fault diagnosis classification.

Fig. 6Pearson correlation coefficient chart

a) HDE Pearson correlation coefficient

b) HRE Pearson correlation coefficient

c) HSE Pearson correlation coefficient

Input the training and testing sets with a length of 100×32 obtained after dimensionality reduction into the least squares support vector machine optimized by particle swarm optimization for training and recognition. The optimization range of the regularization parameters and kernel parameters in the least squares support vector machine is set to [1, 200], the number of iterations is 10, and the maximum and minimum velocities of the particle swarm are 10 and –10. The least squares support vector machine performs 50 internal cross validations during fault classification of samples. The optimal parameter combination and accuracy obtained after training and testing the classification model are shown in Table 7. In order to demonstrate the effectiveness of this method, three entropy extraction methods and three data dimensionality reduction methods were combined to form a total of 12 methods for fault diagnosis. Table 7 also shows the classification accuracy of feature sample sets extracted by other methods in the model.

From the classification accuracy of LSSVM, we can see that for the three different entropy extraction methods, the correlation between the feature samples of different fault types is smaller under the HDE entropy extraction method, so it is easier to distinguish different fault types. From the final accuracy, we can see that after different dimensionality reduction methods, HDE method has higher accuracy than the other two entropy extraction methods; for different dimensionality reduction methods, the accuracy of SELF method is higher than that of the other three methods, which also proves the superiority of the semi-supervised local Fisher discrimination method proposed in this paper.

In order to further demonstrate the effectiveness of the method proposed in this paper, the deep learning classifier of deep extreme learning machine (DELM) is selected to compare again. This method effectively combines autoencoder with extreme learning machine to form a deep neural network structure. It combines the advantages of autoencoder and extreme learning machine, and makes use of the feature learning ability of autoencoder and the fast training characteristics of extreme learning machine, it can better extract the advanced features of data in depth structure, so as to improve the modeling ability of complex data, and can overcome some shortcomings of traditional deep neural network training, such as long-time training, gradient disappearance and so on. Thus, the training efficiency and the generalization ability of the model are improved.

Fig. 7t-SNE visualization diagram

a) t-SNE graph after LLE

b) t-SNE graph after KPCA-LDA

c) t-SNE graph after PCA-LDA

d) t-SNE graph after SELF

Similarly, the reduced samples are input into the depth extreme learning machine optimized by the optimization algorithm for training and recognition. The optimization range of weight parameters of depth extreme learning machine is set to [–10, 10], the number of iterations is 50, and the number of iterative populations is 20. The deep extreme learning machine performs 50 internal cross-validation in the process of sample fault classification. The diagnostic accuracy of 12 fault diagnosis methods in the optimized DELM model is shown in Table 8. As can be seen from Table 8, the accuracy of three different entropy extraction methods after dimensionality reduction by SELF method has reached 100 %, indicating the advantage of SELF method in data dimensionality reduction. For the same dimensionality reduction method, the accuracy of HDE entropy extraction method is also higher than that of other entropy extraction methods, which shows the accuracy of the proposed method. Overall, the optimized DELM method has higher classification accuracy than LSSVM method, and can better classify samples.

4.2. Southeast University (SEU) dataset

Due to the fact that the Case Western Reserve University bearing dataset achieved high accuracy in extracting feature vectors using three entropy extraction methods and then reducing the dimensionality through SELF, both HRE-SELF and HDE-SELF achieved 100 % accuracy. Therefore, in order to better demonstrate the accuracy of the method proposed in this article, further experiments will be conducted on the dataset of Southeast University. The Southeast University (SEU) dataset is a gearbox dataset provided by Southeast University [19], and the fault test bench is shown in Fig. 8. This dataset consists of two sub datasets, including the bearing dataset and the gear dataset. The bearing dataset simulates five bearing operating states under two different operating conditions, namely 20 Hz (1200 r/min) - unloaded 0 V (0 N/m) and 30 Hz (1800 r/min) - loaded 2 V (7.32 N/m). The specific fault types are shown in Table 9, and the fault wave-forms are shown in Fig. 9.

Fig. 8Southeast University bearing test bench

Fig. 9Fault waveform diagram

a) 20 Hz-0 V fault signal

b) 30 Hz-2 V fault signal

Similar to the previous text, ten different types of fault signals are divided into training and testing sets, preprocessed through VMD decomposition, and their HDE entropy values are extracted. The HDE parameter settings are the same as the previous text, but will not be further elaborated here. Further perform dimensionality reduction on the extracted 100×64 entropy samples. Taking 5 types of fault samples with a speed load condition of 20 Hz-0 V as an example, Fig. 10 show the t-SNE scatter plots of four dimensionality reduction methods under HDE method.

From Fig. 10(a) and 10(b), it can be seen that the feature samples after dimensionality reduction using LLE and KPCA-LDA methods still exhibit hybridization, and various types of fault samples are relatively scattered without good discrimination; The PCA-LDA method in Fig. 10(c) performs much better in reducing the dimensionality of the samples compared to the previous two methods, but there is still a situation of mixed sample points in the dimensionality reduction of the first three types of fault samples; From Fig. 10(d), it can be seen that the sample points of each fault type in the data after dimensionality reduction using the SELF method are mapped to separate regions, with good discrimination.

Input the dimensionality reduced data samples into the least squares support vector machine optimized by particle swarm optimization for training and testing. Table 10 shows the optimal parameter combinations and classification accuracy of 12 methods. From the previous analysis and the data in the table, it can be seen that the HDE entropy extraction method has a smaller correlation between the extracted samples compared to the other two entropy extraction methods, which makes its accuracy higher than the other two entropy values after the same method of dimensionality reduction, verifying the previous analysis; From the previous analysis, it can be seen that among the four dimensionality reduction methods, the discrimination between datasets under SELF dimensionality reduction is higher, which is also reflected by the final classification accuracy. The accuracy of the HDE combined with SELF dimensionality reduction method reached 98.2 %, verifying the effectiveness of the HDE-SELF method proposed in this paper.

Fig. 10t-SNE visualization diagram

a) t-SNE graph after LLE

b) t-SNE graph after KPCA-LDA

c) t-SNE graph after PCA-LDA

d) t-SNE graph after SELF

As mentioned above, the DELM method is used to classify SEU fault samples. As can be seen from Table 11, the final fault diagnosis accuracy of the samples processed by HDE and SELF is the highest, reaching 99.6 %. The other two entropy extraction methods also achieve high accuracy in the case of SELF processing, which once again proves the effectiveness of the SELF method. For different entropy extraction methods, the accuracy of HDE method is also higher than that of the other two methods. Generally speaking, the data processed by the HDE-SELF method proposed in this paper has high classification accuracy, which also shows the effectiveness of the method.