1 Introduction

Rolling element bearings are critical components in rotating machines and their states in operation directly affect industrial productivity and workers' safety. Fault diagnosis in the early stage could provide maintenance guidance and avoid huge economic losses so a fast and precise method is necessary. Traditionally, the framework of fault diagnosis for bearings consists of five parts: data acquisition, data preprocessing, feature extraction, dimensionality reduction, and fault classification^[1-3].

Prominent features are critical in fault diagnosis, so multiple domain features have been investigated, including root mean square (RMS), kurtosis, and skewness in time domain and frequency domain, and short-time Fourier transform (STFT)^[4] in time-frequency domain^[5]. To deal with nonlinear dynamic characteristics of bearing fault signals, several entropy-based approaches have been proposed, such as approximate entropy (ApEn)^[6-7], sample entropy (SaEn)^[8-9], fuzzy entropy^[10-11], multiscale entropy (MSE)^[12-13], and permutation entropy (PE)^[14-15]. PE can measure the complexity of time series through comparing adjacent values, which is simple, immune to noise, and suitable for online monitoring. Based on the PE method, multi-scale permutation entropy (MPE) was developed by Aziz and Arif^[16] and was employed to estimate the complexity of the signals in different scales. MPE excels at its stability and robustness and performs better than PE in the application of bearing fault diagnosis^[17-21].

Traditional approaches for dimensionality reduction are principal component analysis (PCA)^[1] and modifications like locally linear embedding (LLE)^[2], Laplacian score (LS)^[3], and non-negative matrix factorization (NMF)^[4]. With the selected features, the fault diagnosis should be performed with a classifier, which consists of two types based on the vibration data, i.e., statistical models^{[1, 3]} and data-driven models^{[2, 11, 13, 15, 20-22]}. The two classifiers have their own advantages for bearing fault classification but are with deficiencies that they are complex and difficult to implement.

Since the essence of fault diagnosis is the problem of pattern recognition, distance would be an essential factor for the discrimination of different feature vectors. Hence, the nearest neighbor principle should be considered and the fault type could be identified with the smallest distance. For example, dynamic time wrapping^[23] has been applied in bearing fault diagnosis^[24], which measures the similarity of vectors with different length, but is with high computation complexity at the same time. In addition, not all the distances (e.g., the Euclidean distance) could differentiate the pattern. Based on the above analysis, MSD^[25] is applied in this paper for bearing fault classification, which considers the shape of feature vectors and is simple to implement because no tuning parameters need to be set up or optimized compared with the abovementioned classifiers.

The rest of this paper is organized as follows: a brief review of the theoretical background concerning MPE and MSD is provided in Section 2; Section 3 illustrates the procedure of MPE-MSD, and two experiments considering different categories are designed and implemented to verify effectiveness of the proposed method in Section 4; finally, conclusions are drawn in Section 5.

2 Theoretical Background 2.1 Multi-scale Permutation Entropy

PE^[14] was firstly used for detecting the dynamic change of time series data and to overcome former entropy method limitations such as the requirement of long datasets and high computational cost. The definition and implement procedure of PE are described as follows.

Considering a time series {x(k), k=1, 2, …, N}, the m dimensional vector at time i can thus be constructed as

$ \begin{array}{l} x_i^m = \{ x(i),x(i + \tau ), \cdots ,x(i + (m - 1)\tau )\} ,\\ \;\;\;\;\;\;\;i = 1,2, \cdots ,N - (m - 1)\tau \end{array} $

(1)

where x_i^m is a re-constructed series, m represents the embedding dimension, and τ is the time delay. Let π_j=(r₀, r₁, …, r_m-1), then x_i^m has a permutation π_j if it satisfies

$ x\left( {i + {r_0}\tau } \right) \le x\left( {i + {r_1}\tau } \right) \le \cdots \le x\left( {i + {r_{m - 1}}\tau } \right) $

(2)

where 0≤r_n≤m-1, and r_n-1 < r_n when

$ x\left( {t + {r_{n - 1}}\tau } \right) = x\left( {t + {r_n}\tau } \right) $

For each permutation π_j, 1≤j≤m!, the relative frequency can be defined as

$ p\left( {{\mathit{\boldsymbol{\pi }}_j}} \right) = \frac{{\# \left\{ {i|i \le N - (m - 1)\tau ,x_i^m\;{\rm{has type }}{\mathit{\boldsymbol{\pi }}_j}} \right\}}}{{N - (m - 1)\tau }} $

(3)

where # represents the number of x_i^m belonging to the type π_j. Then PE with m dimension is defined as

$ {H_{{\rm{PE}}}}(m) = - \sum\limits_j^{m!} p \left( {{\mathit{\boldsymbol{\pi }}_j}} \right)\ln \left[ {p\left( {{\mathit{\boldsymbol{\pi }}_j}} \right)} \right] $

(4)

The MPE algorithm based on PE comprises the following two steps:

1) Considering a time series {x(k), k=1, 2, …, N}, as illustrated in Fig. 1, the coarse-grained time series y_j^(s) with scale factor s can be constructed as

$ y_j^{(s)} = \frac{1}{s}\sum\limits_{i = (j - 1)s + 1}^{js} {{x_i}} \;\;\;1 \le j \le \frac{N}{s} $

(5)

2) The PEs of each coarse-grained time series are calculated based on Eqs. (1)-(5) and then expressed as the function of the scale factor s

