An optimized variational mode decomposition for extracting weak feature of viscoelastic sandwich cylindrical structures

A viscoelastic sandwich cylindrical structure (VSCS) is typically composed of a viscoelastic layer confined between two identical elastic and stiff layers [1–4]. Due to its viscoelastic and low-density properties, viscoelastic sandwich materials can play important roles in suppressing vibration, reducing noise and sealing, etc, almost without increasing the weight. Therefore, VSCSs are increasingly being applied in aviation, navigation and weapons, etc [5]. However, the mechanical performance of the sandwich material is susceptible to aging characteristics, temperature and load factors, etc, in service. The vibration of the VSCS will inevitably increase under external loads due to its performance degradation. What is more serious is that relative rotation or slippage between layers may occur [6], which will lead to the failure of the structure. Therefore, it is highly necessary to detect the features of interlayer slipping as early as possible to avoid accidents. Vibration signal analysis is the most commonly used approach of feature extraction due to its easy measurement, nondestructive nature and high correlation with structural dynamics [7–10]. However, the measured vibration signals from the VSCS are complex and nonstationary in nature, and the feature components are often submerged in background noise, such as excitation frequency components, meaning precise feature extraction in the VSCS faces great difficulties.

For the past few years, many methods have been developed to extract the weak characteristic components of nonstationary signals. For example, due to its multi-resolution property, wavelet transform can analyze nonstationary signals [11, 12]. However, its performance heavily relies on the choice of the mother wavelet function. The wavelet basis function does not always match the fault feature very well [13]. Huang et al proposed empirical mode decomposition (EMD), which adaptively decomposes the signals into a series of intrinsic mode functions (IMF), and the equivalent basis functions are entirely determined by the signal itself [14]. However, because EMD is sensitive to noise, mode mixing (where multiple components are decomposed into an IMF) frequently happens. Later, its extended versions, such as the ensemble EMD (EEMD) and the complementary EEMD (CEEMD), were proposed by Huang et al [15, 16]. EEMD can alleviate mode mixing to some degree by adding white noise to the original signal many times and averaging the decomposed results [17]. However, the reconstructed signal is easily contaminated by the residual noise [18]. The CEEMD, in which white noise is added in pairs to the analyzed signal with plus and minus signs performs better in eliminating the contamination of the residual noise. However, EMD-based methods are empirical algorithms and sensitive to noise.

To overcome the limitations of EMD-based methods, Dragomiretskiy and Zosso recently proposed the variational mode decomposition (VMD) algorithm, which can adaptively decompose a multi-component signal into a number of narrowband modes (called a variational mode function or variational mode component, VMF) in the framework of variational theory [19]. VMFs and their corresponding center frequencies are collectively derived by solving a constrained problem which minimizes the summation of the bandwidth (SBW) of all VMFs, while maintaining the data fidelity of the reconstructed signal. These VMFs can reproduce the original signal either exactly or in a least-squares sense. Due to its excellent decomposition capability, VMD outperforms EMD and the wavelet transform in extracting feature components [20, 21].

However, as a parameterized signal processing method, the decomposition performance of the VMD heavily relies on its parameters, among which the mode number K and the balancing parameter of the data-fidelity constraint $\alpha$ are two crucial parameters [22, 23]. The mode number K defines the number of VMFs, and the parameter $\alpha$ weighs the data fidelity of the reconstructed signal and the SBW of all the VMFs. Since the VMD cannot adaptively select the two parameters, they need to be artificially defined first. Nevertheless, inappropriate parameters not only increase the computational cost but also cause fuzzy results or information loss. In order to overcome the confusion of parameter selection, Wang et al discussed the filter bank property of the VMD algorithm through simulations, and suggested empirical values of parameters for extracting the oscillation components, removing the trend component and detecting the impact components, respectively [24]. Taking the accuracy of the center frequencies of VMFs into consideration, Upadhyay et al selected K to be 2 and a smaller $\alpha$ to detect the instantaneous components of the speech signals [25]. In order to extract the weak fault feature of the locomotive wheel bearing in the strong noise background [26] and the diagnosis turbine rotor system fault [27], both Li et al and Liu et al adopted the correlation coefficient to determine the mode number K so as to overcome the problem of information loss or over decomposition. Shi et al proposed an adaptive algorithm to choose K and $\alpha$ optimally through the setting threshold [28]. The above literature selected the parameters K and $\alpha$ mainly based on empirical values or the setting threshold, which is helpful for extracting powerful components instead of weak ones. Tang et al [29] and Yi et al [30] selected the parameters K and $\alpha$ based on particle swarm optimization (PSO) for extracting the characteristic components of the rolling bearing. However, it is difficult for the PSO to construct an appropriate fitness function. On the other hand, if the feature components are much weaker than the background noise, it is impossible for the VMD to separate them simultaneously [31]. However, weak feature components are always closely related to the status of VSCSs. Therefore, it is necessary to optimize the parameters and to explore new techniques to improve the capability of the VMD to extract weak feature components.

