In underwater drilling and blasting engineering, the blasting vibration signal is mixed with a mass of noises due to the complexity of monitoring environment, the error of monitoring sensors and the reflection of propagation medium. In order to accurately obtain the characteristics of vibration signal, a novel denoising model is established. The complete ensemble empirical mode decomposition with adaptive noise is used to decompose the original signal, and the objective function of the filtering algorithm is used to obtain the optimal denoising signal. The results indicate that the model can not only successfully remove the high-frequency noise but also has no effect on the low-frequency signal components, which verifies the reliability and validity of the denoising model. 在水下钻孔爆破工程中,由于监测环境的复杂性、监测传感器的误差和传播介质的反射,爆破振动信号中混杂着大量的噪声。为了准确获取振动信号的特征,本文建立了一个新颖的去噪模型。利用带有自适应噪声的完全集合经验模态分解法对原始信号进行分解,并利用滤波算法的目标函数获得最优去噪信号。结果表明,该模型不仅能成功去除高频噪声,而且对低频信号成分没有影响,验证了去噪模型的可靠性和有效性。
Keywords Underwater drilling and blasting *\cdot Vibration signal *\cdot Denoising model *\cdot CEEMDAN 关键词 水下钻孔和爆破 *\cdot 振动信号 *\cdot 去噪模型 *\cdot CEEMDAN
List of Symbols 符号列表
IMF Intrinsic mode function IMF 固有模式功能 x(t)quadx(t) \quad Original signal x(t)quadx(t) \quad 原始信号 epsi_(0)quad\varepsilon_{0} \quad Noise coefficient epsi_(0)quad\varepsilon_{0} \quad 噪声系数 r(t)quadr(t) \quad Residual component r(t)quadr(t) \quad 剩余部分 n_(j)(t)quadn_{j}(t) \quad White noises n_(j)(t)quadn_{j}(t) \quad 白噪声 x_(m)quad mx_{m} \quad m th sample point of original signal x_(m)quad mx_{m} \quad m 原始信号的第 1 个采样点 tilde(x)_(m)quad M\tilde{x}_{m} \quad M th sample point of denoised signal 去噪信号的 tilde(x)_(m)quad M\tilde{x}_{m} \quad M 个采样点 MSE_(f)quad\operatorname{MSE}_{f} \quad Mean square error MSE_(f)quad\operatorname{MSE}_{f} \quad 均方误差 u(x),v(x),f(x)quadu(x), v(x), f(x) \quad Smooth curves u(x),v(x),f(x)quadu(x), v(x), f(x) \quad 平滑曲线 h quadh \quad Sampling interval h quadh \quad 采样间隔
SMSE _(f)quad_{f} \quad Mean square error of smoothness SMSE _(f)quad_{f} \quad 平滑度的均方误差 F quadF \quad Objective function F quadF \quad 目标函数 muquad\mu \quad Weight coefficient muquad\mu \quad 重量系数 xiquad\xi \quad Signal-to-noise ratio xiquad\xi \quad 信噪比 epsiquad\varepsilon \quad Mean absolute error epsiquad\varepsilon \quad 平均绝对误差
1 Introduction 1 引言
In the constructions of wading engineering, such as harbors, wharfs and channels, underwater drilling and blasting plays a critical role in the excavation of rock mass. In such engineering, merely 20-30%20-30 \% energy released by the explosion of explosive is used to break rock mass, while the main part adversely affects the surrounding environment [1-3]. Blasting vibration is one of the main harmful effects, which will adversely affect the buildings, bridges and slopes near the blasting area. 在港口、码头、航道等涉水工程建设中,水下钻孔爆破对岩体开挖起着至关重要的作用。在这类工程中,仅仅利用炸药爆炸释放的 20-30%20-30 \% 能量来破碎岩体,而主要部分则对周围环境造成不利影响[1-3]。爆破振动是主要的有害影响之一,会对爆破区域附近的建筑物、桥梁和斜坡造成不利影响。
Blasting vibration signal is the basis for analyzing blasting vibration effects, which reflects the dynamic response characteristics of structures. Noise is an inevitable problem for signals in the process of vibration monitoring due to the complexity of monitoring environment, the error of monitoring sensors, the reflection of propagation medium and the interference of magnetic field [4]. The high-frequency noise will distort blasting vibration signals and conceal the real information of signal, which result in the inaccurate description of vibration characteristics. Therefore, blasting vibration signals must be denoised. 爆破振动信号是分析爆破振动效应的基础,它反映了结构的动态响应特性。在振动监测过程中,由于监测环境的复杂性、监测传感器的误差、传播介质的反射以及磁场的干扰,噪声是信号不可避免的问题[4]。高频噪声会使爆破振动信号失真,掩盖信号的真实信息,导致振动特征描述不准确。因此,必须对爆破振动信号进行去噪处理。
Blasting vibration signal is a sort of non-stationary random signal, and the denoising methods for the signal mainly contain the wavelet decomposition technology [5-7] and the empirical mode decomposition (EMD) technology [8-10]. Wavelet decomposition technology is characterized by multi-resolution analysis and excellent time-frequency 爆破振动信号是一种非稳态随机信号,其去噪方法主要包括小波分解技术[5-7]和经验模态分解(EMD)技术[8-10]。小波分解技术的特点是多分辨率分析和出色的时频分析能力。
Springer 施普林格
localization capability. The real signal and noises can be separated according to different characteristics of wavelet coefficients. Zhang et al. [11] applied the wavelet threshold denoising method to process the blast vibration signal, and results demonstrated that noise could be removed by Rigrsure threshold based on Stein unbiased likelihood estimation. Xie et al. [12] and Xia et al. [13], respectively, proposed different types of wavelet threshold denoising algorithms and applied the algorithms to process the measured vibration signals. However, wavelet denoising methods have to choose basis functions, which inevitably generates defects as Fourier transform. 定位能力根据小波系数的不同特性,可以分离真实信号和噪声。Zhang 等[11]应用小波阈值去噪方法处理爆炸振动信号,结果表明基于 Stein 无偏似然估计的 Rigrsure 阈值可以去除噪声。Xie 等人[12]和 Xia 等人[13]分别提出了不同类型的小波阈值去噪算法,并将算法应用于处理实测振动信号。然而,小波去噪方法必须选择基函数,这不可避免地会产生傅里叶变换的缺陷。
EMD is characterized by multi-resolution and excellent adaptability to non-stationary signals, which is widely applied to process vibration signals. Zhao et al. [14] employed the ensemble empirical mode decomposition (EEMD) to denoise blasting vibration signals and compared it with wavelet denoising method. Yuan et al. [15] applied EMD denoising method and Hilbert transform to analyze the blasting vibration signals in lead-zinc mine. EMD is adaptive and need not to select basis functions, but the method will result in mode mixing when processing abnormal noise signals [16]. EMD具有多分辨率、对非平稳信号适应性强等特点,被广泛应用于振动信号的处理。Zhao 等[14]采用集合经验模态分解(EEMD)对爆破振动信号进行去噪,并与小波去噪方法进行了比较。Yuan等人[15]应用EMD去噪方法和希尔伯特变换分析了铅锌矿的爆破振动信号。EMD 是自适应的,无需选择基函数,但该方法在处理异常噪声信号时会导致模态混合[16]。
In order to accurately analyze the vibration characteristics of underwater drilling and blasting, it is necessary to reduce the noise of the measured signals. The existing denoising methods still have defects and are not effective in more complex underwater environment. Therefore, it is necessary to establish a denoising method suitable for vibration signal induced by underwater drilling and blasting. 为了准确分析水下钻孔和爆破的振动特性,有必要降低测量信号的噪声。现有的去噪方法仍存在缺陷,在较为复杂的水下环境中效果不佳。因此,有必要建立一种适合水下钻孔和爆破引起的振动信号的去噪方法。
The complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) can not only reduce the phenomenon of mode mixing but also accurately reconstruct the original signal [17,18][17,18]. The optimal smooth denoising algorithm can get better filtering signals. This paper proposes a novel denoising model to process vibration signals induced by underwater drilling and blasting, and the model is consisted of above two parts. CEEMDAN is used to decompose the original signal, and the objective function of the filtering algorithm is used to obtain the optimal denoised signal. The novel denoising model can effectively improve the denoising effect and preserve the authenticity and integrity of waveform better. 具有自适应噪声的完全集合经验模态分解(CEEMDAN)不仅能减少模态混合现象,还能准确地重建原始信号 [17,18][17,18] 。最优平滑去噪算法可以获得更好的滤波信号。本文提出了一种处理水下钻孔爆破引起的振动信号的新型去噪模型,该模型由以上两部分组成。利用 CEEMDAN 对原始信号进行分解,并利用滤波算法的目标函数获得最优去噪信号。该新型去噪模型能有效提高去噪效果,更好地保持波形的真实性和完整性。
2 CEEMDAN Theory 2 CEEMDAN 理论
2.1 Empirical Mode Decomposition 2.1 经验模式分解
Empirical mode decomposition (EMD) is an adaptive signal decomposition method for processing the nonlinear and non-stationary signals proposed by Huang et al. [19]. Based on the time-scale characteristics of signals, the 经验模态分解(EMD)是 Huang 等人[19]提出的一种处理非线性和非平稳信号的自适应信号分解方法。根据信号的时间尺度特征,EMD
multi-component signals can be decomposed into a series of intrinsic mode function (IMF) components and residual components, and the IMF components are arranged in the order of instantaneous frequencies from high to low. The method possesses favorable adaptability, completeness and orthogonality. However, there is mode mixing problem in processing signals containing discontinuities, impulses and noises [20]. 多分量信号可被分解为一系列本征模态函数(IMF)分量和残差分量,IMF分量按瞬时频率从高到低的顺序排列。该方法具有良好的适应性、完整性和正交性。然而,在处理含有不连续、脉冲和噪声的信号时,存在模式混合问题[20]。
To solve the above problems, Huang et al. [21] added the even-distributed white noises to the original signal and decomposed the aggregate signals into IMFs by EMD method for several times, and the final IMF was obtained by the averaged IMFs. This method is called the ensemble empirical mode decomposition (EEMD), which can weaken the mode mixing phenomenon of EMD algorithm to a certain extent. The algorithm tries to eliminate the influence of white noise on decomposition results through integration averaging, but it cannot be entirely suppressed. The reconstruction error is highly dependent on the number of integration times. 为了解决上述问题,Huang 等人[21]在原始信号中加入偶数分布的白噪声,用 EMD 方法将集合信号多次分解为 IMF,通过平均 IMF 得到最终的 IMF。这种方法被称为集合经验模态分解(EEMD),可以在一定程度上削弱 EMD 算法的模态混合现象。该算法试图通过积分平均消除白噪声对分解结果的影响,但并不能完全抑制白噪声。重构误差与积分次数有很大关系。
2.2 CEEMDAN
Torres et al. [17] proposed the complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) based on EEMD. Several adaptive white noises are added to the original signal during each EMD decomposition. Even if the integration times are limited, the reconstruction error is almost 0 ; that is, the reconstructed signal is almost identical to the original signal. The main steps of CEEMDAN are as follows. Torres 等人[17] 在 EEMD 的基础上提出了带有自适应噪声的完整集合经验模式分解(CEEMDAN)。在每次 EMD 分解过程中,都会向原始信号添加若干自适应白噪声。即使积分时间有限,重建误差也几乎为 0;也就是说,重建后的信号与原始信号几乎相同。CEEMDAN 的主要步骤如下。
Add the white noises n_(j)(t)n_{j}(t) with different amplitudes to the original signal x(t)x(t). Then, the signal can be expressed as x(t)+epsi_(0)n_(j)(t)x(t)+\varepsilon_{0} n_{j}(t), and epsi_(0)\varepsilon_{0} is the noise coefficient. EMD is used to carry out II times of decomposition, and the first IMF can be obtained by the integrated averaging. The first IMF and the first residual component are shown as follows. 在原始信号 x(t)x(t) 中加入不同振幅的白噪声 n_(j)(t)n_{j}(t) 。然后,信号可以表示为 x(t)+epsi_(0)n_(j)(t)x(t)+\varepsilon_{0} n_{j}(t) , epsi_(0)\varepsilon_{0} 是噪声系数。利用 EMD 进行 II 次分解,通过积分平均可以得到第一个 IMF。第一个 IMF 和第一个残差分量如下所示。 IMF_(1)(t)=(1)/(I)sum_(i=1)^(I)IMF_(i1)(t)\operatorname{IMF}_{1}(t)=\frac{1}{I} \sum_{i=1}^{I} \operatorname{IMF}_{i 1}(t) r_(1)(t)=x(t)-IMF_(1)(t)r_{1}(t)=x(t)-\mathrm{IMF}_{1}(t). EMD_(j)(*)\mathrm{EMD}_{j}(\cdot) is the jj th modal component generated by the EMD algorithm. After the signal r_(1)(t)+epsi_(1)*EMD_(1)(n_(j)(t))r_{1}(t)+\varepsilon_{1} \cdot \operatorname{EMD}_{1}\left(n_{j}(t)\right) is further decomposed for II times, the second IMF can be obtained as follows. EMD_(j)(*)\mathrm{EMD}_{j}(\cdot) 是 EMD 算法生成的 jj 第三模态分量。在对信号 r_(1)(t)+epsi_(1)*EMD_(1)(n_(j)(t))r_{1}(t)+\varepsilon_{1} \cdot \operatorname{EMD}_{1}\left(n_{j}(t)\right) 进一步分解 II 次后,可得到第二个 IMF 如下。 IMF_(2)(t)=(1)/(I)sum_(i=1)^(I)EMD_(1)[r_(1)(t)+epsi_(1)EMD_(1)(n_(i)(t))]\operatorname{IMF}_{2}(t)=\frac{1}{I} \sum_{i=1}^{I} \operatorname{EMD}_{1}\left[r_{1}(t)+\varepsilon_{1} \operatorname{EMD}_{1}\left(n_{i}(t)\right)\right]
Calculate the residual components in the sequence of 2, 3,dots,k3, \ldots, k, which is shown as follows. 计算 2、 3,dots,k3, \ldots, k 序列中的残差分量,如下所示。
Extract the first IMF from the signal r_(1)(t)+epsi_(1)*EMD_(1)(n_(j)(t))r_{1}(t)+\varepsilon_{1} \cdot \operatorname{EMD}_{1}\left(n_{j}(t)\right), and the IMF_(k+1)\mathrm{IMF}_{k+1} is obtained. 从信号 r_(1)(t)+epsi_(1)*EMD_(1)(n_(j)(t))r_{1}(t)+\varepsilon_{1} \cdot \operatorname{EMD}_{1}\left(n_{j}(t)\right) 中提取第一个 IMF,得到 IMF_(k+1)\mathrm{IMF}_{k+1} 。 IMF_(k+1)(t)=(1)/(I)sum_(i=1)^(I)EMD_(k)[r_(k)(t)+epsi_(k)EMD_(k)(n_(k)(t))]\operatorname{IMF}_{k+1}(t)=\frac{1}{I} \sum_{i=1}^{I} \operatorname{EMD}_{k}\left[r_{k}(t)+\varepsilon_{k} \operatorname{EMD}_{k}\left(n_{k}(t)\right)\right]
Until the residual signal can no longer be decomposed, all the IMFs can be obtained, and the final residual component is shown as follows. 直到残差信号无法再被分解为止,所有的 IMF 都可以得到,最终的残差分量如下所示。 r(t)=x(t)-sum_(k=1)^(K)IMF_(k)(t)r(t)=x(t)-\sum_{k=1}^{K} \mathrm{IMF}_{k}(t)
The original signal x(t)x(t) can be expressed as follows. 原始信号 x(t)x(t) 可以表示如下。 x(t)=r(t)+sum_(k=1)^(K)IMF_(k)(t)x(t)=r(t)+\sum_{k=1}^{K} \operatorname{IMF}_{k}(t)
CEEMDAN makes full use of the noise-assisted analysis and can reconstruct the original signal accurately and completely. And adjusting the noise coefficient epsi_(k)\varepsilon_{k}, the noises with different ratios are selected during each decomposition to calculate IMFs. CEEMDAN 充分利用噪声辅助分析技术,可以准确、完整地重建原始信号。通过调整噪声系数 epsi_(k)\varepsilon_{k} ,在每次分解时选择不同比例的噪声来计算 IMF。
3 Optimal Smooth Denoising Model 3 最佳平滑去噪模型
A low-pass filtering algorithm can be obtained by IMFs with CEEMDAN decomposition. The optimal smooth denoising model is established by combining the smoothness and deviation of the filtering algorithm [22]. 低通滤波算法可以通过 IMF 与 CEEMDAN 分解得到。结合滤波算法的平滑度和偏差,可以建立最佳平滑去噪模型[22]。
3.1 Deviation of the Filtering Algorithm 3.1 过滤算法的偏差
After the original signal x(t)x(t) is decomposed by CEEMDAN, the low-pass filtering algorithm can be established, as shown in Eq. (8). 原始信号 x(t)x(t) 经 CEEMDAN 分解后,可建立低通滤波算法,如式(8)所示。 tilde(x)_(k)(t)=x(t)-sum_(1)^(k)IMF_(k)1 <= k <= K\tilde{x}_{k}(t)=x(t)-\sum_{1}^{k} \operatorname{IMF}_{k} 1 \leq k \leq K
Then, the mean square error between the denoised signal and the original signal is shown as follows. 然后,去噪信号与原始信号之间的均方误差如下所示。 MSE_(f)=sqrt((sum_(m=1)^(M)( tilde(x)_(m)-x_(m))^(2))/(M))\mathrm{MSE}_{f}=\sqrt{\frac{\sum_{m=1}^{M}\left(\tilde{x}_{m}-x_{m}\right)^{2}}{M}}
where MM is the number of sample points; x_(m)x_{m} is the value of the mm th sample point of original signal; and tilde(x)_(m)\tilde{x}_{m} is the value of the mm th sample point of denoised signal. MSE_(f)\mathrm{MSE}_{f} reflects the closeness degree between filtered data and original 其中, MM 为采样点的个数; x_(m)x_{m} 为原始信号的第 mm 个采样点的值; tilde(x)_(m)\tilde{x}_{m} 为去噪信号的第 mm 个采样点的值。 MSE_(f)\mathrm{MSE}_{f} 反映了滤波数据与原始数据的接近程度。
data, which is used to evaluate the quality of the filtering algorithm. 数据,用于评估过滤算法的质量。
3.2 Smoothness of the Filtering Algorithm 3.2 滤波算法的平滑性
Smooth curve is that the function x(t)x(t) has a first continuous derivative in the defined interval. For all the points in a certain curve, the left and right derivatives exist and are equal. For a smooth curve formed by more than two curves, the derivatives at the junction of the curves are equal. 平滑曲线是指函数 x(t)x(t) 在定义的区间内具有一阶连续导数。对于某条曲线上的所有点,左导数和右导数都存在且相等。对于由两条以上曲线构成的光滑曲线,曲线交点处的导数相等。
Suppose that the smooth curve formed by two curves u(x)u(x) and v(x)v(x) joins at the point x_(0)x_{0} and the curves are continuous and derivable in the point x_(0)x_{0}, the curves must satisfy Eq. (10). 假设由两条曲线 u(x)u(x) 和 v(x)v(x) 形成的光滑曲线在点 x_(0)x_{0} 相接,且曲线连续并可在点 x_(0)x_{0} 中导出,则曲线必须满足公式 (10)。 {:[u(x_(0)),=v(x_(0))],[u^(')(x_(0)),=v^(')(x_(0))]}\left.\begin{array}{rl}u\left(x_{0}\right) & =v\left(x_{0}\right) \\ u^{\prime}\left(x_{0}\right) & =v^{\prime}\left(x_{0}\right)\end{array}\right\}
The smooth curve should satisfy that the curvatures of u(x)u(x) and v(x)v(x) at the point x_(0)x_{0} are equal. 平滑曲线应满足 u(x)u(x) 和 v(x)v(x) 在 x_(0)x_{0} 点的曲率相等。 {:[K_(u)|_(x=x_(0))=K_(v)|_(x=x_(0))],[(|u^('')(x_(0))|)/({1+[u^(')(x_(0))]^(2)}^(3//2))=(|v^('')(x_(0))|)/({1+[v^(')(x_(0))]^(2)}^(3//2))]:}\begin{aligned} & \left.K_{u}\right|_{x=x_{0}}=\left.K_{v}\right|_{x=x_{0}} \\ \frac{\left|u^{\prime \prime}\left(x_{0}\right)\right|}{\left\{1+\left[u^{\prime}\left(x_{0}\right)\right]^{2}\right\}^{3 / 2}}= & \frac{\left|v^{\prime \prime}\left(x_{0}\right)\right|}{\left\{1+\left[v^{\prime}\left(x_{0}\right)\right]^{2}\right\}^{3 / 2}}\end{aligned}
Substitute Eq. (10) into Eq. (11). 将公式 (10) 代入公式 (11)。 u^('')(x_(0))=v^('')(x_(0))u^{\prime \prime}\left(x_{0}\right)=v^{\prime \prime}\left(x_{0}\right)
By the definition of derivative, uu " (x)(x) and v"(x)v "(x) can be represented as follows. 根据导数的定义, uu " (x)(x) 和 v"(x)v "(x) 可以表示如下。 u^('')(x_(0))=lim_(h rarr0)(u(x_(0)-2h)-2u(x_(0)-h)+u(x_(0)))/(h^(2))u^{\prime \prime}\left(x_{0}\right)=\lim _{h \rightarrow 0} \frac{u\left(x_{0}-2 h\right)-2 u\left(x_{0}-h\right)+u\left(x_{0}\right)}{h^{2}}