Compound Fault Diagnosis of a Wind Turbine Gearbox Based on MOMEDA and Parallel Parameter Optimized Resonant Sparse Decomposition

Abstract

Wind turbines usually operate in harsh environments. The gearbox, the key component of the transmission chain in wind turbines, can easily be affected by multiple factors during the operation process and develop compound faults. Different types of faults can occur, coupled with each other and staggered interference. Thus, a challenge is to extract the fault characteristics from the composite fault signal to improve the reliability and the accuracy of compound fault diagnosis. To address the above problems, we propose a compound fault diagnosis method for wind turbine gearboxes based on multipoint optimal minimum entropy deconvolution adjusted (MOMEDA) and parallel parameter optimized resonant sparse decomposition (RSSD). Firstly, the MOMEDA is applied to the preprocess, setting the deconvolution period with different fault frequency types to eliminate the interference of the transmission path and environmental noise, while decoupling and separating the different types of single faults. Then, the RSSD method with parallel parameter optimization is applied for decomposing the preprocessed signal to obtain the low resonance components, further suppressing the interference components and enhancing the periodic fault characteristics. Finally, envelope demodulation of the enhanced signal is applied to extract the fault features and identify the different fault types. The effectiveness of the proposed method was verified using the actual data from the wind turbine gearbox. In addition, a comparison with some existing methods demonstrates the superiority of this method for decoupling composite fault characteristics.

1. Introduction

Wind turbines, gas turbines and other advanced equipment are used widely in modern industry. The gearbox, a key component in these devices, is prone to failure when running under severe operating conditions such as heavy loads, large temperature differences, corrosive media, and alternating loads [1,2]. As the structure of machinery and equipment tends to be large and complex, the actual operation of the weak single fault will also have a chain reaction with the components of the transmission chain, resulting in various faults occurring in successive cascades to form a compound fault. Therefore, it is important to accurately and reliably extract fault features from vibration signals to achieve a composite fault diagnosis. In addition, this current problem is rarely covered in existing research, which mostly focuses on the compound failure caused by the local failure of the bearing and does not consider the compound failure caused by the bearing failure through the transmission path. The abovementioned situation is the starting point for this study and the basis for this paper.

Feature extraction aims at extracting feature information from the vibration signal to describe the operational status of the mechanical equipment. Over the past two decades, many scholars have explored the rotating machinery fault diagnosis field and introduced many diagnostic theories and methods. For example, empirical mode decomposition (EMD) [3,4], wavelet transform (WT) [5,6], variational mode decomposition fault (VMD) [7,8], and so on. Although these methods and their combination perform well on single faults, there are some limitations: for example, the existence of modal mixing in EMD [9,10], the diagnosis effectiveness of wavelet transform, which depends on the constant quality factor, and the choice of wavelet basis [11,12]. The lack of adaptiveness to complex signals and good separation and decoupling performance makes it impossible to effectively extract fault components from complex, random, and variable vibration signals.

The vibration signal collected by the sensor is regarded as a convolutional mixture of different excitation sources, fault sources and transmission channels, so the recovery of the fault signal can be considered a deconvolution process. In 2007, Endo et al. [13] successfully applied the minimum entropy deconvolution (MED) algorithm to the fault detection of rotating motors for the first time. MED can only highlight a small number of single pulses with high amplitude and cannot extract periodic pulses. In turn, the weak components in the composite fault cannot be extracted, while the filter after the MED iterative solution is not the optimal filter. Subsequently, McDonald et al. [14] developed a method named the maximum correlated kurtosis deconvolution (MCKD) algorithm, using the maximum value of correlated kurtosis as the iteration termination condition of the optimal filter. MCKD extracts periodic pulses while suppressing noise interference. He et al. [15] used MCKD to significantly enhance the periodic pulse component of the fault signal, which allowed the fault characteristics to be more prominent. Yang et al. [16] utilized MCKD as a pre-processing operation to highlight the continuous pulse component of the bearing. The noise interference is reduced while effectively improving the ability to represent the fault characteristics. However, the performance of the MCKD algorithm depends on three parameters: filter length, fault period, and shift number [17]. In addition, the fault period also needs to be rounded after resampling if it is not an integer. To overcome the shortcomings of MED and MCKD, McDonald and Zhao proposed Multipoint optimal minimum entropy deconvolution adjusted (MOMEDA) [18]. MOMEDA determines the period of fault occurrence by multipoint cliff values. At the same time, a non-iterative approach to obtaining the optimal filter enables the analysis of a non-integer number of fault cycles without resampling. Ma and Feng [19] redesigned the objective function of the MOMEDA algorithm based on the planetary bearing vibration signal characteristics and verified the effectiveness of the proposed method by numerical simulation and experimental analysis. Wang et al. [20] improved the ability of the MOMEDA algorithm to capture fault features by constructing an autoregressive mean shift model to improve noise immunity. Xiang et al. [21] combined MOMEDA and 1.5-dimensional Teager kurtosis spectrum analysis to effectively achieve feature extraction of composite bearing faults. Due to the coupling relationship between different fault characteristics, MOMEDA alone is not immune and thus leads to diagnostic failure when dealing with multiple faults.

