# Gearbox Fault Diagnosis in a Wind Turbine Using Single Sensor Based Blind Source Separation

Authors: Yuning Qian; Ruqiang Yan
Journal: Journal of Sensors
 (2016)

---

## Abstract

This paper presents a single sensor based blind source separation approach, namely, the wavelet-assisted stationary subspace analysis (WSSA), for gearbox fault diagnosis in a wind turbine. Continuous wavelet transform (CWT) is used as a preprocessing tool to decompose a single sensor measurement data into a set of wavelet coefficients to meet the multidimensional requirement of the stationary subspace analysis (SSA). The SSA is a blind source separation technique that can separate the multidimensional signals into stationary and nonstationary source components without the need for independency and prior information of the source signals. After that, the separated nonstationary source component with the maximum kurtosis value is analyzed by the enveloping spectral analysis to identify potential fault-related characteristic frequencies. Case studies performed on a wind turbine gearbox test system verify the effectiveness of the WSSA approach and indicate that it outperforms independent component analysis (ICA) and empirical mode decomposition (EMD), as well as the spectral-kurtosis-based enveloping, for wind turbine gearbox fault diagnosis.


---

## Body

## 1. Introduction

Wind energy is one of the renewable energy sources, which have received considerable attention around the world. As critical equipment for wind energy development, the wind turbine is being widely used due to its technological maturity and good infrastructure construction [1, 2]. However, with the high demand for energy efficiency, the size of the wind turbine is increasing over the years, leading to high operation costs and risk of failures. To avoid expensive maintenance costs and economic losses caused by wind turbine failures, researches on condition monitoring and fault diagnosis of the wind turbine have been carried out. For example, Kotzalas and Doll [3] investigated the failure principle of various wind turbines. Ciang et al. [1], García Márquez et al. [4], and Liu et al. [5] reviewed the existing condition-monitoring strategies and relative fault diagnosis methods for wind turbines. Chen et al. [6] summarized main failure modes of the wind turbines and corresponding signal characteristics. It is known that the wind turbine failures often occur in the generator, blade, gearbox, and electrical system. As a key component for wind turbines, the gearbox is vulnerable to damage and the gearbox failure has a severe influence on working status of the whole system. As a result, effective diagnosis of gearbox failures has become a focused research trend. Generally, the gear mesh signal is strong while the gearbox fault-related signal is often weak and transient, causing difficulty to separate such fault-related signatures from gear mesh signals.Over the past, some filtering methods [7, 8], for example, spectral-kurtosis-based and autoregressive model-based methods, were used to extract the fault-related signal components, but the filtering results are highly related to the filter parameter selection and may be sensitive to noise. Li et al. [9] developed a noise-controlled technique based on stochastic resonance method to enhance the fault signal by adjusting the input signal and the noise level. Combet and Gelman [10] and Heyns et al. [11] both utilized the time synchronous averaging technique to remove the gear meshing frequency with the help of an independent encoder signal to resample the measured vibration signals. Different from the above methods, independent component analysis (ICA) is a blind source separation approach, which can separate the raw signal into several independent sources without any prior information or reference signal. For example, both He et al. [12] and Wang et al. [13] used the ICA to extract the fault-related features from gearbox vibration signals. However, the ICA approach needs to assume that each component to be separated is independent and does not take distribution changes of the signals into account. Compared with ICA, the recently developed stationary subspace analysis (SSA) [14] is also a blind source separation approach. The SSA can not only take the signal distribution into consideration [15], but also decompose a multidimensional time series into stationary and nonstationary source components without the independency assumption of these source signals. Since the gearbox fault-related signal often exhibits nonstationary behavior due to complex working condition of the wind turbines, it can be naturally characterized by the nonstationary components resulting from the SSA. Therefore, the SSA presents a good alternative for gearbox fault diagnosis. The applications of the SSA have been seen in change-point detection of the multidimensional time series [16], EEG signal processing [15, 17], geophysical data analysis [18], and image processing [19]. It should be noted that there is a multidimensional requirement when implementing the SSA algorithm. For real-time gearbox fault diagnosis of the wind turbine, although there are several sensors attached on the system, we often have only one sensor near the fault component to collect a single-dimensional signal that is highly related to the fault feature. To address such a problem, this paper presents a single sensor based blind source separation approach for gearbox fault diagnosis by integrating continuous wavelet transform with stationary subspace analysis, which is called the wavelet-assisted stationary subspace analysis (WSSA). The proposed approach first uses CWT to decompose a one-dimensional signal into multiscale wavelet coefficients, which can be considered as a multidimensional signal. Then, the SSA is applied to separating it into a set of stationary and nonstationary source components. It should be emphasized that the key of the proposed approach is to utilize the inherent multiscale analysis capability and the redundancy resulting from the CWT decomposition for dimension extension of a single sensor signal [20] and then select the SSA as a blind source separation approach for both the stationary and nonstationary signal components separation. Besides, this paper also proposes to use modified principal component analysis and runs test methods to select the number of wavelet scales and the number of nonstationary source components, respectively, for further improving the effectiveness of the WSSA approach.The organization of the rest of the paper is as follows. Section2 provides the theoretical background of the CWT and the SSA. The parameter selection methods for the SSA and the framework of the WSSA-based gearbox fault diagnosis approach are introduced in Section 3. Section 4 gives the experimental results, together with some discussions. Conclusions are finally stated in Section 5.


## 2. Wavelet-Assisted Stationary Subspace Analysis

### 2.1. Continuous Wavelet Transform

The CWT of a signalx(t) is defined as the inner product of the signal and a selected wavelet function, which is expressed as(1)Wfxa,b=ψa,bt,xt=a-1/2∫xtψa,b∗tdt,where ψa,bt=a-1/2ψ((t-b)/a) is the scaled and translated wavelet function, a is the scale factor, and b denotes the time location. From this equation, the CWT at a certain scale a can be considered as a continuous correlation operation between the original signal and the wavelet function by changing the time location b. The result of the CWT is a series of wavelet coefficients with the same length as the original signal, which can express the similarity between the signal and the wavelet function at a given scale. It should be noted that the selection of the wavelet function significantly affects the performance of the CWT. Previous study for wavelet function selection using a quantitative measure (i.e., energy-to-Shannon entropy ratio) [22] has shown that the Morlet wavelet is effective for mechanical fault feature extraction and it has also been successfully applied in many cases [23–25]. Hence, the Morlet wavelet is chosen as the wavelet function in this study.


