Abstract
The defect characteristics of rolling bearing are difficult to excavate at the incipient injury phase; in order to effectively solve this issue, an original strategy fusing recursive singular spectrum decomposition (RSSD) with optimized cyclostationary blind deconvolution (OCYCBD) is put forward to achieve fault characteristic enhanced detection. In this diagnosis strategy, the data-driven RSSD method without predetermined component number is proposed. In addition, a new morphological difference operation entropy (MDOE) indicator, which takes advantage of morphological transformation and Shannon entropy, is developed for confirming the influencing parameters of cyclostationary blind deconvolution (CYCBD). During the process of fault detection, RSSD is firstly adopted to preprocess the original signal, and the most sensitive singular spectrum component (SSC) is selected by the envelope spectrum peak (ESP) indicator. Then, the grid search algorithm is adopted to precisely confirm the optimal parameters and OCYCBD is further performed as a postprocessing technology on the most sensitive component to suppress the residual interferences and amplify the fault signatures. Finally, the enhanced fault detection of rolling bearing is able to achieve by analyzing the envelope spectrum of deconvolution signal. The feasibility of the proposed strategy is verified by the simulated and the measured signals, respectively, and its superiority is also demonstrated through several comparison methods. The results manifest this novel strategy has praisable advantages on weak characteristic extraction and intensification.
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. Without proper detection and maintenance, bearing flaw may lead to wind turbine nonplanned shutdown or even result in catastrophic accident. Therefore, the incipient flaw detection of rolling bearing is of great significance to guarantee stable operation of wind turbine [1, 2].
In general, the vibration signal of wind turbine subjects to heavy background noises and multicomponent interferences [3]. The signal decomposition technology, such as local mean decomposition (LMD), wavelet decomposition, and empirical mode decomposition (EMD), can be applied to suppress noises and dig up the sensitive component containing more characteristic signatures. And applications of these methods in the field of fault identification have already got a certain amount of laudable effectiveness in recent decades [4–6]. Nevertheless, some inherent weaknesses of these methods impair their performances on fault characteristic extraction. The classic wavelet decomposition does not possess adaptive processing ability because mother wavelet function and decomposition level need to be predefined in advance [7]. The strict mathematical derivations are lacking in LMD and EMD, which also suffer from the noise sensibility, the endpoint divergence, and the aliasing effect problems, and these problems may cause the obtained signal components losing specific meanings [8]. As a novel data-driven decomposition approach, singular spectrum decomposition (SSD) can adaptively separate the complicated nonlinear and nonstationary signal into a series of simple components [9]. It obviously outperforms EMD and LMD with regard to the fine division property for high-frequency component. Owing to this excellent superiority, SSD attracts a lot of attention and has been successfully used for bearing fault signal processing [10, 11]. The segmentation process of SSD is adaptive while the component number needs to be artificially selected according to different engineering signals, which has the greatest influence on the analysis results and hinders its application in automatic fault diagnosis. Incorrect decomposition or memory overflow often occurs if component number is selected inappropriately. Until now, little effort has been made to solve these problems. Thus, to reasonably process the engineering signal and avoid above issues, the recursive singular spectrum decomposition (RSSD) is proposed in this paper.
In the early injury stage of rolling bearing, the energies of strike vibrations generated by local defect are usually very weak. In addition, the complicated transmission path between the defect location and the installed sensor affects the response peculiarity of collected signal. Thus, the fault signatures are usually submerged by the external noises and the unknown interference sources. Owing to these reasons, the direct application of RSSD without signal postprocessing is a great challenge to achieve satisfying diagnosis effect. In view of the spreading effect of vibration transmission path of the complicated wind turbine system, the blind deconvolution operation is considered to recover the fault source signal. As the premier deconvolution technology, minimum entropy deconvolution (MED) has been successfully utilized for bearing defect identification [12], but its performance is easy to be affected because it preferably recovers a large random impact rather than the periodic impacts [13]. Subsequently, maximum correlated kurtosis deconvolution (MCKD) [14] is further developed to overcome the drawback of MED, and satisfactory diagnosis results have been achieved in some cases [15, 16]. However, some inherent disadvantages also limit its engineering application. For instance, the correlated kurtosis indicator is proposed empirically without researching its statistical properties. In addition, the resample operation needs to be carried out if the deconvolution period is not an integer, which is easy to cause misdiagnosis [17]. Recently, a novel cyclostationary blind deconvolution (CYCBD) approach is investigated by defining the higher-order cyclostationarity maximization indicator [18]. Nevertheless, its performance depends mostly on the cyclic frequency and the filter length parameters, which need to be artificially set in advance through prior knowledge and the operator’s experiences, and the processing result can be entirely different due to the imprecise selection of these two parameters. Thus, the key issue of CYCBD application is how to set the influencing parameters appropriately for the best performance. To acquire the optimal deconvolution result automatically, the grid search algorithm [19] is fused with CYCBD and an optimized cyclostationary blind deconvolution (OCYCBD) method is put forward in this paper. The particular trait of OCYCBD is beneficial to fault characteristic extraction while its performance may be weakened in the presence of strong background interferences. For the sake of reducing the heavy noises, a preprocessing procedure is necessary before carrying out OCYCBD. Fortunately, the excellent specialty of RSSD is able to meet this requirement, and its processing results will be conducive to the influencing parameter optimization progress of CYCBD.
In the view of above statements, the proposed RSSD and OCYCBD methods can be naturally combined to exert their complementary advantages for dealing with the weak injury identification problem of rolling bearing, and then an enhanced fault detection strategy based on RSSD-OCYCBD is put forward in this paper. Within this strategy, the proposed RSSD method, which can overcome the drawback of component number artificial selection in the original SSD algorithm, is firstly utilized as a preprocessing approach to separate the most sensitive SSC automatically. Then, a two-phase optimization scheme based on RSSD and grid search algorithm is put forward to solve the problem of precise influencing parameter selection in the original CYCBD method. In addition, a new MDOE indicator combining the advantages of morphological difference operation and Shannon entropy is investigated to evaluate the richness of periodic impact components, and it is used to guide influencing parameter optimization process. Finally, the presented OCYCBD method is considered as a postprocessing approach to automatically acquire the optimal deconvolution result from the most sensitive SSC. The superior performances of this proposed RSSD-OCYCBD detection strategy on noise interference suppression and fault characteristic amplification are given expectations. And both the simulated signal and the actual collected signals are used for verification. The paper is structured as follows. Section 2 illustrates the basic theory of SSD and the proposed RSSD method. The background of CYCBD, the novel morphological difference operation entropy (MDOE) indicator, and the OCYCBD method are, respectively, provided in Section 3. Section 4 introduces the detailed implementation procedures of proposed diagnosis strategy. The simulated signal is constructed and analyzed in Section 5. The analysis results and the comparison verifications of the experimental signal and the engineering case are, respectively, presented in Section 6 and Section 7. And the conclusions are summarized in Section 8.
2. Recursive Singular Spectrum Decomposition
2.1. Basic Theory of SSD
As a novel kind of data-driven nonlinear nonstationary signal processing approach, singular spectrum decomposition is put forward based on the traditional singular spectrum analysis algorithm, which is a nonparametric spectral estimation technology for processing and predicting the time domain signal. The embedding dimension and the grouping process, which have significant influences on the output results, need to be confirmed according to the experiences in singular spectrum analysis algorithm. However, these key issues are confirmed automatically in SSD algorithm, which realizes the data-driven decomposition and reconstruction of single component from high frequency to low frequency. And the essence of this method is to separate the given nonlinear nonstationary signal into the sum of several singular spectrum components (SSCs) and the residual term adaptively by iterative operation. The detailed expositions about this complicated method can be found in Reference [9]. Due to the limit of the space, we just give a brief introduction to this algorithm in this paper; the brief flowchart of SSD is shown in Figure 1 and the corresponding procedures are as follows.

