Weak Fault Feature Extraction and Enhancement of Wind Turbine Bearing Based on OCYCBD and SVDD

Abstract

The fault feature of wind turbine bearing is usually very weak in the early injury stage, in order to accurately identify the defect location, an original approach based on optimized cyclostationary blind deconvolution (OCYCBD) and singular value decomposition denoising (SVDD) is put forward to extract and enhance the fault feature effectively. In this diagnosis method, the fast spectral coherence is fused with the equal step size search strategy for the cyclic frequency parameter and the filter length parameter optimization, and a new frequency weighted energy entropy (FWEE) indicator which combining the advantages of the frequency weighted energy operator (FWEO) and the Shannon entropy, is developed for deconvolution signal evaluation during parameter optimization process. In addition, a novel singular value order determination approach based on fitting error minimum principle is utilized by SVDD to enhance the fault feature. During the process of defect identification, OCYCBD with the optimal parameters is firstly used to recover the informative source from the collected vibration signal. FWEO is further utilized to highlight the potential impulsive characteristics, and the instantaneous energy signal of deconvolution result can be acquired. The whole interferences contained in the instantaneous energy signal can’t be removed due to the weak fault signature and the severe background noise. Then, SVDD is applied to purify the instantaneous energy signal of deconvolution signal, by which the residual interference component is eliminated and the fault feature is strengthened immensely. Finally, frequency domain analysis is performed on the denoised instantaneous energy signal, and the defect location identification of wind turbine bearing can be achieved through analyzing the obvious spectral lines in the obtained enhanced energy spectrum. The collected signals from the experimental platform and the engineering field are both utilized to verify the feasibility of proposed method, and its superiority is further demonstrated through comparing with several well known diagnosis methods. The results indicate this novel method has distinct advantage on bearing weak feature extraction and enhancement.

1. Introduction

As the joint of wind turbine, rolling bearing is indispensable and important component during wind turbine operation. It is well known that local defect occurring on the bearing is one of the primary reasons for wind turbine failure. In the wind power field, without proper detection and maintenance, bearing flaw may lead to non-planned shutdown or even result in catastrophic accident. Therefore, the incipient fault detection of rolling bearing is of great significant to ensure stable operation of wind turbine [1,2,3].

In the early injury stage of wind turbine bearing, the energy of strike vibration generated by local defect is usually very weak. In additional, the complicated transmission path between the defect location and the installed sensor would influence the response peculiarity of collected signal. And the amplitude of impulsive feature would be further weakened. Thus, the impulsive feature is often masked by the external noises and the unknown vibration sources [4]. Owing to these reasons, it is a great challenge to identify the incipient defect of wind turbine bearing. Aim at this problem, a lot of diagnosis methods have been utilized for processing the bearing incipient defect signal, such as wigner-ville distribution (WVD), wavelet transform (WT), empirical mode decomposition (EMD), and spectral kurtosis (SK). And applications of these methods have already got a certain amount of effectiveness [5]. Nevertheless, some inherent weaknesses of these methods also impair their performances on fault feature extraction. For instance, the cross-term interference in WVD is difficult to avoid when processing mixed signal [6]. WT doesn’t possess adaptive signal decomposition ability because mother wavelet function and decomposition level need to be predefined in advance [7]. EMD has the unresolved problems of noise sensibility, endpoint divergence and aliasing effect, which may cause the obtained intrinsic mode function components losing specific meanings [8]. As for SK, the obtained parameters for filter design may be unreasonable due to random impact interference, and the pass-band of constructed filter may not cover the whole resonant region by reason of fixed tiling pattern [9].

Considering the spreading effect of complicated transmission path of unknown time-variant system, the blind deconvolution operation is investigated to recover the impulsive features from the acquired signal. As the premier deconvolution technology, minimum entropy deconvolution (MED) has been successfully utilized for bearing defect identification [10], but its performance is weakened due to it preferably recovers a large random impact rather than the periodic impacts [11]. Subsequently, maximum correlated kurtosis deconvolution (MCKD) [12] technology is further developed to overcome the drawback of MED, and satisfactory diagnosis results have been achieved in some cases [13,14]. However, some inherent disadvantages also limit its engineering application. For instance, the correlated kurtosis indicator is proposed empirically without researching its statistical properties. In additional, the resample operation needs to be carried out when the deconvolution period isn’t integer, which is easy to cause misdiagnosis. Recently, a novel cyclostationary blind deconvolution (CYCBD) approach, which is able to eliminate the influence of vibration transmission path and separate the fault source effectively, is investigated by defining a simple and effective indicator of higher-order cyclostationarity maximization [15]. Nevertheless, the performance of CYCBD is mostly depended on the cyclic frequency and the filter length parameters, which also need to be artificially selected in advance. Thus, the important problem of applying CYCBD is how to set these two parameters appropriately for the best performance. In order to acquire the optimal deconvolution result automatically, the fast spectral coherence is fused with the equal step size search strategy for parameter optimization, and an optimized cyclostationary blind deconvolution (OCYCBD) method is put forward in this paper, whose superior performance on weak fault feature extraction is given expectation.

Teager Kaiser energy operator is a non-stationary signal processing technique, which has been widely utilized for tracking and enhancing the transient impact due to favorable time resolution and high computational efficiency [16]. However, it is susceptible to the background noises and become invalid under the condition of severe interferences [17]. Recently, the frequency weighted energy operator (FWEO) based on envelope-derivative operation is investigated as an alternative technology to track the instantaneous energy of a given signal [18]. It has been verified that this technology is able to amplify the impact feature without any preprocessing and exhibits better performance and higher robustness in the presence of background noises and multiple interferences [19]. As an effective signal component separation approach, singular value decomposition (SVD) is suitable for bearing fault signal denoising due to its superiorities of easy implementation and few control variables. The Hankel matrix containing the fault signature and the noise interference can be decomposed into different subspaces corresponding to the singular value sequence. And the signal denoising operation can be implemented by simply reconstructing the matrix with a number of larger singular values. This admirable trait of SVD is a benefit to enhancing fault feature during bearing injury judgment [20]. Inspired by the respective advantages of the frequency weighted energy operator and the singular value decomposition, the instantaneous energy signal obtained by FWEO operation is further purified by singular value decomposition denoising (SVDD), whose performance on noise interference suppression and fault feature enhancement is promising.

