Time-frequency Feature Extraction Method of the Multi-Source Shock Signal Based on Improved VMD and Bilateral Adaptive Laplace Wavelet

Vibration signals have the characteristics of multi-source strong shock coupling and strong noise interference owing to the complex structure of reciprocating machinery. Therefore, it is difficult to extract, analyze, and diagnose mechanical fault features. To accurately extract sensitive features from the strong noise interference and unsteady monitoring signals of reciprocating machinery, a study on the time-frequency feature extraction method of multi-source shock signals is conducted. Combining the characteristics of reciprocating mechanical vibration signals, a targeted optimization method considering the variational modal decomposition (VMD) mode number and second penalty factor is proposed, which completed the adaptive decomposition of coupled signals. Aiming at the bilateral asymmetric attenuation characteristics of reciprocating mechanical shock signals, a new bilateral adaptive Laplace wavelet (BALW) is established. A search strategy for wavelet local parameters of multi-shock signals is proposed using the harmony search (HS) method. A multi-source shock simulation signal is established, and actual data on the valve fault are obtained through diesel engine fault experiments. The fault recognition rate of the intake and exhaust valve clearance is above 90% and the extraction accuracy of the shock start position is improved by 10°.


Introduction
High-end mechanical equipment, such as high-power gas turbines, diesel engines, high-pressure reciprocating compressors, and large centrifugal pumps, are the core power equipment of nuclear power, petrochemicals, shipbuilding, and other industries. The development of fault monitoring, diagnosis, and health management technology for key equipment has always attracted the attention of researchers [1,2]. The feature extraction of equipment operating parameters was the basis for establishing effective fault diagnosis and state evaluation models [3,4]. Extracting sensitive features from non-stationary signals with strong noise interference is a difficult problem in early fault diagnosis [5].
Researchers have done a lot of research on feature extraction of vibration signals of rotating machinery. Rotating machinery, e.g., rotors, bearings, and gears, has a close correlation with characteristic frequencies, such as rotation frequency, passing frequency, and meshing frequency, in their dynamic characteristics under failure. Therefore, the focus of related research is on effectively extracting the frequency domain characteristics of early weak faults from the signals of non-stationary conditions [6][7][8]. In reciprocating machinery, its working process has reciprocating periodicity. The For example, diesel engine valve wear fault results in changes to the valve opening and closing shock angle and peak value [9]. Similarly, a reciprocating compressor connecting rod small end bearing wear fault results in a crosshead pin force reversal angle shock [10]. Therefore, it is important to identify the shock source from the complex interference signals and extract the time-frequency domain characteristics of the shocks. Representative methods for decomposition processing and feature extraction of various signals include wavelet transform (WT), empirical mode decomposition (EMD), ensemble empirical mode decomposition (EEMD), variational modal decomposition (VMD), etc. Recently, new methods based on neural networks have been proposed, including autoencoders (AEs) and convolutional neural networks (CNNs).
Research has been conducted on improving these methods, regarding the above-mentioned feature extraction problems. Based on the existing defects of signal decomposition methods in nonlinear signal analysis and the need to set parameters for WT, EEMD, etc., Pan et al. [11] proposed a symplectic geometry mode decomposition (SGMD) method. Li et al. [12] proposed a bandwidth-based method to select the best envelope interpolation method to improve the EMD frequency band aliasing problem. Focusing on the problem of setting the number of decomposition modes and the second penalty factor in the VMD, Zhang et al. [13] combined a genetic algorithm and nonlinear programming to adaptively optimize VMD parameters and apply them to bearing fault diagnosis. Regarding the problem of the fault characteristics of rolling bearings being affected by working conditions, Li et al. [14] proposed a knowledge-mapping adversarial domain adaptive method, which was applied to the bearing feature transfer of a CNN. The existing research mainly focused on the optimization of decomposition methods and the selection of parameters, which are mostly applicable to rotating equipment, such as rolling bearings. There have been few studies on multi-source shock recognition and adaptive time-domain feature extraction of reciprocating machinery. And related researches have mostly focused on extracting the fault frequency features of rotating machinery. There are still research gaps in the field of time-domain shock feature extraction and fault diagnosis of reciprocating machinery faults. And early fault feature extraction and recognition were difficult problems that plagued equipment health management. This paper proposes a time-frequency domain feature extraction method for multi-source shock signals based on improved VMD and bilateral adaptive Laplace wavelet (BALW) methods. First, the VMD method was improved and the vibration signals were adaptively decomposed to obtain signals of different frequency bands. Then a new BALW was then constructed to extract the shock and its characteristics in the time and angle domains. Recently, studies on VMD have mostly focused on parameter optimization. Zhao et al. [15] used envelope nesting and a multi-objective function to optimize the two VMD parameters and applied them to the fault diagnosis of rolling bearings. Tan et al. [16] proposed a fitness function based on mode mutual information to optimize VMD parameters for fault feature extraction of rolling bearings.
The multi-source shocks of the reciprocating mechanical vibration signals had certain differences in both the time and frequency domains. In view of this characteristic, this study optimizes the decomposition mode number based on the difference between adjacent decomposition modes. Simultaneously, because the shock frequency band required a higher signal-to-noise ratio (SNR) and the noise component required a larger entropy value, this study combined the normalized SNR of the shock frequency band component and the normalized power spectral entropy value of the main noise frequency band to optimize the second penalty factor.
In terms of shock recognition and feature extraction of the rotating machinery, time-domain feature extraction methods are the current research focus. Based on the vibration characteristics of the second-order underdamped system, the researchers used a Laplace wavelet (LW) to filter and extract the characteristics of the shock [17,18]. Regarding the unilateral limitation of the LW, to obtain a wavelet filter with linear phase-frequency response characteristics, researchers proposed an antisymmetric real Laplace wavelet (ARLW). Feng et al. [19] proposed a differential evolution (DE) optimization method and an ARLW filter-based method to extract the impulsive features buried in noisy vibration signals. Wang et al. [20] proposed a Bayesian inference method based on the smoothness index, guided to determine the joint posterior probability distributions of ARLW parameters, which was used to identify different bearing faults. However, the multi-source vibration response of reciprocating machinery did not completely conform to the unilateral attenuation characteristics [21,22]. Moreover, owing to the sensor characteristics, the collected vibration signals had a certain degree of distortion [23]. Therefore, ARLW is not fully applicable to the asymmetric vibration and shock signals of the reciprocating machinery. Further, two different damping ratio parameters are used in this study to establish a new BALW, which better matches the shock vibration shape and frequency. In the iterative search of the wavelet parameter interval, the harmony search algorithm [24,25] establishes a BALW parameter search strategy, which solves the automatic search problem of wavelet parameters of different timedomain shock components. Finally, the proposed method is tested using diesel engine vibration signals.
The contributions of the new method demonstrated in this paper can be summarized as follows: ① Combining the characteristics of reciprocating mechanical signals, a parameter optimization method based on VMD is proposed. The simulation signal and actual diesel engine vibration signal processing results show that the new method effectively extracted the shock components of each frequency band. ② The proposed BALW has a better matching performance to the actual shock signal and results in better accuracy when extracting the timedomain features of the shock signal, such as the shock start position and shock frequency characteristics. ③ The multi-shock signal wavelet parameter search strategy based on the HS algorithm significantly improves the speed and effect. ④ According to the fault data, the frequency and time-domain characteristics of the shock signal are extracted, which improves the ability to capture the local features of the early fault shocks. This paper is structured as follows. Section 2 introduces the optimized VMD algorithm and the principle of the BALW. Section 3 introduces the application effect of the method proposed in this study on the simulated signal. Section 4 describes the experimental equipment and experimental data application effects. Finally, the conclusions are presented in Section 5.