### 2.2. Stationary Subspace Analysis (SSA)

The SSA [14] is a blind source separation technique that can separate the stationary source components from the nonstationary source components in a multidimensional signal. Neither independency nor prior information is required for these source components. In the SSA algorithm, the observed signal x(t) with D-dimension is assumed to be generated by a linear mixture of d stationary sources (s-sources) ss(t)=[s1(t),…,sd(t)]T and dn   (dn=D-d) nonstationary (n-sources) sources sn(t)=[sd+1(t),…,sD(t)]T, and x(t) can be expressed as(2)xt=Ast=AsAnsstsnt,where A is an unknown invertible mixing matrix. The spaces spanned by the columns of As and An are called s-space and n-space, respectively. Particularly, different from the ICA algorithm, there is no independency assumption on the sources ss(t) and sn(t). The aim of the SSA algorithm is to find a linear transformation(3)A^-1=P^sP^nthat can separate the stationary sources ss(t) from nonstationary sources sn(t). This can be expressed as(4)s^sts^nt=A^-1xt=P^sxtP^nxt,where s^s(t) and s^n(t) are the estimated stationary and nonstationary sources and P^s and P^n are called s-projection and n-projection. The SSA algorithm uses weak stationarity condition as an optimization criterion to recover the sources as stationary as possible [16]. The weak stationarity condition can be expressed as the mean and covariance of a time series are constant over time. Since the SSA employs an optimization criterion to recover the sources as stationary as possible, it can identify the true ss(t), not the true sn(t). The sn(t) is often identified by maximizing the nonstationarity of the estimated nonstationary sources. The flowchart of the SSA is shown in Figure 1 and the specific procedures are illustrated as follows.Figure 1
Flowchart of the SSA algorithm.(I) The observed signal is divided intoN  (N≥(D-ds)/2+2) equal time epochs and the mean value μ^i and covariance Σ^i for theith epoch are calculated. For any selected P^s, the mean value μ^s,i=P^sμi and covariance matrix Σ^s,i=P^sΣ^i(P^s)T of the estimated stationary source in theith epoch can be obtained. Then the probability distribution of the estimated stationary source in theith epoch can be obtained using normal distribution as Norm(μ^s,i,Σ^s,i).(II) Kullback-Leibler (KL) divergenceDKL [26] is utilized to measure the difference between the distribution of the estimated stationary source and the standard normal distribution in each epoch. Given that ps,i(x) and qs,i(x) are probability distribution functions of two distributions Norm(μ^s,i,Σ^s,i) and Norm(0,I), the KL divergence between the two distributions in theith epoch can be defined as(5)DKLNormμ^s,i,Σ^s,i∥Norm0,I=∑x=1nips,ixlogps,ixqs,ix,where I is unit matrix and ni is the point number in theith epoch. Then the objective function is set as the sum of the KL divergence in each epoch as shown in the following equation:(6)fP^s=∑i=1NDKLNormμ^s,i,Σ^s,i∥Norm0,I=∑i=1NDKLNormP^sμ^i,P^sΣ^iP^sT∥Norm0,I=∑i=1N∑x=1nips,ixlogps,ixqs,ix.It should be noted that P^s is the only variable in (6). Therefore, the optimal s-projection P^s can be obtained by minimizing the objective function shown in (7). Then, the optimal estimated stationary source s^s(t) can be obtained using (4): (7)LP^s=minfP^s.(III) Similar to the procedure described in (I) and (II), the objective function for nonstationary source components is formulated in (8). By maximizing the objective function [16] as shown in (9), the n-projection P^n is obtained and the estimated nonstationary source s^n(t) can be calculated by (4): (8)gP^n=∑i=1NDKLNormμ^n,i,Σ^n,i∥Norm0,IgP^n=∑i=1NDKLNormP^nμ^i,P^nΣ^iP^nT∥Norm0,IgP^n=∑i=1N∑x=1nipn,ixlogpn,ixqn,ix,(9)LP^n=maxgP^n.From the procedures presented above, it can be seen that, as compared to the ICA, the SSA takes the distribution change of the raw signal into account and allows for dependency between the source components, which may better satisfy real-world situations. Combing the CWT with the SSA, a single sensor based blind source separation approach (i.e., WSSA) can now be developed.


### 2.3. Simulation Evaluation

To evaluate the performance of the WSSA for gearbox fault diagnosis, a simulation study is first conducted, where a test signal simulating the gearbox signal is formulated as follows [13]: (10)x1t=e-400∗t1sin⁡2πf1t,t1=modt,1fB,x2t=sin2πf2t,x3t=0.8∗sin2πf3t,x4t=sin2πf4t,nt=0.16∗randnn,1,n=lengtht,xt=x1t+x2t+x3t+x4t+nt.There are four source signals in the simulated signal as shown in Figure 2. x1t represents the outer raceway defect-related signal of a bearing modulated by the resonance signal of the system, where f1 is the resonance frequency and fB is the bearing defect characteristic frequency. The 3-stage gear meshing signals are simulated by x2t, x3t, and x4t, respectively. In addition, nt represents the noise. The values of these frequency components are set as fB = 33 Hz, f1 = 3600 Hz, f2 = 420 Hz, f3 = 160 Hz, and f4 = 50 Hz, respectively. The waveform of the simulated signal and its frequency spectrum are shown in Figure 3, where the three-stage gear meshing frequencies f2, f3, and f4 can be seen clearly, but the bearing defect characteristic frequency cannot be identified.Figure 2
Source components of simulation signal.Figure 3
The simulated signal and its frequency spectrum.
(a)
(b)The simulated signal is then decomposed using the CWT, where the Morlet wavelet is chosen as the wavelet function. Since there are four main source components in the raw signal, the wavelet coefficients at four scales are chosen as input for the SSA decomposition and they are separated into 3 stationary sources and 1 nonstationary source as shown in Figure4. It can be seen that the shape of the nonstationary source is similar to the wave shape of the x1(t). At last, the Hilbert transform is used to extract the envelope of the nonstationary source and spectral analysis is subsequently applied to obtaining the envelope spectrum as shown in Figure 5(b), where the bearing defect characteristic frequency fB is clearly extracted. Although the simulated signal is not very complex when comparing with the signals in real cases, the simulation results verify the effectiveness of the WSSA for gearbox failure signal analysis.Figure 4
SSA results of simulation signal.Figure 5
Envelope analysis for nonstationary source with maximum kurtosis.
(a)
(b)It should be noted that when the WSSA is applied to real gearbox vibration signals, the number of source components contained in the vibration signals is not knowna priori; thus it is necessary to determine how many wavelet scales (D value) should be chosen and which scales need to be selected for the WSSA. In addition, how to choose a proper number of stationary or nonstationary sources is another key issue. These two issues will be addressed in the following section.


