A Rolling Bearing Fault Feature Extraction Algorithm Based on IPOA-VMD and MOMEDA

Since the rolling bearing fault signal captured by a vibration sensor contains a large amount of background noise, fault features cannot be accurately extracted. To address this problem, a rolling bearing fault feature extraction algorithm based on improved pelican optimization algorithm (IPOA)–variable modal decomposition (VMD) and multipoint optimal minimum entropy deconvolution adjustment (MOMEDA) methods is proposed. Firstly, the pelican optimization algorithm (POA) was improved using a reverse learning strategy for dimensional-by-dimensional lens imaging and circle mapping, and the optimization performance of IPOA was verified. Secondly, the kurtosis-square envelope Gini coefficient criterion was used to select the optimal modal components from the decomposed components of the signal, and MOMEDA was used to process the optimal modal components in order to obtain the optimal deconvolution signal. Finally, the Teager energy operator (TEO) was employed to demodulate and analyze the optimally deconvoluted signal in order to enhance the transient shock component of the original fault signal. The effectiveness of the proposed method was verified using simulated and actual signals. The results showed that the proposed method can accurately extract failure characteristics in the presence of strong background noise interference.


Introduction
Rolling bearings play a crucial role in rotating mechanical equipment.Their health status directly affects the operational performance of mechanical equipment and even determines whether production processes can be performed [1].Rolling bearing fault diagnosis mainly contains two aspects.On the one hand, the original signal is preprocessed, and the time-frequency domain's fault features are extracted via signal processing methods.On the other hand, the extracted time-frequency domain features are used for fault identification via deep learning methods [2].However, due to the harsh operating environment of rolling bearings and the mutual friction of the transmission system, the vibration signal contains large amounts of ambient noise, and the fault characteristics are frequently drowned out and can be challenging to discern.Therefore, the effective extraction of fault characteristics of rolling bearings in the presence of strong background noise interference can greatly improve the training efficiency and recognition accuracy of deep learning fault diagnosis algorithms, which is of immense practical significance for the normal operation of mechanical equipment [3,4].
Rolling bearing fault signals exhibit nonlinear and non-smooth characteristics, limiting the effectiveness of traditional signal processing methods for noise reduction analysis [5,6].As such, advanced signal processing methods are required to filter out noise components from the original vibration signal [7,8].Commonly used methods for vibration signal denoising include wavelet transform (WT), empirical mode decomposition (EMD), and ensemble EMD (EEMD) [9].WT can characterize local signals but lacks self-adaptability and requires the manual selection of wavelet basis functions [10].EMD decomposes the signal into a series of intrinsic mode components (IMFs) that characterize the signal; however, it is susceptible to modal aliasing and endpoint effects during the signal decomposition process [11,12].EEMD distinguishes high-and low-frequency signals by using a mean square error criterion but suffers from long decomposition times and signal residue [13].Although the aforementioned methods achieve noise reduction in rolling bearing fault signals to a certain extent, they fail to eliminate modal aliasing and boundary effects encountered in the signal decomposition process.Dragomiretskiy et al. [14] proposed a variational modal decomposition (VMD) algorithm to achieve adaptive signal decomposition by solving the optimal solution of variational modes.This method effectively avoids boundary effects and modal aliasing in the main frequency component, and it is more conducive to the processing of complex signals.However, the decomposition effect of the VMD algorithm is mainly determined by the penalty factor and the number of modal layers, and the improper selection of parameters can result in over-decomposition and spurious components of the signal [15,16].Therefore, to ensure accurate parameter settings in VMD and avoid subjective errors, it is particularly important to adaptively select the optimal core parameters of VMD for subsequent feature extraction.Inspired by the application of swarm intelligence algorithms in the field of parameter optimization, this paper utilizes metaheuristic algorithms to adaptively determine the core parameters of VMD.
In recent years, scholars have applied swarm intelligence optimization algorithms to signal processing, which can not only adaptively establish the core parameters of signal processing algorithms but also improve the signal processing performance of the corresponding algorithms.Luo et al. [17] proposed an improved differential search (DS) optimization algorithm for the adaptive optimization of the core parameters of VMD and combined it with stochastic resonance theory to extract fault features.Zhang et al. [18] combined the grasshopper optimization algorithm (GOA) and the maximum weighted kurtosis index criterion to optimize the core parameters of VMD and extracted rolling bearing fault features.Yang et al. [19] used the marine predator optimization algorithm (MPA) to adaptively obtain the optimal parameters of the VMD and combined it with the fully variational denoising maximum second-order cyclic steady-state blind convolution (TVD-CYCBD) model for fault feature extraction.Ding et al. [20] proposed a VMD parameter optimization algorithm based on gene mutation particle swarm optimization (GMPSO).This method used GMPSO to obtain the optimal parameter combination of the VMD algorithm and then performed envelope spectrum analyses on the optimal modal components, finally extracting the fault features.Wang et al. [21] used an Archimedean optimization algorithm (AOA) to search for the optimal number of decomposition layers and penalty factor of the VMD and to find the IMF components that are the most sensitive to fault features.Although the above swarm intelligence algorithm can optimize the parameters of VMD to some extent, they are prone to fall into local optimization at the later stages of the iteration.Mei et al. [22] improved the pelican optimization algorithm (POA) using chaotic mapping to optimize the random forest model.However, this method only initialized the pelican population, which did not significantly improve the algorithm's optimization seeking ability.Therefore, in order to improve the convergence speed and accuracy of POA and enhance its ability with respect to global optimization searching, this paper proposes a new improvement strategy to adaptively adjust the position of pelican individuals in POA via circle mapping and reverse learning strategies for dimensional-by-dimensional lens imaging.
After the preliminary noise reduction in the original signal using the signal decomposition algorithm, methods for extracting the periodic shock signal from the signal components are also one of the key steps in fault feature extraction.To address this problem, Endo et al. [23] proposed a minimum entropy deconvolution (MED) algorithm.However, this method can only obtain a single pulse signal, while the rolling bearing fault signal is generally a periodic pulse signal.McDonald et al. [24] proposed a maximum correlation kurtosis deconvolution (MCKD) algorithm by combining correlation kurtosis and the deconvolution algorithm.However, the model of the algorithm is complex, and the accuracy of the model is determined by multiple parameters.McDonald et al. [25] proposed a multipoint optimal minimum entropy deconvolution-adjusted (MOMEDA) algorithm and employed it for the fault diagnosis of rolling bearings [26][27][28].MOMEDA eliminates the need for preset failure cycles, and the model is highly generalizable.MOMEDA overcomes the limitations of MED and MCKD, which require constant iterations to obtain the optimal filter, and it can effectively enhance and extract the periodic shock components in vibration signals.
To effectively and precisely extract failure characteristics in the presence of strong background noise, we propose applying the IPOA-VMD and MOMEDA algorithms to the fault detection of rolling bearings.Firstly, the POA is improved using circle mapping and a reverse learning strategy for dimensional-by-dimensional lens imaging, and the core parameters of VMD are adaptively optimized using IPOA.Secondly, the kurtosissquare envelope Gini coefficient (K-SEGI) criterion is used to select the optimal modal components from the decomposed components of the signal, and MOMEDA is used to further denoise the optimal modal components to obtain the optimal deconvolution signal.Finally, the Teager energy operator (TEO) is employed in demodulation, and the optimally deconvoluted signal is analyzed, achieving the accurate extraction of fault features.We verified the effectiveness of the proposed method for failure characteristic extraction in the presence of strong background noise interference using simulated and actual signals.