VMD Parameter Optimization Method
Diesel engine vibration includes the shock response of various components, and there is aliasing of multi-band signals. The fault-sensitive shock signals must be separated from them. Therefore, this study used the widely used VMD method to decompose diesel engine vibration signals. In the VMD method, it is necessary to set the modal number K and the second penalty factor α in advance; however, it is difficult to set these two parameters in practice [26,27]. This paper proposed a VMD parameter optimization method based on the characteristics of large shock vibration signals.

Optimization of the Second Penalty Factor
The second penalty factor, α , affects the signal reconstruction accuracy. It is a significant advantage that the decomposed shock signal contains less noise interference.
The shock components often exist in the mid-band of the signal so the initial decomposition mode number, K = 3, is decomposed into three parts: low-frequency IMF1, intermediate-frequency IMF2, and high-frequency IMF3. The SNR is calculated for the decomposed signal of the intermediate frequency; a higher SNR is desired. The SNR calculation formula is as follows: where the RMS is the root mean square value, IMF2 is the decomposed mid-band signal, and IMF2_noise is the decomposed non-shock component of the mid-band.
For the high-frequency part, the main component is noise, and the entropy value of the noise component is relatively large; therefore, a high-frequency component, IMF3, with a higher power spectrum entropy value must be obtained. The power spectrum of each VMD component is calculated as follows: where F IMF f is the discrete signal Fourier transform of the component, and N is the data length.
The power spectrum entropy value is: where q i is the proportion of the i-th power spectrum value in the entire spectrum. The optimization target of the second penalty factor α is the maximum value of the sum of the SNR of the normalized intermediate frequency component and the power spectrum entropy of the normalized high-frequency component. The optimization objective function is as follows: where nor(SNR(α)) is the normalized SNR vector and nor(H P (α)) is the normalized power spectrum entropy vector.

