# An Improved VMD-Based Denoising Method for Time Domain Load Signal Combining Wavelet with Singular Spectrum Analysis

Authors: Jingjing Fu; Fuyou Cai; Yonghao Guo; Hongda Liu; Wentie Niu
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Measured load data play a crucial role in the fatigue durability analysis of mechanical structures. However, in the process of signal acquisition, time domain load signals are easily contaminated by noise. In this paper, a signal denoising method based on variational mode decomposition (VMD), wavelet threshold denoising (WTD), and singular spectrum analysis (SSA) is proposed. Firstly, a simple criterion based on mutual information entropy (MIE) is designed to select the proper mode number for VMD. Detrended fluctuation analysis (DFA) is adopted to obtain the noise level of the noisy signal, which can optimize the selection of MIE threshold. Meanwhile, the noisy signal is adaptively decomposed into band-limited intrinsic mode functions (BLIMFs) by using VMD. In addition, weighted-permutation entropy (WPE) is applied to divide the BLIMFs into signal-dominant BLIMFs and noise-dominant BLIMFs. Then, the signal-dominant BLIMFs are reconstructed with the noise-dominant BLIMFs processed by WTD. Finally, SSA is implemented for the reconstructed signal. Experimental results of synthetic signals demonstrate that the presented method outperforms the conventional digital signal denoising methods and the related methods proposed recently. Effectiveness of the proposed method is verified through experiments of the measured load signals.


---

## Body

## 1. Introduction

Fatigue failure is one of the main causes of mechanical failure. Thus, it is of great significance to analyze the fatigue durability of the mechanical structure [1]. The basis of durability design of agricultural vehicles is to obtain time domain load signals, which can reflect the real working condition. Usually, measured signals are nonstationary and contain a lot of noise due to the complex and changeable working environment, which makes it meaningless to apply measured signals directly to durability analysis. Therefore, it is of great significance to effectively remove the noise in time domain load signals for the durability analysis of agricultural machinery.There are many general signal denoising methods, whose basic idea is to extract and reconstruct useful signal components through some criteria [2, 3]. Common signal decomposition methods include wavelet transform (WT), empirical mode decomposition (EMD), and variational mode decomposition (VMD). WT can decompose the signal into different frequencies through multiscale analysis and have good time-frequency localization characteristics. Donoho and Johnstone [4] proposed a threshold denoising method based on WT and achieved good results. Therefore, it has been applied in many fields, such as machinery [5], biomedicine [6], and chemistry [7]. However, wavelet threshold denoising (WTD) has poor denoising performance on the signal with low signal-to-noise ratio. Because the measured signal usually contains a lot of background noise, the performance of WTD is not ideal [8]. EMD was put forward by Huang et al. [9] in 1998. Compared with WT, it can adaptively decompose the signal into a series of intrinsic mode functions (IMFs) without much prior information. Therefore, it has been widely used in signal denoising [10–12] and analysis [13, 14]. However, the defects of EMD, such as mode aliasing [15] and end effect [16], restrict its wide application. For this reason, scholars have proposed improved methods of EMD, such as EEMD [17] and CEEMD [18]. However, similar to EMD, the improved methods based on EMD are all iterative algorithms, which will inevitably produce the error accumulation phenomenon [19], and lack solid mathematical foundation [20].In 2014, Dragomiretskiy and Zosso [21] proposed the variational mode decomposition (VMD), which is a noniterative and adaptive signal processing method. Compared with the method based on EMD, it has a solid mathematical theoretical foundation, and it can avoid the error accumulation phenomenon. VMD-based denoising methods have been used in many fields, such as underwater acoustic signal [22], ship-radiated noise [23], seismic signal [24], and laser radar [25]. It is worth noting that the key problem of signal denoising with VMD is the determination of the mode number K and relevant modes. If we use trial and error to get the optimal mode number K, it will need dozens of operations and waste a lot of time. Therefore, the parameter values are usually determined based on experience and convenience, which seriously affects the denoising performance of VMD. For this reason, researchers have proposed many algorithms to select the mode number. For example, Liu et al. [26] determined the mode number by detrend fluctuation analysis (DFA) through a large number of simulation experiments and finally achieved a good denoising performance. Lei et al. [27] proposed a hybrid algorithm combining VMD and WTD and selected the mode number for VMD through the decomposition results of EMD. In addition, the proposed denoising algorithm is shown to be able to effectively recognize different cast speeds. Li et al. [28] proposed an improved VMD (IVMD), where the mode number is determined by a frequency-aided method. IVMD as a novel algorithm combining SE was first proposed for feature extraction of S-RN signals. At present, the main methods of selecting relevant modes are correlation coefficient [22], permutation entropy [29], approximate entropy [30], Hausdorff distance [31], Euclidean distance [32], etc. However, pure VMD denoising can not only remove the high-frequency noise but also reduce the effective high-frequency information. Some scholars have proposed to use the wavelet threshold method to remove noise components from high-frequency components and reconstruct the signal [27]. However, WTD cannot completely remove the high-frequency noise, and with the change of operation conditions, the low-frequency signal will be adulterated with low-frequency interharmonics [33], so the performance of the abovementioned denoising method is not ideal. Singular spectrum analysis (SSA) proposed by Broomhead and King [34] can decompose signals into a series of independent and meaningful principal components [35]. Therefore, SSA may be an effective method to identify useful signals from noise [36].In this paper, DFA is used to get the noise level of the noisy signal, and mutual information entropy (MIE) is used to determine the mode number for adaptive VMD denoising. Meanwhile, the advantages of WTD and SSA are effectively absorbed. Firstly, the noise level of the noisy signal is obtained by using DFA [37], and the mode number for VMD is obtained by using MIE. At the same time, the noisy signal is decomposed into a series of BLIMFs by VMD. Next, weighted-permutation entropy (WPE) is adopted to distinguish between signal-dominated BLIMFs and noise-dominated BLIMFs. Then, the signal-dominant BLIMFs are reconstructed with the noise-dominant BLIMFs, which have been denoised by WTD. Finally, SSA is implemented for the reconstructed signal. Simulation and experimental results show the effectiveness and superiority of the proposed method.The rest of this paper is organized as follows. The related theories, including VMD, WPE, WTD, and SSA, are introduced in Section2. The proposed denoising method is presented in Section 3. In Section 4, the proposed method is compared with other existing methods, and the superiority of the proposed method is proved. In Section 5, the proposed method is applied to measured load signals. Finally, conclusions are drawn in Section 6.


