# Fault Diagnosis of Spindle Device in Hoist Using Variational Mode Decomposition and Statistical Features

Authors: Jun Gu; Yuxing Peng; Hao Lu; Shuang Cao; Bobo Cao
Journal: Shock and Vibration
 (2020)

---

## Abstract

By analyzing nonlinear and nonstationary vibration signals from the spindle device of the mine hoist, it is a challenge to overcome the difficulty of fault feature extraction and accurately identify the fault of rotor-bearing system. In response to this problem, this paper proposes a new approach based on variational mode decomposition (VMD), SVM, and statistical characteristics such as variance contribution rate (VCR), energy entropy (EE), and permutation entropy (PE). Comparisons have gone to evaluate the performance of rolling bearing defect by using EMD (Empirical Mode Decomposition), MEEMD (Modified Ensemble EMD), BP (Back Propagation) network, single or multiple statistical characteristics, and different motor loads. The experiment was carried out on the mechanical failure simulator of the main shaft device of the hoist, which verified the reliability and effectiveness of the method. The results show that the diagnosis method is suitable for feature extraction of bearing fault signals, with the highest diagnosis accuracy. It can provide a good practical reference for the fault diagnosis of mechanical equipment of the hoist spindle device and has certain practical value.


---

## Body

## 1. Introduction

As the “throat” of the mine, the hoist is an important mechanical device connecting the underground and the ground. It is used to drive minerals, people, and important equipment by winding wire ropes. The spindle device is the key component of mine hoist which is responsible for winding steel wire rope, bearing various impact loads, and transferring work. Once the spindle device has a fault defect, it will affect the performance of the whole lifting system to a certain extent and also bring potential safety hazards to the mine production. Therefore, according to the special working condition of mine production, it is very necessary to study the fault diagnosis method of shaft device of mine hoist. The spindle device of the hoist belongs to a typical rotor-bearing rotating system. There are many types of faults in this system, mainly including rotor imbalance, misalignment, and rolling bearing faults [1]. According to incomplete statistics, 70% of the failures of rotating machinery are caused by vibration [2], among which 30% are caused by bearings, 90% by inner and outer rings, and only 10% by rolling body and cage.The failure of rolling bearing mainly occurs in the form of normal state, outer ring damage, inner ring damage, rolling body damage, and comprehensive situation. During the working process of the bearing, it will constantly contact with other parts, producing nonstationary and nonlinear vibration signals. Vibration signals usually contain a large number of mechanical health-related data [3]. Because it is easily affected by the noise in the surrounding environment, the fault characteristic signal is submerged in the vibration signal. How to extract fault characteristic signals effectively is an important step in bearing fault diagnosis. So far, the processing tools for vibration signal analysis used in mechanical fault diagnosis can be divided into three categories: time-domain analysis, frequency-domain analysis, and time-frequency analysis [4]. Time-domain analysis is the simplest analysis method, which mainly analyzes the kurtosis, root mean square, average value, waveform factor, and other statistical indicators extracted directly from the vibration signal time series to diagnose the fault [5]. The frequency-domain method is usually used to diagnose mechanical faults by observing the hidden fault characteristic frequency in the vibration signal [6]. Time-domain analysis and frequency-domain analysis can only extract some useful information from vibration signal because of the loss of information in other fields. In actual situations, the spindle device will generate nonstationary vibration signals. Time-domain and frequency-domain techniques are not suitable for nonstationary or time-varying signals because they are amplitude and frequency-modulated signals [7].Therefore, the time-frequency analysis technology of nonlinear and nonstationary data decomposition has become the focus of fault diagnosis research in recent years. Traditional methods for dealing with nonlinear and nonstationary sequences, such as short-time Fourier transform (STFT) [8], EMD [9], EEMD [10], wavelet transform (WT) [11], CEEMD [12], and MEEMD [13], have been used to study nonstationary signals. These techniques have been used for bearing fault diagnosis.The traditional STFT has the problem of selecting the width and type of window function. WT has the problem of selecting the mother wavelet and the number of layers. The wavelet method cannot decompose the signal adaptively in the time-frequency domain, and its disadvantage is pointed out in reference [14]. Hilbert-Huang transform combines EMD with Hilbert transform to conduct signal decomposition analysis [15]. The EMD method has been used by many researchers for rolling bearing and other failure analyses [16–18]. However, the traditional EMD method has some problems such as pattern mixing, endpoint effect, and complex calculation. In the past decade, researchers have adopted improved methods to remedy some of the problems in EMD, and proposed methods such as EEMD, CEEMD, and MEEMD, which have shown good results. But it still needs to be revised to meet the strict definition of the IMF.Recently, Dragomiretskiy and Zosso [19] introduced an adaptive signal decomposition algorithm called variational mode decomposition (VMD). By searching the optimal solution of variational modulus iteratively, VMD can find the frequency center and bandwidth of each component adaptively. In the past few years, researchers have found that VMD provides better results than other technologies. Aiming at the rolling bearing faults of multistage centrifugal pumps, Zhang et al. [20] established and studied the fault signal calculation models of rolling bearing defects at different positions based on the methods of VMD, FFT, and envelope analysis. Experiments show that VMD can accurately extract the main mode of bearing fault signals and is better than EMD in terms of bearing defect feature extraction and the method can successfully diagnose rolling bearing failures. Ding et al. [21] proposed a gear fault diagnosis method based on genetic variation particle swarm optimization VMD and probabilistic neural network algorithm for decomposition number K and penalty factor. The superiority of the proposed method in gear fault diagnosis is verified by comparing with other methods. These applications demonstrate that the VMD method can effectively deal with complex signals. However, at present, the application of VMD to mechanical fault vibration signal is still rare.In the past two decades, many intelligent technology methods, such as artificial neural network (ANN) [18], generative adversarial network (GAN) [22], support vector machine (SVM) [23], K-nearest neighbor (KNN) [24], and deep neural network (DNN) [25], have been developed. Researchers often use classifiers combined with signal decomposition methods to improve the accuracy of fault diagnosis. Based on the characteristics of the sample data and the application of SVM in fault diagnosis of rotating machinery, SVM is selected as the fault classifier in this paper.The main contributions of this work are as follows:(1) a novel method of intelligent fault diagnosis based on the combination of WPT as a powerful tool for removing the noise influence on original vibration signals, VMD as a decomposition algorithm in order to extract features and SVM as a fault classifier is proposed; (2) up to now, the VMD has been applied in fault diagnosis widely. But the VMD is rarely conducted on fault diagnosis for the spindle device of hoist; (3) the single and combined feature vectors are used to train and evaluate the SVM fault diagnosis model. The effectiveness of SVM is discussed and compared with other classifiers; (4) the classifier can be trained effectively and real-time bearing state detection can be realized with limited training data set; (5) the proposed method in the paper is applied to rolling bearing fault diagnosis under different motor loads and the results show that the proposed method performs good robustness.The rest of the paper is organized as follows. In the second section, the theoretical foundations are presented. The flow of the rolling bearing fault diagnosis method is introduced in the third section. Experimental analysis is carried out in the fourth section. The experimental results are discussed in the fifth section. The conclusion is in the sixth section.


## 2. Theoretical Foundation of Signal Processing and Feature Extraction

### 2.1. Wavelet Packet-Based Denoising

In practical engineering applications, the signals generated by mechanical equipment are nonstationary signals composed of useful signals and noise signals. When performing noise processing, signal decomposition, threshold quantization processing, and signal reconstruction are required. In signal decomposition, wavelet transform and wavelet packet transform have different decomposition characteristics.


#### 2.1.1. Wavelet Decomposition

