Research on Fault Feature Extraction Method of Rolling Bearing Based on SSA–VMD–MCKD

: In response to the problem that nonlinear and non-stationary rolling bearing fault signals are easily disturbed by noise, which leads to the difﬁculty of fault feature extraction, to take full advantage of the superiority of variational mode decomposition (VMD) in noise reduction, and of maximum correlation kurtosis deconvolution (MCKD) in highlighting continuous pulses masked by noise, a method based on sparrow search algorithm (SSA), VMD, and MCKD is proposed, namely, SSA–VM–MCKD, for rolling bearing faint fault extraction. To improve the feature extraction effect, the method uses the inverse of the peak factor squared of the envelope spectrum as the ﬁtness function, and the parameters to be determined in both algorithms are searched adaptively by SSA. Firstly, the parameter-optimized VMD is used to decompose the fault signal to obtain the intrinsic mode function (IMF) components, from which the optimal mode component is selected, and then the optimal component signal is deconvoluted by the parameter-optimized MCKD to enhance the periodic fault pulses in the optimal component signal, and ﬁnally extracts the rolling bearing fault characteristic frequency by envelope demodulation. Experiments on simulated signals and measured data show that the method can adaptively determine the parameters in VMD and MCKD, enhance the fault impact components in the signals, and effectively extract the fault characteristic frequencies of rolling bearings, with a success rate up to 100%, providing a new idea for rolling bearing fault feature extraction.


Introduction
Rolling bearings, as the main components of large industrial rotating machinery, play an important role in improving the efficiency of machinery and equipment, and their health condition directly affects the smooth operation of machinery and equipment [1]. However, rolling bearings often work in high temperature, high pressure, and a complex mechanical environment, and are highly susceptible to failure. Since the shock generated by the early failure of rolling bearings is very weak and easily disturbed by system noise, coupled with the complex vibration transmission path of rolling bearings, the extraction of rolling bearing fault characteristics is very difficult [2,3]. Therefore, it is of great significance to process the original vibration signals of rolling bearings and extract the fault characteristics accurately and effectively for the determination of their working conditions, the prevention of unexpected accidents, and the improvement of equipment operation efficiency [4].
The periodic force shocks generated by rolling bearings excite high-frequency resonance and convolution between the bearing system components, creating a bearing periodic shock signal [5]. To effectively deal with nonlinear and nonstationary rolling bearing fault signals and accurately extract fault features from the processed signals, domestic and foreign experts and scholars have proposed many effective methods. Among them, empirical mode decomposition (EMD), local mean decomposition (LMD), wavelet transform, etc., (1) Section 2 first introduces the basic principles of VMD, MCKD and SSA to optimize the VMD and MCKD parameters.

SSA-VMD-MCKD Method
Firstly, the basic principles and processes of optimizing VMD and MCKD parameters by VMD, MCKD and SSA are introduced, respectively.

Principle and Process of VMD
By improving Wiener filtering and Hilbert transform, VMD applies the variational problem-solving process in signal decomposition, and solves the bandwidth and center frequency of the modal function. It is a new non-recursive modal decomposition algorithm. It can effectively highlight the fluctuation characteristics of the signal in each frequencydomain and has excellent signal processing performance for rolling bearing faults [22]. The VMD algorithm can be divided into two main parts: the first part is the construction of the variational problem, and the second part is the solution of the variational problem. The construction of the variational problem can be described by the following Equation (1): (1) where {u k } = {u 1 , u 2 , . . . , u K } represents the K IMF components decomposed by the input signal f (t): u k (t) = A k (t) cos[ϕ k (t)]; A k (t) and ϕ k (t) are the instantaneous amplitude and phase angle of u k (t), respectively; {ω k } = {ω 1 , ω 2 , . . . , ω K } is the center frequency of the IMF set, ω k (t) = dϕ k (t)/dt; K is the number of decompositions; ∂ t [·] represents the partial derivative of the time parameter t; δ(t) is the impulse function; j is the imaginary unit; "*" indicates convolution; and the Hilbert-transformed one-sided spectrum of the modal components is δ(t) + j πt * u k (t) [23]. The second part is the solution of the variational problem. Here, a Lagrange multiplier λ that can increase the strictness of constraints and a quadratic penalty factor α that can reduce interference from Gaussian noise are introduced. The constrained variational problem in Equation (1) is turned into an unconstrained variational problem, and then solved [24]: By alternating the direction multiplier algorithm, one can iteratively update u n+1 k , ω n+1 k , and λ n+1 to obtain the "saddle point", that is, the optimal solution of Equation (1). The specific process is [25]: Step 1: initialize parametersû 1 k , ω 1 k ,λ 1 , α, and n; Step 2: k = k + 1, until k = K, update u k and ω k according to Formula (3): where Fourier transform is performed on u k and λ;û k andλ can be obtained; Step 3: update λ:λ where τ is the noise tolerance parameter.
Step 4: determine whether the conditions for stopping the iteration are met by Formula (5). If so, stop the iteration; if not, return to step 2: where ε > 0 is the given discrimination accuracy, which is usually set to 1 × 10 −6 .