## 2. Methodology

### 2.1. Brief Description of the VMD Algorithm [21]

The VMD algorithm can decompose any signal into a series of BLIMFs by iteratively searching the optimal solution of the variational model. Each mode is compact around its respective center frequencyωk. In this algorithm, the intrinsic mode function (IMF) is redefined as an amplitude-modulated-frequency-modulated (AM-FM) signal, which can be written as(1)ukt=Aktcosϕkt,where phase ϕkt>0 is a nondecreasing function and Akt is the instantaneous amplitude of ukt.Both the instantaneous amplitudeAkt and the instantaneous frequency ωkt=ϕk′t vary slower than the phase ϕkt so that pseudostationarity can be assumed. In order to search for suitable ukt and ωkt, VMD algorithm is required to solve the following constrained variational problem:(2)minuk,ωk∑k=1K∂tδt+jπt∗ukte−jωkt22s.t.∑k=1Kuk=f,where uk and ωk are modes and their center frequency, respectively, K is the mode number, ∂t represents the gradient with respect to t, δt denotes the unit impulse function, ∗ is the convolution symbol, and f is the original signal to be decomposed.In order to render the problem unconstrained, a quadratic penaltyα and a Lagrangian multiplier λ are taken into consideration. The augmented Lagrangian is expressed as follows:(3)Luk,ωk,λ=α∑K∂tδt+jπt∗ukte−jωkt22+ft−∑Kukt22+λt,ft−∑Kukt.When the alternate direction method of multipliers (ADMM) is used to solve the problem,ukn+1, ωkn+1, and λn+1 are updated until the optimal solution of the variational constrained problem is obtained. The modal component u^kn+1 in the frequency domain can be expressed as(4)u^kn+1w=f^ω−∑i≠ku^iω+λ^ω/21+2αω−ωk2,where u^ is the Fourier transform of u. Equation (4) is the Fourier representation of a Wiener filter with power spectrum prior 1/ω−ωk2. It ensures that the VMD algorithm has a good noise robustness.Similarly, the center frequencyωkn+1 can be expressed as(5)ωkn+1ω=∫0∞ωu^kω2dω∫0∞u^kω2dω.The stopping condition of the iteration is(6)∑k=1Ku^kn+1−u^kn22u^kn22<e,where e is a given value.


### 2.2. A Parameter-Adaptive VMD Method

Compared with EMD and its improved algorithm, VMD has solid theoretical derivation and stronger robustness for data sampling and noise. However, two crucial parameters in the decomposition algorithm need to be resolved, i.e., the decomposition mode numberK and the quadratic penalty α. A relatively moderate value of α is recommended to be 2000 [21].In order to avoid the impact of overbinning or underbinning on the VMD denoising, the appropriate numberK must be selected. MIE [38] can be used to represent the statistical dependence of two random variables, so as to reflect the degree of correlation. The expression is addressed as follows:(7)MIEX,Y=∑i=1r∑j=1spxi,yjlbpxi,yjpxipyj,where X and Y represent different random variables, pxi,yj is the joint probability distribution, pxi and pyi are the marginal probability distributions, r and s are the length of X and Y, respectively, and lb denotes the binary logarithm.The stronger the correlation betweenX and Y, the greater MIE will be obtained. In this paper, The MIEs between BLIMFl l=1,2,…,K and the noisy signal are designed to select the number K. The MIE can be normalized as(8)δl=MIElmaxMIEl,where MIEl is the MIE between BLIMFl l=1,2,…,K and the noisy signal.When the minimum ofδll=1,2,…,K is less than or equal to a specified threshold value δ, it can be considered that the decomposed mode does not contain important information, and the number K will be taken as the optimal mode number of VMD. However, the more noise in the noisy signal, the greater the MIE between the noise-dominant BLIMFs and the noisy signal. Therefore, a single threshold cannot meet all signals with different noise levels. In order to widen the application scope of this method, DFA is introduced to optimize the selection of threshold through a large number of simulation experiments.The scaling exponentα0 [39] provides a quantitative measure of the temporal correlations that exists in the time series. It can also be viewed as an indicator that describes the “roughness” of the original time series: the larger the value of α0, the smoother the time series [39]. In this paper, through a large number of simulation experiments, the scaling exponent α0 is used to represent the noise level of the noisy signal. The appropriate threshold value is selected to determine the mode number of VMD, and the selection expression of the threshold value proposed in this paper is described as follows:(9)δ=0.3,α0<0.6,0.05,0.6≤α0<0.8,0.02,0.8≤α0<1.2,0.015,1.2≤α0<1.5,0.01,1.5≤α0<1.8,0.005,1.8≤α0.


### 2.3. Division of BLIMFs Based on WPE [40]