Wavelet analysis is a time-frequency localized analysis method developed on the basis of Fourier transform, with a fixed window size, but its shape can be changed. It is generally used to process nonstationary signals. When the wavelet transform decomposes the signal, only the low-frequency part is further decomposed, and the high-frequency part, that is, the detailed part of the signal, is no longer decomposed. Therefore, the wavelet transform can well characterize the signal with the low-frequency signal as the main component but cannot decompose and represent signals containing a large amount of detailed information, such as nonstationary mechanical vibration signals, biomedical signals, and seismic signals. The schematic diagram of wavelet decomposition is shown in Figure1, in which three layers of wavelet decomposition are done on the original signal S, A represents the approximate value of the signal, and D represents the detailed value of the signal.Figure 1
The 3-layer wavelet decomposition.The choice of threshold rules (Minimax, Sqtwlolg, Heursure, Rigrsure), decomposed layers (generally 3 to 5 layers), and wavelet basis functions (Haar, Mexican Hat (mexh), Morlets, Gaussian, Daubechies (dbN), Meyer, Symlets (sym), Coiflet (coifN), etc.) will have a certain effect on the noise reduction effect.When the signal continuity and smoothness is good, it is appropriate to use the Sym8 wavelet, mainly because the wavelet basis of Sym8 has good smoothness and continuity. When the signal continuity is weak, the Haar wavelet basis should be used. The db wavelet base is more suitable for the unsteady signals such as structural vibration and microseismic signals.


#### 2.1.2. Wavelet Packet Decomposition

Wavelet packet decomposition is an improvement on wavelet decomposition and is more refined than wavelet decomposition. The wavelet packet transform divides the signal more finely. It can decompose the low-frequency part and the high-frequency part of the signal and can adaptively select the corresponding frequency band according to the characteristics of the analyzed signal to make it match the new frequency spectrum, thereby improving the time-frequency resolution rate. At present, wavelet packet analysis has been widely used in data denoising, vibration information extraction, and other monitoring fields [26, 27].The following are the specific steps of wavelet packet noise reduction: (1) the signal is decomposed byN-layer wavelet packet. The schematic diagram of wavelet packet decomposition is shown in Figure 2, in which three layers of wavelet packet decomposition are done on the original signal S; (2) the minimum cost principle is used to calculate the optimal wavelet packet base; (3) threshold quantization of wavelet packet decomposition coefficient. For each wavelet packet decomposition coefficient, select an appropriate threshold value for threshold quantization of left and right wavelet packet bases; (4) the reconstructed signal is the signal after the optimal wavelet packet noise reduction.Figure 2
The wavelet packet 3-layer decomposition.Based on the characteristics of wavelet packet, wavelet packet decomposition is used to denoise the collected vibration signals.


### 2.2. Variational Mode Decomposition

The essence of the VMD method is an all-modal variant that adaptively decomposes the original signalfx into a modal function uk with a central frequency of wk for K, where K is the predetermined number of components. Assuming that fx is a signal with a multi-frequency component, the VMD algorithm redefines IMF as a signal that is both amplitude- and frequency-modulated, with the expression(1)ukt=Aktcosϕkt,where ϕkt is the phase and Akt is the envelope magnitude of the signal. The constraint variance model can be expressed as(2)minuk,wk∑k=1K∂tδt+jπt∗ukt∗e−jwkt22,where ukt satisfies ∑k=1Kuk=fx, δt is the Dirac delta function, and “∗” is the convolutional sign; uk is the set of components for IMFs; wk is the set of the corresponding central frequency bands for K components; and ∂t denotes the bias.In order to make the constrained-variance problem unconstrained-variance problem, the secondary penalty factorα, which guarantees the accuracy of the reconstruction, and the Lagrange multiplication operator λt, which guarantees the tightness of the constraints, are introduced; its extended Lagrangian expression is as follows:(3)Luk,wk,λ=α∑k=1K∂tδt+jπt∗ukte−jwkt22+ft−∑k=1Kukt22+λt,ft−∑k=1Kukt.On this basis, iterative updates ofukn+1, wkn+1, and λn+1 using the Alternate Direction Method of Multipliers (ADMM) are used to find the “saddle point” of the extended Lagrange expression, i.e., the optimal solution of the constrained-variance model [28]. The VMD algorithm decomposes the raw signal into K modal components as follows:Step 1.
Initialize parameters:u^k1,w^k1,λ^k1,n.Step 2.
Perform an iterative loopn=n+1.Step 3
. The IMFu^kw is updated for all w≥0:(4)u^kn+1w=f^w−∑i≠ku^inw+λ^inw/21+2αw−wkn2.Step 4.
Update the center frequencywk of the modal.(5)w^kn+1=∫0∞wu^kn+1w2dw∫0∞u^kn+1w2dw.Step 5.
Update Lagrange multiplication operatorλ^w.(6)λ^n+1w=λ^nw+τf^w−∑ku^kn+1w.Step 6.
Repeat Steps 2–5 until the inequality converges.(7)∑ku^kn+1−u^kn22u^kn22<ε.Step 7.
At the end,K IMFs are obtained, where n is the number of iterations. In this method, the input parameters play an important role; for example, the balance parameter (α) which is the penalty factor, the fidelity coefficient τ, the number of DC components, w (init) which is the initial frequency, and the threshold ε(ε > 0).


### 2.3. Statistical-Based Feature Extraction

In the research work of this paper, statistical parameters are used as features. The specific expressions and details of statistical characteristics are as follows.


#### 2.3.1. Variance Contribution Rate

The VCR of the IMF component, that is, the VCR of each IMF component, is calculated by calculating the size of the variance of each IMF component, which can be calculated by(8)Mi=Di∑i=1nDi,(9)Di=1N∑k=1NcikΔt2−1N∑k=1NcikΔt2.Mi is the VCR of the IMF component; ci is the IMF component; Di is the variance of the IMF component; N is the number of IMF components. According to statistical theory, the VCR represents the relative importance of each IMF component to the decomposed signal, and the larger the Mi value of the VCR of the IMF, the greater the effect of the IMF component on the decomposed signal and, conversely, the smaller.


#### 2.3.2. Energy Entropy

Entropy is a concept derived from physics, which is a measure of the degree of disorder within the system. The larger the entropy value, the more complicated the system. In a normal working state or a working state with a fault, the energy in the vibration signal is distributed in different frequency bands [29].Decompose the fault signal by VMD, and then calculate the energy of each fault IMF component:E1,E2,E3,…,En, and then the energy distribution in the vibration signal can be obtained. In the case of ignoring the remainder, according to the principle of orthogonality, the total energy of the decomposed signal should be substantially equal to the sum of the energies of the n IMF components. The intrinsic modal energy entropy is defined as follows:(10)HEN=−∑i=1npilogpi.In the formula,pi=Ei/E, which represents the proportion of the energy of the IMF in the total energy.


#### 2.3.3. Permutation Entropy

In order to further increase the distinction between different fault states and facilitate the subsequent classification, the concept of permutation entropy is introduced. Permutation entropy is an average entropy parameter to measure the complexity of one-dimensional time series. Bandt and Pompe [30] proposed to use PE to detect the dynamic changes of random signals in 2002. Compared with LyaPunov index, fractal dimension, and other complexity parameters, it has the advantages of simple calculation, strong robustness, and stable calculation value.The principle of the PE algorithm is as follows: consider a time seriesXi,i=1,2,…,N; the phase space reconstruction is performed using the embedding dimension m and the delay time τ to obtain the matrix:(11)x1x1+τ⋯x1+m−1τ⋮⋮⋮⋮xjxj+τ⋯xj+m−1τ⋮⋮⋮⋮xGxG+τ⋯xG+m−1τ.G is the number of reconstructed vectors in the reconstructed phase space, G=N−m−1τ. Rearrange the components xj,xj+τ,…,xj+m−1τ in ascending order as follows: xi+j1−1τ≤⋯≤xi+jm−1τ. In the formula, j1,j2,…,jm is the index of the column where each element in the reconstructed component is located. If there are equal values in the reconstructed component, such as xi+jp−1τ=xi+jq−1τ, then, according to the original order of jp and jq, if jp<jq, there is xi+jp−1τ≤xi+jq−1τ. For any reconstruction vector Xj, a symbol sequence Sl=j1,j2,…,jm reflecting the order of its elements can be obtained, where l=1,2,…,g, and g≤m!. The m-dimensional phase space maps different symbol sequences j1,j2,…,jm with m! in total, and Sl is one of them. The sequences P1,P2,…,Pg are constructed, and Pg is the probability of the occurrence of the symbol sequence. For the PE of the symbol sequence corresponding to the g reconstruction vectors of a time series, it can be defined as(12)PEm=−∑l=1gPllnPl.It is noticed thatPEm reaches the maximum value, lnm! when Pl=1/m!. Normally, lnm! can be used for standardization; namely,(13)PE=PEmlnm!.In the formula, the value range ofPE is0,1. The larger the value of PE, the more random the time series; the smaller the PE value, the more regular the time series.