Principle and Process of MCKD
The MCKD method aims at maximizing the correlation kurtosis and extracts the periodic impulse components of the signal by iteratively selecting a suitable filter for deconvolution operation. The correlation kurtosis expression is [26]:  (6) where M is the displacement number; T = f s / f fault is the impulse signal period; f s and f fault are the sampling frequency and fault characteristic frequency, respectively; m ∈ [0, M], the number of pulses becomes more with the increase of M, thus improving the fault feature extraction capability of the method, but too large an M leads to a decrease in computational accuracy; the periodic impulse signal y n (n = 1, 2, . . . , N) is the output signal, that is, the processed signal; N is the length of the signal, and the expression of y n is [27]: Electronics 2022, 11, 3404 5 of 24 where L is the filter length; f = [ f 1 , f 2 , . . . , f L ] is the filter coefficient; l ∈ [0, L]; and x n is the input signal. After filtering the signal, the objective function of MCKD is shown in Formula (8), which maximizes the correlation kurtosis: The maximum value of CK M (T) in Equation (8) can be obtained by solving Equation (9): The filter coefficient expression in matrix form is further obtained as: where In summary, the specific steps of MCKD are as follows [28]: Step 1: determine parameters such as L, T, and M; Step 2: calculate X 0 X T 0 and X mT of x n ; Step 3: calculate the filtered signal y n ; Step 4: calculate β and ψ m according to y n ; Step 5: update the filter coefficient f; Step 6: if ∆CK M (T) < ε before and after filtering, ε is the set threshold, then stop the iteration; otherwise return to Step 3.

SSA Optimizes Parameters of VMD and MCKD
The performance of VMD and MCKD methods is greatly affected by the parameter values. Among them, the performance of the VMD method is strongly influenced by parameters α and K: when α increases, the bandwidth of each IMF component becomes smaller and the attenuation becomes faster. On the contrary, when α decreases, the bandwidth Electronics 2022, 11, 3404 6 of 24 becomes larger and the attenuation becomes slower. Too large a setting of K will lead to over-decomposition, too small will lead to incomplete decomposition, leading the fault information to be ignored. The performance of the MCKD method is affected by L, T, and M: if L is set too small, the deconvolution accuracy will decrease, and conversely, the number of samples to be calculated in each iteration will increase, thereby prolonging the calculation time. The setting of the parameter T requires an a priori fault period; the fault pulses of the deconvoluted signal increase as the parameter M becomes larger [29]. If the above parameters are set subjectively by humans, the effect of fault feature extraction will deteriorate. Therefore, an optimization algorithm is required to adaptively determine each parameter in the above two methods.
Xue, J.K. proposed an optimization algorithm in 2020, namely, SSA. Compared with the algorithms such as the genetic algorithm and the ant colony algorithm, SSA has a better optimization effect, and has the advantages of high accuracy, strong robustness, fast convergence speed, and wide application range [30]. Therefore, to obtain better fault feature extraction, the parameters of VMD and MCKD are optimally sought using SSA to adaptively obtain the required parameter combinations.
In SSA, it is divided into three parts: explorer, follower, and early warning. The explorer is mainly used to determine the foraging search direction and area. The position update can be expressed as [31]: where X t i,j denotes the position of the i-th sparrow in the j-th dimension at the t-th iteration; exp(·) represents the exponential function with the natural constant e as the base; α ∈ [0, 1] is a random number; M is the maximum number of iterations; R 2 ∈ [0, 1] is the warning value; ST ∈ [0.5, 1] is the safety value; Q is a random number obeying normal distribution; L is a 1 × d-dimensional matrix, with all elements being 1. If R 2 < ST, it means that the explorer can continue to expand the search range. If R 2 ≥ ST, it indicates that the explorer needs to leave the current area and turn to a safe location to forage.
The follower's position update Formula is [32]: where X worst is the worst position on the board; X P is the optimal position of the explorer; A is a 1 × d-dimensional matrix with each element randomly assigned 1 or −1, and A + = A T AA T −1 . If i > N/2, it indicates the need to go to other areas for forage. Otherwise, it means that the follower has the chance to forage. The early-warning person is responsible for guarding against predation, and the location update Formula is [33]: Among them, X best indicates the best position; β is a random number used to control the step size; K ∈ [−1, 1] is a random number used to control the step size and moving direction; f i , f w , and f g denote the individual fitness, the worst fitness, and the optimal fitness, respectively. When using SSA to optimize VMD and MCKD parameters, the fitness function is the inverse of the peak factor squared of the envelope spectrum. The expression of the peak factor is: where X(z)(z = 1, 2, . . . , Z) is the signal envelope spectrum amplitude sequence. The larger the crest factor, the more obvious the fault characteristic. Therefore, the smaller the value of the fitness function, the stronger the periodic impact characteristic and the better the optimization effect.

Fault Feature Extraction Process
To effectively integrate the strengths of the VMD method for processing rolling bearing fault signals, the superiority of the MCKD method for extracting fault features, and to adaptively seek all parameters in the VMD and MCKD methods, a new rolling bearing fault feature extraction method, namely SSA-VMD-MCKD, is proposed. The SSA-VMD-MCKD method flow is shown in Figure 1. Step 2 Step 3 Step 4 Step 5 Initialize SSA-VMD-MCKD parameters According to the flow chart of the SSA-VMD-MCKD method, the specific steps are as follows: Step 1: acquisition of raw vibration signals using sensors on rolling bearings. According to the vibration signal, use SSA to perform adaptive optimization on the parameters of VMD and obtain the required combination   0 0 , K  . According to the literature [34] and combined with the actual situation, the number of sparrow populations in SSA is set According to the flow chart of the SSA-VMD-MCKD method, the specific steps are as follows: Step 1: acquisition of raw vibration signals using sensors on rolling bearings. According to the vibration signal, use SSA to perform adaptive optimization on the parameters of VMD and obtain the required combination [α 0 , K 0 ]. According to the literature [34] and combined with the actual situation, the number of sparrow populations in SSA is set to 30, the maximum number of iterations is 20, the proportion of explorers is 0.2, the safety value is 0.8, the optimization range of parameter α is [100,2000], and the optimization range of parameter K is [3,10].
Step 2: set the parameters in VMD according o the optimal parameter combination obtained in step 1. Perform VMD processing of the fault signal to obtain the required components of each IMF. Then, calculate the peak factor value of the envelope spectrum of each component according to Formula (14) and select the component with the largest value as the required optimal component. Envelope demodulation analysis is performed on the optimal component, the frequency band containing the frequency line with the largest amplitude in the envelope spectrum is the prominent frequency range, and the frequency band is a symmetric interval with the prominent spectral line as the symmetry axis.
Step 3: according to the optimal components obtained in step 2, use SSA to perform adaptive optimization on the parameters of MCKD, and obtain the required parameter combination [L 0 , T 0 , M 0 ] of filter length, impulse signal period, and shift number. The parameter settings in SSA are the same as in step 1. According to the literature [35] and combined with the actual situation, set the optimization range of parameter L to [100, 1000]. After comprehensively considering the accuracy of fault feature extraction and computational efficiency, the optimization range of parameter T is selected as and the optimization range of parameter M is [3,7]. When the fault characteristic frequency is unknown, the frequency band highlighted from step 2 is used to determine the range of the T search in the MCKD parameters. The size of the frequency band should be moderate. The maximum value of the upper boundary is slightly larger than the maximum fault characteristic frequency of the theoretical bearing in the whole transmission system.
Step 4: after setting the parameters in MCKD according to the optimal parameter combination obtained in step 3, use the MCKD method to analyze the optimal component to strengthen the periodic fault pulses in the signal.
Step 5: perform envelope demodulation on the deconvolved signal and obtain the fault characteristic frequency value and its multiplication frequency. The theoretical fault characteristic frequency values of the bearings are compared with the spectral lines with significant peaks in the envelope spectrum to diagnose the fault type.
Through the above steps, effective extraction of rolling bearing fault characteristics and determination of fault types can be achieved.

