Time-Frequency Analysis Based on Improved Variational Mode Decomposition and Teager Energy Operator for Rotor System Fault Diagnosis

A time-frequency analysis method based on improved variational mode decomposition and Teager energy operator (IVMD-TEO) is proposed for fault diagnosis of turbine rotor. Variational mode decomposition (VMD) can decompose a multicomponent signal into a number of band-limited monocomponent signals and can effectively avoid model mixing. To determine the number of monocomponents adaptively, VMD is improved using the correlation coefficient criterion. Firstly, IVMD algorithm is used to decompose a multicomponent signal into a number of monocompositions adaptively. Second, all the monocomponent signals’ instantaneous amplitude and instantaneous frequency are demodulated via TEO, respectively, because TEO has fast and high precision demodulation advantages to monocomponent signal. Finally, the time-frequency representation of original signal is exhibited by superposing the time-frequency representations of all the monocomponents. The analysis results of simulation signal and experimental turbine rotor in rising speed condition demonstrate that the proposed method has perfect multicomponent signal decomposition capacity and good time-frequency expression. It is a promising time-frequency analysis method for rotor fault diagnosis.


Introduction
The fault vibration signal of turbine and other larger-scale rotating machinery is often characterized by multicomponent and nonstationary on account of large span, complex structure, and variable working condition. Time-frequency analysis method can provide nonstationary signal's localization information both in time domain and in frequency domain, so it has been widely used in rotating mechanical fault diagnosis [1][2][3][4][5][6][7][8][9][10]. The key of time-frequency analysis is to find an effective method to decompose the multicomponent signal. Aiming at this problem, scholars have carried on unceasing exploration. Li and Liang used the generalized synchrosqueezing transform to detect and diagnose faults of the gearbox through analyzing the time-frequency representation [2]. Feng et al. proposed the iterative generalized synchrosqueezing transform to diagnose the lab experimental vibration signals of planetary gearboxes under nonstationary conditions [3]. A time-frequency analysis method based on the Vold-Kalman filter and higher-order energy separation was proposed by Feng et al., and it was used to diagnose the planetary gearbox fault under varying shaft speed [4]. The cracked rotor transient response was analyzed and monitored using HHT by Ramesh Babu et al. [5]. Liu et al. diagnosed the rotor rub-impact fault using Teager Huang transform [6], but these methods based on Empirical Mode Decomposition (EMD) had the problems of mode mixing, over envelope, less envelope, and so forth [7,8], which affected the analysis results. Chandra and Sekhar found that Continuous Wavelet Transform (CWT) was more preferred than HHT for noisy data and it was valid to diagnose the shaft misalignment and rotor stator rub-impact faults [8]. Li et al. distinguished the severity degree of rub-impact fault using Empirical Wavelet Transform (EWT) [9]. Tang and Pang extracted the typical fault characteristics of oil-whirl using Hilbert vibration decomposition (HVD) [10].
Dragomiretskiy and Zosso proposed a new multicomponent signal decomposition method called variational mode 2 Mathematical Problems in Engineering decomposition (VMD) [11]. This method can decompose a multicomponent signal into a given number of bandlimited intrinsic mode functions using frequency domain nonrecursive iteration pattern. VMD has high precision and can effectively avoid the problem of mode mixing [12]. It has better performance than EWT to produce the classification feature vectors of power quality disturbances [13]. Wang et al. found that the filtering characteristic of VMD was wavelet packet-like expansion via numerical simulations and used it to analyze the fixed frequency rotor rub-impact fault in time domain [14]. However, how to select the appropriate parameter of VMD, especially the number of decomposed components, is very important in practice. To solve this problem, the IVMD which uses the correlation coefficient criterion to determine the number of components adaptively is applied in this paper. Combining the fast and high precision advantages of TEO to monocomponent signal, an IVMD-TEO time-frequency analysis new method is proposed. The analysis results of simulation signal and experiment turbine rotor fault signal illustrate the effect of this method.