### 2.4. Classifier of SVM

Cortes and Vapnik [31] proposed a support vector machine (SVM) algorithm based on statistical learning theory, which is a learning method for classification and regression and has the ability to handle high input features.When using SVM for classification, the SVM first maps the input vector to a higher feature space and constructs a maximum separation hyperplane in the linear or nonlinear feature space. Two parallel hyperplanes are constructed on each side of the hyperplane separating the data. Then, use the support vector to determine the optimal hyperplane. This hyperplane can be linear or nonlinear [32]. However, in most cases, the mode is not linearly separable, which is related to the kernel function.In this study, through the work of the above steps, we can get many characteristic samples from these collected vibration signals. Then, the extracted samples are divided into two groups: training samples are used to train the support vector machine model, and test samples are used to prove the effectiveness of the proposed method. During the training process, the support vector machine creates a hyperplane and assigns most points of the same class on the same side. Studies have shown that radial basis function (RBF) kernels give better processing results than other kernel functions [33]. Therefore, this study uses SVM with RBF kernel functions.


## 2.1. Wavelet Packet-Based Denoising

In practical engineering applications, the signals generated by mechanical equipment are nonstationary signals composed of useful signals and noise signals. When performing noise processing, signal decomposition, threshold quantization processing, and signal reconstruction are required. In signal decomposition, wavelet transform and wavelet packet transform have different decomposition characteristics.


### 2.1.1. Wavelet Decomposition

Wavelet analysis is a time-frequency localized analysis method developed on the basis of Fourier transform, with a fixed window size, but its shape can be changed. It is generally used to process nonstationary signals. When the wavelet transform decomposes the signal, only the low-frequency part is further decomposed, and the high-frequency part, that is, the detailed part of the signal, is no longer decomposed. Therefore, the wavelet transform can well characterize the signal with the low-frequency signal as the main component but cannot decompose and represent signals containing a large amount of detailed information, such as nonstationary mechanical vibration signals, biomedical signals, and seismic signals. The schematic diagram of wavelet decomposition is shown in Figure1, in which three layers of wavelet decomposition are done on the original signal S, A represents the approximate value of the signal, and D represents the detailed value of the signal.Figure 1
The 3-layer wavelet decomposition.The choice of threshold rules (Minimax, Sqtwlolg, Heursure, Rigrsure), decomposed layers (generally 3 to 5 layers), and wavelet basis functions (Haar, Mexican Hat (mexh), Morlets, Gaussian, Daubechies (dbN), Meyer, Symlets (sym), Coiflet (coifN), etc.) will have a certain effect on the noise reduction effect.When the signal continuity and smoothness is good, it is appropriate to use the Sym8 wavelet, mainly because the wavelet basis of Sym8 has good smoothness and continuity. When the signal continuity is weak, the Haar wavelet basis should be used. The db wavelet base is more suitable for the unsteady signals such as structural vibration and microseismic signals.


### 2.1.2. Wavelet Packet Decomposition

Wavelet packet decomposition is an improvement on wavelet decomposition and is more refined than wavelet decomposition. The wavelet packet transform divides the signal more finely. It can decompose the low-frequency part and the high-frequency part of the signal and can adaptively select the corresponding frequency band according to the characteristics of the analyzed signal to make it match the new frequency spectrum, thereby improving the time-frequency resolution rate. At present, wavelet packet analysis has been widely used in data denoising, vibration information extraction, and other monitoring fields [26, 27].The following are the specific steps of wavelet packet noise reduction: (1) the signal is decomposed byN-layer wavelet packet. The schematic diagram of wavelet packet decomposition is shown in Figure 2, in which three layers of wavelet packet decomposition are done on the original signal S; (2) the minimum cost principle is used to calculate the optimal wavelet packet base; (3) threshold quantization of wavelet packet decomposition coefficient. For each wavelet packet decomposition coefficient, select an appropriate threshold value for threshold quantization of left and right wavelet packet bases; (4) the reconstructed signal is the signal after the optimal wavelet packet noise reduction.Figure 2
The wavelet packet 3-layer decomposition.Based on the characteristics of wavelet packet, wavelet packet decomposition is used to denoise the collected vibration signals.


## 2.1.1. Wavelet Decomposition

Wavelet analysis is a time-frequency localized analysis method developed on the basis of Fourier transform, with a fixed window size, but its shape can be changed. It is generally used to process nonstationary signals. When the wavelet transform decomposes the signal, only the low-frequency part is further decomposed, and the high-frequency part, that is, the detailed part of the signal, is no longer decomposed. Therefore, the wavelet transform can well characterize the signal with the low-frequency signal as the main component but cannot decompose and represent signals containing a large amount of detailed information, such as nonstationary mechanical vibration signals, biomedical signals, and seismic signals. The schematic diagram of wavelet decomposition is shown in Figure1, in which three layers of wavelet decomposition are done on the original signal S, A represents the approximate value of the signal, and D represents the detailed value of the signal.Figure 1
The 3-layer wavelet decomposition.The choice of threshold rules (Minimax, Sqtwlolg, Heursure, Rigrsure), decomposed layers (generally 3 to 5 layers), and wavelet basis functions (Haar, Mexican Hat (mexh), Morlets, Gaussian, Daubechies (dbN), Meyer, Symlets (sym), Coiflet (coifN), etc.) will have a certain effect on the noise reduction effect.When the signal continuity and smoothness is good, it is appropriate to use the Sym8 wavelet, mainly because the wavelet basis of Sym8 has good smoothness and continuity. When the signal continuity is weak, the Haar wavelet basis should be used. The db wavelet base is more suitable for the unsteady signals such as structural vibration and microseismic signals.


## 2.1.2. Wavelet Packet Decomposition

Wavelet packet decomposition is an improvement on wavelet decomposition and is more refined than wavelet decomposition. The wavelet packet transform divides the signal more finely. It can decompose the low-frequency part and the high-frequency part of the signal and can adaptively select the corresponding frequency band according to the characteristics of the analyzed signal to make it match the new frequency spectrum, thereby improving the time-frequency resolution rate. At present, wavelet packet analysis has been widely used in data denoising, vibration information extraction, and other monitoring fields [26, 27].The following are the specific steps of wavelet packet noise reduction: (1) the signal is decomposed byN-layer wavelet packet. The schematic diagram of wavelet packet decomposition is shown in Figure 2, in which three layers of wavelet packet decomposition are done on the original signal S; (2) the minimum cost principle is used to calculate the optimal wavelet packet base; (3) threshold quantization of wavelet packet decomposition coefficient. For each wavelet packet decomposition coefficient, select an appropriate threshold value for threshold quantization of left and right wavelet packet bases; (4) the reconstructed signal is the signal after the optimal wavelet packet noise reduction.Figure 2
The wavelet packet 3-layer decomposition.Based on the characteristics of wavelet packet, wavelet packet decomposition is used to denoise the collected vibration signals.


## 2.2. Variational Mode Decomposition

