# Fault Diagnosis Approach for Rotating Machinery Based on Feature Importance Ranking and Selection

Authors: Zong Yuan; Taotao Zhou; Jie Liu; Changhe Zhang; Yong Liu
Journal: Shock and Vibration
 (2021)

---

## Abstract

The key to fault diagnosis of rotating machinery is to extract fault features effectively and select the appropriate classification algorithm. As a common signal decomposition method, the effect of wavelet packet decomposition (WPD) largely depends on the applicability of the wavelet basis function (WBF). In this paper, a novel fault diagnosis approach for rotating machinery based on feature importance ranking and selection is proposed. Firstly, a two-step principle is proposed to select the most suitable WBF for the vibration signal, based on which an optimized WPD (OWPD) method is proposed to decompose the vibration signal and extract the fault information in the frequency domain. Secondly, FE is utilized to extract fault features of the decomposed subsignals of OWPD. Thirdly, the categorical boosting (CatBoost) algorithm is introduced to rank the fault features by a certain strategy, and the optimal feature set is further utilized to identify and diagnose the fault types. A hybrid dataset of bearing and rotor faults and an actual dataset of the one-stage reduction gearbox are utilized for experimental verification. Experimental results indicate that the proposed approach can achieve higher fault diagnosis accuracy using fewer features under complex working conditions.


---

## Body

## 1. Introduction

At present, the components of industrial rotating machinery equipment are becoming increasingly complex and compact. As a key component of the transmission system, rotating machinery is an important part of modern industrial machinery equipment, including motors, engines, bearings, and gearboxes [1–3]. When rotating machinery is operating under harsh or complex conditions, its key components are extremely prone to failure, which may cause the shutdown of the entire mechanical equipment, and even endanger the safety of surrounding operators [4, 5]. Therefore, it is significant to construct a fault diagnosis scheme for rotating machinery under complex conditions to accurately detect and diagnose its health or fault state.Fault diagnosis based on vibration signal analysis is the main research hotspot at present, among which the most critical step is feature extraction [6, 7]. Based on vibration signal analysis, the existing methods mainly extract fault features from the time domain, frequency domain, and time-frequency domain [8, 9]. Time-domain features contain root mean square, mean, standard deviation, kurtosis, etc., which may be valid only for certain fault types [10, 11]. The frequency-domain analysis is mostly based on the Fourier transform (FT) [12]. However, these methods are limited by prior knowledge and experience in practical applications due to the nonlinearity and nonstationarity of the original vibration signal, which makes it difficult for them to effectively mine the fault information hidden in the vibration signal [11, 13].As a measure of time-domain uncertainty, entropy-based analysis methods have attracted extensive attention of scholars, which have been widely used in image processing, biological analysis, and other fields [14]. These entropies mainly include sample entropy (SE) [15], approximate entropy (AE) [16], permutation entropy (PE) [17], fuzzy entropy (FE) [18], and symbolic dynamic entropy (SDE) [19]. FE uses a Gaussian function to replace the Heaviside function in SE to measure the similarity between two vectors. In recent studies, these entropies are usually combined with some other processing strategies, e.g., multiscale and improved multiscale techniques, and some signal decomposition methods in the time-frequency domain [19–21]. However, the multiscale analysis may result in the absence of part of the frequency band components, which may lead to the loss of some fault information. In addition, traditional signal decomposition methods also have some disadvantages that cannot be ignored.Research studies based on time-frequency domain analysis have been done for a long time. For example, empirical mode decomposition (EMD) and local mean decomposition (LMD) both are self-adaptive time-frequency decomposition methods [22, 23]. LMD is improved on the basis of EMD to better maintain the local characteristics of the original signal [24]. Variational mode decomposition (VMD) proposed by Dragomiretskiy et al. is a self-adaptive decomposition method that aims to overcome the shortcomings of undershoot, overshoot, mode mixing in EMD [25–27]. However, there are still some defects in VMD, e.g., the massive consumption of computing resources [27, 28]. In addition, these methods are all based on “mode,” and some frequency information will be lost, which means that they do not apply to frequency analysis.Compared with the above methods, the decomposition strategy of wavelet packet decomposition (WPD) is to pass the signal through a series of filters with different central frequencies but the same bandwidth. Therefore, the signal analysis performed by WPD is more refined, especially for the high-frequency components [29]. In [30], an automatic method combining WPD and EMD was proposed to detect the weak defects of rolling bearings. In [31], WPD was combined with PE to extract fault features of rolling bearings. However, how to select the optimal wavelet basis function (WBF) has not been analysed and discussed in these references. The WBF is a group of functions obtained from the expansion and translation, including db wavelets, sym wavelets, and mexh wavelets [32]. Different WBFs are applicable to different analysis objects, and improper selection may affect the accuracy of fault pattern recognition [32, 33]. In this paper, a two-step principle is proposed to select the most suitable WBF for the fault vibration signal. Then, the optimized WPD (OWPD) method is proposed and applied to decompose the vibration signal to obtain its frequency component. In view of the advantage that entropy measure can effectively extract dynamic information of time series, FE is used to extract hidden fault features from decomposed subsignals. Meanwhile, it also has the advantage of being insensitive to background noise and good robustness [18, 20]. Therefore, a novel fault feature exaction method combining OWPD and FE is further proposed in this paper.After feature extraction utilizing OWPD and FE, a classification algorithm with good performance and computational efficiency is needed to give final diagnosis results. In addition, screening redundant features before fault classification can effectively reduce feature dimension and computational burden and further improve classification accuracy [34]. Traditional classification algorithms include support vector machine (SVM) [35], K-nearest neighbor (KNN) [36], artificial neural networks (ANNs) [37], and random forest (RF) [38]. Deep learning (DL) algorithms include convolutional neural network (CNN) [39], autoencoders (AEs) [40], and deep belief network (DBN) [11, 41]. However, these classification algorithms still have some inevitable shortcomings. For example, SVM is not effective for large-scale training samples and sensitive to the selection of parameters and kernel function. RF is easy to overfit in noisy classification or regression problems. The structure and parameters of some DL algorithms, e.g., DBN, are basically determined by human experience, which not only affects the accuracy of diagnostic results but also causes a large amount of computing costs [11, 42]. CatBoost is a new implementation of the gradient boosting decision tree (GBDT) framework [43]. It has the advantages of high efficiency, few parameters, and strong generalization ability and has excellent performance in many machine learning tasks [43–45]. In addition, as an algorithm based on the decision tree, it can obtain the importance of each feature according to the tree model after gradient boosting, and then the valuable features can be effectively selected for model training. Therefore, it is introduced for feature selection to form a feature set that contains the main fault information. To the best of authors’ knowledge, CatBoost algorithm is rarely studied in the field of fault diagnosis of rotating machinery. In this paper, it is introduced not only for fault pattern recognition but also for selecting the optimal features.Finally, the optimization of hyperparameters is also an urgent problem to be solved in the use of the CatBoost algorithm, which usually has a great impact on the performance of the model. The optimization of hyperparameters is to find an acceptable solution for the optimization goal as effectively as possible [46]. Due to the large amount of data and large solution space, the application of traditional solution methods, e.g., grid search and greedy algorithm, has been limited, while intelligent algorithms such as differential evolution have been widely used due to their fast computing efficiency and the ability to obtain global optimal solution [46, 47]. In this paper, Bayesian optimization (BO) algorithm [48] is considered to solve this problem to find the optimal hyperparameters of the CatBoost classifier. It can obtain the global optimal solution through Gaussian process, which has the advantages of high search efficiency and less iteration times, and can be used for the optimization of any black-box function.Based on the above analysis, aiming to solve the defect that traditional feature extraction methods cannot fully explore the deep-level fault features and to improve the performance of fault pattern recognition, a novel fault diagnosis approach based on feature importance ranking and selection is proposed. In summary, the advantages of WPD in signal decomposition, FE in feature extraction, and CatBoost in fault pattern recognition are fully exploited in the proposed approach. The main contributions can be summarized as follows:(1)
A two-step principle is proposed to select the optimal WBF adaptively according to the characteristics of the mechanical vibration signal.(2)
A fault feature extraction method combining OWPD and FE is proposed, where OWPD is utilized to decompose the vibration signal, and FE is further adopted to form the fault feature set.(3)
CatBoost algorithm is introduced not only for fault pattern recognition but also for feature selection to filter redundant features, which helps to reduce model training time and improve the classification accuracy.(4)
BO algorithm is adopted to solve the optimization problem of hyperparameters in CatBoost. On this basis, the BO-CatBoost algorithm is established and applied to the fault diagnosis of rotating machinery.The remainder of this paper is organized as follows. Section2 introduces the theoretical knowledge and methods of the proposed approach. The diagnosis process and the preliminary validation of the proposed fault diagnosis approach with a mechanical fault simulation (MFS) platform dataset are detailed in Section 3. Further experimental verification using another actual dataset of the one-stage reduction gearbox is shown in Section 4. Section 5 contains the conclusions and future research studies.


## 2. Materials and Methods

### 2.1. Optimized Wavelet Packet Decomposition

#### 2.1.1. Wavelet Packet Decomposition