Permutation entropy (PE), proposed by Bandt and Pompe [41], can represent the complexity of time series. It is a nonlinear dynamic method to analyze the complexity of the system. Its algorithm is fast and easy to implement. However, PE simply arranges the components in the phase space vector, which results in the loss of amplitude information in the time series. The weighted-permutation entropy (WPE) proposed by Fadlallah et al. [40] can consider the amplitude information of time series on the basis of PE, which significantly improves the robustness and antinoise performance of PE. WPE can be described as follows:(1)
For a given time seriesztt=1n, the matrix of phase space reconstruction can be expressed as(10)Z=z1z1+τ⋯z1+m−1τz2z2+τ⋯z2+m−1τ⋮⋮⋮⋮zczc+τ⋯zn,whereτ is the time delay, m is the embedding dimension, and c is the number of phase space components, c=n−m−1τ.(2)
The weight of each row ofZ is calculated as(11)wi=1m∑j=1mzi+j−1τ−z¯i2,i=1,2,…,c,wherez¯i is the mean value of each row, z¯i=1/m∑i=1mzi+j−1τ.(3)
Each row ofZ can be arranged in ascending order:(12)zi+k1−1τ≤zi+k2−1τ≤⋯≤zi+kl−1τ≤⋯≤zi+km−1τ,wherekl is the column index of each element in the ith row.(4)
After that, each row is arranged in ascending order, and a group of permutations can be obtained. Since the embedding dimension ism, there will be m! possible permutations. We calculated the appearing number of every permutation na, 1≤a≤m!, and the weighted relative frequency of each permutation can be expressed as(13)pwa=wina∑i=1cωi.(5)
Then, the WPE can be designated as follows:(14)Hwm=−∑i=1m!pwalnpwa,and the WPE of order can be normalized as(15)hw=Hwmlnm!.The range ofhw is 0 to 1. Bandt and Pompe [41] recommended m = 3∼7 and τ = 1. Considering the computational efficiency, m = 3 is selected in this paper. Similar to PE [42], WPE can also be used to distinguish between signal-dominant BLIMFs and noise-dominant BLIMFs.


### 2.4. Wavelet Threshold Denoising Method [4]

The essence of WTD is to filter the signal in the wavelet domain. Firstly, the signal is transformed into the wavelet domain. Then, the noisy wavelet coefficients are threshold processed. Finally, the denoised signal is reconstructed by wavelet coefficients. The primary steps of WTD are as follows:(1)
Decomposing the original signal by wavelet transform (WT) with proper wavelet basis function and decomposition level, and a group of wavelet decomposition coefficients are obtained.(2)
Selecting the appropriate threshold value to process the wavelet decomposition coefficient.(3)
The wavelet coefficients are reconstructed to get the denoised signal.The selection of the threshold has a crucial impact on WTD. There are two common threshold denoising methods: wavelet hard threshold denoising (WHSD) and wavelet soft threshold denoising (WSTD) [4]. After WHSD, the signal is discontinuous, and the reconstructed signal will produce oscillation, while WSTD has better continuity. Therefore, WSTD is adopted in this paper.


### 2.5. Singular Spectrum Analysis

SSA proposed by Broomhead and King [34] is an effective method for analyzing and predicting nonlinear time series, which is suitable for nonlinear and nonstationary signal processing. For a given one-dimensional time series qtt=1N, after selecting the appropriate window length L, L∈2,N, the trajectory matrix can be obtained by using the data embedding method, which can be expressed as(16)Q=q1q2⋯qMq2q3⋯qM+1⋮⋮⋮⋮qLqL+1⋯qN,where M=N−L+1.Then, the singular value decomposition (SVD) ofQQTcan be written as(17)Q=∑i=1dQi=∑i=1dλiUiViT,where d is the number of nonzero singular values of QQT, Qi is the elementary matrices, λ1,λ2,…,λd is the singular values of QQT arranged in descending order, and Ui and Vi are the ith left and right eigen vectors of QQT, respectively.Based on the grouping principle, the set of indices1,2,…,d are partitioned into G irrelevant subsets I1, I2, …, IG. Let Ii=i1,i2,…,ic be the ith subset, the subset matrix and the trajectory matrix can be expressed as(18)QIi=Qi1+Qi2+⋯+Qic,Q=QI1+QI2+⋯+QIG.By applying the diagonal averaging method toQIi=qijL×M, the bth element of signal reconstruction components RCIi is given by(19)RCb=1b∑m=1bqm,b−m+1∗,1≤b<L∗,1L∗∑m=1L∗qm,b−m+1∗,L∗≤b<K∗,1N−b+1∑m=b−K∗+1N−K∗+1qm,b−m+1∗,K∗≤b<N,where L∗=minL,M, K∗=maxL,M, and(20)qij∗=qij,L<M,qji,M≤L.It is worth noting that the window lengthL should be carefully selected because it directly affects the decomposition. However, there is no universal rule for the selection of the window length L. Through a large number of simulation experiments, we find that when the window length Lis larger thanN/10, the denoising performance of the proposed method is good enough. Considering the computational efficiency, L=N/10 is adopted in this paper.Furthermore, the grouping principle is very important in SSA. In this paper, the noisy signal is divided into signal-dominant and noise-dominant components by singular entropy increment [36], and its computational equation is given by(21)Ii=−λi∑j=1dλjlnλi∑j=1dλj,i=1,2,…,d.


## 2.1. Brief Description of the VMD Algorithm [21]

