A Novel Empirical Variational Mode Decomposition for Early Fault Feature Extraction

Early fault features of large-scale and low-speed mechanical equipment with heavy duty are weak and exhibit strong non-stationary characteristics. The adaptive extraction and identification of highly relevant important features from such signals has attracted significant attention. In this study, a novel empirical variational mode decomposition and exact Teager energy operator are proposed to explore valuable information. To highlight the fault impact signal representation, we use the exact energy operator to enhance the weak-impact components in the early fault signal. The proposed binary mechanism effectively distinguishes irrelevant features based on the adaptive decomposition parameter construction strategy. Therefore, interference features are easily removed from similar mixed signals, and the independent mode features are determined. The experimental results of the simulation and collected data are compared with those obtained with existing signal decomposition methods, and the superiority of the proposed method, owing to its better modal distinction and less time consumption, is verified.


I. INTRODUCTION
Large-scale machinery under low-speed conditions (e.g., blast furnace rotating distributor, continuous cast steel ladle return turntable, converter rotating support mechanism, and other equipment) is a special type of heavy-duty rotating machinery. It is characterized by a complex transmission structure, large bearing capacity, low working speed, and typical intermittent operation of the transmission mechanism. Bearings and gears are crucial components of transmission machinery [1]. The rolling bearing system is constantly in the transitional process of starting and braking, and its dynamic characteristics exhibit nonlinearity in the transitional process, making conventional vibration monitoring difficult to conduct and diagnosis technology difficult to operate. An urgent problem for metallurgical enterprises is monitoring The associate editor coordinating the review of this manuscript and approving it for publication was Gerard-Andre Capolino. the working state of low-speed heavy-duty machinery, obtaining the early micro-fault information of bearing, and predicting the safe working time of the equipment [2]. The vibration signals of the early failures of such drive bearing are extremely sparse and contain only a few effective signal components, such as frequency division, frequency doubling, micro-impact, and even a single micro-impact component. In particular, a single micro-impact component describes the early fault characteristics of equipment more effectively. In addition, the equipment runs at a low speed (less than 100 rpm) under heavy load for a long time, and the microimpact signal of early failure is severely weak. The existing signal processing technology makes it difficult to extract fault features directly.
In the past decade, non-stationary signal processing methods for engineering applications have attracted considerable research attention. The short-time Fourier transform (STFT) [3], [4], [5] is a commonly used signal time-frequency (TF) domain analysis method. STFT has a variety of window functions, and it is difficult to select the correct window. Moreover, its time-frequency resolution is fixed. Khalil et al. [6], [7] proposed a fast Fourier transform (FFT) to obtain fault frequency signatures. A wavelet transform (WT) [8], [9] overcomes the shortcoming of the invariable length of the STFT transform window and has the characteristic of multiresolution analysis. Zhang et al. [10] proposed an improved empirical wavelet transform to analyze a bearing signal with a complicated spectrum. Yuan and Tian [11] proposed an intelligent industrial process monitoring and fault diagnosis method based on the discrete WT and deep learning. Chen et al. [12] proposed a WT to solve the problem of the early fault diagnosis of rolling bearings under strong background noise. Wang et al. [13] proposed a sparsityguided empirical wavelet transform to automatically establish the Fourier segments required in empirical wavelet transform for the fault diagnosis of rolling element bearings. Its signal resolution performance is highly dependent on the wavelet basis function and scale parameters; however, it is difficult to match properly. The Wigner-Vile distribution is a Fourier transform of an instantaneous signal [14], [15] and has many outstanding properties. Sur et al. [16] proposed a smoothed pseudo-Wigner-Ville distribution method for extracting the instantaneous frequency of a highly nonstationary polynomial frequency-modulated fringe signal. Khare et al. [17] proposed a smoothed pseudo-Wigner-Ville distribution technique to obtain the scalogram, spectrogram, and time-frequency representation plots of the EEG. A tunable Q-factor wavelet transform was proposed to diagnose bearing faults in high-speed rails [18]. However, it cannot ensure non-negativity, particularly for multicomponent signals, and severe cross-term interference will occur. Although many methods have been proposed to inhibit cross-terms, eliminating them completely is difficult. Ensemble empirical mode decomposition (EEMD) is adaptive and suitable for analyzing nonlinear and non-stationary signals [19], [20], [21], [22]. This method has been successfully applied for the analysis of mechanical vibration signals. Chen et al. [23] proposed complementary ensemble empirical mode decomposition (CEEMD) to diagnose the fault of a gearbox. CEEMD and singular value decomposition (SVD) have been proposed for bearing fault diagnosis [24]. Jiang et al. [25] proposed the CEEMD with adaptive noise (CEEMDAN)-based permutation entropy as a sensitive feature for spiral bevel gear fault identification. An improved CEEMDAD was proposed to diagnose the fault of a high-speed train bogie [26]. CEEM-DAN decomposition was performed on the denoised signal to obtain multiple groups of intrinsic mode functions (IMF) [27]. In [28] and [29], CEEMD and CEEMDAN were used to extract the fault features of bearings. The S-transform (ST) [30], [31], [32], [33] combines the advantages of the short-time Fourier transform and WT and has good timefrequency analysis ability. However, it still exhibits shortcomings in practical applications, such as its inability to gather time-frequency energy and the roughness of its spectrum. Although many attempts have been made to improve these problems, its essence is still a special case of the Gaussian window adopted by the STFT. Higher-order statistical analysis (HOSA) can effectively suppress the Gaussian noise components in the signal [34], [35], [36], [37], [38]. However, it is ineffective for non-Gaussian noise. In addition, it is difficult to apply this method to practical engineering because it requires numerous calculations. Although many of these improvements have significantly reduced the computational complexity of HOSA, it remains unsuitable for practical use. Local mean decomposition [39], [40], [41], [42], [43], [44] is typically used to extract the feature information in different frequency bands of the original signal. However, it has some disadvantages such as end effects and mode mixing. Although many researchers have made improvements to the LMD based on their experience in dealing with EMD problems, a slight effect has been achieved in practice. Intrinsic time-scale decomposition (ITD) [45], [46], [47], [48] has some advantages over EMD and LMD. The definition of the baseline in this method is based on a linear transformation, which results in the distortion of the signal waveform from the second component. Several scholars have attempted to solve this problem. However, it is difficult to eliminate. VMD is a new adaptive signal processing method proposed in 2014 [49]. Its excellent characteristics enable it to solve inherent problems in EEMD, LMD, and ITD, and it is widely used in mechanical fault diagnosis [50], [51], [52], [53]. However, it has an inherent defect, which is that the choice of its key parameters affects its performance. Consequently, VMD has undergone many improvements. Nazari and Sakhaei [54] proposed a consecutive VMD that does not need to consider the initial value of K ; however, it requires an artificially set convergence value. Although the K value of the VMD has a significant impact on its performance, the impact of other parameters on VMD performance is equally important. Chen et al. [55] proposed a self-correcting VMD method, where K and α are automatically updated by the energy ratio and orthogonality of the IMFs in the frequency domain, respectively. This method requires a stop condition provided by the signal-to-noise level. However, the noise level of the actual vibration signal is unavailable, and it must be clarified whether the orthogonality between the modes is linear or nonlinear. Therefore, the practical engineering applications of this method should be verified further. Xu and Hu [56] proposed a grey Wolf optimization algorithm for VMD to search for the optimal combination value of [K, α] of VMD. It uses the minimum mutual information between adjacent IMF components as a fitness function; however, this fitness function is not rigorous. A grasshopper optimization algorithm was used to improve the VMD [57]. It assumes the weighted form of the kurtosis of the IMF component as a fitness function. However, this single fitness function is still not sufficiently rigorous. Liu et al. [58] selected VMD parameters based on kurtosis. However, the search method relies on personal experience and intuition of the search mechanism and lacks a mathematical basis. An identification method [59] was proposed to set the critical parameters K and α of the VMD according to the number of spectral peaks and minimum frequency of the signal FFT. However, FFT cannot effectively handle non-stationary and nonlinear signals. Shi and Yang [60] used the central frequency observation method to set the K and τ values based on residual indicators. However, the VMD's performance is closely related to its key parameters. It is unreasonable to set the K and τ values separately using independent methods and indices, without considering the joint action of the parameters. In addition, the influence of the key parameter α was omitted. Furthermore, the actual vibration signal is complex and variable, and this observation method is unsuitable for engineering applications. The VMD optimization method based on particle swarm optimization (PSO) algorithm has been applied successively in mechanical bearing fault diagnosis [61], [62]. By searching for the optimal values of K and α, the performance of the traditional VMD was improved, and a higher decomposition accuracy was achieved. However, the classical PSO algorithm has some inherent disadvantages such as the ability to easily fall into a local optimal solution. A new variable-dimensional composite chaotic PIO method [63] and a whale optimization algorithm [64] were used to search for the optimal combination value [K , α] of the VMD. The proposed framework exhibits excellent performance; however, the complexity of the algorithm and the feasibility of practical applications should also be verified. In [65], a fusion indicator was designed as an adaptation function to search for the optimal combination value of [K , α] for the VMD. However, it does not consider the impact of the other key parameters of the VMD on its performance. The methods mentioned above select the parameters [K , α] of the VMD with certain results. However, challenges persist, such as performing a feasibility analysis of intelligent optimization algorithms. In addition, the problem of the shared influence of other key parameters of the VMD has not been sufficiently considered (i.e. τ ), and neither has the influence of [K , α]. On the other hand, the computational complexity of parameter optimization using intelligent optimization algorithms is also a non-trivial problem, especially in more union parameters selection. An effective approach is to automatically adjust the parameters according to the VMD decomposition mechanism and depending on the characteristics of the signal itself, and the common influence of joint parameters are considered at the same time. Inspired by the above discussion, this study proposes a novel empirical VMD for binary-tree-based multiparameter selection.
The early fault characteristic signals of large machinery with low speeds and heavy loads are considerably sparse and weak. It is difficult to extract their fault features directly using modern signal processing methods. There is little research on the topic. The Teager energy operator (TEO) is an effective tool for enhancing the transient properties of signals and detecting shock components with the advantage of low computational complexity. With an optimal VMD algorithm, TEO was applied to extract the early fault features of large low-speed heavy-duty mechanical equipment and achieve the desired results [63]. Moreover, the TEO is an energy operator that will change the spectral distribution of the original signal. Fusion of the signal mass was proposed to partly improve the TEO method in [66]. However, this does not change the nature of the TEO as an energy operator. In this study, a TEO based on the enhanced factor γ , termed the exact TEO (ETEO), is proposed. Its output is T Energy = γ ·A 2 · 2 , which is the exact form of the TEO. Moreover, the TEO is based on the energy forms of x 2 i − x i+1 · x i−1 and ∀i ∈ Z , that is, the output value is not negative. However, the vibration signals generated by mechanical faults are positive and negative. The classical TEO can enhance the instantaneous amplitude of the micro-impact components in fault signals and change the positive and negative characteristics, and frequency distribution of the signals. Therefore, the output form of TEO can be obtained as T Energy = (±)γ · A 2 · 2 by discriminating between positive and negative signals. This enhances the amplitude of the impact signal more effectively and retains the frequency distribution of the signal x i . Therefore, in this study, we propose an early fault feature extraction method that combines the advantages of empirical variational mode decomposition (EVMD) and ETEO. Finally, the effectiveness of the proposed method is verified through an experimental analysis of the simulation signals and the measured bearing fault signals. The main innovations can be summarized in the following aspects: 1) The binary-tree mechanism and the mutual information variant are introduced to identify the key information of the signal and determine the number of valid modes; 2) In the signal hierarchical decomposition model, the penalty parameter α of VMD algorithm is not preset as a constant, but is dynamically set by the information entropy of the signal property; 3) Considering the combined effect of the update step parameter τ and other parameters on the VMD decomposition performance, the signal standard deviation is introduced to update τ ; and 4) As an effective improved model of the original TEO, the proposed ETEO is used for impact feature extraction promotion. The algorithm is fully adaptive and has low complexity because it fully considers the joint effect of the key parameters [K , α, τ ] of the VMD.
The remainder of this paper is organized as follows. Section II introduces the basic methods and principles used in this study. Section III presents the proposed methods and principles in detail. Sections IV and V verify the effectiveness and superiority of the proposed method through simulated signal analysis and real measurement signals. Finally, conclusions are presented in Section VI.