Generally, FT has been widely used in traditional vibration signal analysis [49]. However, only the frequency-domain information is retained in FT, while the time-domain information is completely lost, which makes it unsuitable for the analysis of nonstationary time-varying signals. Wavelet transform (WT) can provide information in both frequency and time domains to overcome the deficiency of FT [50]. However, only the low-frequency coefficients will be decomposed again in the WT method, which will cause the problem of missing high-frequency information. WPD was proposed to address this deficiency, where the information in both low-frequency band and high-frequency band is completely preserved [51]. The schematic diagram of a three-layer WPD is shown in Figure 1, and the theory is described as follows.Figure 1
Schematic diagram of the three-layer wavelet packet decomposition.Letj denote the decomposition layer, n denote the frequency factor n=0,1,2,…,2j−1, and ψ(n) and ϕ(n) represent the wavelet function and scale function, respectively. Given φ0(n) = ϕ(n) and φ1(n) = ψ(n), the wavelet packet φi(n) (i = 0, 1, 2, …) can be defined as(1)φ2in=2∑k∈Zhkφi2n−k,φ2i+1n=2∑k∈Zgkφi2n−k,where k is the shift factor, Z is the integer set, hk denotes the low-pass filter, gk denotes the high-pass filter, and hk, gk is a couple of quadruple mirror filters that satisfies gk=−1khk.For a given time seriesx(n), let xjpnp=0,1,2,…,2j−1 denote its subsignal, which can be represented as a linear combination of the corresponding wavelet function of the wavelet packet:(2)xjpn=∑k∈Zdjpkφ2j−1+p2n−k,where djpk denotes the p-th wavelet packet coefficient of the j-th layer, and it can be obtained by the inner product between x(n) and φ2j−1+p2n−k, namely,(3)djpk≤xn,φ2j−1+p2n−k≥∫−∞+∞xnφ2j−1+p2n−kdn.An approximationxj(n) of x(n) with layer j equates the sum of all the subsignals:(4)xjn=∑pxjpn.


#### 2.1.2. Optimized Wavelet Packet Decomposition

The main idea of the OWPD method is to automatically select the most suitable WBF for vibration signal analysis of rotating machinery according to the proposed two-step principle on the basis of WPD. In general, the main characteristics to be considered in selecting the WBF include orthogonality, compactness, symmetry, and vanishing moment. Considering the characteristics of different wavelet families and vibration signals of rotating machinery, the coif wavelets, db wavelets, and sym wavelets are selected as the candidate WBFs. Figure2 shows the proposed two-step principle for selecting the optimal WBF, detailed as follows.Step 1: select the candidate WBFs preliminarily from the same wavelet family according to the principle of maximum energy-to-Shannon entropy ratio (METSE) [52].(1)
Calculate the energy valueE(n) of the n-th node:(5)En=∑i=1mCn,i2,wherei and m are the serial number and total number of discrete points in the n-th node, respectively, and Cn,i is the coefficient corresponding to the discrete point.(2)
The Shannon entropy of then-th node is defined as(6)Sentropyn=−∑i=1mpilog2pi,wherepi is the energy probability distribution of the wavelet coefficients, defined as(7)pi=Cn,i2En.(3)
The ratio of the total energy and the total Shannon entropy of thej-layer WPD is represented as ζ, namely,(8)ζ=∑n=12jEn∑n=12jSentropyn.According to equations (5)–(8), the candidate WBF with the largest ζ value can be selected from the same wavelet family.Step 2: select the optimal WBF further from different wavelet families according to the principle of similarity measure.Figure 2
The proposed two-step principle of selecting the optimal WBF.Firstly, the candidate WBFs from different wavelet families mentioned above are applied to implement WPD, respectively. Then, the signal is reconstructed using the coefficients with the nodes of the last layer. Finally, a standardized Euclidean distance is used to measure the similarity between the original signalxi and reconstructed signal yii=1,2,…,N:(9)d=∑i=1Nxi−yisi2,where si is the standard deviation between xi and yi. The smaller the value of d is, the closer the reconstructed signal is to the original signal, and the corresponding WBF is more suitable for signal analysis.


### 2.2. Fuzzy Entropy

Given anN-dimensional time series μ1,μ2,…,μN, the phase space dimension and the similarity tolerance are defined as m (m ≦ N − 2) and r, respectively. Then, the phase space can be reconstructed as(10)Xi=μi,μi+1,…,μi+m−1−μ0i,i=1,2,…,N−m+1,where(11)μ0i=1m∑j=0m−1μi+j.The fuzzy membership functionA(x) is introduced as(12)Ax=1,x=0,exp−ln2xr2,x>0.Fori=1,2,…,N−m+1, calculate(13)Aijm=exp−ln2⋅dijmr2,j=1,2,…,N−m+1,j≠i,where dijm is the maximum absolute distance between window vector X(i) and X(j), that is,(14)dijm=dXi,Xj=maxp=1,2,…,mμi+p−1−μ0i−μj+p−1−μ0j.The function Φm is defined as(15)Φmr=1N−m+1∑i=1N−m+11N−m∑j=1,j≠iN−m+1Aijm.Therefore, the FE value of the original time series can be calculated as(16)FEm,n,r,N=lnΦmn,r−lnΦm+1n,r.


### 2.3. Categorical Boosting

#### 2.3.1. CatBoost for Classification

In the GBDT algorithm, lots of decision trees are combined to produce a model of high accuracy, and the progress can be written as(17)yx=∑t=1Tftx,θt,where x is the feature vector, T is the number of trees, θtt=1,2,…,T is a learned parameter, and ft (x, θt) represents the decision trees that are learned.Given training samplesD=xk,yk1n, where n is the number of samples, xkk=1,2,…,n is the sample data, and yk is the true sample label. To learn the model introduced in equation (17), the following objective function needs to be minimized:(18)Οft=∑i=1nLyk,y^k+∑t=1TΩft,where y^k is the predicted sample label, L is the loss function that represents the difference between yk and y^k, and Ω is the regular function that is used to punish the complexity of ft, defined as(19)Ω=αq+12βω2,where α is the penalty parameter that controls the number of leaf nodes, q is the number of leaf nodes, β is the regularization parameter, and ω is the weight coefficient.Letg denote the negative gradient of the loss function, and the objective function is minimized in the direction of g, namely,(20)g=−∂Lyk,y^k∂y^k.Traditional GBDT algorithms generally have the problem of prediction offset, which affects the generalization ability of the model. To overcome this defect, CatBoost was proposed with two notable improvements [43]: (1) the ordered boosting strategy was adopted to obtain the unbiased estimation of the gradient and slow down the prediction offset; (2) the oblivious tree was used as the basic learner to increase the reliability of the model and speed up the prediction. In addition, to better deal with categorical features, the greedy target-based statistics strategy was improved by adding prior terms in CatBoost algorithm, which can be summarized as three main steps: (1) all the sample datasets are randomly arranged; (2) samples with the same category are selected, and the average label of similar samples is calculated; and (3) features of each sample are digitized by adding the prior term and its corresponding weight coefficient. The improved greedy target-based statistics strategy can be expressed as(21)x^ki=∑j=1nxji=xki⋅yi+a⋅P∑j=1nxji=xki+a,where xki represents the i-th category feature of the k-th sample, x^ki represents the corresponding numerical feature, P represents the increased prior value, and a represents the weight coefficient (a > 0). The addition of prior values can effectively reduce the noise caused by low-frequency features and avoid the overfitting phenomenon.


#### 2.3.2. CatBoost for Feature Selection

The growth strategies of decision trees are different in different GBDT algorithms. XGBoost uses level-wise strategy, which has the disadvantage of inefficiency [53]. CatBoost uses the symmetric tree strategy to optimize the computation of the leaf value to prevent the model from overfitting. In the case of the basic learner of CatBoost is the tree model, the feature coefficient or importance can be obtained according to a certain evaluation index after training the model, e.g., the change of loss function or prediction values.For a given feature setF=f1,f2,…,fN, the feature importance of fii=1,2,…N in the trained CatBoost model is calculated by(22)featurefi=∑trees,leafSfiν1−avr2⋅c1+ν2−avr2⋅c2,and(23)avr=ν1⋅c1+ν2⋅c2c1+c2,where S denotes the different paths to the leaf nodes in the decision tree, c1 and c2 denote the total weight coefficient in the left and right leaves, respectively, and υ1 and υ2 denote the formula value in the left and right leaves, respectively.


#### 2.3.3. BO Algorithm for Optimizing Hyperparameters

In machine learning algorithms, the hyperparameters usually have a strong influence on the performance of the model. There are some inevitable shortcomings in traditional methods of parameter tuning. For example, greedy algorithm can only obtain the local optimal solution, and the uncertainty and nonconvexity of grid search tend to miss global optimality. Different with these methods, BO algorithm can obtain the global optimal solution through Gaussian process, which is considered to find the optimal hyperparameters for the CatBoost model.The basic thought of BO algorithm is that, for the given data and optimal termination condition (usually, the number of iterations or the expected value of the objective function), Bayesian theory is used to estimate the posterior distribution. The distribution and the information from the previous sampling point are used to select the hyperparameters of the later sampling until that the value of the objective function reaches the maximum globally. Here, the objective function is defined as the maximum of the classification accuracy:(24)OFμ=argmax1K∑k=1KfCB,μ,Dt,Dv,where CB denotes the CatBoost classifier, μ=μ1,μ2,…,μn is the hyperparameters, Dt and Dv represent the training and validation set divided by the K-fold cross-validation, respectively, and f (CB, μ, Dt, Dv) is the classification accuracy.


#### 2.3.4. The Proposed Fault Diagnosis Approach

