# Acoustics Source Identification of Diesel Engines Based on Variational Mode Decomposition, Fast Independent Component Analysis, and Hilbert Transformation

Authors: Xieqi Chen; Jianhui Lin; Hang Jin; Andy Tan; Lei Yan
Journal: Shock and Vibration
 (2021)

---

## Abstract

Diesel engines are widely used in railway systems, particularly in freight trains. Despite their high efficiency in energy conversion, they usually generate high levels of acoustics pollution during operation. In order to mitigate this problem, a series of active/passive acoustics control methods are used to reduce noise. Most of these methods are only effective if the prior knowledge of sources is given. In other words, it is essential to recognize the acoustics source. Variational mode decomposition (VMD) is a signal processing method that enhances the signal corrupted by background noise. However, the decomposed results of VMD depend on their mode parameter and penalty parameter. Therefore, an evaluation method based on system modal parameters (natural frequency and damping ratio) is proposed to select the mode parameter, and the penalty parameter can be selected from the power spectra of signals. In order to increase the accuracy of decomposition for diesel engines and find out the sources of acoustics, a method combining VMD, fast independent component analysis, and Hilbert transformation (VMD-FastICA-HT) is proposed for the separation and identification of different sources for diesel engines. The optimization results indicate that when the penalty parameter value is 1.5 to 16 times the maximum signal amplitude, better decomposition results can be achieved. Therefore, the separated independent acoustics are more accurate in source identification. Furthermore, both simulation data and in situ operational data of diesel engines for vehicles are used to validate the effectiveness of the proposed method.


---

## Body

## 1. Introduction

Diesel engines have been widely used in railways. Diesel locomotives are designed without pantograph-catenary systems that are least affected by environmental weather, and they are powered by on-board diesel engines. Therefore, they play an irreplaceable role in railway freight and passenger transportation [1]. However, with the increase of power for internal combustion engine, the acoustics of locomotives will also increase. Most locomotives meet the requirements of acoustics standard at low speed but exceed the limitations at high speed [2]. Engines are the main sources of acoustics. Specifically, the main acoustics sources from the engine can be classified as combustion acoustics, mechanical acoustics, and exhaust acoustics. The structural acoustics of diesel engines are related to the structural modal frequency [3]. Gas acoustics can be treated as the transfer function inside and outside the cylinder, and the exhaust acoustics frequency component is determined by the ignition timing [1].In order to obtain the acoustics characteristics and the contributions of independent acoustics sources, it is necessary to separate the measured mixed acoustics [4]. However, compared with the vibration response, the acoustics signal is more pluralistic. Therefore, robust signal processing technology is needed to realize the state of detection and fault diagnosis for different systems [5].In recent years, acoustics source identification techniques are widely used in diesel engines of freight trains for acoustic reduction. The techniques, including passive [6] and active [7] acoustics control, have been applied to acoustic reduction. However, these techniques can only be effective if the characteristics of sources are known. The near sound intensity and the near-field acoustics holography technology are used to identify the acoustics and vibration of chain drive system for diesel engines [8]. The wavelet partial coherence analysis is also applied to identify the sources of acoustic signals for diesel locomotives [2]. Blacodon and Lewy [9] adopt the three-sensor coherence method to identify the acoustics source of the turboshaft engine. These methods can identify the acoustics sources according to the coherence values calculated from the signals measured by multiple sensors. Mixed transfer path analysis (TPA) can also be used to predict the internal acoustics caused by vibration of power [10], thus inhibiting the whole vehicle acoustics. Zhang et al. [11] have identified combustion acoustics of gasoline engines, mechanical meshing acoustics of air compressors in driving wheels, and electromagnetic acoustics of generators by independent component analysis (ICA). When multichannel method is used to identify the acoustics source of diesel engines, a large number of sensors are needed to measure multichannel signals simultaneously. This means that the traditional blind source separation methods cannot solve the problem of undetermined separation. In engineering testings, it is usually difficult to determine the exact number of source signals. Therefore, this method is not viable in practical engineering applications. Several researchers have denoted some signal processing methods combining ICA and correlation values. Li et al. [12] apply and verify the independent component analysis to the identification of acoustics sources for engines. Based on this method, they also propose a combined denoising algorithm, based on the ICA-CEEMD wavelet threshold, to solve the mixture of different acoustics [13]. Moreover, Wang et al. [14] combine adaptive wavelet threshold (AWT) denoising, ensemble empirical mode decomposition (EEMD), and correlation dimension (CD) to realize diesel engine fault diagnosis. However, due to the mode mixing problems of EMD method, it takes a lot of time to suppress the aliasing caused by the addition of White Gaussian Noise (WGN). Furthermore, it is also necessary to judge the pseudocomponents of EMD. Therefore, there are still several problems that need to be solved in the application of EMD-based blind source separation method [12–14] in acoustics source identification.VMD is a relatively new method that promises to replace EMD. It decomposes the multicomponent signals adaptively into multiple quasi-orthogonal intrinsic mode functions. It has been proved that VMD is superior to EMD in tone detection, separation, and acoustics robustness [15]. Therefore, the VMD-based blind source separation method is a feasible acoustics source identification method. Yao et al. [16] notice that VMD can be combined with robust independent component analysis to identify the acoustics sources of 6-cylinder diesel engines, and the results show that the method has better performance than the EMD-based methods. It can be found that the mode parameter K and the penalty parameter α mentioned in VMD have great influences on the final results, where the analysis signals are random [15, 17]. When VMD is used to identify engine acoustics sources, the result of acoustics source separation is not satisfied if the value of K is not given in advance. Recently, some progress has been made in the study of VMD parameter selection. For example, Li et al. [18] have proposed a method based on center frequency to select penalty parameters adaptively. Wang et al. [19] determine the optimal mode parameter K according to the ratio of energy in the component signal to total energy. Liu et al. [20] propose an automatic selection method based on kurtosis. Zhang et al. [21] construct the measurement index according to the correlation coefficient of kurtosis index and use the grasshopper optimization algorithm (GOA) to search for the best VMD decomposition parameters. Most of methods, including the above methods and other similar methods [19–21], can select the mode parameter K via the evaluation criteria, lacking the consideration of physical definition for systems. The selection of suitable K values under physical considerations is an essential issue. Therefore, it is a feasible method to select the appropriate K value from the steady-state points of both natural mode and damping ratio for a time-varying system.In this paper, a new parameter evaluation criterion of mode parameterK and penalty parameter K is proposed. The steady-state poles between natural frequency and damping ratio of the time-varying system are used as the evaluation standard for the selection of parameter K, and the penalty parameter α can be selected from the energy of signal spectra. In other words, the core of this evaluation method is that the system modal parameters are used as the basis for selecting parameter K during the decomposition of VMD. According to this evaluation method, both penalty and mode parameters are selected in VMD to enhance the structural vibration signals of diesel engines corrupted by background noise. In order to increase the accuracy of decomposition for diesel engines, inspired by the ICA-CEEMD wavelet threshold method, combined with VMD, fast independent component analysis algorithm, and Hilbert transformation, a new method known as VMD-FastICA-HT is proposed to separate and identify the combustion, mechanical, and exhaust acoustics in diesel engines.This paper is arranged as follows: Section2 briefly provides the background of VMD and FastICA; the process steps of the proposed VMD-FastICA-HT method are introduced. In Section 3, the method is used to analyze the numerical simulation signal and verify the selection of mode and penalty parameters based on the proposed system modal parameter evaluation. Moreover, the process of acoustics source identification is discussed. Section 4 introduces the real test by using a CAT-C18 diesel engine. Additionally, the collected in situ data is processed by the proposed VMD-FastICA-HT method to separate and identify the engine acoustics sources.


## 2. Theoretical Background

### 2.1. Variational Mode Decomposition