Simulation Signal Analysis
When measuring the data, due to the influence of factors such as the accuracy of the sensor and the operating environment (temperature, humidity), the collected vibration signal usually contains a certain degree of error. To verify the feasibility and effectiveness of the method proposed in the paper, the impact signal generated by the bearing inner ring failure is simulated using the rolling bearing failure model, and strong Gaussian white noise is added to simulate the faint failure of the actual working condition of the bearing covered by ambient noise. The original simulation signal is shown in Equation (15) [36]: where attenuation coefficient C = 800; resonance frequency f n = 4000 Hz; amplitude A 0 = 0.5; rotation frequency f r = 25 Hz; sampling frequency f s = 12,800 Hz; fault characteristic frequency f fault = 120 Hz; the number of analysis points is 8192; τ k means obeying a small fluctuation with a mean of zero and a standard deviation of 0.5% f r ; n(t) represents the Gaussian white noise component; and the signal-to-noise ratio of the noisecontaining signal is set to −16 dB.
The time-domain waveforms of the original signal and of the simulated signal after adding noise are shown in Figure 2. As can be seen from Figure 2, the original signal contains obvious periodic shock components. After adding Gaussian white noise to the original signal, the periodic shock signal in Figure 2a is completely drowned out by the noise, and the original periodic impact component cannot be distinguished.
The frequency-domain waveform and envelope spectrum of the noisy signal are shown in Figure 3. From Figure 3a, it can be concluded that the spectrum of the noise-containing signal has no regularity and the frequency components are disordered, making it difficult to find fault information. From the envelope spectrum of Figure 3b, it is difficult to find the prominent frequencies, and it is impossible to identify the fault information effectively. In sum- As can be seen from Figure 2, the original signal contains obvious periodic shock components. After adding Gaussian white noise to the original signal, the periodic shock signal in Figure 2a is completely drowned out by the noise, and the original periodic impact component cannot be distinguished.
The frequency-domain waveform and envelope spectrum of the noisy signal are shown in Figure 3. As can be seen from Figure 2, the original signal contains obvious periodic shock components. After adding Gaussian white noise to the original signal, the periodic shock signal in Figure 2a is completely drowned out by the noise, and the original periodic impact component cannot be distinguished.
The frequency-domain waveform and envelope spectrum of the noisy signal are shown in Figure 3. From Figure 3a, it can be concluded that the spectrum of the noise-containing signal has no regularity and the frequency components are disordered, making it difficult to find fault information. From the envelope spectrum of Figure 3b, it is difficult to find the prominent frequencies, and it is impossible to identify the fault information effectively. In sum- From Figure 3a, it can be concluded that the spectrum of the noise-containing signal has no regularity and the frequency components are disordered, making it difficult to find fault information. From the envelope spectrum of Figure 3b, it is difficult to find the prominent frequencies, and it is impossible to identify the fault information effectively. In summary, the traditional spectrum analysis method has failed.
The fault features are extracted using the SSA-VMD-MCKD method proposed in this paper. The curve of the variation of the value of the fitness function with the number of iterations is obtained by adaptive optimization seeking of the parameters of VMD through SSA as shown in Figure 4. iterations is obtained by adaptive optimization seeking of the parameters of VMD through SSA as shown in Figure 4.  As can be seen from Figure 5, the VMD algorithm is used to decompose the simulated signal to obtain a total of nine IMF components, with the center frequencies of IMF1~IMF9 The SSA-optimized VMD converges in the third generation and the value of the fitness function searched is 5.5 × 10 −3 , which is a dimensionless index. It can be concluded that the optimal combination of parameters [α 0 , K 0 ] for VMD are α 0 = 1423 and K 0 = 9, respectively. After setting the parameters in VMD according to the results obtained from the optimization search, the IMF components obtained by using VMD to process the signal containing noise are shown in Figure 5.
As can be seen from Figure 5, the VMD algorithm is used to decompose the simulated signal to obtain a total of nine IMF components, with the center frequencies of IMF1~IMF9 gradually increasing. From the waveform and frequency change of IMF components, it can be inferred that the forward component is the real signal component, and the higher frequency at the back may be the high-frequency noise component in the signal. The envelope spectral peak factor values of each IMF component in the above figure are calculated according to Equation (14), and the results obtained are shown in Figure 6.
The component with the largest value of the peak factor of the envelope spectrum is selected as the desired optimal component. The magnitude plot of each component in Figure 6 shows that the value of IMF4 is the largest by comparison; therefore, IMF4 is selected as the optimal component. Using the envelope spectrum to analyze this component, the results obtained are shown in Figure 7.
It can be observed that the amplitude of the corresponding spectral line is largest when the horizontal coordinate is taken as 120 Hz, and it can be inferred that this frequency is the most likely frequency close to the fault characteristics. However, in the bearing fault signal envelope spectrum, only f fault and its multiples can be observed to ensure effective extraction of fault features, and the 120 Hz multiples component does not appear in the above figure; therefore, it cannot be determined that 120 Hz is the fault frequency. In order to accurately extract the fault characteristics, further analysis should be performed. According to the prominent spectral line with the largest amplitude in the above figure, the frequency band containing 120 Hz is selected as the prominent frequency range, and the frequency band [90, 150] is selected as the frequency range for solving T after comprehensive consideration.
The SSA-optimized VMD converges in the third generation and the value of the fitness function searched is 3 5.5 10   , which is a dimensionless index. It can be concluded that the optimal combination of parameters   0 0 , K  for VMD are 0 1423   and 0 9 K  , respectively. After setting the parameters in VMD according to the results obtained from the optimization search, the IMF components obtained by using VMD to process the signal containing noise are shown in Figure 5. As can be seen from Figure 5, the VMD algorithm is used to decompose the simulated signal to obtain a total of nine IMF components, with the center frequencies of IMF1~IMF9 gradually increasing. From the waveform and frequency change of IMF components, it can be inferred that the forward component is the real signal component, and the higher  (14), and the results obtained are shown in Figure 6. The component with the largest value of the peak factor of the envelope spectrum is selected as the desired optimal component. The magnitude plot of each component in Figure 6 shows that the value of IMF4 is the largest by comparison; therefore, IMF4 is selected as the optimal component. Using the envelope spectrum to analyze this component, the results obtained are shown in Figure 7.  (14), and the results obtained are shown in Figure 6. The component with the largest value of the peak factor of the envelope spectrum is selected as the desired optimal component. The magnitude plot of each component in Figure 6 shows that the value of IMF4 is the largest by comparison; therefore, IMF4 is selected as the optimal component. Using the envelope spectrum to analyze this component, the results obtained are shown in Figure 7. It can be observed that the amplitude of the corresponding spectral line is largest when the horizontal coordinate is taken as 120 Hz, and it can be inferred that this frequency is the most likely frequency close to the fault characteristics. However, in the bearing fault signal envelope spectrum, only f and its multiples can be observed to ensure To highlight the superiority of parameter-optimized VMD, the noise-containing signal is processed by CEEMDAN to obtain 15 IMF components. The peak factor value of the envelope spectrum of each component is calculated separately, the IMF12 with the largest index is taken as the optimal component for which the envelope demodulation analysis is performed, and the results obtained are shown in Figure 8. To highlight the superiority of parameter-optimized VMD, the noise-containing signal is processed by CEEMDAN to obtain 15 IMF components. The peak factor value of the envelope spectrum of each component is calculated separately, the IMF12 with the largest index is taken as the optimal component for which the envelope demodulation analysis is performed, and the results obtained are shown in Figure 8. As can be seen from Figure 8, the highlighted frequency range does not include the fault characteristic frequencies. Therefore, the superiority of using SSA to optimize VMD parameters is verified.
The parameters of MCKD are adaptively searched for using SSA Figure 9 shows the curve of the value of the fitness function of the optimized MCKD algorithm with the number of iterations of population evolution; the vertical coordinates are dimensionless indicators. As can be seen from Figure 8, the highlighted frequency range does not include the fault characteristic frequencies. Therefore, the superiority of using SSA to optimize VMD parameters is verified.
The parameters of MCKD are adaptively searched for using SSA Figure 9 shows the curve of the value of the fitness function of the optimized MCKD algorithm with the number of iterations of population evolution; the vertical coordinates are dimensionless indicators. To highlight the superiority of parameter-optimized VMD, the noise-containing signal is processed by CEEMDAN to obtain 15 IMF components. The peak factor value of the envelope spectrum of each component is calculated separately, the IMF12 with the largest index is taken as the optimal component for which the envelope demodulation analysis is performed, and the results obtained are shown in Figure 8. As can be seen from Figure 8, the highlighted frequency range does not include the fault characteristic frequencies. Therefore, the superiority of using SSA to optimize VMD parameters is verified.
The parameters of MCKD are adaptively searched for using SSA Figure 9 shows the curve of the value of the fitness function of the optimized MCKD algorithm with the number of iterations of population evolution; the vertical coordinates are dimensionless indicators.  From the above figure, it can be seen that SSA converges in the ninth generation and the fitness function takes the minimum value 5 × 10 −3 . The optimal combination of parameters is obtained as L 0 = 740, T 0 = 99 and M 0 = 5. After setting the parameters in MCKD according to the optimal parameter combination obtained, use MCKD to analyze the optimal component, and obtain the processed time-domain plot and envelope spectrum, as shown in Figure 10. The optimal component is processed by MCKD to highlight continuous shock pulses that are drowned by noise and to improve the correlation kurtosis value of the original signal. The extraction of fault features can be accomplished by envelope demodulation. After completing the above steps, the fault shock component is clearly observable, and the signal shock component is significantly increased in the time-domain waveform of Figure  10a. The spectral lines of both the fault characteristic frequency and its multiples are clearly visible in the deconvoluted envelope spectrum of Figure 10b, and the fault type can thus be identified as an inner-loop fault, indicating that the characteristic frequency is accurately extracted, thus verifying the correctness of the present method.
The reliability of the results obtained by using SSA to optimize the MCKD parameters was verified by modifying the three parameters of MCKD. These were changed in accordance with 10% float from the original parameter value, and the filter length, impulse signal period, and shift number were changed to 666, 90, and 4, respectively. The results are shown in Figure 11. Figure 11. The envelope spectrum after modifying the parameters for MCKD processing.
It is more difficult to observe the fault characteristic frequency in Figure 11. There are many interferences spectral lines around it and the multiplication frequency is not obvi- The optimal component is processed by MCKD to highlight continuous shock pulses that are drowned by noise and to improve the correlation kurtosis value of the original signal. The extraction of fault features can be accomplished by envelope demodulation. After completing the above steps, the fault shock component is clearly observable, and the signal shock component is significantly increased in the time-domain waveform of Figure 10a. The spectral lines of both the fault characteristic frequency and its multiples are clearly visible in the deconvoluted envelope spectrum of Figure 10b, and the fault type can thus be identified as an inner-loop fault, indicating that the characteristic frequency is accurately extracted, thus verifying the correctness of the present method.
The reliability of the results obtained by using SSA to optimize the MCKD parameters was verified by modifying the three parameters of MCKD. These were changed in accordance with 10% float from the original parameter value, and the filter length, impulse signal period, and shift number were changed to 666, 90, and 4, respectively. The results are shown in Figure 11. The optimal component is processed by MCKD to highlight continuous shock pulses that are drowned by noise and to improve the correlation kurtosis value of the original signal. The extraction of fault features can be accomplished by envelope demodulation. After completing the above steps, the fault shock component is clearly observable, and the signal shock component is significantly increased in the time-domain waveform of Figure  10a. The spectral lines of both the fault characteristic frequency and its multiples are clearly visible in the deconvoluted envelope spectrum of Figure 10b, and the fault type can thus be identified as an inner-loop fault, indicating that the characteristic frequency is accurately extracted, thus verifying the correctness of the present method.
The reliability of the results obtained by using SSA to optimize the MCKD parameters was verified by modifying the three parameters of MCKD. These were changed in accordance with 10% float from the original parameter value, and the filter length, impulse signal period, and shift number were changed to 666, 90, and 4, respectively. The results are shown in Figure 11. Figure 11. The envelope spectrum after modifying the parameters for MCKD processing.
It is more difficult to observe the fault characteristic frequency in Figure 11. There are many interferences spectral lines around it and the multiplication frequency is not obvi-  Figure 11. The envelope spectrum after modifying the parameters for MCKD processing.
It is more difficult to observe the fault characteristic frequency in Figure 11. There are many interferences spectral lines around it and the multiplication frequency is not obvious enough, which makes it impossible to extract the fault characteristics effectively, indicating that the results of optimizing MCKD parameters using SSA in this paper are reliable.
To verify the necessity of combining the optimized VMD algorithm and the optimized MCKD algorithm proposed in this paper, the simulated signal is analyzed without the optimized VMD preprocessing, based on the SSA optimized MCKD for the simulated signal. After optimization, the filter length, impulse signal period, and shift number are 848, 107, and 4, respectively, and the results obtained are shown in Figure 12.
Electronics 2022, 11, x FOR PEER REVIEW 15 of 26 ous enough, which makes it impossible to extract the fault characteristics effectively, indicating that the results of optimizing MCKD parameters using SSA in this paper are reliable.
To verify the necessity of combining the optimized VMD algorithm and the optimized MCKD algorithm proposed in this paper, the simulated signal is analyzed without the optimized VMD preprocessing, based on the SSA optimized MCKD for the simulated signal. After optimization, the filter length, impulse signal period, and shift number are 848, 107, and 4, respectively, and the results obtained are shown in Figure 12. Obviously, the direct use of MCKD with parameters optimized by SSA is not satisfactory for processing the rolling bearing fault signal, and there are more interference components in the envelope spectrum, which makes it difficult to clearly distinguish the fault feature frequencies and their multiples, resulting in difficulty in extracting the fault features. This proves the effectiveness of the fault feature extraction method proposed in this paper, and the necessity of using the parameter-seeking VMD to process the signal before MCKD analysis.

