Motor fault diagnosis based on composite multi-scale weighted reverse slope entropy and neighborhood preserving embedding

2.1. Reverse slope entropy

Standardization of vibration sequences $\mathbf{X}$ is performed as:

in which ${\mu }_{\mathbf{X}}$ and ${\sigma }_{\mathbf{X}}$ respectively are the mean and the standard deviation of vibration sequences $\mathbf{X}$, and then a new time series $\mathbf{x}=\left\{x\left(i\right),i=1,2,\cdots ,N\right\}$, is obtained according to Eq. (1) to ensure the following slope entropy computation. Then the subsequence is obtained as following:

in which $j=1,2,\dots ,N-\left(m-1\right)\tau$, $m$ represents the embedding dimension and $\tau$ represents the delay time. Consider the difference in adjacent amplitudes $x_i-x_{i-1}$, and then divide it with the number of segments $q=\text{5}$ as shown in Fig. 1, namely.

1) When $x_i>x_{i-1}+\gamma$, the identifier was set as +2.

2) When $x_i>x_{i-1}+\delta$ and $x_i\le x_{i-1}+\gamma$, the angle is below 45° and above a small area near 0°, set the identifier as +1.

3) When $\left|x_i-x_{i-1}\right|\le \delta$, in a small area near 0°, the identifier was set as 0.

4) When $x_i<x_{i-1}-\delta$ and $x_i\ge x_{i-1}-\gamma$, the angle is above –45° and below the small area near 0°, set the identifier as -1.

5) When $x_i<x_{i-1}-\gamma$, the identifier was set as –2.

Fig. 1Slope entropy with five division levels

There are a total of $q^{m-1}$ types of permutation patterns that may appear after the symbolization of $m$ -dimensional vectors. Let $n$ be the number of all pattern types, ${\pi }_l$ be the $l$th pattern and $\mathrm{\Pi }={\left\{{\pi }_l\right\}}_1^{q^{m-1}}$, then the probability of each pattern is as follows:

in which $1_A\left(u\right)$ is the indicator function of set $\mathbf{A}$. When $u\in \mathbf{A}$, $1_A\left(u\right)=1$; when $u\notin \mathbf{A}$, $1_A\left(u\right)=0$. Substitute Eq. (3) into Eq. (4), and then the slope entropy of the vibration sequence $\mathbf{x}$ at the specified embedding dimension $m$ and delay parameter $\tau$ is obtained as:

Considering the volatility of vibration sequences under the same structure, an improved algorithm is proposed as shown in Fig. 2 [17]. In this algorithm, the same structured data is weighted based on its volatility, and Eq. (3) is replaced with:

When $w_j=C$, $\forall j\le N$ and $C\ge 0$, the above equation degenerates to Eq. (3). Substitute the statistical results of each pattern obtained from Eq. (5) into Eq. (4) and the weighted slope entropy is derived.

In addition, Bandt proposed a reverse permutation entropy [18] to measure the distance between useful information and white noise. On this basis, this paper proposes the weighted reverse slope entropy:

Fig. 2Volatility of data with the same structure [17]

2.2. Composite multi-scale weighted reverse slope entropy

Composite coarse-grained the sequence $\left\{x\left(n\right),n=1,2,\dots ,N\right\}$ to obtain the sequence ${\mathbf{y}}_k^{\left(s\right)}=\left(y_{k,j}^{\left(s\right)}\right)$, that is:

where, $s$ is the scale factor, and $s_{\mathrm{m}\mathrm{a}\mathrm{x}}$ is the maximum scale factor. Under scale factor $s$, calculate the weighted reverse slope entropy ${WRSlE}_k^{\left(s\right)}\left(\mathbf{x},m\right)$ of $s$ coarse-grained sequences $y_{k,j}^{\left(s\right)}$ and take the mean to obtain the composite weighted reverse slope entropy under scale $s$:

Combine the composite weighted reverse slope entropy at multiple scales and the composite multi-scale weighted reverse slope entropy as a vector $\mathbf{C}$ is obtained finally.

2.3. Neighbourhood preserving embedding

Assuming the high-dimensional state space $\mathbf{C}=\left\{c_i\in {\mathbb{R}}^D,i=1,2,\cdots ,D\right\}$, where $D=s_{\mathrm{m}\mathrm{a}\mathrm{x}}$ in this paper, the goal of the NPE algorithm is to find the optimal projection matrix $\mathbf{A}=\left\{{\mathbf{a}}_1,{\mathbf{a}}_2,\cdots ,{\mathbf{a}}_d,\right\}\in {\mathbb{R}}^{D\times d}$, $\left(d<D\right)$ and map high-dimensional data to relatively low dimensional feature spaces $\mathbf{Y}=\left\{y_i\in {\mathbb{R}}^d,i=1,2,\cdots ,N\right\}$ through dimensionality reduction [15]. The detailed algorithm is as follows:

(1) Determine the nearest neighbor points. Use the $k$-nearest neighbor method to find the $k$ sample points closest to the sample point $c_i$, in which Euclidean distance is often employed as a measure between the sample point $c_i$ and the $k$ sample points in the original manifold and then the neighborhood points are selected. Each sample point is only reconstructed from the relevant nearest neighbor points, and when $c_j$ is not a neighborhood of $c_i$, the connection weight $w_{ij}=0$.

(2) Reconstruct the weight matrix. Reconstruct of each data point c_i linearly from neighboring points, and calculate of the reconstruction weight matrix $\mathbf{W}$ for this sample point. Minimize the following reconstruction error:

where, ${\sum }_{j=1}^k{w}_{ij}=1$.

(3) Calculate the feature mapping matrix. Maintain the linear reconstruction weight of the sample points unchanged, reconstruct the original data sample, and minimize the reconstruction error, i.e.:

where, $\mathbf{M}={\left(\mathbf{I}-\mathbf{W}\right)}^T\left(\mathbf{I}-\mathbf{W}\right)$, $\mathbf{I}$ is the identity diagonal matrix. Constraints are:

By using the Lagrangian method, the objective function and constraint conditions are merged and derived, and the minimization problem of the objective function can be transformed into solving the generalized eigenvalue problem as shown in Eq. (11):

Since $\mathbf{C}\mathbf{M}{\mathbf{C}}^T$ and $\mathbf{C}{\mathbf{C}}^T$ are semi positive definite matrices, the eigenvectors corresponding to the smallest $d$ eigenvalues ${\lambda }_1\le {\lambda }_2\le \cdots \le {\lambda }_d$ form the projection matrix $\mathbf{A}$:

Then a low dimensional representation of $\mathbf{C}$ in the low dimensional space ${\mathbb{R}}^d$ is obtained, and $\mathbf{Y}$ is the reduced dimensional data matrix.

2.4. Motor fault diagnosis flowchart

The motor fault diagnosis process based on CMWRSlE and NPE is as follows:

1) Collect vibration signals of motor normal and different fault types, and set the data length $N=\text{2048}$ to extract composite multi-scale weighted reverse slope entropy.

2) Normalize each segment of data and then extract CMWRSlE, with parameters set as: embedding dimension $m=$3, maximum scale factor $s_{max}=$30, symbol segmentation parameters $\gamma =$1, $\delta =$0.001.

3) The number of neighboring points in NPE is $k=$5, and the first 5 principal components (PCs) after dimensionality reduction are taken as training features.

4) Set the number of training samples and testing samples respectively for each type of motor fault, substitute the training samples into the HPA for classification training while the testing samples for testing, and define the classification accuracy as:

where, $n_R$ and $n_T$ respectively represents the total number of correctly classified samples and the total number of test samples.

3.1. Motor rotor fault data from Wuhan University in China