The VMD method abandons the process of iterative screening of the extrema of EMD method, and it introduces signal decomposition into variational model to solve the decomposition problem. In order to find the optimal solution, signal decomposition is realized by seeking for a constrained variational model. In the process of VMD decomposition, the center frequency and bandwidth of each BIMF component are alternately iterated and adaptively decomposed to the appropriate signal frequency band to obtainK-th narrowband components.(1)minukωk∑k=1K∂tδt+jπt∗uktexp−jωkt22,where uk and ωk are shorthand notations for the set of all modes and their center frequencies, and δt is the impulse unit function. In order to solve the constrained variational problem, a quadratic penalty term and the Lagrange multipliers are introduced to the augmented Lagrange method ζ, as shown in (2). The solution to the original minimization problem is now transformed to find the saddle point of augmented Lagrangian by a sequence of iterative suboptimization which is called alternate direction method of multipliers (ADMM).(2)ζuk,ωk,λ=α∑k=1K∂tδt+jπt∗uktexp−jωkt22+ft−∑k=1Kukt22+λt,ft−∑k=1Kukt.Additionally,α is the penalty parameter, λt are the Lagrange multipliers, and the detailed process of VMD has been fully described in [22]. The original signal xt has been decomposed into K-th BIMF components, ukt, and residual signals, εt, with lower energy.(3)xt=∑k=1Kukt+εt.The VMD achieves signal decomposition by iteratively searching for the optimal solution of the variational model in variational problems. VMD keeps the advantages of effective suppression on cross-termslnAkt. Therefore, it has robustness, and it overcomes the shortages of EMD algorithm, including mode aliasing, pseudocomponents, and vulnerability to acoustics interference.


### 2.2. System Modal Parameters

In order to effectively evaluate the mode parameter of VMD, the concept of modal parameters is introduced into VMD. The natural frequency and the damping ratio of modal parameters are calculated. Assuming a continuous signal in the time domain,ut=atcosϕt, where at represents the amplitude of signals and ϕt is the phase angle. The conjugate signal u^t can be obtained by the Hilbert transformation. The constructed analytic signal ht formed by ut and u^t can be shown as (4)ht=ut+ju^t=atexpjθt.Besides,(5)at=u2t+u^2t1/2,θt=arctanu^tut.Therefore, the instantaneous frequency can be defined as the derivative of phase function:(6)ω=dθtdt=1a2uu^′−u^u′,ft=12πdθtdt.For small damping systems, the damping ratio and the natural frequency [23] can be calculated according to the attenuation amplitude of modal response obtained through the following relationship between modal response and time t.(7)lnAkt=−ξkωn,kt+c,ωn,k=dωktt.The estimated modal parameters of the natural frequenciesωn,k and the damping ratio ξk are fitted by using the linear least-squares method.


### 2.3. Fast Independent Component Analysis

Independent component analysis (ICA) is a novel data analysis method that has been used in recent years [24]. It has different applications in several data analysis problems. ICA is developed for cocktail party problems [25]. It involves separating the mixed signals produced by two or more speakers talking simultaneously. The ICA model can be briefly described as(8)xt=At∗st+vt,where st is the unknown signals, which represent the vibration and the corrupting vibration signals. The sources transmitted and detected by n-th sensors. xt are the observed signals, and At is the linear random mixing matrix that represents the propagation from source to sensor.vt is the added noise. It is assumed that the source signal St can be transformed by random matrix A to X(t). Therefore, Xt=ASt and yt=WXt=WASt. Additionally, WA=I and yt=St. Therefore, yt and St can be similar by iterative selection and optimization of the object function to restore the observed signal to the source signal.In this paper, the FastICA algorithm based on negative entropy is adopted. The principle of calculation is to find a direction through system learning, so thatW can obtain a weighted vector wi and map as yi=wiTx with the greatest Gaussianity. The formula taken by negative entropy is shown as (9), where G is an arbitrary quadratic function, yi is a random variable with zero mean and unit variance, and v is a Gauss random variable with zero mean and unit variance. Substituting yi=wiTx into the negative entropy approximation gives (9)Jyi∝EGyi−EGv2,Jwi∝EGwiTz−EGv2.Equation (17) is the object function of FastICA based on negative entropy. FastICA algorithm solves the optimal solution of row vector by maximizing the above objective function. In order to maximize the negative entropy Jwi in (9), the mean value of random variables is zero and the variance is 1. Therefore, (9), evolves into EGwiTz and reaches the maximum when wi=1. Lagrange multiplier method can be obtained: (10)Lwi=EGwiTx−βwi.The first derivative function of (10) is (11)L′wi=∂Lwi∂wi=ExgwiTx−βwi.The quadratic function of the derivation for (10) is (12)L″wi=∂L′wi∂wi=ExxTg′wiTx−βI,where g′ is the derivative of g.(13)ExxTg′wiTx≈ExxT⋅Eg′wiTx=Eg′wiTxI.Solving equationL′wi=0 by Newton’s iteration method and determining the approximate solutions, we get(14)wik+1=wik−ExgwiTx−βwikEg′wiTx−β.Considering the conditions of normalization weight vectors,wi=1 and β=EwiTxgwiTx, using them as the input into (14), and multiplying both sides by β−Eg′wiTx, we obtain the iteration equation as(15)wik+1=ExgwiTx−Eg′wiTxwik.


### 2.4. The Proposed VMD-FastICA-HT Method

The parameter evaluation and VMD-FastICA-HT method proposed in this paper for the process of acoustics source identification is shown in the flowchart in Figure1; the core of this method consists of three parts: (1)
Select the penalty parameterα from the energy of signal spectra.(2)
Select the mode parameterK from the steady-state correlation of system modal parameters (i.e., natural frequency and damping ratio).(3)
Apply FastICA to the band-limited intrinsic mode function (BIMF) of all extended measurement channels, and use the Hilbert transform to estimate and recognize the source signals.Figure 1
The process of acoustics source identification.


## 2.1. Variational Mode Decomposition

The VMD method abandons the process of iterative screening of the extrema of EMD method, and it introduces signal decomposition into variational model to solve the decomposition problem. In order to find the optimal solution, signal decomposition is realized by seeking for a constrained variational model. In the process of VMD decomposition, the center frequency and bandwidth of each BIMF component are alternately iterated and adaptively decomposed to the appropriate signal frequency band to obtainK-th narrowband components.(1)minukωk∑k=1K∂tδt+jπt∗uktexp−jωkt22,where uk and ωk are shorthand notations for the set of all modes and their center frequencies, and δt is the impulse unit function. In order to solve the constrained variational problem, a quadratic penalty term and the Lagrange multipliers are introduced to the augmented Lagrange method ζ, as shown in (2). The solution to the original minimization problem is now transformed to find the saddle point of augmented Lagrangian by a sequence of iterative suboptimization which is called alternate direction method of multipliers (ADMM).(2)ζuk,ωk,λ=α∑k=1K∂tδt+jπt∗uktexp−jωkt22+ft−∑k=1Kukt22+λt,ft−∑k=1Kukt.Additionally,α is the penalty parameter, λt are the Lagrange multipliers, and the detailed process of VMD has been fully described in [22]. The original signal xt has been decomposed into K-th BIMF components, ukt, and residual signals, εt, with lower energy.(3)xt=∑k=1Kukt+εt.The VMD achieves signal decomposition by iteratively searching for the optimal solution of the variational model in variational problems. VMD keeps the advantages of effective suppression on cross-termslnAkt. Therefore, it has robustness, and it overcomes the shortages of EMD algorithm, including mode aliasing, pseudocomponents, and vulnerability to acoustics interference.


## 2.2. System Modal Parameters

In order to effectively evaluate the mode parameter of VMD, the concept of modal parameters is introduced into VMD. The natural frequency and the damping ratio of modal parameters are calculated. Assuming a continuous signal in the time domain,ut=atcosϕt, where at represents the amplitude of signals and ϕt is the phase angle. The conjugate signal u^t can be obtained by the Hilbert transformation. The constructed analytic signal ht formed by ut and u^t can be shown as (4)ht=ut+ju^t=atexpjθt.Besides,(5)at=u2t+u^2t1/2,θt=arctanu^tut.Therefore, the instantaneous frequency can be defined as the derivative of phase function:(6)ω=dθtdt=1a2uu^′−u^u′,ft=12πdθtdt.For small damping systems, the damping ratio and the natural frequency [23] can be calculated according to the attenuation amplitude of modal response obtained through the following relationship between modal response and time t.(7)lnAkt=−ξkωn,kt+c,ωn,k=dωktt.The estimated modal parameters of the natural frequenciesωn,k and the damping ratio ξk are fitted by using the linear least-squares method.


