A Fault Feature Extraction Method Based on Improved VMD Multi-Scale Dispersion Entropy and TVD-CYCBD

In modern industry, due to the poor working environment and the complex working conditions of mechanical equipment, the characteristics of the impact signals caused by faults are often submerged in strong background signals and noises. Therefore, it is difficult to effectivelyextract the fault features. In this paper, a fault feature extraction method based on improved VMD multi-scale dispersion entropy and TVD-CYCBD is proposed. First, the marine predator algorithm (MPA) is used to optimize the modal components and penalty factors in VMD. Second, the optimized VMD is used to model and decompose the fault signal, and then the optimal signal components are filtered according to the combined weight index criteria. Third, TVD is used to denoise the optimal signal components. Finally, CYCBD filters the de-noised signal and then envelope demodulation analysis is carried out. Through the simulation signal experiment and the actual fault signal experiment, the results verified that multiple frequency doubling peaks can be seen from the envelope spectrum, and there is little interference near the peak, which shows the good performance of the method.


Introduction
With the continuous acceleration of the industrialization process in the world, modern industry is gradually developing towards the direction of large-scale, complex, highspeed, and automatic production equipment. Especially since the concept of Industry 4.0 was presented, the state data generated by the operation of mechanical equipment has increased [1,2]. As the most basic and essential mechanical components, rolling bearings are widely used in modern aviation, aerospace, navigation, machine tools, etc. Rolling bearing runs under poor working conditions for a long time, so it is easy to have various faults. It may cause excessive energy consumption, cause a bad mechanical equipment vibration, and affect the performance of the relevant instruments on the equipment. In the worst case, it will lead to the shutdown of mechanical equipment, causing substantial economic losses and even endangering the life safety of the relevant staff. Therefore, the fault diagnosis of this component can not only ensure the smooth and healthy operation of mechanical equipment but also help prevent major accidents. Generally, due to the limitation of the working environment of bearing, it is impossible to diagnose directly. In the process of a practical application, the vibration signal extracted by observing the running state of the bearing contains a lot of information. Meanwhile, the requirements for equipment and staff skills are low in the process of signal acquisition. Therefore, vibration signal analysis has become one of the most widelyused fault detection methods. Because the bearingvibration signal under fault conditions shows prominent non-stationary, nonlinear, and weak fault characteristics [3], it is important and difficult for experts and scholars in this field to enhance fault information from the noise background.
(1) The strategy for optimizing VMD by MPA is proposed, and the initialization parameters of the algorithm are optimized. It enables the signal to be decomposed with a high quality and eliminates adverse effects, such as mode aliasing. (2) A combined weight screening criterion that balances the advantages and disadvantages of the two indicators is constructed. On the basis of this, the evaluation of the IMF signal components and the selection of the best signal components are completed. Meanwhile, the signal is subsequently processed by TVD noise reduction to reduce the noise interference. The arrangement of this paper is as follows: Section 2 describes the basic theory of improved VMD, multi-scale dispersion entropy, TVD, and CYCBD algorithms. Section 3 introduces the specific steps and flow chart of the improved VMD multi-scale dispersion entropy and TVD-CYCBD. Section 4 validates the feasibility of the proposed method. Section 5 shows the practical application performance of the proposed method. Section 6 draws a conclusion.

Improved VMD
The main parameters of the VMD algorithminclude the modal components and penalty factors. The modal component ensures the appropriateness and accuracy of the number of decomposition modes. The penalty factor is related to the accuracy of the signal reconstruction. The selection of the two factors plays an essential role in the decomposition effect of the VMD. When the modal component and penalty factor are too large, it is easy to cause modal aliasing. Otherwise, it will cause the loss of useful information. The marine predator algorithm (MPA) proposed in recent years has the advantages of a simple structure, flexible algorithm, easy implementation, and the ability to coordinate exploration and development. Compared with other algorithms, it has a better optimization performance. Therefore, this paper introduced MPA to optimize VMD. Then, the optimization results are entered into VMD for modeling.

VMD
VMD is an adaptive signal processing method proposed by Dragomiretskiy. It decomposes the signal into finite component signals. The component signal is the limited bandwidth of a specific center frequency, which satisfies that the sum of all the component signals is equal to the original signal. The signal is solved iteratively to minimize the sum of the bandwidth. It is defined as follows: where A k (t) is the instantaneous amplitude. φ k (t) is a phase. ω k (t) is the instantaneous phase.
(1) The unilateral spectrum of each modal component is calculated by Hilbert.
(2) The center frequency is estimated and the spectrum of each mode is modulated to the corresponding fundamental frequency band by exponential correction e −jω k t .
The constraint model is built on the basis of the above formula, and then: where {ω k } = {ω 1 , ω 2 , ω 3 , · · ·ω k } is the k-th central frequency of the signal after VMD. t is the time. δ(t) is the Dirac function. j is an imaginary unit. * is a convolution operation. 2 2 is the L2 norm. x(t) is the original signal. To calculate the above model, the penalty factor and Lagrangian multiplication operator are introduced in the above equation. Then, the problem to be solved is transformed into an unconstrained variational problem and the following formula can be obtained:

MPA Method
The main inspiration for MPA [30] comes from the foraging strategy of marine predators, the Lévy motion and Brownian motion, and the optimal encounter strategy of the interaction between predators and prey. The mathematical description of the MPA optimization process is as follows: (1) Initialization phase. To commence the optimization process, the MPA algorithm will randomly initialize the prey location in the search space. The expression formula is as follows: where X max and X min are the upper and lower bounds of the search space.
(2) MPA optimization stage. At the initial stage, under the condition of a high-speed ratio, the predator does not move. When Iter < 1 3 Max_Iter, the mathematical description of the MPA optimization process is as follows: where stepsice is the moving step. Elite i is the elite matrix. Prey i is the prey matrix. P is a constant. P is taken as 0.5. R and is the uniform random vector. R B is the random vector. Iter, Max_Iter is the current iteration number and the maximum iteration number.
In the middle of the iteration, the optimization process is transformed from the exploration to the development when the speed of the predator and prey is the same. Therefore, half will be used for exploration and the other half for development. The mathematical descriptions of the development and exploration are as follows: where R L is a random vector with Lévy distribution. CF = (1 − Iter/Max_Iter) (2·Iter/Max_Iter) . At the end of the iteration, the predator moves based on the way of Lévy walking and its position changes as follows: (3) Fish aggregating devices (FADs) effect and eddy the current effect. This strategy can enable MPA to overcome the premature convergence problem in the optimization process and escape from the local extremum.