$ {\rm{MPE}}(x,s,m,\tau ) = {\rm{PE}}\left( {y_j^{(s)},m,\tau } \right) $

(6)

Four parameters should be specified before using MPE, including N, m, τ, and s. In this paper, the parameters are set as m=4, τ=1, s=20^[21]. The normalized MPE is expressed as

$ \begin{array}{*{20}{c}} {{\rm{NorMPE}}(x,t,m,\tau ) = {\rm{MPE}}(x,t,m,\tau )/}\\ {\sum\limits_{t = 1}^s {\rm{MPE}} (x,t,m,\tau )} \end{array} $

(7)

2.2 Morphology Similarity Distance

Considering vectors as objects, the morphology similarity distance is proposed to estimate similarity. The morphology similarity distance between two n-dimensionality vectors L_i=(l_i1, l_i2, …, l_in) and L_j=(l_j1, l_j2, …, l_jn) is defined as^[25]

$ {D_{{\rm{MSD}}}} = {D_{{\rm{Euclid}}}} \times \left( {2 - {\rm{ASD}}/{D_{{\rm{Manhattan}}}}} \right) $

(8)

where D_Euclid is the Euclid distance, D_Manhattan is the Manhattan distance, ASD is the absolute sum of the differences as follows:

$ \begin{array}{l} {D_{{\rm{MSD}}}}\left( {{L_i},{L_j}} \right) = \sqrt {\sum\limits_{k = 1}^n {{{\left( {{l_{ik}} - {l_{jk}}} \right)}^2}} } \times \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {2 - \frac{{\left| {\sum\limits_{k = 1}^n {\left( {{l_{ik}} - {l_{jk}}} \right)} } \right|}}{{\sum\limits_{k = 1}^n {\left| {{l_{ik}} - {l_{jk}}} \right|} }}} \right) \end{array} $

(9)

The traditional index used for similarity is the Minkowski metric

$ d\left( {j,k} \right) = {\left( {\sum\limits_{i = 1}^n {{{\left| {{x_{ji}} - {x_{ki}}} \right|}^r}} } \right)^{1/r}} $

(10)

The traditional distances Euclidean (r=2) and Manhattan distances (r=1) can be regarded as special cases of Minkowski distance. However, they both neglect the difference between vectors and fail to reflect shape similarities such as translation, compression, and stretch.

For example, there are four vectors as follows:

$ {\mathit{\boldsymbol{x}}_0} = (5,5,5),{\mathit{\boldsymbol{x}}_1} = (4,4,4) $

$ {\mathit{\boldsymbol{x}}_2} = (5,4,3),{\mathit{\boldsymbol{x}}_3} = (5,4,7) $

Let vector x₀ be the reference vector. The vector differences, Manhattan distance, Euclidean distance, and MSD from x₁, x₂, and x₃ to x₀ are computed and listed in Table 1.

3 Procedure of Fault Diagnosis

A procedure based on MPE and MSD is established for the fault classification of rolling element bearings, and the steps are as follows:

1) Collect vibration signals of healthy and different types of fault bearings;

2) Calculate MPE values of the original signals with different fault types and construct normalized vectors as reference samples and testing samples;

3) Assign the number of the reference samples to each fault type, compute the MSD between every testing sample and all reference samples, and then average the MSD values with the specified fault;

4) Based on the nearest neighbor principle, choose the minimum value among the MSD means from all fault types, and the corresponding fault type can be identified.

The flowchart of the proposed method is shown in Fig. 2.

4 Application 4.1 Case 1: Artificially Seeded Damage Bearings

The bearing data were obtained from Bearing Data Centre of Case Western Reserve University^[26] and the bearing experiment system is shown in Fig. 3. During the experiment, the drive end bearing 6205-2RS JEM SKF was investigated and defective bearings were seeded with single point faults using electro-discharge machining. Vibration signals from the accelerometer were placed at the 12 o'clock position at the drive end of the motor housing with the sampling frequency 48 kHz, including normal, ball fault (BF), inner race fault (IRF), and outer race fault (ORF) (at the 6 o'clock position).

To verify the effectiveness of the proposed method, three groups of tests at defect size 0.007 inch and loading 1 HP were investigated, as depicted in Table 2, including different fault types (Group A) with waveforms shown in Fig. 4(a), different defective sizes (Group B) for damage severity, and a combination of fault types and defective sizes (Group C) for complex conditions. In every condition, non-overlapping segments with length N=4 096 were extracted, and MPE values were computed for each segment.

