Utilisation of Ensemble Empirical Mode Decomposition in Conjunction with Cyclostationary Technique for Wind Turbine Gearbox Fault Detection

: In this paper the application of cyclostationary signal processing in conjunction with Ensemble Empirical Mode Decomposition (EEMD) technique, on the fault diagnostics of wind turbine gearboxes is investigated and has been highlighted. It is shown that the EEMD technique together with cyclostationary analysis can be used to detect the damage in complex and non-linear systems such as wind turbine gearbox, where the vibration signals are modulated with carrier frequencies and are superimposed. In these situations when multiple faults alongside noisy environment are present together, the faults are not easily detectable by conventional signal processing techniques such as FFT and RMS.


Introduction
With the wind energy global installed capacity reaching 651 GW in 2019 [1], effective remote condition monitoring and predictive maintenance strategies for wind turbines are becoming more significant attracting increased attention. Wind turbines are typically expected to have a design lifespan of 20-25 years. However, they rarely meet this target without major overhauls, particularly in their drivetrain [2]. Gearbox faults and failures take considerable time to repair. Hence, they result in major downtime and loss of production capacity leading to increased maintenance costs for wind farm operators [3][4][5][6].
Direct drive wind turbine designs remove the need for a gearbox. However, the power converters of direct drive wind turbines are among the most frequent failing components which cumulatively result in excessive downtime [7,8]. Thus, the failure rate of power electronics in direct drive wind turbines greatly exceeds the failure rate of the gearbox in geared wind turbines [8]. Nonetheless, due to the advantages and disadvantages of each design, technology forecasts suggest that both geared and direct drive wind turbines will continue to be used in the future [9].
Unlike simple gearboxes used in conventional steady-state machinery, the vibration signals recorded from the gearbox of a wind turbine can be contaminated resulting in erroneous interpretation of the data acquired. Wind turbine gearboxes have very complex designs, comprising many different gears and bearings resulting in the vibration of different parts of the gearbox to superimpose in the recorded vibration signal with multiple frequency and amplitude modulations occurring. Hence, it is very challenging to detect a developing fault in its early stages using conventional signal processing techniques such as Fast Fourier Transform or RMS [10]. Moreover, once a fault is detected it is extremely important from a maintenance point of view to evaluate the severity. Identifying a defect alone is not sufficient for a wind farm operator in order to decide if maintenance is required. In addition, if it is deemed that maintenance is required it is very important to identify when it should be carried out without risking loss of production and expensive outages. Gearboxes are rarely available as spare parts and therefore it can take several months before a spare gearbox can become available. Therefore, it is critical for wind farm operators to pinpoint the exact gearbox component affected by the fault and quantify its severity in order to accurately plan maintenance without the risk of loss of production.
The cyclostationary signal processing technique is a very powerful tool for fault detection in rotating machinery. Moreover, it can be utilised to evaluate the severity of developing faults in rolling elements and bearings. However, as it has been suggested by Dong and Chen [11], this technique is not very successful in complex systems when multiple gears and bearings are present.
Here we have investigated the suitability of combining two advanced signal processing techniques for the purpose of fault detection in the gearbox of wind turbines. These techniques are known as Cyclic Spectral Analysis and Ensemble Empirical Mode Decomposition (EEMD), which will be explained in more details in the next section. It should be stressed, that wind turbines operate in variable loading conditions and therefore, the accurate evaluation of the severity of faults is very challenging [12]. Moreover, false indications of non-existent faults are not uncommon when using conventional signal processing techniques.