The VMD algorithm can decompose any signal into a series of BLIMFs by iteratively searching the optimal solution of the variational model. Each mode is compact around its respective center frequencyωk. In this algorithm, the intrinsic mode function (IMF) is redefined as an amplitude-modulated-frequency-modulated (AM-FM) signal, which can be written as(1)ukt=Aktcosϕkt,where phase ϕkt>0 is a nondecreasing function and Akt is the instantaneous amplitude of ukt.Both the instantaneous amplitudeAkt and the instantaneous frequency ωkt=ϕk′t vary slower than the phase ϕkt so that pseudostationarity can be assumed. In order to search for suitable ukt and ωkt, VMD algorithm is required to solve the following constrained variational problem:(2)minuk,ωk∑k=1K∂tδt+jπt∗ukte−jωkt22s.t.∑k=1Kuk=f,where uk and ωk are modes and their center frequency, respectively, K is the mode number, ∂t represents the gradient with respect to t, δt denotes the unit impulse function, ∗ is the convolution symbol, and f is the original signal to be decomposed.In order to render the problem unconstrained, a quadratic penaltyα and a Lagrangian multiplier λ are taken into consideration. The augmented Lagrangian is expressed as follows:(3)Luk,ωk,λ=α∑K∂tδt+jπt∗ukte−jωkt22+ft−∑Kukt22+λt,ft−∑Kukt.When the alternate direction method of multipliers (ADMM) is used to solve the problem,ukn+1, ωkn+1, and λn+1 are updated until the optimal solution of the variational constrained problem is obtained. The modal component u^kn+1 in the frequency domain can be expressed as(4)u^kn+1w=f^ω−∑i≠ku^iω+λ^ω/21+2αω−ωk2,where u^ is the Fourier transform of u. Equation (4) is the Fourier representation of a Wiener filter with power spectrum prior 1/ω−ωk2. It ensures that the VMD algorithm has a good noise robustness.Similarly, the center frequencyωkn+1 can be expressed as(5)ωkn+1ω=∫0∞ωu^kω2dω∫0∞u^kω2dω.The stopping condition of the iteration is(6)∑k=1Ku^kn+1−u^kn22u^kn22<e,where e is a given value.


## 2.2. A Parameter-Adaptive VMD Method

Compared with EMD and its improved algorithm, VMD has solid theoretical derivation and stronger robustness for data sampling and noise. However, two crucial parameters in the decomposition algorithm need to be resolved, i.e., the decomposition mode numberK and the quadratic penalty α. A relatively moderate value of α is recommended to be 2000 [21].In order to avoid the impact of overbinning or underbinning on the VMD denoising, the appropriate numberK must be selected. MIE [38] can be used to represent the statistical dependence of two random variables, so as to reflect the degree of correlation. The expression is addressed as follows:(7)MIEX,Y=∑i=1r∑j=1spxi,yjlbpxi,yjpxipyj,where X and Y represent different random variables, pxi,yj is the joint probability distribution, pxi and pyi are the marginal probability distributions, r and s are the length of X and Y, respectively, and lb denotes the binary logarithm.The stronger the correlation betweenX and Y, the greater MIE will be obtained. In this paper, The MIEs between BLIMFl l=1,2,…,K and the noisy signal are designed to select the number K. The MIE can be normalized as(8)δl=MIElmaxMIEl,where MIEl is the MIE between BLIMFl l=1,2,…,K and the noisy signal.When the minimum ofδll=1,2,…,K is less than or equal to a specified threshold value δ, it can be considered that the decomposed mode does not contain important information, and the number K will be taken as the optimal mode number of VMD. However, the more noise in the noisy signal, the greater the MIE between the noise-dominant BLIMFs and the noisy signal. Therefore, a single threshold cannot meet all signals with different noise levels. In order to widen the application scope of this method, DFA is introduced to optimize the selection of threshold through a large number of simulation experiments.The scaling exponentα0 [39] provides a quantitative measure of the temporal correlations that exists in the time series. It can also be viewed as an indicator that describes the “roughness” of the original time series: the larger the value of α0, the smoother the time series [39]. In this paper, through a large number of simulation experiments, the scaling exponent α0 is used to represent the noise level of the noisy signal. The appropriate threshold value is selected to determine the mode number of VMD, and the selection expression of the threshold value proposed in this paper is described as follows:(9)δ=0.3,α0<0.6,0.05,0.6≤α0<0.8,0.02,0.8≤α0<1.2,0.015,1.2≤α0<1.5,0.01,1.5≤α0<1.8,0.005,1.8≤α0.


## 2.3. Division of BLIMFs Based on WPE [40]

Permutation entropy (PE), proposed by Bandt and Pompe [41], can represent the complexity of time series. It is a nonlinear dynamic method to analyze the complexity of the system. Its algorithm is fast and easy to implement. However, PE simply arranges the components in the phase space vector, which results in the loss of amplitude information in the time series. The weighted-permutation entropy (WPE) proposed by Fadlallah et al. [40] can consider the amplitude information of time series on the basis of PE, which significantly improves the robustness and antinoise performance of PE. WPE can be described as follows:(1)
For a given time seriesztt=1n, the matrix of phase space reconstruction can be expressed as(10)Z=z1z1+τ⋯z1+m−1τz2z2+τ⋯z2+m−1τ⋮⋮⋮⋮zczc+τ⋯zn,whereτ is the time delay, m is the embedding dimension, and c is the number of phase space components, c=n−m−1τ.(2)
The weight of each row ofZ is calculated as(11)wi=1m∑j=1mzi+j−1τ−z¯i2,i=1,2,…,c,wherez¯i is the mean value of each row, z¯i=1/m∑i=1mzi+j−1τ.(3)
Each row ofZ can be arranged in ascending order:(12)zi+k1−1τ≤zi+k2−1τ≤⋯≤zi+kl−1τ≤⋯≤zi+km−1τ,wherekl is the column index of each element in the ith row.(4)
After that, each row is arranged in ascending order, and a group of permutations can be obtained. Since the embedding dimension ism, there will be m! possible permutations. We calculated the appearing number of every permutation na, 1≤a≤m!, and the weighted relative frequency of each permutation can be expressed as(13)pwa=wina∑i=1cωi.(5)
Then, the WPE can be designated as follows:(14)Hwm=−∑i=1m!pwalnpwa,and the WPE of order can be normalized as(15)hw=Hwmlnm!.The range ofhw is 0 to 1. Bandt and Pompe [41] recommended m = 3∼7 and τ = 1. Considering the computational efficiency, m = 3 is selected in this paper. Similar to PE [42], WPE can also be used to distinguish between signal-dominant BLIMFs and noise-dominant BLIMFs.


## 2.4. Wavelet Threshold Denoising Method [4]