Total Variation Denoising
TVD [24] is a nonlinear fast noise reduction method without iteration, which achieves the smooth processing of discrete signals.
Suppose that the signal x with N points is shown as follows: The first-order differential of the matrix (N − 1 × N) is as follows: The second-order differential of the matrix (N − 2 ×N) is as follows: The expression of the l p norm (P ≥ 1) is as follows: In special cases, when p = 1, Equation (16) will become the following form: When p = 2, Equation (16) becomes the following form: Assume that the signal y(n) containing noise is: The noise reduction of this algorithm can be summarized into the following optimization problems: argmin

CYCBD Method
The CYCBD maximizes the second-order cyclostationarity (ICS 2 ) by solving the eigenvalue to extract the required fault signal from the complex source signal. The convolution process is expressed as follows: where s is the original signal. x is the observation signal. h is the inverse filter. The above formula can be expressed as follows: where L is the length of discrete signal s. N is the length of the inverse filter h. Therefore, ICS 2 can be expressed in the form of generalized Rayleigh entropy: where H is the conjugate transpose of the matrix. R XWX is the weighted correlation matrix. R XX is the correlation matrix. The weighting matrix W is as follows.
where W is the weighting matrix. p[|s| 2 ] is the matrix containing the fault characteristic cycle frequency. K is the number of samples. T represents the failure cycle. The cycle frequency is the frequency related to the signal pulse, which is related to mechanical equipment failures such as bearings. Formula (25) is the core of CYCBD. The weighting matrix W can be calculated by the inverse filter h and the observation signal x. Then, using the characteristics of generalized Rayleigh entropy, it can be concluded that the optimal inverse filter h is equivalent to the maximum eigenvector corresponding to the maximum generalized eigenvalue of the weighted correlation matrix R XWX and the correlation matrix R XX in the Formula (24). Then, the calculation result is converted into Formula (21) and the deconvolution signal corresponds to the maximum second-order cyclostationarity, and the fault frequency can be extracted from the observation signal x.

Parameter Optimization of VMD Based on MPA
In the MPA algorithm optimizing VMD, the selection of the fitness function directly affects the optimization effect. Once the bearing is damaged during operation, it will produce a pulse impact. This method can be used to characterize the irregularity of the time series. The greater the entropy of dispersion, the higher the degree of irregularity.
Therefore, the concept of dispersion entropy is used to measure the irregularity of the intrinsic mode function (IMF) time series data obtained by VMD and took the average value of the dispersion entropy of the signal component as the fitness function of MPA. Then, through the global optimization of the mean value of the minimum dispersion entropy, the optimal decomposition level k and the penalty factor αare obtained. The flow chart of MPA-VMD is shown in Figure 1. The specific steps are as follows: (1) Initialize the parameters of the MPA, including the population size, maximum iterations, FADs, and the parameter range to be optimized for VMD. (2) Set the fitness function value, then obtain the initial position of the prey and calculate the fitness value. (3) Iterative optimization is performed from the initial, middle, and late stages of the iteration. Comparing the fitness values in the iteration process can detect whether the optimal fitness value can be found and then update the predator's position. Then, calculate the fitness of the new location and evaluate the impact of the FADs or eddy current effect on the fitness value. Calculate the optimal predator position according to the prey position and behavior and, finally, determine and store the current optimal position. (4) Judge whether the following relationship is true: iter ≥ Max_iter. If not, continue to repeat steps 2 and 3 to find the optimal solution. Otherwise, the calculation is terminated and the set of optimal parameters is the output. (5) The optimal parameters obtained from the MPA algorithm optimization are assigned to VMD to build the algorithm model.

IMF Component Selection Criteria
After VMD, the IMF components obtained may have false components or components with a small correlation with the original signal. If the practical signal components are not filtered, the subsequent noise reduction and fault feature extraction will be somewhat disturbed. Therefore, this paper established a combined weight-effective signal component screening rule that combines kurtosis and cross-correlation coefficients.

IMF Component Selection Criteria
After VMD, the IMF components obtained may have false components or components with a small correlation with the original signal. If the practical signal components are not filtered, the subsequent noise reduction and fault feature extraction will be somewhat disturbed. Therefore, this paper established a combined weight-effective signal component screening rule that combines kurtosis and cross-correlation coefficients.
The cross-correlation coefficient can represent the interdependence between two signal amplitudes. The greater the cross-correlation coefficient between the IMF signal components and the original signal, the stronger the degree of correlation. Kurtosis is an index used to judge the Gaussian performance of vibration signals, which can be used to know whether each IMF signal component carries shock and fault components from the side. Kurtosis reflects the impact of the impact components.
Although kurtosis is sensitive to impact components, proper signals are vulnerable to noise interference. Although the correlation coefficient can measure the correlation between signals, it is easily affected by the total number of samples. By combining the advantages and disadvantages of the above two evaluation indicators, the calculation method based on the combination weight value is as follows: where φ and ϕ are the weights of kurtosis and the correlation coefficient, respectively, and φ + ϕ = 1 is satisfied. In combination with the parameter selection rules in the literature [31], this paper has selected φ = 0.4, ϕ = 0.6 through the experiments.

Implementation Process
The implementation process of the proposed method is as follows. The specific implementation flowchart is shown in Figure 2. advantages and disadvantages of the above two evaluation indicators, the calculation method based on the combination weight value is as follows: where  and  are the weights of kurtosis and the correlation coefficient, respectively, and 1  += is satisfied. In combination with the parameter selection rules in the literature [31], this paper has selected  = 0.4,  = 0.6 through the experiments.