Cyclic Spectral Analysis
The concept of cyclostationary signals has been around for almost 40 years [13,14]. In recent years, the theory of cyclostationary signals has been increasingly used in the field of mechanical fault detection. Vibration signals extracted from bearings and rotating machines are non-stationary signals containing some hidden periodicities in their background. To extract the periodicity many signal processing tools are available. The cyclostationary signal analysis technique is one of the most powerful techniques from those currently available [15]. The statistics of cyclostationary signals have periodicity with respect to time [10]. Currently, second-order cyclostationarity which is known as cyclic correlation [10] has an important role in practical application. Despite the strengths of this technique and its popularity in signal processing of rotating machinery, it has not attracted much attention in the field of wind turbine condition monitoring and signal processing. To describe this technique we start with a random process.
In general, a random process will have a time-varying autocorrelation: where E{. . . } denotes the mathematical expectation operator, and τ is the time delay. If the autocorrelation is considered periodic with a period of T 0 , the ensemble average can be estimated with a time average as below: The autocorrelation function in Equation (1) can be written in form of Fourier series due to its periodicity where α = m/T 0 and m ∈ Z. Combining the above equation with the Equation (2), the Fourier coefficients can be written as below where R α x (τ) is known as cyclic autocorrelation function with α being the cyclic frequency. The . . . symbol represents the time averaging operator. Next, to obtain the cyclic spectrum, the Fourier transform can be applied to the cyclic autocorrelation function with respect to the time delay τ, yielding that is known as cyclic spectrum or the spectral correlation function. In this work the Welch's averaged periodogram technique has been used, that due to its high computational efficiency is one of the most common estimators for cyclic spectrum [16].
The cyclic coherence that measures the strength of the correlation between spectral components distanced by cyclic frequencies can be calculated using [16] C α

Empirical Mode Decomposition
The Empirical Mode Decomposition (EMD) is the fundamental part of a technique known as Hilbert-Huang Transform that is an empirical data analysis technique, developed in late 1990s by Huang et al. [17][18][19]. Unlike other techniques such as Fourier or Wavelet Transforms, that use prior knowledge or fixed basis to decompose the signal, EMD derives its basis adaptively from the data itself and does not rely on any prior knowledge [20]. This technique is based on an assumption that any data is consisting of a set of simple intrinsic modes of oscillations [21].
The EMD decomposition process, decomposes the original signal into a set of signals known as Intrinsic Mode Functions (IMFs). This is accomplished by a novel process known as sifting, which is repeatedly applied to the signal until it converges on criteria that define an IMF [20]. The EMD sifting process extracts the fastest varying component first and continues to the next IMF. The process produces a finite number of IMFs before converging itself.
The extracted IMFs represent simple oscillatory functions with different time-scale. Each IMF by definition will have the same number of extrema and zero-crossings or differ at most by one. Furthermore, the mean value of the envelopes of local maxima and local minima of IMFs are zero at any given point [21]. The technique to calculating the IMFs are briefly explained here but for more details and derivations on the technique one can refer to [17,21].
The first step to calculate the IMF is to find all the local maxima points of the signal and connect them using cubic spline to obtain the upper envelope of the signal. The same procedure is applied to the local minima points to obtain the lower envelope. Next, the mean value of the two envelopes is calculated. The first component of the signal is calculated by subtracting the obtained mean from the original signal where x(t) is the original signal and m 1 is the calculated mean. If the calculated component h 1 meets the IMF criteria, then it will be the first IMF. However if the calculated component does not meet the requirement of an IMF, the procedure is repeated again by replacing the original signal x(t) with h 1 . This cycle is repeated until the calculated component meets the IMF criteria. It is then marked as the first IMF (c 1 ) and subtracted from the original signal.
x new (t) = x(t) − c 1 (t) (8) After this step, the original signal x(t) is replaced with the new signal x new (t) and the whole procedure is repeated to calculate the next IMFs. The remaining signal after calculating the last IMF is called the residual signal r(t) and can be used to reconstruct the original signal using the calculated IMFs using where n represent the number of calculated IMFs. This process is summarised in Figure 1.
Find the local exermas of h i(k−1)