## 2.1. Continuous Wavelet Transform

The CWT of a signalx(t) is defined as the inner product of the signal and a selected wavelet function, which is expressed as(1)Wfxa,b=ψa,bt,xt=a-1/2∫xtψa,b∗tdt,where ψa,bt=a-1/2ψ((t-b)/a) is the scaled and translated wavelet function, a is the scale factor, and b denotes the time location. From this equation, the CWT at a certain scale a can be considered as a continuous correlation operation between the original signal and the wavelet function by changing the time location b. The result of the CWT is a series of wavelet coefficients with the same length as the original signal, which can express the similarity between the signal and the wavelet function at a given scale. It should be noted that the selection of the wavelet function significantly affects the performance of the CWT. Previous study for wavelet function selection using a quantitative measure (i.e., energy-to-Shannon entropy ratio) [22] has shown that the Morlet wavelet is effective for mechanical fault feature extraction and it has also been successfully applied in many cases [23–25]. Hence, the Morlet wavelet is chosen as the wavelet function in this study.


## 2.2. Stationary Subspace Analysis (SSA)

The SSA [14] is a blind source separation technique that can separate the stationary source components from the nonstationary source components in a multidimensional signal. Neither independency nor prior information is required for these source components. In the SSA algorithm, the observed signal x(t) with D-dimension is assumed to be generated by a linear mixture of d stationary sources (s-sources) ss(t)=[s1(t),…,sd(t)]T and dn   (dn=D-d) nonstationary (n-sources) sources sn(t)=[sd+1(t),…,sD(t)]T, and x(t) can be expressed as(2)xt=Ast=AsAnsstsnt,where A is an unknown invertible mixing matrix. The spaces spanned by the columns of As and An are called s-space and n-space, respectively. Particularly, different from the ICA algorithm, there is no independency assumption on the sources ss(t) and sn(t). The aim of the SSA algorithm is to find a linear transformation(3)A^-1=P^sP^nthat can separate the stationary sources ss(t) from nonstationary sources sn(t). This can be expressed as(4)s^sts^nt=A^-1xt=P^sxtP^nxt,where s^s(t) and s^n(t) are the estimated stationary and nonstationary sources and P^s and P^n are called s-projection and n-projection. The SSA algorithm uses weak stationarity condition as an optimization criterion to recover the sources as stationary as possible [16]. The weak stationarity condition can be expressed as the mean and covariance of a time series are constant over time. Since the SSA employs an optimization criterion to recover the sources as stationary as possible, it can identify the true ss(t), not the true sn(t). The sn(t) is often identified by maximizing the nonstationarity of the estimated nonstationary sources. The flowchart of the SSA is shown in Figure 1 and the specific procedures are illustrated as follows.Figure 1
Flowchart of the SSA algorithm.(I) The observed signal is divided intoN  (N≥(D-ds)/2+2) equal time epochs and the mean value μ^i and covariance Σ^i for theith epoch are calculated. For any selected P^s, the mean value μ^s,i=P^sμi and covariance matrix Σ^s,i=P^sΣ^i(P^s)T of the estimated stationary source in theith epoch can be obtained. Then the probability distribution of the estimated stationary source in theith epoch can be obtained using normal distribution as Norm(μ^s,i,Σ^s,i).(II) Kullback-Leibler (KL) divergenceDKL [26] is utilized to measure the difference between the distribution of the estimated stationary source and the standard normal distribution in each epoch. Given that ps,i(x) and qs,i(x) are probability distribution functions of two distributions Norm(μ^s,i,Σ^s,i) and Norm(0,I), the KL divergence between the two distributions in theith epoch can be defined as(5)DKLNormμ^s,i,Σ^s,i∥Norm0,I=∑x=1nips,ixlogps,ixqs,ix,where I is unit matrix and ni is the point number in theith epoch. Then the objective function is set as the sum of the KL divergence in each epoch as shown in the following equation:(6)fP^s=∑i=1NDKLNormμ^s,i,Σ^s,i∥Norm0,I=∑i=1NDKLNormP^sμ^i,P^sΣ^iP^sT∥Norm0,I=∑i=1N∑x=1nips,ixlogps,ixqs,ix.It should be noted that P^s is the only variable in (6). Therefore, the optimal s-projection P^s can be obtained by minimizing the objective function shown in (7). Then, the optimal estimated stationary source s^s(t) can be obtained using (4): (7)LP^s=minfP^s.(III) Similar to the procedure described in (I) and (II), the objective function for nonstationary source components is formulated in (8). By maximizing the objective function [16] as shown in (9), the n-projection P^n is obtained and the estimated nonstationary source s^n(t) can be calculated by (4): (8)gP^n=∑i=1NDKLNormμ^n,i,Σ^n,i∥Norm0,IgP^n=∑i=1NDKLNormP^nμ^i,P^nΣ^iP^nT∥Norm0,IgP^n=∑i=1N∑x=1nipn,ixlogpn,ixqn,ix,(9)LP^n=maxgP^n.From the procedures presented above, it can be seen that, as compared to the ICA, the SSA takes the distribution change of the raw signal into account and allows for dependency between the source components, which may better satisfy real-world situations. Combing the CWT with the SSA, a single sensor based blind source separation approach (i.e., WSSA) can now be developed.


