Research on a Novel Improved Adaptive Variational Mode Decomposition Method in Rotor Fault Diagnosis

: Variational mode decomposition (VMD) with a non-recursive and narrow-band ﬁltering nature is a promising time-frequency analysis tool, which can deal e ﬀ ectively with a non-stationary and complicated compound signal. Nevertheless, the factitious parameter setting in VMD is closely related to its decomposability. Moreover, VMD has a certain endpoint e ﬀ ect phenomenon. Hence, to overcome these drawbacks, this paper presents a novel time-frequency analysis algorithm termed as improved adaptive variational mode decomposition (IAVMD) for rotor fault diagnosis. First, a waveform matching extension is employed to preprocess the left and right boundaries of the raw compound signal instead of mirroring the extreme extension. Then, a grey wolf optimization (GWO) algorithm is employed to determine the inside parameters ( ˆ α , K ) of VMD, where the minimization of the mean of weighted sparseness kurtosis (WSK) is regarded as the optimized target. Meanwhile, VMD with the optimized parameters is used to decompose the preprocessed signal into several mono-component signals. Finally, a Teager energy operator (TEO) with a favorable demodulation performance is conducted to e ﬃ ciently estimate the instantaneous characteristics of each mono-component signal, which is aimed at obtaining the ultimate time-frequency representation (TFR). The e ﬃ cacy of the presented approach is veriﬁed by applying the simulated data and experimental rotor vibration data. The results indicate that our approach can provide a precise diagnosis result, and it exhibits the patterns of time-varying frequency more explicitly than some existing congeneric methods do (e.g., local mean decomposition (LMD), empirical mode decomposition (EMD) and wavelet transform (WT) ). analyzing multi-component signals.


Introduction
As everyone knows, most mechanical fault signals have the characteristics of being nonlinear and non-stationary; it is difficult to accurately reveal the fault feature frequency applying directly fast Fourier transform (FFT) in data processing [1]. Hence, a host of studies have focused on how to improve the ability of fault information excavation, mainly classified as three types, namely the time domain, frequency domain, and time-frequency domain [2][3][4]. Among these approaches, the time-frequency analysis (TFA) tool has gained increasing interest for capturing intrinsic feature information from the waveform matching extension to obtain the extended signal. After that, GWO is introduced to determine adaptively the inside parameters of VMD. Moreover, we employ the parameter-optimized VMD to divide the extended signal into several IMF ingredients and truncate the extended segment to obtain the final decomposition results. Finally, TEO demodulation for all IMFs is conducted to achieve the corresponding instantaneous frequency (IF), and the IF of all IMFs is painted together to obtain our desired time-frequency representation (TFR). One obvious merit of the proposed IAVMD is that it does not depend on any parameterization, which makes this method more adaptive to actual situations. Crucially, the study of simulated and actual laboratory data validates the effectiveness of our approach.
The remainder of this paper is as follows. Section 2 reviews the basic theory of VMD. Section 3 introduces the theory of the waveform matching extension method and shows the parameter optimization process of VMD based on GWO. Meanwhile, the overall realization process of the proposed IAVMD is illustrated in Section 3. Sections 4 and 5 prove the effectiveness of the presented approach through a simulation analysis and experimental examples, respectively. Section 6 draws some conclusions and provides some future works.

VMD
Briefly speaking, the main idea of VMD is to availably solve the constrained variational issue (see Equation (1)): where u k indicates the k-th mode, and ω k means the mid-frequency corresponding to the k-th mode. First, the penalty factorα and Lagrangian multiplier λ(t) are added into Equation (1), so the variational problem (i.e., the augmented Lagrangian) is rewritten as [55]: Second, several sub-signals (see Equation (3)) can be obtained by utilizing the alternate direction method of multipliers:û n+1 k (ω) =f (ω) − i kûi (ω) +λ where the center-frequencies ω k of each mode are updated by: At last, determine if the stop conditions (see Equation (5)) are met, and if it meets the conditions, stop VMD running and output the entire decomposition results: According to Ref. [55], the threshold ε is usually set as 10 −6 . Additional details of VMD can be found in the original article [14].

Waveform Matching Extension
In the waveform matching extension, we first look for all available similar waveforms, then pick out the best matching waveform within them. Last, we connect the waveband before/after the best matching waveform to both left and right ends of the raw waveform to finish the task of boundary processing. Figure 1 shows the flowchart of the waveform matching extension. For one compound waveform x(t), the waveform matching extension procedure is described below [56].
(1) Find all extreme points and boundary points of the original signal x(t) through the peak searching method, and determine the characteristic waveform. Take Figure 2 for an example, M i (i = 1, 2, 3, · · · ) is the maximum value of x(t), corresponding to the time tm i , N i is the minimum value of x(t), corresponding to the time tn i , and S 1 is the left boundary point of the signal x(t). We regard the triangular waveform S 1 − M 1 − N 1 of the original signal x(t) as the characteristic waveform to search the optimum matching waveform.
(2) Seek for all matching waveforms resembling the characteristic waveform S 1 − M 1 − N 1 on the basis of the time ts i corresponding to the starting value S i of the matching waveforms, where the time ts i of the horizontal axis is achieved by the linear interpolation, which is expressed as: (3) Estimate the matching error generated in searching all similar waveforms according to the following Equation (7): where M i+1 − M 2 represents the trend term of the similar waveforms, which can reveal the relative extremum position of all similar waveforms.
(4) Compare the matching error, determine the smallest error, and consider the similar waveform S 1 − M 1 − N 1 containing the smallest error as the best matching waveform. (5) Use the data before the best matching waveform to extend the left ends of the raw waveform x(t).
(6) Apply the same way to extend the right ends of the raw waveform x(t), and finally obtain the extended waveform. This indicates that the waveform matching extension procedure of the original signal x(t) is finished in this stage. That is, the left and right boundaries of the original signal x(t) have been processed thoroughly.
Here, two representative examples for a signal with an abrupt change are utilized to describe intuitively the procedure of waveform matching extension. Without loss of generality, the number of data points extended in both left and right ends are selected as 100 points. Figure 3a,b shows the processing results achieved when applying the waveform matching extension for a cycle shock pulse signal and a modulated compound signal, respectively. Obviously, the waveform extended in both the left and right ends is in line with the variation trend of the raw signal, which indicates that the waveform matching extension is feasible in analyzing the modulated compound signal.