In 2011, Selesnick proposed resonance-based sparse signal decomposition (RSSD) [22]. It is different from the traditional signal decomposition methods based on band division. The quality factor $Q$ and redundancy degree $r$ are flexibly selected to determine the basis function, which can effectively extract the periodic pulse characteristics in the fault signal without losing important information by waveform distortion. It solves the difficulty of separating fault features in traditional methods due to the similarity of decomposition frequencies. RSSD obtains a bank of base functions for high- and low-quality factors by the Tunable Q-Factor Wavelet Transform (TQWT) [23] according to the differences in oscillation properties of the signals (i.e., differences in the quality factor $Q$). Morphological Component Analysis (MCA) [24] is employed to decouple the signal into high and low resonance components with different quality factors. Thus, it realizes the separation effectively of the different elements in the unsteady signal. Traditional resonant sparse decomposition methods in which the quality factor $Q$ and the redundancy factor $r$ are defined by human selection [25,26] make the parameters more subjective and contingent. It will lead to the corresponding basis functions not being optimally matched to the transient shock components and harmonic elements, further reducing the decoupling effect. To address this issue, a parallel parameter optimized RSSD approach is proposed in this article. The compound indicator $KHE$ serves as the objective function, such that the algorithm adaptively determines the values of the quality factor $Q$ and the redundancy degree $r$ based on the signal properties.

Motivated by the above discussion, a novel approach to wind turbine gearbox composite fault diagnosis based on MOMEDA and parallel parameter optimized RSSD is proposed in this article. Firstly, MOMEDA is used to deconvolute the signal to extract multiple periodic faults in the composite fault vibration signal. While achieving decoupling separation of the multi-fault, it effectively eliminates the influence of the transmission channels and external excitation sources. Secondly, the RSSD with parallel parameter optimized constructs the wavelet basis function bank to match the fault characteristics, which enables the interference components to be efficiently suppressed in the decoupled fault signal and the weak fault pulse to be enhanced. The superiority and effectiveness of the proposed approach are verified using the measured signals of the wind turbine gearbox.

The main structure of this article is composed as follows: Section 2 describes the theories of MOMEDA. Section 3 introduces the parallel parameter optimized RSSD based on WOA. Section 4 presents the detailed steps for extracting multi-fault characteristics with the proposed method. In Section 5, the actual composite fault signals are collected from wind turbine gearboxes to verify the proposed method and introduce a quantitative index to evaluate the performance between our approach and two comparison approaches to prove the superiority of the method in this article. Finally, the conclusions of this study are summarized in Section 6.

2. Multipoint Optimal Minimum Entropy Deconvolution Adjusted

MOMEDA utilizes a target vector to define the position and weight of the pulse sequence to be solved and applies multi-point kurtosis values to determine the period of fault occurrence. Multiple pulse target identification and deconvolution algorithms are implemented at determined locations to obtain continuous periodic pulse components.

When a rotating machine fails, the impulse signal $x$ is modulated to $s$ by the system transmission path response $h$, and together with the noise $q$ is collected by the sensor to form the vibration signal $y$. The process is expressed as:

$$y=s+q=h\ast x+q$$

(1)

The impulse signal x is recovered from the vibration signal y by an optimal filter. The MOMEDA solving process for the optimal filter can transform into a search for the maximum value of the multipoint D-parameter, using the multipoint D-parameter to reflect the shock characteristics of the filtered signal, and the related expressions are defined as follows:

$$\mathrm{MDN}(\overrightarrow{y},\overrightarrow{t})=\frac{1}{\Vert \overrightarrow{t}\Vert }\frac{{\overrightarrow{t}}^T\overrightarrow{y}}{\Vert \overrightarrow{y}\Vert }$$

(2)

$$\mathrm{MOMEDA}:\underset{\overrightarrow{f}}{\max }MDN(\overrightarrow{y},\overrightarrow{t})=\underset{\overrightarrow{f}}{\max }\frac{{\overrightarrow{t}}^T\overrightarrow{y}}{\Vert \overrightarrow{y}\Vert }$$

(3)

where the target vector $\overrightarrow{t}$ is a constant that defines the position and weight of the pulse sequence to be deconvoluted; $\overrightarrow{f}$ denotes the filter coefficients.

The extreme value of Equation (3) is obtained by taking the derivative of the filter coefficient ($\overrightarrow{f}=f_1,f_2,\dots ,f_L$).

$$\frac{d}{d\overrightarrow{f}}(\frac{{\overrightarrow{t}}^T\overrightarrow{y}}{\Vert \overrightarrow{y}\Vert })=\frac{d}{d\overrightarrow{f}}\frac{t_1{y}_1}{\Vert \overrightarrow{y}\Vert }+\frac{d}{d\overrightarrow{f}}\frac{t_2{y}_2}{\Vert \overrightarrow{y}\Vert }+\dots +\frac{d}{d\overrightarrow{f}}\frac{t_{N-L}{y}_{N-L}}{\Vert \overrightarrow{y}\Vert }$$

(4)