## 2.3. Simulation Evaluation

To evaluate the performance of the WSSA for gearbox fault diagnosis, a simulation study is first conducted, where a test signal simulating the gearbox signal is formulated as follows [13]: (10)x1t=e-400∗t1sin⁡2πf1t,t1=modt,1fB,x2t=sin2πf2t,x3t=0.8∗sin2πf3t,x4t=sin2πf4t,nt=0.16∗randnn,1,n=lengtht,xt=x1t+x2t+x3t+x4t+nt.There are four source signals in the simulated signal as shown in Figure 2. x1t represents the outer raceway defect-related signal of a bearing modulated by the resonance signal of the system, where f1 is the resonance frequency and fB is the bearing defect characteristic frequency. The 3-stage gear meshing signals are simulated by x2t, x3t, and x4t, respectively. In addition, nt represents the noise. The values of these frequency components are set as fB = 33 Hz, f1 = 3600 Hz, f2 = 420 Hz, f3 = 160 Hz, and f4 = 50 Hz, respectively. The waveform of the simulated signal and its frequency spectrum are shown in Figure 3, where the three-stage gear meshing frequencies f2, f3, and f4 can be seen clearly, but the bearing defect characteristic frequency cannot be identified.Figure 2
Source components of simulation signal.Figure 3
The simulated signal and its frequency spectrum.
(a)
(b)The simulated signal is then decomposed using the CWT, where the Morlet wavelet is chosen as the wavelet function. Since there are four main source components in the raw signal, the wavelet coefficients at four scales are chosen as input for the SSA decomposition and they are separated into 3 stationary sources and 1 nonstationary source as shown in Figure4. It can be seen that the shape of the nonstationary source is similar to the wave shape of the x1(t). At last, the Hilbert transform is used to extract the envelope of the nonstationary source and spectral analysis is subsequently applied to obtaining the envelope spectrum as shown in Figure 5(b), where the bearing defect characteristic frequency fB is clearly extracted. Although the simulated signal is not very complex when comparing with the signals in real cases, the simulation results verify the effectiveness of the WSSA for gearbox failure signal analysis.Figure 4
SSA results of simulation signal.Figure 5
Envelope analysis for nonstationary source with maximum kurtosis.
(a)
(b)It should be noted that when the WSSA is applied to real gearbox vibration signals, the number of source components contained in the vibration signals is not knowna priori; thus it is necessary to determine how many wavelet scales (D value) should be chosen and which scales need to be selected for the WSSA. In addition, how to choose a proper number of stationary or nonstationary sources is another key issue. These two issues will be addressed in the following section.


## 3. Framework of Gearbox Fault Diagnosis Based on WSSA

### 3.1. DimensionD Selection Using mPCA

Generally speaking, the result of the CWT is a high-dimensional wavelet coefficient matrix. The considerable signal dimension severely affects the algorithm efficiency of the WSSA. Hence, dimension reduction is necessary for the high-dimensional wavelet coefficient matrix. As a popular dimension reduction method, principal component analysis (PCA) has been widely used to transform a set of possibly correlated variables into a set of linearly uncorrelated principal components which can better reflect the characteristics of the system [27]. Specifically, given aN-dimension series x(t)=[x1(t),x2(t),…,xN(t)], PCA can transform x(t) into a set of principal components f(t)=[f1(t),f2(t),…,fN(t)]. Each fi(t) is a linear superposition of {x1(t),x2(t),…,xN(t)} as(11)f1t=b11x1t+b12x2t+⋯+b1NxNt,f2t=b21x1t+b22x2t+⋯+b2NxNt,⋮fNt=bN1x1t+bN2x2t+⋯+bNNxNt,where bij (i,j=1 to N) represents the principal component coefficient and each principal component corresponds to a variance contribution σi. Equation (11) can be rewritten using matrix style as(12)fTt=BxTt,where(13)B=b11b12⋯b1Nb21b21⋯b2N⋮⋮⋮⋮bN1bN2⋯bNN.Traditionally, the principal components with high cumulative variance contribution can be selected as low-dimensional features of the raw high-dimensional series. For example, if ∑i=1K-1σi<σth and ∑i=1Kσi≥σth (σth represents cumulative variance contribution), the K principal components f1(t) to fK(t) can be chosen as low-dimensional features. However, there will be a problem when the traditional PCA is used to choose the outstanding wavelet scales. The wavelet coefficient series in each scale represents the feature of the raw signal under a certain frequency-band window, but the chosen principal components by the PCA are a linear superposition of these wavelet coefficients in each scale, leading to multifrequency information mixture in each principal component, which may increase separation difficulty of the SSA. On the other hand, for each principal component, the value of component coefficient bij can reflect the contribution rate of each xj(t). The larger the coefficient bij is, the more important the wavelet series xj(t) will be. Taking all chosen principal components into account, a measure called variable score can be calculated using the following method: based on cumulative variance contribution σth, the number D  (D<N) of outstanding principal components is obtained. Then the coefficients of the first D principal components are processed by (14) to calculate the variable score F=[F1F2⋯FD]:(14)FT=BTΔT,where Δ=σ1σ2⋯σN. This equation can be further expanded as (15)F1=b11σ1+b21σ2+⋯+bN1σN,F2=b12σ1+b22σ2+⋯+bN2σN,⋮FD=b1Dσ1+b2Dσ2+⋯+bNDσN,where Fi stands for theith variable score. It should be noted that each Fi represents the contribution rate of wavelet coefficients in each scale xi(t). As a result, the D series of wavelet coefficients with larger score F, corresponding to D wavelet scales, are chosen from theN-dimensional wavelet series as a D-dimension feature for the WSSA algorithm.


### 3.2. Dimensiond Selection Using Runs Test