An overview of the proposed fault diagnosis approach for rotating machinery based on feature importance ranking and selection is shown in Figure3. The specific steps are as follows:Step 1: the original vibration signals are acquired by accelerometers and the data acquisition system.Step 2: the optimal WBF is selected according to the proposed two-step principle, and then the original vibration signals are decomposed by OWPD.Step 3: the FE values of the decomposed subsignals are calculated to form the fault feature setF.Step 4: CatBoost algorithm is utilized to obtain the importance of each feature inF by a certain strategy. According to the ranking result of feature importance, the candidate features are selected in sequence and combined with the corresponding labels to form dataset S.Step 5: datasetS is divided into two parts according to a certain proportion: training set and test set.Step 6: the training set is used to train the CatBoost classifier, and BO algorithm is adopted to optimize the main hyperparameters.Step 7: the test set is fed into the trained CatBoost classifier to output the diagnostic results.Figure 3
The flowchart of the proposed approach.


## 2.1. Optimized Wavelet Packet Decomposition

### 2.1.1. Wavelet Packet Decomposition

Generally, FT has been widely used in traditional vibration signal analysis [49]. However, only the frequency-domain information is retained in FT, while the time-domain information is completely lost, which makes it unsuitable for the analysis of nonstationary time-varying signals. Wavelet transform (WT) can provide information in both frequency and time domains to overcome the deficiency of FT [50]. However, only the low-frequency coefficients will be decomposed again in the WT method, which will cause the problem of missing high-frequency information. WPD was proposed to address this deficiency, where the information in both low-frequency band and high-frequency band is completely preserved [51]. The schematic diagram of a three-layer WPD is shown in Figure 1, and the theory is described as follows.Figure 1
Schematic diagram of the three-layer wavelet packet decomposition.Letj denote the decomposition layer, n denote the frequency factor n=0,1,2,…,2j−1, and ψ(n) and ϕ(n) represent the wavelet function and scale function, respectively. Given φ0(n) = ϕ(n) and φ1(n) = ψ(n), the wavelet packet φi(n) (i = 0, 1, 2, …) can be defined as(1)φ2in=2∑k∈Zhkφi2n−k,φ2i+1n=2∑k∈Zgkφi2n−k,where k is the shift factor, Z is the integer set, hk denotes the low-pass filter, gk denotes the high-pass filter, and hk, gk is a couple of quadruple mirror filters that satisfies gk=−1khk.For a given time seriesx(n), let xjpnp=0,1,2,…,2j−1 denote its subsignal, which can be represented as a linear combination of the corresponding wavelet function of the wavelet packet:(2)xjpn=∑k∈Zdjpkφ2j−1+p2n−k,where djpk denotes the p-th wavelet packet coefficient of the j-th layer, and it can be obtained by the inner product between x(n) and φ2j−1+p2n−k, namely,(3)djpk≤xn,φ2j−1+p2n−k≥∫−∞+∞xnφ2j−1+p2n−kdn.An approximationxj(n) of x(n) with layer j equates the sum of all the subsignals:(4)xjn=∑pxjpn.


### 2.1.2. Optimized Wavelet Packet Decomposition

The main idea of the OWPD method is to automatically select the most suitable WBF for vibration signal analysis of rotating machinery according to the proposed two-step principle on the basis of WPD. In general, the main characteristics to be considered in selecting the WBF include orthogonality, compactness, symmetry, and vanishing moment. Considering the characteristics of different wavelet families and vibration signals of rotating machinery, the coif wavelets, db wavelets, and sym wavelets are selected as the candidate WBFs. Figure2 shows the proposed two-step principle for selecting the optimal WBF, detailed as follows.Step 1: select the candidate WBFs preliminarily from the same wavelet family according to the principle of maximum energy-to-Shannon entropy ratio (METSE) [52].(1)
Calculate the energy valueE(n) of the n-th node:(5)En=∑i=1mCn,i2,wherei and m are the serial number and total number of discrete points in the n-th node, respectively, and Cn,i is the coefficient corresponding to the discrete point.(2)
The Shannon entropy of then-th node is defined as(6)Sentropyn=−∑i=1mpilog2pi,wherepi is the energy probability distribution of the wavelet coefficients, defined as(7)pi=Cn,i2En.(3)
The ratio of the total energy and the total Shannon entropy of thej-layer WPD is represented as ζ, namely,(8)ζ=∑n=12jEn∑n=12jSentropyn.According to equations (5)–(8), the candidate WBF with the largest ζ value can be selected from the same wavelet family.Step 2: select the optimal WBF further from different wavelet families according to the principle of similarity measure.Figure 2
The proposed two-step principle of selecting the optimal WBF.Firstly, the candidate WBFs from different wavelet families mentioned above are applied to implement WPD, respectively. Then, the signal is reconstructed using the coefficients with the nodes of the last layer. Finally, a standardized Euclidean distance is used to measure the similarity between the original signalxi and reconstructed signal yii=1,2,…,N:(9)d=∑i=1Nxi−yisi2,where si is the standard deviation between xi and yi. The smaller the value of d is, the closer the reconstructed signal is to the original signal, and the corresponding WBF is more suitable for signal analysis.


## 2.1.1. Wavelet Packet Decomposition

Generally, FT has been widely used in traditional vibration signal analysis [49]. However, only the frequency-domain information is retained in FT, while the time-domain information is completely lost, which makes it unsuitable for the analysis of nonstationary time-varying signals. Wavelet transform (WT) can provide information in both frequency and time domains to overcome the deficiency of FT [50]. However, only the low-frequency coefficients will be decomposed again in the WT method, which will cause the problem of missing high-frequency information. WPD was proposed to address this deficiency, where the information in both low-frequency band and high-frequency band is completely preserved [51]. The schematic diagram of a three-layer WPD is shown in Figure 1, and the theory is described as follows.Figure 1
Schematic diagram of the three-layer wavelet packet decomposition.Letj denote the decomposition layer, n denote the frequency factor n=0,1,2,…,2j−1, and ψ(n) and ϕ(n) represent the wavelet function and scale function, respectively. Given φ0(n) = ϕ(n) and φ1(n) = ψ(n), the wavelet packet φi(n) (i = 0, 1, 2, …) can be defined as(1)φ2in=2∑k∈Zhkφi2n−k,φ2i+1n=2∑k∈Zgkφi2n−k,where k is the shift factor, Z is the integer set, hk denotes the low-pass filter, gk denotes the high-pass filter, and hk, gk is a couple of quadruple mirror filters that satisfies gk=−1khk.For a given time seriesx(n), let xjpnp=0,1,2,…,2j−1 denote its subsignal, which can be represented as a linear combination of the corresponding wavelet function of the wavelet packet:(2)xjpn=∑k∈Zdjpkφ2j−1+p2n−k,where djpk denotes the p-th wavelet packet coefficient of the j-th layer, and it can be obtained by the inner product between x(n) and φ2j−1+p2n−k, namely,(3)djpk≤xn,φ2j−1+p2n−k≥∫−∞+∞xnφ2j−1+p2n−kdn.An approximationxj(n) of x(n) with layer j equates the sum of all the subsignals:(4)xjn=∑pxjpn.


## 2.1.2. Optimized Wavelet Packet Decomposition

The main idea of the OWPD method is to automatically select the most suitable WBF for vibration signal analysis of rotating machinery according to the proposed two-step principle on the basis of WPD. In general, the main characteristics to be considered in selecting the WBF include orthogonality, compactness, symmetry, and vanishing moment. Considering the characteristics of different wavelet families and vibration signals of rotating machinery, the coif wavelets, db wavelets, and sym wavelets are selected as the candidate WBFs. Figure2 shows the proposed two-step principle for selecting the optimal WBF, detailed as follows.Step 1: select the candidate WBFs preliminarily from the same wavelet family according to the principle of maximum energy-to-Shannon entropy ratio (METSE) [52].(1)
Calculate the energy valueE(n) of the n-th node:(5)En=∑i=1mCn,i2,wherei and m are the serial number and total number of discrete points in the n-th node, respectively, and Cn,i is the coefficient corresponding to the discrete point.(2)
The Shannon entropy of then-th node is defined as(6)Sentropyn=−∑i=1mpilog2pi,wherepi is the energy probability distribution of the wavelet coefficients, defined as(7)pi=Cn,i2En.(3)
The ratio of the total energy and the total Shannon entropy of thej-layer WPD is represented as ζ, namely,(8)ζ=∑n=12jEn∑n=12jSentropyn.According to equations (5)–(8), the candidate WBF with the largest ζ value can be selected from the same wavelet family.Step 2: select the optimal WBF further from different wavelet families according to the principle of similarity measure.Figure 2
The proposed two-step principle of selecting the optimal WBF.Firstly, the candidate WBFs from different wavelet families mentioned above are applied to implement WPD, respectively. Then, the signal is reconstructed using the coefficients with the nodes of the last layer. Finally, a standardized Euclidean distance is used to measure the similarity between the original signalxi and reconstructed signal yii=1,2,…,N:(9)d=∑i=1Nxi−yisi2,where si is the standard deviation between xi and yi. The smaller the value of d is, the closer the reconstructed signal is to the original signal, and the corresponding WBF is more suitable for signal analysis.


## 2.2. Fuzzy Entropy