IVMD-TEO Time-Frequency
Analysis Method 2.1. Theory of VMD. For a multicomposition real valued signal , VMD assumes that is composed of a given number of subsignals (modes). Each mode is regarded as an amplitude-modulated and frequency-modulated (AM-FM) signal and has mostly compact frequency around a center pulsation. These modes are called intrinsic mode functions (IMFs). The constrained variational problem is the following: (1) To assess the bandwidth of the modes, the following scheme is proposed. (1) Compute the associated analytic signal by means of the Hilbert transform to obtain a unilateral frequency spectrum for each mode. (2) Shift the frequency spectrum of each mode to the baseband by mixing with an exponential tuned to the estimated center frequency. (3) Estimate the bandwidth through the 1 Gaussian smoothness of the demodulated signal, that is, the squared 2 -norm of the gradient [11].
To solve the constrained variational problem, the augmented Lagrangian is introduced and the nonconstrained variational problem is gotten by where denotes the balancing parameter of the fidelity constraint. The saddle point of (2) is the optimal solution of original problem, which can be gained using alternate direction method of multipliers (ADMM) [15]. All the modes can be obtained from (3) in the frequency domain through updating each mode and its center frequency constantly: (3) +1 is regarded as the Wiener filtering result of the residuê( ) − ∑ ̸ =̂( ), which makes the VMD algorithm much more robust to sampling and noise. The new center frequency +1 is put at the center of gravity of the corresponding mode's power spectrum, which can be updated by 2.2. Improved VMD. It is known from VMD algorithm that the number of decomposed components must be set in advance. However, the value is usually difficult to be accurately set for different conditions. In order to solve this problem, the correlation coefficient criterion is introduced to control the number adaptively in this paper, named IVMD.
The correlation coefficient is a quantitative measure of correlation and dependence. It shows the statistical relationships between two or more random variables or observed data values. It has been used in the chromosomes' similarity analysis [16], the relationship analysis between cerebral venous thrombosis and oral contraceptives in adult women [17], and the mechanical vibration signal similarity analysis [18,19]. In this paper, Pearson's correlation coefficient is used and can be measured by where ( , ) is the correlation coefficient of and , cov( , ) is the covariance between and , and and are standard deviations of and .
The correlation coefficient between the decomposed residue of VMD and original signal is computed by (5). When the absolute value is lower than the threshold (0.05 in this paper), it is thought that there is not important information in the residue and the original signal is decomposed completely and appropriately. The iterative decomposition process is finished, and the number can be determined adaptively. Engineering   3 The implementation steps of IVMD are as follows:

Teager Energy Operator.
To express time-frequency representation of each mode, the demodulation analysis is necessary to each IMF. TEO is a kind of nonlinear local differential operator without involving integral transform. It can improve the estimation of the instantaneous spectrum characteristics of the vibration data under certain conditions [7]. It has been widely used in bearing fault diagnosis in recent years [20][21][22].
For the AM-FM signal ( ) = ( ) cos[ ( )], the corresponding TEO is defined as where( ) and( ) are its first and second derivatives, respectively. In the discrete time case, the time derivatives in (6) can be approximated to TEO offers excellent time resolution and demodulation speed because only three sample points are required for the energy computation at each time instant. As a result, the instantaneous change in amplitude and frequency of signal can be quickly tracked. The AM and FM information of vibration signals can be extracted using, respectively [23],

