Abstract
Variational mode decomposition (VMD) has been applied in the field of rolling bearing fault diagnosis because of its good ability of frequency segmentation. Mode number and quadratic penalty term have a significant influence on the decomposition result of VMD. At present, the commonly used method is to determine these two parameters adaptively through intelligent optimization algorithm, namely, the parameter-adaptive VMD (PAVMD) method. The key of the PAVMD method is the setting of an objective function, and the traditional PAVMD method is prone to overdecomposition or underdecomposition. To solve these problems, an improved parameter-adaptive VMD (IPAVMD) method is proposed. A new objective function, the maximum average envelope kurtosis (MAEK), is proposed in this paper. The new objective function fully considers the equivalent filtering characteristics of VMD, and squared envelope kurtosis has good antinoise performance. In the optimization method, this paper uses an improved particle swarm optimization (PSO) algorithm. The MAEK and PSO can make sure the IPAVMD method reaches the best complete decomposition of the signal without an underdecomposition or overdecomposition problem. Through the analysis of simulation data and experimental data, the performance of the IPAVMD and the traditional PAVMD is compared. The comparison results show that the proposed IPAVMD has better performance and stronger robustness than the traditional method and is suitable for both single-fault and multiple-fault cases of rolling bearings. The research results have certain theoretical significance and application value for improving the fault diagnosis effect of rolling bearings.
1. Introduction
Rolling bearing is widely used in rotating machinery, which is one of the easily damaged parts [1]. Its running condition affects the stability of the whole rotating machinery system. As a key component of the transmission system, it is of great practical significance for the safe operation of the whole rotating machinery whether the fault and fault type can be detected in time.
In bearing fault diagnosis, the collected vibration signals are generally nonlinear and nonstationary, and the fault feature information contained in them is often submerged by strong background noise [2]. How to efficiently and accurately extract useful fault features from complex vibration signals has become the focus and difficulty of fault diagnosis.
Empirical mode decomposition (EMD) [3] is an adaptive signal processing method, which can decompose nonstationary signals into modal components of different frequency bands and is widely used in signal processing [4]. However, EMD has some problems such as end effect and mode aliasing [5]. Ensemble empirical mode decomposition (EEMD) is an improved algorithm of EMD, which can effectively suppress the mode aliasing phenomenon [6]. However, due to the increase of the number of iterations, the operation efficiency is reduced, and the endpoint effect problem is not solved. Local mean decomposition (LMD) [7] improves the problem of overenvelope or underenvelope in EMD, but its essence is the same as that of EMD and EEMD [8], which belongs to recursive mode decomposition, and ultimately cannot avoid the end effect and mode aliasing phenomenon. At the same time, these recursive decomposition methods lack strict mathematical theory support. For the abovementioned problems, Dragomiretskiy and Zosso [9] proposed a nonrecursive mode decomposition method according to the binding variational problem, namely, Variational Mode Decomposition (VMD). The algorithm can segment the signal flexibly in the frequency domain and can effectively mine the submerged feature information in the signal. Once put forward, it was favored by many scholars [10–14].
The researchers found that the signal processing effect of VMD was influenced by its parameter setting. Li et al. [15] determined the optimal mode number adaptively by peak search and similarity principle. Guided by a series of indexes such as permutation entropy and kurtosis, Lian et al. [16] used iterative search to determine the mode number . In the application of VMD, the abovementioned two references only considered the influence of mode number and ignored the influence of penalty factor on modal bandwidth [17–19]. In [20], and were determined by the central frequency observation method, which depended on experience and had poor adaptability. At present, one of the popular methods is to adaptively determine the optimal parameter combination of VMD [, ] by using an intelligent optimization algorithm. In [21], simulated annealing (SA) algorithm was used to optimize VMD with kurtosis as fitness function for early chatter identification in the milling process. In [22], particle swarm optimization (PSO) algorithm was used to optimize VMD with maximum entropy as the fitness function to realize liquid pipelines leakage detection. In the field of rolling bearing fault diagnosis, Zhang et al. [23] introduced grasshopper optimization algorithm (GOA) to optimize VMD with the weighted kurtosis as the fitness function. Nassef et al. [24] adaptively selected the best parameters of VMD by using sailfish optimization (SFO) algorithm and Gini index as the fitness function. Gai et al. [25] proposed an improved parameter-adaptive VMD method based on the hybrid gray wolf optimizer (HGWO) and envelope entropy. The abovementioned three studies have achieved good results, but this kind of method is only for single-fault feature extraction, making a component to achieve global optimization, and cannot achieve bearing compound fault analysis. One step closer, Miao et al. [26] proposed a new index, namely, ensemble kurtosis (EK). The grasshopper optimization algorithm (GOA) was used to optimize the VMD with the minimum mean ensemble kurtosis as the objective function, so that the optimized VMD was suitable for single-fault and compound-fault diagnosis.
To sum up, we can find that VMD based on optimization algorithm can adaptively determine and , but the optimization result depends heavily on the objective function. Therefore, constructing an appropriate objective function is the key to the parameter-adaptive VMD method. The authors find that the ensemble kurtosis is still affected by cyclic stationary noise, and the parameter-adaptive VMD proposed by Miao et al. [26] has the problem of overdecomposition or underdecomposition. Further research finds that the ensemble kurtosis is greatly affected by random noise, while the squared envelope kurtosis has strong antinoise performance. Therefore, a new objective function, the maximum average envelope kurtosis (MAEK), is proposed. An improved parameter-adaptive variational mode decomposition method is established by using particle swarm optimization algorithm. The proposed method has no overdecomposition or underdecomposition problem, and its robustness and generalization are outstanding. It is suitable for both single-fault diagnosis and compound-fault diagnosis of rolling bearings. The correctness and efficiency of the new method are verified by simulation and experiment.
The rest of this paper is summarized as follows. In Section 2, we briefly introduce the basic principles of VMD and PSO, clarify the original intention of establishing the new objective function, and introduce the detailed steps of rolling bearing fault diagnosis proposed in this paper. In Section 3, the correctness of IPAVMD is verified through the analysis of single-fault and multiple‐fault simulation signals of bearings. The improved method is compared with the traditional PAVMD proposed in [26], which proves the superiority of the improved method. In Section 4, we use two actual bearing signals to check the performance of the new method and verify the effectiveness of the proposed method. Finally, the conclusions are drawn in Section 5.
2. Improved Parameter-Adaptive Variational Mode Decomposition (IPAVMD)
2.1. VMD Method
VMD is a completely nonrecursive signal decomposition mode, which realizes the adaptive decomposition of the signal by iteratively searching the optimal solution in the variational models. An actual nonstationary signal can be decomposed into discrete modal components , and the center frequency and bandwidth of each component can be determined. VMD algorithm consists of two steps: constructing the variational problems and solving the variational problems. The expression of the constrained variational model is shown as follows:where denotes the component signal, represents the center frequency, and is the mode number.
Quadratic penalty factor and Lagrange multiplier are introduced to transform equation (1) into an unconstrained optimization model, and the expression is as follows:
For equation (2), the alternate direction method of multipliers (ADMM) is used to find the optimal solution. Each mode and the corresponding center frequency are updated by the following equations, respectively:
When and are updated, Lagrange multiplier is also updated by
The VMD iteration process stops until the relative error is less than the convergence tolerance .
2.2. Particle Swarm Optimization (PSO)
Particle swarm optimization (PSO) has been a popular algorithm for several decades because of its characteristics of evolutionary computation and swarm intelligence and good global optimization ability. Many scholars [22, 27, 28] have proposed different improved particle swarm optimization algorithms to improve convergence speed and accuracy. Zhang et al. [28] successfully applied an improved particle swarm optimization algorithm to VMD parameter optimization. This paper will also use it for parameter optimization.
The specific optimization steps are as follows:(1)Random generation of particle swarm locations and velocities.(2)To evaluate the fitness of each particle, the position and fitness of the particle are stored in the individual extremum , and the global extremum is generated at the same time.(3)The particles velocity and position are updated according to where is the inertia weight, , , , are, respectively, the current iteration number and the maximum iteration, and are the learning factors, and and are the random number between . The specific initial parameter values are shown in Table 1.(4)All and are compared, and is updated.(5)determined if the termination condition is met; the search results are output if it is met; otherwise, the iteration is continued.
2.3. The New Objective Function: Maximum Average Envelope Kurtosis (MAEK)
Kurtosis is very sensitive to fault shock signals and is often used as an evaluation index for the diagnosis of surface damage of rolling bearings [29, 30]. Yang et al. [31] pointed out that the bearing fault signal has impulsive characteristics in time domain and cyclic stationarity characteristics in the frequency domain, which should be considered at the same time. If only one of them is considered, it will be interfered by impulsive noise or cyclic stationarity noise, which will affect the diagnosis effect. In [26], the time-domain kurtosis was used to evaluate the impact property of the signal, the frequency-domain kurtosis was used to evaluate the cyclic stability of the signal, and the ensemble kurtosis (EK) index was proposed. It is defined as follows:where denotes the ensemble kurtosis, is the kurtosis of the signal, is the envelope spectrum kurtosis, denotes the envelope spectrum of the signal, and is the sampling number of the envelope spectrum.
In [26], grasshopper optimization algorithm (GOA) was used to optimize the parameters of VMD. The objective function as shown in equation (9) is established. The optimized VMD is used to decompose the signal, and the component with the maximum ensemble kurtosis value is selected, which is considered to contain the bearing fault feature information.where denotes the ensemble kurtosis of each mode and is the mode number.
It is undeniable that EK index is more suitable for evaluating bearing fault impulses. However, the authors find that the parameter-adaptive VMD method based on EK has the problem of overdecomposition or underdecomposition when processing the signal. The reason is that EK index is seriously affected by random noise. Bearing fault signals generally include the cyclic stationary interference part caused by shaft rotation and gear engagement, the fault impulses part, and the random impulses part caused by electromagnetic interference. We analyzed the change of the EK value of these three parts under different signal-to-noise ratios (SNRs), as shown in Figures 1(a)–1(c), and the simulation signals used are described in detail in Section 3.1. It can be intuitively found that the EK values of these three parts vary greatly under different SNRs. Meanwhile, the objective function established by equation (9) is to average the EK values of all components obtained by VMD decomposition. The essence of VMD is multiple adaptive Wiener filters [32]. In the process of searching and , the SNR of the signals in each filter varies greatly, which makes the EK value of the abovementioned three parts fluctuate greatly, resulting in that the determined mode number may not be optimal.
(a)
(b)
(c)
(d)
(e)
(f)
Further study found that the squared envelope kurtosis (SEK) of the signal is less affected by the random noise. Similarly, we analyze the changes of SEK values of the abovementioned three parts of signals under different SNRs, as shown in Figures 1(d)–1(f). It can be clearly found that the SEK value fluctuates little under different SNRs. Therefore, SEK has stronger noise resistance and better stability than EK. In this paper, a new indicator, average envelope kurtosis (AEK), is proposed in this paper. AEK refers to the average value of the squared envelope kurtosis of each modal component obtained by VMD under specific parameters and . The maximum average envelope kurtosis (MAEK) is taken as the objective function, and the optimal parameter combination of VMD is obtained through PSO search, as shown in equation (10). The improved parameter-adaptive VMD method can decompose the cyclic stationary interference part, fault pulse part, and random pulse part at one time, and the problem of overdecomposition or underdecomposition does not occur easily.where and are optimal parameters and is the squared envelope kurtosis value of each component, which is calculated as follows [33]:where is the squared envelope of each component, is the mean value of , is the standard deviation of , and represents the expected value.
2.4. Proposed Method
The detailed steps of the rolling bearing fault diagnosis method proposed in this paper are as follows: Step 1: the range of penalty factor is set as [100, 5000], and the range of mode number is set as . Step 2: parameters in PSO are set according to Table 1. The MAEK is taken as the optimization target to search for the best parameter combination . Step 3: the optimized VMD is applied to decompose the vibration signal and generate the time-domain waveform of each component and the corresponding squared envelope spectrum. Step 4: the envelope spectrum of each component is observed and compared with the theoretical fault frequency, and the fault type is obtained.
The flowchart of the proposed method is illustrated in Figure 2.
In order to demonstrate the correctness and generalization ability of the IPAVMD method, a series of bearing signals, including single fault and compound fault, are used for simulation and experimental verification. The performance of the IPAVMD method and PAVMD method is compared. The comparison results show that IPAVMD can realize the complete decomposition of the signal and does not easily show the underdecomposition or overdecomposition phenomenon.
3. Simulation
3.1. Simulated Signals of Case 1: The Single-Fault Signal
In order to verify the correctness of the method proposed in this paper, the single-fault simulation signal of rolling bearing used in [31, 34, 35] is used for analysis, and the simulation model is as follows:
The five parts in equation (12) represent the harmonic vibration of the rotor, the cyclic stationary interference of the gear meshing, the repetitive transients caused by the bearing fault, high-amplitude random impulse, and Gaussian distributed white noise, respectively, where represents the rotating frequency, means the number of teeth of the gear, represents the fault period, and denotes the roller random slips of [−0.02 , 0.02 ]. The impulse response function is defined as follows:
All the parameters used in the simulated signals of case 1 are listed in Table 2. The sampling frequency and sampling duration are 10000 Hz and 1 s, respectively. The waveforms of the repetitive transients, high-amplitude random impulse, and pure signal are demonstrated in Figures 3(a)–3(c). Gaussian white noise with an SNR of −0 dB is added to the pure signal to generate a fault simulation signal. Its time-domain waveform, frequency spectrum, and squared envelope spectrum are shown in Figures 4(a)–4(c), respectively. No fault characteristic information can be found in the frequency spectrum and envelope spectrum; only the frequency of rotation and gear meshing can be found.
(a)
(b)
(c)
(a)
(b)
(c)
Using the method proposed in this paper to analyze this data, the PSO is used to optimize the VMD. Figure 5 shows the curve of the maximum average envelope kurtosis with the population evolution algebra, which converges in the fourth generation, and the best parameter combination is . According to the optimization results, the parameters in VMD are set, and then, the fault simulation signal is decomposed. The time-domain waveform and the squared envelope spectrum of each component are shown in Figure 6. It can be clearly found from the envelope spectrum that u2 is the bearing fault information. We can see from Figure 6 that the cyclic stationary interference signal and random impulse signal are also well extracted. Figure 7 shows the frequency spectrum of each component. It can be found that the center frequency interval of each mode is far away and there is no overdecomposition phenomenon. Through the abovementioned analysis, the correctness of the proposed method is proved, and this method can realize the complete decomposition of bearing fault signal.
(a)
(b)
(c)
(d)
(e)
(f)
To illustrate the advantages of this method, PAVMD [26] is used to analyze the simulated signal. Figure 8 shows the optimization process, and according to the results of PSO, we set the optimal parameter combination of VMD as . The final decomposition result is shown in Figure 9, and we can easily find that the component u2 is the bearing fault information. Miao et al. [26] pointed out that the component with the largest EK value after decomposition is the component containing fault information. However, it can be seen from Figure 8(b) that the maximum EK value is the component u1. It is proved that EK index is easily disturbed by cyclic stationary noise and cannot point to fault component accurately.
(a)
(b)
(a)
(b)
(c)
(d)
3.2. Simulated Signals of Case 2: The Compound Fault Signal
In order to illustrate the generality of IPAVMD, it is applied to the signal analysis of rolling bearing compound fault. In this paper, the simulation signal of case 1 is recorded as the outer-race fault simulation signal, and the inner-race fault impulses signal is added to form the compound fault signal. The formula of inner-race fault impulses is as follows [26]:where , inner-race fault period , natural frequency = 3600, denotes the roller random slips of [−0.01 , 0.01 ], and initial phase .
The waveforms of the outer-race fault impulses, inner-race fault impulses, high-amplitude random impulse, and pure signal are demonstrated in Figures 10(a)–10(d). Gaussian white noise with an SNR of −0 dB is added to the pure signal to generate a compound fault simulation signal. Its time-domain waveform, frequency spectrum, and squared envelope spectrum are shown in Figures 11(a)–11(c), respectively. No fault characteristic information can be found in the frequency spectrum and envelope spectrum; only the frequency of rotation and gear meshing can be found.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
The new method proposed in this paper is used to analyze the compound fault signal. The PSO is used to optimize the VMD. Figure 12 shows the curve of the maximum average envelope kurtosis with the population evolution algebra, which converges in the fifth generation, and the best parameter combination is . According to the optimization results, the parameters in VMD are set, and then, the fault simulation signal is decomposed. The time-domain waveform and the squared envelope spectrum of each component are shown in Figure 13. It can be clearly found from the envelope spectrum that u2 is the outer-race fault information and u3 is the inner-race fault information. At the same time, it can be seen from Figure 13 that the cyclic stationary interference signal and random impulse signal are also well extracted. Figure 14 shows the frequency spectrum of each component. It can be found that the center frequency interval of each mode is far away and there is no overdecomposition phenomenon. The method proposed in this paper can decompose all kinds of components in bearing compound fault at one time. Through the abovementioned analysis, the correctness and generality of the proposed method are proved.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Similarly, the PAVMD [26] method is used to analyze the simulation signal. Figure 15 shows the optimization process, and according to the results of PSO, we set the optimal parameter combination of VMD as . The final decomposition result is shown in Figure 16. We can clearly find that u1 is the cyclic stationary interference signal, u2 is the inner-race fault information, and there is no outer-race fault information. PAVMD has the problem of underdecomposition, and it cannot separate each single fault in the compound fault at one time. Miao et al. [26] pointed out that the component with the largest EK value after decomposition is the component containing fault information. However, it can be seen from Figure 15(b) that the component u1 with the largest value of EK is not u2. Therefore, it is proved that EK index is easy to be disturbed by cyclic stationary noise and cannot point to fault component accurately.
(a)
(b)
(a)
(b)
(c)
(d)
4. Experimental Verification
4.1. Experimental Signals of Case 1: The Single-Fault Signal
In order to verify the effectiveness of the new method, QPZZ-II rotating machinery vibration and fault simulation experimental bench is used to verify. The test bearing model is 6205-2RS, and a 0.2 mm diameter pitting fault was machined on the outer ring raceway by EDM. The test bench is shown in Figure 17. The main parameters of the test bearing are shown in Table 3.
In the experiment, the rotational speed is1478 r/min, so the rotating frequency = 24.63 Hz, and the sampling frequency is 25.6 kHz. According to the geometric parameters of test bearing, the outer-race characteristic frequency of the bearing can be calculated, BPFO = 88.31 Hz. The acceleration sensor is applied to collect vibration data, and its time-domain waveform, frequency spectrum, and squared envelope spectrum are shown in Figure 18. Due to the noise interference, the outer-race fault characteristic frequency spectrum line is not prominent in the envelope spectrum, and there are many interference lines.
(a)
(b)
(c)
The method proposed in this paper is used to analyze this vibration signal. The PSO is used to optimize the VMD. Figure 19 shows the curve of the maximum average envelope kurtosis with the population evolution algebra, which converges in the fourth generation, and the best parameter combination is . According to the optimization results, the parameters in VMD are set, and then, the vibration signal is decomposed. The time-domain waveform and the squared envelope spectrum of each component are shown in Figure 20. It can be clearly found from the envelope spectrum that u1 is the rotating frequency information and u2 is the outer-race fault information. Figure 21 shows the frequency spectrum of each component. It can be found that the center frequency interval of each mode is far away and there is no overdecomposition phenomenon. Through the abovementioned analysis, it can be seen that the effect of the proposed method on signal separation is quite outstanding.
(a)
(b)
(c)
(d)
In order to illustrate the advantages of the proposed method, as a comparison, PAVMD [26] is used to analyze this experimental signal. Figure 22 shows the optimization process, and according to the results of PSO, we set the optimal parameter combination of VMD as . The final decomposition result is shown in Figure 23. We can clearly find that u4 is the outer-race fault information. However, the components u1–u3 all contain rotating frequency information, indicating that there is overdecomposition in this decomposition. Miao et al. [26] pointed out that the component with the largest EK value after decomposition is the component containing fault information. However, it can be seen from Figure 22(b) that the component u1 with the largest value of EK is not u4. Through the test analysis, it is easy to see that EK index is easily disturbed by cyclic stationary noise and cannot accurately point to the fault component.
(a)
(b)
4.2. Experimental Signals of Case 2: The Compound Fault Signal
The complex working environment and strong background noise of freight trains make it difficult to extract fault features of wheel set bearings, and the compound fault diagnosis is more difficult.
The test bench is shown in Figure 24(a). The RD2 wheel set is installed on the test bench. The main structural parameters of the bearing are shown in Table 4. The bearing used in the test is the maintenance bearing of a freight train after long service, and there are many kinds of damage and wear inside it. There is an obvious spalling fault on the inner ring raceway, as shown in Figure 24(b), while there are only few slight indentations on the outer ring raceway, as shown in Figure 24(c).
(a)
(b)
(c)
The rotating frequency = 7.75 Hz. According to the geometric parameters of test bearing, the fault characteristic frequencies of the bearing can be calculated, outer-race characteristic frequency BPFO = 66.75 Hz and inner-race fault characteristic frequency BPFI = 88.25 Hz. The acceleration sensor is applied to collect vibration data, and its time-domain waveform, frequency spectrum, and squared envelope spectrum are shown in Figure 25. Because of the interference of background noise, no obvious fault shock can be found in the time-domain signal, while only the inner ring fault base spectrum line is prominent in the squared envelope spectrum, and the other spectrum lines are messy. Therefore, the bearing can only be judged to have an inner ring fault.
(a)
(b)
(c)
Using the method proposed in this paper to analyze this data, the PSO is used to optimize the VMD. Figure 26 shows the curve of the maximum average envelope kurtosis with the population evolution algebra, which converges in the second generation, and the best parameter combination is . According to the optimization results, the parameters in VMD are set, and then, the vibration signal is decomposed. The time-domain waveform and the squared envelope spectrum of each component are shown in Figure 27. It can be clearly found from the envelope spectrum that u2 is the inner ring fault information and u3 is the outer ring fault information, and the two single faults are clearly separated.
(a)
(b)
(c)
(d)
(e)
(f)
To better reflect the effect of fault feature extraction, the time-domain waveform and squared envelope spectrum of u2 are separately drawn, as shown in Figures 28(a) and 28(b). The continuous impact information can be seen clearly in Figure 28(a). In Figure 28(b), the rotating frequency and its double frequency can be clearly found. The first six harmonics of the inner-race fault frequency BPFI as well as the modulation sidebands can also be clearly found. Similarly, the time-domain waveform and envelope spectrum of u3 are drawn, as shown in Figure 29. The characteristic frequency information of the outer-race fault can be clearly found in Figure 29(b). Figure 30 shows the frequency spectrum of each component. It can be found that the center frequency interval of each mode is far away and there is no overdecomposition phenomenon. The method presented in this paper can automatically separate the single-channel compound fault signals at one time, and there is no overdecomposition or underdecomposition problem. Through the abovementioned experimental analysis, the correctness of the proposed method is further verified.
(a)
(b)
(a)
(b)
In order to reflect the superiority of the improved method, the traditional PAVMD [26] is used to analyze this experimental signal. Figure 31 shows the optimization process, and according to the results of PSO, we set the optimal parameter combination of VMD as . The final decomposition result is shown in Figure 32, and we can easily find that the component u2 is the inner-race fault feature information. The component u2 has the largest EK value, as shown in Figure 31(b). It is not difficult to find that the traditional PAVMD method can only find the fault information of the inner ring, but not the weak fault information of the outer ring, so there is an underdecomposition problem.
(a)
(b)
(a)
(b)
(c)
(d)
5. Conclusions
(1)Parameter-adaptive VMD method can adaptively determine the two key parameters and , and its core is to construct the appropriate objective function. Different objective functions will lead to different optimization results. The traditional parameter-adaptive VMD method based on average ensemble kurtosis is prone to underdecomposition or overdecomposition in signal processing. The reason is that the ensemble kurtosis value fluctuates greatly under different SNRs, and the ensemble kurtosis index is greatly affected by random noise. In order to overcome this shortcoming, the squared envelope kurtosis index with strong noise resistance is introduced, and a new objective function, the maximum average envelope kurtosis, is constructed. Compared with the traditional method, the performance of the parameter-adaptive VMD method guided by the new objective function has been greatly improved. Through the analysis of two groups of bearing fault simulation signals and two groups of test signals, the correctness and superiority of this method are verified.(2)The improved parameter-adaptive VMD method proposed in this paper can decompose the cyclic stationary component, fault impact component, and high-amplitude random pulse component of the bearing fault signal at one time to achieve the complete decomposition of the signal. This method is not only suitable for single fault of rolling bearings but also for multiple-fault cases. It has strong universality and robustness and has certain engineering application value.Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare no conflicts of interest.
Acknowledgments
The present work was supported by the National Natural Science Foundation of China (nos. 11790282, 12032017, 11802184, 11902205, and 12002221), S&T Program of Hebei (20310803D), Natural Science Foundation of Hebei Province (no. A2020210028), and Shijiazhuang Tiedao University Graduate Innovation Funding Project (no. YC2021087).