Study on Denoising Method of Vibration Signal Induced by Tunnel Portal Blasting Based on WOA-VMD Algorithm

Abstract

Because of the impact of the complex environment of tunnel portals, the measured blasting vibration signals in a tunnel portal contains a lot of high-frequency noise. To achieve effective noise reduction, a novel method of noise reduction for blasting vibration signals based on the whale optimization algorithm (WOA) optimized with the variational mode decomposition (VMD) algorithm was proposed. The WOA algorithm is used to optimize globally for the mode number K and penalty factor α of VMD for measured signals and to determine the optimal parameters of [K, α], and to obtain the intrinsic mode function (IMF). Multi-scale permutation entropy (MPE) was used to identify and remove noise components in IMF, and then the reserved IMF was reconstructed to achieve a denoised signal. The method is applied to the blasting vibration analysis of the Xiali tunnel of the Jixin expressway in the Henan Province. Results indicate that the novel method can acquire the optimal decomposition mode number and identify the high frequency noise. Its denoising effect is better than the conventional VMD algorithm and the complete ensemble empirical mode decomposition with adaptive ability (CEEMDAN), which verifies the self-adaptivity and effectiveness of the WOA-VMD denoising method.

1. Introduction

Under long-term weathering, the surrounding rock in a tunnel portal section is relatively broken and has poor stability and blasting has adverse effects on the surrounding rock stability. A detailed analysis of the harmful effects of blasting vibration is crucial to avoid the disaster of engineering instability [1,2,3]. Due to the complex environmental influence of the tunnel–slope composition, the signal obtained by the field blasting vibration test contains a lot of noise, which changes the real signal amplitude and spectrum characteristics. The measured blasting vibration signal belongs to the non-stable time domain signal. It is the most effective denoising method to obtain the mode component signal of different frequency domains and then judge the noise component [4,5].

The algorithms commonly used for the decomposition of blasting vibration signals include the empirical mode decomposition (EMD) [6], local mean decomposition (LMD) [7], and variational mode decomposition (VMD) [8]. EMD and its improved algorithm can produce more ideal decomposition signals, but there are endpoint effect and modal stacking problems. Shao et al. [9] put forward a signal noise reduction model of complete ensemble empirical mode decomposition (CEEMD) in which low-pass filtering is constructed to effectively remove the mixed noise of the tunnel blasting signal. Peng et al. [10] proposed a smooth denoising model of complete ensemble empirical mode decomposition with adaptive ability (CEEMDAN) to reduce the noise in underwater drilling and blasting vibration signal, and achieve the great denoising effect. The LMD algorithm can adaptively decompose the signal into the product functions of several frequency modulations and envelope signals, which has good results for noise reduction and features the extraction of vibration signals, which is better than the EMD algorithm [7,11]. The methods of LMD, EMD, and their improved algorithms do not need to select base equations and have self-adaptability, but they will lead to mode mixing when dealing with noise signals.

Variational mode decomposition algorithm is a non-recursive signal decomposition with extensive adaptability and good robustness, which can solve endpoint distortion and mode confusion [12]. Peng et al. [13] and Fu et al. [14] studied the noise reduction and feature extraction of the blasting vibration signal of the VMD algorithm and verified its effectiveness. However, the mode number K of this algorithm is difficult to determine, which may lead to large manual intervention errors. The whale optimization algorithm (WOA) is a new type of intelligent optimization algorithm, and the optimization process has the advantages of few parameters, high accuracy, and fast convergence. In this paper, WOA is used to optimize VMD parameters, and is applied with practical engineering methods to verify the noise reduction effect.

2. WOA and VMD Algorithms

2.1. VMD Algorithm

Variational mode decomposition algorithm is an adaptive decomposition algorithm based on one-dimensional Hilbert transform and frequency mixing [15]. Its essence is to solve the variational problem, by iteratively calculating the center frequency and bandwidth of each mode function, which completes the signal frequency domain decomposition. The constrained variational equation is constructed as follows.