Variational Mode Decomposition (VMD)
The constrained variational model constructed by the VMD is described as follows [29]: where {u k } refers to each IMF component that is derived from the VMD decomposition, and f denotes the original signal.The {ω k } represents the center frequency associated with each IMF component, δ(t) represents the unit pulse function, and * represents convolutional operations.
Using an augmented Lagrange function, we can transform a restricted variational problem into an unrestricted variational problem.The expression is as follows: where λ is the Lagrange multiplier, and α is the punishment factor.The alternating direction method of multipliers (ADMMs) is employed to address the problem of solving the unconstrained variational model.The update formulas for {ω k }, {u k }, and λ, as well as the iteration termination conditions, are as follows: Sensors 2023, 23, 8620 2 where û n+1 k , λn+1 , and f are, respectively, denote the Fourier transform corresponding to u n+1 k , λ n+1 , and f ; ε is the convergence accuracy.Cyclic updating is performed using the aforementioned steps until the termination condition is met, finally yielding the K IMF components.

Improved Pelican Optimization Algorithm (IPOA)
POA achieves optimal searches by simulating the hunting process of pelican populations in terrestrial organisms [30].Compared to commonly used optimization algorithms [31][32][33], POA has advantages, such as fast convergence speed, resilience to local optima, and strong approximation ability to find optimal solutions [34].POA comprises two main phases: the exploration phase and the development phase.However, the problem of decreasing population diversity occurs at the later stage of iterations.Therefore, in this paper, we proposed the use of circle mapping and reverse learning strategies for dimensional-by-dimensional lens imaging to enhance population diversity.First, circle mapping is employed to initialize the positions of the pelican population so that each pelican is evenly distributed across the entire search space.Next, the position of each pelican is optimized using the reverse learning strategy for dimensional-by-dimensional lens imaging to enhance the diversity of the pelican population during the later iteration process of the POA in order to improve the convergence rate of the algorithm and reduce the risk of trapping in local optimality.

Circle Mapping Strategy
Chaotic mapping exhibits properties such as randomness, ergodicity, and regularity, which can be utilized to improve the optimization performance of POA.Commonly used chaotic mappings in the field of optimization include logistic mappings, tent mappings, and circle mappings [35].As observed in the distribution of the 1500 sequence values generated using four different methods, as shown in Figure 1, circle mapping generated a more uniform and stable distribution of chaotic sequence values between 0 and 1 compared to logistic mapping, tent mapping, and ordinary random numbers.Therefore, in this study, circle mapping was used to initialize the positions of the pelican population.Circle mapping can be defined as follows: where num i is the i-th chaotic sequence number, and mod is the remainder operation.The expression for IPOA population initialization is as follows: where num ij is the value of the chaotic sequence generated via circle mapping on the dimension of the corresponding pelican individual.
where i num is the i-th chaotic sequence number, and mod is the remainder operation.The expression for IPOA population initialization is as follows: where ij num is the value of the chaotic sequence generated via circle mapping on the dimension of the corresponding pelican individual.

Reverse Learning Strategy for Dimensional-by-Dimensional Lens Imaging
Two methods are mainly used to solve the local optimum problem: (1) maintaining the current optimal position and expanding the search area and (2) abandoning the current optimal position and searching in a new area.In this study, the first method was used, inspired by lens imaging [36], and the reverse learning strategy for dimensional-by-

Reverse Learning Strategy for Dimensional-by-Dimensional Lens Imaging
Two methods are mainly used to solve the local optimum problem: (1) maintaining the current optimal position and expanding the search area and (2) abandoning the current optimal position and searching in a new area.In this study, the first method was used, inspired by lens imaging [36], and the reverse learning strategy for dimensional-by-dimensional lens imaging was used to facilitate the POA in escaping from the local optimal region.Figure 2 illustrates the schematic of the reverse learning strategy for lens imaging [36].
For a space with a search range of [a j , b j ] for feasible solutions, the position of the optimal individual X j best in the j-th dimension represents the projection of an object p with a height of h on the x-axis.A convex lens is placed at base point o.An object p creates an inverted real image p' with height h' on the other side of the convex lens.At this point, the projection of p' on the x-axis is represented as X j best .The following expression is based on the principle of lens imaging: Sensors 2023, 23, 8620