Before the SSA algorithm is executed, we need to determine the number of stationary sourcesdfrom the D-dimensional features. If d is chosen, dn can be calculated as D-d. In [16], this parameter selection issue is discussed and a likelihood ratio test statistic is constructed to choose d using the distribution test. In fact, the essence of this method is to identify the stationarity of the source signals obtained by the SSA. Hence, this d selection problem can be transferred to a stationarity identification problem. The D-dimensional signals can be decomposed by the SSA using different d from 1 to D-1, and the dimension of the most stationary sources is the most suitabled value. Therefore, a novel d selection method, called runs test, has been utilized for the SSA in this study. The runs test method [28] is a stationarity identification method which has been used in various fields, such as wireless networks [29], financial time series analysis [30], and simulation run length determination [31]. In this study, thed selection method using runs test is designed as follows.(a) Based on theD-dimensional wavelet coefficient series, d=i  (i=1 to D) is chosen to implement the SSA and obtain thed stationary sources s^s(t) and D-d nonstationary sources s^n(t).(b) For each stationary source involved in ad-dimensional stationary source, the points larger than mean value will be marked as 1 and the rest will be marked as 0. The number of 1 and 0 can be marked as n1 and n0, respectively. The runsr can be defined as the total number of consecutive 1’s and 0’s in one series. For example, given two series L1=10011 and L2=11000, the runs r1=3 and r2=2, respectively. Based on the definition, the number of runs r can be calculated.(c) The total number of runs in the signal can provide information regarding whether it is stationary or not. A signal that contains very few runs is probably stationary, while a signal that contains many runs is more likely to be nonstationary. Considering that different sizes of different signals may affect runs calculation, the runs valuer cannot be compared directly. Since the runs valuercan be treated as an approximate Gaussian distribution, the mean value and variance ofrcan be acquired as (16) [30] and a statistic Z [32], standardization of r, can be obtained to identify the stationarity of the series as (17). The larger the |Z| is, the more nonstationary the series will be. As a result, the mean value Zd of |Z| for d stationary sources is obtained from (18). Consider (16)E=2n0n1n0+n1+1,D=2n0n12n0n1-n0-n1n0+n12n0+n1-1,(17)Z=r-ED=r-2n0n1/n0+n1-12n0n12n0n1-n0-n1/n0+n12n0+n1-1,(18)Zd=1d∑i=1dZi.(d)d from 1 to D-1 is chosen for the SSA algorithm, and then (D-1)  Zd values are obtained based on procedures (b) and (c). Thedvalue, corresponding to the smallest Zd, will be selected as the final number of stationary sources.


### 3.3. Framework of Gearbox Fault Diagnosis Approach

By integrating the proposed methods, the framework of gearbox fault diagnosis approach is presented in Figure6. The gearbox vibration signal is firstly decomposed by the CWT to obtain the multiscale wavelet coefficients. Secondly, the modified PCA is used to select the coefficients in certain scales and obtain a low-dimensional coefficient series. Thirdly, the runs test method is utilized to decide how many stationary sources should be separated from the wavelet coefficient series. Fourthly, the SSA algorithm is applied to separating the wavelet coefficient series into stationary sources and nonstationary sources. At last, the nonstationary source with the maximum kurtosis value is selected and processed using enveloping spectral analysis. The obtained envelope spectrum is used to detect the defect frequency of the gearbox.Figure 6
Flowchart of gearbox fault diagnosis approach.


## 3.1. DimensionD Selection Using mPCA

Generally speaking, the result of the CWT is a high-dimensional wavelet coefficient matrix. The considerable signal dimension severely affects the algorithm efficiency of the WSSA. Hence, dimension reduction is necessary for the high-dimensional wavelet coefficient matrix. As a popular dimension reduction method, principal component analysis (PCA) has been widely used to transform a set of possibly correlated variables into a set of linearly uncorrelated principal components which can better reflect the characteristics of the system [27]. Specifically, given aN-dimension series x(t)=[x1(t),x2(t),…,xN(t)], PCA can transform x(t) into a set of principal components f(t)=[f1(t),f2(t),…,fN(t)]. Each fi(t) is a linear superposition of {x1(t),x2(t),…,xN(t)} as(11)f1t=b11x1t+b12x2t+⋯+b1NxNt,f2t=b21x1t+b22x2t+⋯+b2NxNt,⋮fNt=bN1x1t+bN2x2t+⋯+bNNxNt,where bij (i,j=1 to N) represents the principal component coefficient and each principal component corresponds to a variance contribution σi. Equation (11) can be rewritten using matrix style as(12)fTt=BxTt,where(13)B=b11b12⋯b1Nb21b21⋯b2N⋮⋮⋮⋮bN1bN2⋯bNN.Traditionally, the principal components with high cumulative variance contribution can be selected as low-dimensional features of the raw high-dimensional series. For example, if ∑i=1K-1σi<σth and ∑i=1Kσi≥σth (σth represents cumulative variance contribution), the K principal components f1(t) to fK(t) can be chosen as low-dimensional features. However, there will be a problem when the traditional PCA is used to choose the outstanding wavelet scales. The wavelet coefficient series in each scale represents the feature of the raw signal under a certain frequency-band window, but the chosen principal components by the PCA are a linear superposition of these wavelet coefficients in each scale, leading to multifrequency information mixture in each principal component, which may increase separation difficulty of the SSA. On the other hand, for each principal component, the value of component coefficient bij can reflect the contribution rate of each xj(t). The larger the coefficient bij is, the more important the wavelet series xj(t) will be. Taking all chosen principal components into account, a measure called variable score can be calculated using the following method: based on cumulative variance contribution σth, the number D  (D<N) of outstanding principal components is obtained. Then the coefficients of the first D principal components are processed by (14) to calculate the variable score F=[F1F2⋯FD]:(14)FT=BTΔT,where Δ=σ1σ2⋯σN. This equation can be further expanded as (15)F1=b11σ1+b21σ2+⋯+bN1σN,F2=b12σ1+b22σ2+⋯+bN2σN,⋮FD=b1Dσ1+b2Dσ2+⋯+bNDσN,where Fi stands for theith variable score. It should be noted that each Fi represents the contribution rate of wavelet coefficients in each scale xi(t). As a result, the D series of wavelet coefficients with larger score F, corresponding to D wavelet scales, are chosen from theN-dimensional wavelet series as a D-dimension feature for the WSSA algorithm.