Implementation Process
The implementation process of the proposed method is as follows. The specific implementation flowchart is shown in Figure 2.
(1) Optimize VMD parameters. Initialize the parameters of the MPA algorithm and determine the optimal k and α by optimizing VMD through MPA. (2) The best parameters are substituted into the VMD algorithm model to complete the modeling. Then, the VMD is carried out for the collected signals and several signal components are obtained. The K-C index of each IMF is calculated according to the screening criteria based on the K-C combination weight. Then, we can evaluate the quality of the IMF signal components. Finally, select the IMF signal component with the largest K-C value. (3) Input the filtered IMF signal components into TVD for noise reduction to reduce noise interference in valuable signals. (4) The fault frequency of the original signal is calculated according to the fault frequency formula and the appropriate cycle frequency is set. Then, the signal is filtered by CYCBD.   (1) Optimize VMD parameters. Initialize the parameters of the MPA algorithm and determine the optimal k and α by optimizing VMD through MPA. (2) The best parameters are substituted into the VMD algorithm model to complete the modeling. Then, the VMD is carried out for the collected signals and several signal components are obtained. The K-C index of each IMF is calculated according to the screening criteria based on the K-C combination weight. Then, we can evaluate the quality of the IMF signal components. Finally, select the IMF signal component with the largest K-C value. (3) Input the filtered IMF signal components into TVD for noise reduction to reduce noise interference in valuable signals. (4) The fault frequency of the original signal is calculated according to the fault frequency formula and the appropriate cycle frequency is set. Then, the signal is filtered by CYCBD. (5) The signal filtered by CYCBD is analyzed by Hilbert envelope.

Simulation Verification
To determine whether the proposed method based on improved VMD multi-scale dispersion entropy and TVD-CYCBD can extract the fault pulse impact component, this research constructs the rolling bearing vibration simulation signal: where y 0 is the displacement constant and its value is 5. f n is the carrier frequency and its value is 3000 Hz. ξ is the damping coefficient and its value is 0.1. f s is the sampling frequency and its value is 20 KHz. t indicates the sampling time. T = 0.01 s. The sampling points are N = 4096. The fault frequency is f 0 = 100 Hz. To simulate the bearing fault according to the actual situation, this experiment adds noise to the simulation signal s(t), and its signal-to-noise ratio (SNR) is −5 dB. The waveform of the simulation signal is shown in Figure 3. The time and frequency domain diagram of the simulation signal after adding noise are shown in Figure 3a