where, $\frac{d}{d\overrightarrow{f}}\frac{t_k{y}_k}{\Vert \overrightarrow{y}\Vert }=\Vert \overrightarrow{y}{\Vert }^{-1}{t}_k{\overrightarrow{M}}_k-\Vert \overrightarrow{y}{\Vert }^{-3}{t}_k{y}_k{X}_0\overrightarrow{y}$, ${\overrightarrow{M}}_k=\left[\begin{array}{c}{x}_{k+L-1}\\ x_{k+L-2}\\ ⋮\\ x_k\end{array}\right]$. Hence, Equation (4) can be written in the following form:

$$\frac{d}{d\overrightarrow{f}}(\frac{{\overrightarrow{t}}^T\overrightarrow{y}}{\Vert \overrightarrow{y}\Vert })=\Vert \overrightarrow{y}{\Vert }^{-1}(t_1{\overrightarrow{M}}_1+t_2{\overrightarrow{M}}_2+\dots +t_{N-L}{\overrightarrow{M}}_{N-L})-\Vert \overrightarrow{y}{\Vert }^{-3}{\overrightarrow{t}}^T\overrightarrow{y}{X}_0\overrightarrow{y}$$

(5)

It can be further simplified as follows:

$$t_1{\overrightarrow{M}}_1+t_2{\overrightarrow{M}}_2+\dots +t_{N-L}{\overrightarrow{M}}_{N-L}=X_0\overrightarrow{t}$$

(6)

Solving the extreme value with Equation (5) to zero, we can obtain:

$$\begin{array}{c}\Vert \overrightarrow{y}{\Vert }^{-1}{X}_0\overrightarrow{t}-\Vert \overrightarrow{y}{\Vert }^{-3}{\overrightarrow{t}}^T\overrightarrow{y}{X}_0\overrightarrow{y}=0\\ \frac{{\overrightarrow{t}}^T\overrightarrow{y}}{\Vert \overrightarrow{y}{\Vert }^2}{X}_0\overrightarrow{y}=X_0\overrightarrow{t}\end{array}$$

(7)

Owing to $\overrightarrow{y}={X_0}^T\overrightarrow{t}$ and the assumed existence of ${(X_0{X_0}^T)}^{-1}$, that:

$$\frac{{\overrightarrow{t}}^T\overrightarrow{y}}{\Vert \overrightarrow{y}{\Vert }^2}\overrightarrow{f}={(X_0{X_0}^T)}^{-1}{X}_0\overrightarrow{t}$$

(8)

Since the multiples of $\overrightarrow{f}$ are also the solutions of Equation (9). So, the multiples of $\overrightarrow{f}={(X_0{X_0}^T)}^{-1}{X}_0\overrightarrow{t}$ are also the solutions of MOMEDA. The advantage of this solution is to avoid iterative operations, Thus, the MOMEDA filter and output solution can be summarized as follows:

$$\overrightarrow{f}={(X_0{X_0}^T)}^{-1}{X}_0\overrightarrow{t}$$

(9)

$$X_0={\left[\begin{array}{cccccc}{x}_L& x_{L+1}& x_{L+2}& \cdots & \cdots & x_N\\ x_{L-1}& x_L& x_{L+1}& \cdots & \cdots & x_{N-1}\\ x_{L-2}& x_{L-1}& x_L& \cdots & \cdots & x_{N-2}\\ ⋮& ⋮& ⋮& \ddots & \cdots & ⋮\\ x_1& x_2& x_3& \cdots & \cdots & x_{N-L+1}\end{array}\right]}_{L \mathrm{by} N-L+1}$$

(10)

$$\overrightarrow{y}={X_0}^T\overrightarrow{f}$$

(11)

When performing target location and fault detection with MOMEDA, a chain of impulses with fault period $T$ as the step is added to the target vector.

$$\begin{array}{c}{t}_n=P_n(T)={\delta }_{round(T)}+{\delta }_{round(2T)}+{\delta }_{round(3T)}+\dots ,\\ \overrightarrow{t}=\overrightarrow{P}(T)\end{array}$$

(12)

where, $\delta$ denotes the pulse of sample $n$. The non-integer period $T$ should be rounded to the integer value closest to the fault pulse.

MOMEDA introduces Multipoint Kurtosis (MK) as a measure of fault characteristics based on kurtosis. When the output result $\overrightarrow{y}$ matches the multiple of $\overrightarrow{t}$, we obtain the standardized MK equation as follows:

$$\mathrm{Multipoint} \mathrm{Kurtosis}=\frac{{({\displaystyle {\sum }_{n=1}^{N-L}}{t}_n^2)}^2}{{\displaystyle {\sum }_{n=1}^{N-L}}{t}_n^8}\frac{{\displaystyle {\sum }_{n=1}^{N-L}}{(t_n{y}_n)}^4}{{\displaystyle {\sum }_{n=1}^{N-L}}{(y_n^2)}^2}$$

(13)

3. The Parallel Parameter Optimized RSSD Base on WOA

3.1. Tunable Q-Factor Wavelet Transform