2.1.1. Trajectory Matrix Construction
Assume s (n) is a multicomponent signal with length of N data points; a M × N trajectory matrix can be constructed as under the condition of embedding dimension M (1 < M < N), where si = {s (i), …, s (N), s (1), …, s (i − 1)} stands for the ith row of S and i = 1, 2, …, M. An instance of constructed trajectory matrix with M = 3 for sequence s (n) = {1, 2, 3, 4, 5} is shown below to illustrate the embedding process.where the left side 3 × 3 matrix corresponds to the trajectory matrix of singular spectrum analysis (SSA) algorithm. For the purposes of intensifying the signal oscillation content and providing the useful property for energy decrease of residual term, the matrix is reconstructed in SSD as
2.1.2. Embedding Dimension Adaptive Confirmation
Firstly, the power spectrum of residual term needs to be calculated, and the residual term rj (n) in the jth iteration can be obtained aswhere rh (n) denotes the sum of former h components.
Then, the primary frequency fmax with largest amplitude in the power spectrum is estimated. In the first iteration, a considerable trend is deemed as the residual term in case the normalized frequency fmax/fs is smaller than the threshold (fs denotes the sample frequency) and the embedding dimension M is selected as floor (N/3). Otherwise, M is set to 1.2 × (fs/fmax) in the jth iteration (j > 1).
2.1.3. SSC Reconstruction
The jth SSC (n) is successively obtained by iterative operation. For the first iteration, only the first left eigenvector , the right eigenvector , and the residual amount of the decomposed signal are adopted to get (n) in case a considerable trend is detected, and (n) is reconstructed by diagonal averaging operation of matrix . Otherwise, for the jth iteration (j > 1), a subset Ij (Ij = {i1, i2, …, ip}) is constructed by choosing the characteristic group whose left eigenvector has outstanding peak in the frequency range [fmax − B, fmax + B] (B denotes the half bandwidth of outstanding peak in the power spectrum of residual term) and contributes most to the peak energy. Then, SSC (n) is reconstructed by diagonal averaging operation of matrix SIj = Si1 + Si2 + ⋯ + Sip.
2.1.4. Iterative Stopping Criterion Setting
A new residual term is obtained as after estimating each SSC . And r(j+1) (n) refers to the input of next iteration j + 1. Then, normalized mean square error (NMSE) between original signal and residual term is calculated as
The iteration will shut off if NMSE is less than the threshold; otherwise, this process continues. Then, the multicomponent signal s (n) can be finally separated into several components and a residual term adaptively:where K denotes the predetermined SSC number and represents the kth component.
2.2. The RSSD Method
The component number K needs to be predefined in SSD, and the obtained results depend closely on this parameter. The useful information will be submerged if K is chosen too small, while excessive decomposition or memory overflow may occur if this parameter is set too large. Nevertheless, the component number selection is a difficult barrier for lack of prior knowledge about the original signal, which is usually very complex in actual engineering application. Based on the decomposition principle of SSD, we can find SSC is a kind of band-limited signal with peak frequency as the center. In order to achieve automatic signal decomposition during bearing fault diagnosis, a novel recursive singular spectrum decomposition (RSSD) method is carried out in this section. In this method, the high-frequency SSC is separated from the original signal one by one through SSD operation with K = 2 until the nearness of the peak frequencies of obtained two components is less than the given threshold. Figure 2 shows the flowchart of RSSD, and the steps are as follows:(1)U (n) and Ch (n) initialization In this procedure, U (n) is firstly replaced with the original signal s (n), and the high-frequency component Ch (n) is initiated by a zero sequence.(2)High frequency component separation Before each SSD operation, the previously obtained high-frequency component Ch (n) is removed from U (n), and U (n) is updated as(3)SSD operation The updated U (n) is decomposed into the low-frequency component Cl (n) and the high-frequency component Ch (n) by SSD operation with K = 2.(4)Nearness calculation The stopping condition depends on the nearness of peak frequencies of the separated low-frequency component and high-frequency component, and the nearness indicator fd can be calculated as where fl and fh correspond to the peak frequencies of Cl (n) and Ch (n), respectively.(5)Decomposed SSC acquisition In case the nearness fd is larger than the given threshold , then the flow continues and returns to step (2), and the high-frequency component Ch (n) is regarded as a new component Ci (n). Otherwise, the excessive decomposition is considered to emerge if fd < and the separated two components are deemed to be a corporate component and combined as Ci (n) = Ch (n) + Cl (n). Meanwhile, the whole flow is over and the output are a series of components Ci (n), which are the sum of the last SSC decomposition results.