Measured Signal Analysis
Based on the successful verification of the simulated signals, the feasibility and effectiveness of the SSA-VMD-MCKD rolling bearing fault feature extraction method in practical engineering applications was further verified using the CWRU data set and the XJTU-SY data set [37].

Analysis of Vibration Data of CWRU Rolling Bearings
The effectiveness of applying the SSA-VMD-MCKD method to rolling bearing fault feature extraction was further validated using a publicly available bearing vibration signal data set from CWRU. The experiment uses deep groove ball bearings, model SKF6205. The failure of the bearing is single point damage with EDM, and the vibration acceleration signal of the bearing is collected using an acceleration sensor [38]. The specific data used is the inner ring fault vibration signal with motor load of 0, approximate speed of 1797 r/min, fault diameter of 0.1778 mm, rotation frequency 1797 60 29.95Hz r f   , and sampling frequency 12000Hz s f  . According to the parameters provided by the data set, the inner ring fault characteristic frequency obtained is 5.4152 162.1852Hz i r f f   [39]. The first 4096 sampling points of the data set are selected as the fault signal test data; the time-domain waveform and envelope spectrum of the experimental signal are shown in Figure 13. Obviously, the direct use of MCKD with parameters optimized by SSA is not satisfactory for processing the rolling bearing fault signal, and there are more interference components in the envelope spectrum, which makes it difficult to clearly distinguish the fault feature frequencies and their multiples, resulting in difficulty in extracting the fault features. This proves the effectiveness of the fault feature extraction method proposed in this paper, and the necessity of using the parameter-seeking VMD to process the signal before MCKD analysis.

