Wind turbine gearbox fault diagnosis via adaptive IMFogram

Abstract

Wind turbines usually operate in mountainous, offshore and other field environments. As the wind speed is constantly changing, the gearboxes of wind turbines operate under variable operating conditions, resulting in a non-stationery and time-varying signal characteristic. For existing time–frequency analysis (TFA) methods, they undergo low time–frequency energy concentration as well as noise interference during the fault feature extraction. IMFogram, a recently proposed time–frequency representation (TFR) method, calculates local frequency and amplitude information of signals based on fast iterative filtering (FIF). Moreover, it is able to remove uncertainty on the time–frequency plane, allowing a long enough time window to reduce the local boundary effects associated with the signal decomposition methods. Therefore, it has potential advantages for TFA. However, the actual performance of IMFogram calculation largely depends on the decomposition effect provided by FIF and the selection of model parameters. With the purpose of obtaining a more distinct and focused TFR of the wind turbine vibration signal, a novel TFA method, namely Adaptive IMFogram (AIMFogram), is proposed in this paper. Firstly, a range of intrinsic mode functions (IMFs) are gained after the decomposition of vibration signals using an improved variational nonlinear chirp mode decomposition, which is treated as the input to AIMFogram to solve the unsatisfied result by traditional FIF under strong noise. Secondly, Renyi entropy-based particle swarm optimization is conducted to perform adaptive parameter optimization of the AIMFogram. Later, using the optimal parameter, the proposed AIMFogram achieves a high-resolution TFR of complex multi-component vibration signals. To verify whether the proposed AIMFogram is effective or not, the method was successfully employed to the fault diagnosis of the rolling bearing in an experimental bench and the gearbox in a 850-kW wind turbine.

Keywords

Adaptive IMFogram variational nonlinear chirp mode decomposition particle swarm optimization wind turbine fault diagnosis

Introduction

With quick depletion of global energy resources, more and more countries are turning to wind energy as their main source of renewable energy. According to statistics, global wind power generation has been growing at an average annual rate of 15.9% since 2009, and the Global Wind Energy Report 2023 states that new grid-connected capacity of global wind power is expected to reach 680 GW in the next 5 years.^1,2 However, as wind turbines often stand in remote mountainous and coastal areas, they are subjected to high winds and other impacts from the outside environment, resulting in a dramatic increase in the failure rate of turbine components under severe operating conditions. The occurrence of a fault will have a significant impact on the normal function of the equipment which is particularly critical for expensive large-scale wind turbines since a small error could also lead to substantial economic losses. Therefore, it has great engineering importance to realize fault feature recognition of wind turbine gearbox.^3,4

The operating conditions of planetary gearboxes in wind turbine drive trains usually vary due to the changing wind speed, resulting in significant non-stationary feature of vibration signals. For these signals, the spectral characteristics are time-varying. Traditional analysis methods in time domain or frequency domain alone cannot fully describe the non-stationary characteristics, while time–frequency analysis (TFA) provides information on the variation of signal spectrum with time, proving itself to be an effective tool for non-stationary signal processing.^5,6

In the course of past research, Short-time Fourier Transform (STFT) and the Continuous Wavelet Transform are frequently employed for TFA. However, due to the limitations imposed by the Heisenberg uncertainty principle, it is difficult for signals to simultaneously achieve optimal resolution in both time and frequency domains.^7,8 These disadvantages of the traditional TFA methods reduce the sensitivity of the diagnostic system to less obvious faults, for instance, early weak failures and faults under noisy environments.⁹ To improve TFA methods’ ability of detecting failures under complicated environments, a number of advanced methods have been presented over the past years. Daubechies et al.¹⁰ proposes the Synchrosqueezing Transform (SST), which utilizes a Synchrosqueezing algorithm to compress time–frequency coefficients to their respective instantaneous frequency locations. This method improves the resolution of the time–frequency plane by concentrating the time–frequency curves. Empirical Mode Decomposition¹¹ (EMD) is a signal decomposition method that is adaptive in nature. It decomposes a signal into multiple Intrinsic Mode Functions (IMFs) that represent its local characteristics of the signal based on its time scale. The IMFs containing fault components exhibit more pronounced transient features.¹² Hilbert-Huang Transform¹³ (HHT) is based on the EMD and then applies the Hilbert transform to compute the analytic signal associated with each IMF, thereby obtaining the instantaneous frequency curve in the time–frequency domain. However, EMD suffers from modal confusion, endpoint effects and difficult determination for stopping conditions.^14,15 For the sake of overcoming these drawbacks, variational mode decomposition (VMD) decomposes the signal into a variable division mode, which essentially employs multiple sets of adaptive Wiener filters to realize adaptive segmentation.¹⁶ Comparing with the EMD, VMD provides greater noise robustness and achieves weaker endpoint effects.¹⁷ However, there exists some difficulties in the decomposition when choosing the appropriate hyperparameters, including the number of decomposition patterns, the initial center frequency and the balance parameters.¹⁸ In addition, VMD is limited to processing signals in narrowband modes and does not work well for the analytical processing of modulated signals in broadband modes, especially those containing crossover modes. Hence, variational nonlinear chirp mode decomposition (VNCMD) is designed to convert broadband signals into narrowband signals using demodulation techniques.¹⁹ Moreover, it is able to extract all signal patterns simultaneously. Nevertheless, similar to VMD, it suffers from the same problem of requiring a priori knowledge to determine parameters such as the initial center frequency.^20,21

Recently, Antonio et al. proposed a TFA method called IMFogram.^22,23 The method first calculates the decomposition of the signal by fast iterative filtering (FIF), and then analyses each decomposed IMF using IMFogram. It identifies both local frequency and local amplitude information for each IMF at the same time and ultimately produces a time–frequency representation (TFR) with accurate resolution. Although IMFogram provides a clearer and more focused TFR, it still has some shortcomings. Firstly, FIF is a fast implementation of Iterative Filtering (IF) on the basis of the Fast Fourier Transform (FFT) framework and is a moving average obtained by means of a double convolution filter. However, this filter is not smooth enough and its IMFs fluctuates spuriously under the disturbance of strong oscillations. Secondly, the practical effectiveness of IMFogram for TFA of multi-component signals is mainly influenced by the choice of several parameters, such as the number of instantaneous frequencies and the length of the time sliding window.