To ensure that enough useful information is contained in the signal component, the band bandwidth of signal component should be equal to at least three times the fault characteristic frequency as indicated in Reference [20]. Among all the characteristic frequencies of rolling bearing, the characteristic frequency of inner ring fault is the largest, and thus the threshold is set to three times this value in this paper.
3. Optimized Cyclostationary Blind Deconvolution
3.1. Theoretical Background of CYCBD
The purpose of blind deconvolution is to retrieve the original fault excitation source s0 from the observed signal x by constructing an inverse filter:where represents the frequency response function of an unknown system, h denotes the inverse filter, s refers to the estimated fault excitation source, and indicates the convolution operation. The convolution operation for discrete signal can be described as the matrix form:where N and L are the lengths of s and h, respectively.
Considering the cyclic stationary trait of fault signal of rotating machinery, a novel blind deconvolution method called CYCBD based on second-order cyclostationarity indicator is investigated by Buzzoni et al., and its specific principles can be found in Reference [18]. A larger indicator value means the periodic impact behavior of deconvolution signal is more prominent. Thus, the target of CYCBD is to get the maximal indicator value. In the deconvolution process driven by cyclostationary maximization, the cyclic frequency is defined aswhere Ts is the cycle related to impact occurrence period.
The mathematical expression of the novel second-order cyclostationarity indicator iswith
Equations (13) and (14) can be further expressed as following matrix form:where , , , and k is the sample index.
Based on above equations, equation (12) can be expressed as
The signal which contains the periodic impact corresponding to all the cyclic frequencies of interest ka can be described as
Substituting equations (9) and (17) into (16), the expression is transformed aswhere the expression of weighting matrix is
Equation (18) is the core of the CYCBD method. By solving equation (18) through equation (19), the ultimate inverse filter h can be calculated. After substituting this calculation result into equation (8), the deconvolution signal s, i.e., the estimated fault excitation source corresponding to the maximum second-order cyclostationarity, can be extracted from the observed signal x according to the cyclic frequency.
3.2. Morphological Difference Operation Entropy
Generally, the traditional kurtosis indicator is adopted to evaluate the impact characteristic of time domain waveform. However, this indicator is very susceptible to accidental impact, which may lead to erroneous assessment [21]. Then, a new dimensionless indicator called morphological difference operation entropy (MDOE) is developed here to avoid the drawback of kurtosis indicator.
Assume x (n) is an original discrete signal with length of N points; (m) is the structural element with length of M (M ≤ N). The morphological dilation and erosion operations of x (n) using (m) are as follows [22]:where and , respectively, refer to the dilation and the erosion operators.
The opening and the closing operations can be defined as below according to the sequential operation of morphological erosion and dilation.where and , respectively, denote the opening and the closing operators. The opening operation suppresses positive impacts and retains negative ones, while the closing operation can restrain negative impacts and preserve positive ones.
The difference operation of closing and opening operations is further described as
As reported in Reference [23], the morphological difference operation has a favorable quality of strengthening the positive and the negative impacts simultaneously, and this trait has been widely applied for impulsive signature enhancement. In another aspect, Shannon entropy has the capacity of describing the uniformity of a sequence [24]. Combining the respective specialties of morphological difference operation and Shannon entropy, a novel morphological difference operation entropy (MDOE) indicator is developed as
For a given signal, the structuring element type and the morphological scale should be selected in advance before calculating its MDOE indicator. In order to guarantee the high computational efficiency and retain the impact characteristic, the simple flat structuring element with small morphological scale [0 0 0 0] is considered according to the research in Reference [25]. To clearly illustrate the feasibility of this indicator in reflecting the periodic impact characteristic of given signal, four kinds of simulated signals are constructed and shown in Figure 3(a). The sample frequency and signal length are, respectively, set to 12800 Hz and 0.1 second. The first signal S1 is created according to Reference [26], which represents the periodic impacts of bearing outer ring defect without any noise. It should be pointed out the amplitudes of the impact components are equal. Three random impacts (marked by red lines) with larger amplitudes are added to S1 to generate the second signal S2, which stands for the accidental impact interferences. The third signal of superimposed harmonics is created using equation . The fourth signal S4 denotes the white noises with standard deviation of 0.4. Figures 3(b) and 3(c) illustrate the corresponding MDOE and kurtosis indicators of four simulated signals. Compared with the other three signals, the ordered degree of periodic impacts in S1 is highest and the corresponding MDOE indicator is the minimum. The regularity of periodic impacts in S2 is affected by the several added random impacts, and thus its MDOE value increases a small value. S3 and S4 do not present periodic impact characteristic, and thus the regularities of these two signals are unobvious and they have larger indicator values.