Obtain upper and lower envelopes
Calculate the mean m i(k−1) is h ik an IMF?
Stop no yes yes no Figure 1. The procedure to obtain Intrinsic Mode Functions (IMFs) using Empirical Mode Decomposition (EMD) technique, adapted from [22]. Refer to the text for more details.
Although this technique is very effective in decomposing non-linear and non-stationary signals, it suffers from certain problems such as end effect or mode mixing (aliasing) which was noted by Huang himself [19]. One solution to avoid the mode mixing problem is to use a newly proposed technique known as Ensemble EMD (EEMD) [23]. This technique uses noise assisted analysis to solve the problem of mode mixing in the EMD algorithm.
This can be done by adding White Gaussian noise to the signal. Consider adding noise to the proposed signal for M times and each time decomposing the signal using EMD technique. Each new signal can be presented as where w i (t) represents the i-th added white noise to the original signal x(t) and x i (t) is the generated noisy signal. Using Equation (9) we can reconstruct the i-th noisy signal as where c ij is the j-th calculated IMF of noisy signal x i (t) and r i (t) is the residual signal. Now we can ensemble the corresponding IMF over noisy signals to obtain the resulting IMFs where M is the total number of times that the white noise was added to the signal. The accuracy of the result obtained using EEMD is highly dependent on the number of ensemble and the amplitude of the added noise (ε). Using a well-established statistical rule such as Equation (13), the standard deviation of error (ε n ) in the final IMFs can be calculated [23].

Measurement Setup
The evaluation of the proposed technique was done using three different testing setups. Initially, a lightly contaminated signal recorded from a simple rolling element bearing was used to evaluate the performance of the cyclostationary analysis, as well as its combination with EEMD.
Next, a more complex scenario was considered with a signal recorded from an input shaft of a 3-stage gearbox in a mock-up wind turbine. Finally, the data recorded from the gearbox of a real wind turbine, during full load operation, was used to evaluate the performance of the technique. Each of these three setups is detailed below.

Rolling Element Test Rig
A simple customised test rig shown in Figure 2 was used for bearing test at the University of Birmingham. The test rig consists of an AC motor rotating a shaft via pulley and belt. A controlled load was applied to the rolling element under investigation. The speed of the bearing was monitored using a tachometer. The sensor which was used to record the vibration signal was a high frequency 712F accelerometer by Wilcoxon Corporation with a nominal range of 3 Hz to 25 KHz.

Tidal Turbine Gearbox Test Rig
The test rig used for these measurements has been developed by TWI Ltd. in the UK and the University of Brunel in an effort to replicate a geared tidal turbine with the ability to be submerged into the water. The test rig consists of a 3-phase 750 Watt AC motor controlled by a Variable Frequency Driver (VFD) capable of running at 1500 rpm, a 3-stage spur and helical parallel gearbox with a ratio of approximately 1:69 and 3 blades at the main shaft turning at approximately 22 rpm. The tidal turbine gearbox test rig setup is shown in Figure 3. This test rig is instrumented with two accelerometers in X and Y directions mounted on the input shaft of the gearbox (i.e., slow speed shaft) over the bearing to record vibration data of the gear box. The sensors used are AC 150-2C accelerometer by Connection Technology Centre Inc with a nominal range of 1 Hz to 10 kHz.

Wind Turbine Gearbox
Vibration data were recorded from an operational industrial wind turbine at full load using an INGESYS system manufactured by INGETEAM Service, Spain. The monitored wind turbine had a rated power of 850 kW with a 3-stage planetary gearbox and a ratio of approximately 1:62. The wind turbine was known to be due for maintenance and it had some gear teeth wear and bearing damage at the planetary stage. The INGESYS monitoring system was installed in an industrial cabinet in the nacelle of the wind turbine. It consists of the power supply, communications and vibration acquisition system. During testing it was connected to the ground data recording unit using fibre optic.
In order to measure vibrations in the drive-train, 8 piezoelectric accelerometers were used. Two of them were specifically selected to measure low vibration frequencies in the main shaft. The two types of the accelerometer sensors that were used are PRE-1010-MS and PRE-1030-MS by IMI Sensors capable of covering the range of 0.5 Hz to 10 kHz and 2 Hz to 10 kHz, respectively. These sensors were installed on both radial and axial directions on the main shaft, planetary stage, high speed shaft and the generator. Some of the accelerometers installed on the wind turbine gearbox are shown in Figure 4.

