Abstract
To solve the difficulty in correctly identifying a compound fault of rolling bearing, a method combining variational mode decomposition (VMD) and harmonic fusion vector bispectrum (HFVB) is proposed. Firstly, to achieve adaptive decomposition of signals, the characteristic ability of envelope entropy to represent signal sparsity is utilized. By employing the minimum envelope entropy as the fitness function for the sparrow search algorithm (SSA), the decomposition levels and penalty factors of VMD are adaptively determined. Secondly, root-mean-square value is treated as fault feature index to self-adaptively choose from intrinsic mode function (IMF) which can embody fault features of bearings. Thirdly, to further highlight fault features, HFVB is used to blend chosen IMFs. Finally, faults of bearings were recognized with spectrum of signals. To verify the effectiveness of proposed method, an analysis was given to vibration signals in different situations, and SSA-VMD-HFVB was compared with classical method. The results demonstrate that the proposed SSA-VMD-HFVB method facilitates the adaptive decomposition of VMD, enabling the selection and effective integration of sensitive fault component signals. This approach enhances the accuracy of complex fault diagnosing in rolling bearings.
Keywords
Introduction
Rolling bearings are extensively employed in modern mechanical equipment, and their health directly impacts the stability and lifespan of the equipment. When a rolling bearing malfunctions, its vibration signal encapsulates rich diagnostic information. Therefore, determining the operational status and type of fault in bearings based on vibration signals is one of the primary research directions in bearing fault diagnosis. 1 The vibration signal during a bearing fault usually contains considerable noise and numerous frequency components, showcasing characteristics such as non-linearity, non-stationarity, weak periodicity, and a low signal-to-noise ratio.2–4 Signals decompose algorithms, addressing non-linearity and non-stationarity, are a typical class of methods for fault feature extraction including local feature scale decomposition,5,6 intrinsic time scale decomposition, 7 and local mean decomposition. Huang et al. proposed the empirical mode decomposition (EMD) method, which offers stabilization of complex signals and decomposes vibration signals into multiple intrinsic mode functions (IMF). However, this method has inherent issues like severe mode mixing. 8 The variational mode decomposition (VMD) can effectively mitigate mode mixing and end effects, but its decomposition performance is dependent on the number of decomposition layers and the selection of the penalty factor. 9 Addressing VMD parameter determination, the most typical approach integrates a signal evaluation index, which embodies fault information, with swarm intelligence optimization algorithms.10–12 Zhang et al. 10 used bacterial foraging algorithm to optimize VMD, and optimized VMD parameters using minimum dispersion entropy as the fitness function. Liu et al. 11 used the minimum fuzzy entropy and kurtosis reciprocal sum as the objective function of the whale algorithm to determine the VMD parameters. Among the signal evaluation indices, envelope entropy effectively evaluates the sparsity characteristics of signals. When an IMF contains more noise and less feature information, its envelope entropy value tends to be high. 13 Therefore, this paper employs the minimum envelope entropy as the fitness function for the sparrow optimization algorithm to adaptively determine VMD parameters.
After signal decomposition, the accurate selection of fault-sensitive component signals is essential for effective fault feature extraction and is key to accurate fault identification. Vibration signals decomposed by VMD yield several IMFs of different modes, with each containing varying amounts of sensitive fault information. When selecting these fault-sensitive component signals, the decision is often based on signal evaluation indices. 14 Reference 15 introduces a method for weak fault recognition of rolling bearing based on a combination of variational mode decomposition and singular value entropy. This method uses the mean square deviation-duclidean distance index for selecting the IMF rich in fault information. Reference 16 employs the time-domain cross-correlation criterion for component signal filtering and uses the envelope spectrum of reconstruction signal for fault diagnosis. The kurtosis index is sensitive to impact characteristics and can effectively detect periodic impact in signals. It is a commonly used indicator parameter when selecting component signals. 17 However, its sensitivity to noise is noted drawbacks. 18 In signal evaluation metrics, root mean square (RMS) is a dimensional characteristic parameter reflecting the overall energy magnitude of vibration signals. When a rolling bearing develops a fault, it’s usually accompanied by intensified impulse impacts, leading to significant changes of signal amplitude; RMS is particularly sensitive to bearing faults. 19 Mi et al. 20 employed mean, variance, and RMS values to design a novel threshold for anomaly detection. This method exhibits heightened sensitivity to anomalies indicative of bearing degradation, outperforming classic detection methods with a lower error rate. Jaypee 21 utilized local mean decomposition (LMD) to analyze non-stationary and non-linear signals of bearing fault, extracting six statistical indicators (RMS, mean, variance, RMA, ama, and kurtosis) for classifier input.
Based on the selected fault-sensitive component signals, typical signal processing methods for fault feature extraction include: (1) Direct analysis of the component signals 22 ; (2) Weighted reconstruction of component signals based on specific weight coefficients, followed by analysis of the reconstructed signals; (3) Feature extraction after merging component signals using fusion algorithms. Direct analysis of component signals often involves lots of signals, potentially containing many components irrelevant to faults, thereby complicating fault identification. 23 The accurate determination of weight coefficients is a key issue when reconstructing signals based on weight coefficients. 24 The vector bispectrum is a bispectral analysis method based on full vector spectrum information fusion. This method can effectively merge the dual-layer information of rotating machinery, offering a more comprehensive and accurate reflection of non-linear fault characteristics embedded in signals.25,26 Yu et al. 27 tackled the challenge of extracting rubbing fault characteristics by integrating complexity parameters, harmonic fusion vector bi-spectrum (HFVB), and intrinsic time-scale decomposition (ITD). This approach, leveraging the squared demodulation spectrum of reconstructed signals, effectively and accurately delineates rubbing fault features, facilitating precise fault identification.
To effectively discern equipment faults, this study employs the root mean square (RMS) value for its sensitivity to signal impact characteristics, using RMS as a fault feature indicator. This allows for the adaptive selection of IMFs that accurately reflect fault characteristics of rolling bearing. Leveraging the advantage of HFVB in extracting nonlinear features from signals, the two optimally selected IMFs are fused based on HFVB method. Then fault information is extracted through spectral analysis of the fused signals, facilitating the determination of composite fault types of bearings.
Algorithm introduction
VMD based on sparrow search algorithm and envelope entropy
The VMD algorithm possesses favorable sparsity properties. It overcomes the end-effects and mode mixing issues inherent in the EMD method and is grounded in robust mathematical theory. VMD is primarily divided into the construction and resolution of the variational problem
28
:
The alternating direction multiplier method (ADMM) is utilized for successive updates to K IMFs and their central frequencies. The iterative formulas for modal components and central frequencies involve:
VMD is primary iterative process concludes upon fulfilling condition (5), targeting a predefined convergence threshold. This entails continuous updates to {uk}, {wk}, and λ in equations (3) and (4) until condition (5) is satisfied, securing the optimal solution.
The alternating direction multiplier method is employed to alternately update and derive the K IMFs and their corresponding central frequencies. As can be seen from equation (2): both K and α influence the decomposition performance of VMD. The selection of K and α is particularly crucial for the outcome of the VMD. In this paper, based on Reference 15, envelope entropy is used as the fitness function for the sparrow search algorithm to optimize and determine the values of K and α. A detailed introduction to the sparrow search algorithm can be found in Reference 15.
Proposed method
Adaptive selection of sensitive fault component signals based on RMS index
The RMS value is a time-domain index known for its discernibility and robustness. In physics, the RMS value is frequently used to analyze noise, reflecting the severity of the signal’s vibration through the effective value of the physical quantity. Let’s assume the original signal is Adaptive selection of component signals.
Signal fusion based on vector bicoherence
The essence of vector bicoherence is to introduce the information fusion analysis method of full vector spectrum into bicoherence analysis. In this study, the full vector spectrum is used to fuse two sensitive fault component signals selected, so that the nonlinear features contained in the signals are more fully and accurately manifested.
Let the signals
According to the full vector rotary energy spectrum evaluation method, the fused rotary energy of each harmonic is the sum of the energy of the fusion intensity based on the main oscillating vector
For
Diagnosis process of compound faults of bearings
To precisely recognize a compound fault of bearings, the paper takes minimum envelope entropy as the fitness function of sparrow search algorithm
15
to self-adaptively determine the number layers and penalty factor of VMD. RMS and vector bispectrum were used to self-adaptively choose and combine sensitive fault component signals. Features of compound faults of bearings were extracted and faults were identified based on the frequency spectrum of blended signals. The specific flow of proposed method is shown in Figure 2. Process diagram.
Detail: Step 1: Initialize the number of populations of sparrow algorithm; maximum iteration parameter is [O, M] (O = 10, M = 10); determine the optimization range of VMD parameter [K, α], (k Step 2: Initialize the position vector of sparrow within value range and subject original signals to VMD according to position vector; with minimum envelope entropy as the fitness function of sparrow algorithm, self-adaptively determine [ K, α ]; Step 3: Substitute self-adaptively determined [ K, α ] into VMD and self-adaptively decompose vibration signals to get k IMFs; Step 4: Calculate self-correlation functions of IMFs and normalize each correlation function; Step 5: Calculate the RMS of correlation function after normalization; Step 6: Choose the IMFs as sensitive fault components, which correspond to 2 maximum RMS; Step 7: Blend 2 IMFs with the vector bispectrum; Step 8: Take frequency spectrum analysis of blended signals to extract fault feature frequencies of bearings.
To verify the accuracy and advantages of proposed method, a comparative analysis has been made among the method of directly analyzing component signals, common method (SSA-VMD-weighted reconstruction) and the proposed method (SSA-VMD-HFVB). To maintain the consistency of comparison, the control method is based on the same data, decomposition and denoising algorithms used by the proposed method.
Verification and analysis
Introduction to bearing failures
The data in this article is sourced from the rotor-rolling bearing experimental rig shown in Figure 3. This tester consists of a rotating shaft, gearbox, rolling bearings, adjustable-speed motor, rotor disk, and bearing mount. The tester can simulate complex failures in rolling bearings. During the experiment, vibration acceleration sensors are placed in the vertical (CH3) and horizontal (CH2, CH4) directions of the bearing mount; an eddy current sensor is used to measure the rotational speed (via the CH1 channel). The bearings are damaged using spark erosion, and the geometric parameters of the faulty bearings are shown in Table 1. The types of complex bearing failures addressed in this article are illustrated in Figure 4. The bearing failure types corresponding to Figure 4(a)–(c) are, respectively: combined inner ring + outer ring + rolling element failure, combined outer ring + rolling element failure, and combined inner ring + rolling element failure. Rotor-rolling bearing experimental rig. Geometric dimensions and parameters of faulty bearings. Compound faults of rolling bearings. (a) Inner ring + outer ring + rolling element combined failure, (b) Outer ring + rolling element combined failure, (c) Inner ring + rolling element combined failure.

Bearing fault characteristic frequencies
In the paper,
Rotational frequency:
Inner ring fault characteristic frequency:
Outer ring fault characteristic frequency:
Rolling element fault characteristic frequency:
Cage fault characteristic frequency:
Bearing fault data
Bearing characteristic frequency.
Validation and analysis of proposed method - Case 1
Initially, as an example, a bearing with a compound fault in the inner ring, outer ring, and rolling element is selected randomly, as detected by the horizontal sensor (CH2) (Case 1 in Table 2). The results of analyzing this fault using proposed method are shown in Figure 5. Figure 5(a1) and (a2) displays the time domain and its spectrum of the fault signal. Using envelope entropy as the fitness function for the sparrow algorithm, the VMD parameters are adaptively determined, as shown in Figure 5(a3). From Figure 5(a3), it can be observed that the minimal envelope entropy value of 12.65 appears after 5 iterations, determining the parameter. Then the VMD is applied to Figure 5(a1) based on this parameter, with results presented in Figure 5(a4). Figure 5(a5) shows the normalized time domain graph of the autocorrelation function for each component signal from Figure 5(a4). The RMS values of the signals from Figure 5(a5) are depicted in Figure 5(a6). From Figure 5(a6), it’s evident that the RMS values corresponding to IMF2 and IMF3 are the highest. Thus, the adaptively chosen sensitive fault component signals are IMF2 and IMF3. The time domains of the two signals obtained from the bispectrum fusion of the selected IMF2 and IMF3 are shown in Figure 5(a7) and (a8), while their spectrum are in Figure 5(a9) and a(10). Result of proposed methods -case 1. (a1) Fault signal. (a2) Spectrum of (a1). (a3) Iterative graph of envelope entropy. (a4) IMFs of (a2) after VMD. (a5) Normalized autocorrelation function of (a4). (a6) Root mean square value of (a5). (a7) Time domain of fusion signals 1. (a8) Time domain of fusion signals 2. (a9) Spectrum of (a7). (a10) Spectrum of (a8).
Analyzing Figure 5(a9) and (a10), the spectrum of the signal obtained using proposed method exhibits the following characteristics: (1) A prominent frequency at 854.5 Hz, which corresponds to eight times the inner ring fault characteristic frequency subtract twice the rotational frequency: ((854.5 + 2f
r
)/8 = 113.19). (2) Prominent frequency components at 920.4 Hz and 2954 Hz, corresponding to 14 times (920.4/14 = 65.74) and 45 times of outer ring fault characteristics frequency (2954/45 = 65.64), respectively. (3) A significant frequency at 971.7 Hz, corresponding to 22 times the rolling element fault characteristic frequency (971.7/22 = 44.17). (4) Noise components have been effectively suppressed and make the fault characteristic frequencies more pronounced.
From points (1)-(4), it can be inferred that proposed method, while effectively suppressing noise, extract characteristic frequencies that directly correspond to the fault types of bearing, thereby accurately diagnosing compound faults of bearing. Additionally, while Figure 5(a9) and (a10) have amplitude differences, their characteristic frequencies are very similar. As such, only one typical fused signal will be analyzed in subsequent discussions.
Comparative validation and analysis under other schemes - Case 1
To further validate the superiority of the proposed method, two typical comparative methods are adopted for fault feature extraction of bearing. For comparative validation, the data in this section is consistent with Case 1.
Comparison scheme 1-direct analysis of component signals
The obtained component signals’ autocorrelation function was directly analyzed (the difference from proposed method being the omission of bispectrum fusion). The results are presented in Figure 6. Figure 6 displays the spectrum of each IMF from Figure 5(a5). Spectrum of (a5)- contrast schemes 1-case 1.
From the analysis of Figure 6, the following observations can be made regarding the spectrum of the obtained signal when the autocorrelation functions of the component signals are directly analyzed: (1) A frequency component at 229.5 Hz is present, which corresponds to the difference between five times the rolling element fault characteristic frequency and once fault characteristic frequency of cage, calculated as ((229.5-fc)/5 = 44.02). There’s a prominent frequency component at 290.5 Hz, which corresponds to the sum of six times the rolling element fault characteristic frequency and once the rotational frequency, given by ((290.5-f
r
)/6 = 44.15). Additionally, there’s a frequency component at 971.7 Hz, which corresponds to 22 times the rolling element fault characteristic frequency (971.7/22 = 44.17). (2) Frequency components at 920.4 Hz and 2954 Hz are observed. These frequencies respectively correspond to 14 times (920.4/14) = 65.74 and 45 times (2954/45 = 65.64) of outer ring fault characteristic frequency. (3) No frequency components corresponding to the inner ring fault characteristic frequency were detected, leading to an inaccurate determination of the combined fault types of bearing. (4) Compared to the proposed method, there are more noise components present.
In summary, by directly analyzing the component signals derived from the original signal (without bispectrum fusion), even though bearing fault detection is achieved, the extracted fault characteristic frequencies of bearings are not comprehensive. The determination of combined fault types of bearings is not accurate, and there are relatively more noise components.
Comparison scheme 2-weighted reconstruction analysis
In this section, we consider the classical method of reconstructing component signals based on weight coefficients and extracting fault characteristics from the resulting reconstructed signal as the second comparison scheme. Due to space constraints, we only analyze the method of weighted reconstruction of the obtained component signals based on equal weights. The results are depicted in Figure 7. Figure 7(a) and 7(b) show the time domain and its spectrum of the reconstructed signal obtained from IMF2 and IMF3 based on Figure 5(a4). Result graph of comparison Scheme 2. (a) Time domain plot of reconstructed signals. (b) Spectrum of reconstructing the signals directly.
Analyzing Figure 7(b) and comparing it with Figure 5 (the method proposed in this paper), we can identify the following typical frequency components in Figure 7(b): (1) 1809 Hz, corresponding to 16 times characteristic frequency of inner ring fault of the bearing (1809/16 = 113.06). (2) 1841 Hz, representing 28 times characteristic frequency of outer ring fault of the bearing (1841/28 = 65.75). (3) 1943 Hz, corresponding to 44 times characteristic frequency of rolling element fault of the bearing (1943/44 = 44.15).
It’s evident that the comparison scheme based on weight coefficient reconstruction of component signals can also extract fault characteristic frequencies of bearings. However, compared to the method proposed in this paper, there is a strong noise component, and the characteristic frequencies are relatively less prominent. This is primarily because the bispectrum analysis method, while ensuring that the nonlinear features contained in the signal are more comprehensively and accurately reflected, significantly reduces the influence of noise on the signal.
Effectiveness analysis of the proposed method under different situations-case 2-4
In this section, we examine the sensitivity of the method proposed in this study to installation positions of sensors and to composite fault types of bearings.
Sensitivity analysis of installation positions of sensors - Case 2
Here, we analyze the sensitivity of the proposed method to the installation positions of sensors. We select the typical signals obtained from the vertically oriented sensor (CH3) that were collected simultaneously with the Case 1. Since the data was collected at the same time as in Case 1, the fault types of bearing, rotation frequencies, and characteristic frequencies of bearing faults are consistent with those in Case 1.
Figure 8(b1) and 8(b2) display the time domain and its spectrum of the fault signal, respectively. Figure 8(b3) represents the iteration curve of the fitness of the sparrow optimization algorithm. From Figure 8(b3), we can see that the minimum envelope entropy value of 12.68 appeared after six iterations, thus [K, α] = [5,3000]. By substituting [K, α] into the VMD for fault signal decomposition, the results are shown in Figure 8(b4). Figure 8(b5) presents the normalized autocorrelation function of Figure 8(b4). Figure 8(b6) illustrates the root mean square value of Figure 8(b5). From Figure 8(b6), we can observe that the RMS values of IMF3 and IMF5 are the largest, and thus IMF3 and IMF5 are selected as the optimal component signals. By integrating IMF3 and IMF5 based on the bispectrum, the time domain and spectrum of the resultant signal are illustrated in Figure 8(b7) and (b8), respectively. Results of the proposed method in Case 2. (b1) Fault signal. (b2) Spectrum of (b1). (b3) Iterative graph of envelope entropy. (b4) IMFs of (b2) after VMD. (b5) Normalized autocorrelation function of (b4). (b6) Root mean square value of (b5). (b7) Time domain plot of fusion signals. (b8) Spectrum of (b7).
From the analysis of Figure 8(b8), we can identify the following prominent frequency components: (1) A frequency of 1696 Hz, which corresponds to 15 times characteristic frequency of the inner ring fault (1696/15 = 113.07). (2) A frequency of 1777 Hz, corresponding to 40 times characteristic frequency of the rolling element fault (1777/40 = 44.425). (3) A frequency of 3209 Hz, which represents 49 times characteristic frequency of the outer ring fault (3209/49 = 65.489).
It’s evident that, although the data is sourced from accelerometers installed in different positions, the method proposed is still able to accurately detect the fault characteristic frequencies that fully correspond with the fault types of bearing, thereby achieving precise identification of the composite bearing faults. This implies that the proposed method is not sensitive to the sensor’s installation position.
Sensitivity to composite fault types - cases 3-4
To verify the effectiveness of the method proposed in this study (SSA-VMD-HFVB) under different composite fault types, this section analyzes data under different composite fault types of bearings. Randomly chosen data represent compound faults involving the outer ring combined with the rolling element, and the inner ring combined with the rolling element. Given that our method is not sensitive to installation positions of sensors, we arbitrarily use data from the horizontal (CH2) direction sensor for this analysis. The characteristic frequencies for the outer ring, inner ring, rolling element, and cage are calculated according to Equation 15–19, as shown in cases 3-4 of Table 2.
Results for the combined outer ring and rolling element fault are presented in Figure 9(c1)–(c7); results for the combined inner ring and rolling element fault are illustrated in Figure 9(d1)–(d7). Figure 9(c1)–(c2) and 9(d1) and (d2) show the time domain and spectrum of the fault signals, respectively. Figure 9(c3) and (d3) display the iteration curves of the fitness function of the sparrow optimization algorithm. From Figure 9(c3), it can be seen that after 5 iterations of SSA, a minimum envelope entropy value of 12.68 emerges, resulting in [K,α] = [4,381.96]. From Figure 9(d3), after 6 iterations of SSA, a minimum local envelope entropy value of 12.48 is detected, yielding optimal parameters [K,α] = [7,961.89]. These parameter sets are then fed into VMD for signal decomposition, as shown in Figure 9(c4) and (d4). The normalized autocorrelation functions for Figure 9(c4) and (d4) are depicted in Figure 9(c5) and (d5), respectively. The root mean square values for each IMF component are presented in Figure 9(c6) and (d6). From Figure 9(c6), it’s observed that IMF2 and IMF3 have the highest root mean square values. Therefore, when inner ring and rolling element faults occur simultaneously, IMF3 and IMF5 are selected as the optimal signals and fused based on the bispectrum. In Figure 9(d6), IMF3 and IMF7 have the largest root mean square values, so IMF3 and IMF7 are selected for fusion. The resulting fused spectra are illustrated in Figure 9(c7) and (d7). Results of the proposed method in case 3-4. (c1) Fault signal. (d1) Fault signal. (c2) Spectrum of (c1). (d2) Spectrum of (d1). (c3) Iterative graph of envelope entropy. (d3) Iterative graph of envelope entropy. (c4) IMFs of (c2) after VMD. (d4) IMFs of (d2) after VMD. (c5) Normalized autocorrelation function of (c4). (d5) Normalized autocorrelation function of (d4). (c6) Root mean square value of (c5). (d6) Root mean square value of (d5). (c7) Spectrum of fusion signals. (d7) Spectrum of fusion signals.
From Figure 9(c7), examining the spectrum of the compound fault of the rolling element and outer ring, the following features are observed: (1) A prominent frequency at 1250 Hz, corresponding to 28 times the characteristic frequency of the rolling element fault (1250/28 = 44.64). (2) A pronounced frequency at 1311 Hz, corresponding to 20 times the characteristic frequency of the outer ring fault (1311/20 = 65.55). (3) A notable frequency at 2548 Hz, corresponding to 39 times the characteristic frequency of the outer ring fault (2548/39 = 65.33). (4) Significant noise reduction, revealing that the extracted characteristic frequencies perfectly match the composite fault types of the bearing.
From Figure 9(d7), considering the spectrum when the bearing has a combined rolling element and inner ring fault, we note the following features: (1) A distinct frequency at 1234 Hz, which corresponds to 29 times the characteristic frequency of the rolling element fault (1234/29 = 42.55). (2) A salient frequency at 1310 Hz, equating to 12 times the characteristic frequency of the inner ring fault (1310/12 = 109.2). (3) Effective noise suppression, with extracted characteristic frequencies precisely aligning with the composite fault types of the bearing.
It is observed that, regardless of the bearing composite fault types, the proposed method effectively suppresses noise and extracts characteristic frequencies that correspond precisely to the fault types of bearings, ensuring accurate judgment of composite fault types of bearings.
Conclusion
To precisely judge the type of compound faults of bearings, the paper took the minimum envelope entropy as the fitness function of sparrow optimization algorithm and embodied the self-adaptive decomposition of VMD. According to root-mean-square value, sensitive fault component signals were effectively chosen; vector bispectrum was used to control the noise and further highlight fault information. The following conclusions can be drawn. (1) Root-mean-square value can be used as feature index of fault to help with the correct choice of sensitive component signals of fault. (2) Vector bispectrum embodies the effective blending of sensitive component signals of fault. Compared with the method of directly analyzing component signals and the one of weighting reconstruction according to weight coefficient, the proposed method can effectively suppress the noise and further highlight the fault characteristic frequency of bearings. (3) The proposed method is insensitive to the installation positions of sensors and the types of compound faults. No matter where the sensors are installed, horizontally or vertically, in different types of compound faults, the proposed method can still effectively extract fault information and precisely identify the type of compound faults of bearings. The proposed method was only verified based on experimental signals. The fault signals will be weaker and more complex in actual engineering. Therefore, the proposed method will require more pre-treatments in actual engineering application.
Footnotes
Acknowledgments
This work was supported by National Natural Science Foundation of China [grant number: 51605309], Natural Science Foundation of Liaoning Province [grant number: 2022-MS-299], Aeronautical Science Foundation of China [grant number: 20230033054001] and Department of Education of Liaoning Province [grant number: LJKMZ20220529].
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: Supported by National Natural Science Foundation of China [grant number: 51605309], Natural Science Foundation of Liaoning Province [grant number: 2022-MS-299], Aeronautical Science Foundation of China [grant number: 20230033054001] and Department of Education of Liaoning Province [grant number: LJKMZ20220529].