The essence of the VMD method is an all-modal variant that adaptively decomposes the original signalfx into a modal function uk with a central frequency of wk for K, where K is the predetermined number of components. Assuming that fx is a signal with a multi-frequency component, the VMD algorithm redefines IMF as a signal that is both amplitude- and frequency-modulated, with the expression(1)ukt=Aktcosϕkt,where ϕkt is the phase and Akt is the envelope magnitude of the signal. The constraint variance model can be expressed as(2)minuk,wk∑k=1K∂tδt+jπt∗ukt∗e−jwkt22,where ukt satisfies ∑k=1Kuk=fx, δt is the Dirac delta function, and “∗” is the convolutional sign; uk is the set of components for IMFs; wk is the set of the corresponding central frequency bands for K components; and ∂t denotes the bias.In order to make the constrained-variance problem unconstrained-variance problem, the secondary penalty factorα, which guarantees the accuracy of the reconstruction, and the Lagrange multiplication operator λt, which guarantees the tightness of the constraints, are introduced; its extended Lagrangian expression is as follows:(3)Luk,wk,λ=α∑k=1K∂tδt+jπt∗ukte−jwkt22+ft−∑k=1Kukt22+λt,ft−∑k=1Kukt.On this basis, iterative updates ofukn+1, wkn+1, and λn+1 using the Alternate Direction Method of Multipliers (ADMM) are used to find the “saddle point” of the extended Lagrange expression, i.e., the optimal solution of the constrained-variance model [28]. The VMD algorithm decomposes the raw signal into K modal components as follows:Step 1.
Initialize parameters:u^k1,w^k1,λ^k1,n.Step 2.
Perform an iterative loopn=n+1.Step 3
. The IMFu^kw is updated for all w≥0:(4)u^kn+1w=f^w−∑i≠ku^inw+λ^inw/21+2αw−wkn2.Step 4.
Update the center frequencywk of the modal.(5)w^kn+1=∫0∞wu^kn+1w2dw∫0∞u^kn+1w2dw.Step 5.
Update Lagrange multiplication operatorλ^w.(6)λ^n+1w=λ^nw+τf^w−∑ku^kn+1w.Step 6.
Repeat Steps 2–5 until the inequality converges.(7)∑ku^kn+1−u^kn22u^kn22<ε.Step 7.
At the end,K IMFs are obtained, where n is the number of iterations. In this method, the input parameters play an important role; for example, the balance parameter (α) which is the penalty factor, the fidelity coefficient τ, the number of DC components, w (init) which is the initial frequency, and the threshold ε(ε > 0).


## 2.3. Statistical-Based Feature Extraction

In the research work of this paper, statistical parameters are used as features. The specific expressions and details of statistical characteristics are as follows.


### 2.3.1. Variance Contribution Rate

The VCR of the IMF component, that is, the VCR of each IMF component, is calculated by calculating the size of the variance of each IMF component, which can be calculated by(8)Mi=Di∑i=1nDi,(9)Di=1N∑k=1NcikΔt2−1N∑k=1NcikΔt2.Mi is the VCR of the IMF component; ci is the IMF component; Di is the variance of the IMF component; N is the number of IMF components. According to statistical theory, the VCR represents the relative importance of each IMF component to the decomposed signal, and the larger the Mi value of the VCR of the IMF, the greater the effect of the IMF component on the decomposed signal and, conversely, the smaller.


### 2.3.2. Energy Entropy

Entropy is a concept derived from physics, which is a measure of the degree of disorder within the system. The larger the entropy value, the more complicated the system. In a normal working state or a working state with a fault, the energy in the vibration signal is distributed in different frequency bands [29].Decompose the fault signal by VMD, and then calculate the energy of each fault IMF component:E1,E2,E3,…,En, and then the energy distribution in the vibration signal can be obtained. In the case of ignoring the remainder, according to the principle of orthogonality, the total energy of the decomposed signal should be substantially equal to the sum of the energies of the n IMF components. The intrinsic modal energy entropy is defined as follows:(10)HEN=−∑i=1npilogpi.In the formula,pi=Ei/E, which represents the proportion of the energy of the IMF in the total energy.


### 2.3.3. Permutation Entropy

In order to further increase the distinction between different fault states and facilitate the subsequent classification, the concept of permutation entropy is introduced. Permutation entropy is an average entropy parameter to measure the complexity of one-dimensional time series. Bandt and Pompe [30] proposed to use PE to detect the dynamic changes of random signals in 2002. Compared with LyaPunov index, fractal dimension, and other complexity parameters, it has the advantages of simple calculation, strong robustness, and stable calculation value.The principle of the PE algorithm is as follows: consider a time seriesXi,i=1,2,…,N; the phase space reconstruction is performed using the embedding dimension m and the delay time τ to obtain the matrix:(11)x1x1+τ⋯x1+m−1τ⋮⋮⋮⋮xjxj+τ⋯xj+m−1τ⋮⋮⋮⋮xGxG+τ⋯xG+m−1τ.G is the number of reconstructed vectors in the reconstructed phase space, G=N−m−1τ. Rearrange the components xj,xj+τ,…,xj+m−1τ in ascending order as follows: xi+j1−1τ≤⋯≤xi+jm−1τ. In the formula, j1,j2,…,jm is the index of the column where each element in the reconstructed component is located. If there are equal values in the reconstructed component, such as xi+jp−1τ=xi+jq−1τ, then, according to the original order of jp and jq, if jp<jq, there is xi+jp−1τ≤xi+jq−1τ. For any reconstruction vector Xj, a symbol sequence Sl=j1,j2,…,jm reflecting the order of its elements can be obtained, where l=1,2,…,g, and g≤m!. The m-dimensional phase space maps different symbol sequences j1,j2,…,jm with m! in total, and Sl is one of them. The sequences P1,P2,…,Pg are constructed, and Pg is the probability of the occurrence of the symbol sequence. For the PE of the symbol sequence corresponding to the g reconstruction vectors of a time series, it can be defined as(12)PEm=−∑l=1gPllnPl.It is noticed thatPEm reaches the maximum value, lnm! when Pl=1/m!. Normally, lnm! can be used for standardization; namely,(13)PE=PEmlnm!.In the formula, the value range ofPE is0,1. The larger the value of PE, the more random the time series; the smaller the PE value, the more regular the time series.


## 2.3.1. Variance Contribution Rate

The VCR of the IMF component, that is, the VCR of each IMF component, is calculated by calculating the size of the variance of each IMF component, which can be calculated by(8)Mi=Di∑i=1nDi,(9)Di=1N∑k=1NcikΔt2−1N∑k=1NcikΔt2.Mi is the VCR of the IMF component; ci is the IMF component; Di is the variance of the IMF component; N is the number of IMF components. According to statistical theory, the VCR represents the relative importance of each IMF component to the decomposed signal, and the larger the Mi value of the VCR of the IMF, the greater the effect of the IMF component on the decomposed signal and, conversely, the smaller.


## 2.3.2. Energy Entropy

Entropy is a concept derived from physics, which is a measure of the degree of disorder within the system. The larger the entropy value, the more complicated the system. In a normal working state or a working state with a fault, the energy in the vibration signal is distributed in different frequency bands [29].Decompose the fault signal by VMD, and then calculate the energy of each fault IMF component:E1,E2,E3,…,En, and then the energy distribution in the vibration signal can be obtained. In the case of ignoring the remainder, according to the principle of orthogonality, the total energy of the decomposed signal should be substantially equal to the sum of the energies of the n IMF components. The intrinsic modal energy entropy is defined as follows:(10)HEN=−∑i=1npilogpi.In the formula,pi=Ei/E, which represents the proportion of the energy of the IMF in the total energy.


## 2.3.3. Permutation Entropy