II. THEORETICAL FOUNDATION A. TEAGER ENERGY OPERATOR
For a continuous discrete time series x(n), the mathematical definition of the TEO of the signal is TEO is a type of nonlinear difference operator that has outstanding adaptability and can recognize the transient components of signals. Moreover, the algorithm is simpler than the Hilbert demodulation and can increase the impact components of rolling bearing faults [67].

B. VARIATIONAL MODE DECOMPOSITION(VMD)
VMD is a non-recursive technique that decomposes an actual signal x(t) into K -independent modes. Its operating principle combines the Hilbert transform and Wiener filter to obtain a set of characteristic components u k (k = 1, 2, . . . , K ) with a limited bandwidth [53]. To solve the variational model, the constrained problem must be transformed into a nonconstrained problem, and an augmented Lagrange expression is introduced: where α is the penalty factor and λ is the Lagrange factor.

C. INFORMATION ENTROPY
Information entropy is a measure of the uncertainty in the amount of information and indicates the average uncertainty of a signal [26]. In a signal, the probability of occurrence of specific information reflects the amount of information that it contains. The complexity of the information contained in a signal can be measured using an information entropy index. Currently, information entropy is widely used in mechanical fault diagnosis and medical diagnosis, and good research results have been obtained. This is defined by the following equation: The random variable was denoted by X . The values of the random variables are {x 1 , x 2 , · · · , x n }, where p(x i ) denotes the probability that the event x i occurs with p(x i ) = 1. We define the information on the event x i as the negative logarithm of its probability of occurrence, denoted as I (x i ).
H(X) is the average information of random variable X; that is, the information entropy of X where P(x i ) denotes the prior probability of event x i and P(x i ) = 1. The base of the log in the formula is related to the information entropy unit. It is widely used with a base of 2, which has a unit of bits, and when using a base of e with a unit of Nat.