Based on the above analysis, we can generalize that a higher ordered degree means a smaller MDOE value. Thus, this indicator is able to effectively distinguish these four signals because it can quantitatively assess the periodic characteristics of impact components, while the traditional kurtosis indicator cannot achieve this target due to its random impact sensibility. Comparing Figures 3(b) and 3(c), we can find that kurtosis value increases from 21.906 to 37.428 when the periodic impacts are disturbed by accidental impacts, and the relative augment reaches to 70.86%. Obviously, the kurtosis indicator is very sensitive to a small number of accidental impacts. On the other hand, MDOE indicator changes a tiny value, from 5.683 to 5.827, and the relative augment is only 2.53%. Therefore, the developed MDOE indicator is extremely robust to the accidental interferences.
Furthermore, different degrees of white noises are added to S1 to acquire the different SNR simulated signals. Figure 4 shows the curve of MDOE indicators corresponding to different SNRs. By comparing the indicator values under different SNR conditions, the effectiveness of this indicator in evaluating the richness of periodic impacts is further verified. It is well known the periodic impact characteristic will be more obvious if SNR is higher. As shown in the figure, the indicator value tends to decrease when the SNR improves, which means smaller indicator represents more abundant periodic impacts. This is because the impact regularity of simulated signal decreases in the presence of noise interferences. Consequently, the developed MDOE indicator is suitable for fairly assessing the deconvolution signal of CYCBD and can be applied as the guidance indicator during parameter optimization process.