of 29
Let h/h = n.Accordingly, Equation ( 9) can be transformed as follows: When n = 1, we obtain the following: As observed in Equation (11), when n = 1, a fixed reverse solution is obtained.Thus, the dynamically varying inverse solution is obtained by adjusting the value of n.First, the optimal individual position of the pelican is updated using the proposed optimization strategy, mapping the positions of each dimension into space to obtain the reverse position.Next, the previous position's fitness value is compared with the fitness value after reverse learning.If the fitness value after reverse learning is better than that of the previous position, the reverse position is selected to replace the previous position.Otherwise, the original position is retained for the next generation.The reverse learning strategy for dimensional-by-dimensional lens imaging not only avoids interference between different dimensions but also expands the search range of the algorithm.dimensional lens imaging was used to facilitate the POA in escaping from the local optimal region.Figure 2 illustrates the schematic of the reverse learning strategy for lens imaging [36].For a space with a search range of [aj, bj] for feasible solutions, the position of the optimal individual j b est X in the j-th dimension represents the projection of an object p with a height of h on the x-axis.A convex lens is placed at base point o.An object p creates an inverted real image p' with height h' on the other side of the convex lens.At this point, the projection of p' on the x-axis is represented as ′ j b est X .The following expression is based on the principle of lens imaging: Let ′ = h h n .Accordingly, Equation ( 9) can be transformed as follows: When n = 1, we obtain the following: As observed in Equation ( 11), when n = 1, a fixed reverse solution is obtained.Thus, the dynamically varying inverse solution is obtained by adjusting the value of n.First, the

IPOA Algorithm Performance Experiment
To verify the optimization performance of IPOA, the swarm intelligence algorithms in the literature [17][18][19][20][21] were compared with IPOA.Multiple benchmark functions were used for testing.The benchmark test function parameters are presented in Table 1.During the testing process, for all the algorithms, the limit for the number of iterations was 100, the initial population was set as 30, and the search space dimension of the population was 30.To avoid the randomness of the experiment, each algorithm was run independently 30 times for each benchmark function.The results of 30 experiments were counted, and the average value Avg and the total number of times to reach the target optimal value MR (the target optimum value was set to 10 −6 ) were calculated.The smaller the value of Avg, the stronger the ability of the optimization algorithm to approximate the optimal solution and the greater the probability of reaching the global optimum.The larger the value of MR, the higher the optimization accuracy; moreover, convergence speeds become faster, and the optimization algorithm exhibits stronger reliability.The iterative curves of different optimization algorithms under single mode and multimodal benchmark test functions are Sensors 2023, 23, 8620 7 of 29 displayed in Figures 3 and 4, respectively.The experimental results of different optimization algorithms are shown in Table 2.
times for each benchmark function.The results of 30 experiments were counted, and the average value Avg and the total number of times to reach the target optimal value MR (the target optimum value was set to 10 −6 ) were calculated.The smaller the value of Avg, the stronger the ability of the optimization algorithm to approximate the optimal solution and the greater the probability of reaching the global optimum.The larger the value of MR, the higher the optimization accuracy; moreover, convergence speeds become faster, and the optimization algorithm exhibits stronger reliability.The iterative curves of different optimization algorithms under single mode and multimodal benchmark test functions are displayed in Figures 3 and 4, respectively.The experimental results of different optimization algorithms are shown in Table 2.  Table 1.Benchmarking function parameters.

Benchmark Function Expression Dimension Value Range Optimal Solution
30 [−5.12,5.12]0 ( ) The above chart clearly shows the test results of different optimization algorithms.As observed in Figures 3 and 4 and Table 2, IPOA converged stably to the optimal solution for single-mode benchmark test functions F1 and F2.It exhibited a better convergence rate, optimization accuracy, and stability compared to the optimization algorithm proposed in the literature [17][18][19][20][21].For multimodal benchmark test functions F3 and F4, IPOA also converged stably to the optimal solution.By employing a reverse learning strategy to avoid falling into local optima in the later stages, the optimization effect was significantly improved.In summary, IPOA exhibits strong stability and robustness.The above chart clearly shows the test results of different optimization algorithms.As observed in Figures 3 and 4 and Table 2, IPOA converged stably to the optimal solution for single-mode benchmark test functions F 1 and F 2 .It exhibited a better convergence rate, optimization accuracy, and stability compared to the optimization algorithm proposed in the literature [17][18][19][20][21].For multimodal benchmark test functions F 3 and F 4 , IPOA also converged stably to the optimal solution.By employing a reverse learning strategy to avoid falling into local optima in the later stages, the optimization effect was significantly improved.In summary, IPOA exhibits strong stability and robustness.

Envelope Spectral Entropy (ESE)
To optimize the core parameters of VMD in IPOA, an appropriate fitness function must be constructed.We introduced the concept of entropy to ensure rational parameter selection.Entropy is an indicator used to reflect the sparsity characteristics of a signal.A smaller entropy value indicates that the sequence contains more meaningful information and is smoother [37].Among information entropy metrics, envelope spectral entropy (ESE) has the characteristics of simple calculation and fewer parameter inputs [38].Therefore, in this paper, ESE was employed as the fitness function for the IPOA-VMD optimization model.The expression for ESE is as follows: where a(i) is the envelope signal of the original signal after Hilbert demodulation, the length of the signal can be denoted as N, and ε(i) is the normalized form of a(i).

IPOA-VMD Optimization Flow
The flow of IPOA for optimizing the core parameters of VMD is illustrated in Figure 5.The concrete steps can be outlined as follows: Sensors 2023, 23, 8620 9 of 29 Step 1: The maximum number of iterations, spatial dimensions, population size, decomposition layers, and the penalty factor for IPOA are set, and the circle mapping strategy is used to initialize the population's position.
Step 2: The VMD decomposition of the original signal yields several IMF components.The ESE value is calculated for each IMF component, and the fitness function for global exploration is determined by selecting the component with the lowest ESE value.
Step 3: After each round of iteration, the ESE value corresponding to each set of parameter combinations is calculated and compared with the current ESE value.If it is less than the current ESE value, the current ESE value will be renewed.
Step 4: Whether the iteration's stop condition has been reached is determined.If the maximum number of iterations is not reached, let t = t + 1.In addition, the pelican population position updated in the previous iteration is used as the initial population position in the next round.Steps 1-5 are repeated until the iteration condition is reached.
Step 5: When the loop iteration ends, the optimal parameter combination will be outputted.