The essence of WTD is to filter the signal in the wavelet domain. Firstly, the signal is transformed into the wavelet domain. Then, the noisy wavelet coefficients are threshold processed. Finally, the denoised signal is reconstructed by wavelet coefficients. The primary steps of WTD are as follows:(1)
Decomposing the original signal by wavelet transform (WT) with proper wavelet basis function and decomposition level, and a group of wavelet decomposition coefficients are obtained.(2)
Selecting the appropriate threshold value to process the wavelet decomposition coefficient.(3)
The wavelet coefficients are reconstructed to get the denoised signal.The selection of the threshold has a crucial impact on WTD. There are two common threshold denoising methods: wavelet hard threshold denoising (WHSD) and wavelet soft threshold denoising (WSTD) [4]. After WHSD, the signal is discontinuous, and the reconstructed signal will produce oscillation, while WSTD has better continuity. Therefore, WSTD is adopted in this paper.


## 2.5. Singular Spectrum Analysis

SSA proposed by Broomhead and King [34] is an effective method for analyzing and predicting nonlinear time series, which is suitable for nonlinear and nonstationary signal processing. For a given one-dimensional time series qtt=1N, after selecting the appropriate window length L, L∈2,N, the trajectory matrix can be obtained by using the data embedding method, which can be expressed as(16)Q=q1q2⋯qMq2q3⋯qM+1⋮⋮⋮⋮qLqL+1⋯qN,where M=N−L+1.Then, the singular value decomposition (SVD) ofQQTcan be written as(17)Q=∑i=1dQi=∑i=1dλiUiViT,where d is the number of nonzero singular values of QQT, Qi is the elementary matrices, λ1,λ2,…,λd is the singular values of QQT arranged in descending order, and Ui and Vi are the ith left and right eigen vectors of QQT, respectively.Based on the grouping principle, the set of indices1,2,…,d are partitioned into G irrelevant subsets I1, I2, …, IG. Let Ii=i1,i2,…,ic be the ith subset, the subset matrix and the trajectory matrix can be expressed as(18)QIi=Qi1+Qi2+⋯+Qic,Q=QI1+QI2+⋯+QIG.By applying the diagonal averaging method toQIi=qijL×M, the bth element of signal reconstruction components RCIi is given by(19)RCb=1b∑m=1bqm,b−m+1∗,1≤b<L∗,1L∗∑m=1L∗qm,b−m+1∗,L∗≤b<K∗,1N−b+1∑m=b−K∗+1N−K∗+1qm,b−m+1∗,K∗≤b<N,where L∗=minL,M, K∗=maxL,M, and(20)qij∗=qij,L<M,qji,M≤L.It is worth noting that the window lengthL should be carefully selected because it directly affects the decomposition. However, there is no universal rule for the selection of the window length L. Through a large number of simulation experiments, we find that when the window length Lis larger thanN/10, the denoising performance of the proposed method is good enough. Considering the computational efficiency, L=N/10 is adopted in this paper.Furthermore, the grouping principle is very important in SSA. In this paper, the noisy signal is divided into signal-dominant and noise-dominant components by singular entropy increment [36], and its computational equation is given by(21)Ii=−λi∑j=1dλjlnλi∑j=1dλj,i=1,2,…,d.


## 3. VMD-WTD-SSA Denoising Procedure

The process of the new efficient denoising method proposed in this paper is shown in Figure1. It consists of four main steps explained as follows:(1)
The noise level of the noisy signal is obtained by using DFA. Then, MIE is used to determine the mode number for adaptive VMD denoising. At the same time, the noisy signal is decomposed into a series of BLIMFs by VMD, whose frequency increases gradually. The noise mainly concentrates in high frequency, and the useful signal is mainly concentrated in the first few BLIMFs.(2)
The WPE of each BLIMF is calculated. If the WPE is less than 0.4, the corresponding BLIMF will be regarded as the signal-dominant component. Otherwise, the corresponding BLIMF will be regarded as the noise-dominant component.(3)
The denoised signal by VMD is obtained by reconstructing the signal-dominant BLIMFs. Then, the denoised signal by VMD-WTD is obtained by reconstructing the signal-dominant BLIMFs and the noise-dominant BLIMFs, which have been denoised by WTD.(4)
The SSA is carried out for the denoised signal by VMD-WTD to obtain the denoised signal by VMD-WTD-SSA.Figure 1
Flow chart of the proposed denoising method.


## 4. Simulation Signal Denoising

### 4.1. Construction of Simulation Signal

In order to verify the performance of the method proposed in this paper, the simulation data are tested. The equation of the noisy signal is as follows:(22)xt=cos2πf1t+0.3cos2πf2t+0.02cos2πf3t,x′t=xt+n′,where f1=5, f2=20, and f3=150 represent the three frequencies of noise-free signal xt, and n′ is the Gaussian white noise with a signal-to-noise ratio of 10 dB. The sampling frequency is 1000 Hz. The noisy and noise-free signals are shown in Figure 2.Figure 2
Simulation signal.


### 4.2. Denoising of Simulation Signal