Therefore, a novel TFA method, called adaptive IMFogram (AIMFogram), is put forward in this paper and is applied to achieve high-resolution TFR of vibration signals of wind turbines. AIMFogram is put forward based on the theoretical framework of the existing IMFogram. Firstly, VNCMD is used to replace the FIF as the input to the AIMFgram to avoid interference from other frequency band oscillations. However, since VNCMD relies on the initial center frequency setting, this paper obtains the initial center frequency of the signal by means of a popular multi-ridge detection (MRD) algorithm.^23,24 Secondly, particle swarm algorithm (PSO) is chosen as the parameter adaptive optimization method for the AIMFogram in this paper. Renyi entropy is also employed as the fitness function to gain the optimal position of the particles, and finally the optimal combination of parameters is selected for the AIMFogram, namely, the number of drawn IMFs and the sliding time window length.

The remaining sections of this paper are organized as follows: in Section “Theoretical descriptions,” the basic theory of IMFogram and VNCMD is introduced, followed by a detailed theoretical description of AIMFogram, and later the process steps of the wind turbine fault diagnosis are described. In Section “Analogue signal simulation analysis,” analogue signal simulation analysis with time-varying features will be employed to verify the superiority of the proposed AIMFogram. In Section “Experimental data analysis,” the proposed method will be applied for vibration signal analysis on a test rig and an 850-kW wind turbine gearbox with rolling bearing fault. These measurement data will further demonstrate the effectiveness of the proposed AIMFogram. Finally, Section “Conclusion” concludes the paper.

Theoretical descriptions

IMFogram

IMFogram is an analogue of a spectrogram which allows fast calculations based on the decomposition results and identification for both local frequency and amplitude information of the signal. Antonio decomposes the signal in the literature²² using FIF, which is an IF based on a fast implementation of the FFT. The basic idea is to decompose a non-linear and non-smooth signal into IMFs by means of a double convolution filter. Next, a stopping criterion is introduced to ensure the termination of the algorithm. Finally, each IMF is analyzed from a time–frequency perspective to obtain a TFR of the IMFogram.

First, the IMF is defined as:

IMF = {(IMF (t_{j}))}_{0}^{n - 1}

(1)

where $t_{j}$ is the time point, $t_{j} \in [t_{0}, t_{n - 1}]$ , Then, the instantaneous amplitude (IA) and instantaneous frequency (IF) of the IMF at $t_{j}$ are $i A^{IMF} (t_{j})$ and $i F^{IMF} (t_{j})$ respectively.

Arrange $g$ as the linear interpolation of the local maxima of $| IMF |$ , and the IA function of $IMF$ is expressed as:

i A^{IMF} (t_{j}) = max {g (t_{j}), | IMF (t_{j}) |}

(2)

where $i A^{IMF} : [t_{0}, t_{n - 1}] \to R$ ; $(t_{j}, i A^{IMF} (t_{j}))$ is the linear interpolation.

Let the zero crossings of $f$ be ${z_{1} < z_{2} < \dots < z_{p}}$ , and define $y_{j} = \frac{1}{z_{j + 1} - z_{j}} (j = 1, 2, \dots, p)$ , $z_{0} = t_{0}$ , $z_{p + 1} = t_{n - 1}$ , $y_{0} = y_{1}$ , $y_{p + 1} = y_{p}$ . Set $f$ as the linear interpolation of $IMF (t_{j})$ on the interval $[t_{0}, t_{n - 1}]$ , then the instantaneous frequency function $i F^{IMF} : [t_{0}, t_{n - 1}] \to R$ is the linear interpolation of ${(z_{j}, 2 y_{j})}_{j = 0}^{p + 1}$ .

The IA and IF of the final IMF at $t_{j}$ are denoted as $i A^{IMF} (t_{j})$ and $i F^{IMF} (t_{j})$ respectively. Next, the local frequencies and amplitudes are obtained from the instantaneous frequencies and amplitudes.

Define the average matrix $T$ of $n \times n$ as:

T_{i, k} = {\begin{matrix} 1 / I, & 1 + (p - 1) I \leq i, k \leq pI, for p = 1, 2, \dots, n / I \\ 0, & otherwise \end{matrix}

(3)

where $I$ is a time window of fixed size and $n$ is usually a multiple of $I$ . $T_{i, k}$ is the $(i, k)$ -entry of $T$ .

In the next step, set $IMF \in R^{n}$ and make $i A^{IMF} \in R^{n}$ as IA and $i F^{IMF} \in R^{n}$ as IF respectively. Their average values over a time window of length $I$ are the local amplitude $L A^{IMF}$ and local frequency $L F^{IMF}$ respectively, which is denoted as:

L A^{IMF} = T (i A^{IMF}), L F^{IMF} = T (i F^{IMF})

(4)

Rewrite the parameters based on the above theoretical derivation, the IMFogram matrix $A$ can be summarized as the following equation:

A ({LF}_{I_{j}}^{(k)}, j) = \bar{A} ({LF}_{I_{j}}^{(k)}, j) + {LA}_{I_{j}}^{(k)}, \forall j = 1, 2, \dots, N

(5)

where $\bar{A}$ is the previously calculated IMFogram matrix with all zero entries as originally set, $I_{j}$ is a sliding time window, $N$ is the number of IMFs, and ${LA}_{I_{j}}^{(k)}$ represents the local amplitude of the k-th IMF on $I_{j}$ . Then, calculate the weighted average of the envelopes of the interpolated $| IM F_{k} |$ local maxima over the sliding time window $I_{j}$ and use the amplitude of $IM F_{k}$ itself over a sliding time window $I_{j}$ as a weight. Furthermore, ${LF}_{I_{j}}^{(k)}$ denotes the local frequency of the k-th IMF on $I_{j}$ and a weighted average of the interpolated envelope of the point set ${(z_{j}^{(k)}, 2 y_{j}^{(k)})}_{j}$ over the time interval $I_{j}$ .

The IMFogram produces TFR results from time domain, local period and local amplitude information, and makes its TFR results more localized by adding windows to the signal. If watching closely, one can see that the TFR results look like they are made up of many rectangular windows. Thus, it is possible to represent the decomposed results accurately on the time–frequency plane.

Proposed AIMFogram algorithm

The actual effectiveness of the original IMFogram depends on the results of the FIF decomposition, but two problems still exist in the vibration signal decomposition. Firstly, to make the result of the decomposition more reflective of the local signal characteristics, the filter function must be of finite length. However, the waveform of a filter function with finite length is not smooth enough in general. Secondly, under the influence of noise, there will still be a noisy component in the multiple IMF components of the FIF decomposition, which will directly lead to a time–frequency curve of noise in the time–frequency diagram obtained by the IMFogram. In addition, the IMFogram-calculated results of the TFR depend on the correct selection of parameters. As a result, this paper will improve on the IMFogram and propose an adaptive IMFogram (AIMFogram). At first, we replace FIF with VNCMD, which means using a demodulation technique to convert a broadband FM signal into a narrowband signal and then decomposing the signal into multiple nonlinear chirp mode components by means of a variational model that is more suitable for the vibration signal analysis under speed changing conditions. However, VNCMD requires the setting of an initial center frequency (ICF), which demands a priori knowledge. Therefore, this paper proposes to use MRD to extract the coarser time–frequency distribution ridges of the STFT as the ICF.^24,25 Then, by introducing a PSO algorithm and combining it with the Renyi entropy model, the AIMFogram is made to achieve adaptive optimization of the parameters.

Improved VNCMD

In practical scenarios, signals often consist of multiple nonlinear frequency-modulated modes accompanied by Gaussian white noise. The mathematical model $x (t)$ can be expressed as follows:

x (t) = \sum_{k = 1}^{K} x_{k} (t) = \sum_{k = 1}^{K} a_{k} (t) \cos (2 π \int_{0}^{t} f_{k} (s) d s + ϕ_{k}) + n (t)

(6)

where $K$ represents the number of frequency-modulated components in the signal; $a_{k} (t)$ and $f_{k} (t)$ respectively represent the IA and IF of the K-th mode; $ϕ_{k}$ represents the initial phase; $n (t)$ represents the noise. Then, Equation (7) can be rewritten in demodulated form below:

x (t) = \sum_{k = 1}^{K} u_{k} (t) \cos (2 π \int_{0}^{t} {\bar{f}}_{k} (s) d s) + v_{k} (t) \sin (2 π \int_{0}^{t} {\bar{f}}_{k} (s) d s) + n (t)

(7)

where $u_{k} (t)$ and $v_{k} (t)$ denote the two demodulated signals; ${\bar{f}}_{k} (s)$ is the frequency function of the individual demodulation operators.

Since the data analyzed are usually sampled in discrete form, when the sampling interval is $t = t_{0}, \dots, t_{N - 1}$ , the resulting constrained problem is expressed as follows:

{\begin{matrix} min_{| u_{k} | \cdot | v_{k} | \cdot | f_{k} |} = {Γ (w) + \sum_{k} (‖ Ω u_{k} ‖_{2}^{2} + ‖ Ω v_{k} ‖_{2}^{2})} \\ s . t . w = x - \sum_{k} (C_{k} u_{k} + D_{k} v_{k}) \end{matrix}

(8)

where $Ω$ is a second-order difference operator used to eliminate boundary effects. $Γ$ denotes the penalty function. $w$ is the noise variable. $C$ and $D$ are two diagonal matrices.

To address the issue in Equation (8), the Lagrangian multiplier $λ'$ is introduced, and the corresponding extended Lagrangian is expressed as follows:

\begin{matrix} L_{a} ({u_{k}} {v_{k}} {f_{k}}, w, λ') = Γ (w) + \sum_{k} (‖ Ω u_{k} ‖_{2}^{2} + ‖ Ω v_{k} ‖_{2}^{2}) \\ + \frac{α}{2} ‖ w + \sum_{k} (C_{k} u_{k} + D_{k} v_{k}) ‖ \\ - \frac{1}{2 α} ‖ λ' ‖_{2}^{2} \end{matrix}

(9)

where $α$ is a penalty parameter. It is worth noting that a change in $α$ changes the bandwidth of the time–frequency filter. In this paper, to avoid excessive fluctuations in the instantaneous frequency curve, a smaller bandwidth parameter is therefore set, with a range of values between $1 \times 10^{- 11} and 1 \times 10^{- 14}$ .

Further, alternating direction method of multipliers algorithm is used for instantaneous frequency optimization as well as modal component solution. The specific steps are followed below:

Step 1, update both demodulation signals:

\begin{matrix} u_{k}^{a + 1} = {(\frac{2}{α} Ω^{T} Ω + C_{k}^{T} C_{k})}^{- 1} C_{k}^{T} \\ (x - \sum_{m \neq k} C_{m} u_{m} - \sum_{m} D_{m} v_{m} - w - \frac{1}{α} λ') \end{matrix}

(10)

\begin{matrix} v_{k}^{a + 1} = {(\frac{2}{α} Ω^{T} Ω + D_{k}^{T} D_{k})}^{- 1} \\ D_{k}^{T} (x - \sum_{m \neq k} C_{m} u_{m} - \sum_{m} D_{m} v_{m} - w - \frac{1}{α} λ') \end{matrix}

(11)

Step 2, the instantaneous frequency $f_{k}$ is updated according to the results obtained in Equations (10) and (11):

Δ f_{k}^{a + 1} = {(\frac{2}{μ} Ω^{T} Ω + I)}^{- 1} (- \frac{1}{2 π} \frac{d}{d t} (\arctan (\frac{v_{i}^{k + 1} (t)}{u_{i}^{k + 1} (t)})))

(12)

f_{k}^{a + 1} = f_{k}^{a} + η Δ f_{k}^{a + 1}

(13)

where $μ$ is a parameter used to control the increment of the instantaneous frequency. $I$ is the unit matrix. $η \in (0, 1)$ is the scaling factor used to stabilize the algorithm.

Step 3, update and reconfigure the signal $x$ :

x_{k}^{a + 1} = C_{k} u_{k}^{a + 1} + D_{k} v_{k}^{a + 1}

(14)

Due to the limitations of the algorithm, VNCMD requires an initial IF that is artificially preset and the estimated frequency band of the pattern depends on the initial IF. For the sake of solving the problem of initial IF for each mode component, a popular multi-ridge detection algorithm of MRD^24,25 is proposed to perform initial IF extraction on the time–frequency results of STFT. By predefining the number of modes K, the calculation formula for this algorithm is as follows:

\begin{matrix} E (ϕ) = \sum_{k = 1}^{K} - \int_{- \infty}^{+ \infty} {| TFR (t, ϕ_{k} (t)) |}^{2} d t \\ + \int_{- \infty}^{+ \infty} (λ \cdot ϕ_{k}^{'} {(t)}^{2} + β \cdot ϕ_{k} {(t)}^{2}) d t \end{matrix}

(15)

where $\sum_{k = 1}^{K} (t, ϕ_{k} (t))$ is an estimate of the IF curve in the time–frequency domain, $λ$ and $β$ are two parameters to adjust the level of regularization.

Adaptive optimization of parameters for IMFogram

If IMFogram is to gain a TFR of the energy concentration of vibration signal, it depends primarily on the number $N$ of IMFs plotted and the length of the sliding time window $I_{j}$ . As $N$ gets larger, the IMFogram obtains a time–frequency diagram that contains more local IMFs and exhibits better local detail but less coherence. And the larger the $I_{j}$ , the better the coherence of the time–frequency ridges, but the poorer the local detail. Different parameters have an impact on the final time–frequency resolution, which in turn increases the estimation error of instantaneous frequency and prevents the time–frequency picture from showing the internal relationship between signal frequency and time. To resolve the issue, the PSO algorithm is chosen as the parameter optimization method for IMFogram in this paper. This algorithm is a population intelligence optimization algorithm with good global optimization capability for parameters, which avoids the situation of individual optimization falling into a local optimum solution.^26,27

Each particle in the PSO has two properties: velocity and position. The two properties of particle $i$ at the d-th at $k + 1$ iterations are denoted by $v_{id}^{k + 1}$ and $x_{id}^{k + 1}$ respectively. When using the PSO, the issue to be solved is modelled as the fitness function. The historical best position of each particle in the particle swarm is denoted by $p_{i}$ . At the same time, the whole population has a global optimal position $p_{g}$ . The formulas for updating particle positions and velocities are as follows:

v_{id}^{k + 1} = ω v_{id}^{k} + c_{1} r_{1} (p_{id}^{k} - x_{id}^{k}) + c_{2} r_{2} (p_{gd}^{k} - x_{id}^{k})

(16)

x_{id}^{k + 1} = x_{id}^{k} + v_{id}^{k + 1}

(17)

In the Equation, k represents the iteration number; i represents the particle’s index; d represents the dimensionality of the particle; $c_{1}$ , $c_{2}$ are learning factors; $r_{1}$ , $r_{2}$ are uniformly distributed pseudo-random numbers.

With the purpose of finding the best position for the particle, an adaptation function is required to be given. Therefore, calculate the fitness function value one time after every iteration and update the position and velocity of particle by comparing the function values. The inertia weights of the PSO only guarantee the general direction of the swarm search but cannot precisely control each iteration of the algorithm. The Renyi entropy, as an extension of the information entropy, reflects the features of the signal in terms of its time–frequency distribution.²⁸ Compared to information entropy, Renyi entropy is more sensitive to the subtle signal changes, making it easier to identify small changes in the signal. Therefore, this paper optimizes the PSO algorithm based on the Renyi entropy model, and obtains characteristics information of the particle swarm in the search process. The inertia weights of each iteration are dynamically adjusted based on the entropy difference, so as to gain the optimal parameter combination for the IMFogram algorithm, that is, the number $N$ of drawn IMFs and the sliding time window length $I_{j}$ .

Renyi entropy is calculated as

H_{α} = \frac{1}{1 - α} \log (\sum_{i = 1}^{M} p {(x_{ki})}^{α})

(18)

where $α \geq 0$ and $α \neq 1$ . $H_{α}$ has a maximum value when $α = 0$ , at which Renyi entropy has a high sensitivity to endpoint changes. The Renyi entropy is equivalent to the information entropy when $α \to 1$ and is less sensitive to changes in endpoints. $M$ is the total number of particle populations. $p (x_{ki})$ is the weight of the i-th particle in the k-th iteration of the particle population.

The entropy value is not always in a straight line down, there can be an increase. The entropy difference represents the absolute magnitude of the difference between the minimum entropy value and the entropy value after the latest iteration of optimization search. The formula is as follows:

d (t) = | H (t_{min}) - H (t) |

(19)

The specific steps of the parametric adaptive optimization IMFogram algorithm are described below:

(1) Initialize the optimization parameters $N$ and $I_{j}$ in the IMFogram, generate a 2 × m two-dimensional space vector in random as the starting position and initialize the particle velocity randomly.

(2) Initialize the PSO parameters and set the Renyi entropy to the fitness function.

(3) Aiming at every particle, compare its current fitness with the fitness at its best position encountered so far. If the fitness improves, update $p_{id}$ .

(4) Compare the best fitness of each particle with the fitness at the best position encountered by the swarm. If the former is better, update the best minimum of individual and population history $p_{gd}$ .

(5) Adjust the speed and position of the particle and repeat steps (1)–(4). If the maximum number of iterations is reached or the entropy difference becomes zero, the particle search ends and optimum position of the particle is output.

(6) The result of experiencing the best positioned particle is used as the global optimal weight, that is, the value after decoding of the best position encountered by the swarm, which is used as the structural parameter of the IMFogram.

Procedure of the proposed AIMFogram

Under variable speed conditions, the obtained fault characteristic frequencies after time–frequency transformation are time-varying, and traditional fault diagnosis methods have difficulty in identifying their changing instantaneous frequency. Therefore, based on the existing IMFogram theoretical framework, this paper proposes an enhanced variational nonlinear chirp mode decomposition, which achieves the acquisition of the IMF after signal decomposition without prior knowledge. A PSO algorithm based on Renyi entropy is then used to gain the best parameter combination for IMFogram, thus achieving the optimal TFR of IMFogram. The proposed fault diagnosis method follows the procedure below:

(1) First, use an improved VNCMD for the mode decomposition of the wind turbine vibration signal to obtain a number of component IMFs.

(2) Then, using the proposed AIMFogram method, an optimal TFR of the multi-component signal of wind turbine vibration is achieved by determining the parameter combination through PSO and Renyi entropy.

(3) Finally, by calculating the fault characteristic frequencies and extracting the relevant instantaneous frequency curves from the time–frequency plots, the identification of rolling bearing fault types in variable speed operating conditions of the wind turbine is achieved. The flow chart of the proposed AIMFogram applied to fault diagnosis is shown as the Figure 1 below.

Figure 1.

General diagram for fault diagnosis using AIMFogram.

Analogue signal simulation analysis

To validate the advantages of the AIMFogram algorithm, a set of multivariate simulated signals with frequency variations over time are utilized for analysis. The simulated signal is configured as follows:

\begin{matrix} s_{1} (t) = 0.2 \sin (2 π (10 t^{2} + 25 t) \\ s_{2} (t) = 0.2 \sin (2 π (70 t - 0.2 \sin (10 π t))) \\ s (t) = s_{1} (t) + s_{2} (t) + n (t) \end{matrix}

(20)

where $s_{1} (t)$ and $s_{2} (t)$ denote two different frequency-modulated components; $n (t)$ is the noise component. The instantaneous frequencies of the two frequency-modulated components are $i F_{1} = 20 t + 25$ and $i F_{2} = 70 - 2 π \cos (10 π t)$ , respectively. To simulate the noise environment of rolling bearings, Gaussian noise with a signal-to-noise ratio (SNR) of $5 dB$ is applied. Additionally, the analogue signal has a length of $L = 2000$ and a sampling frequency of 2000 Hz.

The time domain waveform of the analogue signal is drawn in Figure 2(a) and the corresponding frequency spectrum is plotted in Figure 2(b). From the spectrum, the signal components at 25–45 and $70 Hz$ can be detected. Unfortunately, the time-varying feature still cannot be inspected. The ideal TFR of an analogue signal is shown in Figure 3. Generally, if processed results of the studied TFA method are the closest to the ideal IF curve, then the method has the best TFR.

Figure 2.

Time domain diagram and frequency spectrum of analogue signal. (a) Time domain diagram of analogue signal and (b) frequency spectrum of analogue signal.

Figure 3.

Ideal IF curve for analogue signal.

For the signal as STFT, the window function is chosen as Gaussian window function with a window width of 300, and the computations displayed in Figure 4(a). Obviously, the STFT results are blurred at time–frequency plane. The signal calculation results using SST are expressed as Figure 4(c), where to some extent energy in the time–frequency plane is concentrated, but the IF curve still suffers from energy divergence as can be seen from the local plot in Figure 4(d). Subsequently, Figure 4(e) and (f) demonstrate the TFR and local plots of the IMFogram by FIF after its decomposition. It can be found that energy dispersion in the time–frequency plane has been solved, but due to signal modulation and the influence of noise, the result obtained after FIF decomposition has a large deviation from the ideal IF curve, and there still presents noise. Figure 4(g) and (h) plot the TFR of the above signals using AIMFogram. The initial population size of the particle swarm is 40, the inertia weight $w$ ranges from 0.4 to 0.9, the maximum number of iterations is 20, and the final $N$ and $I_{j}$ attain the best parameter combination 410 and 80. Figure 5 displays the variation of entropy difference with the times of iterations. After the AIMFogram has been decomposed by improved VNCMD to resolve the signal modulation phenomenon and obtain better decomposition results, the best TFR is obtained using the PSO-optimized IMFogram. A comparison of the above methods shows that the AIMFogram produces a more distinct and focused time–frequency energy ridge in the time–frequency diagram.

Figure 4.

TFR of the simulated signal by different methods. (a) Result provided by STFT, (b) partial view of STFT result, (c) result provided by SST, (d) partial view of SST result, (e) result provided by IMFogram, (f) partial view of IMFogram result, (g) result provided by AIMFogram, and (h) partial view of AIMFogram result.

Figure 5.

Variation of entropy difference with the number of iterations.

To demonstrate the superiority of the VNCMD method over the FIF algorithm in signal decomposition, this study conducts an analysis using two performance metrics: correlation coefficient and mean square error (MSE). The correlation coefficient is a metric used to assess the correlation between the decomposed synthesized signal and the original noise-free signal. The value of correlation coefficient closer to 1 always indicates a stronger correlation between the two signals. Commonly, the smaller the MSE value is, the stronger the multi-component signal decomposition ability of the employed algorithm is. Therefore, we investigate the performance of function indicators under different noise environments, as shown in Figure 6. In Figure 6(a), VNCMD clearly outperforms FIF. In Figure 6(b), for relatively high SNR values, VNCMD and FIF display similar MSE values, but as the SNR drops below 10 dB, the MSE of FIF increases dramatically. Therefore, VNCMD exhibits superior decomposition performance compared to FIF.

Figure 6.

The variation of the performance metrics with SNR input. (a) Correlation coefficient, and (b) MSE.

To assess the performance of the aforementioned methods in terms of TFR, this study quantitatively analyzes the concentration of energy in the IF curves of different TFA methods using Renyi entropy values.²⁹ Commonly, the lower Renyi entropy value indicates a more concentrated TFR. Figure 7 shows the Renyi entropy values for four different methods in different noise scenarios. As can be seen from the Figure 7, the researched AIMFogram has the lowest Renyi entropy, showing that it comprehensively outperformed the other methods in TFA.

Figure 7.

Renyi entropy for different methods with different noise.

Experimental data analysis

Application to rolling bearings with outer ring failures on test bench

To further verify whether the proposed method in practical applications is effective or not, experimental analysis was carried out using test bench data. A physical diagram and a sketch of the structure of the test rig is shown in Figure 8. During the experiment, the motor gradually accelerated from a stationary state to 1920 rpm in 10 s. The vibration acceleration sensor was radially placed on the external bearing housing for measurement. We collected the acceleration process and selected the final 1.6 s of the acceleration period for analysis, with a sampling frequency of 2560 Hz. Table 1 lists the parameters of the test bench and the model of the faulty rolling bearing in detail. The vibration signal and its frequency spectrum acquired from the test bench are shown in Figure 8.

Figure 8.

Physical drawing and structural sketch of the experimental bench. 1- Servo motors, 2- Couplings, 3- Bearing seat, 4- Magnetic powder brake, 5- PCB acceleration sensors, 6- External bearing seat, 7- Gearboxes.

Table 1.

The parameters of the test rig gearbox and rolling bearing.

The gear ratio of the gearbox	Type of rolling bearing	Number of rolling elements of rolling bearings	Type of rolling bearing failure
2/3	6020	$z = 8$	Outer race failure

The feature frequency of the bearing failure could not be clearly determined in the spectrum of Figure 9(b); thus the results gained by STFT of the acquired signal are demonstrated in Figure 10(a) and Fig. 10 (b) shows its local enlargement. From this it is clear that due to the inevitable noise and other interference factors in actual measurement process, the time–frequency energy ridge clearly diverges in its local magnification, and its instantaneous frequency cannot be accurately obtained. The analysis results provided by SST are presented in Figure 10(c), with its local enlargement in Figure 10(d), from which the components in the signal can be seen apparently. However, the IF trajectory still exhibits a slight energy divergence. The IMFogram analysis results are shown in Figure 10(e), from which it can be found that the interference of rotating frequency and noise lead to deviations in the components obtained after FIF decomposition, such as triplets of the fault characteristic frequency. There is also information of noise, which causes interference in the next fault diagnosis. Figure 10(f) shows its local magnification; it is obvious that the time–frequency ridge energy is more concentrated compared to those of STFT and SST, but other components exist.

Figure 9.

Vibration signal and its frequency spectrum from the test bench. (a) Vibration signal from the test bench and (b) frequency spectrum of the vibration signal.

Figure 10.

TFR of the vibration signal for four different methods. (a) Result provided by STFT, (b) partial map of STFT result, (c) result provided by SST, (d) partial map of SST result, (e) result provided by IMFogram, (f) partial map of IMFogram result, (g) result provided by AIMFogram, and (h) partial map of AIMFogram result.

The TFA results generated by AIMFogram are presented in Figure 10(g) and (h). The initial population size of the particle swarm is 60, the inertia weight $w$ is taken to range between 0.4 and 0.9, the maximum number of iterations is 20, and the final $N$ and $I_{j}$ obtain optimal parameter combinations of 480 and 260. Figure 11 demonstrates the result of the variation of entropy difference with the number of iterations. The result illustrates that the AIMFogram is better at capturing the fault feature components in the presence of noise interference or other components. It outperforms other methods in both concentration of IF energy ridge lines and clarity of the time–frequency plane, thus achieving a time–frequency representation closer to the ideal.

Figure 11.

Variation of entropy difference with the number of iterations.

Theoretically, the fault characteristic frequencies of the rolling bearing, namely, the ball passing frequency of outer race $f_{o}$ , inner race $f_{i}$ , cage $f_{c}$ and ball $f_{b}$ ^30–33 can be defined as:

{\begin{matrix} f_{o} = \frac{1}{2} z (1 - \frac{d}{D} \cos α) f_{r} \\ f_{i} = \frac{1}{2} z (1 + \frac{d}{D} \cos α) f_{r} \\ f_{c} = \frac{1}{2} (1 - \frac{d}{D} \cos α) f_{r} \\ f_{b} = \frac{D}{2 d} {1 - {(\frac{d}{D} \cos α)}^{2}} f_{r} \end{matrix}

(21)

where $z$ is the number of the balls, $d$ , $D$ and $α$ denote the ball diameter, the pitch diameter and the contact angle, respectively, and $f_{r}$ is the rotating frequency of the shaft.

According to the rolling bearing of test bench type 6206 and the fluctuation range of among 17.9–21.3 Hz, we used the rolling bearing failure characteristic frequency formulae (21) to calculate the failure characteristic frequency of the inner race $f_{i}$ , outer race $f_{o}$ , cage $f_{c}$ , ball $f_{b}$ (as shown in Table 2). By combining the TFA results in Figure 10(g) with the characteristic frequency ranges in the Table 2, it could be conclusively determined that the bearing is experiencing an outer race fault.

Table 2.

Failure characteristic frequency of rolling bearing in test bench.

Number	Failure characteristics	Characteristic frequency ranges (Hz)
1	Inner race $f_{i}$	86–102.4
2	Outer race $f_{o}$	57.3–68.3
3	Cage $f_{c}$	6.8–8.5
4	Ball $f_{b}$	33–39.3

Application to wind turbine gearbox fault identification

The actual measured vibration signal is sourced from an 850-kW wind turbine unit located at a wind farm operated by Goldwind Sci & Tech Co., Ltd. (Urumqi, China). It was collected from the sensor measuring point of the high-speed shaft front bearing of the gearbox during the operation of a specific wind turbine. This dataset is utilized to further investigate the effectiveness of the proposed method in practical operation of mechanical equipment. The model of the wind turbine gearbox is FLENDER. The structure diagram and the measurement point are shown in Figure 12. The measurement point is the front bearing of the high-speed gear shaft that has 25 rolling elements in total. The vibration signal is sampled at a frequency of 2560 Hz. The generator speed is measured by pulse voltage method.³⁴ The collected vibration signal time domain diagram and frequency spectrum are presented in Figure 13.

Figure 12.

Structure and measurement point of wind turbine. (a) Wind turbine structure diagram and (b) measurement point.

Figure 13.

Vibration signal and its frequency spectrum from the wind turbine. (a) Vibration signals from the wind turbine and (b) frequency spectrum of vibration signal.

As the gearbox is operated at variable speed, it is not possible to identify the fault characteristics in the frequency spectrum from Figure 13(b), therefore a TFA is required to obtain more information about the fault under variable speed. In the search for fault characteristic frequency components, the interest frequency band range is set to $0 - 800 Hz$ , where all the fault feature frequencies and their harmonics dominate the spectrum. The fault diagnosis scheme for the rolling bearing of the wind turbine is the same as described in Section “Application to rolling bearings with outer ring failures on test bench.” By utilizing the known parameters of the rolling bearing, measuring the rotational speed of the high-speed gear shaft of the wind turbine, and using the fault characteristic frequency formula (21), the characteristic frequency corresponding to inner ring fault of rolling bearing can be determined.

Four different algorithms were employed to analyze and calculate the wind turbine gearbox vibration signal and the results were demonstrated in Figure 14. By comparison, it is clear that the TFR of STFT and SST results have different degrees of energy divergence. Due to the strong oscillations encountered in the FIF decomposition for the original IMFogram, the time–frequency characteristic curve experienced periodic spurious fluctuations as presented in Figure 14(e) and (f). Instead, the proposed AIMFogram avoids this situation, as shown in Figure 14(g) and (h). The initial population size of the particle swarm is 50, the inertia weight $w$ is taken to range between 0.4 and 0.9, the maximum number of iterations is 20, and the final $N$ and $I_{j}$ obtain optimal parameter combinations of 560 and 450. Figure 15 displays the results of the variation of entropy difference with the number of iterations. Obviously, the TFR provided by AIMFogram offers the clearest visualization, with highly concentrated energy ridges in the TF plane, hence the fault characteristic frequencies $f_{i}$ , the two-fold fault characteristic frequency $2 f_{i}$ and the frequencies modulated by the rotating frequency can be captured effectively.

Figure 14.

TFR of the wind turbine gearbox vibration signal for four different methods.

Figure 15.

Variation of entropy difference with the number of iterations.

By querying the measurement records of the generator’s rotational speed, which was approximately n = 1250 rpm, the rotating frequency was calculated to be $f_{r} = n / 60 \approx 20.8 Hz$ . Additionally, it was known that the rolling bearing had 25 rolling elements at the measurement point. Similar to Section “Application to rolling bearings with outer ring failures on test bench,” we calculated $f_{i}$ , $f_{o}$ , $f_{c}$ , $f_{b}$ using the fault characteristic frequency formulae, and the results were presented in Table 3. By examining the TFA results of the wind turbine signal (as shown in Figure 14(g)), we could initially observe the fault frequency $f_{i}$ and the two-fold fault frequency $2 f_{i}$ . Additionally, we could ensure that the sidebands correspond to the modulation of the inner race fault frequency by the rotational frequency, denoted as $f_{c} + f_{r}$ , $f_{c} + 2 f_{r}$ , $f_{c} - f_{r}$ , $2 f_{c} + f_{r}$ , $2 f_{c} + 2 f_{r}$ , $2 f_{c} - f_{r}$ . According to the basic theory of rolling bearing fault diagnosis,³⁵ if the above fault characteristics frequencies and side frequency band appear in the TFA results, we can judge that the inner ring fault characteristics have occurred. Subsequently, the dismantling inspection of the gearbox reveals an inner ring fault in the rolling bearing, which is exactly consistent with the diagnostic result presented in this paper, and that further confirms the effectiveness of the proposed method in the application of wind turbine gearbox fault diagnosis.

Table 3.

Failure characteristic frequency of bearing in wind turbine gearbox.

Number	Failure characteristics	Characteristic frequencies (Hz)
1	Inner race $f_{i}$	310
2	Outer race $f_{o}$	208
3	Cage $f_{c}$	8
4	Ball $f_{b}$	93

Conclusion

In terms of the wind turbine gearbox under changing operating environment, this paper proposed the AIMFogram to achieve an accurate TFR of multi-components vibration signals aiming at the blurred TFR and hard extraction of signal fault characteristics. The main points of this paper are summarized as following: (1) Original FIF is replaced by an improved VNCMD, and a multi-ridge line detection method of MRD provides the ICF for the VNCMD, enabling adaptive decomposition of the signal without prior knowledge, which is more conducive to the vibration signal analysis under varying velocity. (2) PSO algorithm is employed to optimize the parameters of the structural element scale of the existing IMFogram, and the Renyi entropy is set as the fitness function so as to obtain the best scale parameters. (3) Furthermore, the proposed AIMFogram in this paper has successfully achieved fault diagnosis of a rolling bearing outer ring in a test bench and fault characterization of a gearbox in a 850-kW wind turbine. All the results above fully demonstrate the feasibility of proposed AIMFogram in this paper for rotating machinery fault diagnosis under varying speeds.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work is supported by Hubei Province Key Research and Development Plan (2021BAA194) and Guangxi Key Research and Development Plan (Guike AB21075009). Numerical calculation is supported by High-Performance Computing Center of Wuhan University of Science and Technology.

ORCID iD

Cancan Yi

References

Global Wind Energy Council (GWEC). GWEC Global Wind Report 2019. Bonn, Germany: Global Wind Energy Council, 2017.

Sadorsky

. Wind energy for sustainable development: driving factors and future outlook. J Cleaner Prod 2021; 289: 125779.

Bhardwaj

Teixeira

Soares

. Reliability prediction of an offshore wind turbine gearbox. Renewable Energy 2019; 141: 693–706.

Liu

Zhang

. A review of failure modes, condition monitoring and fault diagnosis methods for large-scale wind turbine bearings. Measurement 2020; 149: 107002.

Chen

, et al. Adaptive multiple second-order synchrosqueezing wavelet transform and its application in wind turbine gearbox fault diagnosis. Meas Sci Technol 2021; 33(1): 015110.

Feng

. A novel adaptive bandwidth selection method for Vold–Kalman filtering and its application in wind turbine planetary gearbox diagnostics. Struct Health Monit 2023; 22(2): 1027–1048.

Gröchenig

. Uncertainty principles for time-frequency representations. In: Feichtinger

Strohmer

(eds) Advances in Gabor Analysis Boston: Birkhäuser, 2003, pp.11–30.

. A concentrated time–frequency analysis tool for bearing fault diagnosis. IEEE Trans Instrum Meas 2019; 69(2): 371–381.

Zagubień

Wolniewicz

. The impact of supporting tower on wind turbine noise emission. Appl Acoust 2019; 155: 260–270.

10.

Daubechies

. Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool. Appl Comput Harmonic Anal 2011; 30(2): 243–261.

11.

Lei

Lin

, et al. A review on empirical mode decomposition in fault diagnosis of rotating machinery. Mech Syst Signal Process 2013; 35(1–2): 108–126.

12.

Wang

Shao

. Fault feature extraction of rotating machinery using a reweighted complete ensemble empirical mode decomposition with adaptive noise and demodulation analysis. Mech Syst Signal Process 2020; 138: 106545.

13.

Osman

Wang

. A new Hilbert-Huang transform technique for fault detection in rolling element bearings. In: Moamar

Sayed-Mouchaweh

(eds) Predictive maintenance in dynamic systems: advanced methods, decision support tools and real-world applications. Switzerand, AG: Switzerand AG. 2019, pp. 207–230.

14.

Zheng

Cheng

. Partly ensemble empirical mode decomposition: an improved noise-assisted method for eliminating mode mixing. Signal Process 2014; 96: 362–374.

15.

Zhang

Ren

, et al. Rolling bearing fault diagnosis utilizing variational mode decomposition based fractal dimension estimation method. Measurement 2021; 181(1–2): 109614.

16.

Dragomiretskiy

Zosso

. Variational mode decomposition. IEEE Trans Signal Process 2013; 62(3): 531–544.

17.

Rehman

Aftab

. Multivariate variational mode decomposition. IEEE Trans Signal Process 2019; 67(23): 6039–6052.

18.

Yao

Xiang

Qian

, et al. Noise source identification of diesel engine based on variational mode decomposition and robust independent component analysis. Appl Acoust 2017; 116: 184–194.

19.

Chen

Dong

Peng

, et al. Nonlinear chirp mode decomposition: a variational method. IEEE Trans Signal Process 2017; 65(22): 6024–6037.

20.

. Nonlinear chirp component decomposition: a method based on elastic network regression. IEEE Trans Instrum Meas 2021; 70: 1–13.

21.

Wei

Wang

Peng

, et al. Variational nonlinear component decomposition for fault diagnosis of planetary gearboxes under variable speed conditions. Mech Syst Signal Process 2022; 162: 108016.

22.

Cicone

Zhou

. New theoretical insights in the decomposition and time-frequency representation of nonstationary signals: the IMFogram algorithm. arXiv preprint arXiv:2205.15702, 2022.

23.

Barbe

Cicone

, et al. Time-frequency representation of nonstationary signals: the IMFogram. arXiv preprint arXiv:2011.14209, 2020.

24.

Carmona

Hwang

Torrésani

. Multiridge detection and time-frequency reconstruction. IEEE Trans Signal Process 1999; 47(2): 480–492.

25.

. Synchroextracting transform. IEEE Trans Ind Electron 2017; 64(10): 8042–8054.

26.

Medoued

Lebaroud

Laifa

, et al. Classification of induction machine faults using time frequency representation and particle swarm optimization. J Electr Eng Technol 2014; 9(1): 170–177.

27.

Dang

. A fault diagnosis scheme for rolling bearing based on particle swarm optimization in variational mode decomposition. Shock Vibr 2016; 2016, pp. 1–10.

28.

Çomak

. A modified particle swarm optimization algorithm using Renyi entropy-based clustering. Neural Comput Appl 2016; 27: 1381–1390.

29.

Popescu

Aiordachioaie

Culea-Florescu

. Basic tools for vibration analysis with applications to predictive maintenance of rotating machines: an overview. Int J Adv Manuf Technol 2022; 118: 2883–2899.

30.

Randall

Antoni

Chobsaard

. The relationship between spectral correlation and envelope analysis in the diagnostics of bearing faults and other cyclostationary machine signals. Mech Syst Signal Process 2001; 15(5): 945–962.

31.

Chen

Peng

, et al. High-accuracy fault feature extraction for rolling bearings under time-varying speed conditions using an iterative envelope-tracking filter. J Sound Vibr 2019; 448: 211–229.

32.

Randall

Antoni

. Rolling element bearing diagnostics—a tutorial. Mech Syst Signal Process 2011; 25(2): 485–520.

33.

Hemmati

Orfali

Gadala

. Roller bearing acoustic signature extraction by wavelet packet transform, applications in fault detection and size estimation. Appl Acoust 2016; 104: 101–118.

34.

Natili

Castellani

Astolfi

, et al. Video-tachometer methodology for wind turbine rotor speed measurement. Sensors 2020; 20(24): 7314.

35.

Dolenc

Boškoski

Juričić

. Distributed bearing fault diagnosis based on vibration analysis. Mech Syst Signal Process 2016; 66: 521–532.