## 2.3. Fast Independent Component Analysis

Independent component analysis (ICA) is a novel data analysis method that has been used in recent years [24]. It has different applications in several data analysis problems. ICA is developed for cocktail party problems [25]. It involves separating the mixed signals produced by two or more speakers talking simultaneously. The ICA model can be briefly described as(8)xt=At∗st+vt,where st is the unknown signals, which represent the vibration and the corrupting vibration signals. The sources transmitted and detected by n-th sensors. xt are the observed signals, and At is the linear random mixing matrix that represents the propagation from source to sensor.vt is the added noise. It is assumed that the source signal St can be transformed by random matrix A to X(t). Therefore, Xt=ASt and yt=WXt=WASt. Additionally, WA=I and yt=St. Therefore, yt and St can be similar by iterative selection and optimization of the object function to restore the observed signal to the source signal.In this paper, the FastICA algorithm based on negative entropy is adopted. The principle of calculation is to find a direction through system learning, so thatW can obtain a weighted vector wi and map as yi=wiTx with the greatest Gaussianity. The formula taken by negative entropy is shown as (9), where G is an arbitrary quadratic function, yi is a random variable with zero mean and unit variance, and v is a Gauss random variable with zero mean and unit variance. Substituting yi=wiTx into the negative entropy approximation gives (9)Jyi∝EGyi−EGv2,Jwi∝EGwiTz−EGv2.Equation (17) is the object function of FastICA based on negative entropy. FastICA algorithm solves the optimal solution of row vector by maximizing the above objective function. In order to maximize the negative entropy Jwi in (9), the mean value of random variables is zero and the variance is 1. Therefore, (9), evolves into EGwiTz and reaches the maximum when wi=1. Lagrange multiplier method can be obtained: (10)Lwi=EGwiTx−βwi.The first derivative function of (10) is (11)L′wi=∂Lwi∂wi=ExgwiTx−βwi.The quadratic function of the derivation for (10) is (12)L″wi=∂L′wi∂wi=ExxTg′wiTx−βI,where g′ is the derivative of g.(13)ExxTg′wiTx≈ExxT⋅Eg′wiTx=Eg′wiTxI.Solving equationL′wi=0 by Newton’s iteration method and determining the approximate solutions, we get(14)wik+1=wik−ExgwiTx−βwikEg′wiTx−β.Considering the conditions of normalization weight vectors,wi=1 and β=EwiTxgwiTx, using them as the input into (14), and multiplying both sides by β−Eg′wiTx, we obtain the iteration equation as(15)wik+1=ExgwiTx−Eg′wiTxwik.


## 2.4. The Proposed VMD-FastICA-HT Method

The parameter evaluation and VMD-FastICA-HT method proposed in this paper for the process of acoustics source identification is shown in the flowchart in Figure1; the core of this method consists of three parts: (1)
Select the penalty parameterα from the energy of signal spectra.(2)
Select the mode parameterK from the steady-state correlation of system modal parameters (i.e., natural frequency and damping ratio).(3)
Apply FastICA to the band-limited intrinsic mode function (BIMF) of all extended measurement channels, and use the Hilbert transform to estimate and recognize the source signals.Figure 1
The process of acoustics source identification.


## 3. Simulation Studies

### 3.1. Simulation Signals Construction

Before selecting the appropriate VMD parameters, the effects of penalty parameters and mode parameters will be discussed. The core of VMD is the assumption that different modal signals have limited bandwidth. Therefore, the number of modes has a significant impact on the decomposition accuracy. Without loss of generality, the influence of penalty parameters is discussed when mode parameterK is fixed. Three simulated signals are constructed for validation, as shown in Figure 2. By setting the time interval as 0.001 s and the signal length as 1000, the generated source signal waveforms are as follows: (16)s1t=sin120t,s2t=2sin20tsin500t,s3t=1+0.5sin30tsin1500t,ft=s1t+s2t+s3t+η.Figure 2
Time-domain comparison of BIMF after VMD.
(a)(b)(c)Because of the high signal-to-noise ratio, the background noise can be neglected in this scenario. The convex function of the optimal Lagrange multiplier and the quadratic penalty parameter can be set to zero. The number of source signals is 3.The simulated signalft is processed in VMD, where the penalty parameter is set to 1000 with K=3. The estimated BIMFs and the three simulated source signals are shown in Figure 2. In terms of the magnitude and phases, the first two BIMF components recover the source signal, whereas the magnitude of BIMF3 is different from that of the source signal S3.The Fourier spectra for BIMFs and the source signals are shown in Figure3. Since the simulated source signals are sinusoidal, the main frequency obtained in the BIMF is consistent with that of the sinusoidal signal. Additionally, it matches the three peaks of the simulated signals in the frequency domain. It can be proved that the simulated signal can be well reconstructed from the decomposition. In addition, the estimated low-frequency source signal is almost perfectly restored after VMD decomposition, and the center-frequency signal can also be restored satisfactorily.Figure 3
Comparison of BIMF frequency between source signal and VMD decomposition.


### 3.2. The Evaluation Criterion of VMD Penalty Parameter Selection

In general, it can be concluded that selecting the mode parameters corresponding to the number of source signals is a good choice, when the signal-to-noise ratio (SNR) is high. However, the penalty parametera needs further study, and the error of two norms ek is used to compare the coincidence of each BIMF component with its corresponding source signal xk.(17)ek=xkt−ukt2xkt2.Apart from this, asα changes, the frequency difference eΔf between the center frequency of BIMF ωk and the source signal xk can be represented by (18)eΔf=ωk−xkxk.It can be seen from Figure4 that the center frequency of initialization is uniformly distributed within the analysis frequency range. According to the decomposition, the smaller the two-norm error is, the higher the degree of reduction is. It is not difficult to see that when the error of center frequency eΔf is less than 1%, the center frequency is closer to the third harmonic frequency with the increase of α. For medium- and high-frequency energy, the change of two norms is small, whereas the two norms’ error is much smaller for low-energy high-frequency signals with the increase of α. For high-frequency low-energy signals, the penalty parameter needs to be larger. Therefore, when the penalty parameter is too small, it will lead to mode aliasing. Meanwhile, if the penalty parameter is too large, it will lead to frequency separation or overlap.Figure 4
The performances of VMD decomposition results under different penalty parameters.
(a)(b)VMD can essentially be treated as a Wiener filter, with the penalty parameterα being related to the signal and noise energy. For BIMF, the penalty parameter α should be large enough near the center frequency, and the energy in a narrow band can be rapidly attenuated to achieve low-pass filtering. In order to verify that the parameter α is related to energy and eliminate accidental factors, 1000 simulated samples of test signals have been constructed based on (16) and the signal-to-noise ratio is adjusted to a value between 15 dB and 40 dB. Under the condition of K = 3, VMD decomposition is carried out for each test signal with different penalty parameters α. The box diagram in Figure 5 shows the size of the required α, when the frequency error is less than 1%.Figure 5
The box diagram of the penalty parameterα when the frequency error is less than 1%.According to the statistics in Figure5, it is found that all the center frequencies coincide with the third harmonic frequencies. The approximate optimal decomposition results can be obtained. When the penalty parameters are greater than 300, 75% of the results can be guaranteed and the frequency error is less than 1%.Better decomposition results are obtained, with the increase of penalty parameters. 25% of center frequencies begin to change when the penalty parameter is less than 8000. Therefore, 1.5 times of upper quartile for lower limit of penalty parameter is selected as the bottom limit (red line). The upper limit is the lower quartile of maximum penalty parameter (blue line). Therefore, the optimal penalty parameter intervalα∈800,8000 of the third harmonic is obtained. Based on the above results, this paper gives the selection formula of penalty parameter α based on the energy of spectra.(19)α=Bmaxℱft,where ℱft denotes the conjugate multiplication of Fourier transform for the signal ft. According to the Fourier transform shown in Figure 3, the conjugate multiplication results can be calculated and α∈800,8000. B can be obtained around 1.5–16. Consequently, when VMD is applied, the penalty parameter α can be set by calculating the magnitude of energy in the frequency domain after the Fourier transform.