Optimization of the Decomposition Mode Number
During the VMD process of the vibration signal, the frequency band of the signal decomposition becomes finer as the value of K increases and the difference between adjacent frequency bands becomes smaller. To obtain the best signal decomposition effect, the difference between adjacent frequency band signals should be maximized. Therefore, this study is based on the principle of a minimum average Pearson correlation coefficient between According to the number of modes, K, the average Pearson correlation coefficient ρ is calculated as follows:

A New BALW
Different types of shocks had different shock response performances, mainly due to the different mechanisms and materials. For example, the ignition process of a diesel engine is a combustion process in which a singlepoint or multiple-point flame spreads to the entire area and the vibration signal response has the characteristics of first rising and then decaying. Figure 1 shows the vibration signal of diesel engine cylinder head and its partial enlarged view. The ignition shock vibration does not conform to the characteristics of LW unilateral attenuation. Meanwhile, the signal decomposition will cause the leakage of impact energy and change its waveform.
Therefore, the commonly used LW derived from the Laplace inverse transformation of the second-order underdamped system was not completely suitable for the reciprocating mechanical shock signal. In this study, a BALW was constructed to meet the needs of matching various impact signals.

Principle of the BALW
The LW is a complex exponential wavelet with unilateral oscillation attenuation, which conforms to the characteristics of unilateral attenuation of actual shock signals. The complex LW is expressed as: where the parameter vector γ = ξ , f determines the characteristics of the wavelet. f ∈ R + is the frequency, which determines the LW oscillation frequency. ξ ∈ [0, 1) is the damping ratio, which determines the decay rate of the LW.
The anti-symmetric real LW is treated as antisymmetric on the basis of the LW: The BALW is proposed on the basis of the ARLW, which has the characteristics of bilateral asymmetric attenuation. When the unilateral damping value is large, the waveform characteristics are close to the unilateral damping type impact. When the bilateral damping values are both within a reasonable range, the waveform characteristics conform to the impact characteristics of first increasing and then attenuating. Therefore, the BALW has a better impact matching ability. The BALW is expressed as: where ξ 1 , ξ 2 ∈ [0, 1) is the damping ratio.
In the actual signal processing, the real part of the BALW is used, Re ψ γ :