## 3.2. Dimensiond Selection Using Runs Test

Before the SSA algorithm is executed, we need to determine the number of stationary sourcesdfrom the D-dimensional features. If d is chosen, dn can be calculated as D-d. In [16], this parameter selection issue is discussed and a likelihood ratio test statistic is constructed to choose d using the distribution test. In fact, the essence of this method is to identify the stationarity of the source signals obtained by the SSA. Hence, this d selection problem can be transferred to a stationarity identification problem. The D-dimensional signals can be decomposed by the SSA using different d from 1 to D-1, and the dimension of the most stationary sources is the most suitabled value. Therefore, a novel d selection method, called runs test, has been utilized for the SSA in this study. The runs test method [28] is a stationarity identification method which has been used in various fields, such as wireless networks [29], financial time series analysis [30], and simulation run length determination [31]. In this study, thed selection method using runs test is designed as follows.(a) Based on theD-dimensional wavelet coefficient series, d=i  (i=1 to D) is chosen to implement the SSA and obtain thed stationary sources s^s(t) and D-d nonstationary sources s^n(t).(b) For each stationary source involved in ad-dimensional stationary source, the points larger than mean value will be marked as 1 and the rest will be marked as 0. The number of 1 and 0 can be marked as n1 and n0, respectively. The runsr can be defined as the total number of consecutive 1’s and 0’s in one series. For example, given two series L1=10011 and L2=11000, the runs r1=3 and r2=2, respectively. Based on the definition, the number of runs r can be calculated.(c) The total number of runs in the signal can provide information regarding whether it is stationary or not. A signal that contains very few runs is probably stationary, while a signal that contains many runs is more likely to be nonstationary. Considering that different sizes of different signals may affect runs calculation, the runs valuer cannot be compared directly. Since the runs valuercan be treated as an approximate Gaussian distribution, the mean value and variance ofrcan be acquired as (16) [30] and a statistic Z [32], standardization of r, can be obtained to identify the stationarity of the series as (17). The larger the |Z| is, the more nonstationary the series will be. As a result, the mean value Zd of |Z| for d stationary sources is obtained from (18). Consider (16)E=2n0n1n0+n1+1,D=2n0n12n0n1-n0-n1n0+n12n0+n1-1,(17)Z=r-ED=r-2n0n1/n0+n1-12n0n12n0n1-n0-n1/n0+n12n0+n1-1,(18)Zd=1d∑i=1dZi.(d)d from 1 to D-1 is chosen for the SSA algorithm, and then (D-1)  Zd values are obtained based on procedures (b) and (c). Thedvalue, corresponding to the smallest Zd, will be selected as the final number of stationary sources.


## 3.3. Framework of Gearbox Fault Diagnosis Approach

By integrating the proposed methods, the framework of gearbox fault diagnosis approach is presented in Figure6. The gearbox vibration signal is firstly decomposed by the CWT to obtain the multiscale wavelet coefficients. Secondly, the modified PCA is used to select the coefficients in certain scales and obtain a low-dimensional coefficient series. Thirdly, the runs test method is utilized to decide how many stationary sources should be separated from the wavelet coefficient series. Fourthly, the SSA algorithm is applied to separating the wavelet coefficient series into stationary sources and nonstationary sources. At last, the nonstationary source with the maximum kurtosis value is selected and processed using enveloping spectral analysis. The obtained envelope spectrum is used to detect the defect frequency of the gearbox.Figure 6
Flowchart of gearbox fault diagnosis approach.


## 4. Experimental Study

In order to verify the effectiveness of the WSSA-based gearbox fault diagnosis approach, vibration signals collected from a real wind turbine gearbox are analyzed. The experimental data is contributed by a project called Gearbox Reliability Collaborative (GRC) hosted by National Renewable Energy Laboratory (NREL) [33]. The test turbine (Figure 7) is a stall-controlled, three-bladed, upwind turbine with a rated power of 750 kW. The gearbox under test is composed of one low speed (LS) planetary stage and two parallel stages and has an overall ratio of 1 : 81.491. The structure of this gearbox is shown in Figure 8. Accelerometers are mounted on several locations outside of the gearbox to collect vibration signals with 40 kHz sampling rate.Figure 7
Test wind turbine.Figure 8
Structure of test gearbox.


### 4.1. Case Study I

The first case study analyzes a vibration signal measured by the accelerometer AN6 (model: IMI622B01), which is close to bearing D (model: SKF NU-2220-ECM, near the high speed shaft). In this study, the rotation speed of the high speed shaft (denoted as HS_shaft) is 1200 rpm and the electric power is 25%. The waveform and frequency spectrum of the vibration signal are shown in Figure9. It can be seen that the rotation frequency of HS_shaft (fHS = 20 Hz) is visible in Figure 9(b), while no defect-related frequency for bearing D is identified.Figure 9
Waveform and frequency spectrum of AN6 signals: (a) waveform; (b) frequency spectrum.
(a)
(b)The WSSA algorithm is then used to process the same signal. After the CWT decomposition, wavelet coefficients of scale 1–32 are obtained and mPCA is used to find four outstanding scales, which are identified as scales 15, 18, 19, and 20. Then the runs test method is utilized to compute the statisticZ1=192.89, Z2=265.10, and Z3=217.06, among which Z1 is the smallest. This means one stationary source and three nonstationary sources can be determined through the SSA. The results for extracted stationary and nonstationary source components are shown in Figure 10, among which the nonstationary source 3 with the maximum kurtosis value is selected and processed using enveloping spectral analysis. From the envelope spectrum shown in Figure 11, the inner raceway defect frequency of the bearing D (fBPFI = 50 Hz) can be clearly identified. In fact, the bearing D indeed suffered from the inner raceway defect as shown in Figure 12. To make a comparison, the wavelet coefficients at the scale with the maximum kurtosis are directly used for enveloping spectral analysis, and the defect frequency 50 Hz cannot be identified from the results shown in Figure 13. In addition, the ICA method is used to process the same vibration signal and conduct the same diagnostic procedure as the SSA does. As shown in Figure 14, although the defect frequency can also be identified, it is not so prominent as compared to the result obtained using the SSA. Empirical mode decomposition (EMD) technique is also used to decompose the signal and the intrinsic mode function (IMF) with maximum kurtosis is chosen for enveloping spectral analysis. From the result shown in Figure 15, the defect frequency is nearly submerged in the surrounding frequency components. Furthermore, the spectral kurtosis is applied to the vibration signals to find the proper central frequency 2500 Hz and bandwidth 312 Hz (shown in Figure 16), and then a band-pass filter with the chosen central frequency and bandwidth is used to filter the vibration signal. From the envelope spectrum of the filtered signal shown in Figure 17, the magnitude of the defect frequency 50 Hz is not large enough as compared to other adjacent frequency components. This case study verifies the effectiveness of the proposed WSSA for wind turbine gearbox fault diagnosis.Figure 10
SSA results of AN6 signals.Figure 11
Envelope analysis of nonstationary source with maximum kurtosis: (a) waveform; (b) frequency spectrum.
(a)
(b)Figure 12
Inner race defect of bearing D [21].Figure 13
CWT and maximum kurtosis comparison results of AN6 signals: (a) waveform; (b) frequency spectrum.
(a)
(b)Figure 14
ICA comparison results of AN6 signals: (a) waveform; (b) frequency spectrum.
(a)
(b)Figure 15
EMD comparison results of AN6 signals: (a) waveform; (b) frequency spectrum.
(a)
(b)Figure 16
Spectral kurtosis of AN6 signals.Figure 17
Filtered results based on spectral kurtosis of AN6 signals: (a) filtered signal; (b) frequency spectrum.
(a)
(b)