Measured Signal Analysis
Based on the successful verification of the simulated signals, the feasibility and effectiveness of the SSA-VMD-MCKD rolling bearing fault feature extraction method in practical engineering applications was further verified using the CWRU data set and the XJTU-SY data set [37].

Analysis of Vibration Data of CWRU Rolling Bearings
The effectiveness of applying the SSA-VMD-MCKD method to rolling bearing fault feature extraction was further validated using a publicly available bearing vibration signal data set from CWRU. The experiment uses deep groove ball bearings, model SKF6205. The failure of the bearing is single point damage with EDM, and the vibration acceleration signal of the bearing is collected using an acceleration sensor [38]. The specific data used is the inner ring fault vibration signal with motor load of 0, approximate speed of 1797 r/min, fault diameter of 0.1778 mm, rotation frequency f r = 1797/60 = 29.95 Hz, and sampling frequency f s = 12,000 Hz. According to the parameters provided by the data set, the inner ring fault characteristic frequency obtained is f i = 5.4152 f r = 162.1852 Hz [39].
The first 4096 sampling points of the data set are selected as the fault signal test data; the time-domain waveform and envelope spectrum of the experimental signal are shown in Figure 13. Insignificant periodic shock components can be observed in the time-domain waveform of the fault signal, which is due to the unavoidable noise interference during data acquisition. Some of the shock components have been covered by noise. In addition, further envelope demodulation analysis of the signal, as shown in Figure 13b, did not reveal prominent frequency components.
The SSA-VMD-MCKD method was used to analyze the experimental signals. First, the parameters of VMD are adaptively optimized by SSA. Figure 14 shows the curve of the change in the value of the fitness function with the number of iterations of population evolution. The SSA optimization converges in the third generation, the minimum value of the fitness function is taken to be 3 8.5 10   , and the optimal combination of parameters obtained is 0 644   and 0 8 K  , respectively. After setting the parameters in VMD according to the combination of the parameters obtained from the optimization search, VMD is performed on the fault signal and eight IMF components are obtained, as shown in Figure  15.  Insignificant periodic shock components can be observed in the time-domain waveform of the fault signal, which is due to the unavoidable noise interference during data acquisition. Some of the shock components have been covered by noise. In addition, further envelope demodulation analysis of the signal, as shown in Figure 13b, did not reveal prominent frequency components.
The SSA-VMD-MCKD method was used to analyze the experimental signals. First, the parameters of VMD are adaptively optimized by SSA. Figure 14 shows the curve of the change in the value of the fitness function with the number of iterations of population evolution. Insignificant periodic shock components can be observed in the time-domain waveform of the fault signal, which is due to the unavoidable noise interference during data acquisition. Some of the shock components have been covered by noise. In addition, further envelope demodulation analysis of the signal, as shown in Figure 13b, did not reveal prominent frequency components.
The SSA-VMD-MCKD method was used to analyze the experimental signals. First, the parameters of VMD are adaptively optimized by SSA. Figure 14 shows the curve of the change in the value of the fitness function with the number of iterations of population evolution. The SSA optimization converges in the third generation, the minimum value of the fitness function is taken to be 3 8.5 10   , and the optimal combination of parameters obtained is 0 644   and 0 8 K  , respectively. After setting the parameters in VMD according to the combination of the parameters obtained from the optimization search, VMD is performed on the fault signal and eight IMF components are obtained, as shown in Figure  15.  The SSA optimization converges in the third generation, the minimum value of the fitness function is taken to be 8.5 × 10 −3 , and the optimal combination of parameters obtained is α 0 = 644 and K 0 = 8, respectively. After setting the parameters in VMD according to the combination of the parameters obtained from the optimization search, VMD is performed on the fault signal and eight IMF components are obtained, as shown in Figure 15. The value of the envelope spectral peak factor is calculated for each IMF component in the above figure, and the result obtained is shown in Figure 16. Through the plot of each component amplitude of each component in Figure 16, it can be seen that IMF4 has the largest value of the envelope spectrum peak factor; therefore, IMF4 is selected as the optimal component. Envelope spectrum analysis is performed for IMF4, and the results are shown in Figure 17. The value of the envelope spectral peak factor is calculated for each IMF component in the above figure, and the result obtained is shown in Figure 16. The value of the envelope spectral peak factor is calculated for each IMF component in the above figure, and the result obtained is shown in Figure 16. Through the plot of each component amplitude of each component in Figure 16, it can be seen that IMF4 has the largest value of the envelope spectrum peak factor; therefore, IMF4 is selected as the optimal component. Envelope spectrum analysis is performed for IMF4, and the results are shown in Figure 17. Through the plot of each component amplitude of each component in Figure 16, it can be seen that IMF4 has the largest value of the envelope spectrum peak factor; therefore, IMF4 is selected as the optimal component. Envelope spectrum analysis is performed for IMF4, and the results are shown in Figure 17.  Figure 18. It can be seen from Figure 18 that when the number of iterations reaches 19, the fitness function takes the minimum value 3 9 10   . After setting the parameters in MCKD according to the optimal parameter combination obtained, use MCKD to perform deconvolution analysis on the optimal components and obtain the analyzed time-domain waveform and envelope spectrum, as shown in Figure 19.   Figure 18.  Figure 18. It can be seen from Figure 18 that when the number of iterations reaches 19, the fitness function takes the minimum value 3 9 10   . After setting the parameters in MCKD according to the optimal parameter combination obtained, use MCKD to perform deconvolution analysis on the optimal components and obtain the analyzed time-domain waveform and envelope spectrum, as shown in Figure 19.  It can be seen from Figure 18 that when the number of iterations reaches 19, the fitness function takes the minimum value 9 × 10 −3 . After setting the parameters in MCKD according to the optimal parameter combination obtained, use MCKD to perform deconvolution analysis on the optimal components and obtain the analyzed time-domain waveform and envelope spectrum, as shown in Figure 19.
From the time-domain waveform in Figure 19a, the shock component can be clearly observed. The spectral lines of the fault characteristic frequency and its multiples are clearly visible in the deconvoluted envelope spectrum, and the fault type of the fault signal can thus be identified as an inner-loop fault, indicating that the characteristic frequency is accurately extracted, thus verifying the effectiveness of the method proposed in this paper. From the time-domain waveform in Figure 19a, the shock component can be clearly observed. The spectral lines of the fault characteristic frequency and its multiples are clearly visible in the deconvoluted envelope spectrum, and the fault type of the fault signal can thus be identified as an inner-loop fault, indicating that the characteristic frequency is accurately extracted, thus verifying the effectiveness of the method proposed in this paper.
To verify the necessity of the method proposed in this paper to process the fault signal using VMD before MCKD analysis, the fault signal is directly analyzed using MCKD with parameters optimized by SSA. The results obtained are shown in Figure 20. Although the fault signal can be processed directly using MCKD with parameters optimized by SSA, and the fault feature frequencies and their multipliers can be obtained, the results are obviously not as clear and those obtained using the method proposed in this paper. This demonstrates the effectiveness of the fault feature extraction method proposed in this paper and the necessity of using VMD after parameter optimization to process the signal before MCKD analysis. To verify the necessity of the method proposed in this paper to process the fault signal using VMD before MCKD analysis, the fault signal is directly analyzed using MCKD with parameters optimized by SSA. The results obtained are shown in Figure 20. From the time-domain waveform in Figure 19a, the shock component can be clearly observed. The spectral lines of the fault characteristic frequency and its multiples are clearly visible in the deconvoluted envelope spectrum, and the fault type of the fault signal can thus be identified as an inner-loop fault, indicating that the characteristic frequency is accurately extracted, thus verifying the effectiveness of the method proposed in this paper.
To verify the necessity of the method proposed in this paper to process the fault signal using VMD before MCKD analysis, the fault signal is directly analyzed using MCKD with parameters optimized by SSA. The results obtained are shown in Figure 20. Although the fault signal can be processed directly using MCKD with parameters optimized by SSA, and the fault feature frequencies and their multipliers can be obtained, the results are obviously not as clear and those obtained using the method proposed in this paper. This demonstrates the effectiveness of the fault feature extraction method proposed in this paper and the necessity of using VMD after parameter optimization to process the signal before MCKD analysis. Although the fault signal can be processed directly using MCKD with parameters optimized by SSA, and the fault feature frequencies and their multipliers can be obtained, the results are obviously not as clear and those obtained using the method proposed in this paper. This demonstrates the effectiveness of the fault feature extraction method proposed in this paper and the necessity of using VMD after parameter optimization to process the signal before MCKD analysis.