Given anN-dimensional time series μ1,μ2,…,μN, the phase space dimension and the similarity tolerance are defined as m (m ≦ N − 2) and r, respectively. Then, the phase space can be reconstructed as(10)Xi=μi,μi+1,…,μi+m−1−μ0i,i=1,2,…,N−m+1,where(11)μ0i=1m∑j=0m−1μi+j.The fuzzy membership functionA(x) is introduced as(12)Ax=1,x=0,exp−ln2xr2,x>0.Fori=1,2,…,N−m+1, calculate(13)Aijm=exp−ln2⋅dijmr2,j=1,2,…,N−m+1,j≠i,where dijm is the maximum absolute distance between window vector X(i) and X(j), that is,(14)dijm=dXi,Xj=maxp=1,2,…,mμi+p−1−μ0i−μj+p−1−μ0j.The function Φm is defined as(15)Φmr=1N−m+1∑i=1N−m+11N−m∑j=1,j≠iN−m+1Aijm.Therefore, the FE value of the original time series can be calculated as(16)FEm,n,r,N=lnΦmn,r−lnΦm+1n,r.


## 2.3. Categorical Boosting

### 2.3.1. CatBoost for Classification

In the GBDT algorithm, lots of decision trees are combined to produce a model of high accuracy, and the progress can be written as(17)yx=∑t=1Tftx,θt,where x is the feature vector, T is the number of trees, θtt=1,2,…,T is a learned parameter, and ft (x, θt) represents the decision trees that are learned.Given training samplesD=xk,yk1n, where n is the number of samples, xkk=1,2,…,n is the sample data, and yk is the true sample label. To learn the model introduced in equation (17), the following objective function needs to be minimized:(18)Οft=∑i=1nLyk,y^k+∑t=1TΩft,where y^k is the predicted sample label, L is the loss function that represents the difference between yk and y^k, and Ω is the regular function that is used to punish the complexity of ft, defined as(19)Ω=αq+12βω2,where α is the penalty parameter that controls the number of leaf nodes, q is the number of leaf nodes, β is the regularization parameter, and ω is the weight coefficient.Letg denote the negative gradient of the loss function, and the objective function is minimized in the direction of g, namely,(20)g=−∂Lyk,y^k∂y^k.Traditional GBDT algorithms generally have the problem of prediction offset, which affects the generalization ability of the model. To overcome this defect, CatBoost was proposed with two notable improvements [43]: (1) the ordered boosting strategy was adopted to obtain the unbiased estimation of the gradient and slow down the prediction offset; (2) the oblivious tree was used as the basic learner to increase the reliability of the model and speed up the prediction. In addition, to better deal with categorical features, the greedy target-based statistics strategy was improved by adding prior terms in CatBoost algorithm, which can be summarized as three main steps: (1) all the sample datasets are randomly arranged; (2) samples with the same category are selected, and the average label of similar samples is calculated; and (3) features of each sample are digitized by adding the prior term and its corresponding weight coefficient. The improved greedy target-based statistics strategy can be expressed as(21)x^ki=∑j=1nxji=xki⋅yi+a⋅P∑j=1nxji=xki+a,where xki represents the i-th category feature of the k-th sample, x^ki represents the corresponding numerical feature, P represents the increased prior value, and a represents the weight coefficient (a > 0). The addition of prior values can effectively reduce the noise caused by low-frequency features and avoid the overfitting phenomenon.


### 2.3.2. CatBoost for Feature Selection

The growth strategies of decision trees are different in different GBDT algorithms. XGBoost uses level-wise strategy, which has the disadvantage of inefficiency [53]. CatBoost uses the symmetric tree strategy to optimize the computation of the leaf value to prevent the model from overfitting. In the case of the basic learner of CatBoost is the tree model, the feature coefficient or importance can be obtained according to a certain evaluation index after training the model, e.g., the change of loss function or prediction values.For a given feature setF=f1,f2,…,fN, the feature importance of fii=1,2,…N in the trained CatBoost model is calculated by(22)featurefi=∑trees,leafSfiν1−avr2⋅c1+ν2−avr2⋅c2,and(23)avr=ν1⋅c1+ν2⋅c2c1+c2,where S denotes the different paths to the leaf nodes in the decision tree, c1 and c2 denote the total weight coefficient in the left and right leaves, respectively, and υ1 and υ2 denote the formula value in the left and right leaves, respectively.


### 2.3.3. BO Algorithm for Optimizing Hyperparameters

In machine learning algorithms, the hyperparameters usually have a strong influence on the performance of the model. There are some inevitable shortcomings in traditional methods of parameter tuning. For example, greedy algorithm can only obtain the local optimal solution, and the uncertainty and nonconvexity of grid search tend to miss global optimality. Different with these methods, BO algorithm can obtain the global optimal solution through Gaussian process, which is considered to find the optimal hyperparameters for the CatBoost model.The basic thought of BO algorithm is that, for the given data and optimal termination condition (usually, the number of iterations or the expected value of the objective function), Bayesian theory is used to estimate the posterior distribution. The distribution and the information from the previous sampling point are used to select the hyperparameters of the later sampling until that the value of the objective function reaches the maximum globally. Here, the objective function is defined as the maximum of the classification accuracy:(24)OFμ=argmax1K∑k=1KfCB,μ,Dt,Dv,where CB denotes the CatBoost classifier, μ=μ1,μ2,…,μn is the hyperparameters, Dt and Dv represent the training and validation set divided by the K-fold cross-validation, respectively, and f (CB, μ, Dt, Dv) is the classification accuracy.


### 2.3.4. The Proposed Fault Diagnosis Approach

An overview of the proposed fault diagnosis approach for rotating machinery based on feature importance ranking and selection is shown in Figure3. The specific steps are as follows:Step 1: the original vibration signals are acquired by accelerometers and the data acquisition system.Step 2: the optimal WBF is selected according to the proposed two-step principle, and then the original vibration signals are decomposed by OWPD.Step 3: the FE values of the decomposed subsignals are calculated to form the fault feature setF.Step 4: CatBoost algorithm is utilized to obtain the importance of each feature inF by a certain strategy. According to the ranking result of feature importance, the candidate features are selected in sequence and combined with the corresponding labels to form dataset S.Step 5: datasetS is divided into two parts according to a certain proportion: training set and test set.Step 6: the training set is used to train the CatBoost classifier, and BO algorithm is adopted to optimize the main hyperparameters.Step 7: the test set is fed into the trained CatBoost classifier to output the diagnostic results.Figure 3
The flowchart of the proposed approach.


## 2.3.1. CatBoost for Classification

In the GBDT algorithm, lots of decision trees are combined to produce a model of high accuracy, and the progress can be written as(17)yx=∑t=1Tftx,θt,where x is the feature vector, T is the number of trees, θtt=1,2,…,T is a learned parameter, and ft (x, θt) represents the decision trees that are learned.Given training samplesD=xk,yk1n, where n is the number of samples, xkk=1,2,…,n is the sample data, and yk is the true sample label. To learn the model introduced in equation (17), the following objective function needs to be minimized:(18)Οft=∑i=1nLyk,y^k+∑t=1TΩft,where y^k is the predicted sample label, L is the loss function that represents the difference between yk and y^k, and Ω is the regular function that is used to punish the complexity of ft, defined as(19)Ω=αq+12βω2,where α is the penalty parameter that controls the number of leaf nodes, q is the number of leaf nodes, β is the regularization parameter, and ω is the weight coefficient.Letg denote the negative gradient of the loss function, and the objective function is minimized in the direction of g, namely,(20)g=−∂Lyk,y^k∂y^k.Traditional GBDT algorithms generally have the problem of prediction offset, which affects the generalization ability of the model. To overcome this defect, CatBoost was proposed with two notable improvements [43]: (1) the ordered boosting strategy was adopted to obtain the unbiased estimation of the gradient and slow down the prediction offset; (2) the oblivious tree was used as the basic learner to increase the reliability of the model and speed up the prediction. In addition, to better deal with categorical features, the greedy target-based statistics strategy was improved by adding prior terms in CatBoost algorithm, which can be summarized as three main steps: (1) all the sample datasets are randomly arranged; (2) samples with the same category are selected, and the average label of similar samples is calculated; and (3) features of each sample are digitized by adding the prior term and its corresponding weight coefficient. The improved greedy target-based statistics strategy can be expressed as(21)x^ki=∑j=1nxji=xki⋅yi+a⋅P∑j=1nxji=xki+a,where xki represents the i-th category feature of the k-th sample, x^ki represents the corresponding numerical feature, P represents the increased prior value, and a represents the weight coefficient (a > 0). The addition of prior values can effectively reduce the noise caused by low-frequency features and avoid the overfitting phenomenon.


## 2.3.2. CatBoost for Feature Selection

The growth strategies of decision trees are different in different GBDT algorithms. XGBoost uses level-wise strategy, which has the disadvantage of inefficiency [53]. CatBoost uses the symmetric tree strategy to optimize the computation of the leaf value to prevent the model from overfitting. In the case of the basic learner of CatBoost is the tree model, the feature coefficient or importance can be obtained according to a certain evaluation index after training the model, e.g., the change of loss function or prediction values.For a given feature setF=f1,f2,…,fN, the feature importance of fii=1,2,…N in the trained CatBoost model is calculated by(22)featurefi=∑trees,leafSfiν1−avr2⋅c1+ν2−avr2⋅c2,and(23)avr=ν1⋅c1+ν2⋅c2c1+c2,where S denotes the different paths to the leaf nodes in the decision tree, c1 and c2 denote the total weight coefficient in the left and right leaves, respectively, and υ1 and υ2 denote the formula value in the left and right leaves, respectively.