### 4.2. Case Study II

The second case study analyzes a vibration signal measured by the accelerometer AN3 (model: IMI 626B02), which is close to bearing H (model: FAG SL181856E). In this study, the rotation speed of the high speed shaft (HS_shaft) is 1800 rpm and the electric power is 25%. The bearing H is used to support the main shaft (rotation speed: 22.09 rpm). The waveform and frequency spectrum of the vibration signal are shown in Figure18, where no defect-related frequency for bearing H is identified.Figure 18
Waveform and frequency spectrum of AN3 signals: (a) waveform; (b) frequency spectrum.
(a)
(b)The vibration signal is then processed by the proposed WSSA approach. Following the same procedure as shown in case study I, wavelet coefficients of the vibration signal at scales 1–32 are obtained by the CWT. Then the coefficients at the three scales 23, 24, and 25 are selected by the mPCA, and the runs test is performed to calculate the statisticZ1=409.68 and Z2=350.80, among which Z2 is smaller than Z1. This means two stationary sources and one nonstationary source can be extracted from the vibration signal through the SSA, and the result is shown in Figure 19. The enveloping spectral analysis of the nonstationary source is then conducted and the resulting envelope spectrum is shown in Figure 20. In addition to the intermediate speed shaft (denoted as IS_shaft) rotation frequency fIS = 7.5 Hz, the outer raceway defect frequency of the bearing H (fBPFO = 8.6 Hz) can also be clearly identified. This is consistent with the physical examination for the bearing H as shown in Figure 21. For comparison study, the methods using the CWT with maximum kurtosis-based scale selection, ICA, EMD, and the band-pass filter based on spectral kurtosis are also used to process the same vibration signal and the results are shown in Figures 22–26, respectively. From these figures, it can be seen that the rotation frequency of the IS_shaft is highlighted, while the outer raceway defect frequency of the bearing H is almost submerged in the spectrum and cannot be identified. The comparison results further confirm the effectiveness of the WSSA for wind turbine gearbox fault diagnosis.Figure 19
SSA result of AN3 signals.Figure 20
Envelope analysis of nonstationary source with maximum kurtosis: (a) waveform; (b) frequency spectrum.
(a)
(b)Figure 21
Outer race defect of the bearing H [21].Figure 22
CWT and maximum kurtosis comparison results of AN3 signals: (a) waveform; (b) frequency spectrum.
(a)
(b)Figure 23
ICA comparison results of AN3 signals: (a) waveform; (b) frequency spectrum.
(a)
(b)Figure 24
EMD comparison results of AN3 signals: (a) waveform; (b) frequency spectrum.
(a)
(b)Figure 25
Spectral kurtosis of AN3 signals.Figure 26
Filtered results based on spectral kurtosis of AN3 signals: (a) filtered signal; (b) frequency spectrum.
(a)
(b)


## 4.1. Case Study I