### 3.3. The Evaluation Criterion of VMD Mode ParameterK

#### 3.3.1. The Selection of Mode ParameterK Based on Structural System Parameters

In order to further discuss the selection of mode parameterK based on system modal parameters in VMD, a simple model is used to numerically verify the effectiveness and accuracy of the proposed method in system identification. In this paper, a four-degree-of-freedom numerical model has been created, as shown in Figure 6.Figure 6
The 4-DOF numerical model.The mass matrix of the model isM:(20)M=m1m2m3m4,where m1=m2=m3=m4=4.9895kg. The model stiffness matrix is K:(21)K=k1+k2−k2−k2k2+k3−k3−k3k3+k4−k4−k4k4.The time-invariant parameters in the time-varying model still exist:k2=1576N/m, k3=1226N/m, and k4=1051N/m. When it is treated as the linear time constant, the parameter k1=1401N/m. When it is set to time-varying, the parameter k1=14011+0.3sin4πtN/m. In order to simplify the model, the proportional damping is used to verify C=aM+bk and set a=0.16,b=3.9188×10−4. Assuming the force ft is applied to m3, we have(22)ft=500N,t=0,0,t≠0.Both of initial displacement and velocity are set to zero. The time domains of free vibration response and the corresponding power spectral density (PSD) ofx1t,x3t in the measurement system are shown in Figure 7.Figure 7
The comparison of the time domain and power spectral density in time-varying and time-invariant signals of the system.
(a)(b)(c)(d)Compared with the time-invariant models, the time-varying models have more frequency responses. Therefore, when using traditional time domain analysis or frequency domain analysis, it is difficult to accurately reflect the nonlinearity of time-varying structures.By using the penalty selection method described in the previous section and taking time-varying observed signals as inputs, the penalty parameter is firstly determined, and then the BIMF components under different mode parametersK are obtained from VMD. According to (7), the fitted natural frequency and damping are determined to adjust the appropriate frequencies in time-varying models. After these stages, the comparison between the natural frequency and damping ratio under different BIMFs components can be illustrated by forming poles shown in Figure 8. The parameter K is selected from the steady-state pole within the changes of frequency and damping ratio.Figure 8
The steady-state poles diagram of mode parameterK in the time-varying model. (a) DOF1. (b) DOF4.
(a)(b)In general, the natural frequencies of BIMF components are not consistent with the change of mode parameterK. We define each natural frequency as a pole, and three types of poles are introduced. The pole that first appeared is the original pole. As the decomposition layer increases, only the changes of natural frequency for the original pole are smaller than the natural frequency threshold. These poles are referred to as frequency poles. Moreover, the changes of both damping ratio and natural frequency in the original pole or frequency pole are smaller than the thresholds. These poles are called steady-state poles. The variation thresholds of natural frequency and damping ratio are set as 1% and 5%.According to these definitions and the values of thresholds, the steady-state poles are found from the lowest decomposition level of a row of poles. If the damping ratio is stable compared with the next level pole, the natural frequency and damping ratio of the pole are selected as the final identification results. It can be seen from Figure8 that several columns of steady-state poles are formed. With the increase of VMD mode parameter, the overdecomposition of frequency is unavoidable. According to the selection principle of steady-state poles, the minimum mode parameter K used for forming steady-state poles, its natural frequency, and damping ratio are selected as the optimal decomposition results. Moreover, the BIMFs calculated by this evaluation method might come from different mode parameter K. According to Figure 8(a), BIMF1, BIMF2, BIMF3, and BIMF4 can be chosen from K = 11, 6, 2, 3. Based on these BIMF components, the time-varying theoretical instantaneous frequency is obtained. It can be seen from Figure 9(a) that the identified instantaneous frequency is similar to the theoretical value. However, the higher frequency parts do not satisfy the Nyquist frequency in Hilbert transform, and the frequency boundaries of the decomposed components are inconsistent, due to endpoints effect.Figure 9
Identification of time-frequency spectrum and theoretical frequency by the numerical model in VMD and EMD methods. (a) VMD identification of time-frequency spectrum. (b) EMD identification of time-frequency spectrum.
(a)(b)In order to verify the effectiveness of VMD, EMD is used as a comparison to decompose the group of signal results and identify its instantaneous frequency in Figure9(b). According to the instantaneous frequency determined by EMD, only IMF2 and IMF3 are relatively close to theoretical values. Two instantaneous frequencies are decomposed into the same IMF component, due to the mode mixing. Additionally, the identified frequencies are higher than the theoretical values. Therefore, the EMD method cannot accurately identify these two instantaneous frequencies. Furthermore, the performance of EMD signal decomposition is easily disturbed by noise, which might lead to the uncertainty of observation results.


### 3.4. The Discussion of Mode ParameterK on FastICA

Considering the limitation of test conditions and the influence of noise, the number of measured mixed signals is less than that of source signals, and the noise should be added to the simulated mixed signals. Therefore, the source signals and observed signals as shown in Figure10 are constructed based on (16) and (8), and the mixing matrix A of FastICA models is defined as (23)A=0.30.20.10.30.40.90.40.80.7.Figure 10
The mixed signalX1 and X2 in time domain and spectrum domain. (a) Mixed signal X1. (b) X1 power spectrum. (c) Mixed signal X2. (d) X2 power spectrum.
(a)(b)(c)(d)Since the number of test signals is less than the number of source signals, only two source signals can be decomposed. In order to further illustrate the best results in FastICA obtained by using VMD parameters, the valuesK = 2, 3, 4, 5, 6, 7, and 8 have been selected in this section. The different mode parameters and the proposed setting method of the parameter are shown in Figure 11. VMD-FastICA is used to decompose two test signals, and the results of simulated and estimated source signal are compared. The estimated source signals whose correlation is the largest with the source signal are counted as shown in Table 1. It can be seen that although the values K = 2 and 6 are related to 95% or more, there are separate signals that are negatively correlated with the source signal. When the parameters K = 3, 5, 7, and 8, only two highly correlated high-separation signals can be found.Figure 11
The comparison results of estimated source signals and original signals with parameterK.
(a)(b)(c)Table 1
The correlation comparison of estimated source signals with mode parameterK.
Correlation factorsS1S2S3K = 2−0.998−0.999−0.999K = 30.974−0.0760.894K = 4 (selected from the proposed evaluation method)0.9980.9960.934K = 5−0.996−0.0130.907K = 6−0.995−0.9890.094K = 7−0.994−0.7090.071K = 80.9880.6950.011WhenK = 4 is selected from the system modal parameters shown in the previous section, the identified correlations are higher than 0.98 in the first two signals. Additionally, the other parameters’ results have the highest correlation with S3. It is shown that the parameter selected from the system modal parameter method has better performance in FastICA.


## 3.1. Simulation Signals Construction