TQWT breaks through the disadvantage of the constant quality factor of the traditional wavelet transform and makes the selection of basis function more flexible by selecting quality factor $Q$ and redundancy factor $r$, which can better match the signals with different vibration properties. TQWT uses a dual-channel filter bank to decompose the signal into multiple scales and obtain the transform coefficients by layer-by-layer decomposition. Thus, it obtains a sparse representation of the high and low resonance components. The TQWT filter bank is clarified in Figure 1. $LPS$ and $HPS$ denote the low-pass scale transform and high-pass scale transform, respectively. $\alpha$ and $\beta$ are the low-pass scale factor and high-pass scale factor.

To achieve perfect reconstruction, the frequency response function $H_0(\omega )$ of the low-pass filter and the frequency response function $H_1(\omega )$ of the high-pass filter are defined as follows:

$$H_0(\omega )=\{\begin{array}{ll}1& \left|\omega \right|\le (1-\beta )\pi \\ \theta \left(\frac{\omega +(\beta -1)\pi }{\alpha +\beta -1}\right)& (1-\beta )\pi <\left|\omega \right|<\alpha \pi \\ 0& \alpha \pi \le \left|\omega \right|\le \pi \end{array}$$

(14)

$$H_1(\omega )=\{\begin{array}{cc}1\hfill & \hfill \left|\omega \right|\le (1-\beta )\pi \hfill \\ \theta \left(\frac{\alpha \pi -\omega }{\alpha +\beta -1}\right)\hfill & \hfill (1-\beta )\pi <\left|\omega \right|<\alpha \pi \hfill \\ 0\hfill & \hfill \alpha \pi \le \left|\omega \right|\le \pi \hfill \end{array}$$

(15)

where, $\theta (\omega )=0.5(1+\cos \omega )\sqrt{2-\cos \omega }$, $\left|\omega \right|\le \pi$ $0<\alpha <1$, $0<\beta \le 1$, $\alpha +\beta >1$. Equation (16) shows the relationship between the scale transformation factor $(\alpha ,\beta )$ and the TQWT parameters $(Q,r)$, and Equation (17) provides the required number of decomposition layers.

$$\beta =\frac{2}{Q+1},\alpha =1-\frac{\beta }{r}$$

(16)

$$J_{\max }=\lfloor \frac{\log (N/4(Q+1))}{\log ((Q+1)/(Q+1-2)/r)}\rfloor$$

(17)

where: $N$ is the signal length and ⎣•⎦ is the downward rounding sign.

The high resonance components $x_1$ and low resonance components $x_2$ are extracted from the original vibration signal, and the relevant expressions are as follows:

$$y=x_1+x_2$$

(18)

Achieving sparse representation of different components of the signal by MCA. The objective function of extracting the different components is expressed as follows:

$$\begin{array}{l}\underset{{\omega }_1,{\omega }_2}{\arg \min } {\displaystyle {\sum }_{j=1}^{j_1+1}{\lambda }_{1,j}}{\Vert {\omega }_{1,j}\Vert }_1+{\displaystyle {\sum }_{j=1}^{j_2+1}{\lambda }_{2,j}}{\Vert {\omega }_{2,j}\Vert }_1\\ x={\mathrm{TQWT}}_1^{-1}({\omega }_1)+{\mathrm{TQWT}}_2^{-1}({\omega }_2)\end{array}$$

(19)

where ${\omega }_{i,j}$ is the sub-band $j$ of $\mathrm{TQWT}(i=1,2)$. ${\lambda }_1$, ${\lambda }_2$ is the regular term parameters.

However, the vibration signal inevitably has interference by background noise in actual works, and the signal separation will be transformed from Equation (19) as follow:

$$y=x_1+x_2+x_{noise}$$

(20)

$$\underset{{\omega }_1,{\omega }_2}{\arg \min }{\Vert y-{\Phi }_1{\omega }_1-{\Phi }_2{\omega }_2\Vert }_2^2+{\displaystyle {\sum }_{j=1}^{j_1+1}{\lambda }_{1,j}}{\Vert {\omega }_{1,j}\Vert }_1+{\displaystyle {\sum }_{j=1}^{j_2+1}{\lambda }_{2,j}}{\Vert {\omega }_{2,j}\Vert }_1$$

(21)

where ${\Phi }_1$ and ${\Phi }_2$ denote the inverse wavelet transform of high Q-factors and low Q-factors, respectively. Then, the minimum value of the above Equation (21) is solved to obtain the high and low resonance components according to Equation (22).

$$x_n={\mathrm{TQWT}}^{-1}({\omega }_n),\hspace{1em}n=1,2$$

(22)

3.2. Whale Optimization Algorithm

Whale Optimization Algorithm (WOA) is a novel intelligent optimization algorithm proposed by Mirjalili in 2016 [27]. The optimal solution is solved via imitating the predatory behavior of whales. The algorithm consists of three processes.

(1) Surrounding the prey

Assuming that the population size is $NP$, and the dimension is $D$, the position of the $i$th whale in the $D$-dimensional space is $X_i=(x_i^1,x_i^2,\dots x_i^D)$, $i=1,2,\dots NP$. The algorithm supposes that the desired problem is the optimal position of the whale, and the other whales will be updated and adjusted oriented to this optimal position. The mathematical model is expressed as follows:

$$X(t+1)=X^{*}(t)-A·D$$

(23)

$$D=\left|C·X^{*}(t)-X(t)\right|$$

(24)

$$A=2a·r_1-a$$

(25)

$$C=2·r$$

(26)

$$a=2(1-t/T_{\max })$$

(27)

where $D$ denotes the distance between the individual position and the optimal position; $A$ is the convergence factor. $X^{*}(t)$ is the optimal position at $t$ iterations, $X(t)$ is the individual position at $t$ iterations, $r_1$ and $r_2$ are random numbers between [0, 1], $T_{\max }$ is the maximum number of iterations, and $a$ decreases linearly from 2 to 0 with the number of iterations.

(2) Bubbling net attack

Shrinkage encircling: The convergence factor $A$ decreases linearly with $a$ to achieve position update. When the value of $A$ is taken between [−1, 1], the position of each individual is between $X(t)$ and $X^{*}(t)$ as a means to achieve prey bracketing.

Spiral updating position: whales spit out bubbles of different sizes for feeding while swimming towards the best position in a spiral posture. The mathematical model is

$$X(t+1)=D^{\prime }·e^{bl}·\cos (2\pi l)+X^{*}(t)$$

(28)

where: $D^{\prime }$ denotes the distance between the individual of the $t$ iterations and the current optimal solution. $l$ is a random number between [−1, 1], $b$ is the spiral shape constant.

During the process of prey encirclement, the whales shrink to surround and spiral forward simultaneously, and the probability of occurrence is 50% for either mode of travel.

$$X(t+1)=\{\begin{array}{ll}{X}^{*}(t)-A·D\hfill & \mathrm{if} p<0.5\hfill \\ D^{\prime }·e^{bl}·\cos (2\pi l)+X^{*}(t)\hfill & \mathrm{if} p\ge 0.5\hfill \end{array}$$

(29)

(3) Searching for prey

The whales will stop approaching the best whale individual in this stage and instead update their position by randomly searching a large area to approach any whale individual. In this case, the value of $A$ is taken as $\left|A\right|>1$. This predation strategy will cause the current individual to deviate from the target prey but will enhance the global search ability and avoid falling into local optimum. The mathematical model as follow:

$$D=\left|C·X_{rand}(t)-X(t)\right|$$

(30)

$$X(t+1)=X_{rand}(t)-A·D$$

(31)

where $X_{rand}$ denotes the location of random individuals in the current population.

3.3. Design Objective Function

Kurtosis is a 4th order statistic that reflects the sharpness of the waveform for the random variable and is sensitive to the impulse component of the signal, which is defined as:

$$K=\frac{E(x-\mu )}{{\sigma }^4}$$

(32)

Information entropy represents the uncertainty of the source information and the randomness of the event occurrence, and its value is only related to the probability distribution of the variables. Suppose a source $X=\{x_1,x_2,\dots ,x_N\}$ consists of a discrete random variable that the probability of occurrence is $p_i=P(x_i)(i=1,2,\cdots N)$ and ${\displaystyle \sum _{i=1}^{N}p=1}$, then the information entropy of the source $X$ is expressed as:

$$H(X)=-{\displaystyle \sum _{i=1}^N{p}_i\ln p_i}$$

(33)

where $\underset{p\to 0}{\lim }p\ln (p)=0$, the more balance the distribution of variables in the source, the greater the value of information entropy. The envelope spectrum is combined with the information entropy, i.e., the envelope spectrum entropy [28]. The fault signal is treated as a signal source after envelope spectrum processing, and the frequency amplitude of each point is regarded as a variable in the signal source. The formula for calculating the envelope spectral entropy is as follows:

$$H_e=-{\displaystyle \sum _{i=1}^N{p}_k\ln p_k}$$

(34)

where $p_k$ represents the envelope spectrum of the vibration signal. The envelope spectrum entropy can measure the uniformity of frequency distribution in the envelope signal and can express the complexity of the signal in the envelope domain.

In this article, we combine kurtosis and envelope spectral entropy, kurtosis as an indicator of time-domain feature can describe the impulsiveness of the signal, and envelope spectral entropy as an indicator of frequency-domain feature can represent the strength of periodic pulses. A composite indicator is constructed to reflect the time-frequency domain properties and the expression is as follows:

$$KHE=\frac{K}{H_e}$$

(35)

This indicator possesses the advantages of kurtosis and envelope spectral entropy, and it can measure the impulsivity and periodicity of the signal at the same time. The more prominent the impulsivity and periodicity of the signal, the larger the value of the indicator.

Overall, the parallel parameter optimized RSSD based on WOA as proposed in this article implements the process as follows:

(1)

The parameters of the algorithm are determined: population size $NP$, population dimension $Dm$, and the maximum number of iterations $T_{\max }$. For resonant sparse decomposition, it is required to find the optimal four parameters: $Q_{\mathrm{H}}$, $Q_{\mathrm{L}}$, $r_{\mathrm{H}}$, $r_{\mathrm{L}}$, so the dimension is set to $Dm$ = 4, the population size set to $NP$ = 30, $T_{\max }$ = 50.

(2)

Population initialization: The optimal parameters should be bounded, and the correlation between quality factors should be as low as possible. The value range of $Q_1$ takes as [8,15], the value range of $Q_2$ takes as [1,3], and the value range of redundancy factors $r_1$ and $r_2$ take as [2,5]. Then, to reduce the calculation, the accuracy of the four parameters is reserved to one single decimal.

(3)

The objective function value is calculated: the composite index constructed by kurtosis and envelope spectral entropy serves as the objective function. The objective function value of an individual is calculated and the current optimal individual is determined.

(4)

The main loop of the algorithm is entered: if $p<0.5$ and $\left|A\right|<1$, the individual updates the current position by Equation (23), otherwise the individual position is updated by Equation (31). When $p\ge 0.5$, the position is updated according to Equation (29).

(5)

Evaluating the whole whale population and iterative optimization until the algorithm converges, it obtains the optimal objective function value $KHE$. Obtain the RSSD parameters after parallel optimization: $Q_{\mathrm{H}}$, $Q_{\mathrm{L}}$, $r_{\mathrm{H}}$, $r_{\mathrm{L}}$.

4. The Procedure of Compound Fault Diagnosis

The method proposed in this article is suitable to separate and extract the compound faults of the gear faults and the bearing faults from the wind turbine gearbox. Firstly, the input vibration signal is pre-processed, the deconvolution period is set according to the fault frequency of the damaged part, and the vibration signal is decoupled into a single fault by MOMEDA. Secondly, the low resonance component is decomposed from the pre-processed signal with optimized RSSD. Finally, the envelope analysis of the low-resonance components is applied to extract the fault characteristic frequency. The flowchart of the method is shown in Figure 2.

5. Application of Proposed Method

5.1. Experiment Introduction

In this paper, 750 kW wind turbine gearbox data provided by The National Renewable Energy Laboratory (NREL) [29] were used to verify the performance of the proposed method. The gearbox first finished run-in in the NREL dynamometer facility and then was sent to the wind plant for testing under actual operating conditions. The gearbox high-speed shaft is operated at 1800 rpm during the test, and the sensors collect vibration acceleration signals with a sampling frequency of 40 kHz. The gearbox consists of a low-speed stage planetary gearbox and a parallel shaft gearbox and has an overall transmission ratio of 1:81.491. Figure 3 gives the internal layout and nomenclature abbreviations of the gearbox.

The gearbox experienced two oil loss events during the actual test, where it damaged its internal bearings and gear components. The damaged gearbox was disassembled for fault analysis, and it was found that the large and small gears of the high-speed gear pair in the gearbox were seriously scuffed, and the inner ring races of the bearings and the two ends of the rolling bodies were overheated. The failure of each component is shown in Figure 4.

The configuration parameters of the gearbox are listed in

Table 1. Parallel shaft gearbox parameters.

Gear Element	No. of Teeth	Root Diameter (mm)	Helix Angle	Face Width (mm)
Intermediate gear	82	678	14R	170
Intermediate pinion	23	174	14L	186
High-speed gear	88	440	14R	110
High-speed pinion	22	100	14L	120

. The experiment takes the HSS pinion failure, the IMS large gear, and the HSS bearing C as the research targets. The failed bearing for this experiment is SKF 3222 J2, and

Table 2. High-speed shaft bearing parameters.

Bearing Pitch diameter/mm	Large End of Rolling diameter/mm	Small End of Rolling diameter/mm	Number of rollers/N	Contact Angle (α/degree)
155.00	24.22	19.03	20	11.63

describes its structural parameters.

Table 3. Characteristic frequency of faulty parts.

HSS Pinion Frequency ${\mathit{f}}_{\mathit{r}1}$	IMS Gear Frequency ${\mathit{f}}_{\mathit{r}2}$	Meshing Frequency ${\mathit{f}}_{\mathit{m}}$	Bearing Inner Ring Fault ${\mathit{f}}_{\mathit{i}}$	Bearing Rollers Fault ${\mathit{f}}_{\mathit{b}}$
30.00	7.50	660.00	345.30	93.51

provides the results of the failure frequency calculations for the damaged parts.

5.2. Experiment Analysis

Figure 5 shows the results of the analysis of the raw vibration signal of the gearbox compound fault. In Figure 5a, the pulse component triggered by faults has been completely submerged in the noise and failed to detect the periodical components. As seen from Figure 5b, the frequency components are mostly concentrated on the meshing frequency $f_m$ and its multiplier $2f_m$ around, and the rotational frequency component $f_{r1}$ is relatively weak. In the amplified plots of the original signal spectrum in Figure 5c,d, the fault characteristic information can see: centered on the meshing frequencies $f_m$ and $2f_m$, with $f_m\pm f_{r1}$, $2f_m\pm f_{r1}$, and $2f_m+2f_{r1}$ as the sidebands. Hence, it can be concluded that the HSS pinion is faulty. In Figure 5c, the fault feature component centered on the meshing frequency $f_m$ and with $f_m\pm f_{r2}$ as the sidebands exist, but the spectrum is non-obvious. Meanwhile, there is no fond of the feature components associated with the IMS gear in Figure 5d. Therefore, it is considered a possibility that the fault is in the IMS gear. The feature components reflect the gear failures and bearing failures uncovered in the envelope spectrum shown in Figure 5e.

The deconvolution period $T_{r1}$ is set with the fault characteristics of the HSS pinion. The RSSD optimized convergence curve is shown in Figure 6. When the number of iterations reaches 14, the algorithm converges to 0.5755. The parameters obtained by optimization are: $Q_{\mathrm{H}}$ = 4.1, $Q_{\mathrm{L}}$ = 1.0, $r_{\mathrm{H}}$ = 5.0, $r_{\mathrm{L}}$ = 3.5. The rest of the optimization process is similar and will not be listed later. The periodic pulses are visible in Figure 7a, and the remarkable fault frequency characteristics with $f_{r1}$ as the interval can be seen in Figure 7b, but $f_{r1}$, $3f_{r1}$, and $4f_{r1}$ fault frequencies are relatively low. The low resonance components obtained by substituting the optimized parameters into RSSD are shown in Figure 7c. we can see that the periodic pulse characteristics are enhanced and the interference is suppressed effectively. From Figure 7d, we can see that the $f_{r1}$–$5f_{r1}$ and $6f_{r1}$–$15f_{r1}$ failure frequencies provide evidence of degradation and decay of the high-speed shaft pinion, thus diagnosing the presence of the HSS pinion faulty.

The deconvolution period $T_{r2}$ is set with the fault characteristics of the IMS gear, and the results as shown in Figure 8. Periodic pulses can be seen in Figure 8a, and the rest of the fault characteristics are not apparent except for $4f_{r2}$ in Figure 8b, but $4f_{r2}$ is easily confused with the HSS pinion fault characteristic frequency $f_{r1}$, leading to misdiagnosis. The optimized parameters are: $Q_{\mathrm{H}}$ = 4.0, $Q_{\mathrm{L}}$ = 1.2, $r_{\mathrm{H}}$ = 5.0, $r_{\mathrm{L}}$ = 3.5, and the low-resonance components obtained by decomposition are shown in Figure 8c,d. The interference components are reduced in the time domain diagram, and the IMS gear fault characteristic frequency $f_{r2}$ and its multiples $2f_{r2}$–$6f_{r2}$ predominate in the envelope spectrum. The resonance decomposition suppresses the interference obviously, and further extracts and enhances the weak fault features.

The deconvolution period $T_i$ is set with the fault characteristic frequency of the bearing inner ring, and the results as shown in Figure 9. From Figure 9b, we can see serious interference around the fault signature frequencies $3f_i$, $4f_i$, and $5f_i$, which is unable to identify the fault signature. The optimized parameters are $Q_{\mathrm{H}}$ = 14.2, $Q_{\mathrm{L}}$ = 2.7, $r_{\mathrm{H}}$ = 3.5, $r_{\mathrm{L}}$ = 3.5, and the low-resonance components obtained by decomposition as shown in Figure 9c,d. The periodicity of the fault signal in the time domain is enhanced, the feature frequency $f_i$ of the bearing inner ring fault and its multiplication frequency $2f_i$–$5f_i$ in the envelope spectrum became more significant, and the interference components become well suppressed.

The deconvolution period $T_b$ is set with the fault characteristic frequency of the bearing rollers, and the results as shown in Figure 10. The envelope spectrum in Figure 10b shows that the fault feature frequencies $2f_b$, $3f_b$, and $4f_b$ are not apparent and the bearing fault characteristics cannot be identified. The optimized parameters are $Q_{\mathrm{H}}$ = 8.0, $Q_{\mathrm{L}}$ = 2.0, $r_{\mathrm{H}}$ = 3.5, $r_{\mathrm{L}}$ = 3.5, and the low-resonance components obtained by decomposition as shown in Figure 10c,d. Though some impulse components are missing in the time domain diagram, the overall periodic characteristics are remarkably improved, and the interference components effectively removed. The high-speed grade bearing rollers fault characteristic frequency $f_b$ and its multiples $2f_b$–$10f_b$ dominate significantly in the envelope spectrum, and the interference component restraining effect is obvious.

In summary, the analysis results show that the proposed method can not only successfully decouple and separate various types of fault characteristics from the composite fault signal, but also has a remarkable suppression effect for interferences, further highlighting the periodicity and impulsiveness of fault characteristics, which helps to improve the accuracy and reliability of wind turbine gearbox composite fault diagnosis.

5.3. Comparative Analysis

The superiority and effectiveness of the proposed method are verified through a comparative analysis using the two methods. The MCKD algorithm in Reference [14] and the improved MCKD algorithm in Reference [16] are used to analyze the composite fault signals.

The two methods mentioned above are used to decouple the HSS pinion fault from the composite fault signal and analyze and process it. The parameters are set as follows: number M = 5 is shifted, length L = 500 is filtered, the above calculation results of the fault frequency are introduced, and the analysis results are shown in Figure 11. From the time domain plot in Figure 11a, the MCKD algorithm can extract only a limited set of pulses because the result after deconvolution is not an iterative optimal solution. The resonant sparse decomposition was used for further processing, the RSSD optimized convergence curve is shown in Figure 12. When the number of iterations reach 5, the algorithm converges to 1.434. The parameters obtained by optimization are: $Q_{\mathrm{H}}$ = 8.3, $Q_{\mathrm{L}}$ = 3, $r_{\mathrm{H}}$ = 5, and $r_{\mathrm{L}}$ = 3.4. From the time domain plot in Figure 11c, though it removes part of the noise interference, the continuous periodic pulse sequence cannot be extracted subject to the algorithm performance of MCKD.

The envelope spectrum analysis shown in Figure 11b, the fault characteristic frequency $f_{r1}$ and its multiples $4f_{r1}$ do not prominent in the envelope spectrum. The fault feature frequencies dominate in the frequency spectrum with the optimized resonance decomposition from Figure 11d. The fault characteristics of the HSS pinion are accurately identified.

The same procedure is applied for the IMS gear, and the parameters obtained with the optimized RSSD are: $Q_{\mathrm{H}}$ = 8.2, $Q_{\mathrm{L}}$ = 2.9, $r_{\mathrm{H}}$ = 5, and $r_{\mathrm{L}}$ = 3.5. Some of the disturbances are removed in the time domain as shown in Figure 13. However, it failed to extract the fault feature frequency of the IMS gear effectively due to the performance of MCKD.

The bearing inner ring fault signals and bearing rollers fault signals are directly processed with the optimized parameters of the method proposed in this paper, so as to imitate the process of artificially selecting parameters in traditional resonance decomposition. The results are shown in Figure 14 and Figure 15. It can be seen from the spectrum that the traditional resonance decomposition has subjectivity and randomness due to the artificial selection of parameters, which means that the basis function cannot effectively match the fault pulse components. The above comparative analysis further illustrates the validity and superiority of the proposed method in this paper.

To evaluate the performance of extracting fault features quantitatively, the fault feature coefficient (FFC) [30] is introduced as a quantitative index to select the fault components. The larger the FFC value, the more adequate the periodic pulse information contained in the fault frequency component. The FFC is defined as follows:

$$FFC=\frac{{\displaystyle \sum {[A(f)]}^2}}{{\displaystyle \sum {[A(f)]}^2}+{\displaystyle \sum {[A(f^{\prime })]}^2}}\times 100\%$$

(36)

where $f$ and $f^{\prime }$ denote the fault component and noise component, respectively, $A(f)$ and $A(f^{\prime })$ denote their amplitudes in the frequency spectrum, respectively. The denominator part is the sum of the envelope spectral amplitude for the time domain signal, and the numerator part is the sum of the envelope spectral amplitude for the fault characteristic frequency and its multiples.

Table 4. Performance indexes of the method in this paper and the comparative methods.

Value	MCKD	MCKD + RSSD	MOMEDA	Proposed Method
HSS pinion fault for FFC	9.34%	15.55%	11.10%	29.76%
IMS gear fault for FFC	27.55%	48.12%	56.24%	65.47%
Bearing inner ring fault for FFC	0.68%	1.35%	1.64%	5.83%
Bearing rollers fault for FFC	1.20%	4.17%	3.40%	11.83%

shows the FFC values of the four methods, and the results demonstrate the superiority of the proposed method in this study.

Table 5. Comparison of the average running time of different methods.

Method	MCKD	MCKD + RSSD	MOMEDA	Proposed Method
Time(s)	23.69	450.48	4.58	280.85

shows the average CPU times of the four methods with 10 tests.

6. Conclusions

A compound fault diagnosis method for wind turbine gearboxes based on MOMEDA and the parallel parameter optimized RSSD is proposed in this study. MOMEDA obtains the deconvolution period based on the fault frequency and obtains periodic continuous pulses in a non-iterative deconvolution manner, thus decoupling and separating the compound fault vibration signals. However, some weak fault features are easily buried in transmission channels and background noise and using MOMEDA alone is not immune to dealing with multiple faults. Therefore, its combination with RSSD for parallel parameter optimization is applied to suppress disturbances and enhance the relevant fault characteristics. The parallel parameter optimized RSSD takes the composite indicator with low resonance components as the objective function and adaptively obtains the best quality factor $Q$ and redundancy $r$ according to the signal properties. The composite indicator measures the periodicity of the signal while measuring impulsivity, thus improving the signal sparse representation performance of RSSD. A bank of base functions matching the fault characteristics with the selected optimal parameters was constructed. RSSD adaptively decomposes the fault signal into high-resonance components and low-resonance components. The feasibility of the proposed method is verified by actual fault signals of wind turbine gearboxes. Compared with MCKD and MCKD-RSSD methods, the proposed method not only possesses excellent performance of decoupling separation and enhancing weak fault characteristics but also clearly and accurately portrays the time-frequency domain characteristics of different fault types, which is more suitable for composite fault diagnosis of wind turbine gearboxes. However, the method also has some drawbacks, such as time-consuming computation time during parameter optimization. In further studies, we will improve this problem to reduce the computation time.