#### Abstract

Early fault diagnosis of rolling element bearing is still a difficult problem. Firstly, in order to effectively extract the fault impulse signal of the bearing, a new enhanced morphological difference operator (EMDO) is constructed by combining two optimal feature extraction-type operators. Next, in the process of processing the test signal, in order to reduce the interference problem caused by strong background noise, the probabilistic principal component analysis (PPCA) method is introduced. In the traditional PPCA method, two important system parameters (decomposition principal component *k* and original variable *n*) are usually set artificially; this will greatly reduce the noise reduction performance of PPCA. To solve this problem, a parameter adaptive PPCA method based on grasshopper optimization algorithm (GOA) is proposed. Finally, combining the advantages of the above algorithms, a comprehensive rolling bearing fault diagnosis method named APPCA-EMDF is proposed in this paper. Experimental comparison results show that the proposed method can effectively diagnose the vibration signals of rolling element bearing.

#### 1. Introduction

Rolling bearings are widely used in the field of rotating machinery. It is vulnerable to damage because it operates at high speed and variable loads. This not only brings difficulties to the safety production of the machinery but also causes huge economic losses to enterprises. When a local fault occurs in the rolling bearing, an impulse signal of periodic vibration will be generated [1, 2]. These fault signals are usually subject to background noise and harmonic interference caused by other equipment. Therefore, the early weak fault extraction of rolling bearings is very necessary [3].

A number of detection techniques have been developed to address this issue, such as empirical mode decomposition (EMD) [4], spectral kurtosis (SK) [5], and variational mode decomposition (VMD) [6]. Then, traditional time-frequency analysis methods tend to denoise, while these algorithms also have their own shortcomings. Different from the above fault diagnosis analysis methods, morphological filter (MF) directly changes the geometry of the original signal through the structural element (SE); at the same time, MF can also handle nonlinear and nonsteady signals [7]. At first, MF was applied to image processing by Matheron and Serra [8], and later it was widely used in mechanical fault diagnosis due to its simple and fast advantages [9–14]. Two important factors affect the performance of MF filter: the selection of SE and morphological operator. For the SE, a large number of literatures show that the length of SE has a great influence on the filtering performance of the signal [15–17]. For the morphological operator, the construction of morphological operators also plays an important role.

Morphological operators are mainly composed of four basic operators: opening, closing, dilation, and erosion. At present, based on these four basic morphological operators, a variety of new morphological operators are cascaded together, and they have been successfully applied to vibration signals. MG is composed of the difference between the dilation and erosion operators. Raj and Murali [18], Li et al. [19], and Hu et al. [20] successfully applied the MG operator to the fault diagnosis of rolling bearings. Dong et al. [21] introduced the AVG operator with the average of opening and closing operators. Hu and Xiang [22] proposed the CMF operator by using the average of the opening-closing and closing-opening operators and demonstrated the denoising performance of CMF. Wang et al. [23] and Meng et al. [24] used the CMF operator to detect bearing failure. DIF is composed of the difference between the closing and opening operators. Li et al. [25, 26] and Zhang et al. [27] successfully applied the DIF operator to the fault diagnosis of mechanical systems.

Previous scholars used the morphological operator with single structure to deal with the fault of rolling bearings. The results show that the morphological operator with a single structure can no longer meet the requirements when dealing with signals with complex frequency components. Therefore, some scholars have studied that the new morphological operator constructed by the product can not only enhance the feature extraction ability but also weaken the interference of random noise. Lv and Yu [13] used basic morphological operators to construct some combination difference operators and chose two optimal operators from them to construct the average combination difference operator (ACDIF). Then, they used the TEK index factor to select the optimal scale of the operator and finally applied the proposed method to fault diagnosis of rolling bearing. Li et al. [15] made a detailed qualitative analysis of some common basic morphological operators. According to the purpose of processing the signal, these morphological operators are divided into two categories: noise reduction-type operators and feature extraction-type operators, and they constructed the MGPO operator based on the superior performance of the two morphological operators. Based on Li’s research, Yan et al. [28] derived a new operator named MHPO, and they successfully applied MHPO to the fault diagnosis of rolling bearings.

Unfortunately, in the above studies, few papers have quantitatively studied the filtering performance of these morphological operators. In this paper, a set of rolling bearing fault simulation signals and four sets of filtering performance indicators are used to quantitatively analyze the feature extraction capability of the feature extraction-type operator and rank their filtering performance. Then, a new enhancement operator named EMDO is constructed.

Although EMDO can extract more impulse signals, due to the interference of background noise, the filtering performance of EMDO will be affected. In order to solve this problem, the probabilistic principal component analysis (PPCA) [29–32] is introduced to make up for the deficiency of EMDO in denoising capability. PPCA first establishes an appropriate probability model for each variable, and it decomposes the principal components and noise of the signal into orthogonal space by increasing the dimension, realizes the decorrelation between principal components and noise, and finally achieves the purpose of noise reduction. In the PPCA model, two important system parameters, the number of decomposed principal components *k* and the original variable *n*, have a great impact on the model’s noise reduction results. At present, existing papers have manually selected the system parameters of PPCA, which greatly affects the noise reduction performance [33].

