基于CEEMDAN-小波包分析的隧道爆破信号去噪方法*

2021-05-27王海龙王海军彭婵媛

爆炸与冲击 2021年5期

王海龙，赵岩，王海军，彭婵媛，仝潇

（1. 中国矿业大学（北京）力学与建筑工程学院，北京 100083；2. 河北省土木工程诊断、改造与抗灾重点实验室，河北张家口 075000；3. 北旺建设集团有限公司，河北承德 067400）

作为铁路、公路隧道掘进最常用的破岩方式，爆破施工带来的结构振动极有可能对在建隧道及周围构筑物产生损害。目前，对爆破振动的分析研究大多是通过对爆破振动信号的处理展开的。隧道施工现场环境复杂，从现场获得的原始数据往往夹杂着各种各样的无用信号（噪声），尤其随着隧道掘进深度的增加，施工现场各个工序交叉进行，噪声与有效信号相互交织，对爆破信号处理分析产生极大的不利影响[1]。因此，对采集的原始隧道爆破振动数据进行降噪预处理显得尤为重要。

目前，应用于爆破振动信号消除干扰分量的方法主要包括小波阈值方法[2]、小波包阈值方法[3]、最小二乘拟合方法[4]、经验模态分解（empirical mode decomposition, EMD）方法[5]、集合经验模态分解(ensemble empirical mode decomposition, EEMD)方法[6]、变分模态分解（variational mode decomposition,VMD）方法[7]、EMD-小波阈值联合方法及EEMD-小波阈值联合方法[8-9]等。其中针对爆破信号中的高频噪声干扰分量，傅里叶分析、小波分析及小波包分析均可应用于去除噪声的领域；同时，EMD、EEMD 以及其与小波分析、小波包分析的联合方法也可以取得更为良好的效果。针对信号中存在的低频趋势项，最小二乘法、EMD 分解、EEMD 分解及VMD 分解等方法都能在保持信号特征的同时较为准确地去除趋势项。

作为爆破振动信号降噪的常用方法，小波降噪首先通过选择合适的小波基函数及分解层数对原始信号进行分解，然后选择合适的阈值函数对小波分解系数进行量化，对于超过阈值的小波分解系数予以保留，对于小于阈值的小波分解系数通过置零操作予以消除；然后，将小波系数进行重构即完成了信号消噪。小波包分析作为对小波算法的改进算法，通过对小波变换过程中忽略的高频部分进行精细分解，提高了小波算法的高频分辨率[10-11]。但是，小波及小波包算法均受制于小波基函数、阈值类型及阈值函数的选择，给复杂的隧道爆破振动信号的降噪处理带来困难。

EMD 分解[12]无需像小波分析一样设定小波基函数及阈值函数，就能将信号分解为一系列本征模态分量（IMF）。EEMD 分解针对EMD 分解中存在的模态混叠现象[13]，对原始信号添加白噪声信号进行分析。许多学者通过经验模态分解去除高频模态分量的方式实现去噪，然而这种“一刀切”式的降噪方式会在降噪的同时消除信号中的主要成分，造成信号失真。

小波分析与经验模态类分解的联合去噪方法可以综合两种方法的优点，在消噪的同时保持原始信号的完整性和特殊性；EMD-小波阈值降噪、EEMD-小波阈值降噪等方法已经应用于爆破振动信号的处理中。然而，EMD-小波阈值降噪、EEMD-小波阈值等降噪方法仍存在不少缺陷，例如：虽然EEMD 分解通过在原始信号中添加高斯白噪声解决了EMD 分解存在的模态混叠的问题，但是EEMD 并未对添加的白噪声信号进行隔离，造成残留噪声由高频过渡到低频，影响降噪效果。