IVMD-TEO Time-Frequency Analysis Method.
Combining the advantages of IVMD and TEO, a time-frequency analysis method named IVMD-TEO is proposed in this paper. First of all, IVMD algorithm is used to decompose a multicomponent signal into IMFs adaptively. Then, the instantaneous amplitude and instantaneous frequency of each IMF are demodulated by TEO. Finally, the timefrequency representation of the original signal is obtained by superposing the time-frequency representations of all the monocomponents, and the fault types are diagnosed by analyzing the fault characteristic frequency components showed in the time-frequency representation diagram. In the process of VMD, the balancing parameter is a parameter to balance the fidelity. The bigger value is suitable to filter low-frequency component (e.g., trend term) and the smaller value is suitable to detect impacts from a time series (or signal) [24], so the smaller value ( = 500) is used in this paper.
The time-frequency analysis process of the proposed IVMD-TEO method for rotor system is shown in Figure 1, and the specific implementation steps are as follows.
(1) Gather original signal data from experimental rotor rig or practical rotor system.
(3) Demodulate the instantaneous amplitude and instantaneous frequency of each IMF using TEO.
The sampling frequency is 2048 Hz, and the analysis points are 2048 points. The waveform and spectrum of ( ) are shown in Figure 2. The spectrum basically shows each component. First, ( ) is decomposed by the IVMD algorithm, and seven IMFs are extracted and shown in    IMF1  IMF2  IMF3  IMF4  IMF5  IMF6  IMF7     Although there is a little interference, all time-frequency representation of each component is clearly segregated.
For comparison, the time-frequency spectrum of simulation signal analyzed by the HHT is shown in Figure 5. It can be seen from Figure 5 that the frequency modulation phenomenon is confused. Two FM components in 1 ( ) are extracted partly, but the low disturbance frequency is emerged. There are frequency components in the vicinity of 200 Hz, but the detailed composition cannot be expressed. The analysis result is bad.

Experiment Rotor System Vibration Signal Analysis
The turbine rotor faults are simulated on the Bently RK-4 test rig ( Figure 6). The test rig is an oil-fluid supporting rotor-bearing system. An electrical motor is attached to the shaft that carries two discs by a flexible coupling. A phase sensor closed to the coupling is fixed to measure the rotating speed. Two eddy current sensors, located near the rub copper thimble, are fixed by holder in the middle of test rig. A safety device is designed to ensure the experiment rotor safety. An  oil film journal bearing is part of an assembly connected to the oil pump system. The bearing is equipped with two proximity eddy current sensors in the horizontal and vertical directions. Rub-impact signal and oil film instability signal are gathered by sensor 1 and sensor 2, respectively using ZonicBook/618E vibration measurement system [6]. The sampling frequency is 5120 Hz. In order to reduce the amount of analysis data, the following analysis frequency is 512 Hz. That is to say, one point in ten is selected in sequence.

Rub-Impact Fault
Analysis. The copper thimble that has a small separation from the shaft is fixed by holder in the middle of test rig and rubs against the high-speed rotating shaft slightly. The rub-impact signal when shaft speed is 1620 r/min∼1740 r/min is analyzed and its waveform and spectrum are shown in Figure 7. The waveform amplitude becomes bigger with time, which only shows the rising rotational speed process. The spectrum only shows that the continuous rotational frequency is 27∼29 Hz. It is difficult to diagnose the slight rub-impact fault accurately. Seven IMFs obtained by IVMD for slight rub-impact signal are shown in Figure 8. In apparent periodicity components, IMF5 is corresponding to the rotational frequency (1 ), IMF3 is corresponding to 2 , IMF2 is corresponding to 3 , and IMF1 is corresponding to 4 . The correlation coefficients of the former six IMFs are bigger than the threshold value. Their time-frequency representations demodulated by TEO are shown in Figure 9. It can be seen that the amplitude of rotational frequency is the maximum and 2 ∼4 frequency components are also extracted by IVMD-TEO. At the same time, the high harmonics frequency fluctuates in the range of 500th to 1000th point. It illustrates that the rubimpact process is instable. All of the frequency characteristics above accord with the rub-impact fault characteristics [25], so the occurrence of rotor slight rub-impact fault can be diagnosed accurately.
In order to compare the analysis results, the timefrequency analysis result for slight rub-impact signal by HHT is shown in Figure 10. It only shows demodulated rotational frequency and cannot show the high harmonics frequency. The slight rub-impact characteristics are not extracted. The analysis effect is not as good as the effect of IVMD-TEO.