Before selecting the appropriate VMD parameters, the effects of penalty parameters and mode parameters will be discussed. The core of VMD is the assumption that different modal signals have limited bandwidth. Therefore, the number of modes has a significant impact on the decomposition accuracy. Without loss of generality, the influence of penalty parameters is discussed when mode parameterK is fixed. Three simulated signals are constructed for validation, as shown in Figure 2. By setting the time interval as 0.001 s and the signal length as 1000, the generated source signal waveforms are as follows: (16)s1t=sin120t,s2t=2sin20tsin500t,s3t=1+0.5sin30tsin1500t,ft=s1t+s2t+s3t+η.Figure 2
Time-domain comparison of BIMF after VMD.
(a)(b)(c)Because of the high signal-to-noise ratio, the background noise can be neglected in this scenario. The convex function of the optimal Lagrange multiplier and the quadratic penalty parameter can be set to zero. The number of source signals is 3.The simulated signalft is processed in VMD, where the penalty parameter is set to 1000 with K=3. The estimated BIMFs and the three simulated source signals are shown in Figure 2. In terms of the magnitude and phases, the first two BIMF components recover the source signal, whereas the magnitude of BIMF3 is different from that of the source signal S3.The Fourier spectra for BIMFs and the source signals are shown in Figure3. Since the simulated source signals are sinusoidal, the main frequency obtained in the BIMF is consistent with that of the sinusoidal signal. Additionally, it matches the three peaks of the simulated signals in the frequency domain. It can be proved that the simulated signal can be well reconstructed from the decomposition. In addition, the estimated low-frequency source signal is almost perfectly restored after VMD decomposition, and the center-frequency signal can also be restored satisfactorily.Figure 3
Comparison of BIMF frequency between source signal and VMD decomposition.


## 3.2. The Evaluation Criterion of VMD Penalty Parameter Selection

In general, it can be concluded that selecting the mode parameters corresponding to the number of source signals is a good choice, when the signal-to-noise ratio (SNR) is high. However, the penalty parametera needs further study, and the error of two norms ek is used to compare the coincidence of each BIMF component with its corresponding source signal xk.(17)ek=xkt−ukt2xkt2.Apart from this, asα changes, the frequency difference eΔf between the center frequency of BIMF ωk and the source signal xk can be represented by (18)eΔf=ωk−xkxk.It can be seen from Figure4 that the center frequency of initialization is uniformly distributed within the analysis frequency range. According to the decomposition, the smaller the two-norm error is, the higher the degree of reduction is. It is not difficult to see that when the error of center frequency eΔf is less than 1%, the center frequency is closer to the third harmonic frequency with the increase of α. For medium- and high-frequency energy, the change of two norms is small, whereas the two norms’ error is much smaller for low-energy high-frequency signals with the increase of α. For high-frequency low-energy signals, the penalty parameter needs to be larger. Therefore, when the penalty parameter is too small, it will lead to mode aliasing. Meanwhile, if the penalty parameter is too large, it will lead to frequency separation or overlap.Figure 4
The performances of VMD decomposition results under different penalty parameters.
(a)(b)VMD can essentially be treated as a Wiener filter, with the penalty parameterα being related to the signal and noise energy. For BIMF, the penalty parameter α should be large enough near the center frequency, and the energy in a narrow band can be rapidly attenuated to achieve low-pass filtering. In order to verify that the parameter α is related to energy and eliminate accidental factors, 1000 simulated samples of test signals have been constructed based on (16) and the signal-to-noise ratio is adjusted to a value between 15 dB and 40 dB. Under the condition of K = 3, VMD decomposition is carried out for each test signal with different penalty parameters α. The box diagram in Figure 5 shows the size of the required α, when the frequency error is less than 1%.Figure 5
The box diagram of the penalty parameterα when the frequency error is less than 1%.According to the statistics in Figure5, it is found that all the center frequencies coincide with the third harmonic frequencies. The approximate optimal decomposition results can be obtained. When the penalty parameters are greater than 300, 75% of the results can be guaranteed and the frequency error is less than 1%.Better decomposition results are obtained, with the increase of penalty parameters. 25% of center frequencies begin to change when the penalty parameter is less than 8000. Therefore, 1.5 times of upper quartile for lower limit of penalty parameter is selected as the bottom limit (red line). The upper limit is the lower quartile of maximum penalty parameter (blue line). Therefore, the optimal penalty parameter intervalα∈800,8000 of the third harmonic is obtained. Based on the above results, this paper gives the selection formula of penalty parameter α based on the energy of spectra.(19)α=Bmaxℱft,where ℱft denotes the conjugate multiplication of Fourier transform for the signal ft. According to the Fourier transform shown in Figure 3, the conjugate multiplication results can be calculated and α∈800,8000. B can be obtained around 1.5–16. Consequently, when VMD is applied, the penalty parameter α can be set by calculating the magnitude of energy in the frequency domain after the Fourier transform.


## 3.3. The Evaluation Criterion of VMD Mode ParameterK

### 3.3.1. The Selection of Mode ParameterK Based on Structural System Parameters

In order to further discuss the selection of mode parameterK based on system modal parameters in VMD, a simple model is used to numerically verify the effectiveness and accuracy of the proposed method in system identification. In this paper, a four-degree-of-freedom numerical model has been created, as shown in Figure 6.Figure 6
The 4-DOF numerical model.The mass matrix of the model isM:(20)M=m1m2m3m4,where m1=m2=m3=m4=4.9895kg. The model stiffness matrix is K:(21)K=k1+k2−k2−k2k2+k3−k3−k3k3+k4−k4−k4k4.The time-invariant parameters in the time-varying model still exist:k2=1576N/m, k3=1226N/m, and k4=1051N/m. When it is treated as the linear time constant, the parameter k1=1401N/m. When it is set to time-varying, the parameter k1=14011+0.3sin4πtN/m. In order to simplify the model, the proportional damping is used to verify C=aM+bk and set a=0.16,b=3.9188×10−4. Assuming the force ft is applied to m3, we have(22)ft=500N,t=0,0,t≠0.Both of initial displacement and velocity are set to zero. The time domains of free vibration response and the corresponding power spectral density (PSD) ofx1t,x3t in the measurement system are shown in Figure 7.Figure 7
The comparison of the time domain and power spectral density in time-varying and time-invariant signals of the system.
(a)(b)(c)(d)Compared with the time-invariant models, the time-varying models have more frequency responses. Therefore, when using traditional time domain analysis or frequency domain analysis, it is difficult to accurately reflect the nonlinearity of time-varying structures.By using the penalty selection method described in the previous section and taking time-varying observed signals as inputs, the penalty parameter is firstly determined, and then the BIMF components under different mode parametersK are obtained from VMD. According to (7), the fitted natural frequency and damping are determined to adjust the appropriate frequencies in time-varying models. After these stages, the comparison between the natural frequency and damping ratio under different BIMFs components can be illustrated by forming poles shown in Figure 8. The parameter K is selected from the steady-state pole within the changes of frequency and damping ratio.Figure 8
The steady-state poles diagram of mode parameterK in the time-varying model. (a) DOF1. (b) DOF4.
(a)(b)In general, the natural frequencies of BIMF components are not consistent with the change of mode parameterK. We define each natural frequency as a pole, and three types of poles are introduced. The pole that first appeared is the original pole. As the decomposition layer increases, only the changes of natural frequency for the original pole are smaller than the natural frequency threshold. These poles are referred to as frequency poles. Moreover, the changes of both damping ratio and natural frequency in the original pole or frequency pole are smaller than the thresholds. These poles are called steady-state poles. The variation thresholds of natural frequency and damping ratio are set as 1% and 5%.According to these definitions and the values of thresholds, the steady-state poles are found from the lowest decomposition level of a row of poles. If the damping ratio is stable compared with the next level pole, the natural frequency and damping ratio of the pole are selected as the final identification results. It can be seen from Figure8 that several columns of steady-state poles are formed. With the increase of VMD mode parameter, the overdecomposition of frequency is unavoidable. According to the selection principle of steady-state poles, the minimum mode parameter K used for forming steady-state poles, its natural frequency, and damping ratio are selected as the optimal decomposition results. Moreover, the BIMFs calculated by this evaluation method might come from different mode parameter K. According to Figure 8(a), BIMF1, BIMF2, BIMF3, and BIMF4 can be chosen from K = 11, 6, 2, 3. Based on these BIMF components, the time-varying theoretical instantaneous frequency is obtained. It can be seen from Figure 9(a) that the identified instantaneous frequency is similar to the theoretical value. However, the higher frequency parts do not satisfy the Nyquist frequency in Hilbert transform, and the frequency boundaries of the decomposed components are inconsistent, due to endpoints effect.Figure 9
Identification of time-frequency spectrum and theoretical frequency by the numerical model in VMD and EMD methods. (a) VMD identification of time-frequency spectrum. (b) EMD identification of time-frequency spectrum.
(a)(b)In order to verify the effectiveness of VMD, EMD is used as a comparison to decompose the group of signal results and identify its instantaneous frequency in Figure9(b). According to the instantaneous frequency determined by EMD, only IMF2 and IMF3 are relatively close to theoretical values. Two instantaneous frequencies are decomposed into the same IMF component, due to the mode mixing. Additionally, the identified frequencies are higher than the theoretical values. Therefore, the EMD method cannot accurately identify these two instantaneous frequencies. Furthermore, the performance of EMD signal decomposition is easily disturbed by noise, which might lead to the uncertainty of observation results.