Simulation Verification
To determine whether the proposed method based on improved VMD multi-scale dispersion entropy and TVD-CYCBD can extract the fault pulse impact component, this research constructs the rolling bearing vibration simulation signal: where y0 is the displacement constant and its value is 5. fn is the carrier frequency and its value is 3000 Hz.  is the damping coefficient and its value is 0.1. fs is the sampling frequency and its value is 20KHz. t indicates the sampling time. T = 0.01s. The sampling points are N = 4096. The fault frequency is f0 = 100 Hz. To simulate the bearing fault according to the actual situation, this experiment adds noise to the simulation signal s(t), and its signal-to-noise ratio (SNR) is −5dB. The waveform of the simulation signal is shown in Figure 3. The time and frequency domain diagram of the simulation signal after adding noise are shown in Figure 3a Firstly, the MPA optimization algorithm is used to optimize the parameters in VMD. The parameter search range in VMD is as follows: the k∊ [3,15] and the α∊[100,5000]. In the MPA algorithm, the population size is 10 and the FADs are 0.2. Meanwhile, the maximum iteration number is 20. In the process of optimization, the fitness function adopted is the average value of dispersion entropy. After parameter optimization, the fitness curve of the MPA obtained is shown in Figure 4. According to the results shown in Figure  4, when the fitness value is optimal, the corresponding k and α are 3 and 3348, respectively. Therefore, the above parameters are substituted into the VMD for modeling and then VMD is performed on the simulation signal added with noise. To compare and analyze the effects of the MPA-VMD method and other methods, this research also conducted experiments based on the EMD and EEMD methods. The decomposition results obtained based on MPA-VMD, EMD, and EEMD are shown in Figures 5 and 6. Firstly, the MPA optimization algorithm is used to optimize the parameters in VMD. The parameter search range in VMD is as follows: the k ∈ [3,15] and the α ∈ [100, 5000]. In the MPA algorithm, the population size is 10 and the FADs are 0.2. Meanwhile, the maximum iteration number is 20. In the process of optimization, the fitness function adopted is the average value of dispersion entropy. After parameter optimization, the fitness curve of the MPA obtained is shown in Figure 4. According to the results shown in Figure 4, when the fitness value is optimal, the corresponding k and α are 3 and 3348, respectively. Therefore, the above parameters are substituted into the VMD for modeling and then VMD is performed on the simulation signal added with noise. To compare and analyze the effects of the MPA-VMD method and other methods, this research also Compared with Figures 5 and 6, it is shown that the number of signal components obtained based on MPA-VMD is the least, and the modal components have a good separability and weak modal aliasing. More signal components are obtained based on the EMD and EEMD methods: 11 signal components and 1 residual component. Meanwhile, it can be seen that some signal components have severe mode aliasing and many false components are generated. Next, the kurtosis and correlation values of the signal components decomposed by the MPA-VMD, EMD, and EEMD methods will be calculated. Each signal component's K-C combined weight index will be calculated according to the calculation formula of the combined weight proposed in this paper. The calculation results based on MPA-VMD, EMD, and EEMD are shown in Tables 1-3   Compared with Figures 5 and 6, it is shown that the number of signal components obtained based on MPA-VMD is the least, and the modal components have a good separability and weak modal aliasing. More signal components are obtained based on the EMD and EEMD methods: 11 signal components and 1 residual component. Meanwhile, it can be seen that some signal components have severe mode aliasing and many false components are generated. Next, the kurtosis and correlation values of the signal components decomposed by the MPA-VMD, EMD, and EEMD methods will be calculated. Each signal component's K-C combined weight index will be calculated according to the calculation formula of the combined weight proposed in this paper. The calculation results based on MPA-VMD, EMD, and EEMD are shown in Tables 1-3.  Compared with Figures 5 and 6, it is shown that the number of signal components obtained based on MPA-VMD is the least, and the modal components have a good separability and weak modal aliasing. More signal components are obtained based on the EMD and EEMD methods: 11 signal components and 1 residual component. Meanwhile, it can be seen that some signal components have severe mode aliasing and many false components are generated. Next, the kurtosis and correlation values of the signal components decomposed by the MPA-VMD, EMD, and EEMD methods will be calculated. Each signal component's K-C combined weight index will be calculated according to the calculation formula of the combined weight proposed in this paper. The calculation results based on MPA-VMD, EMD, and EEMD are shown in Tables 1-3.  By comparing the K-C combined weight index of each signal component in Table 1, it is shown that the result value of IMF1 obtained based on the MPA-VMD method is the largest. The correlation and kurtosis values of this signal component are the largest compared with other signal components. Therefore, IMF1 is selected as the optimal signal component for the reduction in signal noise. Comparing the K-C combined weight index of each signal component in Table 2 shows that the result value of IMF2 obtained based on the EMD method is the largest. Therefore, IMF2 is selected as the optimal signal component for the reduction in the signal noise. Comparing the K-C combined weight index of each signal component in Table 3 shows that the result value of IMF2 obtained based on the EEMD method is the largest. Therefore, according to the comparison of the computer values, the IMF2 will be used as the optimal signal component for the reduction in the signal noise. Furthermore, it can also be seen from Tables 1-3 that using one of the cross-correlation coefficients and kurtosis values to filter will cause a specific interference in the filtering of the signal components.   By comparing the K-C combined weight index of each signal component in Table 1, it is shown that the result value of IMF1 obtained based on the MPA-VMD method is the largest. The correlation and kurtosis values of this signal component are the largest compared with other signal components. Therefore, IMF1 is selected as the optimal signal component for the reduction in signal noise. Comparing the K-C combined weight index of each signal component in Table 2 shows that the result value of IMF2 obtained based on the EMD method is the largest. Therefore, IMF2 is selected as the optimal signal component for the reduction in the signal noise. Comparing the K-C combined weight index of each signal component in Table 3 shows that the result value of IMF2 obtained based on the EEMD method is the largest. Therefore, according to the comparison of the computer values, the IMF2 will be used as the optimal signal component for the reduction in the signal noise. Furthermore, it can also be seen from Tables 1-3 that using one of the cross-correlation coefficients and kurtosis values to filter will cause a specific interference in the filtering of the signal components. Figures 7-9 are the time-domain waveforms of the signal components filtered based on MPA-VMD, EMD, and EEMD methods after noise reduction. Through analysis, the following conclusions can be drawn: after the TVD noise reduction, the noise interference has been reduced to varying degrees. Several evaluation indexes are introduced in this study to compare the noise reduction effects. They are SNR, RMSE, and MAE. Calculating the above index values shows the results in Table 4.  9 are the time-domain waveforms of the signal components filtered based on MPA-VMD, EMD, and EEMD methods after noise reduction. Through analysis, the following conclusions can be drawn: after the TVD noise reduction, the noise interference has been reduced to varying degrees. Several evaluation indexes are introduced in this study to compare the noise reduction effects. They are SNR, RMSE, and MAE. Calculating the above index values shows the results in Table 4.  The data shown in Table 4 are the results of the three methods. It shows that the SNR based on the MPA-VMD-TVD method is 4.225, the correlation coefficient is 0.792, the RMSE is 0.357, and the MAE is 0.179. Compared with EMD-TVD and EEMD-TVD, the SNR is higher and the RMSE is smaller. Therefore, according to the above results, we can come to the conclusion that the method proposed in this paper has a better noise reduction effect. Next, envelope spectrum analysis is performed on the noise-reduced signals using the above three methods. The envelope spectrum obtained by the hilbert transformation is shown in Figures 10-12. Through the analysis of the above three results, it can be seen that the fault characteristic frequency and its multiple frequency can appear in the envelope spectrum obtained based on the EMD-TVD method. However, the amplitude of the fault frequency is relatively low, which brings some interference to the extraction of the fault characteristics. Meanwhile, although the characteristic frequency of the fault and its double frequency, three-time frequency, four-time frequency, and other components can be seen from Figure 10, due to an insufficient signal noise reduction, the interference of other irrelevant frequencies exists near the peak value of the multiple fault frequencies.
From the envelope spectrum obtained on the basis of MPA-VMD-TVD, the fault charac- The data shown in Table 4 are the results of the three methods. It shows that the SNR based on the MPA-VMD-TVD method is 4.225, the correlation coefficient is 0.792, the RMSE is 0.357, and the MAE is 0.179. Compared with EMD-TVD and EEMD-TVD, the SNR is higher and the RMSE is smaller. Therefore, according to the above results, we can come to the conclusion that the method proposed in this paper has a better noise reduction effect. Next, envelope spectrum analysis is performed on the noise-reduced signals using the above three methods. The envelope spectrum obtained by the hilbert transformation is shown in Figures 10-12. Through the analysis of the above three results, it can be seen that the fault characteristic frequency and its multiple frequency can appear in the envelope spectrum obtained based on the EMD-TVD method. However, the amplitude of the fault frequency is relatively low, which brings some interference to the extraction of the fault characteristics. Meanwhile, although the characteristic frequency of the fault and its double frequency, three-time frequency, four-time frequency, and other components can be seen from Figure 10, due to an insufficient signal noise reduction, the interference of other irrelevant frequencies exists near the peak value of the multiple fault frequencies.
From the envelope spectrum obtained on the basis of MPA-VMD-TVD, the fault charac-  The data shown in Table 4 are the results of the three methods. It shows that the SNR based on the MPA-VMD-TVD method is 4.225, the correlation coefficient is 0.792, the RMSE is 0.357, and the MAE is 0.179. Compared with EMD-TVD and EEMD-TVD, the SNR is higher and the RMSE is smaller. Therefore, according to the above results, we can come to the conclusion that the method proposed in this paper has a better noise reduction effect.
Next, envelope spectrum analysis is performed on the noise-reduced signals using the above three methods. The envelope spectrum obtained by the hilbert transformation is shown in Figures 10-12. Through the analysis of the above three results, it can be seen that the fault characteristic frequency and its multiple frequency can appear in the envelope spectrum obtained based on the EMD-TVD method. However, the amplitude of the fault frequency is relatively low, which brings some interference to the extraction of the fault characteristics. Meanwhile, although the characteristic frequency of the fault and its double frequency, three-time frequency, four-time frequency, and other components can be seen from Figure 10, due to an insufficient signal noise reduction, the interference of other irrelevant frequencies exists near the peak value of the multiple fault frequencies. From the envelope spectrum obtained on the basis of MPA-VMD-TVD, the fault characteristic frequency and its double frequency, three-time frequency, four-time frequency, and other components are clearly displayed. Meanwhile, the outside interference near the fault frequency is significantly reduced. Finally, the above three signals after TVD noise reduction are inputted to the CYCBD filter for further filtering to enhance the signal's periodic impact characteristics. Finally, the above three signals after TVD noise reduction are inputted to the CYCBD filter for further filtering to enhance the signal's periodic impact characteristics. teristic frequency and its double frequency, three-time frequency, four-time frequency, and other components are clearly displayed. Meanwhile, the outside interference near the fault frequency is significantly reduced. Finally, the above three signals after TVD noise reduction are inputted to the CYCBD filter for further filtering to enhance the signal's periodic impact characteristics. Finally, the above three signals after TVD noise reduction are inputted to the CYCBD filter for further filtering to enhance the signal's periodic impact characteristics. In this algorithm, the cyclic frequency set is set to [100,200, . . . ,1000]. Then, the signal filtered by the CYCBD is analyzed in the envelope. The time domain diagrams of the three signals filtered by the CYCBD filter are shown in Figures 13-15, and the generated envelope spectrum is shown in Figures 16-18. It can be seen from Figures 13-15 that after filtering, the periodic impact components in the time domain diagrams based on MPA-VMD-TVD, EMD-TVD, and EEMD-TVD methods have been greatly improved. Through the comparative analysis of Figures 16-18, it is shown that the characteristic frequency of the faultand its double to nine times of the fault frequency can be extracted by using the above three methods. At the same time, in general, the amplitude of the fault frequency and multiple frequencies of the fault frequency obtained based on the proposed method are the highest, and there is no outside interference near the fault frequency. However, there are still different degrees of outside interference near the fault frequency obtained by the other two methods. Therefore, the proposed method achieved better results.
In this algorithm, the cyclic frequency set is set to [100,200, …,1000]. Then, the signal filtered by the CYCBD is analyzed in the envelope. The time domain diagrams of the three signals filtered by the CYCBD filter are shown in Figures 13-15, and the generated envelope spectrum is shown in Figures 16-18. It can be seen from Figures 13-15 that after filtering, the periodic impact components in the time domain diagrams based on MPA-VMD-TVD, EMD-TVD, and EEMD-TVD methods have been greatly improved. Through the comparative analysis of Figures 16-18, it is shown that the characteristic frequency of the faultand its double to nine times of the fault frequency can be extracted by using the above three methods. At the same time, in general, the amplitude of the fault frequency and multiple frequencies of the fault frequency obtained based on the proposed method are the highest, and there is no outside interference near the fault frequency. However, there are still different degrees of outside interference near the fault frequency obtained by the other two methods. Therefore, the proposed method achieved better results. In this algorithm, the cyclic frequency set is set to [100,200, …,1000]. Then, the signal filtered by the CYCBD is analyzed in the envelope. The time domain diagrams of the three signals filtered by the CYCBD filter are shown in Figures 13-15, and the generated envelope spectrum is shown in Figures 16-18. It can be seen from Figures 13-15 that after filtering, the periodic impact components in the time domain diagrams based on MPA-VMD-TVD, EMD-TVD, and EEMD-TVD methods have been greatly improved. Through the comparative analysis of Figures 16-18, it is shown that the characteristic frequency of the faultand its double to nine times of the fault frequency can be extracted by using the above three methods. At the same time, in general, the amplitude of the fault frequency and multiple frequencies of the fault frequency obtained based on the proposed method are the highest, and there is no outside interference near the fault frequency. However, there are still different degrees of outside interference near the fault frequency obtained by the other two methods. Therefore, the proposed method achieved better results. In this algorithm, the cyclic frequency set is set to [100,200, …,1000]. Then, the signal filtered by the CYCBD is analyzed in the envelope. The time domain diagrams of the three signals filtered by the CYCBD filter are shown in Figures 13-15, and the generated envelope spectrum is shown in Figures 16-18. It can be seen from Figures 13-15 that after filtering, the periodic impact components in the time domain diagrams based on MPA-VMD-TVD, EMD-TVD, and EEMD-TVD methods have been greatly improved. Through the comparative analysis of Figures 16-18, it is shown that the characteristic frequency of the faultand its double to nine times of the fault frequency can be extracted by using the above three methods. At the same time, in general, the amplitude of the fault frequency and multiple frequencies of the fault frequency obtained based on the proposed method are the highest, and there is no outside interference near the fault frequency. However, there are still different degrees of outside interference near the fault frequency obtained by the other two methods. Therefore, the proposed method achieved better results.