3.3. Parameter Optimization Based on Grid Search Algorithm
As a kind of parametric signal processing method, CYCBD needs to set the cyclic frequency a and the filter length L in advance. And the obtained result will be unsatisfactory if these two parameters are selected unreasonably. Thus, the important problem of applying CYCBD to bearing fault detection is how to choose the influencing parameters appropriately for the best performance. We have known that the cyclic frequency a is closely related to the characteristic frequency. Nevertheless, the desired result is difficult to obtain if this parameter is directly confirmed by referring to the theoretical characteristic frequency. The reason is that the actual characteristic frequency is often inconsistent with the theoretical calculating value due to inevitable roller random slippage during bearing operation. As for filter length selection, there is no reference. Literature [27] has shown the performance of CYCBD is not proportional to the filter length. Whether the selected filter length L is too big or too small, it will not be conducive to fault characteristic extraction. To automatically acquire the optimal influencing parameters, the grid search algorithm is introduced here and fused with CYCBD to achieve this purpose.
The grid search algorithm divides the ranges of independent variables into certain spaces and in turn calculates the objective function of each grid point under the constraint conditions [19]. And the combination of independent variables making the objective function optimal is deemed as the final output. In essence, parameter optimization of CYCBD is a nonlinear programming problem. The developed MDOE indicator can be utilized as the objective function to guide the whole search process. And objective function minimum is considered as the final target in grid search algorithm. Then, the corresponding expressions can be constructed as follows:where a and L are the independent variables, i.e., the cyclic frequency and the filter length, au and ad refer to the upper and the lower boundaries of cyclic frequency, and Lu and Ld stand for the upper and the lower boundaries of filter length.
Then, the constraint conditions can be imposed by setting the boundary ranges of independent variables; it should be pointed out the grid precision and the search efficiency are determined by the boundary ranges and the step sizes. Assume that the step sizes for searching a and L are, respectively, z and q. Then, the search scopes of cyclic frequency and filter length are severally divided into (au − ad)/z + 1 and (Lu − Ld)/q + 1 sections. On the one hand, the search precision will be very high and the calculation burden will be very heavy if the boundary ranges are chosen too large and the step sizes are set too small. On the other hand, the search efficiency will rise while the precise search result is easy to miss in case these parameters are set inversely. The MDOE indicator of each grid point is in turn calculated and the location corresponding to the minimum value is the optimal parameters, i.e., (ao, Lo) = argmin (MDOE). Taking advantage of gird search algorithm, the influencing parameters of CYCBD can be precisely confirmed. Then, the OCYCBD method is put forward to automatically get the optimal deconvolution result, and its characteristic extraction and amplification performances are given expectations.
4. Implementation Procedures of the Proposed RSSD-OCYCBD Strategy
In actual engineering application, the bearing injury signatures are usually inconspicuous and covered by strong external noises in the early stage. For a complicated observed signal of rolling bearing, the SNR can be improved to some extent and the weak fault characteristics are able to be discovered more easily if the original signal can be separated into several simple components. Owing to the excellent ability of decomposing the multicomponent signal, RSSD is regarded as a preprocessing method to deal with the original signal. Nevertheless, a lot of noise interferences still remain in the obtained results because the SNR of original signal is very low. Thus, it is difficult to achieve the ideal fault recognition target. CYCBD is capable of eliminating the spreading effect of unknown transmission path and intensifying the impulsive characteristics. Moreover, compared with the traditional deconvolution algorithm, it is more robust against interference noises. Thus, for the purpose of getting rid of the redundant interferences and enhancing the fault signatures, the OCYCBD method, whose influencing parameters are adaptively selected by grid search algorithm, is further applied as a postprocessing approach. On the basis of these statements, a novel enhanced fault detection strategy based on RSSD-OCYCBD is put forward to improve the accuracy of defect identification for rolling bearing. Figure 5 illustrates its schematic, and the detailed implementation procedures are as follows: Step 1: data acquisition. The original signal is collected using the corresponding data acquisition equipment. Step 2: RSSD preprocessing. The proposed RSSD method in Section 2 is adopted to adaptively decompose the original signal into a number of components. On the one hand, noise restrain can be achieved through this operation. On the other hand, the decomposition results are conducive to determining the search center for cyclic frequency parameter optimization. Step 3: sensitive SSC selection. If the SCC component contains more useful information, then the peak spectral line related to the theoretical defect characteristic frequency will be more prominent in the corresponding envelope spectrum. Therefore, to quantitatively compare the obtain SSC, the envelope spectrum of each component is calculated. And envelope spectrum peak (ESP) index used in Reference [28] is considered as the criterion to select the sensitive SSC; the definition of this index is where eRMS refers to the root mean square value of envelope spectrum and max (A (f)) denotes the maximal amplitude of frequency component in envelope spectrum. For the sake of avoiding the interferences of unrelated spectral lines as far as possible, only the frequency areas around the theoretical fault characteristic frequencies are analyzed, and the corresponding formula expressions are where fi, fo, fe, and fc, respectively, stand for the theoretical bearing fault characteristic frequencies of inner ring, outer ring, roller, and cage. In this paper, the analysis region interval ΔI is set to 5 Hz. The most sensitive SSC, which contains more obvious fault signatures than original signal, can be selected according to the largest ESP index. However, considering that the component also has a lot of noise interferences, the ideal diagnosis effect is difficult to obtain. Then, the following procedures are carried out to eliminate the redundant interferences and enhance the fault signatures. Step 4: search center of cyclic frequency confirmation. The cyclic frequency parameter is difficult to precisely select due to the wide frequency domain scope. Thus, for the purpose of improving the search efficiency and easing the computational burden, the search center of cyclic frequency needs to confirm, and the subsequent grid search algorithm is carried out around this center. The envelope spectrum analysis of the most sensitive SSC is beneficial to the realization of this goal. On the basis of the illustrated contents in step 3, the spectral line corresponding to the largest amplitude among the theoretical fault characteristic frequency regions is determined as the central value of cyclic frequency search range. And the corresponding expression can be described as ac = argmax (A (f)) based on equation (26). Considering the fence effect of envelope spectrum analysis, this central value may not be the optimal cyclic frequency parameter, and thus more elaborate grid search process around this center is carried out in the next step. Step 5: search scope and step size initialization. In order to optimize the parameter a and L, a two-dimensional grid coordinate system is constructed. The search scope of cyclic frequency is initialized as [ad, au] = [floor (ac − Δf), ceil (ac + Δf)] with step size of 0.1, where denotes the frequency resolution of envelope spectrum of sensitive SSC and floor and ceil, respectively, refer to the round down and the round up operations. As for the filter length L, the deconvolution signal will be distorted if it is too large, while the processing effect will be inconspicuous if this parameter is too small. Referring to literature [29], its search scope is set to [Ld, Lu] = [20, 200] with step size of 1. Step 6: optimal influencing parameters selection. Firstly, the cyclic frequency is fixed at the lower boundary a = ad, while the filter length L continuous increases from the lower boundary Ld to the upper boundary Lu with step size of 1. Then, the cyclic frequency changes with step size of 0.1 to the next location a = ad + 0.1, and filter length L still continuously increases from Ld to Lu with step size of 1. This similar process is executed until the cyclic frequency reaches the upper boundary a = au. In each grid point location, the corresponding current influencing parameters are, respectively, substituted into CYCBD to perform deconvolution operation on the most sensitive SSC, and the MDOE indicator of the obtained deconvolution signal is calculated. Then, all of the grid point locations are evaluated according to MDOE values and the parameters corresponding to the minimum indicator are regarded as the optimal outputs. Step 7: OCYCBD postprocessing. The obtained optimal parameters ao and Lo in step 6 are substituted into CYCBD, and OCYCBD is further adopted on the most sensitive SSC to enhance the fault characteristics. Step 8: defect identification. The envelope demodulation operation is carried out on the obtained deconvolution signal, and the incipient injury of rolling bearing is able to be judged by identifying the peak frequencies in the envelope spectrum.

5. Simulated Signal Verification
To manifest the feasibility of proposed enhanced detection strategy, it is assumed that a local flaw occurs on the inner ring of rolling bearing, and a simulated fault signal comprising periodic impact components and Gaussian white noises is constructed as [30]where refers to the convolution operation, k is the positive integer, represents the Dirac delta function, and c = 800 rad/s stands for the structural attenuation factor of bearing system. fr = 30 Hz and fn = 2500 Hz, respectively, represent the rotating frequency and the resonance frequency excited by defect point strike. D = 0.3 refers to the impact amplitude fluctuation caused by rotating frequency modulation, and fault characteristic frequency of inner ring is set to fi = 140 Hz, which is equal to the reciprocal of impact period T. R (b, k) = bT denotes the random fluctuation of expected repetitive period of kth impact due to possible rotating speed variation or roller slippage, which usually accounts for 1%∼5% of T. Then, in this paper, b is set to the random number in the range of [0.01 T, 0.05 T]. Furthermore, the sample rate and sample length are set to fs = 12800 Hz and N = 8192 points, respectively. n (t) denotes the added noises with standard deviation of 1, and the SNR of simulated signal is −13 dB, which can be calculated by the formula in Reference [31].
Figure 6 shows the waveforms as well as the corresponding spectra of simulated inner ring fault signal. After adding the Gaussian white noises, the periodic impacts cannot be found in the simulated signal. Meanwhile, the resonance frequency band is indistinguishable in the frequency spectrum. The peak spectral lines related to characteristic frequency also cannot be found in the envelope spectrum.