Oil Film Instability Analysis.
The oil film instability of the rotor system usually consists of two stages that are the nearly half speed oil-whirl and oil-whip. The monitor and control to initial oil-whirl are significant for turbine rotor to avoid great loss. The initial oil-whirl and the first-order oil-whip spectrum characteristics of the test rig are analyzed subsequently.
In the process of oil film instability, the oil-whirl occurs at the rising speed of 1620 r/min∼1860 r/min, and the waveform and spectrum of the oil-whirl are shown in Figure 11. The rising speed can be seen in the waveform. In the spectrum, the oil-whirl frequency ( wl ) (12∼14 Hz) and 1 (27∼31 Hz) can be seen. The wl value is about half value of 1 and higher amplitude than 1 , but all of them lack the corresponding time information.
Six IMFs analyzed by IVMD to initial oil-whirl signal are shown in Figure 12. In apparent periodicity components, IMF1 is corresponding to wl . IMF2 and IMF3 are corresponding to 1 . The correlation coefficients of the former five IMFs are bigger than the threshold value. Their timefrequency representations demodulated by TEO, respectively, are shown in Figure 13. It can be seen that a little oilwhirl amplitude appears before 1000th point and the strong oil-whirl amplitude accompanied by amplitude modulation phenomenon appears after the 1000th point. The process Mathematical Problems in Engineering shows that the initial oil-whirl is instability process. The timefrequency representation clearly shows the whole process of the occurrence and development of oil-whirl. A phenomenon that a frequency component whose frequency is about 2 before the 1200th point and becomes 1 + wl after the 1200th point appears too. This detailed changing process suggests the high precision decomposition of IVMD-TEO.
For comparison, the time-frequency spectrum of oilwhirl signal analyzed by HHT is shown in Figure 14. The frequency mixing phenomenon is obvious, which cannot clearly show the whole process of the oil-whirl.
When the rotational speed is up to 2760 r/min∼ 3000 r/min, the rotor system test rig shows biggish vibration because of the oil-whip. Waveform and spectrum of the oil-whip signal are shown in Figure 15. The rotational frequency (1 ) (46∼50 Hz) and the oil-whip frequency ( wp ) (22∼24 Hz) can be seen in the spectrum, but they are continuous spectrum and lack corresponding time information. The oil-whip can be diagnosed preliminary.
Eleven IMFs obtained by IVMD to oil-whip signal are shown in Figure 16. The correlation coefficients of the former ten IMFs are bigger than the threshold value. Their timefrequency representations demodulated by TEO, respectively, are shown in Figure 17. Besides the main 1 and wp frequency components, the related combination frequency components wp /2, 1 − wp /2, 1 + wp , 2 − wp /2, 2 , and 2 + wp are also extracted clearly. The oil-whip can be diagnosed accurately and reliably [25]. The analysis results also indicate the decomposition accuracy of the IVMD-TEO method.
To compare the analysis results, the time-frequency spectrum of oil-whip signal analyzed by HHT is shown in Figure 18. The 1 , wp , and wp /2 frequency components are observed only and other frequency components cannot be extracted out. The analysis results are bad to IVMD-TEO method.
Qualitative analysis results in Figures 9, 13, and 17 for rotor faults show better effectiveness of the proposed method than HHT. The quantitative description in main fault   Table 2. For slight rub-impact fault, HHT method cannot extract fault frequency components; however, the proposed method can extract 2 ∼4 fault frequency components and their amplitudes are about 4 m. For oil film instability fault (oil-whirl and oil-whip), most of the main fault frequency amplitudes extracted by proposed method, such as 1 oil-whirl, wl , and wp , are improved within 12.3%∼35.5% comparing with HHT except 1 oil-whip frequency. The new extracted detail fault frequency amplitudes which cannot be extracted by HHT are also shown in Table 2. Both qualitative analysis and quantitative analysis can illustrate the advantage of IVMD-TEO method.

Conclusion
VMD is improved to determine the number of decomposition components adaptively using the correlation coefficient criterion. Integrating the advantages of IVMD and TEO, a time-frequency analysis method named IVMD-TEO is proposed in this paper. The analysis results of experimental turbine rotor signals under variable speed conditions illustrate that the method can show the time-frequency characteristics of slight rub-impact clearly, monitor the process of oil-whirl occurrence and development exactly, and extract the complicated time-frequency components of oil-whip accurately. It enhances the precision and accuracy of rotor fault diagnosis. Comparative analysis results show that the proposed method is more effective than HHT, so this method has important engineering application value for rotor fault diagnosis.