$$\begin{array}{c}\underset{\left\{u_k\right\},\left\{{\omega }_k\right\}}{\min }\left\{{\displaystyle \sum _k{\Vert {\partial }_t\left[\left(\delta (t)+\frac{j}{\pi t}\right)*u_k(t)\right]e^{-j{\omega }_{k}t}\Vert }_2^2}\right\}\\ s.t.{\displaystyle \sum _{k=1}^K{u}_k}=x(t)\end{array}$$

(1)

where x(t) is original signal; {u_k} is the intrinsic mode function (IMF) component; {ω_k} is the center frequency; δ(t) is the impulse function.

In order to solve Equation (1), the penalty coefficient α and Lagrange operator should be introduced to relieve the constraint of equation. Solve the extending Lagrange formula by alternating parameters of $u_k^{n+1}$, ${\omega }_k^{n+1}$, and ${\lambda }^{n+1}$

$$\begin{array}{c}{\widehat{u}}_k^{n+1}(\omega )=\frac{\widehat{f}(\omega )-{\displaystyle \sum _{i\ne k}{\widehat{u}}_i(\omega )}+\frac{\lambda (\omega )}{2}}{1+2\alpha (\omega -{\omega }_k)}\\ {\omega }^{n+1}=\frac{{\displaystyle {\int }_0^{\infty }\omega {\left|{\widehat{u}}_k(\omega )\right|}^{2}d\omega }}{{\displaystyle {\int }_0^{\infty }{\left|{\widehat{u}}_k(\omega )\right|}^{2}d\omega }}\end{array}$$

(2)

The original signal x(t) and the reconstructed signal $\widehat{x}(t)$ can be expressed as follow.

$$x(t)={\displaystyle \sum _{k=1}^K{u}_k}+u_r={\displaystyle \sum _{k=1}^{K}IM{F}_k}+u_r$$

(3)

$$\widehat{x}(t)={\displaystyle \sum _{k=1}^{K}IM{F}_k}$$

(4)

2.2. WOA Algorithm

Whale optimization algorithm simulates hunting behavior of humpback whales, and is applied to solve complex optimization problems [16]. The calculation process mainly includes three parts: surrounding prey, attacking prey, and searching prey.

(1) Surrounding prey

The hunting process is able to determine the optimal prey position and surround it. The mathematical model of this process is shown as follow.

$$\overrightarrow{D}=\left|\overrightarrow{C}{\overrightarrow{X}}^{*}-\overrightarrow{X}\left(t\right)\right|$$

(5)

$$\overrightarrow{X}\left(t+1\right)={\overrightarrow{X}}^{*}\left(t\right)-\overrightarrow{A}\overrightarrow{D}$$

(6)

where t is the number of iterations; ${\overrightarrow{X}}^{*}$ is the vector of prey position; $\overrightarrow{X}$ is the whale position vector; $\overrightarrow{D}$ is the hunting distance; $\overrightarrow{A}$ and $\overrightarrow{C}$ are coefficient vectors which are expressed as:

$$\overrightarrow{A}=2ag\overrightarrow{r}-a$$

(7)

$$\overrightarrow{C}=2g\overrightarrow{r}$$

(8)

where a is a convergence factor, and its value is [0, 2]; $\overrightarrow{r}$ is the random vector of [0, 1].

(2) Attacking prey

Two methods are used to simulate the bubble-net attacking strategy of whales. The mathematical model is shown as follow.

① Shrinking Encirclement

This behavior can be realized by reducing a in Equation (7), and it should be noted that $\overrightarrow{A}$ also decreases as a decreases.

② Updating position

With the spiral motion close to the prey, the mathematical model is shown as follow.

$$\overrightarrow{X}\left(t+1\right)={\overrightarrow{D}}^{\prime }{e}^{bl}\cos \left(2\pi l\right)+{\overrightarrow{X}}^{*}\left(t\right)$$