In order to further increase the distinction between different fault states and facilitate the subsequent classification, the concept of permutation entropy is introduced. Permutation entropy is an average entropy parameter to measure the complexity of one-dimensional time series. Bandt and Pompe [30] proposed to use PE to detect the dynamic changes of random signals in 2002. Compared with LyaPunov index, fractal dimension, and other complexity parameters, it has the advantages of simple calculation, strong robustness, and stable calculation value.The principle of the PE algorithm is as follows: consider a time seriesXi,i=1,2,…,N; the phase space reconstruction is performed using the embedding dimension m and the delay time τ to obtain the matrix:(11)x1x1+τ⋯x1+m−1τ⋮⋮⋮⋮xjxj+τ⋯xj+m−1τ⋮⋮⋮⋮xGxG+τ⋯xG+m−1τ.G is the number of reconstructed vectors in the reconstructed phase space, G=N−m−1τ. Rearrange the components xj,xj+τ,…,xj+m−1τ in ascending order as follows: xi+j1−1τ≤⋯≤xi+jm−1τ. In the formula, j1,j2,…,jm is the index of the column where each element in the reconstructed component is located. If there are equal values in the reconstructed component, such as xi+jp−1τ=xi+jq−1τ, then, according to the original order of jp and jq, if jp<jq, there is xi+jp−1τ≤xi+jq−1τ. For any reconstruction vector Xj, a symbol sequence Sl=j1,j2,…,jm reflecting the order of its elements can be obtained, where l=1,2,…,g, and g≤m!. The m-dimensional phase space maps different symbol sequences j1,j2,…,jm with m! in total, and Sl is one of them. The sequences P1,P2,…,Pg are constructed, and Pg is the probability of the occurrence of the symbol sequence. For the PE of the symbol sequence corresponding to the g reconstruction vectors of a time series, it can be defined as(12)PEm=−∑l=1gPllnPl.It is noticed thatPEm reaches the maximum value, lnm! when Pl=1/m!. Normally, lnm! can be used for standardization; namely,(13)PE=PEmlnm!.In the formula, the value range ofPE is0,1. The larger the value of PE, the more random the time series; the smaller the PE value, the more regular the time series.


## 2.4. Classifier of SVM

Cortes and Vapnik [31] proposed a support vector machine (SVM) algorithm based on statistical learning theory, which is a learning method for classification and regression and has the ability to handle high input features.When using SVM for classification, the SVM first maps the input vector to a higher feature space and constructs a maximum separation hyperplane in the linear or nonlinear feature space. Two parallel hyperplanes are constructed on each side of the hyperplane separating the data. Then, use the support vector to determine the optimal hyperplane. This hyperplane can be linear or nonlinear [32]. However, in most cases, the mode is not linearly separable, which is related to the kernel function.In this study, through the work of the above steps, we can get many characteristic samples from these collected vibration signals. Then, the extracted samples are divided into two groups: training samples are used to train the support vector machine model, and test samples are used to prove the effectiveness of the proposed method. During the training process, the support vector machine creates a hyperplane and assigns most points of the same class on the same side. Studies have shown that radial basis function (RBF) kernels give better processing results than other kernel functions [33]. Therefore, this study uses SVM with RBF kernel functions.


## 3. Proposed Methodology

The flowchart of the rolling bearing fault diagnosis method proposed in this paper is shown in Figure3. As shown in the figure, the method includes three steps: data processing, feature extraction, and model-based fault diagnosis for classification identification.Figure 3
The flowchart of the proposed method.The first step is to use the acceleration sensor data collection system to collect the vibration signals of the healthy rolling bearings of the spindle equipment and the rolling bearings of different damaged parts. And noise reduction processing is performed on the collected signal. In the second step, the VMD is used to adaptively decompose the vibration signal into K IMFs without modal mixing and non-IMF problems. The IMFs of VMD are concentrated around the center frequencyw. Statistical characteristics are extracted from these IMFs through statistical analysis and normalized. Finally, the training samples and test samples are constructed according to the used feature parameter matrix. The third step is to use the training samples to design a classification model for classification prediction and fault diagnosis. Experiments were carried out on the elevator mechanical failure simulator to verify the effectiveness of the method. Experimental results show that the method is robust to diagnose known faults and monitor unknown faults.


## 4. Experiment Analysis

In order to verify the effectiveness of the method of VMD combined with VCR, EE, and PE, this paper first selects the rolling bearing test data of Case Western Reserve University’s Electrical Engineering Laboratory [34] for analysis to prove the effectiveness of the method. In the past decades, the bearing data has become a standard reference used to test the performance of the algorithms proposed by the fault diagnosis researchers [35].Then, this method is applied to the mechanical failure simulation experiment platform of the independently developed hoist spindle device.


### 4.1. Experiment 1

The electrical engineering experiment platform of Case Western Reserve University is shown in Figure4. The left side is a motor, the middle part is a torque sensor, and the right side is a dynamometer. The loading power of the rotating motor is 0 hp, 1 hp, 2 hp, and 3 hp (1 hp = 736 W), and the 2 hp was selected for test analysis first. The test uses EDM technology to process a single point of failure on the bearing. The diameter of the fault point is 0.1778 mm. The bearing speed is 1750 r/min. At the sampling frequency of 12 kHz, the vibration signal data of the driving end of the normal state, the inner ring fault, the outer ring fault (6.00 o’clock position), and the rolling element fault were collected, respectively, and the experimental data was analyzed.Figure 4
The mechanical failure simulation test bench.81920 data points were intercepted for the four states, respectively, and divided into 40 groups, each group having a length of 2048. A total of 160 sets of experimental data samples were obtained in all states. Figure5 is a time-domain diagram of four state data sequences. It is obviously not intuitive to see that the bearing has failed.Figure 5
The time-domain diagram of four states of rolling bearing.


#### 4.1.1. Determination of the Number of Modes

Whether the choice of the decomposition numberK is reasonable will have different effects on the decomposition result. Assuming that the number of IMF components K is less than the number of useful components in the processed signal, it will result in insufficient decomposition; on the contrary, it will make the center frequencies of adjacent IMF components close, resulting in frequency aliasing. The rolling element fault signal is used for VMD decomposition. The center frequency under different K values is shown in Table 1.Table 1
The center frequency of components with differentK values.
KCenter frequency (Hz)363528813366458613332884336655571102283832603430655711012838325934284164It can be seen from the table that when the number of components is 5, the center frequencies of 3260 Hz and 3430 Hz are relatively close, and modal aliasing may occur, so it is more appropriate to select 4 as the number of components. Select the first set of data of the inner circle fault, and use the input parametersα = 2000, DC = 0, K = 4, t = 0, tol = 1e−7, w (init) = 1, ε = 1e−7 to perform VMD decomposition.The time-domain diagram and corresponding spectrum diagram after the data is decomposed by VMD are shown in Figures6 and 7.Figure 6
The time-domain diagram after inner ring fault signal decomposition.Figure 7
The frequency-domain diagram after inner ring fault signal decomposition.


#### 4.1.2. Feature Extraction

According to formulae (8), (10), and (12), extract the Mi, EE(HEN), PE of each IMF component to construct a 160∗12 feature vector matrix. After normalizing the data, the box plot is shown in Figure 8. It can be seen from the figure that the fault features all exhibit good clustering characteristics, which are all between 0 and 1, and there are no abnormal data points.Figure 8
The box diagram of fault data.


#### 4.1.3. Bearing Fault Identification and Classification