## 2.3.3. BO Algorithm for Optimizing Hyperparameters

In machine learning algorithms, the hyperparameters usually have a strong influence on the performance of the model. There are some inevitable shortcomings in traditional methods of parameter tuning. For example, greedy algorithm can only obtain the local optimal solution, and the uncertainty and nonconvexity of grid search tend to miss global optimality. Different with these methods, BO algorithm can obtain the global optimal solution through Gaussian process, which is considered to find the optimal hyperparameters for the CatBoost model.The basic thought of BO algorithm is that, for the given data and optimal termination condition (usually, the number of iterations or the expected value of the objective function), Bayesian theory is used to estimate the posterior distribution. The distribution and the information from the previous sampling point are used to select the hyperparameters of the later sampling until that the value of the objective function reaches the maximum globally. Here, the objective function is defined as the maximum of the classification accuracy:(24)OFμ=argmax1K∑k=1KfCB,μ,Dt,Dv,where CB denotes the CatBoost classifier, μ=μ1,μ2,…,μn is the hyperparameters, Dt and Dv represent the training and validation set divided by the K-fold cross-validation, respectively, and f (CB, μ, Dt, Dv) is the classification accuracy.


## 2.3.4. The Proposed Fault Diagnosis Approach

An overview of the proposed fault diagnosis approach for rotating machinery based on feature importance ranking and selection is shown in Figure3. The specific steps are as follows:Step 1: the original vibration signals are acquired by accelerometers and the data acquisition system.Step 2: the optimal WBF is selected according to the proposed two-step principle, and then the original vibration signals are decomposed by OWPD.Step 3: the FE values of the decomposed subsignals are calculated to form the fault feature setF.Step 4: CatBoost algorithm is utilized to obtain the importance of each feature inF by a certain strategy. According to the ranking result of feature importance, the candidate features are selected in sequence and combined with the corresponding labels to form dataset S.Step 5: datasetS is divided into two parts according to a certain proportion: training set and test set.Step 6: the training set is used to train the CatBoost classifier, and BO algorithm is adopted to optimize the main hyperparameters.Step 7: the test set is fed into the trained CatBoost classifier to output the diagnostic results.Figure 3
The flowchart of the proposed approach.


## 3. Case I: Experimental Verification with the MFS Dataset

### 3.1. Experimental Setup and Data Description

To prove the effectiveness of the proposed fault diagnosis approach, a hybrid dataset of bearing and rotor faults collected by the machinery fault simulator (MFS) platform was used for experimental verification [54]. The experiment setup is shown in Figure 4, the MFS is driven by the AC motor with the speed of 2100 rpm, and the power of it is transmitted to the rotating plate and the drive shaft and through the coupling. The sampling frequency is 6 kHz. Through replacing different components, ten different types of datasets are collected with the data acquisition box, including nine fault types and one normal type, with detailed information shown in Table 1. There are 160 samples for each type, and each sample contains 1000 nonoverlapping data points.Figure 4
Experimental platform of the MFS dataset [54].Table 1
Details of ten types of faults.
Fault classLabelBall fault1Inner race fault2Outer race fault3Combination fault4Normal state5Spindle central bent6Couple bent7Cocked rotor8Unbalanced rotor9Eccentric rotor10To observe the difference between different types of vibration signals, a sample of each fault type is randomly selected to draw the waveform in time and frequency domains, respectively, and the results are shown in Figure5. It can be seen that all fault types are time-varying and frequency-varying signals, which indicates that the original vibration signals are nonstationary. In addition, in view of the frequency domain, most types of fault information are concentrated in the low-frequency band, while the high-frequency band contains less.Figure 5
Time- and frequency-domain waveforms of ten fault types. (a) Time domain. (b) Frequency domain.
(a)(b)Each sample data is standardized by thez-score method. In addition, missing values are detected in advance. If there is a missing value, it is filled with the Lagrange interpolation formula. The design and improvement of the experimental algorithms are implemented by Python 3.7.3 with a computer configured with Intel Core i5-6000hq CPU and 12G RAM.


### 3.2. Research on Feature Extraction

#### 3.2.1. Parameter Settings of OWPD

Different decomposition layers in OWPD will result in different frequency resolutions of subsignals, which affect the accuracy and time consumption of fault diagnosis. If the decomposition layer isl, the frequency resolution of the signal is(25)df=fs2l+1,where fs is the sampling frequency, and here, it is 6 kHz. To make sure that df is greater than 1 Hz, the value of l needs to be less than 12. In addition, the number of features and calculation time will increase with the increase in the number of subbands. Therefore, l is preliminarily selected as 5, and each sample is averagely separated into 32 parts in the frequency domain. The influence of the value of l on the diagnostic results will be analysed in detail in the following experiments.The optimal WBF will be selected according to the two-step principle described in Section2.1.2. To reduce the calculation time, 10 samples are randomly selected under each fault type to form a new dataset, and 5-level WPD is carried out for each sample with different WBFs. In the first step, the total energy-to-Shannon entropy ratio ζ of each sample is calculated, respectively, according to equations (5)–(8), and then the average value of these 100 samples is taken. The results are detailed in Table 2. As can be seen, the WBFs with the largest average ζ value are db7, sym7, and coif3 in the same wavelet family, respectively.Table 2
The average value of the total energy-to-Shannon entropy ratioζ.
WBFAverageζdb14.538db25.349db35.660db45.960db55.986db65.951db76.311db85.986db96.285db106.237sym25.349sym35.660sym45.854sym56.044sym65.968sym76.310sym86.151coif15.362coif25.857coif36.319coif46.057coif56.241————In the second step, the above candidate WBF is used to reconstruct the signal, and the average value of its similarity coefficientd with the original signal is calculated according to equation (9). The results are detailed in Table 3. As can be seen, the original signal is most similar to the reconstructed signal when coif3 WBF is selected. coif WBF has orthogonality and compact support. In addition, compared with db WBF, it has better symmetry. Therefore, it is more effective to extract fault vibration signals of impulsive and nonstationary characteristics.Table 3
The average value of the similarity coefficientd.
WBFdb7sym7coif3Similarity coefficientd (10−12)20.0369.2137.632


#### 3.2.2. Parameter Settings of FE

After determining the decomposition layer and optimal WBF, the original signals are decomposed by OWPD. Then, the FE values of these subsignals are calculated separately. That is to say, the dimension of the feature set is 32. The parameters of FE are set according to [14, 18], as shown in Table 4, where STD represents the standard deviation of each sample. Time consumption of feature extraction is 1.29 s/sample. To visualize the calculation results, a sample from each fault type is randomly selected to carry out 5-layer OWPD using coif3 WBF. Figure 6 shows the FE values of the subsignals of different fault types. It can be seen that the difference of FE values among different fault types is quite obvious, which indicates that the proposed OWPD-FE method can effectively extract fault features. In addition, the FE values of some subsignals for certain fault types are almost 0, which also indicates that there is some feature redundancy.Table 4
Parameter settings of the FE method.
ParameterDescriptionValueΛTime delay1mEmbedding dimension2rSimilarity tolerance0.15∗ STDnGradient of similarity tolerance2Figure 6
The FE values of subsignals decomposed by OWPD.


#### 3.2.3. Feature Visualization

t-SNE is the most commonly used algorithm for data visualization and dimensionality reduction [55]. Here, it is adopted to project the extracted 32-dimensional features into two-dimensional (2D) space for visualization. For each fault type, 20 groups of data after dimensionality reduction by t-SNE are selected to draw a scatter diagram, as shown in Figure 7. It is clear that, in the 2D plane, the features of ten fault types overlap little, and the boundaries among different types can be clearly distinguished, which shows that the OWPD-FE method is effective in extracting fault information of bearings and rotors.Figure 7
2D projection of the extracted fault features using t-SNE.


### 3.3. Feature Selection Using CatBoost

An appropriate feature selection method can reduce the feature redundancy and improve the diagnostic performance of the model. In this paper, CatBoost algorithm is considered for this step. According to Section2.3.2, equations (22) and (23) are used to get the values of importance of the 32 fault features extracted by the OWPD-FE method, according to which the feature importance can be ranked in the descending order. Then, candidate features and their corresponding labels are added sequentially to train and test the model. Finally, through the relationship between the number of features and the classification accuracy, the feature set used to reach the highest classification accuracy can be obtained. Figure 8 shows the normalized calculation results of the importance of 32 features, whose sum is 100.Figure 8
The value of the normalized feature importance.


### 3.4. Diagnosis Results and Analysis