The aim of this paper is to propose a novel method based on the optimized VMD to precisely extract the weak feature of the VSCS. In this method, a new parameter optimization algorithm is proposed to reveal the influences of parameters K and $\alpha$ on the decomposition performance of the VMD, and to adaptively obtain the optimized parameters simultaneously. Then, the original signals of the VSCS are decomposed twice by the VMD with the optimized parameters, obtaining the strong and the weak components step by step. The rest of this paper is organized as follows: the VMD algorithm is briefly reviewed in section 2. In section 3, in order to obtain optimized results, an optimized VMD is proposed. First, a novel index is introduced to measure the performance of the VMD algorithm. Then the parameter optimization algorithm is proposed to select the optimized parameters simultaneously. Based on the optimized VMD, a weak feature extraction method is developed in section 4. In section 5, simulation and experiment are carried out to verify the effectiveness of the proposed method. Its robustness to Gaussian white noise and impulse noise is discussed. Finally conclusions and the next work are summarized in section 6.

Based on the idea that the SBW of all VMFs is much smaller than the bandwidth of a multi-component signal, the VMD can adaptively decompose the multi-component signal into a series of quasi-orthogonal VMFs by solving a constrained optimization problem in the framework of variational theory. The constrained optimization problem views the VMF ${{u}_{k}}$ and the corresponding center frequency ${{\omega }_{k}}$ as variables. VMD is a nonrecursive, intrinsic and adaptive signal decomposition method, which is robust to noise and sampling.

For a 1D signal $x$, given that it consists of VMFs ${{u}_{k}}(k=1,2,\ldots ,K)$, the process of the VMD algorithm is described as follows [19]:

• Step 1,  
  estimate the bandwidth of the kth VMF ${{u}_{k}}(k=1,2,\ldots ,K)$ as follows:

  • (1)  
    For each VMF ${{u}_{k}}$, compute its analytic signal $u_{k}^{a}$ by the Hilbert transform to obtain a unilateral frequency spectrum, i.e. $u_{k}^{a}=(\delta (t)+\frac{i}{\pi t})*u{\hspace{0pt}}_{k}(t)$.
  • (2)  
    With the estimated center frequency ${{\omega }_{k}}$, shift the frequency spectrum of the analytic signal $u_{k}^{a}$ to the baseband to get the demodulated signal $u_{k}^{d}$, i.e. $u_{k}^{d}=[(\delta (t)+\frac{i}{\pi t})*u{\hspace{0pt}}_{k}(t)]{{{\rm e}}^{-{\rm i}{{\omega }_{k}}t}}$.
  • (3)  
    Estimate the bandwidth of the demodulated signal $u_{k}^{d}$ through the ${{H}^{1}}$ Gaussian smoothness (the squared ${{L}^{2}}$-norm of the gradient), i.e. ${{\Delta }_{k}}=\left\Vert {{\partial }_{t}}[(\delta (t)+\frac{i}{\pi t})*u{\hspace{0pt}}_{k}(t)]{{{\rm e}}^{-j{{\omega }_{k}}t}} \right\Vert_{2}^{2}$.
• Step 2,  
  construct the constrained variational problem, where the objective function is the SBW of all VMFs, and the constrained function is that the original signal is reconstructed by all VMFs, that is

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}cccccccccccccccccccc@{}} & \underset{\{{{u}_{k}}\},\{{{\omega }_{k}}\}}{\mathop{{\rm min}}}\,\left\{\sum\limits_{k}{{{\Delta }_{k}}} \right\} \nonumber \\ & {\rm s}.{\rm t}.\quad \quad \sum\limits_{k}{u{\hspace{0pt}}_{k}(t)}=x \nonumber \\ \end{array}\nonumber \end{align} \tag{ 1 }$$

  where $\left\{u{\hspace{0pt}}_{k} \right\}:=\{u{\hspace{0pt}}_{1},\cdots ,u{\hspace{0pt}}_{K}\}$ and ${\{\omega {\hspace{0pt}}_{k} \}:=\{\omega {\hspace{0pt}}_{1},\cdots ,}$ $\omega {\hspace{0pt}}_{K}\}$ are shorthand notations for the set of all VMFs and their center frequencies, respectively. $\delta$ is the Dirac distribution, t is time and $\ast$ denotes convolution.