(a)

(b)

(c)

(d)
Analysis results of simulated signal by proposed strategy are shown in Figure 7. Four components are acquired automatically after RSSD preprocessing, and the waveforms of components are shown in Figure 7(a). The corresponding ESP indexes of four components are calculated as shown in Figure 7(b). Among the obtained decomposition results, the ESP index of SSC2 is the largest value (9.21). Figure 7(c) further shows the frequency spectrum of SSC2, which manifests that this component is exactly located on the setting resonance frequency band. It is indicated that SSC2 is the most sensitive component and contains the most abundant characteristic information. Thus, envelope demodulation is carried out on this component and the obtained envelope spectrum is obtained in Figure 7(d). In this figure, we can identify the spectral peaks of 140.6 Hz, 279.7 Hz, and 700 Hz, and these frequencies suspiciously correspond to the characteristic frequency fi, its second harmonic 2fi, and fifth harmonic 5fi. For reasons of the random fluctuations of periodic impact intervals and the fence effect of envelope spectrum analysis, the identified characteristic frequencies deviate from the theoretical setting values, and this is inevitable in actual engineering application. In order to present the fault signatures more clearly, OCYCBD needs to be applied as a postprocessing method. During the process of CYCBD parameter optimization, the 140.6 Hz peak frequency is regarded as the search central of cyclic frequency based on envelope spectrum analysis result of SSC2. As the frequency resolution of envelope spectrum Δf = fs/N = 1.5625 Hz, the search scope of cyclic frequency is set as [ad, au] = [floor (140.6–1.5625), ceil (140.6 + 1.5625)] = [139, 148] with step size of 0.1. As illustrated in Section 4, the boundary scope of filter length is set as [Ld, Lu] = [20, 200] with step size of 1. Thus, the grid point locations of cyclic frequency and filter length are, respectively, divided into 91 and 181 sections. And the obtained two-dimensional grid coordinate system is shown in Figure 7(e), based on which we can confirm the optimal cyclic frequency ao = 139.8 and the optimal filter length Lo = 50. Then, CYCBD with these two parameters is executed on SSC2; the waveform and the corresponding envelope spectrum of obtained deconvolution signal are shown in Figures 7(f) and 7(g). And the characteristic frequency and its harmonics fi∼4fi are very outstanding in the envelope spectrum.

To verify the necessity of RSSD pretreatment, CYCBD with optimal parameters is directly applied to the original signal, and the obtained results are shown in Figure 8. It can be found that there exists a great gap between Figures 7 and 8 by comparing the waveform and the envelope spectrum of deconvolution signal. Unfortunately, we cannot find any useful information in Figure 8. Without RSSD preprocessing, the deconvolution effect is seriously affected by heavy background noises. Then, it is clearly verified the advantage of fusing RSSD with OCYCBD on noise interference suppression and fault characteristic enhancement.

(a)

(b)
To further illustrate the superiority of gird search algorithm for precise parameter selection, the obtained SSC2 is, respectively, processed by CYCBD under the conditions of the setting parameters deviating from the optimal search values. Firstly, the filter length is fixed at L = 50, and we directly select a = 140.6, which varies within 1% from the optimal cyclic frequency ao = 139.8. Then, the cyclic frequency is fixed at a = 139.8, and the filter length is changed to L = 70. The obtained comparison results are, respectively, illustrated in Figures 9(a) and 9(b), and the characteristic frequency components are disturbed by background noises and interference spectral lines, and thus the obtained results are bad for achieving precise bearing flaw detection. It is also indicated that the deconvolution performance of CYCBD is seriously dependent on the set key parameters and the characteristic extraction effect decreases rapidly due to the changes of key parameters. The obvious difference between Figures 9(a) and 9(b) shows the performance of CYCBD is more easily influenced by the cyclic frequency compared with the filter length; then, this parameter must be set as precise as possible to get the favorable result. Fortunately, we can achieve this purpose with the help of gird search algorithm.

(a)

(b)
6. Experimental Verification
6.1. Introduction of Experimental Platform
The collected vibration signal from experimental platform is utilized for feasibility verification. Figure 10 shows 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, i.e., 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 10(a). The injury was inserted on the roller of experimental bearing, as shown in Figure 10(b). The width and the depth of local defect were, respectively, 0.2 mm and 1.53 mm. The bearing structure parameters are listed in Table 1. The installed PCB acceleration sensors on the bearing seat are shown in Figure 10(c), and the signal acquisition process is shown in Figure 10(d). During the process of bearing running, the sample frequency was set as fs = 12800 Hz and the signal with length of 8192 points is intercepted for analysis. The defect characteristic frequencies of experimental bearing are, 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 Reference [32].

(a)

(b)

(c)

(d)
Figure 11 shows the waveform of experimental signal as well as its spectra. There exists obvious impact phenomenon in the waveform, while the interval between adjacent impacts is not the reciprocal of any fault characteristic frequency. The frequency components are mainly concentrated below 2500 Hz, and there exists broad resonance band in the frequency spectrum. However, the spectral line related to bearing defect cannot 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 characteristic frequency component in the envelope spectrum. Thus, the bearing injury could not be detected and the obtained envelope spectrum is invalid.

(a)

(b)