Component Screening
Selecting the optimal component from multiple IMF components is important for accurately extracting fault features.Metrics such as kurtosis, correlation coefficient, and sample entropy are generally used to screen valid components [18].However, a single screening metric is susceptible to noise interference, resulting in the misidentification of valid components [39].
The Gini coefficient is an index used for measuring the sparsity of a sequence [40], and it has been used in the fault diagnosis of rolling bearings due to its high stability  Selecting the optimal component from multiple IMF components is important for accurately extracting fault features.Metrics such as kurtosis, correlation coefficient, and sample entropy are generally used to screen valid components [18].However, a single screening metric is susceptible to noise interference, resulting in the misidentification of valid components [39].
The Gini coefficient is an index used for measuring the sparsity of a sequence [40], and it has been used in the fault diagnosis of rolling bearings due to its high stability against noise disturbances [41,42].To better screen the effective components, in this study, the squared envelope Gini coefficient (SEGI) was used.The value of SEGI is between 0 and 1; the closer it is to 1, the better the balance of the sequence.The expression for SEGI is as follows: where SE represents the squared envelope of the original signal, and SE 1 is the L1 paradigm of the SE.
Kurtosis is sensitive to shock signals and is a good indicator for detecting periodic shocks.The expression for kurtosis can be described as follows: where x is the expected value, the value of N represents the total number of signal points, and σ denotes the standard deviation.
To utilize the advantages of these two indicators, we proposed a screening criterion called K-SEGI.Due to the different dimensions and value ranges of these two indicators, it is necessary to standardize their amplitudes.First, the amplitude is normalized.Next, the normalized amplitude is exponentially increased based on a base of 2, and the resulting value is used as the final amplitude.The calculation formula for the K-SEGI screening criterion is as follows: where K' and SEGI' are the amplitudes of the two indicators after normalization.

MOMEDA Algorithm
MOMEDA is a weak signal enhancement method and is a non-iterative deconvolution process for obtaining optimal finite impulse response (FIR) filters [25].Assuming that the original vibration signal is x, the following expression holds: where y is the periodic shock signal, h is the transfer function, e is the ambient noise, and * denotes the convolution operation.MOMEDA recovers the periodic impulse signal y by searching for the optimal FIR filter.The process of solving the optimal filter can be translated into finding the maximum value of the multipoint D-paradigm number (MDN): where fault period T can be defined as the ratio between the sampling frequency and the eigenfrequency of the fault, and t represents the target vector that signifies the position and weight of the deconvolution target shock signal.When MDN reaches the maximum value, the following expression holds: By taking the derivative of Equation (19) and setting it to 0, the optimal filter for MOMEDA can be obtained as follows: In addition, Substituting Equation ( 21) into Equation ( 20) yields Let Because y = X 0 T f and assuming (X 0 X 0 T ) −1 exists, we obtain In summary, the optimal filter of MOMEDA can be expressed as follows: Features extraction using only MOMEDA is less effective due to the interference of strong background noise.To address this problem, TEO [43] is used to demodulate and analyze the optimal deconvolution signal obtained after MOMEDA processing to further suppress the interference of noise in order to accurately extract fault characteristics.

Fault Feature Extraction Method Process
The method for extracting fault features is outlined in Figure 6.The method consists of the following steps: Step 1: VMD core parameters are optimized using IPOA, and the original fault signal is adaptively decomposed.
Step 2: The K-SEGI value for each IMF component after signal decomposition is calculated.
Step 3: The modal component with the highest K-SEGI value is chosen as the optimal IMF component.
Step 4: The optimal modal component is further denoised using MOMEDA to obtain the optimal deconvolution signal.
Step 5: The optimal deconvolution signal is analyzed using TEO demodulation to extract the bearing's fault features.

culated.
Step 3: The modal component with the highest K-SEGI value is chosen as the optimal IMF component.
Step 4: The optimal modal component is further denoised using MOMEDA to obtain the optimal deconvolution signal.
Step 5: The optimal deconvolution signal is analyzed using TEO demodulation to extract the bearing's fault features.
where n(t) represents the Gaussian white noise signal, y(t) denotes the simulated signal for bearing faults, h(t) represents the periodic shock signal, t represents the sampling time, T represents the repetition period, A 0 represents the displacement constant, α represents the attenuation constant, f r is the rotation frequency of the bearing, and f n is the intrinsic frequency.
For A 0 = 0.01, α = 700, and t = 0.2 s, the number of sampling points is 2400, intrinsic frequency is f n = 4 kHz, bearing rotation frequency is f r = 30 Hz, repetition period is T = 0.0067 s, characteristic frequency of the bearing's failure is f i = 149.25 Hz, and sampling frequency is f s = 12 kHz.Figure 7 illustrates the basic characteristics of the periodic shock signal.
To create a simulation of high background noise environments, Gaussian white noise with a signal-to-noise ratio (SNR) of −15 dB was introduced to the periodic shock signal.Figure 8 displays the basic characteristics of the simulated signal.As depicted in Figure 8a, the characteristics of the periodic shock signal are heavily contaminated by Gaussian white noise signals.Furthermore, as depicted in Figure 8b, noise submerged the characteristic frequency and its octave in the simulated signal, making it impossible to extract fault features from the simulated signal.
resents the repetition period, A0 represents the displacement constant, α represents the attenuation constant, fr is the rotation frequency of the bearing, and fn is the intrinsic frequency.
For A0 = 0.01, α = 700, and t = 0.2 s, the number of sampling points is 2400, intrinsic frequency is fn = 4 kHz, bearing rotation frequency is fr = 30 Hz, repetition period is T = 0.0067 s, characteristic frequency of the bearing's failure is fi = 149.25 Hz, and sampling frequency is fs = 12 kHz.Figure 7 illustrates the basic characteristics of the periodic shock signal.To create a simulation of high background noise environments, Gaussian white noise with a signal-to-noise ratio (SNR) of −15 dB was introduced to the periodic shock signal.Figure 8 displays the basic characteristics of the simulated signal.As depicted in Figure 8a, the characteristics of the periodic shock signal are heavily contaminated by Gaussian white noise signals.Furthermore, as depicted in Figure 8b, noise submerged the characteristic frequency and its octave in the simulated signal, making it impossible to extract fault features from the simulated signal.