Accelerometer Gearbox
Radial Accelerometer Generator Axial Accelerometer

Rolling Element Analysis
Two identical tapered roller bearing with different induced faults were used in this setup. The first roller bearing was damaged on its both outer races and the second roller bearing was scuffed on the roller surface as shown in Figure 5. Each scuffed roller was placed in one of the bearing races. The recorded signal from these bearings was also expected to be contaminated with some noise originating from the AC motor, pulley, belt and other nearby bearings used in the test-rig itself. The characteristic frequencies of the rolling element used for the experiments using the customised test rig at the University of Birmingham are summarised in Table 1.  In each case, it is expected to observe the presence of the bearing characteristic frequency corresponding to the induced fault. In the case of the first bearing, the Ball Pass Frequency of Outer race (BPFO) is expected to be dominant. For the second bearing, the Ball Spin Frequency (BSF) that indicates the frequency of which the roller is turning can be used to identify the fault. However, as the rollers turn, the scuffed area will hit both the inner and outer race in each revolution. Thus, the observed defect frequency is expected to be twice the BSF frequency.
The Cyclic Spectral Coherence was calculated over 4 s of the vibration data. As it can be seen in Figure 6, the BPFO frequency is clearly visible. Alongside the BPFO frequency, some modulations at the BSF frequency are also present which indicate damage in the rolling elements. This was confirmed by inspection of the rollers after testing. The outer race damage and heavy load on the bearing caused the surface of the rollers to degrade. However, the degradation was not significant.  As it was mentioned earlier, the fault frequency (2 × BSF) was observed in the Cyclic Spectral Coherence. However, unlike the first bearing frequency, it was more spread across the frequency spectrum. This could have been caused due to two damaged rollers or slight change in the bearing speed during the test arising from the load imposed.
As it can be seen in Figure 7, the fundamental train frequency (4.4 Hz) is also very dominant. This was assumed to be due to the interaction between the damaged roller and the cage. However, after inspection, it was found that one of the damaged rollers led to fracture of the cage of the second bearing.

Simplified Gearbox Analysis
A simple gearbox operating at steady-state was studied with imbalance caused in the rotor i.e., slow shaft, and the high speed shaft rotating at 750 rpm. The characteristic frequencies of the gearbox are summarised in Table 2. The Cyclic Spectral Coherence was calculated using 4 s vibration signals. The calculated results are presented in Figure 8. As the induced fault was imbalance, modulations and cyclic behaviour at the gear-mesh (GM) frequency of the first stage of the gearbox (10.5 Hz) were caused. As it can be seen in Figure 8a, the result from the analysis of the original raw signal is less clear in comparison with the bearing test. For comparison purposes, the analysis was also applied to the de-noised signal using the Wavelet de-noising technique. This was carried out using Daubechies-3 wavelet. However, de-noising the signal worsened the Cyclic Spectral Coherence result as it is shown in Figure 8b.
Using the EEMD technique, the signal was decomposed into 8 IMFs and the Cyclic Spectral Coherence of each IMF was calculated. The result of the first two IMFs is presented in Figure 8c,d respectively. As it can be seen, the result is slightly less contaminated compared with the raw and de-noised signal. However, more interestingly, in Figure 8e,f which shows the result on the 6th and 8th IMFs two other frequencies are visible. These frequencies are the first stage GAPF frequency (5.3 Hz) and the second stage rotational frequency (3.5 Hz).
The presence of both first stage gear mesh and GAPF cyclic behaviour clearly indicates the imbalance applied to the first stage. The existence of the second stage rotational frequency could as well indicate the imbalance. Nonetheless, it can also show the normal vibration contaminating the signal.