Wavelet Parameter Search Method
The wavelet parameter set is φ = F × Z 1 × Z 2 × T , a four-dimensional matrix, and hence the step-by-step iterative search process for each parameter has a significant computational cost. Thus, a search algorithm is needed to quickly search for parameters. The harmony search (HS) algorithm [24] is an optimization algorithm inspired by the beautiful harmonizing of musical instruments. It has been widely studied and applied in the field of combinatorial optimization. The main features of the algorithm are fewer setting parameters, strong solution ability of continuous optimization and discrete optimization methods, and a simple concept. This study used the HS algorithm to accelerate the search process of wavelet parameters.
The search objective function is defined as f (x) , where X = [F , Z 1 , Z 2 , T ] . Set the initialized harmony memory bank HM, the HM considering rate (HMCR), the pitch adjusting rate (PAR), the fine-tuning amplitude b w , and the maximum number of iterations MAX. HM and HMCR are used to generate a new solution X new i as follows: where X rand i,j is a random individual generated in the possible range of values.
The PAR parameter is used to adjust the pitch, as follows: The maximum iteration number was selected as the stopping criterion of the analyses in this study, although there are various alternative criterions.

Simulation Signal Generation
In actual vibration signals, strong shocks can excite twoorder or even multi-order vibration modes, and there is coupling and superposition of different shock signals. Hence, single-frequency wavelet matching is not suitable. To reflect the applicability of the multi-impact signal search strategy proposed in this study, a simulation signal s is constructed for verification, where s 1 is a bilaterally attenuated impulse signal, s 2 is the double frequency signal of s 1 at the same time as τ 1 , s 3 is another unilateral attenuation shock signal at time τ 2 , and s 4 is the lowfrequency sinusoidal signal. The waveform diagram and frequency domain diagram of the simulated signal s are shown in Figures 2 and 3, respectively.

VMD Processing of the Simulated Signal
The optimized VMD is then performed on the simulated signal. According to the optimized VMD parameter optimization method proposed in this paper, the searched optimal secondary penalty factor α = 859 and the optimal decomposition modal number K = 4, as shown in Figures 4 and 5. The optimized VMD modal signals and their spectrograms are shown in Figure 6. It can be seen that the optimized VMD effectively decomposed the shock signal components, among which IMF2 is a 40 Hz signal component and IMF3 is an 80 Hz signal component.
Further, the unoptimized VMD is used to process the simulation signal. The parameters are set as α = 2000 , K = 5 [28]. As shown in Figure 7, the simulation signal is  The cosine distance shown in Eq. (18) is used to measure the similarity between the simulated shock signal and the decomposed signal.
As shown in Table 1, The decomposed signals based on optimized VMD processing has smaller cosine distance, which indicates that it is closer to the original signal and has less energy loss.  Table 1 The cosine distance of simulated shock signal and the decomposed signal To further verify the decomposition effect and advantages of the improved VMD method proposed in this paper, it is also compared with other popular signal decomposition methods in recent years. EMD is often used in fault diagnosis, especially in bearing vibration signal decomposition and feature extraction [29,30]. However, there are end effects in the EMD, which affect the decomposition effect. EEMD introduces white noise to improve the end effect [31], but it has the same modal aliasing problem as EMD.
The simulation signal is processed by the EEMD method. As shown in Figure 8, there are 40 Hz and 80 Hz frequency bands in IMF1. IMF3 and IMF5 contain part of the 10 Hz harmonic signal of IMF4. Therefore, EEMD has a strong frequency band aliasing phenomenon in the process of processing the above signals. It shows the advantages of optimized VMD in processing multi-shock signals.