IPOA-VMD Signal Decomposition
The adaptive decomposition of the simulated signal was performed using IPOA-VMD.The initialization parameters for IPOA are presented in Table 3.The optimization iteration curve of IPOA-VMD is shown in Figure 9.

IPOA-VMD Signal Decomposition
The adaptive decomposition of the simulated signal was performed using IPOA-VMD.The initialization parameters for IPOA are presented in Table 3.The optimization iteration curve of IPOA-VMD is shown in Figure 9.
As illustrated in Figure 9, the ESE value decreased with the increase in the number of iterations during the optimization process, and the minimum ESE value occurred in the 15th generation, after which the minimum ESE value remained constant.When the loop iteration ended, the optimal combination of parameters was [4,19,219].The VMD decomposition of the simulated signal by using the optimal parameter combinations yielded four IMF components.The waveforms of these components are displayed in Figure 10a,b, both in the time and frequency domains.

IPOA-VMD Signal Decomposition
The adaptive decomposition of the simulated signal was performed using IPOA-VMD.The initialization parameters for IPOA are presented in Table 3.The optimization iteration curve of IPOA-VMD is shown in Figure 9.As illustrated in Figure 10b, the simulated signals were distributed in different frequency bands after the adaptive decomposition of IPOA-VMD, effectively avoiding the boundary effects and modal aliasing phenomena.To accurately select the optimal modal components from each IMF component, the evaluation parameters of each IMF component were calculated and are presented in Table 4.As evident from Table 4, the K-SEGI value of IMF1 was the largest, indicating that this component contained more shock components and less noise content; thus, it was selected as the optimal modal component.The kurtosis index selected IMF1 as the optimal modal component, while the SEGI index selected IMF4 as the optimal modal component.This indicated that a single screening indicator can easily lead to the misidentification of the optimal modal component, which proved the superiority of the comprehensive screening indicator proposed in this paper.As illustrated in Figure 9, the ESE value decreased with the increase in the number of iterations during the optimization process, and the minimum ESE value occurred in the 15th generation, after which the minimum ESE value remained constant.When the loop iteration ended, the optimal combination of parameters was [4,19,219].The VMD decomposition of the simulated signal by using the optimal parameter combinations yielded four IMF components.The waveforms of these components are displayed in Figure 10a,b, both in the time and frequency domains.As illustrated in Figure 10b, the simulated signals were distributed in different frequency bands after the adaptive decomposition of IPOA-VMD, effectively avoiding the boundary effects and modal aliasing phenomena.To accurately select the optimal modal components from each IMF component, the evaluation parameters of each IMF component were calculated and are presented in Table 4.As evident from Table 4, the K-SEGI value of IMF1 was the largest, indicating that this component contained more shock components and less noise content; thus, it was selected as the optimal modal component.The kurtosis index selected IMF1 as the optimal modal component, while the SEGI index selected IMF4 as the optimal modal component.This indicated that a single screening indicator can easily lead to the misidentification of the optimal modal component, which proved the superiority of the comprehensive screening indicator proposed in this paper.

Fault Feature Extraction
IMF1 was further denoised using MOMEDA by setting the period of deconvolution as T = f s /f r = 80 and the filter length as L = 700.The optimal deconvolution signal timedomain waveform of IMF1 after MOMEDA processing is shown in Figure 11, and the results obtained via demodulation analysis using TEO are displayed in Figure 12.As observed in Figure 11, the periodic shock component of the simulated signal was very obvious, and the noise reduction effect was good.IMF1 was further denoised using MOMEDA by setting the period of deconvolution as T = fs/fr = 80 and the filter length as L = 700.The optimal deconvolution signal timedomain waveform of IMF1 after MOMEDA processing is shown in Figure 11, and the results obtained via demodulation analysis using TEO are displayed in Figure 12.As observed in Figure 11, the periodic shock component of the simulated signal was very obvious, and the noise reduction effect was good.The results analyzed for IMF1 using only TEO demodulation are depicted in Figure 13a, and the results processed using only MOMEDA are depicted in Figure 13b.By The results analyzed for IMF1 using only TEO demodulation are depicted in Figure 13a, and the results processed using only MOMEDA are depicted in Figure 13b.By comparing Figures 12 and 13, it can be observed that the analyses of the signal using only TEO demodulation and processing using only MOMEDA are not satisfactory, the fault characteristics are not obvious, and the peaks are not prominent.This suggested that using TEO or MOMEDA alone did not accurately extract the fault features, and the characteristic frequency is still disturbed by a small amount of noise.In contrast, the fault characteristics frequency of IMF1 obtained via demodulation analysis using MOMEDA combined with TEO were extremely obvious, and peak f i and its octave were very prominent, corresponding to f i -6f i .Therefore, the proposed method enables the precise extraction of fault characteristics.comparing Figures 12 and 13, it can be observed that the analyses of the signal using only TEO demodulation and processing using only MOMEDA are not satisfactory, the fault characteristics are not obvious, and the peaks are not prominent.This suggested that using TEO or MOMEDA alone did not accurately extract the fault features, and the characteristic frequency is still disturbed by a small amount of noise.In contrast, the fault characteristics frequency of IMF1 obtained via demodulation analysis using MOMEDA combined with TEO were extremely obvious, and peak fi and its octave were very prominent, corresponding to fi-6fi.Therefore, the proposed method enables the precise extraction of fault characteristics.