XJTU-SY Bearing Life Cycle Data Analysis
The XJTU-SY bearing life cycle data compiled by Dr. Wang Biao of Lei Yaguo's research group at Xi'an Jiao tong University are used to further verify the effectiveness of the method proposed in this paper [40]. The bearing accelerated life test platform is shown in Figure 21. The test bearing is LDK UER204, the number of balls is eight, and the contact angle is 0 • . An acceleration sensor is used to obtain the bearing vibration signal. The rotating speed is 2100 r/min, sampling frequency f s = 25,600 Hz, rotation frequency f r = 2100/60 = 35 Hz, radial force is 12 kN, outer ring fault characteristic frequency f o = 107.91 Hz, and a fault signal of a length of 4096 was selected for the analysis [41].

XJTU-SY Bearing Life Cycle Data Analysis
The XJTU-SY bearing life cycle data compiled by Dr. Wang Biao of Lei Yaguo's research group at Xi'an Jiao tong University are used to further verify the effectiveness of the method proposed in this paper [40]. The bearing accelerated life test platform is shown in Figure 21. The test bearing is LDK UER204, the number of balls is eight, and the contact angle is 0  . An acceleration sensor is used to obtain the bearing vibration signal. The rotating speed is 2100 r/min, sampling frequency 25  The time-domain waveform and envelope spectrum of the experimental signal are shown in Figure 22.
The XJTU-SY bearing life cycle data compiled by Dr. Wang Biao of Lei Yaguo's research group at Xi'an Jiao tong University are used to further verify the effectiveness of the method proposed in this paper [40]. The bearing accelerated life test platform is shown in Figure 21. The test bearing is LDK UER204, the number of balls is eight, and the contact angle is 0  . An acceleration sensor is used to obtain the bearing vibration signal.  The time-domain waveform hardly observes any shock component, and there is no pattern after the envelope demodulation analysis of the signal.
The experimental signal was analyzed using the method proposed in this paper. The parameters of VMD were optimized by SSA, and the variation curve of the fitness function value with the number of iterations is shown in Figure 23.
The time-domain waveform hardly observes any shock component, and there is no pattern after the envelope demodulation analysis of the signal.
The experimental signal was analyzed using the method proposed in this paper. The parameters of VMD were optimized by SSA, and the variation curve of the fitness function value with the number of iterations is shown in Figure 23. After obtaining the optimal parameter combinations of the number of components and the penalty factor in the VMD method as 0 8 K  and 0 904   , respectively, perform VMD processing on the fault signal, and obtain eight IMF components, as shown in Figure 24. Calculate the peak factor value of the envelope spectrum of each IMF component in the above figure; the result is shown in Figure 25.  Figure 23. The fitness function varies with the number of iterations.
After obtaining the optimal parameter combinations of the number of components and the penalty factor in the VMD method as K 0 = 8 and α 0 = 904, respectively, perform VMD processing on the fault signal, and obtain eight IMF components, as shown in Figure 24.
The time-domain waveform hardly observes any shock component, and there is no pattern after the envelope demodulation analysis of the signal.
The experimental signal was analyzed using the method proposed in this paper. The parameters of VMD were optimized by SSA, and the variation curve of the fitness function value with the number of iterations is shown in Figure 23. After obtaining the optimal parameter combinations of the number of components and the penalty factor in the VMD method as 0 8 K  and 0 904   , respectively, perform VMD processing on the fault signal, and obtain eight IMF components, as shown in Figure 24.    Figure 24. VMD decomposition result of fault signal.
Calculate the peak factor value of the envelope spectrum of each IMF component in the above figure; the result is shown in Figure 25.
It can be seen that IMF1 is the optimal component. In addition, the range of MCKD optimization parameters T is calculated as [207,267]. SSA is used to adaptively optimize the parameters of MCKD, and the required filter length, impulse signal period, and shift number parameters are combined as L 0 = 1111, T 0 = 229, and M 0 = 4. The curve of the fitness function value changing with the number of iterations is shown in Figure 26. It can be seen that IMF1 is the optimal component. In addition, the range of MCKD optimization parameters T is calculated as [207,267]. SSA is used to adaptively optimize the parameters of MCKD, and the required filter length, impulse signal period, and shift number parameters are combined as 0 1111 L  , 0 229 T  , and 0 4 M  . The curve of the fitness function value changing with the number of iterations is shown in Figure 26. After setting the parameters in MCKD according to the optimal parameter combination obtained, the time-domain waveform and envelope spectrum obtained by deconvolution analysis of the optimal components, using MCKD, are shown in Figure 27.   It can be seen that IMF1 is the optimal component. In addition, the range of MCKD optimization parameters T is calculated as [207,267]. SSA is used to adaptively optimize the parameters of MCKD, and the required filter length, impulse signal period, and shift number parameters are combined as 0 1111 L  , 0 229 T  , and 0 4 M  . The curve of the fitness function value changing with the number of iterations is shown in Figure 26. After setting the parameters in MCKD according to the optimal parameter combination obtained, the time-domain waveform and envelope spectrum obtained by deconvolution analysis of the optimal components, using MCKD, are shown in Figure 27.  After setting the parameters in MCKD according to the optimal parameter combination obtained, the time-domain waveform and envelope spectrum obtained by deconvolution analysis of the optimal components, using MCKD, are shown in Figure 27. In Figure 27a, the signal impact component is significantly increased. In Figure 27b, the spectral lines of the fault feature frequency and its multiples are clearly visible, indicating that the rolling bearing fault features are effectively extracted, and the fault type of the fault signal can be determined as outer ring fault; thus, the effectiveness of the fault In Figure 27a, the signal impact component is significantly increased. In Figure 27b, the spectral lines of the fault feature frequency and its multiples are clearly visible, indicating that the rolling bearing fault features are effectively extracted, and the fault type of the fault signal can be determined as outer ring fault; thus, the effectiveness of the fault feature extraction method proposed in this paper is further verified. Figure 28 shows the results of the outer ring fault signal by fast spectral kurtosis (FSK) processing [41].
In Figure 27a, the signal impact component is significantly increased. In Figure 27b, the spectral lines of the fault feature frequency and its multiples are clearly visible, indicating that the rolling bearing fault features are effectively extracted, and the fault type of the fault signal can be determined as outer ring fault; thus, the effectiveness of the fault feature extraction method proposed in this paper is further verified. Figure 28 shows the results of the outer ring fault signal by fast spectral kurtosis (FSK) processing [41]. The bandwidth of the feature signal with the largest kurtosis is 200 Hz and the center frequency is 7300 Hz. The bandwidth is narrow, the noise interference is large, only the transient frequency is found in the filtered signal envelope spectrum, and the fault feature frequency and its multiples cannot be observed, which is not effective. The bandwidth of the feature signal with the largest kurtosis is 200 Hz and the center frequency is 7300 Hz. The bandwidth is narrow, the noise interference is large, only the transient frequency is found in the filtered signal envelope spectrum, and the fault feature frequency and its multiples cannot be observed, which is not effective.