• Step 3,  
  convert the problem (1) into an unconstrained optimization problem by the method of multipliers. The augmented Lagrangian function is expressed as follows:

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle &L(\{{{u}_{k}}\},\{{{\omega }_{k}}\},\lambda):=\alpha \sum\limits_{k}{\left\Vert {{\partial }_{t}}[(\delta (t)+\frac{i}{\pi t})*u{\hspace{0pt}}_{k}(t)]{{{\rm e}}^{-{\rm i}{{\omega }_{k}}t}} \right\Vert}_{2}^{2}\nonumber \\ &+\left\langle \lambda (t),x(t)-\sum\limits_{k}{u{\hspace{0pt}}_{k}(t)} \right\rangle +\left\Vert x(t)-\sum\limits_{k}{u{\hspace{0pt}}_{k}(t)} \right\Vert_{2}^{2}\nonumber \end{align} \tag{ 2 }$$

  where $\alpha$ is the balancing parameter of the data-fidelity constraint and $\lambda (t)$ denotes the Lagrangian multiplier.
• Step 4,  
  solve the unconstrained problem (2) by the alternate direction method of multipliers (ADMM), which iterates the variables ${{u}_{k}}$ and ${{\omega }_{k}}$ alternately in the framework of variational theory. The solutions of problem (1) are given as follows:

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \hat{u}_{k}^{n+1}(\omega):=\frac{\hat{x}(\omega)-\sum\nolimits_{i\ne k}{{{{\hat{u}}}_{i}}(\omega)+\frac{\hat{\lambda }(\omega)}{2}}}{1+2\alpha {{(\omega -\omega _{k}^{n})}^{2}}}\nonumber \end{align} \tag{ 3 }$$

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \omega _{k}^{n+1}:=\frac{\int_{0}^{\infty }{\omega {{\left| \hat{u}_{k}^{{}}(\omega) \right|}^{2}}{\rm d}\omega }}{\int_{0}^{\infty }{{{\left| \hat{u}(\omega) \right|}^{2}}{\rm d}\omega }}\nonumber \end{align} \tag{ 4 }$$

  where ${{\hat{u}}_{k}}$,$\hat{x}$ and $\hat{\lambda }$ denote the Fourier transforms of the VMF ${{u}_{k}}$, the original signal $x$ and the Lagrangian multiplier $\lambda$, respectively. ${{\omega }_{k}}$ is the center of gravity of the corresponding VMF power spectrum, and the superscript n is the number of iterations.

Finally, taking the real part of the inverse Fourier transform of equation (3), the time domain form of the VMF is obtained. More details can be found in [19].

From the VMD algorithm depicted above, it can be seen that it is a parameterized signal decomposition algorithm. The parameters, such as K, $\alpha$ and $\lambda$, must be initialized in advance. Using the method of multipliers to solve the constraint optimization problem (1), the Lagrangian multiplier term (i.e. the middle term on the right-hand side of equation (2)) is to enforce the constraint. When the Lagrangian multiplier $\lambda$ is positive, VMFs can accurately reconstruct the input signal. In this case, the VMD is good at dealing with the noise-free signal or low noise signal. On the other hand, if $\lambda$ is zero, the Lagrangian function is reduced to the penalty function [32, 33], and the input signal is reconstructed by the VMFs in the least-squares sense. In this situation, the VMD can deal with the noisy signal well. Since the practical signals are more or less subjected to noise, $\lambda$ is assumed to be zero in the following. For parameters K and $\alpha$, detailed studies are given in section 3.

Since the parameters play important roles on its decomposition performance, the VMD with the optimized parameters instead of the empirical ones should be used to decompose the signals and effectively extract the feature components. To this end, a novel index is introduced to measure the decomposition performance of the VMD, and a parameter optimization algorithm is proposed to optimize the crucial parameters K and $\alpha$ to obtain better decomposition results.

3.1. Measurement index

From the constraint optimization problem (1), the SBW is closely related to the decomposition results of the VMD algorithm. So the SBW of all the VMFs is introduced to measure the influence of parameters K and $\alpha$ on the decomposition performance of the VMD.

For an analytic signal s(t), its bandwidth is equivalent to the width of the rectangle whose area is equal to that of the power spectrum of s(t), and whose height is the amplitude at the centroid of the power spectrum [19, 34]. The bandwidth of s(t) can be calculated as follows:

The frequency domain energy (E) of s(t)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E=\int_{0}^{\infty }{{{\left| \hat{s}(\,f) \right|}^{2}}}{\rm d}f.\nonumber \end{align} \tag{ 5 }$$

The center frequency (CF) of s(t)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle CF=\frac{1}{E}\int_{0}^{\infty }{f{{\left| \hat{s}(\,f) \right|}^{2}}{\rm d}f}.\nonumber \end{align} \tag{ 6 }$$