D. STANDARD DEVIATION
The variance represents the extent to which the random signal x(t) deviates from its meanx and is the dynamic component of the depicted signal, whose discrete data expression is The standard deviation is a variation of the variance (positive square root), denoted by σ x . This indicates the degree to which the signal deviates from the steady state.

E. LEAST SQUARES MUTUAL INFORMATION
The mutual information (MI) is more accurate than the correlation coefficient method [33]. Most signals, such as typical vibration and shock signals, satisfy the zero-mean property. According to the principle of uncorrelated and orthogonality equivalence between zero-mean random signals, MI can measure the level of similarity between the IMF components and residuals obtained during VMD decomposition, i.e., it can measure the occurrence and degree of modal mixing. MI is defined as follows: As can be observed from the above equation, the logarithmic function in MI reacts significantly to outliers and affects the accuracy of the estimates. Therefore, to overcome this problem, the square loss MI is used in this study to replace the logarithmic function to reduce the interference of outliers and obtain more accurate estimates of the MI. The substitution equation is defined as follows: To avoid calculating the individual probabilities of joint probability p(x,y), edge probability p(x), and edge probability p(y), the least-squares estimation method is introduced to calculate the squared loss MI; instead, the density ratio function that combines them is learned directly and substituted into the following equation, which is equivalent to the squared loss MI: The least squares MI estimator (LSMI) can be obtained as follows: where the regularization parameter λ and basis function ψ contain parameters that can be determined by the optimization algorithm associated with rule J.