The SVM model uses a cross-validation search algorithm to select the gamma functiong and the penalty coefficient c in the kernel function, and select c = 0.0625 and g = 0.0625. There are 160 sets of data samples, of which the normal state number is 1–40, the inner ring fault number is 41–80, the outer ring fault number is 81–120, and the rolling element fault number is 121–160. Meanwhile, define labels on the sample data, respectively; 1 represents the normal state, 2 represents the inner circle failure, 3 represents the outer circle failure, 4 represents rolling element failure. Three-quarters of the data were randomly selected from each state sample and the order was shuffled. The corresponding labels are also scrambled accordingly. The data was input into SVM classifier for training, and the SVM model for fault data classification was obtained. The remaining sample data is used as test data. Use the trained SVM model to classify and predict the disordered test data. The classification results of the tests are shown in Figure 9. It can be seen from the figure that the prediction test set classification results are the same as the actual test set classification. The confusion matrix of the classification results is presented in Figure 10. The columns stand for the actual label of each condition, and the rows stand for the predict label. It shows that the four conditions have the accuracy of 100%.Figure 9
The predicted results with VMD.Figure 10
The confusion matrix of the result with VMD.In order to further verify the effectiveness of the VMD method, the EMD and MEEMD algorithms are used to process the sample data of the four fault states. Meanwhile, single or multiple statistical features are selected as feature parameters and combined with decomposition algorithms. The correlation coefficient method was used to select the IMF components of the first four orders obtained from the decomposition of each sample data and extract the characteristics. The diagnostic process is the same as above.It can be seen from Table2 that, for the combination of VCR, EE, and PE, the accuracy rates of the three have reached 100%. However, for individual characteristics VCR, EE, and PE, VMD has the highest accuracy rate of 100%. For the combination of VCR and EE, VMD also has the highest accuracy at 100%.Table 2
The classification results with SVM (2 hp).
FeaturesAccuracy of classification (%)VMDEMDMEEMDVCR100.0100.097.5EE100.0100.085.0PE100.095.077.5VCR + EE100.095.097.5VCR + EE + PE100.0100.0100.0For comparison, BP neural network is used to classify and recognize the same feature parameters. The learning rate of BP neural network is 0.1, and the number of iterations is 50. The numbers of input layers, hidden layers, and output layers are 4, 5, 4 (VCR, EE, PE), 8, 9, 4 (VCR + EE), and 12, 13, 4 (VCR + EE + PE). The classification and recognition results are shown in Table3.Table 3
The classification results with BP network (2 hp).
FeaturesAccuracy of classification (%)VMDEMDMEEMDVCR100.088.7591.25EE100.098.085.25PE100.094.2577.35VCR + EE100.098.090.5VCR + EE + PE100.0100%97.25Through the comparison of the values in Table2 and 3, for the same feature parameters, the SVM classifier is generally better than the BP neural network classification effect. Meanwhile, test experiments with different motor loads (0 hp and 1 hp) were added, and the classification results are shown in Table 4. The training models of 0 hp and 1 hp were used to test and verify the characteristic parameters of 3 hp, and the classification results are shown in Table 5.Table 4
The classification results of 0 hp and 1 hp (SVM).
FeaturesAccuracy of classification (%)0 hp1 hpVMDEMDMEEMDVMDEMDMEEMDVCR100.097.592.5100.0100.0100.0EE100.097.585.0100.0100.095.0PE97.5100.077.5100.095.087.5VCR + EE100.0100.085.0100.0100.097.5VCR + EE + PE100.0100.0100.0100.0100.0100.0Table 5
The test results of 3 hp under the training model of 0 hp and 1 hp (SVM).
FeaturesClassification accuracy of 3 hp (%)0 hp training model1 hp training modelVMDEMDMEEMDVMDEMDMEEMDVCR71.258043.7594.37576.87535.0EE99.3757534.37599.37576.2535.625PE9591.2558.7596.2593.12545.625VCR + EE76.2575.62536.2599.37576.87538.125VCR + EE + PE98.7598.12538.12599.37599.37536.875Through the comparison of the values in Tables4 and 5, it can be seen that the test classification accuracy of 3 hp reaches 99.375%, indicating that the training classification model has certain robustness. Therefore, based on the VMD method to decompose the signal and extract the characteristics, it is more effective in diagnosing rolling bearing failure.


### 4.2. Experiment 2

The mine hoist is mainly composed of driving motor, spindle device, hydraulic and braking system, electrical control system, and depth indicator. The spindle device, which bears all the dynamic and static loads acting on the spindle device through the wire rope and the transmission torque of the drive motor, is the key rotating body supporting the mine hoist and is mainly composed of the reel, the spindle, the main bearing, the bearing end cover, and so on. The VMD method in the previous section is applied to the mechanical failure simulation experiment platform of the hoist spindle device in Figure11.Figure 11
The fault simulation test bench.The test bench selects the Panasonic model MSMF082 L wire-type motor with rated speed of 3000 r/min and output power of 750 W. The servo motor speed is set at 20 r/min. The model of AC192-1D vibration acceleration sensor produced by American CTC vibration acceleration sensor company is used to record the fault vibration signal of the bearing under different conditions: normal state, inner ring damage, outer ring damage, and rolling element damage. The model of acquisition card is Smacq-USB-4000, and the vibration signal is collected at a sampling frequency of 5000 Hz.In order to simulate the working condition of the lifting load of the hoist, the test bench uses a magnetic powder dynamometer to apply a load to the steel wire rope. The size of the load is adjusted and displayed by the NX6002 speed torque measuring instrument. Under different fault conditions, the torque of the dynamometer is adjusted to 20 N·m. A single point of failure (width 1.5 mm, depth 1.5 mm) was processed on the outer and inner rings of the bearing by electrical discharge machining, and a fine slit (width 1 mm, depth 1 mm) was cut on the roller by wire cutting. The relevant parameters of the bearing are shown in Table6. Figure 12 shows the four failure states of rolling bearings.Table 6
The specifications of the test bearing.
FeatureDetailBearing typeSKF22207CA/W33Roller diameter d (mm)8.2Pitch diameter D (mm)55Number of balls (N)30Contact angleθ°0Figure 12
The object diagram of the roller bearing: (a) outer ring fault; (b) inner ring fault; (c) rolling element fault; (d) normal.
(a)(b)(c)(d)The vibration signal collected by the sensor near the servo motor driver is selected as the analysis signal, and the wavelet packet is used for noise reduction. The 81920 data points were intercepted from the four states of the bearing and divided into 40 groups, with each length 2048. A total of 160 data samples were obtained. Figure13 shows the time-domain diagram of bearing vibration signals in four states.Figure 13
The time-domain diagram of four states of rolling bearing.


## 4.1. Experiment 1

The electrical engineering experiment platform of Case Western Reserve University is shown in Figure4. The left side is a motor, the middle part is a torque sensor, and the right side is a dynamometer. The loading power of the rotating motor is 0 hp, 1 hp, 2 hp, and 3 hp (1 hp = 736 W), and the 2 hp was selected for test analysis first. The test uses EDM technology to process a single point of failure on the bearing. The diameter of the fault point is 0.1778 mm. The bearing speed is 1750 r/min. At the sampling frequency of 12 kHz, the vibration signal data of the driving end of the normal state, the inner ring fault, the outer ring fault (6.00 o’clock position), and the rolling element fault were collected, respectively, and the experimental data was analyzed.Figure 4
The mechanical failure simulation test bench.81920 data points were intercepted for the four states, respectively, and divided into 40 groups, each group having a length of 2048. A total of 160 sets of experimental data samples were obtained in all states. Figure5 is a time-domain diagram of four state data sequences. It is obviously not intuitive to see that the bearing has failed.Figure 5
The time-domain diagram of four states of rolling bearing.


### 4.1.1. Determination of the Number of Modes

Whether the choice of the decomposition numberK is reasonable will have different effects on the decomposition result. Assuming that the number of IMF components K is less than the number of useful components in the processed signal, it will result in insufficient decomposition; on the contrary, it will make the center frequencies of adjacent IMF components close, resulting in frequency aliasing. The rolling element fault signal is used for VMD decomposition. The center frequency under different K values is shown in Table 1.Table 1
The center frequency of components with differentK values.
KCenter frequency (Hz)363528813366458613332884336655571102283832603430655711012838325934284164It can be seen from the table that when the number of components is 5, the center frequencies of 3260 Hz and 3430 Hz are relatively close, and modal aliasing may occur, so it is more appropriate to select 4 as the number of components. Select the first set of data of the inner circle fault, and use the input parametersα = 2000, DC = 0, K = 4, t = 0, tol = 1e−7, w (init) = 1, ε = 1e−7 to perform VMD decomposition.The time-domain diagram and corresponding spectrum diagram after the data is decomposed by VMD are shown in Figures6 and 7.Figure 6
The time-domain diagram after inner ring fault signal decomposition.Figure 7
The frequency-domain diagram after inner ring fault signal decomposition.


### 4.1.2. Feature Extraction