The bandwidth (BW) of s(t)

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta _{f}^{2}=\int_{0}^{\infty }{{{(\,f-CF)}^{2}}\frac{{{\left| \hat{s}(\,f) \right|}^{2}}}{E}}{\rm d}f\nonumber \end{align} \tag{ 7 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle BW=2{{\Delta }_{f}}\nonumber \end{align} \tag{ 8 }$$

where $\hat{s}(\,f)$ is the amplitude spectrum of s(t) and $f$ denotes the frequency.

The bandwidth can measure the degree of signal energy concentration with respect to the center frequency. The smaller the BW, the more energy is concentrated around the center frequency of the signal. The BW from equation (8) is essentially the same as the bandwidth in the VMD according to the ${{H}^{1}}$ Gaussian smoothness of the demodulated signal. However, for the convenience of calculation, the bandwidth from equation (8) is adopted in the following:

Define the SBW as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm SBW}=\sum\limits_{k=1}^{K}{{\rm B}{{{\rm W}}_{k}}}\nonumber \end{align} \tag{ 9 }$$

where ${\rm B}{{{\rm W}}_{k}}$ denotes the bandwidth of the kth VMF ${{u}_{k}}$.

For a multicomponent signal composed of narrowband modes without a frequency spectrum overlap, its bandwidth is much greater than the SBW of all the narrowband modes. In theory, the higher the SNR of a narrowband mode, the smaller the bandwidth.

3.2. Parameter optimization algorithm

The VMD algorithm decomposes the original signal into a series of VMFs based on pre-defined parameters weighing the fidelity of the reconstructed signal and the SBW. Since the parameters K and $\alpha$ affect the decomposition results, they will also influence the SBW. If the parameter K is greater than the actual number of components in the signal, it will inevitably lead to spurious modes (consisting of noise content) or mode splitting (the same component is shared by several VMFs). In contrast, a too small K easily results in mode mixing or mode loss. In these two cases of K, the SBW will get larger. For the parameter $\alpha$, a too small value easily leads to VMFs still containing a large amount of noise, and even the occurrence of mode mixing. Conversely, a too large $\alpha$ will cause mode loss accompanied by spurious modes or mode splitting. In these two situations of $\alpha$, the SBW will also become larger. In short, the inappropriate K or $\alpha$ will result in the decomposition performance of VMD degradation and a larger SBW, and vice versa.

Based on the idea that the SBW is smaller only if optimal K and $\alpha$ are selected, an adaptive parameter optimization algorithm is proposed, based on the index SBW, in order to obtain a better decomposition effect. The flowchart of the proposed algorithm is shown in figure 1, which consists of two nested loops. The inner loop aims to find the local minimum of the SBW with respect to $\alpha$ assuming the K constant, denoted by SBW_r. The outer loop is to find the local minimum of the SBW_r with respect to K, i.e. ${\rm SB}{{{\rm W}}_{{\rm op}}}$. The optimized parameters K and $\alpha$ are the parameters corresponding to ${\rm SB}{{{\rm W}}_{{\rm op}}}$. The detailed process is as follows:

• Step 1,  
  initialize parameters K and $\alpha$, ${{K}_{i}}={{K}_{0}}:1:{{K}_{t}}$, ${{\alpha }_{i}}={{\alpha }_{0}}:\Delta \alpha :{{\alpha }_{t}}$, where $\Delta \alpha$ is the step size, and let $j=1$, $m=1$.
• Step 2,  
  let $K={{K}_{i}}(m)$.
• Step 3,  
  let $\alpha ={{\alpha }_{i}}(\,j)$.
• Step 4,  
  execute the VMD algorithm to get VMFs, i.e. ${{u}_{k}}(k=1,2,\ldots ,K)$.
• Step 5,  
  calculate the SBW by equation (9), denoted by $s(\,j)={\rm SBW}$.
• Step 6,  
  determine whether the SBW is a local minimum with regard to $\alpha$ (except for the endpoints), and if so, go to step 8.
• Step 7,  
  determine whether $\alpha$ is smaller than ${{\alpha }_{t}}$, and if so, let $\alpha =\alpha +\Delta \alpha$, then go to step 4.
• Step 8,  
  record K and $\alpha$, and the local minimum of SBW with respect to $\alpha$, denoted by SBW_r, i.e. ${{K}_{r}}(m)=K,{{\alpha }_{r}}(m)={{\alpha }_{i}}(\,j-1),$$SB{{W}_{r}}(m)=s(\,j-1)$.
• Step 9,  
  determine whether $SB{{W}_{r}}$ is a local minimum with respect to K (except for the endpoints), and if so, go to step 11.
• Step 10,  
  determine whether K is smaller than ${{K}_{t}}$; if so, let $m=m+1$, and go to step 3.
• Step 11,  
  record the minimum of ${\rm SB}{{{\rm W}}_{r}}$ with respect to K, denoted by ${\rm SB}{{{\rm W}}_{{\rm op}}}$, the corresponding parameters K and $\alpha$ are the optimized parameters.

Zoom In Zoom Out Reset image size

The proposed algorithm is not only able to adaptively obtain the optimized parameters K and $\alpha$, but also reveal the change of decomposition performance of VMD with parameters K and $\alpha$, which can overcome the deficiency of the empirical value or thresholds.