Signal Processing Performance of the Proposed Method
In order to prove the reliability of the method proposed in this paper, experiments were conducted by adding noise signals with an SNR of −16 dB, −17 dB, and −18 dB to the pure signal.According to the calculation formula of SNR, the noise power in the above four (including the noise signal with SNR of −15 dB) simulated signals was 31.6 times, 39.8 times, 50.1 times, and 63.1 times the signal power, respectively.Figure 14 shows the spectrum of the components of each simulated signal after IPOA-VMD decomposition, and Figure 15 shows the TEO spectrum of the optimal IMF of each simulated signal after MOMEDA processing.As observed in Figure 14, the IMFs are distributed in different frequency bands without mode aliasing.As observed in Figure 15, the peaks of fault characteristic frequency and its octave are extremely prominent.This proved the reliability of the proposed method in this paper at different SNR levels.During the experiment, when a noise signal with an SNR of −20 dB was added to the pure signal (the noise power was 100 times the signal power), it was observed that the signal processing time was too long and always exceeded the limit value As observed in Figure 15, the peaks of fault characteristic frequency and its octave are extremely prominent.This proved the reliability of the proposed method in this paper at different SNR levels.During the experiment, when a noise signal with an SNR of −20 dB was added to the pure signal (the noise power was 100 times the signal power), it was observed that the signal processing time was too long and always exceeded the limit value of the preset parameters.This indicated that the limits of the signal processing method proposed in this paper were about to be reached.

Measured Signal Analysis
The experimental data were obtained from the vibration acceleration signal collected by a university in the US on the rolling bearing fault experimental platform [44].Figure 16 illustrates the schematic of the experimental platform.The basic parameters of the bearing are presented in Table 5 [44].In this paper, the fault data of the inner race of the bearing at the driving end of four kinds of motors with different actual load power (0 kW, 0.735 kW, 1.47 kW, and 2.205 kW) are selected for experimental analyses.The motor load is generally measured using percentages, which are calculated as the ratio of the actual load power to the rated power of the motor multiplied by 100%.The number of sample points is 2400, the sampling frequency f s2 is 12 kHz, and the rated power of the motor is 1.47 kW.Table 4 lists the parameters, such as the speed, diameter, and motor load of the bearings.The motor rotation frequency f r and the characteristic frequency f i of the inner race fault were calculated as follows: where D represents the pitch diameter of the bear'ng's raceway, Z represents the number of balls, d represents the diameter of the ball, β represents the contact angle of the bearing, and n represents the motor's speed.
Sensors 2023, 23, x FOR PEER REVIEW 20 of 30 power to the rated power of the motor multiplied by 100%.The number of sample points is 2400, the sampling frequency fs2 is 12 kHz, and the rated power of the motor is 1.47 kW.Table 4 lists the parameters, such as the speed, diameter, and motor load of the bearings.The motor rotation frequency fr and the characteristic frequency fi of the inner race fault were calculated as follows: . ( ) where D represents the pitch diameter of the bear'ng's raceway, Z represents the number of balls, d represents the diameter of the ball, β represents the contact angle of the bearing, and n represents the motor's speed.Based on the parameters presented in Table 6, when the motor load was 0% (indicating that the motor was in a no-load state), the fault characteristic frequency was calculated using fi2 = 162.19Hz in Equation (27), and the motor rotational frequency was calculated as fr2 = 29.95Hz.To create a simulation of a high background noise environment, Gaussian white noise with a signal-to-noise ratio (SNR) of −6 dB was introduced to the actual signal.Figure 17 displays the basic characteristics of the inner ring fault signal.Based on the parameters presented in Table 6, when the motor load was 0% (indicating that the motor was in a no-load state), the fault characteristic frequency was calculated using f i2 = 162.19Hz in Equation (27), and the motor rotational frequency was calculated as f r2 = 29.95Hz.To create a simulation of a high background noise environment, Gaussian white noise with a signal-to-noise ratio (SNR) of −6 dB was introduced to the actual signal.Figure 17 displays the basic characteristics of the inner ring fault signal.The adaptive decomposition of inner-ring bearing fault signals was performed using IPOA-VMD.The initialization parameters for IPOA are presented in Table 7.The IPOA-VMD optimization search iteration curve is displayed in Figure 18.As observed in Figure 18, the minimum ESE value occurred in the 11th generation, after which the minimum ESE value remained constant.When the loop iteration ended, the optimal combination of parameters was [3,2536].We applied the optimal parameter combination to perform VMD decomposition on the inner ring fault signal, and three modal components were obtained.The waveforms of these components are displayed in Figure 19a,b, both in the time and frequency domains.As observed in Figure 18, the minimum ESE value occurred in the 11th generation, after which the minimum ESE value remained constant.When the loop iteration ended, the optimal combination of parameters was [3,2536].We applied the optimal parameter combination to perform VMD decomposition on the inner ring fault signal, and three modal components were obtained.The waveforms of these components are displayed in Figure 19a,b, both in the time and frequency domains.The evaluation parameters of each IMF component are presented in Table 8.As evident from Table 8, the K-SEGI value of IMF2 was the largest.Thus, IMF2 was selected as the optimal modal component.As observed in Figure 19a, the fault shock characteristics of IMF2 were more obvious; however, a small amount of noise interference remained.Thus, MOMEDA was used for further noise reduction in IMF2 by setting the period of deconvolution as T = fs2/fr2 = 73.99 and filter length as L = 700.The optimal deconvolution signal timedomain waveform of IMF2 after MOMEDA processing is shown in Figure 20, and the results obtained via performing demodulation analysis using TEO are shown in Figure 21.The evaluation parameters of each IMF component are presented in Table 8.As evident from Table 8, the K-SEGI value of IMF2 was the largest.Thus, IMF2 was selected as the optimal modal component.As observed in Figure 19a, the fault shock characteristics of IMF2 were more obvious; however, a small amount of noise interference remained.Thus, MOMEDA was used for further noise reduction in IMF2 by setting the period of deconvolution as T = f s2 /f r2 = 73.99 and filter length as L = 700.The optimal deconvolution signal time-domain waveform of IMF2 after MOMEDA processing is shown in Figure 20, and the results obtained via performing demodulation analysis using TEO are shown in Figure 21.The evaluation parameters of each IMF component are presented in Table 8.As evident from Table 8, the K-SEGI value of IMF2 was the largest.Thus, IMF2 was selected as the optimal modal component.As observed in Figure 19a, the fault shock characteristics of IMF2 were more obvious; however, a small amount of noise interference remained.Thus, MOMEDA was used for further noise reduction in IMF2 by setting the period of deconvolution as T = fs2/fr2 = 73.99 and filter length as L = 700.The optimal deconvolution signal timedomain waveform of IMF2 after MOMEDA processing is shown in Figure 20, and the results obtained via performing demodulation analysis using TEO are shown in Figure 21.To further illustrate the effectiveness of the method proposed in this paper, the rolling bearing inner ring fault data with motor loads of 50%, 100%, and 150% were selected for experiments.According to Equation ( 27), when the motor load is 50%, the fault characteristic frequency fi3 is 159.91 Hz, and the rotation frequency fr3 is 29.53 Hz.When the motor load is 100%, the fault characteristic frequency fi4 is 157.96Hz, and the rotation frequency fr4 is 29.17 Hz.When the motor load is 150%, the fault characteristic frequency fi5  To further illustrate the effectiveness of the method proposed in this paper, the rolling bearing inner ring fault data with motor loads of 50%, 100%, and 150% were selected for experiments.According to Equation ( 27), when the motor load is 50%, the fault characteristic frequency f i3 is 159.91 Hz, and the rotation frequency f r3 is 29.53 Hz.When the motor load is 100%, the fault characteristic frequency f i4 is 157.96Hz, and the rotation frequency f r4 is 29.17 Hz.When the motor load is 150%, the fault characteristic frequency f i5 is 156.12Hz, and the rotation frequency f r5 is 28.83 Hz.According to the above fault feature extraction process, Gaussian white noise with an SNR of −6 dB was added to the inner race signals of rolling bearings under different loads, and the parameters of the decomposition algorithm were consistent with those shown in Table 7. Figure 22 displays the time-domain waveforms and TEO spectrum of each optimal IMF after MOMEDA processing.As observed in Figure 22, when the motor load is 50%, the actual frequency of fault characteristics is 159.42Hz.When the motor load is 100%, the actual frequency of fault characteristics is 157.52 Hz.When the motor load is 150%, the actual frequency of fault characteristics is 156.21Hz.They are all very close to the theoretically calculated frequency, and the peak and octave frequencies are extremely prominent.This showed that the algorithm proposed in this paper can accurately extract the characteristics of the inner race fault signal under different motor load conditions.This proved the effectiveness and reliability of the proposed method.