(c)
6.2. Experimental Signal Analysis and Result Comparison
The presented strategy is applied for analyzing this experimental signal. Firstly, four components are acquired by RSSD processing, and Figure 12(a) shows their waveforms. Figure 12(b) shows the corresponding ESP indexes. Among the obtained decomposition results, the ESP index of SSC3 is the largest, which means this component is the most sensitive component and contains the most abundant characteristic information. Then, the envelope spectrum of this component is further calculated and shown in Figure 12(c). As marked in this figure, the spectral peaks of 115.6 Hz and 231.3 Hz are suspected to be the roller fault characteristic frequency and its second harmonic. Based on the selection principle of search center of cyclic frequency, which has been detailed expounded in Section 4, the peak frequency 115.6 Hz among the theoretical fault characteristic frequency regions is confirmed as the search central ac of cyclic frequency. Due to the frequency resolution of envelope spectrum Δf = fs/N = 1.5625 Hz, the lower and the upper boundaries of cyclic frequency are ad = 114 and au = 118. And the lower and the upper boundaries of filter length are Ld = 20 and Lu = 200. The corresponding grid search result is shown in Figure 12(d); we can confirm the optimal cyclic frequency ao = 115.8 and the optimal filter length Lo = 140. Then, these two parameters are substituted into CYCBD to further process SSC3. The waveform and the envelope spectrum of deconvolution signal are, respectively, presented in Figures 12(e) and 12(f). Owing to the favorable specialty of OCYCBD, a lot of interference components in the sensitive component are eliminated and the fault related spectral lines fe∼6fe with large amplitudes are able to be presented more obviously. Thus, it is indicated that there exists local flaw in the roller, and this conclusion is consistent to the actual condition. Then, the feasibilities of the proposed strategy in getting rid of the severe background noises and intensifying the incipient fault signatures of bearing is verified.

To verify the superiority of RSSD pretreatment, CYCBD with optimal parameters is directly adopted to process the original experimental signal, and the obtained results are shown in Figure 13. It can be found the analysis results are not as ideal as those in Figure 12. Some interference spectral lines appear in the envelope spectrum, which hinders the characteristic frequency identification. In addition, the CYCBD parameter selection process will be more troublesome without RSSD preprocessing. Owing to the low SNR, the suspectable fault frequency cannot be detected based on envelope spectrum analysis of original signal. Then, the search center and the search range of cyclic frequency cannot be confirmed quickly. Therefore, it will bring inconvenience for the whole diagnosis progress.

(a)

(b)
The current methods of MCKD, MED, MOMEDA [33], EWT [34], and SK [35] are also applied for comparison. The analysis results using MCKD, MED, and MOMEDA are shown, respectively, in Figures 14–16. The influencing parameters of these three methods are confirmed based on the grid search result. As the optimal cyclic frequency ao = 115.8 and the optimal filter length Lo = 140, the deconvolution period of MCKD and MOMEDA is calculated as fs/ao = 110.535, and the filter lengths of these three methods are set to 140. Regrettably, the characteristic frequency information can be found in Figures 14(b), 15(b), and 16(b) due to the SNR of original signal being very low. The waveforms and the corresponding envelope spectra of the separated 6 components obtained by EWT are shown in Figure 17. The characteristic frequency and its harmonics are able to be detected in the envelope spectra of C5 and C6. It is noticeable that the fault signatures in these two envelope spectra are not as obvious as those in Figure 12(f); the numbers of characteristic frequency harmonic are fewer and some unrelated components simultaneously appear. As for SK analysis results, Figure 18(a) illustrates the calculated kurtogram, and the optimal resonance band in level 6 is marked by black oval. And a bandpass filter is designed based on the kurtogram. However, the intervals of periodic impacts of filtered signal in Figure 18(b) are not related to bearing defect, and the characteristic frequency is not able to be extracted from the envelope spectrum in Figure 18(c). Owing to the strong noise interferences, the obtained information for filter design is unreasonable and the passband of the constructed filter does not cover the true resonance region. Therefore, the SK method is invalid for analyzing this experimental signal. Based on above analysis results, it is verified that the capacity of the proposed strategy on weak fault detection is superior to these comparison methods.

(a)

(b)

(a)

(b)

(a)

(b)

(a)

(b)

(a)

(b)

(c)
7. Engineering Case Verification
7.1. Description of Wind Turbine
An engineering case of generator bearing injury of a wind turbine, whose nominal output power is 750 KW, is used to verify the proposed method. Figure 19(a) shows the schematic diagram of the wind turbine, which is composed of an impeller, NGC-FD825 gearbox, and double-fed asynchronous generator. In order to collect the vibration signal of the generator bearing, the RH103 pass-frequency acceleration sensors with sensitivity of 100 mV/g were installed on the front-back bearing housings of the generator through the magnetic base connection, as shown in Figures 19(b) and 19(c). The sampling frequency of the data acquisition system was set as 16384 Hz. The generator front bearing vibration signal with length of 8192 points is analyzed here. And the bearing type is SKF6324, whose structure parameters are listed in Table 2. The speed of the main shaft was = 21.63 rpm, and transmission ratios of the planetary gears, the middle stage gears, and the high speed stage gears in the gearbox are 4.98, 3.81, and 3.7, respectively. Thus, the speed of the generator shaft was = 1519 rpm, and the corresponding rotating frequency was = 25.3 Hz. The theoretical outer race characteristic frequency, the inner race characteristic frequency, the rolling element characteristic frequency, and the cage characteristic frequency of the generator bearing are, respectively, calculated as = 79.21 Hz, = 123.18 Hz, = 110.96 Hz, and = 15.39 Hz.

(a)

(b)

(c)
Figure 20 describes the waveform and the spectra of collected engineering signal. The periodic impact components in the waveform are not abundant. And the kurtosis index of this engineering signal, which is used to evaluate the impulsive phenomenon, is only 3.15. 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 characteristic frequency of generator bearing can be found. And we cannot make any conclusion based on limited information.