III. PROPOSED METHODOLOGY A. EXACT TEAGER ENERGY OPERATOR
In this study, the enhanced factor γ is introduced, and (1) is modified as follows: where A is the amplitude, = ω 2π is the frequency, and ζ is the attenuation coefficient. Let γ = sin 2 2 · e −2ζ n and then modify the above formula to The output is retained as the square of the instantaneous amplitude of the signal multiplied by the instantaneous frequency squared. Here, the above equation is modified to give where T is the enhanced TEO, which is the exact output of T . Moreover, γ = sin 2 2 · e −2ζ n ∈ (0, 1) in this equation. Evidently, the output of the energy operator T is significantly larger than that of the energy operator T . Therefore, T is the augmented form of T.
Moreover, the sin 2 2 term in the enhanced factor γ is related to variable . The solution of can be obtained according to the above equations using the following procedure: Considering x n = e −ζ n ·A·sin( ·n+ϕ), the instantaneous angular frequency is independent of the attenuation coefficient ζ and the magnitude of the initial amplitude A of the signal for a single-DOF vibration system. Let e −ζ n = 1 and A = 1.
We then obtain Therefore, (13) as shown at the bottom of the page.

B. EMPIRICAL VMD
In this study, we propose an empirical VMD algorithm that can effectively resolve the selection of the VMD parameters K , α, and τ and make the decomposition process of VMD an adaptive process. The detailed steps and pseudocode (described in Algorithm 1) are as follows.
Step 2. Initialize the kernel parameters of the Gaussian radial basis kernel function of the LSMI and set the threshold δ of the LSMI estimator and the threshold η of the reconstruction error ρ.
Step 3. Compute the LSMI estimates for the IMF1 and IMF2 components. If LSMI is equal to zero, there is no identical information between the IMF components. If LSMI = 1, then the information between the IMF components is the same. We then determine whether the LSMI is greater than threshold δ. If it is, decomposition is stopped. Conversely, the reconstruction error, ρ, is computed. If ρ > threshold η, the decomposition is terminated. Otherwise, the IMF components IMF1 and IMF2 (obtained by signal x(t)) are decomposed as new signals x 1 (t) and x 2 (t), respectively, and continue to decompose according to Steps 1 and 2.
Step 4. Decompose x 1 (t) to obtain the IMF11and IMF12 components. The LSMI between the two items is calculated. If the LSMI is greater than threshold δ, the reconstruction error ρ is determined. If ρ is less than the threshold, the VMD key parameter α = round (E × (f s /2) × log(K ))), where f s is the sampling frequency of the signal, round(·) is the rounding function, E is the information entropy, and log(K ) is the log base 10 of K . E can measure the complexity of the signal well and adjust the value of α dynamically, and τ = σ + rand. Decompose x 1 (t) again and obtain two components: IMF11 and IMF12. Thereafter, Steps 2 and 3 are followed. If ρ is larger than the threshold, x 1 (t) is named IMF1 as the first label obtained. Similarly, x 2 (t) is decomposed to obtain the components IMF21 and IMF22. If the LSMI between the two components is larger than threshold δ, the reconstruction error ρ is determined. If ρ > the threshold η, reset the key parameters α = round (E × (f s /2) × log(K ))) and τ = σ + rand, and decompose x 2 (t) again to obtain IMF21 and IMF22. Steps 2 and 3 are executed. If ρ is larger than the threshold, x 2 (t) is selected as the second decomposed IMF component IMF2.
Step 5 can be divided into four steps. In the first case, x 1 (t) is decomposed into IMF11 and IMF12, and x 2 (t) is decomposed into IMF21 and IMF22. Thereafter, the LSMI estimates are calculated for each component. If the LSMI between two IMF components is greater than or equal to threshold δ, they are summed as new IMF components. If the LSMI between the two components is smaller than threshold δ, the corresponding IMF components are used as new IMF components. In the second case, x 1 (t) is decomposed to obtain the components IMF11 and IMF12, whereas x 2 (t) is a single IMF component, and the estimates of the LSMI between each component are calculated. If the LSMI between the two IMF components is greater than the threshold δ, the corresponding IMF components are summed as a new IMF. If the LSMI between two IMF components is less than the threshold δ, the corresponding IMF components are used as new IMF components. In the third case, x 1 (t) is a single IMF component and x 2 (t) is decomposed to obtain components IMF21 and IMF22. We then compute the estimated LSMI for each component. If the LSMI between the two IMF components is greater than the threshold δ, the corresponding IMF components are summed up as a new IMF. Otherwise, the corresponding IMF components are used as new IMF items. In the fourth case, x 1 (t) and x 2 (t) are both mono-component signals and the decomposition stops. Specifically, x 1 (t) and x 2 (t) are the final IMF components.
Step 6. Finally, cases 1-3 in step 5 are judged and handled according to steps 1, 2, and 3. When the LSMI between the two IMF components is greater than threshold δ and the reconstruction error ρ > threshold η, the decomposition is terminated to obtain the final independent components.
After the empirical VMD decomposition process, the complex multicomponent and nonstationary signals can be adaptively decomposed into individual mutually orthogonal monocomponents. The process is presented in Algorithm 1.
Because the impact components of early faults are extremely sparse and weak, EVMD is directly used for decomposition, which causes misjudgment of the signal representing early faults as noise or worthless signals, and fault features cannot be extracted. Therefore, before using EVMD to decompose the vibration signals in this study, ETEO was used to effectively strengthen the impact components in the vibration signals to improve the signal-to-noise ratio (SNR). The EVMD was used to efficiently decompose the enhanced signal and obtain a single IMF component. Subsequently, the average value of the ratio of kurtosis to the information entropy of the IMFs was calculated, and the components of the IMFs that were larger than the average value were selected for signal reconstruction. Finally, the EES of the reconstructed signal was calculated to extract the fault features. A flowchart of the proposed method is shown in Fig. 1.