The scaling exponentα0 of the noisy signal is 1.22, and the normalized MIE threshold δ is 0.015 by formula (9). Then, the normalized MIE between each BLIMF and the noisy signal (as defined in equation (8)) is obtained, as shown in Table 1.Table 1
Normalized MIE between each component and the noisy signal.
VMD
MIE
BLIMF1
1
BLIMF2
0.1336
BLIMF3
0.0119
BLIMF4
0.0223
BLIMF5
0.0224
BLIMF6
0.0333
BLIMF7
0.0187
BLIMF8
0.0178
BLIMF9
0.0180
BLIMF10
0.0257It can be seen that the optimal mode numberK is 10, and VMD is carried out for the simulation signal, and the BLIMFs are shown in Figure 3.Figure 3
Decomposition results of the noisy signal. (a) BLIMF1. (b) BLIMF2. (c) BLIMF3. (d) BLIMF4. (e) BLIMF5. (f) BLIMF6. (g) BLIMF7. (h) BLIMF8. (I) BLIMF9. (j) BLIMF10.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(I)
(j)The WPE of each BLIMF is shown in Table2. Obviously, BLIMF1 and BLIMF2 are signal-dominant BLIMFs, and the rest of the BLIMFs are noise-dominant BLIMFs. Then, the BLIMF1 and BLIMF2 are reconstructed to get the denoised signal by VMD, as shown in Figure 4. It can be seen that after VMD denoising, the signal appears with partial distortion and oscillation.Table 2
WPE of BLIMFs.
VMD
WPE
BLIMF1
0.3870
BLIMF2
0.3884
BLIMF3
0.4629
BLIMF4
0.5748
BLIMF5
0.7058
BLIMF6
0.8608
BLIMF7
0.9678
BLIMF8
0.9934
BLIMF9
0.8672
BLIMF10
0.7787Figure 4
Comparison of the noise-free signal and the denoised signal by VMD.WTD is conducted for the noise-dominant BLIMFs. In this paper, the db3 wavelet is selected and the number of decomposition levels is 3. Signal-dominant BLIMFs and noise-dominant BLIMFs after WTD are combined for signal reconstruction. VMD-WTD denoising performance is as shown in Figure5. Obviously, the partial distortion and oscillation still exist, which demonstrates that the signal denoising performance after VMD-WTD denoising is still not ideal.Figure 5
Comparison of the noise-free signal and the denoised signal by VMD-WTD.SSA is implemented for the denoised signal by VMD-WTD, and Figure6 shows singular entropy increments ΔE of the signal. It is found that the information quantity of the signal reaches the saturated state at order 4, so the first 4 elementary matrices are used for reconstruction. After diagonal averaging, the denoised signal by VMD-WTD-SSA is obtained, as shown in Figure 7. It can be seen that the signal waveform is well recovered after denoising. Through the comparison, it can be seen that VMD-WTD-SSA denoising performance is superior to VMD and VMD-WTD denoising performances obviously.Figure 6
Singular entropy increments of the denoised signal by VMD-WTD.Figure 7
Comparison of the noise-free signal and the denoised signal by VMD-WTD-SSA.


### 4.3. Comparison with Other Methods

In this paper, the signal-to-noise ratio (SNR), root mean square error (RMSE), and maximum absolute error (MAE) are used as the evaluation indexes of denoising performance. SNR reflects the performance of the denoising method, and the larger the value, the better the denoising performance; RMSE reflects the similarity between the denoised signal and the noise-free signal, and the smaller the value, the better the denoising performance; MAE reflects the real error between the denoised signal and the noise-free signal, and the smaller the value, the better the denoising performance. The three evaluation indexes are defined as follows:(23)SNR=10log∑i=1ny2t∑i=1nyt−yt2,RMSE=1n∑i=1nyt−y′t2,MAE=1n∑i=1nyt−y′t,where y (t) is the original signal, y′t is the denoised signal, and n is the length of the original signal.Six typical denoising methods are selected to compare with VMD-WTD-SSA. They are wavelet soft threshold denoising (WSTD) [4], EMD hard threshold denoising (EMD-HT) [11], EMD partial reconstruction (EMD-PR) [12], DFA-VMD [26], VMD-WTD [27], and VMD-AE [30]. Denoising results of different methods for the noisy signals with SNRin values varying from −10 dB to 20 dB are shown in Figure 8. It can be seen that compared with the other methods, VMD-WTD-SSA has higher SNRout and lower RMSE and MAE. Therefore, the results show that the VMD-WTD-SSA proposed in this paper has the best denoising performances with different SNRin values. The presented method is better than the other methods on the whole, and its performance on noisy signals with low SNR is more prominent. Besides, as the SNRin values continuously increases, the differences among these methods are narrowed.Figure 8
Denoising results for simulation signals with different SNRin: (a) SNRout, (b) RMSE, and (c) MAE.
(a)
(b)
(c)


## 4.1. Construction of Simulation Signal

In order to verify the performance of the method proposed in this paper, the simulation data are tested. The equation of the noisy signal is as follows:(22)xt=cos2πf1t+0.3cos2πf2t+0.02cos2πf3t,x′t=xt+n′,where f1=5, f2=20, and f3=150 represent the three frequencies of noise-free signal xt, and n′ is the Gaussian white noise with a signal-to-noise ratio of 10 dB. The sampling frequency is 1000 Hz. The noisy and noise-free signals are shown in Figure 2.Figure 2
Simulation signal.


## 4.2. Denoising of Simulation Signal