## 3.3.1. The Selection of Mode ParameterK Based on Structural System Parameters

In order to further discuss the selection of mode parameterK based on system modal parameters in VMD, a simple model is used to numerically verify the effectiveness and accuracy of the proposed method in system identification. In this paper, a four-degree-of-freedom numerical model has been created, as shown in Figure 6.Figure 6
The 4-DOF numerical model.The mass matrix of the model isM:(20)M=m1m2m3m4,where m1=m2=m3=m4=4.9895kg. The model stiffness matrix is K:(21)K=k1+k2−k2−k2k2+k3−k3−k3k3+k4−k4−k4k4.The time-invariant parameters in the time-varying model still exist:k2=1576N/m, k3=1226N/m, and k4=1051N/m. When it is treated as the linear time constant, the parameter k1=1401N/m. When it is set to time-varying, the parameter k1=14011+0.3sin4πtN/m. In order to simplify the model, the proportional damping is used to verify C=aM+bk and set a=0.16,b=3.9188×10−4. Assuming the force ft is applied to m3, we have(22)ft=500N,t=0,0,t≠0.Both of initial displacement and velocity are set to zero. The time domains of free vibration response and the corresponding power spectral density (PSD) ofx1t,x3t in the measurement system are shown in Figure 7.Figure 7
The comparison of the time domain and power spectral density in time-varying and time-invariant signals of the system.
(a)(b)(c)(d)Compared with the time-invariant models, the time-varying models have more frequency responses. Therefore, when using traditional time domain analysis or frequency domain analysis, it is difficult to accurately reflect the nonlinearity of time-varying structures.By using the penalty selection method described in the previous section and taking time-varying observed signals as inputs, the penalty parameter is firstly determined, and then the BIMF components under different mode parametersK are obtained from VMD. According to (7), the fitted natural frequency and damping are determined to adjust the appropriate frequencies in time-varying models. After these stages, the comparison between the natural frequency and damping ratio under different BIMFs components can be illustrated by forming poles shown in Figure 8. The parameter K is selected from the steady-state pole within the changes of frequency and damping ratio.Figure 8
The steady-state poles diagram of mode parameterK in the time-varying model. (a) DOF1. (b) DOF4.
(a)(b)In general, the natural frequencies of BIMF components are not consistent with the change of mode parameterK. We define each natural frequency as a pole, and three types of poles are introduced. The pole that first appeared is the original pole. As the decomposition layer increases, only the changes of natural frequency for the original pole are smaller than the natural frequency threshold. These poles are referred to as frequency poles. Moreover, the changes of both damping ratio and natural frequency in the original pole or frequency pole are smaller than the thresholds. These poles are called steady-state poles. The variation thresholds of natural frequency and damping ratio are set as 1% and 5%.According to these definitions and the values of thresholds, the steady-state poles are found from the lowest decomposition level of a row of poles. If the damping ratio is stable compared with the next level pole, the natural frequency and damping ratio of the pole are selected as the final identification results. It can be seen from Figure8 that several columns of steady-state poles are formed. With the increase of VMD mode parameter, the overdecomposition of frequency is unavoidable. According to the selection principle of steady-state poles, the minimum mode parameter K used for forming steady-state poles, its natural frequency, and damping ratio are selected as the optimal decomposition results. Moreover, the BIMFs calculated by this evaluation method might come from different mode parameter K. According to Figure 8(a), BIMF1, BIMF2, BIMF3, and BIMF4 can be chosen from K = 11, 6, 2, 3. Based on these BIMF components, the time-varying theoretical instantaneous frequency is obtained. It can be seen from Figure 9(a) that the identified instantaneous frequency is similar to the theoretical value. However, the higher frequency parts do not satisfy the Nyquist frequency in Hilbert transform, and the frequency boundaries of the decomposed components are inconsistent, due to endpoints effect.Figure 9
Identification of time-frequency spectrum and theoretical frequency by the numerical model in VMD and EMD methods. (a) VMD identification of time-frequency spectrum. (b) EMD identification of time-frequency spectrum.
(a)(b)In order to verify the effectiveness of VMD, EMD is used as a comparison to decompose the group of signal results and identify its instantaneous frequency in Figure9(b). According to the instantaneous frequency determined by EMD, only IMF2 and IMF3 are relatively close to theoretical values. Two instantaneous frequencies are decomposed into the same IMF component, due to the mode mixing. Additionally, the identified frequencies are higher than the theoretical values. Therefore, the EMD method cannot accurately identify these two instantaneous frequencies. Furthermore, the performance of EMD signal decomposition is easily disturbed by noise, which might lead to the uncertainty of observation results.


## 3.4. The Discussion of Mode ParameterK on FastICA

Considering the limitation of test conditions and the influence of noise, the number of measured mixed signals is less than that of source signals, and the noise should be added to the simulated mixed signals. Therefore, the source signals and observed signals as shown in Figure10 are constructed based on (16) and (8), and the mixing matrix A of FastICA models is defined as (23)A=0.30.20.10.30.40.90.40.80.7.Figure 10
The mixed signalX1 and X2 in time domain and spectrum domain. (a) Mixed signal X1. (b) X1 power spectrum. (c) Mixed signal X2. (d) X2 power spectrum.
(a)(b)(c)(d)Since the number of test signals is less than the number of source signals, only two source signals can be decomposed. In order to further illustrate the best results in FastICA obtained by using VMD parameters, the valuesK = 2, 3, 4, 5, 6, 7, and 8 have been selected in this section. The different mode parameters and the proposed setting method of the parameter are shown in Figure 11. VMD-FastICA is used to decompose two test signals, and the results of simulated and estimated source signal are compared. The estimated source signals whose correlation is the largest with the source signal are counted as shown in Table 1. It can be seen that although the values K = 2 and 6 are related to 95% or more, there are separate signals that are negatively correlated with the source signal. When the parameters K = 3, 5, 7, and 8, only two highly correlated high-separation signals can be found.Figure 11
The comparison results of estimated source signals and original signals with parameterK.
(a)(b)(c)Table 1
The correlation comparison of estimated source signals with mode parameterK.
Correlation factorsS1S2S3K = 2−0.998−0.999−0.999K = 30.974−0.0760.894K = 4 (selected from the proposed evaluation method)0.9980.9960.934K = 5−0.996−0.0130.907K = 6−0.995−0.9890.094K = 7−0.994−0.7090.071K = 80.9880.6950.011WhenK = 4 is selected from the system modal parameters shown in the previous section, the identified correlations are higher than 0.98 in the first two signals. Additionally, the other parameters’ results have the highest correlation with S3. It is shown that the parameter selected from the system modal parameter method has better performance in FastICA.


## 4. In Situ Validation Test

### 4.1. Diesel Engine Acoustics Test