Parameter Optimization of VMD Using GWO
Two important variables (i.e., the penalty factorα and mode number K) are required to define in advance when we employ VMD to a process complex compound signal. Some previous studies have shown that the two variables (α and K) have a great influence on the decomposition performance of VMD. On the one hand, the biggerα represents the smaller bandwidth of modes. Otherwise, the smallerα represents the greater bandwidth of modes. On the other hand, if the mode number K is too big, it easily leads to the over decomposition phenomenon (i.e., information redundancy), while this will cause the under-decomposition phenomenon (i.e., information loss) when the mode number K is too small. Hence, how to pick out the appropriate variables (α, K), is a challenge when applying VMD to deal with complex compound signals. To address this issue, a new nature-inspired algorithm called GWO is introduced to determine and select automatically the combination parameters (α, K) of VMD in this subsection due to its excellent global optimization performance. Obviously, the objective function should be constructed when we use GWO to optimize the parameters of VMD. At present, there are many indicators (e.g., smoothness index, entropy, sparseness, kurtosis, and correlation coefficient) in the analysis of the mechanical vibration signal, but it is not optimal to use any of their single indicators to measure vibration signal characteristics. This indicates that the combination of two or more of their single indicators is valuable for practical vibration detection. Zhang et al. [57] presented a new combination index named weighted kurtosis (WK) as the objective function of their method to optimize the key parameters of VMD, but the proposed weighted kurtosis does not consider some abnormal impacts (i.e., the outlier) with disperse distributions. Hence, to solve this problem, a comprehensive evaluation indicator named weighted sparseness kurtosis (WSK) is formulated in this paper to work as the objective function of the parameter optimization process of VMD, which is defined by: where S is the sparseness of the signal x(n), N is the length of the signal x(n), ρ represents the Pearson correlation coefficient between two signals (i.e., x and y) and satisfies ρ ≤ 1, and E[·] denotes the expectation operator. This indicates that the correlation coefficient ρ amounts to the weight of the sparseness kurtosis indicator. That is, Equation (8) is reputed essentially as the weighted sparseness kurtosis indicator in this paper. Note that in Equation (8) the sparseness is a statistical indicator reflecting the amplitude distribution of the signal, especially for an abnormal impact amplitude with discrete distribution, while kurtosis is very sensitive to the periodic impulse series of vibration signals. Besides, the correlation coefficient can measure the similarity between two vibration signals.
To prove the efficacy of the presented WSK indicator, comparisons between four indicators (i.e., kurtosis, entropy, WK, and WSK) are conducted by analyzing a simulation signal. The simulation signal shown in Figure 4a is composed of the periodic impulse, harmonic signal, and random noise. Figure 4b,c shows the simulation signal containing a single outlier and multiple outlier, respectively. Table 1 gives the comparison results of different indicators for the simulation signal. From the chart we know that the proposed WSK indicator has the smallest change when the outlier of the signal increases gradually. This implies that the proposed WSK indicator is more reliable and robust to the outlier of the signal compared with the single indicator and WK. Hence, through the combination of three indicators (i.e., the sparseness, kurtosis, and correlation coefficient), the WSK indicator is more comprehensive and has more advantages in measuring vibration features. fitness value ) ( i WSK mean is small enough. If that is the case, stop the iteration and output the best wolves best X (i.e., the best parameters best alpha x and best K x of VMD). Otherwise, let t = t+1, back to step (2) to keep working.
In short, this section has two contributions. First, the WSK indicator is proposed and is more reliable and robust to the outlier of the signal compared with the previous WK and three single indicators (i.e., the sparseness, kurtosis, and correlation coefficient). Second, two key parameters ( αˆ, K) of VMD can be determined automatically by using the GWO method, where the mean of the proposed WSK for all modes is regarded as the objective function of the optimization procedure.    Obviously, a larger WSK implies a richer impact feature and a better signal sparseness and decomposition results. Here, to further ensure a global optimal solution, we deem the mean of the WSK of modes acquired by VMD as an objective function of the optimization procedure. Concretely, the goal of the parameter optimization process of VMD using GWO is visually interpreted as an effective search for the minimum value of the opposite of the objective function, as shown in Equation (12): where K denotes the mode number,α denotes the penalty factor, mean(WSK i ) is the mean value of WSK, and WSK i is the WSK corresponding to the i-th mode. GWO has a strict population hierarchy, which is mainly divided into four layers (see Figure 5), where the first layer is α (i.e., leaders/the optimal solution), the second layer is β (i.e., subordinate wolves-advisors/the suboptimal solution), the third layer is δ (i.e., scouts, sentinels, elders, hunters, and caretakers/the third optimal solution), and the last layer is ω (i.e., lowest ranking wolves-scapegoats/other solutions). Moreover, GWO is applied for parameter optimization by implementing several steps (e.g., hunting, searching for prey, encircling prey, and attacking prey). Figure 6 shows the flowchart of the parameter optimization of VMD using the GWO method. The general procedure of VMD parameter optimization is elaborated in detail below: (1) Initialize the grey wolf population X i (i = 1, 2, · · · , m), and set the fitness function and the parameters of the GWO algorithm. Generally speaking, the higher the number of wolves and iterations is, the better the optimization results are, but the longer the computation time is; thus, according to Ref. [42], we define the maximum number of wolves m = 30 and maximum iterations T = 10, which is aimed at achieving a tradeoff between the optimization results and computational efficiency. Besides, due to the pre-optimization, the objective only involves two variables (α, K), and thus each wolf is expressed as X i = (x K , x alpha ), i = 1, 2, · · · , m, where x K and x alpha represent the penalty factor and mode number, respectively. Due to fault feature information of the practical signal usually spread over the first few high-frequency components obtained by VMD, the first two to eight components are usually selected and analyzed. Besides, the feature information of the main components obtained by VMD is distributed over an appropriate bandwidth, which means that the parameterα is supposed to stay within an empirical range and is neither too big nor too small compared to the default of VMD. For all those reasons, here the search ranges of the two parametersα and K are empirically set as [2,8] and [100,2000], respectively.
(2) Calculate the fitness value mean(WSK i ) of each wolf X i (i = 1, 2, · · · , m), look for the minimum fitness value (i.e., min mean(WSK i ) ) between the four agents (i.e., the fitness value of the individual grey wolf α, β, δ, and ω) and record the current best wolves X best .
(3) Update the position of wolves in terms of the gray wolf movement pattern shown in Equation (13): where A is the convergence factor and meets A = 2d · r 1 − d, C is the swing factor and meets C = 2 · r 2 , d is the range control parameter which linearly decays from 2 to 0 over the whole iteration, and r 1 and r 2 are the random numbers between 0 and 1. i . (5) Determine whether the stopping condition is satisfied. Concretely, decide whether the current iterations are lower than the maximum iterations (i.e., t ≤ T) or whether the opposite of the fitness value mean(WSK i ) is small enough. If that is the case, stop the iteration and output the best wolves X best (i.e., the best parameters x best alpha and x best K of VMD). Otherwise, let t = t + 1, back to step (2) to keep working.
In short, this section has two contributions. First, the WSK indicator is proposed and is more reliable and robust to the outlier of the signal compared with the previous WK and three single indicators (i.e., the sparseness, kurtosis, and correlation coefficient). Second, two key parameters (α, K) of VMD can be determined automatically by using the GWO method, where the mean of the proposed WSK for all modes is regarded as the objective function of the optimization procedure.

Teager Energy Operator
For a given continuous-time signal with modulated information (t x  mean the first and second order derivative functions of x(t), respectively. For a discrete signal x(n), TEO is modified into where n means the length of x(n). The instantaneous amplitudes IA and IF of x(n) are capable of being acquired using the discrete energy separation method (see the DESA-2 function in MATLAB), as shown in Equations (16) and (17): represents the instantaneous amplitude, and f(n) represents the instantaneous frequency. One major advantage of TEO compared with HT is that the calculated IA and IF have fewer errors.

Teager Energy Operator
For a given continuous-time signal with modulated information ..
x(t) mean the first and second order derivative functions of x(t), respectively. For a discrete signal x(n), TEO is modified into ψ[x(n)] = x 2 (n) − x(n − 1)x(n + 1), where n means the length of x(n). The instantaneous amplitudes IA and IF of x(n) are capable of being acquired using the discrete energy separation method (see the DESA-2 function in MATLAB), as shown in Equations (16) and (17): where a(n) represents the instantaneous amplitude, and f (n) represents the instantaneous frequency. One major advantage of TEO compared with HT is that the calculated IA and IF have fewer errors.

The Proposed IAVMD Method
As mentioned before, the advantages of waveform matching extension are utilized by this paper for signal pretreatment. The inside parameters (i.e.,α and K) of VMD are determined beforehand by GWO, which can make the algorithm adaptive and practicable. One advantage of TEO is that the edge errors of HT can be avoided in the signal demodulation process. Hence, by merging together several procedures, a novel TFA method called IAVMD is formulated in this paper. Figure 7 shows the flowchart of the proposed method. For a multi-component signal, the realization process of the proposed IAVMD method is as follows: Step 1: Preparations of IAVMD. Use the waveform matching extension illustrated in Figure 1 to stretch the left and right boundary points of the multi-component signal.
Step 2: Parameter optimization automatically. Introduce GWO to adaptively search the optimal inside parameters (α, K) of VMD according to steps (1) to (5) of Section 3.2.
Step 3: Signal decomposition. Apply the parameter optimization VMD method to decompose the extended signal acquired by step 1 into several IMF components shown in Equation (3), and truncate the extended data points of each IMF component to obtain the final decomposition results.
Step 4: Signal demodulation. Perform TEO to estimate the IF and IA of each IMF component, and plot the IF of all IMF components together to get the ultimate time-frequency representations of the original signal. Observe and identify the frequency contents displayed in TFR, and output the fault diagnosis results. beforehand by GWO, which can make the algorithm adaptive and practicable. One advantage of TEO is that the edge errors of HT can be avoided in the signal demodulation process. Hence, by merging together several procedures, a novel TFA method called IAVMD is formulated in this paper. Figure 7 shows the flowchart of the proposed method. For a multi-component signal, the realization process of the proposed IAVMD method is as follows: Step 1: Preparations of IAVMD. Use the waveform matching extension illustrated in Figure 1 to stretch the left and right boundary points of the multi-component signal.
Step 2: Parameter optimization automatically. Introduce GWO to adaptively search the optimal inside parameters ( αˆ, K) of VMD according to steps (1) to (5) of Section 3.2.
Step 3: Signal decomposition. Apply the parameter optimization VMD method to decompose the extended signal acquired by step 1 into several IMF components shown in Equation (3), and truncate the extended data points of each IMF component to obtain the final decomposition results.
Step 4: Signal demodulation. Perform TEO to estimate the IF and IA of each IMF component, and plot the IF of all IMF components together to get the ultimate time-frequency representations of the original signal. Observe and identify the frequency contents displayed in TFR, and output the fault diagnosis results.

Case 1: Sinusoidal Superposition Signal
Considering that the real rotor vibration signal is generally expressed as multiple harmonic superposition signals (e.g., the linear superposition of fundamental frequency, fractional frequency, and ultraharmonics), to confirm the availability of the presented approach we here employ a numerical signal ) (t y consisting of four sine waves (their frequencies are, respectively, 30 Hz, 60 Hz, 100 Hz, and 150 Hz) to simulate the rotor fault signal, which is considered as follows:

Case 1: Sinusoidal Superposition Signal
Considering that the real rotor vibration signal is generally expressed as multiple harmonic superposition signals (e.g., the linear superposition of fundamental frequency, fractional frequency, and ultraharmonics), to confirm the availability of the presented approach we here employ a numerical signal y(t) consisting of four sine waves (their frequencies are, respectively, 30 Hz, 60 Hz, 100 Hz, and 150 Hz) to simulate the rotor fault signal, which is considered as follows: where x 1 (t) is used to simulate the fundamental frequency component of the rotor fault signal, x 2 (t) is used to simulate the double frequency component of the rotor fault signal, x 3 (t) is used to simulate the fractional frequency component of the rotor fault signal, and x 4 (t) is used to simulate the higher harmonic component of the rotor fault signal. For the simulation signal y(t), its sampling frequency and sampling number are set as 3000 Hz and 1500, respectively. Figure 8 draws the numerical signal x(t) and its FFT spectrum.
In the IAVMD method, the numerical signal y(t) is first preprocessed by waveform matching extension, after which the optimal combined parameters (α, K) of VMD obtained using GWO are (1986,4). That is, VMD withα = 1986 and K = 4 is applied to analyze the numerical signal y(t), and the ultimate decomposition results are shown in Figure 9. Note that, except for the input signal and two key parameters (α, K), another four parameters (e.g., time-step of the dual ascent τ, DC part, uniform initialization of omegas, and tolerance of the convergence criterion ε) of VMD have a fixed value, set as 0, 0, 0, and 10 −6 , respectively. Seen from Figure 9, four frequency components of the numerical signal y(t) are separated efficiently. Figure 10a,b plots the two-dimensional and three-dimensional TFR based on IAVMD, respectively. One can clearly see that the time-frequency trajectory obtained by IAVMD is very clear in Figure 10, which indicates that the proposed IAVMD can exactly decompose the sinusoidal superposition signal containing multiple components.    As a comparison, particle swarm optimization-based VMD (OVMD) and three other techniques (i.e., LMD, EMD, and WT) are devoted to disposing the same signal y(t). The optimization range of OVMD is the same as for the proposed method. Here, the optimized parameters of OVMD are respectively determined asα = 1120 and K = 4 by using the PSO algorithm. The decomposition results acquired by OVMD, LMD, EMD, and WT are shown in Figure 11a-d, respectively. Besides, the acquired TFR using four techniques (i.e., OVMD, LMD, EMD, and WT) are shown in Figure 12a-d, respectively. Note that, in the comparison process, to ensure a fair comparison, the TFR of four methods (i.e., OVMD, LMD, EMD, and WT) are obtained by TEO, where WT adopts the 'db1' wavelet function containing a good orthogonality and where three decomposition layers are performed (i.e., four reconstructed wavelet coefficients can be obtained for implementing the time-frequency analysis). As shown in Figure 12a, four sub-components are extracted effectively, but they suffer from certain boundary fluctuations (see the red dotted line). From Figure 12b,c, one can observe that the first two IMFs suffer from serious mode aliasing and that the signal representation in the time-frequency plane is highly dispersed. That is, the IMFs derived from LMD and EMD cannot coincide well with the real components. As can be seen from Figure 12d, many discrete points appear on the time-frequency trajectory. That is, the time-frequency trajectory obtained by WT is very scattered and unfocused, so it is difficult to effectively identify the frequency components. This means that the time-frequency resolution of WT is not as good as those provided in Figure 10a. Consequently, through the comparison of Figure 12a-d, it is obvious that the potential of our presented IAVMD approach in analyzing non-stationary compound signals is preliminarily verified. Besides, the decomposition performance of our presented IAVMD is superior to that of some classical algorithms mentioned in this case.

Case 2: AM-FM Superposition Signal
To further prove the efficacy of the proposed IAVMD method, we analyze and investigate a numerical signal z(t) shown in Equation (19), which is composed of three AM-FM signals and a stochastic noise: q 1 (t) = (0.8 + 0.8 sin(30πt)) cos(800πt + 0.6 sin(24πt)) q 2 (t) = (1 + sin(20πt)) cos(400πt + 0.6 sin(16πt)) q 3 (t) = (0.5 + 0.5 sin(10πt)) cos(200πt + 0.6 sin(8πt)) where three components q 1 (t), q 2 (t), and q 3 (t) are AM-FM signals containing the main frequency of 400 Hz, 200 Hz, and 100 Hz, respectively. The component q 4 (t) is a stochastic noise added in the numerical signal z(t), and its value is generated by the randn() function in MATLAB. The sampling frequency and sampling number of the simulation signal z(t) are 3000 Hz and 3000, respectively. Figure 13 shows the simulation signal z(t) and its corresponding FFT spectrum.   The proposed IAVMD method is adopted to analyze the simulation signal z(t). First, waveform the matching extension is used to preprocess the simulation signal z(t), after which the GWO method is employed to determine automatically the VMD parametersα = 1023 and K = 3. Finally, VMD containing the optimal parameters is adopted to analyze the simulation signal z(t), and the ultimate decomposition results are displayed in Figure 14. The corresponding time-frequency representations are shown in Figure 15. Figure 15a,b shows the two-dimensional and three-dimensional TFR based on IAVMD, respectively. It can be clearly seen from Figures 14 and 15 that three components obtained by the proposed method are very close to three real AM-FM components of the original signal, which indicates that the IAVMD method is effective in analyzing multi-component signals with modulation phenomena.      As a contrast, four methods (i.e., OVMD, LMD, EMD, and WT) are employed to deal with the same simulation signal z(t). Note that to consolidate the results of the proposed method two key parameters of OVMD are optimized asα = 751 and K = 3. Besides, the implementation of the TFR of other three methods is the same as that of simulation 1. The TFR obtained by the four methods are plotted in Figure 16a-d, respectively. Seen from Figure 16a, in the first two components obtained by OVMD, there is some deviation from the true component. That is, the first two components have a low resolution in TFR based on OVMD. As shown in Figure 16b,c, the time-frequency trajectory is not clear in the TFR based on EMD and LMD, and their time-frequency aggregation is not as good as in the proposed method. As shown in Figure 16d, the time-frequency trajectories are scattered, and three frequency components of the original signal cannot be identified in TFR based on WT. Therefore, the above comparison results further validate the effectiveness of the proposed method in multi-component signal analysis.

Experiment Platform and Data Description
To further verify the presented IAVMD method, we process the experimental data collected from the rotor-bearing system in the vibration measurement laboratory at North China Electric Power University (NCEPU). The experimental device used in the testing is a Bently RK-4 test bench In short, the contribution of this section is that two simulation cases are adopted to verify the effectiveness of our proposed IAVMD method. Besides, the superiority of the IAVMD method is demonstrated via the comparisons of various similar methods (e.g., OVMD, LMD, EMD, and WT).

Experiment Platform and Data Description
To further verify the presented IAVMD method, we process the experimental data collected from the rotor-bearing system in the vibration measurement laboratory at North China Electric Power University (NCEPU). The experimental device used in the testing is a Bently RK-4 test bench (see Figure 17), which mainly consists of a signal front adapter, speed-adjusting, and a bearing oil pump system. The specification of the data collection equipment is ZonicBook/618E. The eddy current sensor and key phasor transducer are mounted on both sides of the disk to collect the rotor experimental vibration signal and rotational speed signal, respectively. First, plastic rods near the revolving shaft were adopted to produce the slight local rub-impact fault under 1790 r/min (i.e., rotating frequency f 1 = 29.83 Hz). Second, the rotor oil-whirl and oil-whip fault were generated by adjusting the device such as the oil pump, preload scaffold, and motor speed controller when the rotation speed was respectively set at 2600 r/min (i.e., rotating frequency f 2 = 43.33 Hz) and 4500 r/min (i.e., rotating frequency f 3 = 75 Hz). Note that the critical speed of oil-whip is about 1800 r/min (i.e., oil-whip frequency f w = 30 Hz). Specifically, in this experiment, the rotor system operates at above twice the critical speeds, resulting in rotor oil-whip. The sampling frequency and data length during testing were set at 1280 Hz and 1024 points, respectively. Generally speaking, when the rotor shows a rub-impact fault, the expected main fault feature frequencies in the frequency spectrum are the rotating frequency f 1 , its harmonics (e.g., 2 f 1 , 3 f 1 , and 4 f 1 ) and sub-harmonics (e.g., 0.5 f 1 , 1.5 f 1 , and 2.5 f 1 ). Besides, as we know the rotor oil-whirl is also known as half-speed whirl. That is, the expected main fault frequencies in the frequency spectrum are half of the rotating frequency (0.5 f 2 ) and its related harmonics (e.g., f 2 , 1.5 f 2 , 2 f 2 , 2.5 f 2 , and 3 f 2 ) when the rotor shows oil-whirl. Similarly, the expected main fault frequencies in the frequency spectrum are the oil-whip frequency f w , rotating frequency f 3 , and their related harmonics (e.g., 0.5, f 3 − f w , f 3 + f w , 2 f w , + 1.5 f w , and 2 f 3 ) when the rotor shows oil-whip.  Figure 18a,b shows the waveform and FFT spectrum of the rotor early rub-impact signal, respectively. As seen in Figure 18a,b, the rotor slight rub-impact signal is expressed approximatively by a sine wave, the frequency component of 29.5 Hz (approximately equal to the rotating frequency 1 f ) is very clear in the FFT spectrum, but the amplitude at the harmonics is very small. That is to say, it is not easy to judge whether rub-impact fault occurs when directly using the FFT spectrum.

Case 1: Rotor Rub-Impact Fault Detection
Five approaches (i.e., IAVMD, OVMD, LMD, EMD, and WT) are applied to decompose the rub-impact fault signal, respectively. Figure 18c-g respectively depicts the TFR resulting from the five methods. The optimal inside parameters in IAVMD are selected as αˆ = 1177 and K = 4 through using the GWO algorithm. The optimal parameters of OVMD are determined as αˆ = 536 and K = 4 through using the PSO algorithm. As seen in Figure 18c, the fundamental frequency of 29.5 Hz and its various harmonics (59 Hz, 147.5 Hz, and 265.5 Hz) can be identified. In addition, the high frequency is relatively ambiguous due to a resolution problem, but the spectral line is very clear in the low frequency. This indicates that our presented IAVMD approach can detect early characteristics of rotor rub-impact faults.
As shown in Figure 18d, we can easily find the frequency components of 29.5 Hz and its  Figure 18a,b shows the waveform and FFT spectrum of the rotor early rub-impact signal, respectively. As seen in Figure 18a,b, the rotor slight rub-impact signal is expressed approximatively by a sine wave, the frequency component of 29.5 Hz (approximately equal to the rotating frequency f 1 ) is very clear in the FFT spectrum, but the amplitude at the harmonics is very small. That is to say, it is not easy to judge whether rub-impact fault occurs when directly using the FFT spectrum. marked with a red line). In Figure 18e,f, the instantaneous frequencies of 29.5 Hz can be found, but their time-frequency trajectory is obscure. It is obvious from Figure 18g that the TFR obtained by WT can dimly extract the rotating frequency of 29.5 Hz, but the frequency points are very scattered and their time-frequency resolution is clearly worse. The contrastive analysis demonstrated that the efficacy of IAVMD in rotor rub-impact fault detection is superior to that of the other four approaches (i.e., OVMD, LMD, EMD, and WT).  Figure 19a,b shows the waveform and FFT spectrum of the rotor oil-whirl signal, respectively. Figure 19a shows that the rotor oil-whirl signal is represented by a sum of sine waves, and that in the FFT spectrum there is a prominent amplitude at the rotating frequency 2 f of 43.33 Hz. Besides, the amplitude of the oil-whirl frequency of 21.5 Hz (equal to about 0.5 2 f ) is greater than that of the rotating frequency, but the amplitude of the high frequency is not obvious. The oil-whirl signal is analyzed by five methods (i.e., IAVMD, OVMD, LMD, EMD, and WT), respectively. Figure 19c-g provides the TFR derived from the five methods, respectively. In IAVMD, the optimal parameters obtained by the GWO algorithm are αˆ = 1919 and K = 4, respectively. In OVMD, the optimal parameters are set as αˆ = 952 and K = 4. As shown in Figure 19c, the rotating frequency 43.33 Hz and its harmonics 86.67 Hz can be found. Besides, the oil-whirl frequency 21.5 Hz and its related frequency 65 Hz (approximately equal to the sum of the oil-whirl frequency 21.5 Hz and the rotating frequency 43.33 Hz) is also extracted clearly, which means that the rotor shows an oil-whirl fault. It can be seen from Figure 19d that only three frequency components (i.e., 21.5 Hz, 43.33 Hz, and 86.67 Hz) are visible due to the overlap of frequency components, their related frequency 65 Hz cannot be found, and the end effect is obvious compared with the proposed method. In Figure 19e,f, the oil-whip frequency 21.5 Hz and rotating frequency 43.33 Hz can be seen, Five approaches (i.e., IAVMD, OVMD, LMD, EMD, and WT) are applied to decompose the rub-impact fault signal, respectively. Figure 18c-g respectively depicts the TFR resulting from the five methods. The optimal inside parameters in IAVMD are selected asα = 1177 and K = 4 through using the GWO algorithm. The optimal parameters of OVMD are determined asα = 536 and K = 4 through using the PSO algorithm. As seen in Figure 18c, the fundamental frequency of 29.5 Hz and its various harmonics (59 Hz, 147.5 Hz, and 265.5 Hz) can be identified. In addition, the high frequency is relatively ambiguous due to a resolution problem, but the spectral line is very clear in the low frequency. This indicates that our presented IAVMD approach can detect early characteristics of rotor rub-impact faults.

Case 2: Rotor Oil-Whirl Fault Detection
As shown in Figure 18d, we can easily find the frequency components of 29.5 Hz and its second-harmonic (59 Hz), but the results have overlap and there is a serious end effect (see the parts marked with a red line). In Figure 18e,f, the instantaneous frequencies of 29.5 Hz can be found, but their time-frequency trajectory is obscure. It is obvious from Figure 18g that the TFR obtained by WT can dimly extract the rotating frequency of 29.5 Hz, but the frequency points are very scattered and their time-frequency resolution is clearly worse. The contrastive analysis demonstrated that the efficacy of IAVMD in rotor rub-impact fault detection is superior to that of the other four approaches (i.e., OVMD, LMD, EMD, and WT). Figure 19a,b shows the waveform and FFT spectrum of the rotor oil-whirl signal, respectively. Figure 19a shows that the rotor oil-whirl signal is represented by a sum of sine waves, and that in the FFT spectrum there is a prominent amplitude at the rotating frequency f 2 of 43.33 Hz. Besides, the amplitude of the oil-whirl frequency of 21.5 Hz (equal to about 0.5 f 2 ) is greater than that of the rotating frequency, but the amplitude of the high frequency is not obvious.

Case 2: Rotor Oil-Whirl Fault Detection
2 amplitude of the oil-whirl frequency of 21.5 Hz (equal to about 0.5 2 f ) is greater than that of the rotating frequency, but the amplitude of the high frequency is not obvious. The oil-whirl signal is analyzed by five methods (i.e., IAVMD, OVMD, LMD, EMD, and WT), respectively. Figure 19c-g provides the TFR derived from the five methods, respectively. In IAVMD, the optimal parameters obtained by the GWO algorithm are αˆ = 1919 and K = 4, respectively. In OVMD, the optimal parameters are set as αˆ = 952 and K = 4. As shown in Figure 19c, the rotating frequency 43.33 Hz and its harmonics 86.67 Hz can be found. Besides, the oil-whirl frequency 21.5 Hz and its related frequency 65 Hz (approximately equal to the sum of the oil-whirl frequency 21.5 Hz and the rotating frequency 43.33 Hz) is also extracted clearly, which means that the rotor shows an oil-whirl fault. It can be seen from Figure 19d that only three frequency components (i.e., 21.5 Hz, 43.33 Hz, and 86.67 Hz) are visible due to the overlap of frequency components, their related frequency 65 Hz cannot be found, and the end effect is obvious compared with the proposed method. In Figure 19e,f, the oil-whip frequency 21.5 Hz and rotating frequency 43.33 Hz can be seen, but they show the end effect and mode mixing problem. Furthermore, their time-frequency trajectory is not as good as those in Figure 19c. Likewise, as seen in Figure 19g, the fault feature information of the rotor oil-whirl is not obvious, and the points on the frequency components are discrete and remarkably decentralized. The contrastive analysis further confirms that the performance of the proposed IAVMD in detecting rotor oil-whirl is better than that of the other four methods (i.e., OVMD, LMD, EMD, and WT).  The oil-whirl signal is analyzed by five methods (i.e., IAVMD, OVMD, LMD, EMD, and WT), respectively. Figure 19c-g provides the TFR derived from the five methods, respectively. In IAVMD, the optimal parameters obtained by the GWO algorithm areα = 1919 and K = 4, respectively. In OVMD, the optimal parameters are set asα = 952 and K = 4. As shown in Figure 19c, the rotating frequency 43.33 Hz and its harmonics 86.67 Hz can be found. Besides, the oil-whirl frequency 21.5 Hz and its related frequency 65 Hz (approximately equal to the sum of the oil-whirl frequency 21.5 Hz and the rotating frequency 43.33 Hz) is also extracted clearly, which means that the rotor shows an oil-whirl fault. It can be seen from Figure 19d that only three frequency components (i.e., 21.5 Hz, 43.33 Hz, and 86.67 Hz) are visible due to the overlap of frequency components, their related frequency 65 Hz cannot be found, and the end effect is obvious compared with the proposed method. In Figure 19e,f, the oil-whip frequency 21.5 Hz and rotating frequency 43.33 Hz can be seen, but they show the end effect and mode mixing problem. Furthermore, their time-frequency trajectory is not as good as those in Figure 19c. Likewise, as seen in Figure 19g, the fault feature information of the rotor oil-whirl is not obvious, and the points on the frequency components are discrete and remarkably decentralized. The contrastive analysis further confirms that the performance of the proposed IAVMD in detecting rotor oil-whirl is better than that of the other four methods (i.e., OVMD, LMD, EMD, and WT). Figure 20a,b shows the waveform and FFT spectrum of the rotor oil-whip signal with a length of 1024, respectively. As depicted in Figure 20b, the oil-whip frequency f w of 30 Hz and rotating frequency f w of 75 Hz can be clearly seen in the FFT spectrum, but the amplitudes at a high frequency are weak. Besides, the FFT spectrum is short of time-varying frequency contents.

Case 3: Rotor Oil-Whip Fault Detection
Analogously, five methods (i.e., IAVMD, OVMD, LMD, EMD, and WT) are utilized to deal with the oil-whip signal. Figure 20c-g plots, respectively, the TFR obtained by the five methods. The inside parameters (α = 1987 and K = 4) of IAVMD are determined by using GWO. Two key parameters of OVMD are determined asα = 1175 and K = 4 by using PSO. As shown in Figure 20c, the oil-whip frequency (i.e., f w = 30 Hz) and rotating frequency (i.e., f 3 = 75 Hz) are very obvious. Moreover, the associated frequency components 120 Hz (i.e., f 3 + 1.5 f w ) and 150 Hz (i.e., 2 f 3 ) are also extracted, and the aggregation along the time-frequency trajectory is very explicit. This indicates that the oil-whip fault can be detected accurately by using the proposed method. As shown in Figure 20d, only three frequencies (i.e., f w = 30 Hz, f 3 = 75 Hz, and the relevant frequency 120 Hz) can be found, but the harmonic frequency of 150 Hz is almost invisible due to the discrete distribution of high frequency points. Besides, there is an end effect problem in Figure 20d. Seen from Figure 20e,f, the TFR derived from LMD and EMD can also extract the oil-whip frequency f w = 30 Hz, but the mode mixing phenomenon is very serious, and the time-frequency trajectory at higher frequency is discrete and blurry. The oil-whip frequency f w = 30 Hz can be seen in Figure 20g, but the resolution of the WT-based TFR is worse than that in our approach. Hence, according to the above comparison results, compared with the other methods, our approach is more effective for detecting rotor oil-whip. In short, we can know from Figures 18-20 that the proposed method is once again validated as being valuable in analyzing multi-component signals. 20g, but the resolution of the WT-based TFR is worse than that in our approach. Hence, according to the above comparison results, compared with the other methods, our approach is more effective for detecting rotor oil-whip. In short, we can know from Figures 18-20 that the proposed method is once again validated as being valuable in analyzing multi-component signals.

Result and Discussion
As analyzed above, the presented approach is applicable to processing multi-component signals and is more effective than the existing methods (i.e., OVMD, LMD, EMD, and WT) in extracting time-frequency contents and fault features. The advantage of merging several procedures (i.e., waveform matching extension, GWO, VMD, and TEO) is that it makes the resulting algorithm more adaptive to the practical situation and that it can directly extract time-varying characteristics. Besides, the presented approach can improve the precision of signal decomposition. However, the above results are geared to qualitative analysis, which indicates that a quantitative analysis should also be conducted to further show the performance of the proposed method. Here, three evaluation indexes (i.e., the energy error θ , orthogonal index OI , and correlation coefficient shown in

Result and Discussion
As analyzed above, the presented approach is applicable to processing multi-component signals and is more effective than the existing methods (i.e., OVMD, LMD, EMD, and WT) in extracting time-frequency contents and fault features. The advantage of merging several procedures (i.e., waveform matching extension, GWO, VMD, and TEO) is that it makes the resulting algorithm more adaptive to the practical situation and that it can directly extract time-varying characteristics. Besides, the presented approach can improve the precision of signal decomposition. However, the above results are geared to qualitative analysis, which indicates that a quantitative analysis should also be conducted to further show the performance of the proposed method. Here, three evaluation indexes (i.e., the energy error θ, orthogonal index OI, and correlation coefficient shown in Equation (11)) are adopted to compare, on a quantificational level, the decomposability of different methods.
The energy error θ of one signal x(t) can be given by [58]: where RMS x represents the root mean square of x(t), RMS i represents the root mean square of the i-th modes, and n + 1 represents the mode number. According to Ref. [58], the smaller θ value represents a higher decomposition accuracy and a lesser end effect.
The orthogonal index of one signal x(t) can be given by [59]: where N C represents the mode number, N means the data length of modes, C ik (t) and C jk (t) respectively denote the i-th and j-th modes at the sifting step k, and x k and r k denote the raw signal and residual term, respectively. Similarly, according to Ref. [59], the smaller the value of the orthogonal index is, the better the decomposition performance is. Furthermore, it should be explained that we evaluate their performance through calculating the correlation coefficient between the raw signal with the signals reconstructed by different methods. Tables 2-4 list the energy error θ, orthogonal index OI, and correlation coefficient of the decomposition results obtained by different approaches for the simulation data and experimental signal, respectively. It can be known from Tables 2 and 3 that two evaluation indexes (θ and OI) of our proposed IAVMD approach are lower than those in four similar approaches (i.e., OVMD, LMD, EMD, and WT), which implies that our approach is equipped with satisfactory a decomposition performance in analyzing complicated compound signals compared with other similar methods (i.e., OVMD, LMD, EMD, and WT). However, it should be pointed out that the computational load of our approach is greater than those in other methods (i.e., OVMD, LMD, EMD, and WT) because of the waveform matching and parameter optimization procedure. Besides, OVMD also has a high computation time including the time cost of two stages (i.e., parameter optimization and signal decomposition). Hence, with the help of the high-end computing platform and computer technology, the improvement of the calculation efficiency of the proposed approach is regarded as our following work and research direction. Besides, as shown in Table 4, the correlation coefficient of our presented approach is greater than those in four similar approaches, which indicates that our presented approach is equipped with a lesser reconstruction error and a better decomposition performance compared with some typical techniques mentioned in the study.
To compare the decomposability of different methods more intuitively and quantify the time/memory costs of different methods more precisely, we add the percentage differences of the IAVMD method with respect to the others (OVMD, LMD, EMD, and WT). Tables 5 and 6 list the percentage differences between the IAVMD and other approaches in the simulation and experiment signal, respectively. It can be seen from Tables 5 and 6 that the index value (energy error θ and orthogonal index OI) of the IAVMD method is reduced compared with the other approaches (OVMD, LMD, EMD, and WT). The smallest percentage difference of IAVMD is about 0.1%, but the largest percentage difference can reach 9%. This indicates that the proposed IAVMD method has a certain degree of advantage in signal decomposition. However, at present, the significance of quantification of this small degree is still undetermined, which is the focus of our future research. Table 7 lists the memory cost and average time of different methods. Note that the memory cost of different methods is measured based on the memory usage of the MATLAB file of different methods. As shown in Table 7, the memory cost of IAVMD is the highest (above 8 Kbyte), the second biggest memory cost is the OVMD method, and WT has the smallest memory cost (i.e., the average cost is 1.39 Kbyte).
In addition, the average time costs of our method and OVMD are high (above 25 s), whereas the other three methods (LMD, EMD, and WT)) have low time costs which are under 4 s.   In a word, by observing the frequency components of the time-frequency representations, we can find that the IAVMD method can excavate a better and more accurate frequency content compared with other approaches (OVMD, LMD, EMD, and WT). Besides, through the above quantification, we can draw the conclusion that the IAVMD method can improve the decomposition precision of the signal to some extent, but that it also brings about an additional computational cost due to an increased preparation stage. If one wants to extend the proposed IAVMD method to the practical engineering application, its time cost needs to be decreased, which is also the focus of our future study.

Conclusions
In this article, to solve the shortcoming of the subjective parameter selection and end effect phenomenon existing in the original VMD, a novel time-frequency analysis tool (i.e., IAVMD) consisting of waveform matching extension, GWO-based parameter optimization, and TEO demodulation is presented, which can automatically determine the inside parameters (α, K) of VMD and separate a non-stationary compound signal into a superposition of modes. Several simulated data and experimental examples are provided to validate the effectiveness of our approach. Studies confirm that our approach can effectively identify fault features related to rotor early rub-impact, oil-whirl, and oil-whip. Besides, our proposed method outperforms four available approaches (i.e., OVMD, LMD, EMD, and WT) in analyzing the non-stationary compound signal with multi-frequency contents. The novelties and attractions of this paper are as follows: (1) The concept of waveform matching extension is introduced to alleviate the end effect problem in the traditional VMD method.
(2) A novel time-frequency analysis algorithm termed as IAVMD is presented, which can avoid the effective manual parameter selection and improve diagnostic accuracy.
(3) The efficacy of our method is demonstrated through the simulated and experimental data.
Our future direction will focus on extending the proposed method to analyze multi-component signals in practical applications. Moreover, in our future work, fault diagnosis analyses under variable speeds will also be investigated by using the presented approach. Finally, in the follow-up work, the availability of our presented approach will be discussed in detail when strong noise and more interference are appended to the analyzed vibration signal.