Based on the above discussion, for the sake of dealing with the weak defect identification problem of wind turbine bearing, an effective feature extraction and enhancement method combining OCYCBD with SVDD is put forward in this paper. Both the experimental signal and the actual engineering case are applied to verify this method. The chapters of this paper are structured as follows. The theoretical background of CYCBD, the influences of key parameters, and the parameter optimization strategy of OCYCBD are introduced in Section 2. Section 3 is devoted to the description of SVDD approach. Section 4 provides the detailed procedures of proposed diagnosis method. Section 5 presents the analysis results and the comparison results of roller defect experimental signal. In Section 6, an actual engineering case is utilized to further verify the proposed method. And the final conclusions are drawn in Section 7.

2. Optimized Cyclostationary Blind Deconvolution

2.1. Theoretical Background of Cyclostationary Blind Deconvolution

The purpose of blind deconvolution is to separate the input signal source $s_0$ from the collected signal $x$ by constructing an inverse filter:

$$s=x\ast h=(s_0\ast g)\ast h\approx s_0$$

(1)

where $g$ represents the frequency response function of an unknown system, $h$ denotes the inverse filter, $s$ refers to the estimated signal source and $\ast$ indicates the convolution operation. The convolution operation for discrete signal can be described as the following matrix form:

$$s=Xh$$

(2)

$$\left[\begin{array}{c}s[L-1]\\ ⋮\\ s[N-1]\end{array}\right]=\left[\begin{array}{ccc}x[L-1]& \cdots & x[0]\\ ⋮& \ddots & ⋮\\ x[N-1]& \cdots & x[N-L-2]\end{array}\right]\left[\begin{array}{c}h[0]]\\ ⋮\\ h[L-1]\end{array}\right]$$

(3)

where N and L are the length of $s$ and $h$, respectively.

A new deconvolution algorithm based on cyclostationarity maximization is proposed in reference [15]. In this algorithm, the cyclic frequency is defined as:

$$a=1/T_s$$

(4)

where $T_s$ is the cycle related to fault occurrence rate.

The novel indicator called second-order cyclostationarity to drive the deconvolution process is described as follows:

$$IC{S}_2=\frac{{\displaystyle {\sum }_{k>0}|c_s^k{|}^2}}{|c_s^0{|}^2}$$

(5)

with

$$c_s^k=\langle |s{|}^2,{\mathrm{e}}^{\mathrm{j}2\mathsf{\pi }kan}\rangle =\frac{1}{N-L+1}{\displaystyle \sum _{n=L-1}^{N-1}{|s[n]|}^2{\mathrm{e}}^{-\mathrm{j}2\mathsf{\pi }kan}}$$

(6)

$$c_s^0=\frac{{\Vert s\Vert }^2}{N-L+1}$$

(7)

Equations (6) and (7) can be further expressed as following matrix form:

$$c_s^k=\frac{{\mathrm{E}}^{\mathrm{H}}|s{|}^2}{N-L+1}$$

(8)

$$c_s^0=\frac{s^{\mathrm{H}}s}{N-L+1}$$

(9)

where $|s{|}^2={\left[|s[L-1]{|}^2,\cdots ,|s[N-1]{|}^2\right]}^{\mathrm{T}}$, $\mathrm{E}=\left[{\mathrm{e}}_1\cdots {\mathrm{e}}_k\cdots {\mathrm{e}}_K\right]$, ${\mathrm{e}}_k={\left[{\mathrm{e}}^{-\mathrm{j}2\mathsf{\pi }ka(L-1)}\cdots {\mathrm{e}}^{-\mathrm{j}2\mathsf{\pi }ka(N-1)}\right]}^{\mathrm{T}}$, and $k$ denotes the sample index.

Based on above Equations, Equation (5) can be expressed as:

$$IC{S}_2=\frac{|s{|}^{2\mathrm{H}}{\mathrm{EE}}^{\mathrm{H}}|s{|}^2}{|s^{\mathrm{H}}s{|}^2}$$

(10)

The signal which contains the periodic component $|s{|}^2$ corresponding to all the cyclic frequencies of interest $ka$, can be described as:

$$\Re [s]=\frac{1}{N-L+1}{\displaystyle \sum _k{\mathrm{e}}_k({\mathrm{e}}_k^{\mathrm{H}}|s{|}^2)=\frac{{\mathrm{EE}}^{\mathrm{H}}|s{|}^2}{N-L+1}}$$

(11)

Substituting Equations (2) and (11) into Equation (10), and the obtained final outcome is:

$$IC{S}_2=\frac{h^{\mathrm{H}}{X}^{\mathrm{H}}WXh}{h^{\mathrm{H}}{X}^{\mathrm{H}}Xh}=\frac{h^{\mathrm{H}}{R}_{XWX}h}{h^{\mathrm{H}}{R}_{XX}h}$$

(12)

where the expression of weighting matrix $W$ is:

$$W=\mathrm{diag}\left(\frac{\Re [|s{|}^2]}{s^{\mathrm{H}}s}\right)(N-L+1)=\left[\begin{array}{ccc}\ddots & & 0\\ & \Re [|s{|}^2]& \\ 0& & \ddots \end{array}\right]\frac{(N-L+1)}{{\displaystyle {\sum }_{n=L-1}^{N-1}|s{|}^2}}$$

(13)

Equation (12) is the core of CYCBD. By solving Equation (12) through Equation (13), the deconvolution signal exhibiting the maximum second-order cyclostationarity can be extracted according to the cyclic frequency.

2.2. Research on the Influences of Key Parameters

As a kind of parametric signal processing method, CYCBD needs to set the cyclic frequency parameter $a$ and the filter length parameter L in advance. To research the influences of these two parameters on the deconvolution performance, it is assumed that a local defect occurs on the inner ring of rolling bearing, and a simulated fault signal comprising periodic impact components and Gaussian white noises is constructed as follows [21,22]:

$$x(t)={\displaystyle \sum _{i=0}^{M-1}\left[[(1+A\cos (2\pi f_{r}t))\delta (t-iT)]\otimes [e^{-ct}\cos (2\pi f_{n}t)]\right]}+n(t)$$

(14)

where $\otimes$ denotes the convolution operation, $\delta (t)$ represents the Dirac delta function and c = 700 rad/s is the structural attenuation factor of bearing system. $M$ = 60, $f_r$ = 30 Hz and $f_n$= 4000 Hz respectively represent the number of impact, the rotating frequency and the resonant frequency excited by defect point strike. $A$ = 0.3 refers to the impact amplitude fluctuation caused by rotating frequency modulation, and fault feature frequency of inner ring is set to 120 Hz, which is equal to the reciprocal of impact period $T$. Furthermore, the sample rate and sample length are set to 16,000 Hz and 8192 points, respectively. $n(t)$ denotes the added noises with standard deviation of 1, and the SNR of simulated signal is -10dB, which can be calculated by the formula in reference [23].

Figure 1 displays the waveforms as well as the corresponding spectra of simulated fault signal. After adding the Gaussian white noises, the periodic impact components can’t be found in the mixed signal. Meanwhile, the feature frequency spectral peak is indistinguishable not only in the frequency spectrum but also in the envelope spectrum.

Then this simulated fault signal is applied to research the performance of CYCBD with different parameters on noise interference suppression and fault feature extraction. Firstly, the influence of filter length parameter L is discussed when cyclic frequency parameter is fixed at $a$ = 120 according to the fault feature frequency 120Hz. The waveforms and envelope spectra of deconvolution signals gained by CYCBD with filter lengths of 50,100,150 and 200 are shown in Figure 2. In order to quantificationally compare the analysis results, the feature energy ratio (FER) index of envelope spectrum used in reference [24] is taken as the criterion in this paper to evaluate the obtained different deconvolution signals, and the FER index is described as:

$$\mathrm{FER}=\left({\displaystyle \sum _{k=1}^K{F}_k}/{\displaystyle \sum _{i=1}^N{E}_i}\right)\times 100\%$$

(15)

where $F_k$ and K respectively denote the amplitude of kth feature frequency spectral line and the number of spectral line, and we set K = 4 here. $E_i$ and N respectively refer to the amplitude of ith frequency component and the number of frequency component.

In Figure 2, the fault feature frequency and its harmonics are all successfully gained by CYCBD with different filter lengths. But the FER indexes of these four envelope spectra are different, which means there exist gaps among the obtained results, and this conclusion also can be drawn by comparing the waveforms and the feature frequency components. Among these results, it is indicated that the extracted features of the third deconvolution signal, whose FER index is the largest value 22.53%, are the most obvious. And the presented features of the second deconvolution signal, whose FER index is the smallest value 7.32%, are the worst. Then the comparison results show the performance of CYCBD isn’t proportional to the filter length. Whether the selected filter length parameter is too big or too small, it will not be conducive to fault feature extraction.

After that, the influence of cyclic frequency parameter $a$ is further discussed when filter length is set to L = 150. The identical simulated fault signal is processed by CYCBD under the conditions of selected cyclic frequencies deviating from the authentic feature frequency 120Hz, and the obtained results when cyclic frequency $a$ is set to 120, 120.1, 120.5 and 121 are respectively displayed in Figure 3. Based on these results, we can find the first deconvolution signal obtained by CYCBD with $a$ = 120 are closely similar to the second deconvolution signal when $a$ = 120.1. The FER index of the second deconvolution signal stills remain large value 21.89% and the feature frequency components are also effectively discovered from the simulated signal when the cyclic frequency deviation reaches to 0.1Hz (0.0833% of fault feature frequency), which means CYCBD has definite tolerance to inaccurate cyclic frequency. However, Figure 3f,h illustrate the FER index decreases rapidly when the deviation comes up to 0.5Hz (0.4167% of fault feature frequency), and CYCBD fails to discover feature information when the predetermined cyclic frequency $a$ = 121 deviates from the real value up to 0.833%. Then the comparison results show the cyclic frequency must be set as precise as possible to get the favorable result.

Based on the above discussions, it is indicated that the performance of CYCBD on bearing fault signal processing is seriously depended on the setting filter length parameter and the cyclic frequency parameter, and the influence of cyclic frequency is more serious compared with filter length. To effectively acquire the optimal feature extraction result, a novel OCYCBD method is put forward in the subsequent section, in which the optimal influencing parameters is confirmed by fusing the fast spectral coherence with the equal step size search strategy, and the whole optimization process is guided by the proposed frequency weighted energy entropy indicator.

2.3. Frequency Weighted Energy Entropy Indicator

For a given signal $x(t)=A\cos (\omega t+\varphi )$, its analytic signal $P(t)$ can be expressed as [25]:

$$P(t)=x(t)+\mathrm{j}H[x(t)]$$

(16)

where $H[·]$ represents the Hilbert transform.

The instantaneous energy of $x(t)$ is usually represented by the square of the module of the analytic signal $P(t)$:

$$E[x(t)]={|x(t)+\mathrm{j}H[x(t)]|}^2=A^2$$

(17)

Through the above equation, it can be found that only the amplitude information is involved during the calculation process of instantaneous energy of $x(t)$, but the frequency information is ignored. In order to avoid this drawback, the derivative operation is regarded as the weighting filter, and the frequency weighted energy operator (FWEO) is defined as follows:

$$\Gamma [x(t)]={|\dot{x}(t)+\mathrm{j}H[\dot{x}(t)]|}^2={\dot{x}}^2(t)+H{[\dot{x}(t)]}^2=A^2{\omega }^2$$

(18)

where $\dot{x}(t)$ is the first-order derivative of $x(t)$.

By performing FWEO operation on the given signal, the instantaneous energy signal which contains the amplitude information as well as the frequency information can be gained. For the discrete signal $x(n)$ with length of N, its FWEO operation is expressed as:

$$\Gamma [x(n)]=\frac{1}{4}\left[x^2(n+1)+x^2(n-1)+h^2(n+1)+h^2(n-1)\right]+\frac{1}{2}\left[x(n+1)x(n-1)+h(n+1)h(n-1)\right]$$

(19)