Conclusions
This paper takes rolling bearings as the engineering research background and conducts a study on their fault feature extraction methods. Because the existing fault feature extraction methods usually have poor self-adaptation and an unsatisfactory fault feature extraction effect, a more suitable fault feature extraction method for rolling bearings is proposed: SSA-VMD-MCKD. Simulated signals and measured data verify the effectiveness of the proposed method. The following conclusions can be drawn: (1) The SSA-based improved VMD method is capable of adaptively searching to obtain the decomposition modal number and penalty factor, avoiding the interference of human subjective factors on the selection of VMD parameters, achieving effective suppression of noise components, and highlighting transient shocks. (2) The improved MCKD method based on SSA can adaptively search for the optimal combination of filter length, deconvolution signal period, and shift number. The interference of human subjective factors on the selection of MCKD parameters is avoided, and the optimal deconvolution of fault signals is achieved. The SSA-VMD-MCKD method proposed in this paper can provide theoretical and technical support for rolling bearing fault features extraction and has good prospects for engineering applications.
The study of the fault feature extraction method proposed in this paper is only validated with the simulation data and the preset fault data of the test bench, and its generality needs to be further studied in conjunction with the actual equipment. Meanwhile, although the method proposed in this paper has a better fault feature extraction effect compared with other methods, there is a deficiency of long parameter search time, and more advanced intelligent optimization algorithms can be used to improve the method proposed in this paper.  Data Availability Statement: These data were collected through many experiments and are copyrighted by the laboratory and therefore cannot be made public, thank you.

Conflicts of Interest:
The authors declare no conflict of interest.