At first, the dataset is divided into training set and test set, and the ratio of them is set as 3 : 2. That is to say, 96 samples of each fault type are used for model training and the rest of 64 samples for test. In addition, ten-fold cross-validation is conducted on the training set. The main hyperparameters of CatBoost optimized by BO algorithm are shown in Table5.Table 5
The main hyperparameters of CatBoost optimized by BO algorithm.
ParameterValueLoss functionCross-entropyNumber of estimators450L2 regularization3Learning rate0.254Max depth of one tree4Random seed number50Then, to study the effect of the number of features on the classification results, the feature selection process is carried out according to the analysis of Section3.3. The experimental results are detailed in Figure 9. It can be seen that the model training time is positively correlated with the number of features, which is consistent with the actual experience. The average accuracy of ten-fold cross-validation on the training set has reached 100% when 7 features are selected. When the number of features is 22, the test set accuracy is the highest, reaching 99.17%. However, when all 32 features are used, it decreased by 0.21 percent (only 98.96%). Therefore, the classification accuracy of using 22 features is considered as the final diagnosis result. The time consumption of model training in this case is 10.13 s, which is 1.68 s less than that without feature selection. Experimental results show the reliability of the proposed feature selection method and the effectiveness of the classification algorithm.Figure 9
Diagnosis results and time consumption using different number of features.As indicated in Figure10, the confusion matrix of the diagnosis result using 22 features is presented in detail. It is not hard to see that the diagnostic accuracy of all fault types is above 98%, and it reaches 100% for 6 fault types (corresponding category labels are 2, 3, 6, 7, 9, and 10). The experimental results show that the proposed approach can effectively identify the hybrid fault states of the rotor and bearing.Figure 10
Confusion matrix of using 22 features selected by CatBoost.


### 3.5. Comparison of Different Decomposition Layers

In this section, in order to further demonstrate the effectiveness and reasonability of setting decomposition layerl to 5, the influence of the value of l on the diagnostic performance is investigated. Firstly, l needs to be less than 12 according to equation (25). Therefore, l is set to 1 to 8, respectively, to perform OWPD. Then, the parameters of FE are set according to Table 4. The partition of the sample dataset is described in Section 3.4. CatBoost is still applied to feature selection and fault pattern recognition of the 8 datasets, and BO algorithm is used to optimize hyperparameters. The experimental results are shown in Figure 11.Figure 11
The change of diagnostic results with the number of decomposition layers. (a) Classification accuracy. (b) Feature extraction and model training time.
(a)(b)As can be seen from Figure11(a), the classification accuracy generally presents an upward trend with the increase of l. In detail, when l is equal to 1, the average accuracy of the training set and validation set is very low, only 82.02% and 53.75%, respectively, while the test set accuracy only reaches 57.29%. When l is greater than 2, all the training sets’ accuracy reaches 100%. The validation set accuracy reaches the maximum when l is equal to 6 (99.29%). When l is equal to or greater than 5, all the test sets’ accuracy reaches 99.17%. As can be seen from Figure 11(b), the time of feature extraction and model training increases with the increase of l. The feature extraction time increases exponentially with the increase of l. Experimental results show that high-quality features containing effective fault information can be obtained by selecting appropriate l to perform WPD. An inappropriate value of l will result in too low classification accuracy (l less than 4) or too high computational cost (l greater than 6). Therefore, it is reasonable to set l to 5 considering the classification accuracy and time consumption comprehensively.In addition, to illustrate the necessity of feature selection, the diagnostic performance with and without feature selection is compared, and the experimental results are shown in Figure12. As can be seen from Figure 12(a), there is no feature redundancy when l is less than 4. However, the number of redundant features gradually increases when l is greater than 4, which is caused by the increase of frequency bands containing no or less fault information. Accordingly, the implementation of feature selection greatly reduces model training time from the perspective of computational cost. Meanwhile, it can be seen from Figure 12(b) that feature selection can also effectively improve the classification accuracy, especially when l is 5, 7, and 8.Figure 12
Comparison of diagnostic results with or without feature selection. (a) Model training time and the number of features used for classification. (b) Test set accuracy.
(a)(b)


### 3.6. Comparison of Different Classifiers

In this section, to justify the superiority of the proposed BO-CatBoost algorithm and the applicability of the OWPD-FE method combined with other classifiers, SVM, RF, GBDT, and XGBoost are adopted for comparison. The dataset consisting of 22 high-quality fault features described in Section3.4 is input into the above classifiers, respectively, for model training and testing. To obtain the optimal diagnostic performance, the main hyperparameters of these classifiers are all optimized by BO algorithm, as detailed in Table 6. The diagnostic results are shown in Figure 13. It can be seen that the average training set accuracy of RF, GBDT, and CatBoost all reaches 100%, while SVM is the lowest, only 97.88%. In terms of the validation set accuracy, SVM has the lowest average accuracy (only 95.27%), while it is the highest of CatBoost (98.39%). RF has the lowest standard deviation, indicating that it is relatively stable. In terms of the test set accuracy, CatBoost is the highest, followed by GBDT and XGBoost (both 98.54%), while SVM performs the worst. In practical applications, different classifiers are applicable to different task requirements. Experimental results show that the diagnostic performance of GBDT algorithm and its variant (XGBoost) is close to that of CatBoost under the premise of setting appropriate hyperparameters. In addition, it is also demonstrated that the OWPD-FE method can extract high-quality fault features that are easy to identify.Table 6
The main hyperparameters of different classifiers optimized by BO algorithm.
ClassifierParameterValueSVMPenalty factor1799Coefficient of the kernel function0.833Kernel functionrbfDecision function shapeovoRFNumber of estimators928Max depth of one tree3Maximum number of features allowed for a single tree0.1Nf1Minimum number of samples needed to split a node6GBDTNumber of estimators208Learning rate0.115Max depth of one tree2Subsampling ratio0.795XGBoostNumber of estimators516Learning rate0.129Max depth of one tree4Proportion of randomly sampled features in each tree0.584Nf1L1 regularization0.105L2 regularization0.591Subsampling ratio0.1Booster modelgbtreeObjective functionmulti: softmaxNumber of classes101Nf denotes the total number of features.Figure 13
Comparison of classification accuracy of different classifiers.


## 3.1. Experimental Setup and Data Description

To prove the effectiveness of the proposed fault diagnosis approach, a hybrid dataset of bearing and rotor faults collected by the machinery fault simulator (MFS) platform was used for experimental verification [54]. The experiment setup is shown in Figure 4, the MFS is driven by the AC motor with the speed of 2100 rpm, and the power of it is transmitted to the rotating plate and the drive shaft and through the coupling. The sampling frequency is 6 kHz. Through replacing different components, ten different types of datasets are collected with the data acquisition box, including nine fault types and one normal type, with detailed information shown in Table 1. There are 160 samples for each type, and each sample contains 1000 nonoverlapping data points.Figure 4
Experimental platform of the MFS dataset [54].Table 1
Details of ten types of faults.
Fault classLabelBall fault1Inner race fault2Outer race fault3Combination fault4Normal state5Spindle central bent6Couple bent7Cocked rotor8Unbalanced rotor9Eccentric rotor10To observe the difference between different types of vibration signals, a sample of each fault type is randomly selected to draw the waveform in time and frequency domains, respectively, and the results are shown in Figure5. It can be seen that all fault types are time-varying and frequency-varying signals, which indicates that the original vibration signals are nonstationary. In addition, in view of the frequency domain, most types of fault information are concentrated in the low-frequency band, while the high-frequency band contains less.Figure 5
Time- and frequency-domain waveforms of ten fault types. (a) Time domain. (b) Frequency domain.
(a)(b)Each sample data is standardized by thez-score method. In addition, missing values are detected in advance. If there is a missing value, it is filled with the Lagrange interpolation formula. The design and improvement of the experimental algorithms are implemented by Python 3.7.3 with a computer configured with Intel Core i5-6000hq CPU and 12G RAM.


## 3.2. Research on Feature Extraction

### 3.2.1. Parameter Settings of OWPD

Different decomposition layers in OWPD will result in different frequency resolutions of subsignals, which affect the accuracy and time consumption of fault diagnosis. If the decomposition layer isl, the frequency resolution of the signal is(25)df=fs2l+1,where fs is the sampling frequency, and here, it is 6 kHz. To make sure that df is greater than 1 Hz, the value of l needs to be less than 12. In addition, the number of features and calculation time will increase with the increase in the number of subbands. Therefore, l is preliminarily selected as 5, and each sample is averagely separated into 32 parts in the frequency domain. The influence of the value of l on the diagnostic results will be analysed in detail in the following experiments.The optimal WBF will be selected according to the two-step principle described in Section2.1.2. To reduce the calculation time, 10 samples are randomly selected under each fault type to form a new dataset, and 5-level WPD is carried out for each sample with different WBFs. In the first step, the total energy-to-Shannon entropy ratio ζ of each sample is calculated, respectively, according to equations (5)–(8), and then the average value of these 100 samples is taken. The results are detailed in Table 2. As can be seen, the WBFs with the largest average ζ value are db7, sym7, and coif3 in the same wavelet family, respectively.Table 2
The average value of the total energy-to-Shannon entropy ratioζ.
WBFAverageζdb14.538db25.349db35.660db45.960db55.986db65.951db76.311db85.986db96.285db106.237sym25.349sym35.660sym45.854sym56.044sym65.968sym76.310sym86.151coif15.362coif25.857coif36.319coif46.057coif56.241————In the second step, the above candidate WBF is used to reconstruct the signal, and the average value of its similarity coefficientd with the original signal is calculated according to equation (9). The results are detailed in Table 3. As can be seen, the original signal is most similar to the reconstructed signal when coif3 WBF is selected. coif WBF has orthogonality and compact support. In addition, compared with db WBF, it has better symmetry. Therefore, it is more effective to extract fault vibration signals of impulsive and nonstationary characteristics.Table 3
The average value of the similarity coefficientd.
WBFdb7sym7coif3Similarity coefficientd (10−12)20.0369.2137.632


### 3.2.2. Parameter Settings of FE