where $h(n)=H[x(n)]$.

In the OCYCBD method, an evaluation indicator is needed to fairly assess the quality of deconvolution signal during parameter optimization process. The traditional kurtosis indicator has been frequently used to evaluate the impact feature of time domain waveform. However, this indicator will change sharply if an accidental impact appears in the waveform, which may lead to erroneous assessment due to its sensibility to accidental impact. As a nonlinear instantaneous energy tracking approach, frequency weighted energy operator (FWEO) is in a position to strengthen the impact symptom effectively [19]. And Shannon entropy has a good capacity of reflecting the uniformity of a given time series [26]. Then a new frequency indicator called weighted energy entropy (FWEE), which combining the respective specialties of FWEO and Shannon entropy, is proposed to overcome the drawback of traditional kurtosis indicator. And it is utilized to evaluate the deconvolution signal and guide the parameter optimization process. The expression of FWEE indicator of discrete signal $x(n)$ is as follows:

$$\{\begin{array}{l}\mathrm{FWEE}=-{\displaystyle \sum _{n=1}^N{p}_{n}In{p}_n}\\ p_n=\Gamma [x(n)]/{\displaystyle \sum _{n=1}^N\Gamma [x(n)]}\end{array}$$

(20)

Four simulated signals are constructed to demonstrate the feasibility of this indicator in describing the richness of periodic impact in the time series. The sample frequency and signal length are respectively set to 8192Hz and 0.1 s. The first simulated signal $x_1(t)$ is constructed as Equation (14), which presents the periodic impact signal of bearing defect without noises. The second simulated signal $x_2(t)$ denotes the Gauss white noises with standard deviation of 0.125. The third simulated signal $x_3(t)$ is mixture of $x_1(t)$ and $x_2(t)$. The fourth simulated signal $x_4(t)$, which presents the superimposed harmonic signal, is expressed as follows:

$$x_4(t)=\cos (60\pi t)+\sin (120\pi t)$$

(21)

Figure 4 illustrates the waveforms of four different simulated signals and their FWEE indicators. It can be seen the FWEE indicator of $x_1(t)$ is the smallest value 5.897, $x_2(t)$, $x_3(t)$ and $x_4(t)$ have larger indicator values, which means FEWW indicator is able to effectively distinguish these four signals due to it can reflect the subtle changes of time series. If the simulated signal contains more prominent periodic impact component, then the regularity of this signal is more distinct and the FWEE indicator is lower. Otherwise, the value of this indicator is higher.

We further add different degrees of white noises into $x_1(t)$ to acquire the simulated signals with different SNRs. Figure 5 displays the histogram of FWEE indicators corresponding to different SNRs. By comparing the FWEE values under the conditions of different SNRs, the effectiveness of this indicator in evaluating the richness of periodic impact component can be further verified. It is well known that the periodic impact feature will be more obvious if SNR is higher. As shown in the histogram, the indicator value tends to decline when the SNR increases, which means lower indicator represents the periodic impact component is richer. This is owing to the impact regularity of simulated signal decreases in the presence of noise interferences. Thus, FWEE can be applied as the guidance indicator to optimize the key influencing parameters of CYCBD.

2.4. Optimal Parameter Selection Stragegy Guided by FWEE Indicator

In Section 2.2, the influences of filter length and cyclic frequency have been investigated through the simulated fault signal. And the obtained CYCBD result will be unsatisfactory if these two parameters are selected unreasonably. Thus, the important problem of applying CYCBD to bearing fault feature extraction is how to choose the key influencing parameters appropriately for the best performance. We have known that the cyclic frequency parameter is closely related to the fault feature frequency. In terms of this parameter selection, if it is directly confirmed by referring to the theoretical fault characteristic frequency, then desired result will be difficult to obtain. The reason is because the actual feature frequency is often inconsistent with the theoretical calculating value due to the inevitable random slippage of roller during bearing operation. As for filter length selection, there is not any reference. If this key parameter is artificially chosen only through the operator’s experiences, then the performance of CYCBD may be weakened. In order to essentially overcome the inherent deficiency of influencing parameters selection and automatically acquire the optimal deconvolution result, the OCYCBD method is put forward here. In this method, firstly, the cyclic frequency $a$ is precisely confirmed by the fast spectral coherence and the equal step size search strategy, then the filter length $L$ is further reasonably selected to highlight the superiority of CYCBD on the basis of satisfying the prime requirement of $a$. The detailed implementation steps are as follows:

Step 1

Confirm the search center of cyclic frequency. It is a great challenge to determinate the cyclic frequency parameter with a high precision owing to the wide search scope. Thus, for the purpose of improving the search efficiency and easing the computational burden, the search center of cyclic frequency needs to be confirmed firstly, and the whole search process is carried out around this center. Spectral coherence provides a new interpretation of periodic flows of energy across the analysis frequency band and the cyclic frequency band [27]. With the help of its excellent ability in revealing the presence of modulation and describing the cyclostationarity, the fast spectral coherence is applied to accurately estimate the search center, and the relevance theory of fast spectral coherence can refer to literature [28]. The cyclic frequency location corresponding to the maximum energy distribution in the fast spectral coherence is regarded as the search center $a^M$. In order to avoid the interference as far as possible, only the regions [$f-df$, $f+df$] ($f\in \{fi,fo,fe,fc\}$) around the theoretical defect feature frequencies are analyzed, where $fi$, $fo$, $fe$ and $fc$ respectively represent the theoretical defect feature frequencies of inner ring, outer ring, roller, and cage. In this paper, the analysis region interval $df$ is set to 5Hz. However, the precise cyclic frequency parameter is still difficult to determinate due to the influence of cyclic frequency resolution. Then more precise search process around the search center $a^M$ is performed in the following steps.

Step 2

Determinate the search scope and search step size. The search scope of cyclic frequency parameter is initialized as [$a^I$, $a^E$] = [floor($a^M-\Delta f$), ceil($a^M+\Delta f$)] with search step size of 0.05Hz. Where $a^E$ and $a^I$ are the upper boundary and the lower boundary of search scope, $\Delta f$ denotes the cyclic frequency resolution of spectral coherence, floor $(·)$ and ceil $(·)$ respectively refer to the round down operation and the round up operation. As for the filter length $L$, the deconvolution signal will be distorted if it is too large, while the treatment efficiency will be inconspicuous if this parameter is too small. Referring to literature [15], the search scope of filter length is set as [$L^I$, $L^E$] = [20,200] with search step size of 1, where $L^E$ and $L^I$ denote the upper boundary and the lower boundary of search scope.

Step 3

Determinate the optimal cyclic frequency parameter. In the circumstances of different $L$ values, the mean of frequency weighted energy entropy (MFWEE) is calculated to select the optimal cyclic frequency $a^o$. When $a$ = $a^I$, the different filter lengths $L$ = [20,40,60,$\cdots$,200] are respectively substituted into CYCBD, and the corresponding 10 set of deconvolution signals are obtained. Then the FWEE indicator of each deconvolution signal is calculated and the mean of 10 FWEE indicators corresponding to $a^I$ can be acquired. For $a=a^I+0.05k$(k = 0,1,2,$\cdots$, $(a^E-a^I)/0.05$ is the number of step size movement), the MFWEE value corresponding to each cyclic frequency $(a^I+0.05k)$ is in turn calculated based on the similar process. And all MFWEE values corresponding to the whole search scope [floor ($a^M-\Delta f$), ceil ($a^M+\Delta f$)] can be obtained. If the MFWEE value is smaller, then it means the periodic impact signature in the obtained deconvolution signal is more prominent and the contained useful information is richer. Thus, during the process of key parameter optimization, the cyclic frequency corresponding to the MFWEE minimum value is regarded as the optimal parameter $a^o$, then store this optimal parameter in the memory.

Step 4

Select the optimal filter length parameter. In order to further optimize the filter length on the basis of $a$ = $a^o$, $L$ = [20,21,22,$\cdots$,200] are substituted into CYCBD respectively, and the corresponding 181 set of deconvolution signals are obtained, then the FWEE indicator of each deconvolution signal is calculated to evaluate the filter length parameter. For $L=L^I+g$ (g = 0,1,2,$\cdots$, $(L^E-L^I)$ is the number of step size movement), the FWEE indicator corresponding to each filter length $(L^I+g)$ is in turn calculated based on the similar process. And all FWEE indicators corresponding to the whole search scope [20,200] can be obtained. Then confirm the optimal filter length parameter $L^o$ when FWEE indicator is the smallest.

3. Singular Value Decomposition Denoising

3.1. Basic Theory of Singular Value Decomposition

The special property of singular value decomposition (SVD) has attracted a lot of attentions in the field of signal denoising. And this approach is capable of clearing off the redundant interferences by simply reconstructing the matrix with the chosen effective singular values [29].

Assume that $X=[x(1),x(2),\cdots ,x(L)]$ is a given discrete signal, and the constructed Hankel matrix for conducting SVD is as follows:

$$A={\left[\begin{array}{cccc}x(1)& x(2)& \cdots & x(l)\\ x(2)& x(3)& \cdots & x(l+1)\\ ⋮& ⋮& \cdots & ⋮\\ x(L-l+1)& x(L-l+2)& \cdots & x(L)\end{array}\right]}_{m\times l}$$

(22)

where $l=L/2$ and $m=L/2+1$.

The orthogonal decomposition of matrix A is the essence of SVD, it can be described as:

$$A=UE{V}^{\mathrm{T}}$$

(23)

where $E=(diag({\sigma }_1,{\sigma }_2,\cdots ,{\sigma }_p),0)\in R^{m\times l}$ is the diagonal matrix, $U=(u_1,u_2,\cdots ,u_m)\in R^{m\times m}$ and $V=(v_1,v_2,\cdots ,v_l)\in R^{l\times l}$ refer to the orthogonal matrixes, ${\sigma }_1\ge {\sigma }_2\ge \cdots \ge {\sigma }_p>0$ indicates the gained singular value sequence.

Another expression form of matrix $A$ can be described as:

$$A={\sigma }_1{u}_1{v}_1^{\mathrm{T}}+{\sigma }_2{u}_2{v}_2^{\mathrm{T}}+\cdots +{\sigma }_p{u}_p{v}_p^{\mathrm{T}}=A_1+A_2+\cdots +A_p$$

(24)

For the purpose of removing the interference components, the matrix $A^{\mathrm{R}}$ corresponding to the denoised signal needs to reconstruct through the selected former i singular values as follows:

$$A^{\mathrm{R}}={\sigma }_1{u}_1{v}_1^{\mathrm{T}}+{\sigma }_2{u}_2{v}_2^{\mathrm{T}}+\cdots +{\sigma }_i{u}_i{v}_i^{\mathrm{T}}=A_1+A_2+\cdots +A_i$$

(25)

And the purified signal is able to be gained through diagonal mean operation of $A^{\mathrm{R}}$.

3.2. Singular Value Order Determination

The order of singular value needs to be predefined in denoising process. And the result is closely depended on the selected singular value order. However, this problem is a difficult barrier for lack of priori knowledge about the original signal, which is very complex in actual engineering application. The useful feature will be lost if the selected singular value order is too small using the traditional difference spectrum method, while the excessive redundant noises will be remained if the selected order is too large using the median value or the mean value method [30,31]. Usually, after performing SVD on the original signal, there exists an elbow in the obtained singular value curve. In the left side of the elbow, the singular values decrease rapidly, while the singular values in the right side of the elbow decrease slowly. Thus, the elbow of singular value curve should be given more attention. The trend lines of singular value distribution can be drawn using fitting method. And the diacritical location of the available and the useless singular values in the elbow can be determined by the degree of closeness with these two trend lines. Then, in order to select the reasonable order of reconstructed matrix during denoising process, a principle based on fitting error minimum is utilized in this paper. The realization processes of singular value order determination are as follows:

(1)

The singular value sequence $Q=[{\sigma }_1,{\sigma }_2,\cdots ,{\sigma }_p]$ and the order sequence $D=[1,2,\cdots ,p]$ are obtained by applying SVD to the given signal using Equations (22) and (23). The first singular value, which corresponding the trend component of the signal, is obvious larger than the others in the singular value sequence. In order to avoid causing large deviation, then the first singular value is removed in the whole process.

(2)

The initial order r is confirmed.

(3)

According to the order r, the original singular value sequence $Q$ is split into two set of sequences $Q_1=[{\sigma }_2,{\sigma }_3,\cdots ,{\sigma }_r]$ and $Q_2=[{\sigma }_{r+1},{\sigma }_{r+2},\cdots ,{\sigma }_p]$, and the original order sequence $D$ is also split into two set of sequences $D_1=[2,3,\cdots ,r]$ and $D_2=[r+1,r+2,\cdots ,p]$.

(4)

The order sequences $D_1=[2,3,\cdots ,r]$ and $Q_2=[{\sigma }_{r+1},{\sigma }_{r+2},\cdots ,{\sigma }_p]$ are respectively considered as the independent variables and the singular value sequences $Q_1=[{\sigma }_2,{\sigma }_3,\cdots ,{\sigma }_r]$ and $Q_2=[{\sigma }_{r+1},{\sigma }_{r+2},\cdots ,{\sigma }_p]$ are respectively regarded as the dependent variables. Then the least squares quadratic polynomial fitting method is utilized to construct the fitting function $\varsigma {}_1(k)=a_1{k}^2+a_{2}k+a_3$ and $\varsigma {}_2(g)=b_1{g}^2+b_{2}g+b_3$. Where $a_1$, $a_2$, $a_3$, $b_1$, $b_2$ and $b_3$ are the polynomial coefficients, $k=2,3,\cdots ,r$ and $g=r+1,r+2,\cdots ,p$ are the independent variables.

(5)

The fitting error ${\mathrm{FE}}_r$ corresponding to the order r is calculated by the following expression: ${\mathrm{FE}}_r=\sqrt{{\displaystyle \sum _{k=2}^r{({\sigma }_k-{\varsigma }_1(k))}^2}}+\sqrt{{\displaystyle \sum _{g=r+1}^p{({\sigma }_g-{\varsigma }_2(g))}^2}}$.

(6)

The order r is reassigned as r = r + 1.

(7)

The steps (3)–(6) are repetitive executed until r = p.

(8)

The minimum value of fitting error ${\mathrm{FE}}_i$ is confirmed, then the corresponding order i is regarded as the diacritical location of the available and the useless singular values. And the former i singular values except for the first one are applied for purified signal reconstruction.

4. Fault Feature Extraction and Enhancement Method Based on OCYCBD and SVDD

In practical engineering, the bearing injury signature is usually inconspicuous and covered by strong external noises in the early stage. For a given signal of rolling bearing with microlesion, if the impulsive fault source can be recovered from the collected observed signal, then the SNR can be improved drastically and the weak fault feature can be extracted preferably. Owing to the excellent ability of eliminating the spreading effect of unknown transmission path, the above proposed OCYCBD technology is regarded as a pre-processing method to achieve this purpose. And this will facilitate the follow-up feature identification progress. FWEO is capable of tracking and intensifying the impulsive feature. Moreover, compared with the traditional energy operator algorithm, it is more robust against severe interference noises. Thus, FWEO is further conducted to calculate the instantaneous energy signal of deconvolution signal, which can presents the regular impulsive feature better. Nevertheless, some interference components still remain in the instantaneous energy signal because the SNR of collected fault signal is very low, and the defect identification is easily affected by these redundant interference components. Thus, for the purpose of enhancing the fault signature and getting rid of the remained interferences, SVDD is applied as a post-processing approach to purify the instantaneous energy signal of deconvolution signal. On the basis of these statements, a novel feature extraction and enhancement method is put forward to improve the accuracy of defect identification for wind turbine bearing. Figure 6 illustrates the schematic of this method, and the diagnosis procedures are as follows:

(1)

Data acquisition. The original signal is collected using the corresponding data acquire equipments.

(2)

OCYCBD processing. The fast spectral coherence and the equal step size search strategy are combined organically to search for the optimal cycle frequency and the filter length, and the corresponding steps have been elaborate explained in Section 2.4. The obtained optimal parameters are substituted into CYCBD. Then the informative fault source with higher SNR can be recovered from the original collected signal by deconvolution operation.

(3)

Instantaneous energy signal calculation. The FWEO operation is carried out on the obtained deconvolution signal to calculate the instantaneous energy signal. The impulsive signature of the instantaneous energy signal is more outstanding than the original signal, but there still exist the interference components. Then the following procedure is carried out to get rid of the redundant components in the instantaneous energy signal and enhance the defect signature.

(4)

SVDD processing. The SVDD approach is further utilized for instantaneous energy signal denoising. And the effective singular value order to reconstruct the matrix can be determined using the fitting error minimum principle. Because the quadratic polynomial fitting algorithm is utilized to establish the fitting functions, thus the initial order is set to r = 4. And the denoised instantaneous energy signal can be gained through the diagonal mean operation of reconstructed matrix.

(5)

Enhanced energy spectrum analysis. The Fourier transform based spectrum analysis is performed on the purified instantaneous energy signal of deconvolution signal, and the corresponding enhanced energy spectrum of deconvolution signal can be acquired. Then the defect location of wind turbine bearing is able to be judged by analyzing the spectral peak in the enhanced energy spectrum.

5. Experiment Verification

5.1. Introduction of Experimental Platform

The vibration signal collected from experimental platform is utilized for feasibility verification of proposed method. Figure 7 displays the experimental test site. Two SKF6025 bearings were used to support the rotating shaft, which was driven by chain wheel. The rotating speed of main shaft was 1470 rpm, viz. the rotating frequency $fr$ = 24.5 Hz. It can be found the bearing on the right side is normal while the left side bearing is defective in Figure 7a. The injury was inserted on the roller of experimental bearing, as displayed in Figure 7b. The width and the depth of local defect were respective 0.2mm and 1.53mm. The bearing structure parameters are listed in

Table 1. Structure parameters of SKF6205 bearing.

Number of Balls	Diameter of Balls (mm)	Pitch Diameter (mm)	Contact Angle (°)
9	7.94	39.04	0

. The installed PCB acceleration sensors on the bearing seat are shown in Figure 7c and the signal acquisition process is displayed in Figure 7d. During the process of bearing running, the sample frequency was set as $fs$ = 12,800 Hz and the collected signal with length of 8192 points is intercepted for analysis. The defect feature frequencies of inner ring $fi$, outer ring $fo$, roller $fe$, and cage $fc$ of experimental bearing can be respectively calculated as $fi$ = 132.67 Hz, $fo$ = 87.83 Hz, $fe$ = 115.48 Hz and $fc$ = 9.76 Hz based on the structure parameters and the theoretical equations in literature [32].

5.2. Experimental Signal Analysis and Result Comparsion

Figure 8 displays the waveform of experimental signal as well as its spectra. There exists obvious impact phenomenon in the waveform, while the interval between adjacent impacts isn’t the reciprocal of any fault feature frequency. The frequency components are mainly concentrated below 2500Hz and there exists broad resonant band in the frequency spectrum. However, the spectral line related to bearing defect can’t be found in the low frequency region. The traditional envelope spectrum analysis is also executed, but there exist a large number of unknown spectral lines except for feature frequency component in the envelope spectrum.

Then the presented method is applied for analyzing this experimental signal. Firstly, the precise influencing parameters of CYCBD need to be selected. During the process of parameter optimization, the fast spectral coherence is calculated as shown in Figure 9a to determinate the search center of cyclic frequency, from which we can find the energy distribution is largest when cyclic frequency is 114.7Hz (the red block indicated by green arrow), then the search center $a^M$ is confirmed to perform more elaborate search process. As the cyclic frequency resolution of fast spectral coherence $\Delta f$ = 1.572 Hz, then the search scope is set as [$a^I$, $a^E$] = [floor (114.7−1.572), ceil(114.7 + 1.572)] = [113 Hz, 117 Hz] with search step size of 0.05 Hz. The MFWEE values corresponding to different cyclic frequencies within the search scope are illustrated in Figure 9b. And the optimal parameter $a^o$ = 115.7 Hz can be automatically confirmed in terms of MFWEE minimum value. Then the filter length parameter is further optimized on the basis of $a$ = $a^o$, and the search scope is set as [20,200] with search step size of 1. The FWEE indicators of deconvolution signals, which are obtained by CYCBD with different filter lengths and optimal cyclic frequency $a^o$ = 115.7, are shown in 9c. And the optimal filter length parameter $L^o$ = 117 can also be chosen according to FWEE minimum value. Then the optimal parameters $a^o$ = 115.7 and $L^o$ = 117 are substituted into CYCBD to recover the fault source from the original signal, and the acquired deconvolution signal is displayed in Figure 9d. Secondly, FWEO operation is further carried out to gain the instantaneous energy signal of deconvolution result, in which the impulsive phenomenon is amplified to a certain extent. Thirdly, SVDD is performed on the instantaneous energy signal to get rid of the influence of interference components. And the fitting error minimum principle is applied to adaptively select the effective singular value order during denoising process. In order to present the result more clearly, Figure 9e only displays the former 150 singular value orders and the corresponding fitting error values. As marked in this figure, the 24th singular value order is considered as the diacritical location of the available and the useless singular values. Therefore, the former 24 singular values except for the first one are adopted to construct the purified instantaneous energy signal. Finally, after performing spectrum analysis on the denoised instantaneous energy signal, the enhanced energy spectrum of deconvolution signal is obtained. In Figure 9f, the fault related spectral lines $fe~8fe$ with obvious peaks are able to be presented distinctly. Thus, we can determinate that there exists local flaw in the ball, and this conclusion is consistent to the actual condition. Then the feasibility of proposed method in getting rid of the severe background noises and intensifying the weak fault signature of rolling bearing is verified through this experimental signal.

To clearly verify the advantage of fusing OCYCBD with SVDD on noise interference suppression and fault feature enhancement, the energy spectrum of deconvolution signal without denoising operation is displayed in Figure 10. It is indicated that the fault feature frequency components are disturbed by background noises and interference spectral peaks, and this isn’t conducive to achieve precise bearing flaw detection. Owing to the favorable specialty of SVDD, the redundant interference components in the instantaneous energy signal can be removed absolutely, and the fault feature spectral lines in the enhanced energy spectrum can be presented more obvious.

In order to further demonstrate the necessity of equal step size search strategy for precise cyclic frequency selection, the parameter $a$ = 114.7 which is confirmed only by the fast spectral coherence, is substituted into CYCBD. And the filter length is still set to L = 117. The obtained results by CYCBD with these two parameters are displayed in Figure 11. It can be found the analysis results are non-ideal because the periodic impact components in the waveform aren’t abundant and the feature spectral lines in the energy spectrum are vague, which hinder the determination of bearing defect location. It is indicated that the diagnosis results may be entirely different due to the imprecise selection of influencing parameters. Fortunately, we can favorably achieve the aureate fault identification with the help of proposed parameter optimization strategy.

To demonstrate the superiority of proposed method, the well known SK method, MED method and MCKD method are respectively applied for comparison. Figure 12a illu strates the calculated kurtogram, and the optimal resonant band in level 3 is marked by black oval. Based on the obtained information in the kurtogram, the center frequency and the bandwidth of designed filter are 400Hz and 800Hz. Nerveless, only the harmonics of rotating frequency can be seen in Figure 12b. The fault feature frequency isn’t discovered due to the selected filter frequency band is unreasonable. The analysis results using MED and MCKD are shown respectively in Figure 13 and Figure 14. The key parameters of there two methods are confirmed based on the obtained search results in Figure 9b,c. As the optimal cyclic frequency is $a^o$ = 115.7 and the optimal filter length is $L^o$ = 117, thus the deconvolution period of MCKD is calculated as $fs/a^o$ = 110.63, and the filter lengths of MED and MCKD are both set to 117. The interval of periodic impact in Figure 13a isn’t related to bearing defect, and the fault feature frequency of roller isn’t able to be found from the corresponding envelope spectrum in Figure 13b. Thus, the analysis results of MED method losing their practical meanings. It is noticeable that the fault signatures in Figure 14a aren’t as outstanding as those in Figure 9d, the basic fault frequency spectral line and a lot of unrelated components are simultaneously appear in Figure 14b. But the harmonics of feature frequency aren’t detected and the obtained envelope spectrum is also non-ideal. Compared with above analysis results, it is verified the capacity of proposed method on weak fault diagnosis is superior to these widely used methods.

6. Engineering Case Verification

6.1. Description of Wind Turbine

An engineering case of generator bearing defect of wind turbine is applied for demonstrating the proposed method. The sketch of installed acceleration sensors of generator bearings is shown in Figure 15a. And the vibration signal was sampled at $fs$ = 16,384 Hz. Figure 15b illustrates the structure chart of wind turbine, and the type of generator bearings is SKF6324, whose structure parameters are listed in

Table 2. Structure parameters of SKF6324 bearing.

Number of Balls	Diameter of Balls (mm)	Pitch Diameter (mm)	Contact Angle (°)
8	41.275	190	0

. The speeds of the main shaft and the generator rotor were $n_1$ = 21.63rpm and $n_2$ = 1519rpm respectively. Then the rotating frequency of generator rotor was $fr$ = 25.3 Hz. The monitoring system showed the vibration of generator front bearing was abnormal. Thus, the collected engineering signal with length of 8192 points is intercepted and analyzed here. The theoretical fault feature frequencies of inner ring $fi$, outer ring $fo$, roller $fe$, and cage $fc$ of the abnormal generator bearing can be respectively calculated as $fi$ = 123.18 Hz, $fo$ = 79.21 Hz, $fe$ = 55.48 Hz and $fc$ = 15.39 Hz.

Figure 16 describes the waveform and the spectra of collected engineering signal. The periodic impact components in the waveform aren’t abundant. And the kurtosis index of this engineering signal, which is used to evaluate the impulsive phenomenon, is only 2.46. It is indicated that the SNR of engineering signal is very low. The frequency spectrum is very simple and there exist three energy concentration regions. In the envelope spectrum, only several peak spectral lines unrelated to any theoretical fault feature frequency of generator bearing can be found. And we can’t make any conclusion based on these limited information.

6.2. Engineering Signal Analysis and Result Comparsion

The proposed method is then applied to this engineering signal. Figure 17a is the calculated fast spectral coherence, and the search center of cyclic frequency is confirmed as $a^M$ = 80.47 Hz. Then the cyclic frequency is optimized within the scope of [78 Hz,83 Hz] and the search step size is 0.05 Hz. Figure 17b shows the calculated MFWEE values corresponding to different cyclic frequencies, and the optimal cyclic frequency $a^o$ = 79.4 Hz is selected based on the MFWEE minimum value. Within the scope of [20,200], the optimal filter length parameter is further searched with step size of 1. Figure 17c display s the calculated FWEE indicators corresponding to different filter lengths, and we can determinate the optimal filter length $L^o$ = 157 according to the FWEE minimum value. Then the optimal influencing parameters are substituted into CYCBD to perform deconvolution operation on the engineering signal, and the obtained result is displayed in Figure 17d. The periodic impact components in the waveform are very rich, and the kurtosis index of this deconvolution signal is 8.74, which means the SNR of engineering signal has been improved drastically. Then the instantaneous energy signal of deconvolution signal is further calculated by FWEO operation. And SVDD is executed to reduce the noise interference and enhance the fault signature. During the process of signal denoising, the fitting error values corresponding to different singular value orders are obtained. From Figure 17e, we can find that the order 9 accords with the requirement of effective singular value order determination principle, then the former 9 singular values excepted for the first one are adopted to acquire the purified instantaneous energy signal. The corresponding enhanced energy spectrum of deconvolution signal is calculated and displayed in Figure 17f. And the peak spectral lines $fo$~$3fo$ corresponding to outer ring fault are noticeably presented in this figure. Thus, it is highly doubt there is local defect on the outer ring of generator front bearing, and the boarding check results verify this diagnosis conclusion. Thus, the effectiveness of presented method in achieving accurate weak fault detection for wind turbine bearing is demonstrated through this engineering case.

The comparison methods of SK, MED, and MCKD are also applied to this collected engineering signal. The kurtogram is depicted in Figure 18a, and the band-pass filter is constructed in terms of the kurtogram information to process the engineering signal. Figure 18b displays the envelope spectrum of filtered signal, and the feature frequency component $fo$ and its harmonic $2fo$ are visible. Nerveless, this spectrum isn’t as satisfied as that shown in 17f. Figure 19 and Figure 20 respectively illustrate the obtained results using MED and MCKD. Several prominent spectral lines are visible in the corresponding envelope spectra of deconvolution signals. Unfortunately, these spectral lines are both unrelated to outer ring defect of generator bearing, and the analysis results using MED and MCKD are invalid. Based on these results, it is indicated that the comparison methods are unable to tackle this defect identification case satisfactorily. And the superiority of proposed method based on OCYCBD and SVDD on fault feature extraction and enhancement for wind turbine bearing diagnosis is verified once again.

7. Conclusions

A novel fault feature extraction and enhancement method based on OCYCBD and SVDD is presented to diagnosis the local defect of wind turbine bearing in the early injury stage. The contributions of this study are as follows:

(1)

The influences of the cyclic frequency and the filter length on the performance of CYCBD are researched by the simulated fault signal. In addition, the FWEE indicator, which can effectively reflect the richness of periodic impact component, is proposed to evaluate the quality of deconvolution signal during parameter optimization process.

(2)

The OCYCBD method fusing the fast spectral coherence with the equal step size search strategy is put forward to overcome the drawback of CYCBD, by which the optimal deconvolution result can be achieved automatically.

(3)

A novel fitting error minimum principle is used by SVDD to select the effective singular value order, the redundant interferences can be suppressed and the fault feature can be enhanced tremendously through denoising operation.

Both the experimental signal and the actual engineering case are applied to demonstrate this proposed diagnosis method. The analysis results show it can effectively recover the fault source and accurately identify the weak defect of wind turbine bearing. Compared with some widely used diagnosis methods, the comparison results verify the advantage of proposed method on weak fault feature extraction and enhancement. It is expected this method can be applied in a broader field. And our future research target is ready to modify this method for dealing with the compound fault problem of rolling bearing under variable speed condition.