Wind Turbine Gearbox Analysis
The studied gearbox in this part contained 11 different types of bearings.
Including the corresponding frequencies of the gears, more than 50 different frequencies of interest were available to analyse. In the present study only the result corresponding to the planetary stage of the gearbox and its bearing defect are shown.
The vibration signal used to analyse the bearing defect at the planetary stage was obtained from the Radial sensor placed at the planetary stage of the gearbox.
The damaged bearing was a cylindrical roller bearing with its characteristic frequencies listed in Table 3. Given the wind turbine working at full rated load, the planets were turning at around 1.07 Hz. The other frequency of interest in the planetary stage fault is the gear-mesh of the planets which was calculated to be 140.84 Hz. The next nearest calculated frequencies to this frequency were 99.1 Hz and 192.5 Hz belonging to the intermediate and high-speed shaft bearings. As it can be seen in the Figure 9a,b, the gear-mesh of the planets (140.8 Hz) are present in both spectra, indicating possible damage in the planetary gear. It is also evident that the cyclic spectral coherence of the first IMF of the signal is considerably less contaminated with noise originating from other gears. Similarly in the Figure 9c,d, the Ball Passing Frequency of Inner race (BPFI) at 8.8 Hz is distinctly visible. The other closest present frequencies to this frequency are 7 Hz and 9.5 Hz corresponding to the sun gear defect and carrier bearing frequencies. Again, similar to the gear mesh analysis the cyclostationary analysis on the first IMF of the signal shown in Figure 9d shows less artefacts related to the other gears.
For comparison purposes, the same signal recorded from the planetary stage was also analysed using FFT, Hilbert transform and Cepstrum that are commonly used in condition monitoring and fault detection. As it can be seen in Figure 10a,b the gear-mesh frequency is not detectable in the applied Cepstrum, FFT or Hilbert analysis while the cyclic spectral coherence enhanced with EEMD preconditioning was able to detect it. Similarly, Figure 10c,d show the same analysis applied to the signal to detect BPFI frequency. Similar to the gear mesh analysis, the BPFI frequency was not detected using either Cepstrum or FFT analysis. Although, the Hilbert transform analysis shows some peaks around the BPFI frequency, the amplitude of the peaks are not very significant which can lead to underestimation of the severity. Again, using the cyclic spectral coherence with EEMD the BPFI frequency was detected clearly reducing the likelihood of underestimation of the severity of the defect.

Conclusions
The rapid growth in the use of wind turbines combined with their operation under severe loading conditions has increased the need for efficient condition monitoring technologies. Early fault detection is important for wind farm operators in order to improve predictive maintenance strategies, reduce unexpected downtime and associated repair costs. Continuously variable operation of the wind turbines and complexity of the multi-stage gearboxes employed makes it very challenging to detect developing faults particularly during the early stages of evolution.
The cyclostationary technique has been proven to be a very effective tool to detect faults in rotating machinery. However as the complexity of the system increases, the performance of this technique degrades and the results get contaminated with various artefacts coming from different parts of the system.
In this paper the Ensemble Empirical Mode Decomposition technique was used to precondition the vibration signal prior to perform the cyclostationary technique. It has been shown that this can improve the performance of the cyclostationary technique and reduce the contamination in the resulting analysed signals.
The Ensemble Empirical Mode Decomposition technique is a very powerful method that has not been fully explored in the field of wind turbine condition monitoring. The present study has shown the potential of this technique to be employed in conjunction with other techniques to enhance the overall performance of various signal processing techniques in wind turbine gearbox condition monitoring. The ability of this technique to separate the signal into its intrinsic modes, as it was demonstrated in this paper, opens new opportunities in advanced signal processing that are yet to be explored in more depth.