The simulation test bench for motor rotor faults at Wuhan University in China is shown in Fig. 3 [4]. The vibration signals of the motor rotor state include four types: normal, contact friction, imbalance, and misalignment. The rotor imbalance is simulated by screwing a 2 g mass block into the threaded hole at the edge of the mass disk; misalignment is achieved by changing the relative position of the shafts at both ends of the coupling; contact friction fault is achieved by screwing the friction screw into the shell of the friction screw and then contacting the rotating shaft. During the signal acquisition process of normal and various faults, the rotor speed is 1200 r/min, the sampling frequency is 2048 Hz, and the sampling time is 1 s. The rotor state data used in this paper are pre-processed with wavelet threshold denoising, as shown in Fig. 4.

Fig. 3Rotor laboratory bench [4]

45 sets of data were collected for each type of rotor state, and a total of 180 sets of data were acquired finally. The mean and the standard deviation (std) of the entropy values at each scale factor are calculated as shown in Fig. 5(a), while the mean and the std of the first three PCs after NPE dimensionality reduction are shown in Fig. 5(b). As the first three PCs shown in Fig. 6, it is indicated that various faults of the motor rotor can be clearly identified.

Randomly select 10 training samples from each state, which are then dimensionally reduced by NPE and brought into the HPA for training. The rest 35 samples are then dimensionally reduced by NPE and used as testing samples. Randomly select 10 training samples and repeat the above test 100 times. The classification results are shown in Table 1, and the classification accuracy of each test is 100 %.

Fig. 4Vibration signals of rotor faults

Fig. 5Statistics of motor rotor fault features

a) The mean with the corresponding std of the entropy values at each scale factor

b) The mean with the corresponding std of the first three PCs

Fig. 6The first three PCs of state features of motor rotor with NPE

3.2. Motor bearing failure data from Case Western Reserve University in the United States

The artificial fault bearing experimental device at Case Western Reserve University (CWRU) consists of a driving motor, a loading motor, a torque sensor, etc. The test bearing is installed at the driving end, and the deep groove ball bearing model is SKF6205, as shown in Fig. 7 [5]. There are three types of single point damage on bearings processed by electric discharge: rolling element failure (BF), inner ring failure (IRF), and outer ring failure (ORF). The diameter of single point damage for each fault type is shown in Table 2, which represents the degree of damage. During the signal acquisition process of normal and various bearing faults, the sampling frequency is 12 kHz and the sampling time is 1 second. Fig. 8 shows the fault signals of each rolling bearing component at the damage diameter 0.007 inches with the sampling length 2048.

Fig. 7Bearing fault test rig of CWRU [5]

Extract 50 samples of fault feature vectors for normal and various fault types of bearings in Table 2, and a total of 500 samples are obtained. The mean and the standard deviation (std) of the entropy values at each scale factor are calculated as shown in Fig. 9(a), while the mean and the std of the first three PCs of various fault characteristics of motor bearings in Table 2 after NPE dimensionality reduction are shown in Fig. 9(b). As the first three PCs shown in Fig. 10(a), it can be seen that there is a slight overlap between BF21 and ORF14, but the two can be clearly separated in the fourth to sixth PCs as shown in Fig. 10(b).

10 samples of training samples were randomly selected for the normal state and each type of fault of the motor bearings. After NPE dimensionality reduction, they were brought into the HPA for training, while the rest 40 samples were used as testing samples after NPE dimensionality reduction. Randomly select 10 training samples and repeat the above test 20 times, and the classification accuracy of each test is 100 %.

Fig. 8Vibration signals at the damage level 0.007 inch for different fault types of bearings on the condition of load 2 hp and rotating speed 1750 r/min

Fig. 9Statistics of motor bearing fault features

a) The mean with the corresponding std of the entropy values at each scale factor

b) The mean with the corresponding std of the first three PCs

Fig. 10The first six PCs of fault features of motor bearings with NPE

a) The first three PCs

b) The fourth to sixth PCs