隧道爆破振动信号中的噪声分量主要来自于钻孔装药、初期支护及运输渣土等施工工序中产生的高频干扰信号；同时，信号采集过程中，仪器温度变化造成的零点漂移及测试仪器的安装松动也会造成采样信号产生基线偏离及形状不规则。针对上述隧道爆破信号处理中存在的噪声干扰问题，本文引入一种基于CEEMDAN 分解及小波包分析的联合去噪方法：首先，对采集得到的原始信号进行CEEMDAN分解，获取多个本征模态分量（IMF）；然后，利用相关系数筛选出包含噪声的模态分量[14]，并对含噪分量进行小波包阈值降噪处理；最后，将经过处理的模态分量与未处理的分量重构得到纯净的爆破振动信号。通过CEEMDAN 分解来解决添加白噪声由高频传递到低频的问题[15-17]，提高高频成分的频率分辨率，以期更为精细地去除噪声成分。

1 CEEMDAN-小波包分析降噪原理

1.1 CEEMDAN 分解原理阐述

EEMD 与CEEMD 分解算法均是在待分解信号中加入高斯白噪声或者成对正负高斯白噪声来消除EMD 分解的模态混叠问题。然而，这两种算法均未对残留噪声进行隔离，造成添加的白噪声信号总是可以从高频转移到低频。因此，分解得到的模态本征分量中总会残留一定的白噪声信号，影响后续信号分析处理。

CEEMDAN 分解从两个方面解决了上述问题：(1) 加入经EMD 分解后含辅助噪声的IMF 分量，而不是将高斯白噪声信号直接添加在原始信号中；(2) EEMD 分解和CEEMD 分解是将经验模态分解后得到的模态分量进行总体平均，CEEMDAN 分解与之不同，则是在得到的第一阶IMF 后就进行总体平均计算，得到最终的第一阶IMF，然后对残余部分重复进行如上操作，这样便有效地解决了噪声从高频到低频的转移传递问题。接下来阐述CEEMDAN 分解原理。

（2）对产生的N 个模态分量进行加总平均得到：

（3）计算去除第1 个模态分量后的残差：

（4）在r1(t)中加入正负成对高斯白噪声得到新信号，以新信号为载体进行EMD 分解，得到第1 阶模态分量，由此可以得到第2 阶模态分量：

（5）计算去除第二个模态分量的残差：

（6）重复上述步骤，直到获得的残差信号为单调函数，不能继续分解，算法结束。此时得到的本征模态分量数量为K，则原始信号z(t)被分解为：

1.2 小波包降噪原理

小波降噪是通过小波变换对分解得到的小波系数进行分类处理。经小波分解后，有用信号的小波分解系数较大，噪声的小波分解系数较小。因此，通过设置阈值函数进行去噪处理，保留分解系数大于阈值的信号，而通过置零操作消除分解系数小于阈值的信号。小波包分析以小波变换为基础，对小波变换中未处理的高频分量再次细化分解，相比于小波降噪，小波包分析具有更高的频率分辨率，可以进一步消除高频部分存在的噪声余量，提高去噪精度。

可以得到小波包分解结果：

1.3 CEEMDAN-小波包分析联合降噪

图1 为CEEMDAN-小波包分析联合降噪流程图。如图所示，首先，利用CEEMDAN 算法将爆破振动信号分解，得到多个本征模态分量；然后，通过比较各个模态分量与原始信号的相关系数筛选出包含噪声的本征模态分量，并通过各本征模态分量的频谱分布特征及方差贡献率进行校核；最后，利用小波包分析对筛选出含有噪声的模态分量进行阈值去噪，将去噪完成的模态分量与未经处理的模态分量重构得到纯净信号。

图 1 去噪流程Fig. 1 Flow of de-noising

2 仿真信号分析

为验证本文引入的CEEMDAN-小波包阈值降噪方法的可行性，且不失一般性，首先对仿真信号进行降噪分析。近年来，有学者研究发现可以近似地利用正余弦函数模拟爆破振动信号进行分析[19]。因此，本文亦利用正余弦函数叠加得到的仿真信号进行研究。仿真信号表达式如下：