IV. SIMULATION VALIDATION A. SIMULATION SIGNAL
It is not appropriate to use the actual measured signals to qualitatively describe the characteristics of nonlinear and VOLUME 10, 2022
where y(t) denotes a simulated signal. y 1 (t), y 2 (t), and y 3 (t) represent three sinusoidal signals with different center frequencies, and f 1 = 20 Hz, f 2 = 35 Hz, and f 3 = 210 Hz. y 3 (t) is a high-frequency intermittent signal, and y 4 (t) is a periodic pulse-decaying sine signal with a frequency of 8 Hz.
where m(t) is the simulated exponentially decaying signal and C is the decay factor (C = 750). Furthermore, f r is the simulated rotation frequency (fr = 1 Hz), and f n is the resonance frequency (f n = 3000 Hz), where T denotes the average pulse period. v i represents a small fluctuation, a i is the amplitude of the impulse sequence, and a 0 is the initial amplitude.

B. VALIDATION OF THE DECOMPOSITION METHOD
In this section, the EVMD method is tested using simulation signals. The standard deviation σ and information entropy E of the simulated signal y(t) in Fig. 2 are calculated, and the key parameters of the initialized VMD are fixed as K = 2, α = round(E × (f s /2) × log(K ))), τ = σ, ω = 1, and ε = 1e-7. y(t) is decomposed using EVMD, and the results are illustrated in Figs. 3-16. Fig. 3 shows the first VMD decomposition performed on y(t), and two components, IMF1 and IMF2, are obtained and are illustrated in Fig. 3. The LSMI and reconstruction error between IMF1 and IMF2 are calculated, and the results are listed in Table 1. In the table,  LSMI [   decomposition of y(t) do not satisfy the stopping conditions, and the VMD decomposition of components IMF1 and IMF2 is continued, that is, the second VMD decomposition is initiated. First, the IMF1 component in Fig. 4 is decomposed, and the obtained components are shown in Fig. 5. As shown in the figure, the waveform diagrams of IMF11 and IMF12 are different. To quantify the similarity between the two components, the LSMI between the two components is calculated, and the results are shown in Table 1. LSMI [IMF11,IMF12] = 0.0214 for components IMF11 and IMF12 does not satisfy    the stopping condition, and it is necessary to decompose components IMF11 and IMF12 further. Similarly, the IMF2 component in Fig. 6 is decomposed, and the results are shown in Fig. 7. In the figure, the waveform diagrams of components IMF21 and IMF22 are not the same, and the LSMI between the two components is calculated. The results are listed in Table 1. The LSMI of components IMF21 and LSMI [IMF21,IMF22] = 0.0593 do not satisfy the stopping condition. Therefore, components IMF11 and IMF12 are further decomposed. The decomposition results for the third VMD are shown in Figs. 8-9. IMF11 decomposes into IMF111 and IMF112 components. From the calculated specifications listed in Table 1, the LSMI between IMF111 and IMF112 does not satisfy the stopping condition. The decomposition process should therefore be continued. The decomposition of IMF21 yields IMF211 and IMF212 components, VOLUME 10, 2022     as illustrated in Figs. 13        The results reveal that the correlation between IMF2' and IMF3' is large, that is, IMF2' and IMF3' are merged to obtain the new item IMF2. However, the LMSI values between the other components are small (i.e., IMF1' and IMF2'), indicating that the correlation between them is weak. They are separated as independent components, and finally, four new IMF components, namely IMF1', IMF2', IMF3, ', and IMF4, are obtained, and the results are shown in Fig. 16. To further illustrate the quantitative relationship between the new IMF components and y(t), the LSMI between the new IMF components and y(t) is calculated as follows: LSMI [IMF1,y1] = 0.9886, LSMI [IMF2,y2] = 0.9881, LSMI [IMF3,y3] = 0.9958, and LSMI [IMF4,y4] = 0.8658. From the LSMI values between the IMF component and its corresponding simulation signal, it can be observed that, except for the LSMI between component IMF4 and simulation signal x(t), which is slightly less than 1. In particular, there is a certain difference between the information of components IMF4' and y(t); the other IMF components and their corresponding simulation signals are almost the same. These results indicate that the EVMD algorithm is viable and effective.
The above experiments confirmed the effectiveness of EVMD. To further verify its superiority, the traditional methods of VMD, CEEMDAN, LMD, and ITD, which are commonly used for nonlinear and non-stationary signal analyses, are compared. The decomposition results are presented in Figs. 17 and 20. In Fig. 18, the values of the relevant parameters of the VMD are set based on personal experience. Four IMF components are obtained from the results shown in Fig. 18. The IMF1 component is not a mono-component signal, but a superposition of y 1 (t) and y 2 (t) in the original signal, indicating that the component is not completely  decomposed. In addition, by observing the IMF2 component, we learn that it is nearly the same as y 3 (t) in the original signal. Therefore, we can observe IMF3 and IMF4, which are nearly the same component and the same as y 4 (t) in the original signal, and are over-decomposed.  The final decomposition results indicate that the traditional VMD method does not effectively decompose the original signal. Fig. 18 shows the decomposition results of the CEEM-DAN method. Evidently, CEEMDAN adaptively decomposes the signal into 15 IMF components, and the number of components does not match the fraction of components contained in the original signal. Fig. 19 shows the results of the LMD decomposition method, where the LMD adaptively decomposes the signal into four components. However, the waveforms are inconsistent with the actual components of the original signal. Fig. 20 shows the results of the ITD method decomposition. This approach adaptively decomposes the signal into five PR components. Except for the fact that the number of components is inconsistent with the original signal, the waveform of the effective components could not be identified. Therefore, from the experimental results of the above comparison, the conventional signal analysis methods, such as VMD, CEEMDAN, LMD, and ITD, do not effectively extract the original signal x(t); however, the proposed method is superior.
As mentioned previously, some researchers have proposed several intelligent optimization algorithms (i.e., PSO, PIO, and WOA) in [26], [27], [28], [29], and [30]. These methods optimize the performance of VMD, and the premise is to design an appropriate fitness function. Thus, a mixed multicomponent signal is decomposed by parameterized VMD into several single and orthogonal IMF under ideal conditions. These orthogonal IMFs are then superimposed to reconstruct the original signal accurately. However, the following cases may occur in practical applications. For example, the central frequencies of the two signal components are close to each other, resulting in mode mixing and formation of a false IMF. Nevertheless, this component is orthogonal to the other IMF components and can reconstruct the original signal to the greatest extent possible. If the above case occurs, the orthogonality between the IMF components and the precision of the signal reconstruction cannot be effectively used as an evaluation index of the VMD performance nor can it guarantee that the VMD parameters are optimal. To ensure the effectiveness and consistency of the different intelligent algorithms, a fitness function is designed as follows [63]:  Table 2.
The optimal values of [K , α] are determined using intelligent optimization algorithms such as PSO, PIO, WOA, CQPSO, and VDCPIO. In Fig. 21, the number of iterations and convergence values of the different intelligent optimization algorithms vary, and the detailed data are presented in Table 2. In the table, PSO,    to the different intelligent optimization algorithms are shown in Fig. 22. The VMD based on PSO, PIO, and WOA decomposes the original signal x(t) into three IMF components, and the number of IMFs is not consistent with the number of real components. In addition, they cannot effectively decompose the first IMF component synthesized at two similar frequencies. In addition, as shown in Figs. 22(d)-(e), the VMD optimized by CQPSO and VDCPIO decomposes x(t) into four IMFs, and the number of IMFs is the same as that of the original signal. However, they cannot effectively decompose the first IMF synthesized using two similar frequencies.
Finally, the experimental results indicate that when the signal components in the original signal have similar waveforms and frequencies, the VMD easily estimates their center frequencies at the same center frequency; that is, the synthesized signal is incorrectly estimated as one signal component. In addition, utilizing an intelligent optimization algorithm to search for the VMD [K , α] combination value is not exact. This is because the VMD performance is also affected by other parameters. Moreover, only the VMD optimal combined value [K , α] is obtained, which consumes time and is unsuitable for practical engineering applications. Therefore, the results of these experiments indicate that the dichotomous method-based empirical VMD algorithm is advantageous.

C. EVALUATION OF THE ETEO
When the mechanical equipment operates at a low speed (less than 100 RMP) and heavy load for a long time, the micro-impact signal is weak. It is difficult to directly use existing signal-processing technology to extract fault features. The TEO is commonly used to enhance the micro-impact of mechanical vibration signals. However, this changed the directionality of the raw vibration signal. Therefore, the ETEO was proposed to preserve the directionality of the original signal. The signal in (15) was used for analysis to illustrate and verify the effectiveness of the proposed method. Because the actual micro-impact signal of the early fault is extremely weak, the amplitude and frequency of the simulated impact signal shown in Fig. 2 are set to 0.1 V and 8 Hz in this section to simulate the early fault feature of the low-speed and heavy-load mechanical equipment. The mixed signal is shown in Fig. 23(a). Fig. 23(b) is an additive Gaussian white noise of −15 dB added to Fig. 23(a). The envelope spectrum (ES) and envelope energy spectrum (EES) of the mixed signal are calculated in Figs. 24 (a) and 24 (b), respectively. In Figs. 24(a)-24(b), the amplitude of the ES is smaller, and the energy of the EES is more evident. However, their spectral spikes corresponding to the 25, 47, 57, 65, 81, and 105 Hz signals (marked in red) are not the correct multiplications. When the SNR is small, it is difficult to correctly extract the fundamental frequency and its multiplier frequencies information. The mixed signal is then processed using the TEO and ETEO. In Fig. 24, the enhanced signal is an energy signal based on the TEO. In contrast, the enhanced signal based on the ETEO maintains the directionality of the original signal in Fig. 24(b). However, the enhancement effect of the two methods cannot be significantly evaluated. Therefore, the EES of the enhanced signal presented in Fig. 24 is calculated separately, and it can better highlight the spectral lines of the signal, as shown in Fig. 25. First, the ES and EES of the original signal are calculated, respectively, as shown in Fig. 25(a). Comparing the ES with the EES in Fig. 25(a), the amplitude of the ES of the original signal is smaller, and in Figs. 25(a)-25(b), the amplitude of the EES of the original signal is larger. However, their octave spectra that the spectral spikes corresponding to 25,47,57,65,81, and 105 Hz signals (marked in red) are not the correct doubling frequencies. This indicates that it is difficult to correctly extract the fundamental frequency of the impact signal characterizing the fault and its multiplier frequency. Thereafter, the ES of the TEO and ETEO are calculated in Fig. 25(b). Comparing their ES, the TEO can correctly characterize the fundamental and multiple frequencies of the fault frequency in its ETEO, except for the deviation in the fivefold frequency. However, the EES of the ETEO is fully capable of correctly characterizing the fundamental and multiplicity frequencies of the  Fig. 25(d), k0 denotes the kurtosis of the ES of the original signal, and k1 is the kurtosis of the EES. k2 and k3 are the kurtosis of the ES of the TEO and ETEO, respectively, and k4 and k5 represent the kurtosis of the EES of the TEO and ETEO, respectively. The kurtosis value of the EES of the ETEO is the most significant. The above experimental results indicate that the ETEO has a better enhancement effect on the impact component of the signal and maintains the directionality of the signal, as well as shows that the EES has a better characterization of the fault characteristic frequency.

V. EXPERIMENT VERIFICATION
To further illustrate the applicability of the new method in practical applications, a rolling bearing fault diagnosis  platform with low speed and heavy duty is developed. Fig. 26(a) shows the fault test bench, (b) structural diagram of the bearing fault, and (c) and (d) faults in the outer and inner parts of the rolling bearing, respectively. The theoretical parameters of the fault frequencies are listed in Table 3.
Based on the collected data, the ES and EES of the signal are shown in Fig. 27. Fig. 27(a) Figs. 27-29, the EES based on the ETEO is significantly better than that based on the other two strategies. However, the characteristic component is still disturbed by other invalid spectral lines and is not in the entire frequency band. Finally, the EVMD is used to decompose the signal preprocessed by the ETEO, as shown in Fig. 29(a), and the decomposition is shown in Fig. 30(a) Fig. 30(b). By contrast, Fig. 30(c) shows the EES of the same signal. The fundamental frequency (7.1716 Hz), which characterizes the fault characteristic frequency of the inner part and its harmonics (14.3433 and 21.5149 Hz), can be clearly observed, whereas the interference spectral lines are also significantly reduced. These results indicate that the proposed method is more effective and superior than the comparison methods. To further verify the effectiveness and superiority, the extraction of the fault signal of the outer race is conducted,   Finally, EVMD is used to decompose the heightened signal based on ETEO in Fig. 32(a), and the process results are shown in Fig. 33(a). The heightened signal is decomposed     IMF components [IMF1, IMF2] corresponding to the ratios are retained and the new signal is reconstructed, as shown in Fig. 34(b). The envelope energy spectrum of the reconstructed signal is shown in Fig. 34 (c). From the EES in Fig. 34(c), the fundamental frequency (4.8828 Hz), which characterizes the outer fault frequency, as well as its double frequency (9.7656 Hz) and triple frequency (14.6484 Hz), can be clearly observed. The interference items are significantly reduced. Therefore, the above experimental results indicate that the proposed method is more accurate than other methods. Furthermore, it provides a reference diagnosis strategy for early fault diagnosis of large-scale mechanical equipment under low speed and heavy duty.

VI. CONCLUSION
In this study, EVMD based on a binomial tree and ETEO is proposed and combined with EES for the early fault feature detection of large-scale heavy-duty mechanical equipment operating under low speed. The results of simulations and real test experiments proved its effectiveness and superiority. The conclusions are as follows.
1) The TEO can effectively increase the instantaneous amplitude of the impact component in the signal. However, it changes the positive and negative nature of the micro-impact signal amplitude, which in turn changes the frequency distribution of the micro-impact signal. The ETEO is an exact output form of the TEO that maintains the positive and negative of the original signal; that is, it maintains the frequency distribution of the original signal.
2) It is difficult to extract the early fault features of largescale machinery under low-speed and heavy loads. Consequently, VMD tends to identify the early fault characteristics as ''noise,'' and the corresponding shock components cannot be correctly decomposed into the corresponding IMFs. Using the ETEO to enhance the original fault signal before decomposing it with VMD can effectively prevent the VMD from identifying it as ''noise.'' 3) VMD parameters must be set in advance. The choice of relevant parameter values significantly affects the decomposition performance. Although the optimal values of the key parameters of VMD can be selected automatically using an intelligent optimization algorithm, the decomposition time is considerably long, which increases the difficulty of practical application. In addition, most intelligent optimization algorithms search only for the optimal value of the [K , α] combination of the VMD and ignore the effects of the other parameters on the performance of the VMD. The EVMD has a fixed mode parameter of K = 2. The embedded parameters α and τ are dynamically set by the signal to be decomposed to α = round (E × (f s /2) × log(K ))) and τ = standard deviation(σ ), respectively. The experimental results indicate that the method significantly reduces the computing time of VMD and has a higher accuracy of decomposition, which is suitable for the real-time demand of practical engineering. 4) Large-scale and low-speed mechanical equipment with heavy duty are characterized by complex transmission structures, large load capacities, and low operating speeds. They work as a typical intermittent transmission mechanism, and the transmission mechanism is constantly in transition between starting and braking.
In future research, we will focus on verifying the effectiveness of our proposed method. In particular, it can be used for fault diagnosis of rotating machinery under noise interference. The decomposition accuracy can be further improved by studying the optimization of the theoretical model for parameter selection.
BO XU received the M.Sc. degree in control science and engineering and the Ph.D. degree in industrial engineering and engineering management from the Wuhan University of Science and Technology, Wuhan, China, in 2012 and 2019, respectively. He is currently engaged in postdoctoral research at the Engineering Research Center for Metallurgical Automation and Measurement Technology, Wuhan University of Science and Technology. His current research interests include fault diagnosis, prognostic health management, machine learning-based pattern recognition, and intelligent control.
HUIPENG LI received the B.E. degree in information engineering and the M.Sc. degree in control science and engineering from the Wuhan Institute of Technology, Wuhan, China, in 2004 and 2007, respectively. He is currently pursuing the Ph.D. degree with the Department of Information Science and Engineering, Wuhan University of Science and Technology, Wuhan. His current research interests include data driven fault diagnosis, intelligent control systems, and fault-tolerant control of real-time systems. VOLUME 10, 2022