To ease understanding of the proposed algorithm, a synthesis signal shown in equation (10) is designed, which is composed of one fundamental frequency component, three harmonic components and white noise. The signal is shown in figure 2.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle y(t)=&\,3{\rm cos}(60\pi t)+4{\rm cos}(120\pi t+1)\nonumber \\ &+2{\rm cos}(180\pi t+1.5)+{\rm cos}(240\pi t-0.5)+n(t).\nonumber \end{align} \tag{ 10 }$$

Zoom In Zoom Out Reset image size

Define $K=2:1:6$ and $\alpha =300:200:100\,000$. The results of the parameter optimization algorithm are shown in figure 3, where (a)–(d) are the changing of SBW with $\alpha$ when K is 2, 3, 4 and 5, respectively. The results imply that SBW is influenced by parameters K and $\alpha$. On the one hand, for each fixed K, there is a local minimum of SBW with respect to $\alpha$, denoted by ${\rm SB}{{{\rm W}}_{r}}$, such as A (5100, 0.013 44), B (5100, 0.018 16), C (59 500, 0.010 42) and D (900, 0.031 49) shown in figure 3. This also indicates that the proposed algorithm was able to find the optimized $\alpha$ directly once K was known. On the other hand, for different K, the corresponding ${\rm SB}{{{\rm W}}_{r}}$ is also different. Once the local minimum of ${\rm SB}{{{\rm W}}_{r}}$ with regard to K, namely ${\rm SB}{{{\rm W}}_{{\rm op}}}$ (0.010 42), is achieved, the proposed algorithm terminates. The corresponding parameters K and $\alpha$, i.e. (4, 59 500), are viewed as the resulting optimized parameters. The corresponding decomposition results are shown in figure 4. VMF1–VMF4 are the estimates of 1×  (30 Hz), 4×  (120 Hz), 3×  (90 Hz) and 2×  (60 Hz) components in the original signal, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Besides obtaining the optimized parameters, the proposed algorithm can reveal the characteristics of the influence of parameters on the decomposition performance of the VMD. On the one hand, it is obvious that the SBW does not monotonically change with $\alpha$ for each constant K, as shown in figure 3. When $\alpha$ is smaller, the SBW decreases with the increase of $\alpha$. However, a too large $\alpha$ will result in a larger SBW. On the other hand, the change of SBW in its nature reflects the different decomposition performance of the VMD. A larger SBW indicates the poor decomposition performance of the VMD. For example, with $\alpha =100,K=4$ the decomposition results lead to a larger SBW, owing to mode mixing in VMF2 and VMF4 being a spurious mode, as shown in figure 5. With $\alpha =60\,000,K=4$ the decomposition results shown in figure 6 result in a larger SBW, owing to the 2×  component losing and VMF4 being a spurious mode. The VMD may decompose the original signal rather well between the interval where the SBW is relatively small, for example the interval ${{C}_{0}}-C$ in figure 3(c). However, at the beginning of the interval, the noise content contained in the decomposition results is relatively large. Therefore, the parameters corresponding to the local minimum of the SBW are looked on as the optimized parameters.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The above analysis can also be illustrated by the convergence of the center frequency of each VMF versus $\alpha$ shown in figure 7(a), from which it can be seen that except for the first and last cases, the center frequency of each VMF approximates a real value, respectively. The observations agree with the results from figure 3(c). In addition, the number of iterations versus $\alpha$ for K  =  4 is shown in figure 7(b). It can be seen that the convergence speed of the optimization algorithm fluctuates with $\alpha$. Although the convergence rate of the optimization algorithm with the optimized parameters is not the fastest, the number of iterations is not large at all.

Zoom In Zoom Out Reset image size

There are a few points to be noted when the parameter optimization algorithm is used. First the initial ranges and step sizes of the parameters $\alpha$ and K. As for $\alpha$, the lower limit is empirically larger than ten and its upper limit can be defined arbitrarily to be a larger value. Its step size is not strict. A larger step saves computational time and a smaller step size helps to improve its resulting accuracy. For K, its range can be defined according to the spectrum of the signal. Its lower limit should be larger than one generally; its step size must be one. From the results shown in figure 3, it can be seen that as long as each initial range covers the optimal value, respectively, the optimal parameters will be adaptively chosen and the initial ranges will not influence the optimized results. Second, the parameter $\alpha$ must change from small to large in order to maintain data fidelity while obtaining the local minimum of SBW. Otherwise its initial range is difficult to define, and inappropriate parameters will be selected. In order to illustrate this point, we revise the parameter optimization algorithm not to terminate until $\alpha$ reaches its last value. Define K  =  3, 4, $\alpha =300:200:100\,000$, and the change of SBW with $\alpha$ is given in figure 8. From the figure it is obvious that there is more than one local minimal SBW due to the information loss caused by a larger $\alpha$. Once $\alpha$ changes from large to small, the first local minimal SBW will vary with the upper limit of $\alpha$, and no optimal parameters or optimal decomposition will result.