According to formulae (8), (10), and (12), extract the Mi, EE(HEN), PE of each IMF component to construct a 160∗12 feature vector matrix. After normalizing the data, the box plot is shown in Figure 8. It can be seen from the figure that the fault features all exhibit good clustering characteristics, which are all between 0 and 1, and there are no abnormal data points.Figure 8
The box diagram of fault data.


### 4.1.3. Bearing Fault Identification and Classification

The SVM model uses a cross-validation search algorithm to select the gamma functiong and the penalty coefficient c in the kernel function, and select c = 0.0625 and g = 0.0625. There are 160 sets of data samples, of which the normal state number is 1–40, the inner ring fault number is 41–80, the outer ring fault number is 81–120, and the rolling element fault number is 121–160. Meanwhile, define labels on the sample data, respectively; 1 represents the normal state, 2 represents the inner circle failure, 3 represents the outer circle failure, 4 represents rolling element failure. Three-quarters of the data were randomly selected from each state sample and the order was shuffled. The corresponding labels are also scrambled accordingly. The data was input into SVM classifier for training, and the SVM model for fault data classification was obtained. The remaining sample data is used as test data. Use the trained SVM model to classify and predict the disordered test data. The classification results of the tests are shown in Figure 9. It can be seen from the figure that the prediction test set classification results are the same as the actual test set classification. The confusion matrix of the classification results is presented in Figure 10. The columns stand for the actual label of each condition, and the rows stand for the predict label. It shows that the four conditions have the accuracy of 100%.Figure 9
The predicted results with VMD.Figure 10
The confusion matrix of the result with VMD.In order to further verify the effectiveness of the VMD method, the EMD and MEEMD algorithms are used to process the sample data of the four fault states. Meanwhile, single or multiple statistical features are selected as feature parameters and combined with decomposition algorithms. The correlation coefficient method was used to select the IMF components of the first four orders obtained from the decomposition of each sample data and extract the characteristics. The diagnostic process is the same as above.It can be seen from Table2 that, for the combination of VCR, EE, and PE, the accuracy rates of the three have reached 100%. However, for individual characteristics VCR, EE, and PE, VMD has the highest accuracy rate of 100%. For the combination of VCR and EE, VMD also has the highest accuracy at 100%.Table 2
The classification results with SVM (2 hp).
FeaturesAccuracy of classification (%)VMDEMDMEEMDVCR100.0100.097.5EE100.0100.085.0PE100.095.077.5VCR + EE100.095.097.5VCR + EE + PE100.0100.0100.0For comparison, BP neural network is used to classify and recognize the same feature parameters. The learning rate of BP neural network is 0.1, and the number of iterations is 50. The numbers of input layers, hidden layers, and output layers are 4, 5, 4 (VCR, EE, PE), 8, 9, 4 (VCR + EE), and 12, 13, 4 (VCR + EE + PE). The classification and recognition results are shown in Table3.Table 3
The classification results with BP network (2 hp).
FeaturesAccuracy of classification (%)VMDEMDMEEMDVCR100.088.7591.25EE100.098.085.25PE100.094.2577.35VCR + EE100.098.090.5VCR + EE + PE100.0100%97.25Through the comparison of the values in Table2 and 3, for the same feature parameters, the SVM classifier is generally better than the BP neural network classification effect. Meanwhile, test experiments with different motor loads (0 hp and 1 hp) were added, and the classification results are shown in Table 4. The training models of 0 hp and 1 hp were used to test and verify the characteristic parameters of 3 hp, and the classification results are shown in Table 5.Table 4
The classification results of 0 hp and 1 hp (SVM).
FeaturesAccuracy of classification (%)0 hp1 hpVMDEMDMEEMDVMDEMDMEEMDVCR100.097.592.5100.0100.0100.0EE100.097.585.0100.0100.095.0PE97.5100.077.5100.095.087.5VCR + EE100.0100.085.0100.0100.097.5VCR + EE + PE100.0100.0100.0100.0100.0100.0Table 5
The test results of 3 hp under the training model of 0 hp and 1 hp (SVM).
FeaturesClassification accuracy of 3 hp (%)0 hp training model1 hp training modelVMDEMDMEEMDVMDEMDMEEMDVCR71.258043.7594.37576.87535.0EE99.3757534.37599.37576.2535.625PE9591.2558.7596.2593.12545.625VCR + EE76.2575.62536.2599.37576.87538.125VCR + EE + PE98.7598.12538.12599.37599.37536.875Through the comparison of the values in Tables4 and 5, it can be seen that the test classification accuracy of 3 hp reaches 99.375%, indicating that the training classification model has certain robustness. Therefore, based on the VMD method to decompose the signal and extract the characteristics, it is more effective in diagnosing rolling bearing failure.


## 4.1.1. Determination of the Number of Modes

Whether the choice of the decomposition numberK is reasonable will have different effects on the decomposition result. Assuming that the number of IMF components K is less than the number of useful components in the processed signal, it will result in insufficient decomposition; on the contrary, it will make the center frequencies of adjacent IMF components close, resulting in frequency aliasing. The rolling element fault signal is used for VMD decomposition. The center frequency under different K values is shown in Table 1.Table 1
The center frequency of components with differentK values.
KCenter frequency (Hz)363528813366458613332884336655571102283832603430655711012838325934284164It can be seen from the table that when the number of components is 5, the center frequencies of 3260 Hz and 3430 Hz are relatively close, and modal aliasing may occur, so it is more appropriate to select 4 as the number of components. Select the first set of data of the inner circle fault, and use the input parametersα = 2000, DC = 0, K = 4, t = 0, tol = 1e−7, w (init) = 1, ε = 1e−7 to perform VMD decomposition.The time-domain diagram and corresponding spectrum diagram after the data is decomposed by VMD are shown in Figures6 and 7.Figure 6
The time-domain diagram after inner ring fault signal decomposition.Figure 7
The frequency-domain diagram after inner ring fault signal decomposition.


## 4.1.2. Feature Extraction

According to formulae (8), (10), and (12), extract the Mi, EE(HEN), PE of each IMF component to construct a 160∗12 feature vector matrix. After normalizing the data, the box plot is shown in Figure 8. It can be seen from the figure that the fault features all exhibit good clustering characteristics, which are all between 0 and 1, and there are no abnormal data points.Figure 8
The box diagram of fault data.


## 4.1.3. Bearing Fault Identification and Classification