Search Strategy for Wavelet Parameters of the Multi-shock Signals
A group of vibration signals often contains multi-shock components. During the global search process for the optimal wavelet parameters, there are multiple shocks corresponding to the optimal local parameter values. However, the harmony search algorithm searched for the global optimal solution and only one of the shock components could be found. Therefore, the search algorithm lacked the ability to search for wavelet parameters for multi-shock vibration signals. Aiming at the above problems, a wavelet parameter harmony search strategy for multi-shock signals was proposed.
(1) Perform VMD processing on the original signal to obtain the frequency band signal of the shock component, and determine the search frequency interval of each frequency band signal. (2) Calculate the peak position of each shock frequency band signal separately to determine the time point parameter τ . The time interval T ∈ [τ − δ, τ + δ] can be set to capture the peak position more accurately owing to noise interference. (3) Calculate the best wavelet parameters based on the cosine distance. Set the wavelet envelope Env cutoff point threshold to 0.01 (the threshold size can be selected and set as required) and zero the signal    between the wavelet cut-off points of the original signal. The signal envelope formula is as follows: (4) Repeat (2) and (3) until the cosine distance is less than the threshold of 0.2 so as to complete the wavelet matching parameter search of all shock components.

Multi-shock Recognition of the Simulated Signal
The wavelet parameter search strategy of multi-shock signals is used to match the wavelet parameters of the two frequency band signals of IMF1 and IMF2. First, we determine the frequency search interval The damping ratio at the rising end of the search waveform of the s 2 signal differed from the actual value by 0.01, while the remaining parameters matched. This is caused by the loss of the signal decomposition. The damping ratio of the falling end of the search of the s 3 signal matches and the damping ratio of the rising end reaches 0.21, which is consistent with the original decomposed signal in terms of shape and frequency. Meanwhile, the calculation time between different search methods is compared. As shown in Table 2, the HS algorithm search time is 3.15 s, the iterative search time is 361.23 s, and the differential evolution algorithm search time is 3.53 s. The search time of HS algorithm is better than that of other methods.
As shown in Figures 9, 10, 11, 12, in terms of vibration waveforms, the decomposed modal signals and the parameter matching signals have a high degree of similarity, indicating that the method proposed in this paper has a good shock parameter extraction result. The specific matching parameters are listed in Table 3.
We analyzed and compared the search results and search time of the iterative search, differential evolution algorithm, and harmony search algorithm for the s 1 signal. All three methods show consistent search results.

Test Bench Introduction and Fault Setting
We conduct a confirmatory experiment on a TBD234V12 direct-injection diesel engine. The common abnormal   valve clearance fault of diesel engines is selected as the analysis object. According to the transmission path of vibration, the vibration generated by the valve finally reflects the surface of the diesel engine cylinder head. Hence, the upper surface of the diesel engine cylinder head is taken as the vibration measuring point (as shown in Figure 13). This diesel engine is a 12-cylinder diesel engine and its cylinders are divided into columns A and B. The vibration signals and pulse signals are sampled by the DATA acquisition (DAQ) system. The DAQ card has a 24-bit ADC resolution, with a maximum sampling rate of 102.4 K/s for each channel and a maximum of 32 analog inputs. Under normal circumstances, the intake valve clearance is 0.3 mm and the exhaust valve clearance is 0.5 mm. Taking the intake valve as the research object, the B4 cylinder is used in an intake clearance fault simulation experiment. The intake valve clearance is set to 0.3 mm and 0.9 mm and the exhaust valve clearance is set to 0.5 mm and 1.1 mm. The valve clearance setting process is shown in Figure 14. When the diesel engine is running at 1500 r/min, 1800 r/min, and 2100 r/min, the cylinder head vibration signals are collected.