After determining the decomposition layer and optimal WBF, the original signals are decomposed by OWPD. Then, the FE values of these subsignals are calculated separately. That is to say, the dimension of the feature set is 32. The parameters of FE are set according to [14, 18], as shown in Table 4, where STD represents the standard deviation of each sample. Time consumption of feature extraction is 1.29 s/sample. To visualize the calculation results, a sample from each fault type is randomly selected to carry out 5-layer OWPD using coif3 WBF. Figure 6 shows the FE values of the subsignals of different fault types. It can be seen that the difference of FE values among different fault types is quite obvious, which indicates that the proposed OWPD-FE method can effectively extract fault features. In addition, the FE values of some subsignals for certain fault types are almost 0, which also indicates that there is some feature redundancy.Table 4
Parameter settings of the FE method.
ParameterDescriptionValueΛTime delay1mEmbedding dimension2rSimilarity tolerance0.15∗ STDnGradient of similarity tolerance2Figure 6
The FE values of subsignals decomposed by OWPD.


### 3.2.3. Feature Visualization

t-SNE is the most commonly used algorithm for data visualization and dimensionality reduction [55]. Here, it is adopted to project the extracted 32-dimensional features into two-dimensional (2D) space for visualization. For each fault type, 20 groups of data after dimensionality reduction by t-SNE are selected to draw a scatter diagram, as shown in Figure 7. It is clear that, in the 2D plane, the features of ten fault types overlap little, and the boundaries among different types can be clearly distinguished, which shows that the OWPD-FE method is effective in extracting fault information of bearings and rotors.Figure 7
2D projection of the extracted fault features using t-SNE.


## 3.2.1. Parameter Settings of OWPD

Different decomposition layers in OWPD will result in different frequency resolutions of subsignals, which affect the accuracy and time consumption of fault diagnosis. If the decomposition layer isl, the frequency resolution of the signal is(25)df=fs2l+1,where fs is the sampling frequency, and here, it is 6 kHz. To make sure that df is greater than 1 Hz, the value of l needs to be less than 12. In addition, the number of features and calculation time will increase with the increase in the number of subbands. Therefore, l is preliminarily selected as 5, and each sample is averagely separated into 32 parts in the frequency domain. The influence of the value of l on the diagnostic results will be analysed in detail in the following experiments.The optimal WBF will be selected according to the two-step principle described in Section2.1.2. To reduce the calculation time, 10 samples are randomly selected under each fault type to form a new dataset, and 5-level WPD is carried out for each sample with different WBFs. In the first step, the total energy-to-Shannon entropy ratio ζ of each sample is calculated, respectively, according to equations (5)–(8), and then the average value of these 100 samples is taken. The results are detailed in Table 2. As can be seen, the WBFs with the largest average ζ value are db7, sym7, and coif3 in the same wavelet family, respectively.Table 2
The average value of the total energy-to-Shannon entropy ratioζ.
WBFAverageζdb14.538db25.349db35.660db45.960db55.986db65.951db76.311db85.986db96.285db106.237sym25.349sym35.660sym45.854sym56.044sym65.968sym76.310sym86.151coif15.362coif25.857coif36.319coif46.057coif56.241————In the second step, the above candidate WBF is used to reconstruct the signal, and the average value of its similarity coefficientd with the original signal is calculated according to equation (9). The results are detailed in Table 3. As can be seen, the original signal is most similar to the reconstructed signal when coif3 WBF is selected. coif WBF has orthogonality and compact support. In addition, compared with db WBF, it has better symmetry. Therefore, it is more effective to extract fault vibration signals of impulsive and nonstationary characteristics.Table 3
The average value of the similarity coefficientd.
WBFdb7sym7coif3Similarity coefficientd (10−12)20.0369.2137.632


## 3.2.2. Parameter Settings of FE

After determining the decomposition layer and optimal WBF, the original signals are decomposed by OWPD. Then, the FE values of these subsignals are calculated separately. That is to say, the dimension of the feature set is 32. The parameters of FE are set according to [14, 18], as shown in Table 4, where STD represents the standard deviation of each sample. Time consumption of feature extraction is 1.29 s/sample. To visualize the calculation results, a sample from each fault type is randomly selected to carry out 5-layer OWPD using coif3 WBF. Figure 6 shows the FE values of the subsignals of different fault types. It can be seen that the difference of FE values among different fault types is quite obvious, which indicates that the proposed OWPD-FE method can effectively extract fault features. In addition, the FE values of some subsignals for certain fault types are almost 0, which also indicates that there is some feature redundancy.Table 4
Parameter settings of the FE method.
ParameterDescriptionValueΛTime delay1mEmbedding dimension2rSimilarity tolerance0.15∗ STDnGradient of similarity tolerance2Figure 6
The FE values of subsignals decomposed by OWPD.


## 3.2.3. Feature Visualization

t-SNE is the most commonly used algorithm for data visualization and dimensionality reduction [55]. Here, it is adopted to project the extracted 32-dimensional features into two-dimensional (2D) space for visualization. For each fault type, 20 groups of data after dimensionality reduction by t-SNE are selected to draw a scatter diagram, as shown in Figure 7. It is clear that, in the 2D plane, the features of ten fault types overlap little, and the boundaries among different types can be clearly distinguished, which shows that the OWPD-FE method is effective in extracting fault information of bearings and rotors.Figure 7
2D projection of the extracted fault features using t-SNE.


## 3.3. Feature Selection Using CatBoost

An appropriate feature selection method can reduce the feature redundancy and improve the diagnostic performance of the model. In this paper, CatBoost algorithm is considered for this step. According to Section2.3.2, equations (22) and (23) are used to get the values of importance of the 32 fault features extracted by the OWPD-FE method, according to which the feature importance can be ranked in the descending order. Then, candidate features and their corresponding labels are added sequentially to train and test the model. Finally, through the relationship between the number of features and the classification accuracy, the feature set used to reach the highest classification accuracy can be obtained. Figure 8 shows the normalized calculation results of the importance of 32 features, whose sum is 100.Figure 8
The value of the normalized feature importance.


## 3.4. Diagnosis Results and Analysis

At first, the dataset is divided into training set and test set, and the ratio of them is set as 3 : 2. That is to say, 96 samples of each fault type are used for model training and the rest of 64 samples for test. In addition, ten-fold cross-validation is conducted on the training set. The main hyperparameters of CatBoost optimized by BO algorithm are shown in Table5.Table 5
The main hyperparameters of CatBoost optimized by BO algorithm.
ParameterValueLoss functionCross-entropyNumber of estimators450L2 regularization3Learning rate0.254Max depth of one tree4Random seed number50Then, to study the effect of the number of features on the classification results, the feature selection process is carried out according to the analysis of Section3.3. The experimental results are detailed in Figure 9. It can be seen that the model training time is positively correlated with the number of features, which is consistent with the actual experience. The average accuracy of ten-fold cross-validation on the training set has reached 100% when 7 features are selected. When the number of features is 22, the test set accuracy is the highest, reaching 99.17%. However, when all 32 features are used, it decreased by 0.21 percent (only 98.96%). Therefore, the classification accuracy of using 22 features is considered as the final diagnosis result. The time consumption of model training in this case is 10.13 s, which is 1.68 s less than that without feature selection. Experimental results show the reliability of the proposed feature selection method and the effectiveness of the classification algorithm.Figure 9
Diagnosis results and time consumption using different number of features.As indicated in Figure10, the confusion matrix of the diagnosis result using 22 features is presented in detail. It is not hard to see that the diagnostic accuracy of all fault types is above 98%, and it reaches 100% for 6 fault types (corresponding category labels are 2, 3, 6, 7, 9, and 10). The experimental results show that the proposed approach can effectively identify the hybrid fault states of the rotor and bearing.Figure 10
Confusion matrix of using 22 features selected by CatBoost.


## 3.5. Comparison of Different Decomposition Layers

In this section, in order to further demonstrate the effectiveness and reasonability of setting decomposition layerl to 5, the influence of the value of l on the diagnostic performance is investigated. Firstly, l needs to be less than 12 according to equation (25). Therefore, l is set to 1 to 8, respectively, to perform OWPD. Then, the parameters of FE are set according to Table 4. The partition of the sample dataset is described in Section 3.4. CatBoost is still applied to feature selection and fault pattern recognition of the 8 datasets, and BO algorithm is used to optimize hyperparameters. The experimental results are shown in Figure 11.Figure 11
The change of diagnostic results with the number of decomposition layers. (a) Classification accuracy. (b) Feature extraction and model training time.
(a)(b)As can be seen from Figure11(a), the classification accuracy generally presents an upward trend with the increase of l. In detail, when l is equal to 1, the average accuracy of the training set and validation set is very low, only 82.02% and 53.75%, respectively, while the test set accuracy only reaches 57.29%. When l is greater than 2, all the training sets’ accuracy reaches 100%. The validation set accuracy reaches the maximum when l is equal to 6 (99.29%). When l is equal to or greater than 5, all the test sets’ accuracy reaches 99.17%. As can be seen from Figure 11(b), the time of feature extraction and model training increases with the increase of l. The feature extraction time increases exponentially with the increase of l. Experimental results show that high-quality features containing effective fault information can be obtained by selecting appropriate l to perform WPD. An inappropriate value of l will result in too low classification accuracy (l less than 4) or too high computational cost (l greater than 6). Therefore, it is reasonable to set l to 5 considering the classification accuracy and time consumption comprehensively.In addition, to illustrate the necessity of feature selection, the diagnostic performance with and without feature selection is compared, and the experimental results are shown in Figure12. As can be seen from Figure 12(a), there is no feature redundancy when l is less than 4. However, the number of redundant features gradually increases when l is greater than 4, which is caused by the increase of frequency bands containing no or less fault information. Accordingly, the implementation of feature selection greatly reduces model training time from the perspective of computational cost. Meanwhile, it can be seen from Figure 12(b) that feature selection can also effectively improve the classification accuracy, especially when l is 5, 7, and 8.Figure 12
Comparison of diagnostic results with or without feature selection. (a) Model training time and the number of features used for classification. (b) Test set accuracy.
(a)(b)