A CAT-C18 diesel engine is mounted on a diesel locomotive. In order to eliminate the influence of other sound sources, the universal joint between the diesel engine and the gearbox and the baffle around the vehicle are removed. The vehicle is parked on the open track with no barrier within 15 m. The in situ test set-up is shown in Figure12.Figure 12
In situ test set-up.In the test process, the outdoor wind speed is recorded as 0.5 m/s and the background acoustics level is 45.6 dBA. The basic parameters of diesel engine are shown in Table2.Table 2
Basic parameters of the diesel engine.
CharacteristicsParametersEngine typeIn-line engineCylinder number4Maximum output power470 kWRated power speed2100 r/minTransmission ratio2.98Firing order1-5-3-6-2-4Maximum torque2137 N·mThe acoustic signal is measured by the DGO-9767-CD Electret Microphone with a frequency response ranging from 20 Hz to 20 kHz. A BBM’s acoustics analysis instrument is used to select a 140 dBA dynamic range for measurement and a sampling frequency of 32768 Hz. The diesel engine runs at a speed of 750 rpm without load. In this case, the diesel engine operates stably, making the acoustic source relatively easy to identify. The measurement positions are shown in Figure13.Figure 13
The measuring points of the diesel engine.The acoustic signals measured on both sides of the diesel engine sound level meter are shown in Figure14. Additionally, the energy spectra are shown in Figure 15. Under 750 RPM running conditions, the vibration signal of diesel engine has very obvious periodicity. For a six-cylinder four-stroke engine, the crankshaft rotates twice and each cylinder ignites once. Every 120°, three times at each turning point, each ignition is counted as an impulse, and the vibration frequency of the diesel engine is w=r∗3/60=37.25Hz. For a four-stroke engine, the exhaust acoustics are determined by the ignition frequency. Exhaust acoustics contain not only ignition frequency but also its harmonics. Figure 15 clearly shows the 37.5 Hz and 75 Hz which are associated with exhaust acoustics and match the vibration frequency calculated above. The gas acoustics and the structural acoustics are not obvious; this requires a signal processing method to separate and identify the acoustics sources.Figure 14
The measured acoustics signals of the diesel engine. (a) The sound pressure in measuring point 1. (b) The sound pressure in measuring point 2.
(a)(b)Figure 15
The power spectrum of diesel engine acoustics signals. (a) The power spectrum in measuring point 1. (b) The power spectrum in measuring point 2.
(a)(b)


### 4.2. Acoustics Source Separation and Identification in Real-Life Operating Diesel Engine

The first step is to use the VMD algorithm to decompose the acoustic signal. However, in order to determine the VMD decomposition parameters, the penalty parameterα can be found according to (19) and the power spectral density in Figure 14. Additionally, the steady-state poles of diesel engine acoustics are shown in Figure 16. The steady-state poles can be considered as frequency. The damping ratio does not change with the change of mode parameter K. The frequency pole means that the natural frequency is not changed, whereas the damping ratio is changed. Moreover, the original pole can be defined as the frequency and the damping ratio calculated from the measured signals.Figure 16
750 RPM diesel engine steady-state poles diagram. (a) Test point 1. (b) Test point 2.
(a)(b)According to Figure16(a), when K = 16 and 17, a steady-state pole with unchanged natural frequency and damping ratio is formed. As shown in Figure 16(b), the steady-state poles have been formed at K = 17. After the optimal modal parameters are determined, BIMF can be obtained by the variational mode decomposition algorithm. The variational modal component and the original acoustic signal form a new signal group. The FastICA algorithm is used to extract the independent components. Each obtained independent component corresponds to an acoustics source of diesel engine. The objective of the test is to separate combustion, structural, and exhaust acoustics. Combined with the Hilbert transform method and the necessary prior knowledge of diesel engine, the calculation results are further analyzed.IC1, IC2, and IC3 are components obtained by the VMD- FastICA method. The results of IC1 are shown in Figure17. It can be seen from Figure 17(b) that the frequency of IC1 component is concentrated around 36 Hz. The largest amplitude appears in 130° and 610° in Figure 17(a). Since the sensor is placed near cylinders 1 and 6, the energy is slightly larger than that of the other cylinders. Low-frequency acoustics are generated due to the ignition of each strong “blow” event. Due to the phase difference of sound wave transmission between cylinders 5 and 2, the maximum energy in the time spectrum is delayed by 30° in Figure 17(c). As mentioned before, the exhaust acoustic frequency component is determined by the ignition timing [1], and IC1 can be recognized as exhaust acoustics.Figure 17
The results of the IC1 component separated by the VMD-FastICA method. (a) Time-domain waveform of IC1. (b) Fourier spectrum of IC1. (c) Hilbert spectrum of IC1. (d) Continuous wavelet transform of IC1.
(a)(b)(c)(d)However, the IC2 component frequency is concentrated around 580 Hz. Figure18(a) shows the largest amplitude that appears at around 240° and 480°. The internal combustion engine is placed in the locomotive. Although some sound absorption treatment has been carried out, cylinder 1 is closer to the wall and has stronger reverberation. It is assumed that cylinder 1 is at 120°, and the ignition sequence should be 4-1-5-3-6-2. According to the previous experiment, the frequency acoustics can be attributed to the combustion acoustics caused by the change of pressure in the cylinder. These acoustics are added and mixed up with the mechanical acoustics of engine when they are transmitted to the outside. Because of the attenuation during the transmission process in IC2, the amplitude of IC1 is larger than that of IC2. Moreover, since the fundamental frequency of combustion acoustics usually does not exceed 250 Hz [26], it can be inferred that the fundamental frequency ratio of combustion efficiency acoustics is submerged in IC1, due to the masking effect. It is very difficult to identify. IC2 is recognized as the main frequency acoustics of combustion acoustics.Figure 18
The results of the IC2 component separated by the VMD-FastICA method. (a) Time-domain waveform of IC2. (b) Fourier spectrum of IC2. (c) Hilbert spectrum of IC2. (d) Continuous wavelet transform of IC2.
(a)(b)(c)(d)As can be seen from Figure19, the frequency of IC3 component is concentrated around 2400 Hz. These high-frequency acoustics are usually caused by the mechanical operation of various components for diesel engine, and they can also be treated as the results of piston acoustics caused by the piston pin moving up and down the cylinder bore of engine. In the process, it generates knock acoustics which transmits through the solid and engine shell and eventually radiates through the surface of the structure [27]. It is a kind of high-frequency acoustics with obvious low-frequency periodic modulation. Therefore, IC3 can be recognized as mechanical acoustics.Figure 19
The results of the IC3 component separated by the VMD-FastICA method. (a) Time domain waveform of IC3. (b) Fourier spectrum of IC3. (c) Hilbert spectrum of IC3. (d) Continuous wavelet transform of IC3.
(a)(b)(c)(d)The results of the continuous wavelet transform (CWT) are shown in Figures17(d), 18(d), and 19(d). Compared with the Hilbert method, CWT uses multiscale analysis. Because of the Heisenberg uncertainty principle, increasing the time resolution will reduce the frequency resolution, and vice versa. This will cause the low-frequency resolution and unclear identification of CWT method. Therefore, it can be considered that the VMD-FastICA-HT method has better performance than the VMD-FastICA-CWT method for a diesel engine. The CWT cannot clearly identify the frequency due to low-frequency resolution over a relatively long period time.


## 4.1. Diesel Engine Acoustics Test