In order to obtain the optimal parameters of PPCA, this paper chooses grasshopper optimization algorithm (GOA) to optimize the system parameters of PPCA. GOA is a mathematical optimization algorithm proposed by Saremi et al. [34] in 2017, which is used to simulate the behavior of grasshoppers in nature. Compared with other optimization methods, GOA has the characteristics of fast convergence and difficulty in generating local optimal solutions [34]. Finally, this paper combines the advantages of the above two methods to propose a new comprehensive fault diagnosis method for rolling bearing fault diagnosis.

The remainder of this article is organized as follows. Section 2 briefly introduces mathematical morphological theory and proposes the EMDO operator. Section 3 describes the APPCA method for noise reduction as well as the construction process of the APPCA-EMDF algorithm. Experimental verification of the proposed method is illustrated in Section 4. Conclusions are drawn in Section 5.

#### 2. Basic Theory of Mathematical Morphology

##### 2.1. Morphological Filter

Assume that the original signal *f* (*n*) is a group of one-dimensional discrete array *F* = (0, 1, 2, …, *N* − 1). The structure element (*m*) is defined as another one-dimensional discrete array *G* = (0, 1, 2, …, *M* − 1) and satisfies *N* ≥ *M*. Then, the four basic morphological operators can be expressed as follows.

Dilation:

Erosion:

Opening:

Closing:where , , ○, and ● represent the dilation, erosion, opening, and closing operators, respectively. Besides, several main morphological filters are introduced as follows.

Average of closing and opening (AVG) [21]:

Morphological gradient (MG) [18–20]:

Difference filter (DIF) [25–27]:

The opening-closing *F*_{OC} and closing-opening *F*_{CO} filters are defined as

A combination morphological filter with *F*_{OC} and *F*_{CO} (CMF) [23,24] is defined as

*F*_{CO} and *F*_{OC} gradient operation (*G*_{CO&OC}) [35]:

Black top-hat (BTH) operation and positive BTH (PBTH) operation [27]:

White top-hat (WTH) operation and negative WTH (NWTH) operation [22]:

In order to verify the filtering performance of the morphological operator, the SE needs to be selected. Some scholars [17, 19, 35] have proved that only the length of SE has a great impact on the filtering results. Therefore, in order to improve the calculation efficiency, this paper chooses a flat SE with zero height. The relationship between the length of SE and the scale is *L* = *λ* + 2, and the selection range of scale is usually from 1 to *f*_{s}/*f*_{o}, where *f*_{s} and *f*_{o} represent sampling frequency and fault frequency, respectively.

##### 2.2. Property of Morphology Operators

According to the purpose of processing fault signals, morphological operators can generally be divided into two categories: noise reduction-type operations and feature extraction-type operations [15, 17]. The noise reduction-type operations tend to denoise the signal, while the feature extraction-type operations are more inclined to extract fault features. The specific classification results of the 15 morphological operators introduced in Section 2.1 are shown in Table 1.

In order to extract weak fault signals from strong background noise, this paper mainly studies the filtering performance of feature extraction-type operators. Li et al. [25] introduced the characteristics of some basic morphological operators in detail, but they did not quantitatively study the fault feature extraction capabilities of morphological operators under different SNRs. Therefore, this paper quantitatively evaluates the filtering performance of the feature extraction-type operators. In order to analyze their engineering applicability to rolling bearing faults, a set of bearing fault simulation signal models is established [36]. It is as follows:

The fault model of the inner race of a bearing is composed of four parts: the inner race fault signal , the random impulse interference signal *y*_{2}(*t*), the harmonic interference signal , and the white Gaussian noise signal *δ*(*t*). Here, a local fault is generated in the inner race of the bearing. The fault frequency of the signal *f*_{o} = 16 Hz, the number of random impulses of the signal is 3, and the resonant frequency is 300 Hz. The sampling frequency of the signal is 1024 Hz, and the number of sampling points is 1024.

Feature energy factor (FEF) [17, 28] and amplitudes at the 3th harmonics of the bearing fault frequency (3thA) [37,38] have been proved to be an effective index factor for evaluating the feature extraction capability of morphological operators. As the statistical measurement index, permutation entropy (PE) [25, 39] and envelope spectrum sparsity (ESS) [40] represent the randomness and sparsity of signals, which can effectively reflect the intensity of background noise in the test signal, and they are a set of effective indicator factors for measuring signal denoising performance. Therefore, in this paper, four representative evaluation factors are selected to prove the filtering performance of these feature extraction-type operators. The four evaluation indicators are shown in Table 2.

For signal *y*(*t*), flat SE with a height of zero is selected and the scale is 5. When the SNR changes from −10 dB to 0 dB, the variation of four evaluation index factors with the SNR is shown in Figure 1.

**(a)**

**(b)**

**(c)**

**(d)**

For FEF factor, the larger the FEF, the stronger the operator’s fault feature extraction ability. In formula *F*1 of Table 2, *M* = 4 and *f* = 200 Hz. It can be clearly seen from Figure 1(a) that the FEF of the 7 operators increases with the increase of SNR. MG received the largest FEF, followed by DIF. WTH and NWTH obtain the minimum FEF. For PE factor, if the signal is more regular, then the PE is smaller. In Figure 1(b), the PE index of 7 operators does not increase significantly with the change of SNR. The PE index of MG is the global minimum, followed by *G*_{co_oc} and DIF. It shows that the filtering performance of these three operators is better than the other four morphological operators. The greater the ESS, the stronger the denoising performance of the operator. As shown in Figure 1(c), when SNR changes from −10 dB to −6 dB, the ESS of *G*_{co_oc} is the best. When SNR changes from −6 dB to 0 dB, MG has the best ESS. The ESS index of DIF is stable between 5 and 6, and its ESS ranks the third place among all operators. BTH, PBTH, NWTH, and WTH operators obtain almost the same ESS, indicating their poor filtering performance. The 3thA is a key index to evaluate the feature extraction ability of the operator. The 3thA of these operators with SNR is shown in Figure 1(d). In Figure 1(d), when SNR increased from −10 dB to 0 dB, MG received the largest 3thA than other operators, while the 3th A of DIF was slightly lower than the MG, ranking second. Among all operators, the 3thA index factor of *G*_{co_oc} was the lowest.

It can be seen from the above analysis that the feature extraction ability of BTH, WTH, PBTH, and NWTH is poor, while the MG operator and DIF operator show superior filtering performance.

Research through reference [17, 19] shows that for the simulation signal *y*(*t*) of equation (16), the optimal scale of these operators is 1–10. In addition, if the scale of the SE is too large, some useful information will be removed. In order to further prove the filtering performance of the feature extraction-type operator at different scales, when the SNR of signal *y*(*t*) is −3, the evaluation indicator of the operator varying with the scale is shown in Figure 2.

**(a)**

**(b)**

**(c)**

**(d)**

In Figure 2(a), when the scale is between 1 and 8, MG received the largest FEF than other operators. When the scale is between 9 and 10, DIF obtains the best FEE. Therefore, from FEF factor analysis, it is concluded that MG and DIF show excellent feature extraction performance. In Figure 2(b), the PE value of the operators decreases and becomes stable with the increase of scale. Compared with other operators, MG obtains the smallest PE in the whole scale, followed by *G*_{co_oc}, and DIF ranks third. The PE values of PBTH, NWTH, BTH, and WTH are larger, indicating that their filtering performance is not as good as the MG, DIF, and *G*_{co_oc} operators. The results of ESS index factor of these operators changing with scale are shown in Figure 2(c). In Figure 2(c), the change trend of the operator’s ESS with the scale is opposite to PE. The same conclusion can be obtained from the analysis results. In Figure 2(d), the variation trend of 3thA index of these operators is similar to FEF. MG and DIF still get the largest value at full scale.

Through the comprehensive analysis of the evaluation indexes of the above operators, it can be concluded that the rank of feature extraction ability of 7 operators is as follows: MG > DIF > *G*_{co_oc} > BTH = PBTH > NWTH = WTH

Based on the above quantitative research on filtering performance of feature extraction operator, it can be concluded that MG and DIF have superior fault feature extraction ability. Therefore, a new morphological differential operator named EMDO is constructed based on the excellent performance of the two operators.

The newly proposed EMDO operator has two advantages: (1) if the fault information is detected by both MG and DIF, then EMDO will amplify the fault information; (2) because Gaussian noise shows randomness, when Gaussian noise is detected by only one of them, EMDO will output a smaller amplitude.

##### 2.3. Comparison of Morphology Operators

When SNR of signal *y*(*t*) is −6 dB, its time-domain diagram and frequency-domain diagram are shown in Figure 3, respectively. In Figure 3(b), the fault frequency *f*_{o} can no longer be detected due to the interference of strong background noise, so the simple envelope analysis can no longer be used to detect early weak fault of the rolling bearings.

**(a)**

**(b)**

In order to further verify the superiority of the proposed operator in this paper, five other morphological operators (MG [19], DIF [25], MGPO [15], MHPO [28], and ACDIF [13]) were selected to compare with the proposed operator in fault feature extraction capabilities. Since the FEF index factor is sensitive to the characteristic information of the signal, it can directly reflect the operator’s ability to extract fault characteristic information. This paper chooses FEF to optimize the scales of the operators. In formula *F*1 of Table 2*M* = 4 and *f* = 200 Hz. The FEF values of above morphological operators changing with the scale are shown in Figure 4.

It can be seen from Figure 4 that before reaching the optimal scale, the FEF of each operator increases with the scale, and when it reaches the optimal scale, it decreases with the increase of the scale. The optimal scale distribution of these operators is between 1 and 10. When the optimal scale is 4, EMDO obtains the highest FEF of 33.67%, which is larger than other operators. The effectiveness and superiority of EMDO have been proven. The FEFs of MG and DIF were significantly lower than EMDO. This shows that the newly constructed operator is enhanced in feature information extraction capability. The maximum FEFs of MG, DIF, MGPO, MHPO, and ACDIF operators are 23.37%, 22.61%, 26.32%, 25.14% and 23.36%, respectively. The FEF values of the above six operators are far greater than the envelope spectrum analysis. The results of processing the signal *y*(*t*) in Figure 3(a) by these operators under their optimal scale are shown in Figure 5.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

It can be seen from Figure 5 that all six operators can extract the fault characteristic frequency *f*_{o}, but 2*f*_{o} extracted by MG, DIF, and ACDIF is weak. Due to the interference of background noise, the frequency components in the range of 0–200 Hz are seriously polluted. Although MGPO can extract the fault frequencies *f*_{o}, the amplitude of the fault frequency is relatively weak, and 4*f*_{o} cannot be detected. MHPO operator can detect 4*f*_{o}; however, there are many interference frequency components in the range of 0–200 Hz. Especially when the background noise further increases, MHPO is not as robust as EMDO. As shown in Figure 5(f), the operator proposed in this paper can clearly detect the fault frequencies *f*_{o}, 2*f*_{o}, 3*f*_{o}, and 4*f*_{o}. Combined with the above analysis, it is further proved that the EMDO is superior to other morphological operators in the extraction ability of fault features.

#### 3. The Method of APPCA-EMDF

It can be seen from the above analysis that the advantages of feature extraction-type operators lie in the ability of fault feature extraction, but their denoising performance is insufficient. In order to make up the deficiency of feature extraction-type operator in denoising ability, this paper presents an adaptive PPCA denoising method.

##### 3.1. PPCA

In the PPCA model [29–32], it is assumed that the main fault information of the original signal is stored in the principal component space, and the remaining noise information and linear correlation information are saved in the remaining subspace. The PPCA model for the *n*-dimensional original variable *X* is as follows:where is a matrix (it is generated by a one-dimensional signal where is the test failure data and is the mean of ) of *n* × *m*, *n* is the dimension of the original variable *X*, and *m* is the length of the signal. is a matrix of *n* × *k*, and the limiting condition is *k* < *n*, where *k* is the number of main components. is the principal component matrix of *k* × *m*. **E** is an isotropic Gaussian noise matrix; **u** and **E** both satisfy Gaussian distribution, , where **I** is the identity matrix and *σ*^{2} is the variance of the noise variable. Therefore, **X** follows a Gaussian distribution .

The probability distribution of **u** is as follows:

The prior probability distribution of **x** under the variable **u** is as follows:

According to equations (19) and (20), the probability distribution of the original data *x* is as follows:where is the covariance matrix, which is determined by the parameters **P** and *σ*^{2}.

According to equations (19)–(21), once the values of the parameters **P** and *σ*^{2} are determined, the probability model is obtained. The above parameters can be estimated using the expectation maximization (EM) algorithm [41], and the iterative formula derived is as follows:where **S** is the covariance matrix of the original data and is the trace of the matrix. The parameters **P** and *σ*^{2} are obtained by iterative convergence of equations (22) and (23). Finally, a PPCA model is established.

When the PPCA model is established, the solution method of the reduced-dimensional data is shown below.

From equation (24), each principal component signal (data after dimensionality reduction) is a projection of the original variable data **X** on the corresponding principal component vector **p**_{i}. Although PPCA is a better noise reduction algorithm, the system parameters of PPCA (the number of principal components *k* and the original variable *n*) need to be determined by personal experience, which will seriously affect its noise reduction performance. In PPCA, the parameter *n* determines the dimension of the construction matrix, and the parameter *k* determines the number of eigenvalues after dimensionality reduction.

In order to test the noise reduction performance after signal processing, Deng et al. [40] proposed a new and effective signal noise reduction index SIESS by combining ESS and PE. Because the kurtosis *K*_{r} is suitable for early fault diagnosis of rolling bearings, this paper proposed a new dimensionless comprehensive evaluation index KSP to detect the noise reduction performance of the signal.where the expressions of ESS and PE are shown in Table 2. The larger the KSP, the better noise reduction performance of the signal. Therefore, the maxKSP can indicate the optimal noise reduction result of PPCA. For the signal *y*(*t*) in Figure 3(a), the distribution of KSP under different system parameters (*n*, *k*) is shown in Figure 6.

From Figure 6, the distribution of KSP is relatively discrete and random with the change of *n* and *k*. When *n* = 19 and *k* = 8, the maximum value of KSP is 14.98. The noise reduction effect of PPCA is affected by the unreasonable setting of parameters *n* and *k*. Therefore, in order to achieve the best noise reduction effect of PPCA, this article chooses GOA to adaptively optimize the system parameters of PPCA. The GOA [34] is described as follows.

##### 3.2. GOA

GOA is a mathematical optimization algorithm that simulates the behavior of grasshopper swarms in nature. There are two distinct main features of the grasshopper swarm behavior: (1) during the larval stage, the grasshoppers move slowly and the number of steps is small; the grasshopper behavior in adulthood shows a long distance of movement and agility; (2) the process by which grasshopper swarms look for food sources can be divided into two trends: exploration and exploitation. During exploration, the search swarms are encouraged to move suddenly, and in exploitation, they are encouraged to move locally. Grasshopper swarms naturally accomplish both functions and seek goals. Therefore, the mathematical model used to simulate the behavior of the grasshopper swarm is as follows:where *X*_{i} is the position of the *i*th grasshopper, *S*_{i} represents social interaction, *G*_{i} defines the gravity force on the *i*th grasshopper, and *A*_{i} represents wind advection.where indicates the distance between the *i*th and *j*th grasshopper (the interval of *d*_{ij} is usually 1 to 4 in the GOA), shows a unit vector from the *i*th grasshopper to the *j*th grasshopper, *s*(•) is used to represent the function of social forces between grasshoppers, represents the intensity of attraction, *l* indicates attractive length scale, and and *l* are two important parameters of the GOA (they are described in [34]). Usually, *f* and *l* are equal to 1.5 and 0.5.where represents the gravitational constant and is a unity vector towards the center of Earth.where *u* indicates a constant drift and represents a unity vector of the wind direction.

Substituting *S*, *G*, and *A* into equation (26):where the expression of *s*(•) is as shown in equation (28) and *N* is the number of grasshoppers.

However, the above mathematical model cannot be directly used to solve the optimization problem because the grasshoppers reach the comfort zone too fast and the swarms cannot converge in the designated area. Therefore, the following modified equation is used to solve the above problem [34].where ub_{d} is the upper boundary of *D*-dimensional space, lb_{d} is the lower boundary of *D*-dimensional space, indicates the optimal target found so far in the *D*th dimension, and *c* represents a decreasing coefficient; its equation is as follows:where *c*_{max} represents the maximum value of the decreasing coefficient *c*, *c*_{min} represents the minimum value of the decreasing coefficient *c*, *l* indicates the current number of iterations, and *L* indicates the total number of iterations. In this paper, *c*_{min} and *c*_{max} are taken as 0.00001 and 1, respectively.

The mathematical model of the GOA is described above. The specific optimization process is as follows. Firstly, initialize the swarm and randomly generate the swarms *X*_{i}, initialize system parameters (e.g. *c*_{max}, *c*_{min}, *l*, and *L*), calculate the fitness of each search agent, and let *T* be equal to the optimal fitness. Secondly, using equation (32) to update the position of the current target search agent, in order to avoid the optimization model falling into the local optimal solution, the decreasing coefficient *c* is updated each time by using equation (33), and the distance between grasshoppers is normalized in each iteration. Finally, when the iteration condition is satisfied, the global optimal solution is obtained. A more detailed introduction to the GOA is given in [34].

##### 3.3. APPCA-EMDF

The proposed parameter adaptive PPCA method takes the maxKSP value as the objective function and optimizes the parameters through the GOA. The specific objective function is shown in the following equation:where fitness represents the objective function of the optimization, KSP_{i} (*i* = 1, 2, …, *γ*) is the KSP index after PPCA decomposition and reconstruction of the signal, and *γ* = (*n*, *k*) represents the optimized system parameters of PPCA. This article proposes an early weak fault detection method for rolling bearings based on a combination of adaptive PPCA and enhanced morphological differential filtering.

The specific flowchart is shown in Figure 7.

The specific working steps are as follows: Step 1: input the vibration signal *y*(*t*), set the optimization range of PPCA parameters (*n*, *k*), and determine the initial parameters of GOA (the maximum iteration steps *L* and the number of population *N*) [34]. Step 2: in the process of each iteration, the signal *y*(*t*) is decomposed and the fitness of the signal is calculated after each decomposition. The maximum fitness is retained. Step 3: determine the termination condition of the program. If *l* ≥ *L*, return yes; otherwise, return no; then *l* = *l* + 1, and continue iteration until the condition is satisfied [34]. Step 4: save the optimal system parameters (*n*, *k*) and obtain the maximum objective function value. PPCA applies optimal system parameters to the signal for noise reduction. Step 5: the morphological analysis scale range is determined and the SE is determined. Step 6: calculate the FEF value at each scale and determine the optimal scale *λ* for EMDF. Step 7: EMDF uses the optimal scale *λ* to perform the morphological calculation. Step 8: envelope analysis finally realizes fault diagnosis.

#### 4. Simulation

##### 4.1. Case 1: The Traditional Bearing Fault Model

The method proposed in this paper is used to process the simulation signal *y*(*t*) in Figure 3(a). The range of optimization parameters (*n*, *k*) of APPCA is selected according to equation (34). The parameters of GOA are selected as follows: maximum number of iterations *L* = 10 and search agents *N* = 20. The GOA convergence curve of the maxKSP value obtained after each iteration is shown in Figure 8(a). In Figure 8(a), when the number of iterations reaches 3, the curve starts to converge. The optimal *n*, *k*, and maxKSP values obtained by GOA are 16, 8, and 14.981, respectively. PPCA processes the signal *y*(*t*) with the optimal *n* and *k*. The noise reduction result is shown in Figure 8(b). Comparing Figure 8(b) with Figure 3(a), it can be clearly found that the background noise is effectively suppressed. Next, the calculated full-scale FEF value of the denoising signal is shown in Figure 9(a). In Figure 9(a), it can be seen that the optimal scale obtained by EMDF is 4. The signal of Figure 8(b) is processed with EMDF, and the time-domain and frequency-domain diagrams of the calculation results are shown in Figures 9(b) and 9(c), respectively. Comparing the calculated results with Figures 5(a)–5(f), it can be found that the APPCA-EMDF method can clearly detect the bearing fault characteristic frequency of 16 Hz and its harmonic frequencies of 32 Hz, 48 Hz, and 64 Hz. At the same time, the background noise is effectively suppressed, and the feature extraction capability of the fault is greatly enhanced.

**(a)**

**(b)**

**(a)**

**(b)**

**(c)**

##### 4.2. Case 2: Complex Fault Model of Bearings

In order to further verify the anti-interference ability and weak fault extraction ability of the algorithm proposed in this paper, a new set of rolling bearing simulation fault model is established, whose expression is as follows:

The fault model of the bearing is composed of five parts: the fault signal *y*_{1}(*t*), the random pulse interference signal *y*_{2}(*t*), the harmonic interference signals *y*_{3}(*t*) and *y*_{4}(*t*), and the white Gaussian noise signal *δ*(*t*). In the model, *A* represents the amplitude of the fault signal, *T* represents the time interval between the two shock signals, and *τ* represents a random variable that is usually used to simulate random sliding of the roller bearing. *B* is a variable representing the amplitude of the random impulse signal. *C* represents the amplitude of the harmonic interference signal, *f* represents the frequency of the harmonic signal, and *φ* represents the phase angle. *s*(*t*) represents the impulse response function of the mechanical system and can be described as follows:where *M* represents the amplitude, *a* is the decay factor, *f*_{n} is the resonance frequency, and *φ* is the phase angle.

A simulated signal is established based on equation (35). Here, a local fault is generated in the inner race of the bearing. The fault frequency of the signal *f*_{o} = 16 Hz, the amplitude *A* = 4, the decay factor *a* is set to 100, the natural frequency *f*_{n} = 200 Hz, and the phase is zero. The number of random impulses of the signal is 10, and the resonant frequency is 300 Hz. In harmonic interference, the harmonic frequencies of the signals *y*_{3}(*t*) are 30 Hz and 40 Hz, with an amplitude of 1.2 and 1.1, respectively. The harmonic frequencies of the signals *y*_{4}(*t*) are 30 Hz and 50 Hz, with an amplitude of 1 and 1.2, respectively. The SNR of white Gaussian noise is −10 dB. The sampling frequency of the signal is 2048 Hz, and the number of sampling points is 2048. The components of the simulated signal are shown in Figure 10. The obtained simulated signal *y*′(*t*) and its envelope spectrum are shown in Figure 11. It can be found in Figure 11(b) that due to the serious interference of harmonic frequency and background noise, only 10 Hz, 20 Hz, and 50 Hz interference frequencies can be identified, making it difficult to detect the bearing fault frequency.

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

The method proposed in this paper is used to process the simulation signal *y*′(*t*). The obtained GOA optimization curve and the time-domain diagram after noise reduction are shown in Figures 12(a) and 12(b), respectively. In Figure 12(a), when iterating to 4, the objective function value converges, the maximum KSP value obtained is 12.945, and the optimized PPCA parameters *n* = 18 and *k* = 7. Next, feature extraction is carried out for the denoised signal. The processed results are shown in Figures 13(a)–13(c). In Figure 13(c), the fault characteristic frequencies 16 Hz, 32 Hz, and 48 Hz still can be clearly detected. This shows that the algorithm proposed in this paper has strong anti-interference ability and weak fault feature extraction ability.

**(a)**

**(b)**

**(a)**

**(b)**

**(c)**

For comparison, the three algorithms (ACDIF [13], VMD [6], and SK [5]) were used to process the simulation signal in Figure 11(a), and the processing results are shown in Figures 14–16. In Figure 14(c), ACDIF can only detect a weak fault frequency 16 Hz, and the harmonic interference frequency 10 Hz is clearly present. In Figure 15(c), although VMD can identify the fault frequency of 16 Hz, the frequency components of 32 Hz and 48 Hz have weak amplitude, and they are seriously disturbed. It can be seen in Figure 16(c) that the SK method cannot effectively detect the fault characteristics of the bearing if the background noise and other frequency interference components are serious.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

#### 5. Experimental Validations

##### 5.1. Case 1: Bearing Fault Diagnosis of Wind Turbine

Wind energy is used as a clean energy source, and the number of wind turbines in the world is gradually increasing. Due to the long-term service of wind turbines in extreme environments, the transmission system of wind turbines is easily damaged. Once the wind turbine is damaged, it will bring huge economic losses to the wind turbine operator. Therefore, fault monitoring of wind turbines is particularly important.

By monitoring the vibration of multiple wind turbines in Tuoshan Wind Farm, the abnormal vibration amplitude of the front bearing of No. 16 wind turbine can be found, so the fault monitoring is carried out. The structural model of the wind turbine and its sensor arrangement are shown in Figure 17(a). The test system monitors 8 positions in the cabin (test points 1–8 in Figure 17(a)). During the test, abnormal vibration signal of measuring point 2 was found, and the field measured photos are shown in Figure 17(b). The tested generator front-end bearing model is SKF NU1030. The specific structural geometric parameters of the bearing are shown in Table 3. The rated power of the wind turbine is 1500 kW, the blade speed is stable at 22.3 rpm, the motor speed is 1740 rpm, and the rotation frequency of the generator *f*_{r} = 29 Hz. The sampling frequency of the test system is 16384 Hz, and the total length of the test signal is 81920. The fault frequencies of the inner ring, outer ring, cage, and rolling element of the rolling bearing calculated from Table 4 are *f*_{i} = 384.7 Hz, *f*_{o} = 315.9 Hz, *f*_{c} = 13.11 Hz, and *f*_{B} = 97.1 Hz.

**(a)**

**(b)**

The time-domain diagram and the FFT diagram of the collected generator front-end vibration signals are shown in Figures 18(a) and 18(b), respectively. In Figure 18(a), impulse signals with different energy amplitudes can be seen, but due to the interference of background noise, the periodic vibration impulse signals of the bearings are not very obvious, so further processing and analysis are needed. In Figure 18(b), the vibration energy of the test signal is mainly concentrated in the middle frequency band. In the low frequency band, the weak outer ring failure frequency of the rolling bearing can only be identified at 315.6 Hz.

**(a)**

**(b)**

Therefore, the method proposed in this paper is used to process the experimental signal to detect the bearing fault at the front of the generator. Firstly, the experimental signal was denoised using PPCA, and the GOA was used to adaptively optimize the system parameters (*n*, *k*) of PPCA. The obtained objective function convergence curve is shown in Figure 19(a). In Figure 19(a), when iterating to 3, the objective function converges, the maxKSP obtained is 1.6684, and the optimal system parameters of PPCA are *n* = 20 and *k* = 8. The optimal parameters are used to denoise the experimental signals, and the results are shown in Figure 19(b). Compared with Figure 18(a), it can be found that the signal denoised by APPCA presents relatively clean impulse vibration. It can be seen in Figure 19(c) that there is an impulse vibration signal with a period interval of 3.2 ms, which is fully consistent with the characteristic frequency of the outer ring fault of the rolling bearing of 315.6 Hz. Finally, in order to further enhance the fault information, the EMDF is used to process the denoised signal. In formula *F*1 of Table 2, *M* = 3 and *f* = 1000 Hz. As can be seen from Figure 20(a), when the scale is 3, EMDF obtains the largest FEF factor. The time-domain and frequency-domain results of the denoised signals processed using the optimal scale are shown in Figures 20(b) and 20(c), respectively. In Figure 20(c), in addition to the rotation frequency of 29 Hz, the fault frequency of bearing outer ring of 315.8 Hz and its harmonic frequencies of 630 Hz and 946.9 Hz can be clearly detected. The analysis results show that the method proposed in this paper can both detect the bearing faults in the front end of the generator and enhance the fault information.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

For comparison, the three algorithms (ACDIF, VMD, and SK) were still used to process the experimental signal in Figure 18(a), and the processing results are shown in Figures 21–23. In the ACDIF method, the optimal scale selection of the operator is shown in Figure 21(a), and the processed frequency-domain results are shown in Figure 21(c). In Figure 21(c), the rotation frequency of 29 Hz can be detected, but the amplitude of the fault frequency *f*_{o} of the bearing outer ring is relatively weak, and the effect of fault feature extraction is not obvious. In the VMD method, the penalty factor *α* is 2000, and the decomposition mode *k* is 5. According to the maximum weight kurtosis index, the most sensitive mode selected is u5, and the selection result is shown in Figure 22(a). The final processing result is shown in Figure 22(c). The fault frequency *f*_{o} and its harmonic frequency 2*f*_{o} of the bearing outer ring can be detected, but the amplitude of the fault characteristic frequency is relatively weak, and the fault feature extraction capability is not strong. In Figure 23(c), from the results obtained by processing the experimental signals with the SK method, the fault frequency *f*_{o} is severely contaminated due to interference from the background noise and can hardly be detected.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

##### 5.2. Case 2: Early Fault Diagnosis of Bearings

In order to verify the early fault feature extraction ability of the proposed algorithm, the bearing life cycle data collected by the University of Cincinnati’s Intelligent Maintenance System (IMS) [42] were applied. Bearing test benches, sensors, and their schematic diagrams are shown in Figures 24(a) and 24(b), respectively. During the experiment, a total of three sets of experimental data were tested. In this paper, the second set of experimental data is selected for processing. This group of data consists of 984 files. Each group of data is collected every 10 minutes. The specific geometric parameters of the tested bearing are shown in Table 5. During the experiment, the motor speed was 2000 rpm and the sampling frequency was 20 kHz. According to Table 5, the fault frequency of bearing outer ring (*f*_{BPFO}) is 236.4 Hz. The root mean square (RMS) value of the full life cycle of the bearing is shown in Figure 25(a). As can be seen in Figure 25(a), the damage form of the bearing is divided into three stages: early, middle, and late stages, corresponding to the A, B, and C points in the Figure 25(a). Their corresponding time-domain diagrams are shown in Figures 25(b)–25(d), respectively. It can be seen from Figure 25(b) that the impulse amplitude of bearing early fault is relatively weak.

**(a)**

**(b)**

**(a)**

**(b)**

**(c)**

**(d)**

In this paper, the time-domain signal in Figure 25(b) is processed. Its FFT spectrum and envelope spectrum are shown in Figures 26(a) and 26(b), respectively. In Figure 26(b), only the weak fault frequency *f*_{FBPO} can be identified, and the bearing fault information is seriously disturbed by background noise.

**(a)**

**(b)**

The method proposed in this paper is used to process the experimental data. First, the GOA optimization curve and the noise reduction results are shown in Figures 27(a) and 27(b), respectively. In Figure 27(a), when iterating to 4, the objective function value converges, the maximum KSP value obtained is 17.059, and the optimized PPCA parameters *n* = 15 and *k* = 6. By comparing Figure 27(b) with Figure 25(b), it can be found that significant bearing fault impulse signals can be detected after noise reduction processing by PPCA. Next, feature extraction is performed on the denoised signal. The processed results are shown in Figures 28(a) to 28(c). In Figure 28(c), the fault characteristic frequency (*f*_{BPFO}) of the bearing outer ring can be clearly detected, and its corresponding high-order characteristic frequency (2*f*_{BPFO}, 3*f*_{BPFO}, and 4*f*_{BPFO}) can also be detected. The fault feature extraction effect is obvious.

**(a)**

**(b)**

**(a)**

**(b)**

**(c)**

ACDIF, VMD, and SK were used to process the same experimental signals for comparison. The experimental results obtained are shown in Figure 29–Figure 31, respectively. In Figure 29(b), although the fault frequency (*f*_{BPFO}, 2*f*_{BPFO}, and 3*f*_{BPFO}) of the bearing outer ring can be detected, the fault characteristics are more seriously polluted by noise and there are many interference components. It can be seen from Figures 30(b) and 31(b) that the fault feature extraction capabilities of the VMD method and the SK method are still relatively weak.

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

Through the comparison of the above experiments and engineering applications, it can be clearly found that the method proposed in this paper can effectively extract the fault characteristic information of rolling bearings. This method not only has good noise reduction performance but also can play a role in enhancing the early weak faults and has good filtering performance.

#### 6. Discussion

The previous analysis results can qualitatively prove that the method proposed is effective. In order to further compare the superiority of the APPCA-EMDF in this paper with the existing decomposition methods (ACDIF, VMD, and SK), some quantitative indicators are applied, such as kurtosis and FEF indicators. The larger the kurtosis and FEF index, the stronger the feature extraction ability of the signal. At the same time, the calculation efficiency of the above methods is also compared. The test was performed on a computer configured with Inter (R) Core (TM) i7-6700HQ, CPU 2.6 Hz, and RAM 8 GB. The comparison results are shown in Tables 6–8, respectively.

From the comparison results of Table 6–Table 8, it can be seen that although the SK method has the highest computational efficiency, its fault feature extraction ability is the weakest. The APPCA-EMDF method in this paper obtains the largest kurtosis value and FEF value, which shows that its fault feature extraction ability is better than other decomposition methods (ACDIF, VMD, and SK). However, the method proposed in this article is less computationally efficient. The reason for this phenomenon is that GOA is used to optimize PPCA system parameters, and a lot of calculation costs are caused in the process of iterative optimization. Therefore, we need to further improve the computational efficiency of the algorithm in the future.

In future work, we can try to use the multiscale morphological analysis method because morphological operators can contain more abundant fault information components at full scale. We can also try to combine the algorithm proposed in this paper with the intelligent recognition algorithm to achieve accurate classification of bearing faults. Finally, we can try to find out how the noise elimination performance of PPCA is under the condition of other interference signals. How to improve the calculation efficiency of the algorithm proposed in this paper is also a problem to be studied in the future.

#### 7. Conclusions

In this paper, a new rolling element bearing fault method named APPCA-EMDF is proposed. By quantitatively analyzing the fault extraction capability of the feature extraction-type operators, a morphological operator named EMDO was proposed. Compared with MG, DIF, MGPO, MHPO, and ACDIF operators, the fault extraction ability of EMDO is superior, and it has an enhanced effect on fault characteristics. Next, in order to solve the interference problem caused by strong background noise, a parameter adaptive PPCA method based on the GOA method is proposed. APPCA can make up for the deficiency of EMDF in denoising capability. Simulation and engineering application results show that the method proposed in this paper is effective in detecting early fault of rolling element bearings. Compared with ACDIF, VMD, and SK methods, the proposed method has certain advantages.

#### Nomenclature

PPCA: | Probabilistic principal component analysis |

EMDF: | Enhanced morphological difference filter |

EMDO: | Enhanced morphological difference operator |

GOA: | Grasshopper optimization algorithm |

MF: | Morphological filter |

SE: | Structural element |

MG: | Morphology gradient |

AVG: | Average of opening and closing operators |

CMF: | Average of the opening-closing and closing-opening operators |

ACDIF: | Average combination difference morphological filter |

TEK: | Teager energy kurtosis |

MGPO: | Morphology gradient product operation |

MHPO: | Morphology-hat product operation |

BTH: | Black top-hat |

PBTH: | Positive black top-hat |

WTH: | White top-hat |

NWTH: | Negative white top-hat |

G_{CO&OC}: | Closing-opening and opening-closing gradient operation |

FEF: | Feature energy factor |

PE: | Permutation entropy |

ESS: | Envelope spectrum sparsity |

3thA: | Amplitudes at the 3rd harmonics of the bearing fault frequency |

VMD: | Variational mode decomposition |

SK: | Spectral kurtosis. |

#### Data Availability

The experimental data were acquired from the Prognostics Center of Excellence (PCoE) through the prognostic data repository donated by the Center of Intelligent Maintenance Systems (IMS), University of Cincinnati (J. Lee, H. Qiu, G. Yu, J. Lin, R.T. Services, Bearing data set, in: U.O. Cincinnati (Ed.), IMS, University of Cincinnati, NASA Ames Prognostics Data Repository, 2007.)

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This study was supported by the National Natural Science Foundation of China under grant nos. 51675350 and 51705337, Natural Science Foundation of Liaoning Province under grant nos. 2019MS245 and 20180551036 and the PhD Start-Up Foundation of Liaoning Province under grant no. 20170520111.