To illustrate specific fault types in the signal, the envelop spectra based on Teager energy operator corresponding to BF, IRF, and ORF are shown in Fig. 4(b). The remarkable amplitudes occurring in the location of the fault frequencies clearly indicated the fault types. The structural parameters of the bearing are listed in Table 3, and the corresponding fault frequencies are computed as follows:

Ball spin frequency (BSF)

$ {f_B} = 0.5{f_r}\left( {1 - {d^2}{{\cos }^2}\alpha /{D^2}} \right)D/d $

(11)

Ball pass frequency on inner race (BPFI)

$ {f_I} = 0.5{f_r}(1 + d\cos \alpha /D)z $

(12)

Ball pass frequency on outer race (BPFO)

$ {f_O} = 0.5{f_r}(1 - d\cos \alpha /D)z $

(13)

where d is the diameter of the rolling element, D is the pitch diameter, α is the contact angle, z is the number of rolling elements, and f_r is the shaft speed.

Then, samples with numbers 10, 20, and 30 were utilized as reference and the rest were used for verification. The classification results are listed in Table 4.

The classification accuracies reached 100%, 100%, and 97.54% in Group A, Group B, and Group C with 10 reference samples respectively, as illustrated in Fig. 5. Moreover, the accuracy was not highly dependent on the number of the reference samples according to the comparison among different sample numbers in Table 4. In Fig. 5(a), for example, when the bearing was normal, the mean MSD between the testing sample and the normal reference sample was the minimum, which was the same in the other three fault conditions. The four types of fault could be clearly identified based on the prominent distance between the minimum mean MSD and the other three MSDs. As illustrated in Fig. 5(b), distances from the normal condition were BF014 < BF007 < BF201 at the first 100 samples, which means the damage severity was non-monotone because of the non-continuous feature, which could not represent the overall tendency. In Fig. 5(c), a more complicated condition was considered, where a ten-sample interval was set to display the result clearly. Compared with the previous work, considering the same four conditions as normal, BF, IRF, and ORF as listed in Table 5, the results in our work are satisfied.

4.2 Case 2: Run-to-Failure Bearings

The data were obtained from Intelligent Maintenance System (IMS), University of Cincinnati. Four bearings were tested at one time on the same shaft in the bearing test rig as shown in Fig. 6. The shaft is 2 000 r/min and a radial load of 6 000 LB is forced on the shaft with a spring mechanism. The data sampling rate is 20 kHz with the data length 20 480 points. Three groups of tests were conducted in the experiment in total and Ref. [28] can provide more details. Table 6 presents four types of fault bearings and the typical waveforms are shown in Fig. 7. One hundred samples were extracted from each bearing fault condition, and 400 samples in whole were generated with length N=20 480. Then MPE values were computed for each sample.

Following the same abovementioned procedure, samples with numbers 10, 20, 30 were utilized as reference, and the rest were used for verification. The classification accuracies were all 100% with 10, 20, 30 samples, and the results are illustrated in Fig. 8. A comparison of previous work considering the same four conditions as normal, BF, IRF, and ORF is listed in Table 7, which shows the advantage of the proposed method.

The whole life of Bearing 1 in Set No. 2 with ORF was employed for the investigation of damage severity. MPE values were computed at every sample time point. The first tenth points were abandoned because they were not stable at the beginning, and the confidence value (CV) was introduced as an index to evaluate the severity. A normalization function combined the Sigmoid function with MSD is introduced as^[31]

$ {\rm{CV}} = 2 - \frac{2}{{1 + {e^{ - {c_0}{\rm{MSD}}}}}} $

(14)

in which the scale parameter

$ {c_0} = - \ln \left( {\frac{{{\rm{C}}{{\rm{V}}_{{\rm{pre}}}}}}{{2 - {\rm{C}}{{\rm{V}}_{{\rm{pre}}}}}}} \right)/{\rm{Mean}}\left( {{\rm{MS}}{{\rm{D}}_{{\rm{normal}}}}} \right) $

(15)

where Mean (MSD_normal) is the average of all MSDs under normal condition, and CV_pre corresponds to the Mean (MSD_normal), which was determined artificially. In this paper, the scale parameter c₀=10.8 and CV_pre=0.99. For comparison, additional features like RMS, kurtosis, and skewness were considered, and the corresponding CV with MPE is illustrated in Fig. 9. The comparison shows that the CV with MPE only was stable and could monotonously represent the damage severity along the time.

5 Conclusions

To fast and accurately diagnose the fault type in bearings, this paper presents a novel bearing fault diagnosis method based on MPE-MSD. From the above research, some conclusions are drawn as follows:

1) MPE could be employed as the representative feature for fault diagnosis and MSD is an efficient and simple approach for classification without the set-up or optimization of the tuning parameters.

2) Two experiments were performed to verify the feasibility of the proposed method. In comparison with previous studies, high accuracies were obtained in fault diagnosis for different conditions without feature selection.

3) A normalization function combining the Sigmoid function with MSD is proposed for continuous damage severity without considering full fault dataset. The degradation path shows stability and monotonicity with MPE, compared with the conditions considering MPE, RMS, kurtosis, and skewness.