The first case study analyzes a vibration signal measured by the accelerometer AN6 (model: IMI622B01), which is close to bearing D (model: SKF NU-2220-ECM, near the high speed shaft). In this study, the rotation speed of the high speed shaft (denoted as HS_shaft) is 1200 rpm and the electric power is 25%. The waveform and frequency spectrum of the vibration signal are shown in Figure9. It can be seen that the rotation frequency of HS_shaft (fHS = 20 Hz) is visible in Figure 9(b), while no defect-related frequency for bearing D is identified.Figure 9
Waveform and frequency spectrum of AN6 signals: (a) waveform; (b) frequency spectrum.
(a)
(b)The WSSA algorithm is then used to process the same signal. After the CWT decomposition, wavelet coefficients of scale 1–32 are obtained and mPCA is used to find four outstanding scales, which are identified as scales 15, 18, 19, and 20. Then the runs test method is utilized to compute the statisticZ1=192.89, Z2=265.10, and Z3=217.06, among which Z1 is the smallest. This means one stationary source and three nonstationary sources can be determined through the SSA. The results for extracted stationary and nonstationary source components are shown in Figure 10, among which the nonstationary source 3 with the maximum kurtosis value is selected and processed using enveloping spectral analysis. From the envelope spectrum shown in Figure 11, the inner raceway defect frequency of the bearing D (fBPFI = 50 Hz) can be clearly identified. In fact, the bearing D indeed suffered from the inner raceway defect as shown in Figure 12. To make a comparison, the wavelet coefficients at the scale with the maximum kurtosis are directly used for enveloping spectral analysis, and the defect frequency 50 Hz cannot be identified from the results shown in Figure 13. In addition, the ICA method is used to process the same vibration signal and conduct the same diagnostic procedure as the SSA does. As shown in Figure 14, although the defect frequency can also be identified, it is not so prominent as compared to the result obtained using the SSA. Empirical mode decomposition (EMD) technique is also used to decompose the signal and the intrinsic mode function (IMF) with maximum kurtosis is chosen for enveloping spectral analysis. From the result shown in Figure 15, the defect frequency is nearly submerged in the surrounding frequency components. Furthermore, the spectral kurtosis is applied to the vibration signals to find the proper central frequency 2500 Hz and bandwidth 312 Hz (shown in Figure 16), and then a band-pass filter with the chosen central frequency and bandwidth is used to filter the vibration signal. From the envelope spectrum of the filtered signal shown in Figure 17, the magnitude of the defect frequency 50 Hz is not large enough as compared to other adjacent frequency components. This case study verifies the effectiveness of the proposed WSSA for wind turbine gearbox fault diagnosis.Figure 10
SSA results of AN6 signals.Figure 11
Envelope analysis of nonstationary source with maximum kurtosis: (a) waveform; (b) frequency spectrum.
(a)
(b)Figure 12
Inner race defect of bearing D [21].Figure 13
CWT and maximum kurtosis comparison results of AN6 signals: (a) waveform; (b) frequency spectrum.
(a)
(b)Figure 14
ICA comparison results of AN6 signals: (a) waveform; (b) frequency spectrum.
(a)
(b)Figure 15
EMD comparison results of AN6 signals: (a) waveform; (b) frequency spectrum.
(a)
(b)Figure 16
Spectral kurtosis of AN6 signals.Figure 17
Filtered results based on spectral kurtosis of AN6 signals: (a) filtered signal; (b) frequency spectrum.
(a)
(b)


## 4.2. Case Study II

The second case study analyzes a vibration signal measured by the accelerometer AN3 (model: IMI 626B02), which is close to bearing H (model: FAG SL181856E). In this study, the rotation speed of the high speed shaft (HS_shaft) is 1800 rpm and the electric power is 25%. The bearing H is used to support the main shaft (rotation speed: 22.09 rpm). The waveform and frequency spectrum of the vibration signal are shown in Figure18, where no defect-related frequency for bearing H is identified.Figure 18
Waveform and frequency spectrum of AN3 signals: (a) waveform; (b) frequency spectrum.
(a)
(b)The vibration signal is then processed by the proposed WSSA approach. Following the same procedure as shown in case study I, wavelet coefficients of the vibration signal at scales 1–32 are obtained by the CWT. Then the coefficients at the three scales 23, 24, and 25 are selected by the mPCA, and the runs test is performed to calculate the statisticZ1=409.68 and Z2=350.80, among which Z2 is smaller than Z1. This means two stationary sources and one nonstationary source can be extracted from the vibration signal through the SSA, and the result is shown in Figure 19. The enveloping spectral analysis of the nonstationary source is then conducted and the resulting envelope spectrum is shown in Figure 20. In addition to the intermediate speed shaft (denoted as IS_shaft) rotation frequency fIS = 7.5 Hz, the outer raceway defect frequency of the bearing H (fBPFO = 8.6 Hz) can also be clearly identified. This is consistent with the physical examination for the bearing H as shown in Figure 21. For comparison study, the methods using the CWT with maximum kurtosis-based scale selection, ICA, EMD, and the band-pass filter based on spectral kurtosis are also used to process the same vibration signal and the results are shown in Figures 22–26, respectively. From these figures, it can be seen that the rotation frequency of the IS_shaft is highlighted, while the outer raceway defect frequency of the bearing H is almost submerged in the spectrum and cannot be identified. The comparison results further confirm the effectiveness of the WSSA for wind turbine gearbox fault diagnosis.Figure 19
SSA result of AN3 signals.Figure 20
Envelope analysis of nonstationary source with maximum kurtosis: (a) waveform; (b) frequency spectrum.
(a)
(b)Figure 21
Outer race defect of the bearing H [21].Figure 22
CWT and maximum kurtosis comparison results of AN3 signals: (a) waveform; (b) frequency spectrum.
(a)
(b)Figure 23
ICA comparison results of AN3 signals: (a) waveform; (b) frequency spectrum.
(a)
(b)Figure 24
EMD comparison results of AN3 signals: (a) waveform; (b) frequency spectrum.
(a)
(b)Figure 25
Spectral kurtosis of AN3 signals.Figure 26
Filtered results based on spectral kurtosis of AN3 signals: (a) filtered signal; (b) frequency spectrum.
(a)
(b)


## 5. Conclusions

This paper presents a single sensor based blind source separation approach (i.e., the WSSA) for wind turbine gearbox fault diagnosis by integrating the CWT with the SSA. The modified PCA (mPCA) and runs test methods are brought in for the presented approach to quantitatively choose the suitable wavelet scales and the number of stationary sources, respectively. The proposed method can use a single sensor signal for SSA separation based on the signal distribution without any independency assumption. Furthermore, the mPCA can effectively reduce the computation load of the SSA by wavelet coefficient dimension reduction and the runs test can ensure an optimal separation mode for the SSA. In the experimental study performed on a real wind turbine gearbox, the proposed WSSA approach can well separate the stationary signal components from the nonstationary signal components mixed in the vibration signal and accurately identify the defect frequency of the bearing in the wind turbine gearbox. In addition, comparison study using CWT with maximum kurtosis, ICA, EMD, and spectral-kurtosis-based approaches indicates that the proposed WSSA approach is more beneficial than these existing methods for wind turbine gearbox diagnosis. However, the influence of the time epoch numberN proposed in Section 2.2 on the SSA result is not considered in this paper which needs to be further investigated in the future study. In addition, it should be noted that the proposed diagnosis approach is tested using the signal collected from constant speed gearbox. By combining this approach with order tracking method, it is able to deal with the signals under variable-speed conditions and can also be applied to monitoring variable-speed wind turbines.



---

*Source: 6971952-2015-11-16.xml*