Comparative Analysis of Different Feature Extraction Methods
To substantiate the excellence of the proposed method in fault feature extraction in the presence of strong background noise interference, the EEMD algorithm, fixed-parameter VMD algorithm (FP-VMD), and empirical wavelet decomposition-independent component analysis (EWT-ICA) algorithm were compared with the proposed method in this paper.The penalty factor in the FP-VMD algorithm was 20,000, and the number of decomposition layers was eight.The EWT-ICA algorithm first performed an EWT decomposition of the signal.Then, several IMF components were selected as observation signals for the ICA algorithm based on the screening index, and other IMF components were used as virtual noise channel signals for the algorithm.Finally, the obtained source signals were enveloped and demodulated to extract the fault features.To increase the credibility of the experiment, the actual signals and screening metrics selected are consistent with this paper.Figure 23 shows the optimal IMF envelope spectra obtained using the EEMD algorithm for rolling bearing inner race fault signals under different loads.Figure 24 shows the optimal IMF envelope spectrum obtained using the FP-VMD algorithm for rolling bearing inner race fault signals under different loads.Figure 25 shows the optimal source signal envelope spectrum obtained using the EWT-ICA algorithm for rolling bearing inner race fault signals under different loads.lated to extract the fault features.To increase the credibility of the experiment, the actual signals and screening metrics selected are consistent with this paper.Figure 23 shows the optimal IMF envelope spectra obtained using the EEMD algorithm for rolling bearing inner race fault signals under different loads.Figure 24 shows the optimal IMF envelope spectrum obtained using the FP-VMD algorithm for rolling bearing inner race fault signals under different loads.Figure 25 shows the optimal source signal envelope spectrum obtained using the EWT-ICA algorithm for rolling bearing inner race fault signals under different loads.As observed in Figures 23-25, the EEMD, PF-VMD, and EWT-ICA algorithms were only able to extract the rotational frequency of the motor for the signals of the inner race of the rolling bearing under different loads, and the characteristic frequency and its multiplicative frequency were swamped by noise.In addition, the EEMD algorithm was only able to extract the rotational frequency of the motor when there was no load.but all the other characteristic frequencies are submerged by noise.In contrast, as observed in

Conclusions
We combine the improved POA-VMD with MOMEDA-TEO to research the rolling bearing fault feature extraction method in the presence of strong background noise environments.The main conclusions of this paper are as follows: As observed in Figures 23-25, the EEMD, PF-VMD, and EWT-ICA algorithms were only able to extract the rotational frequency of the motor for the signals of the inner race of the rolling bearing under different loads, and the characteristic frequency and its multiplicative frequency were swamped by noise.In addition, the EEMD algorithm was only able to extract the rotational frequency of the motor when there was no load.but all the other characteristic frequencies are submerged by noise.In contrast, as observed in Figures 21 and 22, the actual characteristic frequency of the inner race fault approached the theoretical value in an extreme manner, and the fault characteristics were very obvious, enabling the accurate extraction of the characteristic frequency and its multiplicative frequency.The above results indicated that the EEMD, FP-VMD, and EWT-ICA algorithms cannot fully extract fault features in strong background noise environments, thus fully demonstrating the superiority of the method proposed in this paper.

Conclusions
We combine the improved POA-VMD with MOMEDA-TEO to research the rolling bearing fault feature extraction method in the presence of strong background noise environments.The main conclusions of this paper are as follows: (1) The POA is improved to form the IPOA using circle mapping and reverse lens imaging learning strategies, and the IPOA is used to optimize the penalty factor and the number of modal layers in the VMD algorithm to adaptively obtain the optimal parameter combination.This method overcomes the over-decomposition of the signal and modal aliasing problems caused by the improper setting of parameters based on human experience and subjectivity.(2) Based on the advantages of SEGI and kurtosis screening metrics, a new optimal modal component screening metric of K-SEGI is proposed.It can effectively screen the components with the most periodic shocks and the most stable components while avoiding the mis-selection problem of the optimal component caused by using a single screening metric of SEGI and kurtosis.(3) After the preliminary noise reduction in the original signal using the IPOA-VMD decomposition algorithm, the optimal IMF component is further denoised by using MOMEDA, and the impact signal is highlighted better using TEO demodulation analyses.We conducted fault extraction experiments using simulated signals and four measured signals of inner race fault under different loads.The experimental results show that the proposed method can accurately extract the fault frequency and its octave.In addition, compared with EEMD, FP-VMD, and EWT-ICA feature extraction algorithms, the proposed method exhibits obvious advantages with respect to fault feature extraction in strong background noise environments.