针对上述模拟信号，利用CEEMDAN 分解将其分解为一系列本征模态分量，并利用MATLAB 中互相关函数的概念计算得到各个模态分量与原始信号的相关性，并依据相关性绝对值的大小筛选得到需要进行处理的含噪分量。仿真信号z(t)及CEEMDAN 分解得到的本征模态分量（C1～C12）见图2。

图 2 仿真信号及模态分量波形图Fig. 2 Simulation signal and modal component waveform

图 3 仿真信号及降噪处理后的纯净信号Fig. 3 Simulated signal and pure signal after de-noising

图 4 草帽山隧道进口工区[20]Fig. 4 Caomaoshan tunnel entrance area[20]

3 工程实例分析

3.1 工程背景

本文以新建京张高铁草帽山隧道爆破施工为工程背景进行研究。草帽山隧道地形起伏较大，穿越草帽山主脉，穿越第四系上更新统新黄土，粉质黏土地等地层，穿越全风化凝灰岩，具有中等膨胀性，对隧道爆破施工造成较大的安全隐患。本文使用的原始数据来自于此隧道进口段钻爆法施工产生的爆破振动速度信号。图4 为隧道进口工区场地布置图。

3.2 爆破振动速度信号采集

此次爆破信号采集使用中科测控公司研发的TC-4850 爆破测振仪。该爆破测振仪可以同时采集X、Y 及Z 方向的爆破振动速度，同时利用Blasting vibration analysis 的分析软件可以对采集到的振速数据进行初步分析。采集过程中，利用配置好的不锈钢夹片将速度传感器固定于隧道初期支护上。其中X 方向朝向掌子面的掘进方向，Y 方向朝向隧道径向，Z 方向垂直于XY 平面垂直向上。测点布置见图5。

图 5 测点布置[21]Fig. 5 Layout of measuring points

图6 为某次爆破施工采集到的Z 方向爆破振动速度原始信号，从图6 可知，受爆破施工现场复杂环境的影响，爆破振速时程曲线中夹杂着各种毛刺噪声信号。

图 6 爆破振动速度原始信号Fig. 6 Original signal of blasting vibration speed

3.3 CEEMDAN-小波包阈值降噪实现

经CEEMDAN 分解，原始爆破信号被分解为14 个本征模态分量。依据EEMD 分解的取值范围，利用“试错法”得到本次试验中CEEMDAN分解的特征参数为：正负高斯白噪声标准差为0.2，加入噪声的次数为100，允许的最大筛选迭代次数为3 000。

利用式（12）计算得到每个模态分量与原始速度信号的相关系数（ri），计算结果见表1。

式中：xi表示各本征模态分量，y 表示原始爆破信号。

表 1 本征模态分量（IMF）的相关系数Table 1 Correlation coefficients of modal components (IMF)

从表1 中可以看出，各个本征模态分量与原始信号的相关系数差异较大。其中C6、C7、C8、C9的相关系数均大于0.2，可以看作是优势模态分量；而C1～C5及C10～C14的相关系数均小于0.14，初步认定C1～C5及C10～C14这两组模态分量中含有较多的噪声信号。

对C1～C5及C10～C14两组分量进行快速傅里叶分析得到频谱见图7。从图7 可以看出，C1～C3分量的主频主要集中在800～1 000 Hz，C4分量的主频集中在400 ～600 Hz，C5分量的主频主要分布在300 Hz左右；而C10～C14分量的主频率均小于5 Hz。由于本文振动测试仪的采集频率主要范围为5～250 Hz，因此从频谱分布特征来看，C1～C5及C10～C14分量确实含有某些干扰信号。

定义模态分量平方的算术平均值与其均值的平方之差为模态分量的方差（D(j)），通过各个模态分量的方差贡献率（e(j)）[22-23]校核上述选择的可行性。模态分量的方差及方差贡献率计算如下：