Experimental Verification
The experimental data used in this study are from the public data set of CWRU [32]. The fault acquisition equipment is shown in Figure 19 [33] and the bearing structure coefficient is shown in Table 5. The motor speed is 1797 rpm. The sampling frequency is 12 KHz. The model used at the drive end is 6205-2RSJEMSKF, the motor speed is 1797 rpm, and the sampling frequency is 12 KHz.

Experimental Verification
The experimental data used in this study are from the public data set of CWRU [32]. The fault acquisition equipment is shown in Figure 19 [33] and the bearing structure coefficient is shown in Table 5. The motor speed is 1797 rpm. The sampling frequency is 12 KHz. The model used at the drive end is 6205-2RSJEMSKF, the motor speed is 1797 rpm, and the sampling frequency is 12 KHz.

Experimental Verification
The experimental data used in this study are from the public data set of CWRU [32]. The fault acquisition equipment is shown in Figure 19 [33] and the bearing structure coefficient is shown in Table 5. The motor speed is 1797 rpm. The sampling frequency is 12 KHz. The model used at the drive end is 6205-2RSJEMSKF, the motor speed is 1797 rpm, and the sampling frequency is 12 KHz.

Experimental Verification
The experimental data used in this study are from the public data set of CWRU [32]. The fault acquisition equipment is shown in Figure 19 [33] and the bearing structure coefficient is shown in Table 5. The motor speed is 1797 rpm. The sampling frequency is 12 KHz. The model used at the drive end is 6205-2RSJEMSKF, the motor speed is 1797 rpm, and the sampling frequency is 12 KHz.