The scaling exponentα0 of the noisy signal is 1.22, and the normalized MIE threshold δ is 0.015 by formula (9). Then, the normalized MIE between each BLIMF and the noisy signal (as defined in equation (8)) is obtained, as shown in Table 1.Table 1
Normalized MIE between each component and the noisy signal.
VMD
MIE
BLIMF1
1
BLIMF2
0.1336
BLIMF3
0.0119
BLIMF4
0.0223
BLIMF5
0.0224
BLIMF6
0.0333
BLIMF7
0.0187
BLIMF8
0.0178
BLIMF9
0.0180
BLIMF10
0.0257It can be seen that the optimal mode numberK is 10, and VMD is carried out for the simulation signal, and the BLIMFs are shown in Figure 3.Figure 3
Decomposition results of the noisy signal. (a) BLIMF1. (b) BLIMF2. (c) BLIMF3. (d) BLIMF4. (e) BLIMF5. (f) BLIMF6. (g) BLIMF7. (h) BLIMF8. (I) BLIMF9. (j) BLIMF10.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(I)
(j)The WPE of each BLIMF is shown in Table2. Obviously, BLIMF1 and BLIMF2 are signal-dominant BLIMFs, and the rest of the BLIMFs are noise-dominant BLIMFs. Then, the BLIMF1 and BLIMF2 are reconstructed to get the denoised signal by VMD, as shown in Figure 4. It can be seen that after VMD denoising, the signal appears with partial distortion and oscillation.Table 2
WPE of BLIMFs.
VMD
WPE
BLIMF1
0.3870
BLIMF2
0.3884
BLIMF3
0.4629
BLIMF4
0.5748
BLIMF5
0.7058
BLIMF6
0.8608
BLIMF7
0.9678
BLIMF8
0.9934
BLIMF9
0.8672
BLIMF10
0.7787Figure 4
Comparison of the noise-free signal and the denoised signal by VMD.WTD is conducted for the noise-dominant BLIMFs. In this paper, the db3 wavelet is selected and the number of decomposition levels is 3. Signal-dominant BLIMFs and noise-dominant BLIMFs after WTD are combined for signal reconstruction. VMD-WTD denoising performance is as shown in Figure5. Obviously, the partial distortion and oscillation still exist, which demonstrates that the signal denoising performance after VMD-WTD denoising is still not ideal.Figure 5
Comparison of the noise-free signal and the denoised signal by VMD-WTD.SSA is implemented for the denoised signal by VMD-WTD, and Figure6 shows singular entropy increments ΔE of the signal. It is found that the information quantity of the signal reaches the saturated state at order 4, so the first 4 elementary matrices are used for reconstruction. After diagonal averaging, the denoised signal by VMD-WTD-SSA is obtained, as shown in Figure 7. It can be seen that the signal waveform is well recovered after denoising. Through the comparison, it can be seen that VMD-WTD-SSA denoising performance is superior to VMD and VMD-WTD denoising performances obviously.Figure 6
Singular entropy increments of the denoised signal by VMD-WTD.Figure 7
Comparison of the noise-free signal and the denoised signal by VMD-WTD-SSA.


## 4.3. Comparison with Other Methods

In this paper, the signal-to-noise ratio (SNR), root mean square error (RMSE), and maximum absolute error (MAE) are used as the evaluation indexes of denoising performance. SNR reflects the performance of the denoising method, and the larger the value, the better the denoising performance; RMSE reflects the similarity between the denoised signal and the noise-free signal, and the smaller the value, the better the denoising performance; MAE reflects the real error between the denoised signal and the noise-free signal, and the smaller the value, the better the denoising performance. The three evaluation indexes are defined as follows:(23)SNR=10log∑i=1ny2t∑i=1nyt−yt2,RMSE=1n∑i=1nyt−y′t2,MAE=1n∑i=1nyt−y′t,where y (t) is the original signal, y′t is the denoised signal, and n is the length of the original signal.Six typical denoising methods are selected to compare with VMD-WTD-SSA. They are wavelet soft threshold denoising (WSTD) [4], EMD hard threshold denoising (EMD-HT) [11], EMD partial reconstruction (EMD-PR) [12], DFA-VMD [26], VMD-WTD [27], and VMD-AE [30]. Denoising results of different methods for the noisy signals with SNRin values varying from −10 dB to 20 dB are shown in Figure 8. It can be seen that compared with the other methods, VMD-WTD-SSA has higher SNRout and lower RMSE and MAE. Therefore, the results show that the VMD-WTD-SSA proposed in this paper has the best denoising performances with different SNRin values. The presented method is better than the other methods on the whole, and its performance on noisy signals with low SNR is more prominent. Besides, as the SNRin values continuously increases, the differences among these methods are narrowed.Figure 8
Denoising results for simulation signals with different SNRin: (a) SNRout, (b) RMSE, and (c) MAE.
(a)
(b)
(c)


## 5. Application to Load Signals

### 5.1. Acquisition of Load Signals

The signal acquisition system for the load test of a corn combine harvester is mainly composed of sensors (strain gauges and strain rosettes), an HBM data collecting instrument, and a computer, as shown in Figure9.Figure 9
Signal acquisition system.According to existing experience, positions with large stress are determined as the measuring locations, of which the strain signals are measured by sensors. The layout of measuring locations of the frame is shown in Figure10. There are six measuring locations in total.Figure 10
Layout of measurement locations of the frame.The sampling frequency of strain signals is 500 Hz, which can fully meet the engineering test. In order to study the fatigue durability of the frame of a corn combine harvester, it is necessary to obtain the load data of the frame under multiple working conditions, such as cement road, field road, and gravel road. The load data under the field working condition is the most representative, which is used in the following analysis. Besides, working sections of a corn combine harvester under the field condition include straight-line harvesting sections and turning sections.


### 5.2. Denoising of Measured Signals

The method, proposed in this paper, is used to process the measured signal segment at location 5 in the straight-line harvesting section, which is shown in Figure11. The scaling exponent α0 of the noisy signal is 1.49, and the normalized MIE threshold is 0.015 by equation (9). It is determined that the mode number K is 6 by MIE, and VMD is carried out for the signal, as shown in Figure 12.Figure 11
Measured signal in the straight-line harvesting section.Figure 12
Decomposition results of the measured signal in the straight-line harvesting section. (a) BLIMF1. (b) BLIMF2. (c) BLIMF3. (d) BLIMF4. (e) BLIMF5. (f) BLIMF6.
(a)
(b)
(c)
(d)
(e)
(f)Then, the first four BLIMFs are determined as the signal-dominant BLIMFs by WPE. Therefore, the last two BLIMFs after WTD processing are combined with the first four BLIMFs to get the denoised signal by VMD-WTD. After that, SSA is implemented for the denoised signal by VMD-WTD, whose singular entropy incrementsΔE are shown in Figure 13. It is found that the information quantity of the signal reaches the saturated state at order 11, so the first 11 elementary matrices are used for reconstruction. After diagonal averaging, the denoised signal by VMD-WTD-SSA is obtained, which is shown in Figure 14(a).Figure 13
Singular entropy increments of the measured signal in the working section of linear harvesting.Figure 14
Comparison before and after denoising of the measured signal in the working section of linear harvesting: (a) measured signal, (b) spectra diagram of the measured signal, (c) denoised signal, and (d) spectra diagram of the denoised signal.
(a)
(b)
(c)
(d)The measured signals in different sections and their spectra before and after denoising are compared, as shown in Figures14 and 15. It can be seen that the denoised signals by VMD-WTD-SSA not only retain the low-frequency effective components of the original signal but also eliminate the interference of high-frequency noise, which fully shows the superiority of the proposed method. Therefore, the denoising method proposed in this paper has a significant denoising performance on time domain load signals, which lays the foundation for further analysis and processing.Figure 15
Comparison before and after denoising of the measured signal in the working section of turning: (a) measured signal, (b) spectra diagram of the measured signal, (c) denoised signal, and (d) spectra diagram of the denoised signal.
(a)
(b)
(c)
(d)