(9)

$${\overrightarrow{D}}^{\prime }=\left|{\overrightarrow{X}}^{*}-\overrightarrow{X}\left(t\right)\right|$$

(10)

where ${\overrightarrow{D}}^{\prime }$ is the best solution; b is the helical shape parameter; l is the random number of [0, 1].

In order to simulate these two behaviors of humpback whale hunting more accurately, the probability of contraction encirclement and spiral movement is considered to be 50%. Therefore, constructing Equation (11) to update the position of whales:

$$\overrightarrow{X}\left(t+1\right)=\{\begin{array}{cc}{\overrightarrow{X}}^{*}\left(t\right)-\overrightarrow{A}\overrightarrow{D}& p<0.5\\ {\overrightarrow{D}}^{\prime }{e}^{bl}\cos \left(2\pi l\right)+{\overrightarrow{X}}^{*}\left(t\right)& p\ge 0.5\end{array}$$

(11)

where p is the random number.

(3) Searching prey

The search for preys is a random search that uses the relative position of the whale and does not depend on the location of the prey. The mathematical model is

$$\overrightarrow{D}=\left|\overrightarrow{C}{\overrightarrow{X}}_{rand}\left(t\right)-\overrightarrow{X}\left(t\right)\right|$$

(12)

$$\overrightarrow{X}\left(t+1\right)={\overrightarrow{X}}_{rand}\left(t\right)-\overrightarrow{A}\overrightarrow{D}$$

(13)

where ${\overrightarrow{X}}_{rand}\left(t\right)$ is a randomly determined whale location.

According to the above mathematical model, when p > 0.5, the whales attack the prey. Otherwise the whales are in the search or siege phase. |A| ≥ 1 is a random search algorithm to update the location of the whale, and |A| < 1 refers to the stage of surrounding the prey after the maximum number of iterations to stop calculation.

3. WOA-VMD Denoising Method

3.1. VMD Parameter Optimization

In the process of VMD signal decomposition, the mode number K and penalty parameter α must be determined. If the value of K is excessive, the signal is decomposed too much, otherwise, the signal is underdecomposed. If the value of α is excessive, the signal component will be lost, otherwise, the information will be redundant [17]. In the actual analysis, the optimal parameters should be determined according to the characteristics of each signal. At present, the optimal K value is usually determined by observing the frequency variation characteristics of the center [13], but this method is subjective and accidental, and the optimal penalty parameter α cannot be obtained.

In this paper, the WOA algorithm is used to optimize the VMD parameters, and envelope entropy is used as the evaluation index of the optimization effect. When the signal contains more noise, the feature information is less, the entropy value is larger, and vice versa. Firstly, the position vector [K, α] of the whale group is initialized to calculate the fitness of each whale. Then, the formula is determined by judging the convergence size, and the optimal solution is obtained by iterating repeatedly until the end. The process is shown in Figure 1.

3.2. Denoising Method

Multi-scale permutation entropy (MPE) can be used to analyze signal mutation characteristics. The timing signal x = {x₁, x₂..., x_L} carries out multi-scale coarser granulation and reconstruction, and calculates the PE of signals at different scales. The steps are as follows [18]:

① Multi-scale coarsening of temporal signals.

In this step, the time series is divided into different scales and averaged in each segment so as to obtain the coarse-grained time series Y^s_t. The calculation formula is shown in Equations (14) and (15).

$$\begin{array}{cc}{y}_j^s=\frac{1}{s}{\displaystyle \sum _{i=(j-1)s+1}^{js}{x}_i}& 1\le j\le L\end{array}$$

(14)

$$Y_t^s=\left\{y_t^s,y_{t+\tau }^s,\cdots ,y_{t+(m-1)\tau }^s\right\}$$

(15)

where s is the scale factor; y^s_j is the time series; τ is the delay; m is the embedding dimension; Y^s_t is the reconstructed sequence.

② Reconstruction of temporal signals.

Reconstruct the Y^s_t ascending order to obtain the m! permutation. Moreover, the number of occurrences N₁ and probability P^s_l of each permutation were calculated.

$$y_{t+(j_1-1)\tau }^s\le y_{t+(j_2-1)\tau }^s\le \cdots \le y_{t+(j_m-1)\tau }^s$$

(16)

$$P_l^s=\frac{N_l}{n/s-m+1}$$

(17)

③ Calculation of permutation entropy at multiple scales

The permutation entropy is calculated at different scale factor s and normalized to obtain multi-scale permutation entropy.

$$h_P^s=H_P^s/\ln (m!)=\frac{-{\displaystyle \sum _{l=1}^{m!}{P}_l^s\ln }{P}_l^s}{\ln (m!)}$$

(18)

The improved algorithm is applied to decompose the measured signals, and the K components of IMF are obtained. The dominant frequency increases with the increase in K value. MPE values of different scales were calculated for each IMF, and noise components were identified by means. After removing the noise components, the IMF was reconstructed to obtain the denoised signal.

3.3. Evaluation of Noise Reduction Effect

The signal to noise ratio (s) and the root mean square error (ε) before and after noise reduction are taken as evaluation indexes, which can effectively evaluate the effect of the noise reduction algorithm [19]. The evaluation parameters are defined as follows:

(1) Signal-to-noise ratio s

$$s=10\lg \left(\frac{{\displaystyle \sum _{i=1}^L{(x_i)}^2}}{{\displaystyle \sum _{i=1}^L{(x_i-{\tilde{x}}_i)}^2}}\right)$$

(19)

(2) Root mean square error ε

$$\epsilon =\sqrt{\frac{1}{L}{\displaystyle \sum _{i=1}^L{({\tilde{x}}_i-x_i)}^2}}$$

(20)

where ${\tilde{x}}_i$ is the vibration amplitude of the noise reduction signal. s reflects the energy relationship between noise and signal. The larger s is, the fewer signal features will be destroyed during noise reduction. ε reflects the mean noise energy, and the smaller ε is, the better the noise reduction effect is.

4. Engineering Applications

4.1. Vibration Signal Noise Reduction

The Xiali tunnel is one of the key control projects of the Ji-Xin expressway in the Henan province, with a total length of 1020 m. The bedrock in the tunnel site is exposed and the stratigraphic lithology indicates the presence of medium-strong weathered limestone. The portal section of the Xiali tunnel is a steep slope, which is affected by combined joints, and the slope stability is poor. The blasting of the shallow buried section of the tunnel portal will cause great vibrations to the overlying slope.

Smooth blasting was used to excavate the tunnel on site. Detonator delays in the blasting have six sections, Ms1–Ms11. Due to the joint development of the surrounding rock and slope in the tunnel portal, the displacement and blasting vibration were monitored during the construction process. The blasting vibration of the tunnel portal is monitored by the vibration meter. Set the sampling frequency to 4000 Hz and the time to 1 s. A typical signal is selected for noise reduction, and the blasting vibration time history curve is shown in Figure 2. The Y axis represents particle vibration velocity, which is directly obtained by the sensors [20].

According to Figure 2, blasting vibration monitoring sensors are usually arranged on rock outcrop in the shallow buried section of a tunnel portal. The blasting vibration propagation is affected by a complex rock and soil layer structure, and the measured signal has a large amount of noise, which has a greater impact on signal feature identification.

Therefore, the WOA was used to optimize the VMD parameters [K, α], set the population size as 50, the maximum number of iterations as 30, the iteration range of K as [2, 15], and the iteration range of α as [1000, 10,000].

After 30 iterations, the minimum envelope entropy of the fifth iteration is 1.045, and the output optimal parameters are [K, α] = (8, 5600). The IMF component and its center frequency of signal VMD decomposition were obtained, and the vibration energy and percentage of the original signal and each IMF component were calculated. The results are shown in Figure 3 and

Table 1. Mean MPE of IMF.

Signals	Original Signal (s1)	IMF1	IMF2	IMF3	IMF4	IMF5	IMF6	IMF7	IMF8
Center frequency/Hz	25.41	3.25	30.68	77.66	219.31	966.44	1258.43	1512.42	1790.42
Vibration energy/(cm2/s2)	110.1172	31.3148	56.5706	21.0999	0.8941	0.0566	0.0622	0.0572	0.0618
Energy percentage %	100	28.438	51.373	19.161	0.812	0.051	0.056	0.052	0.056

.

According to Figure 3 and Table 1, eight IMF components are obtained from WOA-VMD decomposition, and the center frequencies of IMF1–IMF8 gradually increase. The vibration amplitude and energy of IMF4–IMF8 are low, accounting for less than 1% of the total energy. The center frequency of IMF5–IMF8 is above 900 Hz, which is obviously higher than the frequency of the general blasting vibration signal. Therefore, it can be basically judged that IMF4–IMF8 may be high-frequency noise.

In order to accurately identify the noise components, the mean MPE values of different scales of each IMF were calculated. The trial calculation determined the embedding dimension m = 6, the delay τ = 1, and the scale factor s = 5, and the mean MPE was obtained as shown in

Table 2. Mean value of MPE.

	IMF1	IMF2	IMF3	IMF4	IMF5	IMF6	IMF7	IMF8
MPE	0.4152	0.3774	0.4242	0.6337	0.7403	0.6990	0.6910	0.7186

.

According to the literature [21], the MPE threshold of the blasting vibration signal is 0.6, so IMF4–IMF8 is determined as signal noise. After removing IMF4–IMF8 and reconstructing other IMFs, the blasting vibration signal after noise reduction is obtained, as shown in Figure 4.

The noise reduction process shows that the WOA algorithm can adaptively obtain the VMD decomposition results and avoid the calculation error caused by human factors. Comparing Figure 2 and Figure 5, it can be seen that the waveform of the denoised signal is basically the same as that of the original signal, the noise is significantly eliminated, and the impact on the peak vibration velocity is small.

4.2. Energy Characteristic Analysis

Adaptive optimal kernel (AOK) time-frequency technology is suited to analyze the time-frequency characteristics of blasting vibration signals, which can effectively obtain the detailed characteristics of vibration signals [22,23]. In order to verify the noise reduction effect of the WOA-VMD algorithm, AOK time-frequency technology is used to obtain the signal time spectrum before and after noise reduction, as shown in Figure 5, where X is the energy and Y is the dominant frequency.

As Figure 5 shows, the dominant frequency of the signal before and after noise reduction is 26.4 Hz, the peak energy is 361.3 J and 358.9 J, the vibration duration is about 0.5 s, and the vibration energy is mainly distributed in the range of 0–250 Hz. By comparing Figure 5a,b, we can see the signal spectrum before noise reduction has significant noise points, as shown in the red box in the figures. These noise points are noise components with a frequency above 250 Hz, which are the main cause of signal distortion. After signal noise reduction, the original noise components are effectively eliminated. The dominant signal frequency before and after noise reduction is the same, and the peak energy is only reduced by 2.4 J. It shows that the noise reduction method can eliminate the high frequency noise effectively and has little influence on the low frequency vibration.

5. Comparison Analysis

In order to further analyze the noise reduction effect of the WOA-VMD algorithm, the CEEMDAN algorithm with adaptive ability, the traditional VMD algorithm and WOA-VMD algorithm are respectively applied to decompose three groups of measured signals induced by tunnel blasting in the portal section (as shown in Figure 6), and MPE is used to determine the noise composition. The parameters [K, α] of the three algorithms are shown in

Table 3. Denoised effect index.

Original Signal	Noise Reduction Algorithm	Modal Number K	Penalty Factor α	Signal to Noise Ratio s/dB	Root Mean Square Error ε
s1	CEEMDAN	10	-	14.021	0.095
	VMD	8	4800	18.242	0.057
	WOA-VMD	8	5600	18.322	0.055
s2	CEEMDAN	6	-	13.871	0.101
	VMD	6	4800	16.520	0.070
	WOA-VMD	7	6100	18.325	0.061
s3	CEEMDAN	8	-	10.210	0.125
	VMD	8	5000	14.525	0.081
	WOA-VMD	10	6000	16.658	0.072

, and the parameters of traditional VMD algorithm are all processed manually by the method used in the literature [13]. Figure 6 is the signal comparison before and after noise reduction, and the signal-to-noise ratio s and the root mean square error ε of the denoised signal and the original signal are calculated in Table 3.

According to the waveform curve (Figure 6), the blasting vibration signals of the tunnel portal before and after noise reduction were compared and analyzed. The CEEMDAN, VMD, and WOA-VMD algorithms all achieved a certain degree of noise reduction effect and removed part of the high-frequency noise, among which the WOA-VMD algorithm obtained the best waveform and the least noise. According to Table 3, comparing the noise reduction effect indexes of the three groups of signals, the signal-to-noise ratio s of the WOA-VMD and VMD algorithms is greater than that of the CEEMDAN algorithm, and the root mean square error ε of the WOA-VMD and VMD algorithm is smaller than that of the CEEMDAN algorithm. It shows that the denoising effect of the VMD algorithm is better than that of the CEEMD class decomposition algorithm.

Comparing the noise reduction effect of the WOA-VMD and VMD algorithms, it can be seen that there are differences in the mode number K and penalty factor α determined by the WOA optimization algorithm and manual determination for different signals. For signal s1, the two methods obtain the same mode number, but the penalty factor is slightly different, and the noise reduction effect is not different. However, for signals s2 and s3, the modal number K obtained by WOA-VMD decomposition are 1~2 more modal functions than that obtained by VMD algorithm, the decomposition is more refined, and the denoising effect is better.

From the analysis of the waveform and evaluation index, it is found that the denoising effect of the WOA-VMD algorithm is better than that of other methods, and the high-frequency noise in signal is more effectively removed. At the same time, the algorithm has strong adaptive energy, avoids the error of manual intervention, and has obvious superiority in the noise reduction of the blasting vibration signal at the tunnel portal.

6. Conclusions

Due to the influence of the complex environment at the tunnel portal, the blasting vibration signal of the slope in tunnel portal sections contains a lot of high-frequency noise, which will cover the characteristics of the real signal and is not conducive to further analysis of its propagation and attenuation law. Therefore, the WOA-VMD signal noise reduction method is proposed for the blasting vibration signal processing of tunnel portals.

(1) The measured signals were decomposed by VMD to obtain IMFs with different dominant frequencies, which avoided the problems of endpoint distortion and mode confusion. The WOA algorithm adaptively determines the VMD parameters [K, α], realizes the optimal decomposition, and avoids the randomness of subjective decision.

(2) The noise reduction analysis of the blasting signal at the tunnel portal shows that the accurate decomposition of the WOA-VMD algorithm and the noise recognition of MPE effectively improve the noise reduction effect. The AOK time-frequency energy spectra of vibration signals before and after noise reduction were compared and analyzed. The WOA-VMD noise reduction method removed the high-frequency noise energy and effectively retained the main energy of the signal. The denoised signal can be used to judge whether the vibration caused by blasting is under or over the accepted limits, which can determine the stability of the portal.

(3) The denoising methods of CEEMDAN, VMD, and WOA-VMD all have good fidelity and a certain denoising effect. The noise reduction signal waveform and evaluation index of the WOA-VMD algorithm are better than that of the CEEMDAN and VMD algorithms, which verifies the effectiveness of the novel method, and provides technical support for a detailed analysis of the blasting vibration propagation law and harmful effects at tunnel entrances.