Figure 2 .
Figure 2. Schematic of the reverse learning strategy for lens imaging.

Figure 2 .
Figure 2. Schematic of the reverse learning strategy for lens imaging.

Figure 3 .
Figure 3. Iterative curves of different optimization algorithms under the single modal benchmark function: (a) iteration curve for benchmark test function F1 and (b) iteration curve for benchmark test function F2.

Figure 3 .
Figure 3. Iterative curves of different optimization algorithms under the single modal benchmark function: (a) iteration curve for benchmark test function F 1 and (b) iteration curve for benchmark test function F 2 .

Figure 4 .
Figure 4. Iterative curves of different optimization algorithms under multimodal benchmark functions: (a) iteration curve for benchmark test function F3 and (b) iteration curve for benchmark test function F4.

Figure 4 .
Figure 4. Iterative curves of different optimization algorithms under multimodal benchmark functions: (a) iteration curve for benchmark test function F 3 and (b) iteration curve for benchmark test function F 4 .

Figure 5 .
Figure 5. Flowchart of the IPOA for the optimization of VMD parameters.

Figure 5 .
Figure 5. Flowchart of the IPOA for the optimization of VMD parameters.

Figure 7 .
Figure 7. Basic characteristic diagram of the periodic shock signal: (a) time-domain diagram and (b) Hilbert envelope spectrogram.

Figure 7 .Figure 8 .
Figure 7. Basic characteristic diagram of the periodic shock signal: (a) time-domain diagram and (b) Hilbert envelope spectrogram.Sensors 2023, 23, x FOR PEER REVIEW 14 of 30

Figure 8 .
Figure 8. Basic characteristic diagram of the simulated signal: (a) time-domain diagram and (b) Hilbert envelope spectrogram.

Figure 8 .
Figure 8. Basic characteristic diagram of the simulated signal: (a) time-domain diagram and (b) Hilbert envelope spectrogram.

Figure 10 .
Figure 10.Time-frequency domain diagram of the IMF component after IPOA-VMD decomposition (simulated signal): (a) time domain and (b) frequency domain.

Figure 10 .
Figure 10.Time-frequency domain diagram of the IMF component after IPOA-VMD decomposition (simulated signal): (a) time domain and (b) frequency domain.

Figure 11 .
Figure 11.Time-domain diagram of the processed IMF1 using MOMEDA.

Figure 11 .
Figure 11.Time-domain diagram of the processed IMF1 using MOMEDA.

Figure 11 .
Figure 11.Time-domain diagram of the processed IMF1 using MOMEDA.

Figure 13 .
Figure 13.Results of IMF1 processing by different methods: (a) TEO demodulation analysis and (b) MOMEDA processing.4.1.4.Signal Processing Performance of the Proposed Method In order to prove the reliability of the method proposed in this paper, experiments were conducted by adding noise signals with an SNR of −16 dB, −17 dB, and −18 dB to the pure signal.According to the calculation formula of SNR, the noise power in the above

Figure 13 .
Figure 13.Results of IMF1 processing by different methods: (a) TEO demodulation analysis and (b) MOMEDA processing.

Figure 16 .
Figure 16.Schematic of the rolling bearing fault experimental platform.

Figure 16 .
Figure 16.Schematic of the rolling bearing fault experimental platform.

Figure 17 .
Figure 17.Basic characteristic diagram of inner ring fault signals: (a) time domain and (b) Hilbert envelope spectrum.

Figure 17 .
Figure 17.Basic characteristic diagram of inner ring fault signals: (a) time domain and (b) Hilbert envelope spectrum.

Figure 19 .
Figure 19.Time-and frequency-domain waveforms of each IMF component after IPOA-VMD decomposition (actual signal): (a) time domain and (b) frequency domain.
composition (actual signal): (a) time domain and (b) frequency domain.

Figure 24 .
Figure 24.Optimal IMF envelope spectra of rolling bearing inner race fault signals obtained using the FP-VMD algorithm under different loads: (a) motor load is 0%, (b) motor load is 50%, (c) motor load is 100%, and (d) motor load is 150%.
Figures 21 and 22, the actual characteristic frequency of the inner race fault approached the theoretical value in an extreme manner, and the fault characteristics were very obvious, enabling the accurate extraction of the characteristic frequency and its multiplicative frequency.The above results indicated that the EEMD, FP-VMD, and EWT-ICA algorithms cannot fully extract fault features in strong background noise environments, thus fully demonstrating the superiority of the method proposed in this paper.

Figure 24 .Figure 25 .
Figure 24.Optimal IMF envelope spectra of rolling bearing inner race fault signals obtained using the FP-VMD algorithm under different loads: (a) motor load is 0%, (b) motor load is 50%, (c) motor load is 100%, and (d) motor load is 150%.Sensors 2023, 23, x FOR PEER REVIEW 27 of 30

Figure 25 .
Figure 25.Optimal IMF envelope spectra of rolling bearing inner race fault signals obtained using the EWT-ICA algorithm under different loads: (a) motor load is 0%, (b) motor load is 50%, (c) motor load is 100%, and (d) motor load is 150%.

Table 2 .
Comparison of the test results of different optimization algorithms.

Table 2 .
Comparison of the test results of different optimization algorithms.

Table 4 .
Evaluation parameters of each IMF component (simulated signal).

Table 4 .
Evaluation parameters of each IMF component (simulated signal).

Table 6 .
Fundamental characteristics of the rolling bearing.

Table 6 .
Fundamental characteristics of the rolling bearing.

Table 8 .
Evaluation parameters of each IMF component (actual signal).

Table 8 .
Evaluation parameters of each IMF component (actual signal).