Feature Extraction of Valve Clearance Fault
According to the flywheel gear pulse signals, the vibration signals are intercepted throughout the period, and the time-domain vibration signals are converted to correspond to the crankshaft angle. Figures 15 and 16 show the vibration signals of B4 cylinder under normal and fault conditions. Compared with the normal signal, the intake valve closing shock peak in the fault state is significantly larger (shown in the red box). The shock amplitude of the intake valve closing increases and the starting position is advanced (the dotted frame area near the 480° crankshaft angle). However, the valve closing shock amplitude is affected by the operating conditions of the diesel engine and the valve closing shock start phase is relatively fixed. Thus, the shock start position is a very critical fault feature. From the spectrum analysis (as shown in Figures 17  and 18), the shock frequency band is mainly around According to the VMD parameter optimization method proposed in this paper, VMD ( K = 4, α = 61 , the iterative process of parameter optimization is shown in Figures 19 and 20) processed modal decomposition signals and their spectra are shown in Figure 21.
It is determined that the shock frequency band is mainly concentrated at around 5 kHz, i.e., the IMF2 signal. The search frequency band F ∈ 4 k : 0.1 k : 6 k and the damping ratio range Z ∈ {0.01 : 0.01 : 0.99} are set, and a wavelet parameter search for the valve closing shock is performed according to the method proposed in this paper. As shown in Figures 22, 23, 24, the BALW proposed in this paper has a better waveform matching ability for the closing shock of the intake valve after VMD than the LW and ARLW. The proposed method based on the BALW accurately extracts the starting angle of the valve closing shock. However, the extraction angle of the other two methods differed by approximately 10° from the real situation, which induces significant interference with the actual valve fault monitoring process; this causes fault misjudgment and is thus not suitable for actual situations.
Moreover, the initial angle position and peak characteristics of 360 groups of inlet valve closing shocks are Crest factor  Figures 25 and  26, the feature of initial shock angle has more advantages than the feature of peak value in fault recognition.

Fault Identification Result
The abnormal valve clearance fault signals of diesel engines are used for fault diagnosis. As shown in Figure 27, the main process of this fault diagnosis process is as follows: (1) The vibration signals of the diesel engine cylinder head are processed by optimized VMD. Then, the signals of the shock frequency band are extracted. (2) The shock search algorithm based on BALW and the harmony search algorithm is used for shock search. Then the feature values are extracted to train the Softmax classifier, and finally, the fault classification is realized. The extracted features are shown in Table 4.  Table 4. (4) Finally, compare the two methods in (2) and (3).
The commonly used Softmax classifier is used for fault classification. There are 240 sets of data for the three speeds, a total of 720 sets, of which normal and fault data each account for half. The 480 sets of data are randomly selected as the training set and the remaining 360 sets of data as the test set. The clearance fault of the intake and exhaust valve is then recognized. Table 5 shows the classification results. The fault recognition rate of the intake and exhaust valve clearance of the methods (optimized VMD + BALW) proposed in this paper was above 90%. However, the fault recognition rate of the intake and exhaust valve clearance for the VMD + conventional features method was 60%-70%. The results showed that the shock feature extraction method proposed in this paper realized the effective extraction of the feature and obtained an accurate fault recognition effect.

Conclusions
Based on the characteristics of multi-source shock coupling and strong noise interference in the shell vibration signal of reciprocating machinery, this study improved the parameter selection of the VMD method.
By introducing different bilateral damping ratio parameters to improve the ARLW, a new BALW was established, which had bilateral asymmetric attenuation characteristics. Considering the difficult problem of optimal selection of wavelet local parameters for multi-shock signals, a wavelet local parameter search strategy for multi-shock signals was established via the harmony search algorithm. The search time was 3.15 s, which improved upon the other methods.
The new method proposed in this paper solved the problem of accurately extracting different types of timedomain shock signals, and realized the optimization of wavelet parameters of multi-source shocks and the adaptive extraction of local shock features. A diesel engine test bench was used to obtain abnormal data for simulated valve clearance, the shock feature extraction effect was verified, and the fault recognition accuracy rate exceeded 90%. Moreover, the extraction accuracy of the shock start position was improved by 10°.
In the future, we will combine other methods to further improve the performance of shock extraction.