A CAT-C18 diesel engine is mounted on a diesel locomotive. In order to eliminate the influence of other sound sources, the universal joint between the diesel engine and the gearbox and the baffle around the vehicle are removed. The vehicle is parked on the open track with no barrier within 15 m. The in situ test set-up is shown in Figure12.Figure 12
In situ test set-up.In the test process, the outdoor wind speed is recorded as 0.5 m/s and the background acoustics level is 45.6 dBA. The basic parameters of diesel engine are shown in Table2.Table 2
Basic parameters of the diesel engine.
CharacteristicsParametersEngine typeIn-line engineCylinder number4Maximum output power470 kWRated power speed2100 r/minTransmission ratio2.98Firing order1-5-3-6-2-4Maximum torque2137 N·mThe acoustic signal is measured by the DGO-9767-CD Electret Microphone with a frequency response ranging from 20 Hz to 20 kHz. A BBM’s acoustics analysis instrument is used to select a 140 dBA dynamic range for measurement and a sampling frequency of 32768 Hz. The diesel engine runs at a speed of 750 rpm without load. In this case, the diesel engine operates stably, making the acoustic source relatively easy to identify. The measurement positions are shown in Figure13.Figure 13
The measuring points of the diesel engine.The acoustic signals measured on both sides of the diesel engine sound level meter are shown in Figure14. Additionally, the energy spectra are shown in Figure 15. Under 750 RPM running conditions, the vibration signal of diesel engine has very obvious periodicity. For a six-cylinder four-stroke engine, the crankshaft rotates twice and each cylinder ignites once. Every 120°, three times at each turning point, each ignition is counted as an impulse, and the vibration frequency of the diesel engine is w=r∗3/60=37.25Hz. For a four-stroke engine, the exhaust acoustics are determined by the ignition frequency. Exhaust acoustics contain not only ignition frequency but also its harmonics. Figure 15 clearly shows the 37.5 Hz and 75 Hz which are associated with exhaust acoustics and match the vibration frequency calculated above. The gas acoustics and the structural acoustics are not obvious; this requires a signal processing method to separate and identify the acoustics sources.Figure 14
The measured acoustics signals of the diesel engine. (a) The sound pressure in measuring point 1. (b) The sound pressure in measuring point 2.
(a)(b)Figure 15
The power spectrum of diesel engine acoustics signals. (a) The power spectrum in measuring point 1. (b) The power spectrum in measuring point 2.
(a)(b)


## 4.2. Acoustics Source Separation and Identification in Real-Life Operating Diesel Engine

The first step is to use the VMD algorithm to decompose the acoustic signal. However, in order to determine the VMD decomposition parameters, the penalty parameterα can be found according to (19) and the power spectral density in Figure 14. Additionally, the steady-state poles of diesel engine acoustics are shown in Figure 16. The steady-state poles can be considered as frequency. The damping ratio does not change with the change of mode parameter K. The frequency pole means that the natural frequency is not changed, whereas the damping ratio is changed. Moreover, the original pole can be defined as the frequency and the damping ratio calculated from the measured signals.Figure 16
750 RPM diesel engine steady-state poles diagram. (a) Test point 1. (b) Test point 2.
(a)(b)According to Figure16(a), when K = 16 and 17, a steady-state pole with unchanged natural frequency and damping ratio is formed. As shown in Figure 16(b), the steady-state poles have been formed at K = 17. After the optimal modal parameters are determined, BIMF can be obtained by the variational mode decomposition algorithm. The variational modal component and the original acoustic signal form a new signal group. The FastICA algorithm is used to extract the independent components. Each obtained independent component corresponds to an acoustics source of diesel engine. The objective of the test is to separate combustion, structural, and exhaust acoustics. Combined with the Hilbert transform method and the necessary prior knowledge of diesel engine, the calculation results are further analyzed.IC1, IC2, and IC3 are components obtained by the VMD- FastICA method. The results of IC1 are shown in Figure17. It can be seen from Figure 17(b) that the frequency of IC1 component is concentrated around 36 Hz. The largest amplitude appears in 130° and 610° in Figure 17(a). Since the sensor is placed near cylinders 1 and 6, the energy is slightly larger than that of the other cylinders. Low-frequency acoustics are generated due to the ignition of each strong “blow” event. Due to the phase difference of sound wave transmission between cylinders 5 and 2, the maximum energy in the time spectrum is delayed by 30° in Figure 17(c). As mentioned before, the exhaust acoustic frequency component is determined by the ignition timing [1], and IC1 can be recognized as exhaust acoustics.Figure 17
The results of the IC1 component separated by the VMD-FastICA method. (a) Time-domain waveform of IC1. (b) Fourier spectrum of IC1. (c) Hilbert spectrum of IC1. (d) Continuous wavelet transform of IC1.
(a)(b)(c)(d)However, the IC2 component frequency is concentrated around 580 Hz. Figure18(a) shows the largest amplitude that appears at around 240° and 480°. The internal combustion engine is placed in the locomotive. Although some sound absorption treatment has been carried out, cylinder 1 is closer to the wall and has stronger reverberation. It is assumed that cylinder 1 is at 120°, and the ignition sequence should be 4-1-5-3-6-2. According to the previous experiment, the frequency acoustics can be attributed to the combustion acoustics caused by the change of pressure in the cylinder. These acoustics are added and mixed up with the mechanical acoustics of engine when they are transmitted to the outside. Because of the attenuation during the transmission process in IC2, the amplitude of IC1 is larger than that of IC2. Moreover, since the fundamental frequency of combustion acoustics usually does not exceed 250 Hz [26], it can be inferred that the fundamental frequency ratio of combustion efficiency acoustics is submerged in IC1, due to the masking effect. It is very difficult to identify. IC2 is recognized as the main frequency acoustics of combustion acoustics.Figure 18
The results of the IC2 component separated by the VMD-FastICA method. (a) Time-domain waveform of IC2. (b) Fourier spectrum of IC2. (c) Hilbert spectrum of IC2. (d) Continuous wavelet transform of IC2.
(a)(b)(c)(d)As can be seen from Figure19, the frequency of IC3 component is concentrated around 2400 Hz. These high-frequency acoustics are usually caused by the mechanical operation of various components for diesel engine, and they can also be treated as the results of piston acoustics caused by the piston pin moving up and down the cylinder bore of engine. In the process, it generates knock acoustics which transmits through the solid and engine shell and eventually radiates through the surface of the structure [27]. It is a kind of high-frequency acoustics with obvious low-frequency periodic modulation. Therefore, IC3 can be recognized as mechanical acoustics.Figure 19
The results of the IC3 component separated by the VMD-FastICA method. (a) Time domain waveform of IC3. (b) Fourier spectrum of IC3. (c) Hilbert spectrum of IC3. (d) Continuous wavelet transform of IC3.
(a)(b)(c)(d)The results of the continuous wavelet transform (CWT) are shown in Figures17(d), 18(d), and 19(d). Compared with the Hilbert method, CWT uses multiscale analysis. Because of the Heisenberg uncertainty principle, increasing the time resolution will reduce the frequency resolution, and vice versa. This will cause the low-frequency resolution and unclear identification of CWT method. Therefore, it can be considered that the VMD-FastICA-HT method has better performance than the VMD-FastICA-CWT method for a diesel engine. The CWT cannot clearly identify the frequency due to low-frequency resolution over a relatively long period time.


## 5. Conclusion

During the process of acoustic identification in diesel engines, the different sources can be identified by using the blind source separation and the EMD method. Because of the mode mixing of EMD, several frequency components will be decomposed into one component simultaneously. Not all acoustics of diesel engine can be identified. As a result, the VMD method based on blind source separation can effectively overcome mode mixing. However, the decomposition results of VMD depend on its mode parameter and penalty parameter. Therefore, this paper proposes an evaluation criterion of VMD mode parameter and penalty parameter based on system modal parameters and power spectra of signals. Firstly, the range of penalty parameters is selected based on the power of signals in the frequency domain. Then, the appropriate mode parameterK of VMD can be selected from the steady-state poles according to the percentage errors of both natural frequency and damping ratio. This method can determine the optimal combination of mode parameters and penalty parameters to be analyzed effectively.In addition, the measured acoustic signals contain different sources, including combustion acoustics, mechanical acoustics, and exhaust acoustics. Therefore, the multicomponent signal, processed by the VMD method, should be separated. In this paper, a method combining VMD, FastICA, and Hilbert transform has been applied to these multicomponent acoustic signals to separate and identify these different types of diesel engine sources. The results show that the VMD-FastICA-HT method can accurately separate and recognize combustion acoustics, mechanical acoustics, and exhaust acoustics of a practical diesel engine. In the future research, since the VMD method still has some shortages, including endpoint effects and long calculation time, it might not be suitable to analyze acoustic signals in real time. Additionally, the other nonlinear and nonstationary signals that cannot be solved by VMD need further research in the future. Moreover, due to the complex structure of internal combustion engine, the transmission path of acoustics in internal combustion engine will be deeply explored, so as to put forward reasonable active noise reduction strategies for combination acoustics, mechanical acoustics, and exhaust acoustics.



---

*Source: 8832932-2021-01-20.xml*