计算得到各个模态分量的方差贡献率见表2。由表2 可知，C6～C9模态分量的方差贡献率较大，均大于8，而剩余模态分量的方差贡献率较小，其中贡献率最大的C10也只有1.59。因此可以看出，由方差贡献率得到的结论与相关系数筛选及频谱分析结果保持一致。综上，可以认定C1～C5和C10～C14两组分量包含噪声信号。

将C1～C5及C10～C14重构得到新的分量C15，以C15为新信号进行小波包分析。由于本文中爆破测振仪TC-4 850 设置的最小工作采样频率为5 Hz，根据采样定理[24]，信号采集频率为5 000 Hz，Nyquist 频率为2 500 Hz，因此根据小波包分解算法，利用具有良好的紧支撑性、光滑性及近似对称性的“db8”小波基函数对新信号进行8 层分解，对应的最低频带为0～9.765 625 Hz；然后，利用默认的软阈值函数进行去噪处理。图8 为小波包降噪处理前后C15信号时程，从图8 可以看出，经过小波包降噪重构后的C15分量的噪声信号基本被消除。

图 7 C1～C5 及C10～C14 分量频谱Fig. 7 Spectra of C1-C5 and C10-C14 components

表 2 模态分量（IMF）的方差贡献率Table 2 Variance contribution rate of modal component (IMF)

将经过降噪处理的C15与优势分量C6、C7、C8、C9重构得到最终的纯净信号如图9 所示。由图9 可知，与原始信号相比，经CEEMDAN-小波包降噪后的纯净信号在保证局部波形特征及峰值不变的基础上，基本消除了原始信号中存在的噪声分量。

图 8 小波包降噪处理前后的C15 信号Fig. 8 signal (C15) before and after wavelet packet noise reduction

图 9 CEEMDAN-小波包阈值降噪后的信号Fig. 9 CEEMDAN-wavelet packet threshold signal after noise reduction

3.4 降噪效果评价

爆破振动信号降噪效果常利用信噪比（η）、均方根差（σ）进行衡量，信噪比越大，均方根差越小，去噪效果越好。信噪比、均方根差的计算公式[25-26]如下：

为了评价本文方法的降噪效果，利用小波包降噪、EMD-小波包联合降噪及EEMD-小波包联合降噪三种方法分别对原始爆破信号进行处理。四种方法的降噪指标见表3，降噪后的时程曲线如图10 所示。

对比四种方法的客观降噪指标发现，经验模态分解与小波包分析联合降噪的总体效果均优于小波包降噪。同时，与EMD-小波包联合降噪及EEMD-小波包联合降噪相比，CEEMDAN-小波包降噪方法得到的信噪比（94.080 2）最大，得到的均方根差（2.400 1×10-5）最小。对比图9 及图10 可以看出，与另外三种方法相比，经CEEMDAN-小波包阈值降噪后的纯净信号在保留原始信号完整信息的同时，可以更有效地去除信号中的噪声分量，有助于提高后续对隧道爆破信号分析的准确性。

表 3 去噪效果对比Table 3 Comparison of noise reduction effects

图 10 用于对比的几种方法的降噪效果Fig. 10 Noise reduction effect of several methods for comparison

4 小波包能量谱分析

以原始爆破信号及经降噪处理后的纯净信号为研究对象进行小波包分析，选用“db8”小波基函数分别进行8 层分解，并利用MATLAB对式(17)～(20)进行编程，得到各个子频带的能量及信号总能量，并求得各频带能量百分比如图11 所示。

由图11 中可知，与原始信号相比，CEEMDAN-小波包降噪得到的纯净信号中低频带能量百分比有较大幅度的增加；同时，降噪后信号的高频带能量百分比呈减小的趋势。研究发现[27]，隧道爆破信号的特征信息主要集中在中低频率带。因此，可以得到：经过CEEMDAN-小波包降噪处理的隧道爆破信号在去除噪声分量的同时，信号能量分布向中低频带转移，很好地保留了原始爆破信号中的特征信息，为后续对爆破信号进行深入分析奠定了基础。