The SVM model uses a cross-validation search algorithm to select the gamma functiong and the penalty coefficient c in the kernel function, and select c = 0.0625 and g = 0.0625. There are 160 sets of data samples, of which the normal state number is 1–40, the inner ring fault number is 41–80, the outer ring fault number is 81–120, and the rolling element fault number is 121–160. Meanwhile, define labels on the sample data, respectively; 1 represents the normal state, 2 represents the inner circle failure, 3 represents the outer circle failure, 4 represents rolling element failure. Three-quarters of the data were randomly selected from each state sample and the order was shuffled. The corresponding labels are also scrambled accordingly. The data was input into SVM classifier for training, and the SVM model for fault data classification was obtained. The remaining sample data is used as test data. Use the trained SVM model to classify and predict the disordered test data. The classification results of the tests are shown in Figure 9. It can be seen from the figure that the prediction test set classification results are the same as the actual test set classification. The confusion matrix of the classification results is presented in Figure 10. The columns stand for the actual label of each condition, and the rows stand for the predict label. It shows that the four conditions have the accuracy of 100%.Figure 9
The predicted results with VMD.Figure 10
The confusion matrix of the result with VMD.In order to further verify the effectiveness of the VMD method, the EMD and MEEMD algorithms are used to process the sample data of the four fault states. Meanwhile, single or multiple statistical features are selected as feature parameters and combined with decomposition algorithms. The correlation coefficient method was used to select the IMF components of the first four orders obtained from the decomposition of each sample data and extract the characteristics. The diagnostic process is the same as above.It can be seen from Table2 that, for the combination of VCR, EE, and PE, the accuracy rates of the three have reached 100%. However, for individual characteristics VCR, EE, and PE, VMD has the highest accuracy rate of 100%. For the combination of VCR and EE, VMD also has the highest accuracy at 100%.Table 2
The classification results with SVM (2 hp).
FeaturesAccuracy of classification (%)VMDEMDMEEMDVCR100.0100.097.5EE100.0100.085.0PE100.095.077.5VCR + EE100.095.097.5VCR + EE + PE100.0100.0100.0For comparison, BP neural network is used to classify and recognize the same feature parameters. The learning rate of BP neural network is 0.1, and the number of iterations is 50. The numbers of input layers, hidden layers, and output layers are 4, 5, 4 (VCR, EE, PE), 8, 9, 4 (VCR + EE), and 12, 13, 4 (VCR + EE + PE). The classification and recognition results are shown in Table3.Table 3
The classification results with BP network (2 hp).
FeaturesAccuracy of classification (%)VMDEMDMEEMDVCR100.088.7591.25EE100.098.085.25PE100.094.2577.35VCR + EE100.098.090.5VCR + EE + PE100.0100%97.25Through the comparison of the values in Table2 and 3, for the same feature parameters, the SVM classifier is generally better than the BP neural network classification effect. Meanwhile, test experiments with different motor loads (0 hp and 1 hp) were added, and the classification results are shown in Table 4. The training models of 0 hp and 1 hp were used to test and verify the characteristic parameters of 3 hp, and the classification results are shown in Table 5.Table 4
The classification results of 0 hp and 1 hp (SVM).
FeaturesAccuracy of classification (%)0 hp1 hpVMDEMDMEEMDVMDEMDMEEMDVCR100.097.592.5100.0100.0100.0EE100.097.585.0100.0100.095.0PE97.5100.077.5100.095.087.5VCR + EE100.0100.085.0100.0100.097.5VCR + EE + PE100.0100.0100.0100.0100.0100.0Table 5
The test results of 3 hp under the training model of 0 hp and 1 hp (SVM).
FeaturesClassification accuracy of 3 hp (%)0 hp training model1 hp training modelVMDEMDMEEMDVMDEMDMEEMDVCR71.258043.7594.37576.87535.0EE99.3757534.37599.37576.2535.625PE9591.2558.7596.2593.12545.625VCR + EE76.2575.62536.2599.37576.87538.125VCR + EE + PE98.7598.12538.12599.37599.37536.875Through the comparison of the values in Tables4 and 5, it can be seen that the test classification accuracy of 3 hp reaches 99.375%, indicating that the training classification model has certain robustness. Therefore, based on the VMD method to decompose the signal and extract the characteristics, it is more effective in diagnosing rolling bearing failure.


## 4.2. Experiment 2

The mine hoist is mainly composed of driving motor, spindle device, hydraulic and braking system, electrical control system, and depth indicator. The spindle device, which bears all the dynamic and static loads acting on the spindle device through the wire rope and the transmission torque of the drive motor, is the key rotating body supporting the mine hoist and is mainly composed of the reel, the spindle, the main bearing, the bearing end cover, and so on. The VMD method in the previous section is applied to the mechanical failure simulation experiment platform of the hoist spindle device in Figure11.Figure 11
The fault simulation test bench.The test bench selects the Panasonic model MSMF082 L wire-type motor with rated speed of 3000 r/min and output power of 750 W. The servo motor speed is set at 20 r/min. The model of AC192-1D vibration acceleration sensor produced by American CTC vibration acceleration sensor company is used to record the fault vibration signal of the bearing under different conditions: normal state, inner ring damage, outer ring damage, and rolling element damage. The model of acquisition card is Smacq-USB-4000, and the vibration signal is collected at a sampling frequency of 5000 Hz.In order to simulate the working condition of the lifting load of the hoist, the test bench uses a magnetic powder dynamometer to apply a load to the steel wire rope. The size of the load is adjusted and displayed by the NX6002 speed torque measuring instrument. Under different fault conditions, the torque of the dynamometer is adjusted to 20 N·m. A single point of failure (width 1.5 mm, depth 1.5 mm) was processed on the outer and inner rings of the bearing by electrical discharge machining, and a fine slit (width 1 mm, depth 1 mm) was cut on the roller by wire cutting. The relevant parameters of the bearing are shown in Table6. Figure 12 shows the four failure states of rolling bearings.Table 6
The specifications of the test bearing.
FeatureDetailBearing typeSKF22207CA/W33Roller diameter d (mm)8.2Pitch diameter D (mm)55Number of balls (N)30Contact angleθ°0Figure 12
The object diagram of the roller bearing: (a) outer ring fault; (b) inner ring fault; (c) rolling element fault; (d) normal.
(a)(b)(c)(d)The vibration signal collected by the sensor near the servo motor driver is selected as the analysis signal, and the wavelet packet is used for noise reduction. The 81920 data points were intercepted from the four states of the bearing and divided into 40 groups, with each length 2048. A total of 160 data samples were obtained. Figure13 shows the time-domain diagram of bearing vibration signals in four states.Figure 13
The time-domain diagram of four states of rolling bearing.


## 5. Results and Discussion

The rolling bearing diagnosis process of experiment 2 is referred to experiment 1. Through the center frequency method,K is set to 6. By using VMD, the rolling bearing signals are decomposed using the input parameters α = 2000, DC = 0, K = 6, t = 0, tol = 1e−7, w (init) = 1, e = 1e−7. Taking the data on the 40th line of the rolling element as an example, Figures 14 and 15, respectively, show the time-domain diagram and the spectrum diagram of IMFs.Figure 14
The time-domain diagram after rolling element fault signal decomposition.Figure 15
The frequency-domain diagram after rolling element fault signal decomposition.The IMFs of rolling bearing signals are then used to evaluate features such as VCR, EE, and PE. Single or multiple statistical features are selected as feature parameters and feature sets are constructed. From each state sample data, 75% of the data were randomly selected as the training data.The trained SVM model was used to classify and predict the test data. The prediction accuracy is shown in Table7. When the feature is a combination of VCR, EE, and PE, the classification results of the test are shown in Figure 16. As can be seen from the figure, there are points where the predicted label and the original label do not coincide. The confusion matrix of the classification results is presented in Figure 17. It can be seen that the accuracies of the four conditions achieve 60%, 80%, 60%, and 90%.Table 7
The classification accuracy of the classifier.
FeaturesClassification accuracy (%)VCR47.5EE60PE40VCR + EE70VCR + PE40EE + PE65VCR + EE + PE72.5Figure 16
The comparison of VMD test results.Figure 17
The confusion matrix of the classification results.As can be seen from Table7, the classification accuracy of individual features VCR, EE, and PE is relatively low. For the combination of VCR, EE, and PE, the classification accuracy is up to 72.5% (29/40).


## 6. Conclusions

Bearing damage is an important problem in the system of hoist spindle. In this paper, VMD method is used for fault diagnosis of damaged bearings. Vibration signal acquisition experiments were carried out on rolling bearings damaged in different parts. The VMD algorithm was used to decompose the signal into different IMFs. The contribution rate, energy entropy, and permutation entropy of each IMF were extracted. These important features are used as input parameters to the classifier. For comparison, some feature parameters are used as input parameters alone, and some are used as input parameters in combination. SVM classifier was used for classification.The following conclusions are drawn: (1) the VMD algorithm was applied to the vibration signal analysis of rolling bearing. The center frequency method was used to determine the number of K values. Experiments proved the feasibility of the VMD algorithm for signal decomposition analysis. VMD can be used as an adaptive signal processing tool to reveal the weak fault characteristics of complex vibration signals. (2) Fault feature extraction can be realized by quantifying fault feature information with VCR, EE, and PE, and creating feature vector matrix. Based on SVM, the sample matrix data is effectively diagnosed compared with BP network. (3) The results show that the classification accuracy of SVM classifier is the highest when combining the characteristics of VCR, EE, and PE rather than a single feature.The method proposed in this paper can provide a good practical reference for the fault diagnosis of mechanical equipment of the hoist spindle device.



---

*Source: 8849513-2020-09-23.xml*