(a)

(b)

(c)
7.2. Engineering Signal Analysis and Result Comparison
The proposed strategy is then applied to this engineering signal. Figure 21(a) shows the obtained components by RSSD. And the calculated ESP indexes show the third component SSC3 is the most sensitive component, and its envelope spectrum is shown in Figure 21(c). Compared with the original signal, SSC3 contains more plentiful useful information. Some outstanding impacts appear in the waveform, and a suspectable fault frequency spectral line of 80 Hz is visible in the envelope spectrum. However, in engineering application, the multiple frequency components of characteristic frequency, which contributes to make a more accurate decision, are expected to appear. Based on the confirmation principle of search center of cyclic frequency, the peak frequency 80 Hz is considered as ac. As the frequency resolution of envelope spectrum Δf = (16384/8192) Hz = 2 Hz, the cyclic frequency is optimized within the boundary scope [78, 82] and the filter length is selected within the boundary scope [20, 200]. As shown in Figure 21(d), the optimal outputs are ao = 79.8 and Lo = 170. Then, CYCBD with these two parameters are carried out on SSC3. The periodic impacts in Figure 21(e) are richer than SSC3. It means the SNR of this component has been improved by OCYCBD postprocessing. Figure 21(f) shows the envelope spectrum of deconvolution signal, and peak frequencies fo∼4fo corresponding to outer ring fault are clearly presented. Thus, it is highly doubtful that there is local defect on the outer ring of the generator front bearing, and the boarding check results verify this conclusion. Thus, the effectiveness of presented strategy in achieving accurate weak fault detection for rolling bearing is demonstrated.

The analysis results of CYCBD with optimal parameters on the original engineering signal are shown in Figure 22. It is found that the periodic shocks in Figure 22(a) are not as abundant as those in Figure 21(e). In addition, some interference components exist in Figure 22(b), and the harmonics of characteristic frequency are unobvious. Based on the comparison results, there is no doubt that the RSSD pretreatment is more beneficial for fault identification.

(a)

(b)
The MCKD, MED, MOMEDA, EWT, and SK methods are also performed on the same engineering signal. Figures 23–25, respectively, illustrate the obtained results using MCKD, MED, and MOMEDA. The basic characteristic frequency and its second harmonic are visible in the corresponding envelope spectra of deconvolution signals. The analysis results by EWT are shown in Figure 26; the basic characteristic frequency can be identified in C8 envelope spectrum. Nevertheless, the multiple frequency components are not rich enough and the acquired results are not as satisfactory as that in Figure 21(f). The kurtogram is depicted in Figure 27(a), and the bandpass filter is constructed in terms of kurtogram information to process the engineering signal. Figures 27(b) and 27(c) show the waveform and the envelope spectrum of filtered signal; only a weak basic characteristic frequency fo can be detected, and it is clearly not an ideal result. Based on these results, it is indicated that the capacities of the proposed RSSD-OCYCBD strategy on characteristic extraction and enhancement for rolling bearing diagnosis are better than the comparison methods.

(a)

(b)

(a)

(b)

(a)

(b)

(a)

(b)

(a)

(b)

(c)
8. Conclusions
An original enhanced fault detection strategy called RSSD-OCYCBD is presented to identify the local defect of rolling bearing in the early injury stage. Within this detection strategy, the highlights are summarized as follows:(1)The novel RSSD method is put forward to overcome the drawback of component number artificial selection in the original SSD algorithm, and then the automatic signal decomposition can be achieved. The noise interferences can be suppressed to some extent, and the most sensitive SSC with higher SNR can be separated from the original signal through RSSD processing.(2)The new MDOE indicator combining the respective specialties of morphological difference operation and Shannon entropy is investigated to evaluate the quality of deconvolution signal during parameter optimization process, and the advantage of this indicator in quantitatively reflecting the richness of periodic impact components has been compared and verified by different simulated signals.(3)The two-phase optimization scheme based on RSSD and grid search algorithm is proposed to solve the problem of precise influencing parameter selection in the original CYCBD method. And the optimal deconvolution result can be automatically obtained using the OCYCBD method.
The simulated signal, the experimental signal, and the engineering field signal are, respectively, applied to demonstrate this proposed strategy. And the analysis results show that it can effectively improve the diagnosis effect and accurately identify the weak defect of the wind turbine bearing. Compared with some widely used diagnosis methods, the comparison results verify the advantages of proposed enhanced diagnosis strategy on background noise suppression and weak fault characteristic intensification. It is expected that this strategy can be applied in a broader field. And our future research target is to modify this method for solving the compound fault problem of rolling bearing under variable speed condition.
Abbreviation
CYCBD: | Cyclostationary blind deconvolution |
OCYCBD: | Optimized cyclostationary blind deconvolution |
SSD: | Singular spectrum decomposition |
RSSD: | Recursive singular spectrum decomposition |
SSC: | Singular spectrum component |
MDOE: | Morphological difference operation entropy |
ESP: | Envelope spectrum peak |
EMD: | Empirical mode decomposition |
LMD: | Local mean decomposition |
NMSE: | Normalized mean square error |
MCKD: | Maximum correlated kurtosis deconvolution |
MED: | Minimum entropy deconvolution |
SNR: | Signal-to-noise ratio |
SK: | Spectral kurtosis. |
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (grant no. 51777074), the Natural Science Foundation of Hebei Province, China (grant no. E2019502047), and the Fundamental Research Funds for the Central Universities of China (grant nos. 2018MS124 and 2017MS190).