Zoom In Zoom Out Reset image size

Although the optimized VMD has excellent capability to decompose a multi-component signal, it is intractable for it to extract the weak feature components from strong background noise directly. However, changes of these weak feature components often indicate changes of the operating condition of the structures. Therefore, it is significant for failure detection to be able to accurately extract them from the vibration response signals. In this section, a novel weak feature component extraction method is proposed based on the optimized VMD.

A flowchart of the proposed method is shown in figure 9. According to the research in section 3, the decomposition results of the VMD are strongly influenced by parameters K and $\alpha$. If a smaller K and a larger $\alpha$ are pre-defined, the VMD seriously suffers from mode loss. This so-called shortcoming is made use of in this paper. First, through pre-defining $K=1$ and a larger $\alpha$ in the VMD, the low-frequency strong component is extracted from the original signal $y(t)$, denoted by VMF1. Then the difference between the original signal and VMF1 is defined as the residual signal $r(t)$, which is $r(t)=y(t)-{\rm VMF}1$. Finally, weak feature components are extracted from the residual signal by VMD with the optimized parameters.

Zoom In Zoom Out Reset image size

5.1. Simulation validation

When studying the fault characteristics of the interlayer slipping of the VSCS or validating the dynamical model, the VSCS is often excited under sinusoidal signals. The measured vibrational signal is composed of a powerful exciting frequency component accompanied by weak components. In order to verify the effectiveness of the proposed method, a synthesized signal is designed, which consists of a strong fundamental frequency component, three weak harmonic components and Gauss white noise, as shown in equation (11). The SNR is 35db, the sampling frequency is 2000 Hz and the data length is 10 000. The signal is shown in figure 10.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle y(t)=&\,30{\rm cos}(60\pi t)+0.5{\rm cos}(120\pi t)+0.3{\rm cos}(180\pi t)\nonumber \\ &+0.2{\rm cos}(240\pi t)+n(t).\nonumber \end{align} \tag{ 11 }$$

Zoom In Zoom Out Reset image size

Define $K=1,\alpha ={{10}^{5}}$; the low-frequency strong fundamental frequency component (30 Hz) of the original signal is extracted first. The result is shown in figure 11. Compared with the original signal, it can be seen that the VMF1 is mainly composed of the 1×  (30 Hz) component.

Zoom In Zoom Out Reset image size

The residual signal $r(t)$ mainly consists of three weak components, as shown in figure 12.

Zoom In Zoom Out Reset image size

In order to obtain the optimized decomposition of the three weak components, the residual signal is analyzed by the optimized VMD. Define $K=2:1:4$,$\alpha =100:50:100\,000$. The results of the optimized parameters are shown in table 1. The ${\rm SB}{{{\rm W}}_{{\rm op}}}$ is $0.024\,95$ and the optimized parameters are $K=3,\alpha =5200$. The corresponding decomposition results are shown in figure 13, VMF2–VMF4 are the reconstructed components of the 2×  (60 Hz), 4×  (120 Hz) and 3×  (90 Hz) harmonic components, respectively. It can be seen that the three weak components are precisely separated without mode mixing, etc.

Table 1. The results of the optimized parameters.

K	$\alpha$	${\rm SB}{{{\rm W}}_{r}}$
2	1150	0.07212
3	5200	0.02495
4	4800	0.05044

Zoom In Zoom Out Reset image size

5.2. Robustness of the proposed method

The above results show that the proposed method can effectively extract the weak feature components mingled in low-frequency strong ones by taking full advantage of the characteristics of the effect of the parameters on the decomposition performance of the VMD. However, the actual measured VSCS vibration signals are often contaminated by wide-band background noise. Thus, the fault feature information is usually submerged in the noise environment. In order to prove the robustness of the method to noise, we discuss the feature extraction effect of the proposed method under different noise levels. First, Gauss white noise is added to the original signal making SNR 30 dB, 33 dB, 35 dB, 40 dB, 45 dB and 50 dB in turn. Then, the noisy signal is processed by the proposed method. Finally, the cross correlation coefficient between the VMFs and the original components under the different noise levels is introduced to evaluate their similarity, as defined by equation (12) [27]. The higher the cross correlation coefficient, the more similar the VMF and the original component, and vice versa. If the cross correlation coefficient is one, the VMF is thought of as a reproduction of the original component. From the results shown in figure 14, it can be seen that the higher the SNR, the better the decomposition effect. The obtained strong fundamental component is little affected by the added noise, and the noise mainly influences the extraction effect of the weaker components. After the above-mentioned noises are added to the original signal, the SNRs of the obtained residual signal are about  −3.75 dB, −0.75 dB, 1.25 dB, 6.25 dB, 11.25 dB and 16.25 dB, respectively. Due to the spectral overlap between the weak components and the Gauss white noise, it is hard for the proposed method to extract weak components accurately under strong Gauss white noise conditions. However, when the SNR is no smaller than 35 dB, the cross correlation coefficients of all VMFs are larger than 0.9, and the extraction effect of the weak components is satisfactory.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle C={\sum\limits_{n=1}^{N}{({{u}_{k}}(n)-{{{\bar{u}}}_{k}})({{y}_{k}}(n)-{{{\bar{y}}}_{k}})}}/{\sqrt{\sum\limits_{n=1}^{N}{{{({{u}_{k}}(n)-{{{\bar{u}}}_{k}})}^{2}}\sum\limits_{n=1}^{N}{{{({{y}_{k}}(n)-{{{\bar{y}}}_{k}})}^{2}}}}}}\;\nonumber \end{align} \tag{ 12 }$$