## 5.1. Acquisition of Load Signals

The signal acquisition system for the load test of a corn combine harvester is mainly composed of sensors (strain gauges and strain rosettes), an HBM data collecting instrument, and a computer, as shown in Figure9.Figure 9
Signal acquisition system.According to existing experience, positions with large stress are determined as the measuring locations, of which the strain signals are measured by sensors. The layout of measuring locations of the frame is shown in Figure10. There are six measuring locations in total.Figure 10
Layout of measurement locations of the frame.The sampling frequency of strain signals is 500 Hz, which can fully meet the engineering test. In order to study the fatigue durability of the frame of a corn combine harvester, it is necessary to obtain the load data of the frame under multiple working conditions, such as cement road, field road, and gravel road. The load data under the field working condition is the most representative, which is used in the following analysis. Besides, working sections of a corn combine harvester under the field condition include straight-line harvesting sections and turning sections.


## 5.2. Denoising of Measured Signals

The method, proposed in this paper, is used to process the measured signal segment at location 5 in the straight-line harvesting section, which is shown in Figure11. The scaling exponent α0 of the noisy signal is 1.49, and the normalized MIE threshold is 0.015 by equation (9). It is determined that the mode number K is 6 by MIE, and VMD is carried out for the signal, as shown in Figure 12.Figure 11
Measured signal in the straight-line harvesting section.Figure 12
Decomposition results of the measured signal in the straight-line harvesting section. (a) BLIMF1. (b) BLIMF2. (c) BLIMF3. (d) BLIMF4. (e) BLIMF5. (f) BLIMF6.
(a)
(b)
(c)
(d)
(e)
(f)Then, the first four BLIMFs are determined as the signal-dominant BLIMFs by WPE. Therefore, the last two BLIMFs after WTD processing are combined with the first four BLIMFs to get the denoised signal by VMD-WTD. After that, SSA is implemented for the denoised signal by VMD-WTD, whose singular entropy incrementsΔE are shown in Figure 13. It is found that the information quantity of the signal reaches the saturated state at order 11, so the first 11 elementary matrices are used for reconstruction. After diagonal averaging, the denoised signal by VMD-WTD-SSA is obtained, which is shown in Figure 14(a).Figure 13
Singular entropy increments of the measured signal in the working section of linear harvesting.Figure 14
Comparison before and after denoising of the measured signal in the working section of linear harvesting: (a) measured signal, (b) spectra diagram of the measured signal, (c) denoised signal, and (d) spectra diagram of the denoised signal.
(a)
(b)
(c)
(d)The measured signals in different sections and their spectra before and after denoising are compared, as shown in Figures14 and 15. It can be seen that the denoised signals by VMD-WTD-SSA not only retain the low-frequency effective components of the original signal but also eliminate the interference of high-frequency noise, which fully shows the superiority of the proposed method. Therefore, the denoising method proposed in this paper has a significant denoising performance on time domain load signals, which lays the foundation for further analysis and processing.Figure 15
Comparison before and after denoising of the measured signal in the working section of turning: (a) measured signal, (b) spectra diagram of the measured signal, (c) denoised signal, and (d) spectra diagram of the denoised signal.
(a)
(b)
(c)
(d)


## 6. Conclusion

In order to improve the denoising performance of time domain load signals, an efficient denoising method based on VMD, WTD, and SSA is proposed in this paper. Firstly, a simple criterion including MIE and the scaling exponentα0, which can be obtained by using DFA, is designed to select the proper mode number K for VMD. In addition, the BLIMFs are divided into the noise-dominate and signal-dominate parts by using WPE, which can help identify the useful information from the noisy signal. Furthermore, the presented method takes full advantages of VMD, WTD, and SSA, which can suppress noise and retain useful components of the original signal. The denoising performance of the proposed method is quantitatively evaluated by using the simulation signals, which are composed of different frequency components and Gaussian white noise. Experimental results of the presented simulation signals with SNRin values varying from −10 dB to 20 dB demonstrate that the proposed method outperforms WSTD, EMD-HT, EMD-PR, DFA-VMD, VMD-WTD, and VMD-AE. Moreover, applied to measured load signals, the presented method can not only effectively remove the noise component but also retain the effective component of the original signal, which proves its feasibility.Although the effectiveness of the proposed method has been verified by the denoising results of synthetic and measured signals, its universality is worth investigating further through a large number of simulation signals, which are comprised of various original signals and noise. Besides, the feasibility of the proposed method is qualitatively evaluated by the spectra diagrams of denoised measured load signals. Nevertheless, the influence of the denoising process on the fatigue durability analysis of mechanical structures needs to be further studied. Further work will be carried out to study the universality of the presented approach and apply the denoised measured load signals to the fatigue durability analysis of mechanical structures.



---

*Source: 1485937-2020-04-29.xml*