Analysis of Inner Ring Vibration Signal
The time and frequency domain waveforms of the experimental signal are shown in Figure 20a,b. To test the performance of the method in terms of the anti-noise interference, this experiment adds noise to the original signal, and its SNR is −2dB. The signal's time and frequency domain diagram added with noise are shown in Figure 21a

Analysis of Inner Ring Vibration Signal
The time and frequency domain waveforms of the experimental signal are shown in Figure 20a,b. To test the performance of the method in terms of the anti-noise interference, this experiment adds noise to the original signal, and its SNR is −2dB. The signal's time and frequency domain diagram added with noise are shown in Figure 21a

Analysis of Inner Ring Vibration Signal
The time and frequency domain waveforms of the experimental signal are shown in Figure 20a,b. To test the performance of the method in terms of the anti-noise interference, this experiment adds noise to the original signal, and its SNR is −2dB. The signal's time and frequency domain diagram added with noise are shown in Figure 21a  Firstly, the MPA is used to optimize the decomposition level k and the penalty factor α according to the VMD parameter optimization process proposed in this paper. The fitness function used in MPA optimization is the average value of dispersion entropy. The parameter search range in VMD is as follows: k ∈ [3,15] and α ∈ [100, 3000]. In the MPA algorithm, the population size is 10. Meanwhile, the maximum iteration number is 20 andthe FADs are 0.2. After parameter optimization, the fitness curve of the obtained MPA is shown in Figure 22. According to the results shown in Figure 22, when the fitness value is optimal, the decomposition level k and penalty factor α are 3 and 1916, respectively. Therefore, k = 3 and α = 1916 are taken as the optimal values. Then, the above two parameters are substituted into the VMD algorithm for modeling, and VMD decomposes the bearing inner race fault signal after adding noise. The decomposition result is shown in Figure 23.  Firstly, the MPA is used to optimize the decomposition level k and the penalty factor αaccording to the VMD parameter optimization process proposed in this paper. The fitness function used in MPA optimization is the average value of dispersion entropy. The parameter search range in VMD is as follows: k∊ [3,15] and α∊ [100,3000]. In the MPA algorithm, the population size is 10. Meanwhile, the maximum iteration number is 20 andthe FADs are 0.2. After parameter optimization, the fitness curve of the obtained MPA is shown in Figure 22. According to the results shown in Figure 22, when the fitness value is optimal, the decomposition level k and penalty factor α are 3 and 1916, respectively. Therefore, k = 3 and α = 1916 are taken as the optimal values. Then, the above two parameters are substituted into the VMD algorithm for modeling, and VMD decomposes the bearing inner race fault signal after adding noise. The decomposition result is shown in Figure 23. Next, the kurtosis and correlation values of the signal components decomposed by the above MPA-VMD method are calculated and the K-C index of each IMF are calculated according to the formula of the combined weights. The calculation results of the signal components are shown in Table 6. Comparing the K-C index of each signal component in Table 6 indicates that the K-C index of the IMF2 component obtained based on the MPA-VMD method is the largest. Therefore, IMF2 is selected as the optimal signal component for the next TVD noise reduction operation. Firstly, the MPA is used to optimize the decomposition level k and the penalty factor αaccording to the VMD parameter optimization process proposed in this paper. The fitness function used in MPA optimization is the average value of dispersion entropy. The parameter search range in VMD is as follows: k∊ [3,15] and α∊ [100,3000]. In the MPA algorithm, the population size is 10. Meanwhile, the maximum iteration number is 20 andthe FADs are 0.2. After parameter optimization, the fitness curve of the obtained MPA is shown in Figure 22. According to the results shown in Figure 22, when the fitness value is optimal, the decomposition level k and penalty factor α are 3 and 1916, respectively. Therefore, k = 3 and α = 1916 are taken as the optimal values. Then, the above two parameters are substituted into the VMD algorithm for modeling, and VMD decomposes the bearing inner race fault signal after adding noise. The decomposition result is shown in Figure 23. Next, the kurtosis and correlation values of the signal components decomposed by the above MPA-VMD method are calculated and the K-C index of each IMF are calculated according to the formula of the combined weights. The calculation results of the signal components are shown in Table 6. Comparing the K-C index of each signal component in Table 6 indicates that the K-C index of the IMF2 component obtained based on the MPA-VMD method is the largest. Therefore, IMF2 is selected as the optimal signal component for the next TVD noise reduction operation. Next, the kurtosis and correlation values of the signal components decomposed by the above MPA-VMD method are calculated and the K-C index of each IMF are calculated according to the formula of the combined weights. The calculation results of the signal components are shown in Table 6. Comparing the K-C index of each signal component in Table 6 indicates that the K-C index of the IMF2 component obtained based on the MPA-VMD method is the largest. Therefore, IMF2 is selected as the optimal signal component for the next TVD noise reduction operation.  Figure 24 is the noise reduction waveform of the optimal signal components filtered based on the MPA-VMD method. It shows that the noise interference is reduced to a certain extent by the TVD method. Then, the signal after noise reduction is analyzed by the Hilbert envelope. The following conclusions are drawn from the envelope spectrum in Figure 25: the fault feature frequency and its double frequency appear in the figure. Still, their amplitude is low, complicating the fault characteristics' effective extraction. Next, the noise-reduced signal is input into the CYCBD filter. In this algorithm, the cyclic frequency set is set to [162. 19 Figure 24 is the noise reduction waveform of the optimal signal components filtered based on the MPA-VMD method. It shows that the noise interference is reduced to a certain extent by the TVD method. Then, the signal after noise reduction is analyzed by the Hilbert envelope. The following conclusions are drawn from the envelope spectrum in Figure 25: the fault feature frequency and its double frequency appear in the figure. Still, their amplitude is low, complicating the fault characteristics' effective extraction. Next, the noise-reduced signal is input into the CYCBD filter. In this algorithm, the cyclic frequency set is set to [162. 19, 324.38, 486.57, …,1135.33]. The signal filtered by the CYCBD filter is shown in Figure 26 and the generated envelope spectrum is shown in Figure 27. Observing the waveform shows that the periodic impact component has been significantly enhanced. It is shown in Figure 27 that the frequency characteristic of the fault and its double frequency, three-times frequency, four-times frequency, and five-times frequency components have appeared, and the   Figure 24 is the noise reduction waveform of the optimal signal components filtered based on the MPA-VMD method. It shows that the noise interference is reduced to a certain extent by the TVD method. Then, the signal after noise reduction is analyzed by the Hilbert envelope. The following conclusions are drawn from the envelope spectrum in Figure 25: the fault feature frequency and its double frequency appear in the figure. Still, their amplitude is low, complicating the fault characteristics' effective extraction. Next, the noise-reduced signal is input into the CYCBD filter. In this algorithm, the cyclic frequency set is set to [162. 19, 324.38, 486.57, …,1135.33]. The signal filtered by the CYCBD filter is shown in Figure 26 and the generated envelope spectrum is shown in Figure 27. Observing the waveform shows that the periodic impact component has been significantly enhanced. It is shown in Figure 27 that the frequency characteristic of the fault and its double frequency, three-times frequency, four-times frequency, and five-times frequency components have appeared, and the The signal filtered by the CYCBD filter is shown in Figure 26 and the generated envelope spectrum is shown in Figure 27. Observing the waveform shows that the periodic impact component has been significantly enhanced. It is shown in Figure 27 that the frequency characteristic of the fault and its double frequency, three-times frequency, fourtimes frequency, and five-times frequency components have appeared, and the amplitude is exceptionally high. Therefore, it can be determined that the bearing has inner race failure.

Analysis of Outer Ring Vibration Signal
The time and frequency domain waveforms of the experimental signal are shown in Figures 28a and 29b. Similarly, to verify the method's performance in anti-noise interference, the experiment adds noise to the original signal, and its SNR is−2dB. The time and frequency diagram of the signal after adding noise are shown in Figure 29a,b.
Next, the MPA optimization algorithm is used to optimize the k and α in VMD. Similarly, the parameter search range in VMD is as follows: the decomposition level k ∊ [3,15] and the penalty factor α∊ [100, 3000]. In the MPA algorithm, the population size is 10, the maximum iteration is 20, and the FADs are 0.2. After parameter optimization, the fitness curve of the obtained MPA is shown in Figure 30. It is shown from the change in the results of the fitness curve that when the fitness value is optimal, the corresponding optimal k and α are 3 and 2541, respectively. Therefore, the above two parameters are substituted into the VMD for modeling, and then VMD is carried out for the simulation experiment signal added with excessive noise. The decomposition result is shown in Figure 31. amplitude is exceptionally high. Therefore, it can be determined that the bearing has inner race failure.

Analysis of Outer Ring Vibration Signal
The time and frequency domain waveforms of the experimental signal are shown in Figures 28a and 29b. Similarly, to verify the method's performance in anti-noise interference, the experiment adds noise to the original signal, and its SNR is−2dB. The time and frequency diagram of the signal after adding noise are shown in Figure 29a,b.
Next, the MPA optimization algorithm is used to optimize the k and α in VMD. Similarly, the parameter search range in VMD is as follows: the decomposition level k ∊ [3,15] and the penalty factor α∊ [100, 3000]. In the MPA algorithm, the population size is 10, the maximum iteration is 20, and the FADs are 0.2. After parameter optimization, the fitness curve of the obtained MPA is shown in Figure 30. It is shown from the change in the results of the fitness curve that when the fitness value is optimal, the corresponding optimal k and α are 3 and 2541, respectively. Therefore, the above two parameters are substituted into the VMD for modeling, and then VMD is carried out for the simulation experiment signal added with excessive noise. The decomposition result is shown in Figure 31.

Analysis of Outer Ring Vibration Signal
The time and frequency domain waveforms of the experimental signal are shown in Figures 28a and 29b. Similarly, to verify the method's performance in anti-noise interference, the experiment adds noise to the original signal, and its SNR is−2dB. The time and frequency diagram of the signal after adding noise are shown in Figure 29a  Next, the MPA optimization algorithm is used to optimize the k and α in VMD. Similarly, the parameter search range in VMD is as follows: the decomposition level k ∈ [3,15] and the penalty factor α ∈ [100, 3000]. In the MPA algorithm, the population size is 10, the maximum iteration is 20, and the FADs are 0.2. After parameter optimization, the fitness curve of the obtained MPA is shown in Figure 30. It is shown from the change in the results of the fitness curve that when the fitness value is optimal, the corresponding optimal k and α are 3 and 2541, respectively. Therefore, the above two parameters are substituted into the VMD for modeling, and then VMD is carried out for the simulation experiment signal added with excessive noise. The decomposition result is shown in Figure 31.  Next, the kurtosis and correlation values of the signal components decomposed by the MPA-VMD method are calculated, and the K-C combination weight calculation formula calculates the K-C index of each IMF. The calculation results of the signal components are shown in Table 7. By comparing the K-C index of each signal component in Table 7, it is shown that the K-C index of the IMF2 component obtained based on the MPA-VMD method is the largest. Therefore, the IMF2 is selected as the optimal signal component for the next TVD noise reduction operation.  Figure 31.VMD result.
Next, the kurtosis and correlation values of the signal components decomposed by the MPA-VMD method are calculated, and the K-C combination weight calculation formula calculates the K-C index of each IMF. The calculation results of the signal components are shown in Table 7. By comparing the K-C index of each signal component in Table 7, it is shown that the K-C index of the IMF2 component obtained based on the MPA-VMD method is the largest. Therefore, the IMF2 is selected as the optimal signal component for the next TVD noise reduction operation.  Figure 32 shows the waveform after the noise reduction of the optimal signal components filtered based on the MPA-VMD method. Noise interference is shown to be significantly reduced by the TVD method. Then, a Hilbert envelope analysis is performed on the noise-reduced signal and the results are shown in Figure 33. It can be seen from Figure 33 that the fault characteristic frequency, double frequency, five-times frequency, and six-times frequency components appear in the envelope spectrum. However, their amplitude is low, and some irrelevant interference components appear around them, which makes it difficult to extract the fault features effectively. Next, the noise-reduced signal is input into the CYCBD filter. In this algorithm, the cyclic frequency set is set to    Figure 32 shows the waveform after the noise reduction of the optimal signal components filtered based on the MPA-VMD method. Noise interference is shown to be significantly reduced by the TVD method. Then, a Hilbert envelope analysis is performed on the noise-reduced signal and the results are shown in Figure 33. It can be seen from Figure 33 that the fault characteristic frequency, double frequency, five-times frequency, and six-times frequency components appear in the envelope spectrum. However, their amplitude is low, and some irrelevant interference components appear around them, which makes it difficult to extract the fault features effectively. Next, the noise-reduced signal is input into the CYCBD filter. In this algorithm, the cyclic frequency set is set to [ Figure 31.VMD result.
Next, the kurtosis and correlation values of the signal components decomposed by the MPA-VMD method are calculated, and the K-C combination weight calculation formula calculates the K-C index of each IMF. The calculation results of the signal components are shown in Table 7. By comparing the K-C index of each signal component in Table 7, it is shown that the K-C index of the IMF2 component obtained based on the MPA-VMD method is the largest. Therefore, the IMF2 is selected as the optimal signal component for the next TVD noise reduction operation.  Figure 32 shows the waveform after the noise reduction of the optimal signal components filtered based on the MPA-VMD method. Noise interference is shown to be significantly reduced by the TVD method. Then, a Hilbert envelope analysis is performed on the noise-reduced signal and the results are shown in Figure 33. It can be seen from Figure 33 that the fault characteristic frequency, double frequency, five-times frequency, and six-times frequency components appear in the envelope spectrum. However, their amplitude is low, and some irrelevant interference components appear around them, which makes it difficult to extract the fault features effectively. Next, the noise-reduced signal is input into the CYCBD filter. In this algorithm, the cyclic frequency set is set to [107.36, 214.72, 322.08, …,966.24].  The signal filtered by the CYCBD filter is shown in Figure 34 and the envelope spectrum generated is shown in Figure 35. It is shown that Figure 34 that the periodic impact component has been significantly improved. At the same time, the fault characteristic frequency and its second-to nine-times frequency components have appeared and the amplitude is exceptionally high. Therefore, it can be determined that the bearing has outer ring failure. The signal filtered by the CYCBD filter is shown in Figure 34 and the envelope spectrum generated is shown in Figure 35. It is shown that Figure 34 that the periodic impact component has been significantly improved. At the same time, the fault characteristic frequency and its second-to nine-times frequency components have appeared and the amplitude is exceptionally high. Therefore, it can be determined that the bearing has outer ring failure.  The signal filtered by the CYCBD filter is shown in Figure 34 and the envelope spectrum generated is shown in Figure 35. It is shown that Figure 34 that the periodic impact component has been significantly improved. At the same time, the fault characteristic frequency and its second-to nine-times frequency components have appeared and the amplitude is exceptionally high. Therefore, it can be determined that the bearing has outer ring failure.  The signal filtered by the CYCBD filter is shown in Figure 34 and the envelope spectrum generated is shown in Figure 35. It is shown that Figure 34 that the periodic impact component has been significantly improved. At the same time, the fault characteristic frequency and its second-to nine-times frequency components have appeared and the amplitude is exceptionally high. Therefore, it can be determined that the bearing has outer ring failure.

Conclusions
In this paper, the method combining improved VMD multi-scale dispersion entropy and TVD-CYCBD is used to research the fault feature extraction of rolling bearing signals under noise interference. The conclusions are as follows.
(1) The MPA is used to optimize the decomposition level and penalty factor in the VMD algorithm to find the optimal parameter combination. It can avoid over-decomposition and under-decomposition problems caused by the traditional VMD algorithm and realize the proper decomposition and the reconstruction of signals. This method overcomes the difficulties of mode aliasing and the end effect caused by impact components and noise interference in EMD and EEMD methods. (2) The K-C index is constructed by balancing the advantages and disadvantages of kurtosis and the cross-correlation coefficient, which solves the problems of identifying the optimal signal components. It effectively removes the component signals with a weak correlation with the original signal. Second, the TVD method reduces the noise of the selected optimal component. It effectively reduces the noise interference. Compared with the noise reduction results obtained based on EMD and EEMD methods, the waveform of the signal received by the proposed noise reduction method is more similar to the original signal. It is also optimal in comparing the noise reduction evaluation index results. (3) The fault impact component can be better highlighted by CYCBD filtering on the signal after TVD denoising. In the simulation signal experiment and the experimental verification of the CWRU data, the envelope spectrum generated by the proposed method can successfully extract the fault frequency and its multiple frequencies. At the same time, the result of fault feature extraction is better than that based on EMD and EEMD, which shows the effectiveness of the proposed method.