where ${{u}_{k}}$,${{y}_{k}}$ denotes the kth VMF and the original component, respectively, and N is the data length.

Zoom In Zoom Out Reset image size

Besides white noise contamination, the measured signals are often disturbed by impulse noises in engineering. To discuss the robustness of the proposed method further, the impulse noise is added to the original noise free signal making SNR 25 dB, as shown in figures 15(a)–(d). Then the noisy signal is decomposed by the proposed method. The strong 1 ×  component (30 Hz) is separated perfectly, as shown in figures 15(e) and (f), and the obtained residual signal is effected severely by the impulse noise, as shown in figures 15(g) and (h). The impulse noise can be viewed as the sum of multiple sinusoidal components. There are no spectral overlaps between the weak components and the sinusoidal components, as shown in figure 15(h), therefore the weak components can be extracted from the residual signal precisely by the proposed method. The results are shown in figure 16. The cross relation coefficients of the 1×, 2×, 3×  and 4×  harmonic component and the corresponding reconstructed signals are 1, 0.9963, 0.9969 and 0.9991, respectively. The above analysis results indicate the proposed method is robust to impulse noise.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

5.3. Experiment evacuation of the proposed method

In this section, the effectiveness of the proposed method is verified by the VSCS experiment, and the results are compared with those by the VMD with nonoptimized parameters and the optimized CEEMD method.

5.3.1. Experiment setup and data acquisition.

The experimental system consists of two parts: the vibration controlling subsystem and the data acquisition subsystem. The former controls the vibration of the vibration table and is composed of sensors, the vibration table controller and controller system, etc, as shown in figure 17(a). The latter acquires signals from the VSCS experimental model, which is connected to the vibration table with bolts. The data acquisition subsystem is composed of acceleration sensors and LMS data acquisition systems, etc. The classical VSCS is a three-layer structure, which consists of a constrained layer, a sandwich layer and a base layer, as shown in figure 17(b). The constrained layer and the base layer are metal layers and the sandwich layer is viscoelastic silicone foam. The radial gap between the constrained layer and the base layer is 3 mm, the thickness of the sandwich layer material is 3.5 mm, and the sandwich layer material is compressed between the constrained layer and the base layer by pre-stressed force. In order to investigate the interlayer slipping feature, sinusoidal vibrations of different intensities and frequencies are excited in the VSCS model in the vertical direction. At each given excitation frequency, the intensity is increased step-by-step until interlayer slipping occurs. At each vibrational intensity, the duration of excitation is about 5 min, and the vibration response signals are measured by acceleration sensors 1–5 mounted on the constrained layer, as shown in figure 17(b).

Zoom In Zoom Out Reset image size

The vibration signal from sensor 5 is shown in figure 18, where the excitation frequency is set at 43.75 Hz, the sampling frequency is 2048 Hz and the data length is 8192. It can be seen from the figure that the vibration signal is mainly composed of the strong fundamental frequency (43.75 Hz) component, accompanied by a large number of weak harmonic components.

Zoom In Zoom Out Reset image size

5.3.2. Results and analysis.

In order to extract the weak feature components contained in the vibration signal shown in figure 18, the signal is decomposed by the proposed method in section 4.

First, the low-frequency strong component VMF1 (43.75 Hz) in the original vibration signal is extracted by the VMD with $K=1,\alpha =300\,000$, as shown in figures 19(a) and (b). As expected, the low-frequency strong component is well obtained, and the residual signal $r(t)$, as shown in figures 19(c) and (d), consists of weak harmonic components.

Zoom In Zoom Out Reset image size

Then, the residual signal is processed by the optimized VMD. The parameters are initialized as $K=6:1:10$, $\alpha =400:100:10\,000$. The optimized results are shown in table 2. It can be seen that there is the local minimum of SBW with respect to parameters K and $\alpha$, i.e. ${\rm SB}{{{\rm W}}_{{\rm op}}}=0.07868$. The resulting optimized parameters are ${{K}_{{\rm op}}}=8,{{\alpha }_{{\rm op}}}=3700$.

Table 2. The results of the optimized parameters.