## 3.6. Comparison of Different Classifiers

In this section, to justify the superiority of the proposed BO-CatBoost algorithm and the applicability of the OWPD-FE method combined with other classifiers, SVM, RF, GBDT, and XGBoost are adopted for comparison. The dataset consisting of 22 high-quality fault features described in Section3.4 is input into the above classifiers, respectively, for model training and testing. To obtain the optimal diagnostic performance, the main hyperparameters of these classifiers are all optimized by BO algorithm, as detailed in Table 6. The diagnostic results are shown in Figure 13. It can be seen that the average training set accuracy of RF, GBDT, and CatBoost all reaches 100%, while SVM is the lowest, only 97.88%. In terms of the validation set accuracy, SVM has the lowest average accuracy (only 95.27%), while it is the highest of CatBoost (98.39%). RF has the lowest standard deviation, indicating that it is relatively stable. In terms of the test set accuracy, CatBoost is the highest, followed by GBDT and XGBoost (both 98.54%), while SVM performs the worst. In practical applications, different classifiers are applicable to different task requirements. Experimental results show that the diagnostic performance of GBDT algorithm and its variant (XGBoost) is close to that of CatBoost under the premise of setting appropriate hyperparameters. In addition, it is also demonstrated that the OWPD-FE method can extract high-quality fault features that are easy to identify.Table 6
The main hyperparameters of different classifiers optimized by BO algorithm.
ClassifierParameterValueSVMPenalty factor1799Coefficient of the kernel function0.833Kernel functionrbfDecision function shapeovoRFNumber of estimators928Max depth of one tree3Maximum number of features allowed for a single tree0.1Nf1Minimum number of samples needed to split a node6GBDTNumber of estimators208Learning rate0.115Max depth of one tree2Subsampling ratio0.795XGBoostNumber of estimators516Learning rate0.129Max depth of one tree4Proportion of randomly sampled features in each tree0.584Nf1L1 regularization0.105L2 regularization0.591Subsampling ratio0.1Booster modelgbtreeObjective functionmulti: softmaxNumber of classes101Nf denotes the total number of features.Figure 13
Comparison of classification accuracy of different classifiers.


## 4. Case II: Experiment Verification with the One-Stage Reduction Gearbox Dataset

### 4.1. Experimental Setup and Data Description

Since the working condition of the MFS dataset is relatively simple, this section further verifies the effectiveness of the proposed approach in practical applications through an actual gearbox dataset with more complex working conditions. The experimental platform is composed of a one-stage reduction gearbox, a torque sensor, a servo motor, etc., as shown in Figure14. Four fault types are formed by processing gears with different crack lengths (0, 5, 10, and 15 mm). The sample frequency is 5 kHz. The details of the dataset can be found in [2]. Here, data collected under 20 different working conditions are used, as shown in Table 7. These data constitute 10 different datasets, as shown in Table 8. Datasets D1 to D9 contain a relatively simple number of working conditions. Composed of all 20 working conditions, D10 is the most complex dataset that is closest to the actual working conditions. For a single working condition, there are 40 samples for each fault type. The total number of samples is 3200, each containing 1500 consecutive data points.Figure 14
Experimental platform for the one-stage gearbox [2].Table 7
20 different working conditions of the one-stage gearbox.
Shaft speed (rpm)Load (N·m)02468600W1W2W3W4W5900W6W7W8W9W101200W11W12W13W14W151500W16W17W18W19W20Table 8
Details of 10 datasets used for experimental verification.
LabelWorking conditionsSample sizeD1W1,W2,W3,W4,W5800D2W6,W7,W8,W9,W10800D3W11,W12,W13,W14,W15800D4W16,W17,W18,W19,W20800D5W1,W6,W11,W16640D6W2,W7,W12,W17640D7W3,W8,W13,W18640D8W4,W9,W14,W19640D9W5,W10,W15,W20640D10W1,W2,…,W203200


### 4.2. Diagnosis Results and Analysis

The steps and settings of feature extraction refer to the descriptions in Section3.2, where sym5 is the optimal WBF. The BO-CatBoost algorithm is still used to identify and diagnose fault types. Figure 15 shows the diagnostic results of ten datasets after selecting the optimal feature subset. For all ten datasets, the training set accuracy reaches 100%. For D1 to D9, the test set accuracy is higher than 96.56% and even higher than 99% on D2, D6, D7, and D9. The test set accuracy of D10 under the most complicated conditions is 98.65%, which is 0.22 percent higher than using the default CatBoost hyperparameters. However, in terms of time consumption, it takes 206.09 s to complete the training process with default hyperparameters, while 6.94 s with hyperparameters optimized by BO algorithm, which greatly improves the computational efficiency. Therefore, experimental results in this study show that the proposed approach is also effective for fault diagnosis of the gearbox.Figure 15
The classification results of the 10 datasets used for experimental verification.


## 4.1. Experimental Setup and Data Description

Since the working condition of the MFS dataset is relatively simple, this section further verifies the effectiveness of the proposed approach in practical applications through an actual gearbox dataset with more complex working conditions. The experimental platform is composed of a one-stage reduction gearbox, a torque sensor, a servo motor, etc., as shown in Figure14. Four fault types are formed by processing gears with different crack lengths (0, 5, 10, and 15 mm). The sample frequency is 5 kHz. The details of the dataset can be found in [2]. Here, data collected under 20 different working conditions are used, as shown in Table 7. These data constitute 10 different datasets, as shown in Table 8. Datasets D1 to D9 contain a relatively simple number of working conditions. Composed of all 20 working conditions, D10 is the most complex dataset that is closest to the actual working conditions. For a single working condition, there are 40 samples for each fault type. The total number of samples is 3200, each containing 1500 consecutive data points.Figure 14
Experimental platform for the one-stage gearbox [2].Table 7
20 different working conditions of the one-stage gearbox.
Shaft speed (rpm)Load (N·m)02468600W1W2W3W4W5900W6W7W8W9W101200W11W12W13W14W151500W16W17W18W19W20Table 8
Details of 10 datasets used for experimental verification.
LabelWorking conditionsSample sizeD1W1,W2,W3,W4,W5800D2W6,W7,W8,W9,W10800D3W11,W12,W13,W14,W15800D4W16,W17,W18,W19,W20800D5W1,W6,W11,W16640D6W2,W7,W12,W17640D7W3,W8,W13,W18640D8W4,W9,W14,W19640D9W5,W10,W15,W20640D10W1,W2,…,W203200


## 4.2. Diagnosis Results and Analysis

The steps and settings of feature extraction refer to the descriptions in Section3.2, where sym5 is the optimal WBF. The BO-CatBoost algorithm is still used to identify and diagnose fault types. Figure 15 shows the diagnostic results of ten datasets after selecting the optimal feature subset. For all ten datasets, the training set accuracy reaches 100%. For D1 to D9, the test set accuracy is higher than 96.56% and even higher than 99% on D2, D6, D7, and D9. The test set accuracy of D10 under the most complicated conditions is 98.65%, which is 0.22 percent higher than using the default CatBoost hyperparameters. However, in terms of time consumption, it takes 206.09 s to complete the training process with default hyperparameters, while 6.94 s with hyperparameters optimized by BO algorithm, which greatly improves the computational efficiency. Therefore, experimental results in this study show that the proposed approach is also effective for fault diagnosis of the gearbox.Figure 15
The classification results of the 10 datasets used for experimental verification.


## 5. Conclusions

In this paper, aiming at the fault diagnosis of rotating machinery under complex working conditions, a novel approach based on feature importance ranking and selection is proposed. Firstly, the OWPD method is proposed to decompose the vibration signal, where a two-step principle of selecting the optimal WBF is introduced. On this basis, it is combined with FE to extract hidden and high-quality fault features from the decomposed subsignals. Then, in order to filter out redundant fault features that are not conducive to the diagnosis result, the CatBoost model is constructed and preliminarily applied to calculate the importance of each feature for further feature selection. Moreover, the classification model based on BO-CatBoost algorithm can effectively solve the optimization problem of hyperparameters, which can greatly reduce model training time and improve the diagnosis accuracy. Finally, experimental results on the MFS dataset and the one-stage gearbox dataset under complex working conditions demonstrate the practicability and the generalization performance of the proposed approach, and the classification accuracy reaches 99.17% and above 96.56%, respectively. In addition, the robustness of the proposed approach under different working conditions is also verified by the one-stage gearbox dataset.In the future work, the effect of combining OWPD with other information entropies or some dimensionless time-domain indexes still needs to be discussed. In addition, the vibration data collected by only one acceleration sensor are utilized in this paper, while the fusion of multisensor data may provide more real fault information, which is also worth investigating in the next step.



---

*Source: 8899188-2021-02-04.xml*