K	$\alpha$	${\rm SB}{{{\rm W}}_{r}}$
6	1200	0.084 05
7	700	0.1072
8	3700	0.078 68
9	900	0.1755

Using the optimized parameters, the decomposition results are shown in figure 20. As can be seen from the figure, the multiple weak components in the residual signal are well separated, and the VMF2–VMF9 are mainly estimates of the 2×, 3×, 5×, 8×, 6×, 13×, 14×  and 10×  harmonic components, respectively. Obviously, the multiple weak feature components are well extracted from the strong ones by the proposed method.

Zoom In Zoom Out Reset image size

5.3.3. Comparisons.

In order to highlight the superiority of the proposed method, the original signal shown in figure 18 is analyzed by using the original VMD and the optimized CEEMD, respectively.

5.3.3.1. Comparison with the original VMD.

The original signal is decomposed by the original VMD with nonoptimized parameters, such as $\alpha =1000,K=9$. The results are shown in figure 21. It can be seen that VMF1, VMF2 and VMF7 share the 1×  component, and VMF3 and VMF8 share the 2×  harmonic component, which implies that mode splitting happens. VMF4 and VMF5 are the 3×  and 5×  harmonic components, respectively. VMF6 suffers from mode mixing, and VMF9 is a low-frequency component. Compared with the results of the proposed method, the results of the original VMD not only suffer from mode splitting but are also unable to capture enough weak feature components.

Zoom In Zoom Out Reset image size

5.3.3.2. Comparison with the optimized CEEMD.

The original signal is again decomposed by CEEMD, which is a standard substitute for EEMD [16, 18]. In CEEMD, there are two crucial parameters to be prescribed—that is the amplitude of the added white noise and the ensemble number—which will influence the decomposition [16]. In order to obtain better decomposition results, an index termed the relative root mean square error (RRMSE) was introduced to evaluate the decomposition performances under different noise levels to find the optimal level [18, 35]. The RRMSE expressed as equation (13) is defined as the ratio between the root mean square (RMS) of the error and the RMS of the original signal $x$, where the error is the difference between the original signal $x$ and the IMF ${{c}_{{\rm max}}}$ having the highest correlation with the vibration signal. More details can be found in [18, 35]. The range of noise levels is set [2 0.01], and the decreasing step size is 0.02. The optimal results are shown in figure 22, and the optimal level is 0.2. As far as the ensemble number is concerned, it is set to be 800 (400 pairs) since a larger ensemble number is helpful for reducing the remaining noise in each IMF [16, 18, 35].

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm RRMSE}=\sqrt{{\sum\limits_{k=1}^{N}{{{(x(k)-{{c}_{{\rm max}}}(k))}^{2}}}}/{\sum\limits_{k=1}^{N}{{{x}^{2}}(k)}}\;}\nonumber \end{align} \tag{ 13 }$$

where N is the number of samples in the original signal.

Zoom In Zoom Out Reset image size

A total of 15 IMFs were obtained by CEEMD with optimized parameters. To save space, only the first eight dominant IMFs are displayed, as shown in figure 23. It can be seen from the figure that IMF4 and IMF5 share the 1x component, which indicates that mode splitting occurs. IMF1–IMF3 contain multiple principal components, respectively, which implies that they all suffer from mode mixing. IMF6–IMF8 are meaningless low-frequency components. The results indicate that the optimized CEEMD fails to extract the weak feature components.

Zoom In Zoom Out Reset image size

From the above analysis, it can be concluded that the proposed method is effective. Compared with the optimized CEEMD and the original VMD, the proposed method can precisely extract weak feature components submerged in strong low- frequency ones without mode mixing, etc. The proposed method could be used to analyze other vibration signals under different loads so as to seek the feature components sensitive to interlayer slipping.

In this paper, a new weak feature extraction method based on a parameter optimization algorithm is proposed to solve the problem of crucial parameter selection and to improve the performance of the weak feature extraction of VMD. The results of the simulations and experiments show that

• (1)  
  The parameter optimization algorithm is not only able to adaptively select the optimized parameters K and $\alpha$ to obtain optimal decomposition results, but also reveal the changing characteristic of the decomposition performance of the VMD algorithm with the parameters K and $\alpha$. In addition, when K is known, the parameter optimization algorithm can also be used to select parameter $\alpha$ directly.
• (2)  
  The proposed method can precisely extract the weak feature components submerged in the low-frequency strong ones. Compared with the optimized CEEMD and the original VMD, the proposed method can avoid mode mixing, etc, and outperforms the optimized CEEMD and the original VMD in the weak feature component extraction.

The next work will study the influence of the initialization of the center frequencies on the decomposition performance of the VMD algorithm, and explore the feature information of the interlayer slipping of the VSCS.

The authors would like to express appreciation to the reviewers and editors for their valuable comments to improve the paper. The work is supported by Science Challenge Project and National Natural Science Foundation of China (No. 51775410).