基于自适应VMD-MPE算法的矿山爆破地震波信号降噪方法研究

彭亚雄， 刘广进， 苏 莹， 陈春晖， 刘运思

(1. 湖南科技大学 岩土工程稳定控制与健康监测湖南省重点实验室， 湖南 湘潭 411201；2. 中国地质大学(武汉) 工程学院， 武汉 430074)

爆破工程中采用实测信号对爆破有害效应进行详细分析与评价。然而受矿山环境复杂性、监测设备电磁干扰、岩体介质的反射及折射等原因的影响，大量高频噪声导致信号波形畸变，掩盖了真实信号成分[1-2]。为了掌握爆破地震波波形、强度和频率特征，必须通过信号处理技术对实测信号进行降噪处理。

爆破地震波信号属于典型的非平稳随机信号，目前针对这类信号的降噪算法主要包括小波类算法[3-4]和经验模态分解(empirical mode decomposition, EMD)算法[5-6]。在爆破信号降噪研究领域，熊正明等[7]通过多次使用平移不变小波对爆破信号进行平移并进行阈值降噪，实现了摆脱小波基的信号降噪处理。路亮等[8]提出了利用提升小波包变换的最优基搜索改进算法，有效去除爆破振动信号噪声。黄智刚等[9]提出了基于CEEMDAN-MPE算法的隧道爆破地震波信号降噪方法，不仅能够去除高频噪声，对信号所含主要信息影响较小。Peng等[10]针对水下钻孔爆破振动信号所含噪声特点，提出了CEEMDAN光滑降噪模型，有效降低了实测信号噪声含量。这些降噪算法的提出与工程应用，充分说明了信号处理技术对实测爆破信号降噪处理的有效性。然而，小波算法在降噪过程中小波基函数和分解层次难以确定，使得这类方法的自适应性不强，降噪效果难以保证；EMD及其改算法具有较好的自适性，但其端点效应和模态混叠问题无法有效解决，使得这两类算法在处理爆破地震波信号的鲁棒性较低。

变分模态分解(variational mode decomposition, VMD)是一种非递归的信号分解方法，迭代搜寻变分的最优解，以确定信号各部分频率中心及带宽，完成频域划分和分量分离。该方法避免了端点效应和模态混叠问题[11]，对各类信号具有普遍适应性和鲁棒性，但是存在模态数K难以确定的问题。因此，本文提出了基于能量差参数ξ的自适应VMD分解方法，将分解得到的本征模态函数(intrinsic mode function, IMF)进行多尺度排列熵(multi-scale permutation entropy, MPE)随机性检测，去除噪声IMF达到降噪目的，并应用于实测矿山爆破地震波信号降噪处理。

1 VMD分解原理

VMD算法是在信号频域内利用迭代计算获得变分模态最优解，变换模态函数和中心频率，得到不同带宽的本征模态函数。将信号模态估计变为变分问题进行求解，基本步骤如下[12]。

将原信号x(t)分解为K个中心频率为ωk的模态函数uk，其中K为模态数，需人为设定[13]。模态函数表示为

uk(t)=Ak(t)cos(φk(t))

(1)

式中:uk(t)为第k个分量；Ak(t)和φk(t)分别为瞬时幅值和相位。

对uk(t)进行Hibert变换，得到解析信号和单边频谱ψ；解析信号加入估计中心频率，将模态频谱转换至基频带；计算解析信号梯度的平方L2范数，估计模态信号带宽，成为约束性变分问题，可以求解出x(t)的IMF

(2)

式中:{uk}为信号分解得到的K个IMF分量；{ωk}为各分量对应的中心频率；δ(t)为脉冲函数。

(3)

2 自适应VMD-MPE降噪方法

2.1 自适应VMD算法

VMD算法中必须人为设定模态数K(即IMF分量个数)和惩罚因子α。当模态数K过少，将导致原信号缺失；过多则会产生频率混叠、过分解等问题。当惩罚因子α较小时，IMF带宽较大易产生混叠现象，当α较大则影响收敛速度[14]。

通过观察不同模态数K分解得到的各IMF中心频率，可以选出最合适的模态数K。然而该方法在判断中心频率变化规律时，受到人为因素影响较大[15]。为此，本文提出了基于信号能量的自适应VMD算法。爆破地震波信号能量定义[16]为

(6)

式中:E为地震波信号能量;Es为原信号能量；第k个IMF能量表示为Ek；v(t)为信号序列；t为信号采集长度。

不同模态数K的IMF能量和与原信号能量不同，为了确定最适合模态数K，定义了重构信号与原信号能量差参数ξ的计算式

(7)

能量差参数ξ随着模态数K增加而增加；信号过分解前，ξ的增长幅度很小，而当信号出现分解时，ξ将大幅增长。因此，由K=2,3,…进行信号分解，计算能量差参数ξK；根据文献[17-18]，当|ξK+1-ξK|>0.1时，认为模态数K+1分解的信号处于过分解状态，则确定最合适模态数为K。

在模态数K值确定后，分别计算不同惩罚因子α条件下信号VMD分解的信噪比SNR

(8)

2.2 MPE算法

MPE算法是基于排列熵(permutation entropy, PE)的改进算法，可用于检测信号的随机和动力突变特性。将时间序列进行多尺度粗粒化，再计算PE。具体步骤如下[19]

① 时间序列x={x1,x2,…,xL}的多尺度粗粒化

(9)

(10)

(13)

2.3 自适应VMD-MPE降噪

自适应VMD-MPE降噪方法是对原信号进行自适应VMD分解，计算最适合的惩罚因子α和模态数K，进而获得K个IMF分量，其中IMF1～IMFK的主频逐渐增加。然后对各IMF分量进行MPE随机性检测，利用IMF分量的MPE均值判断是否为噪声成分。对信号主成分IMF分量进行重构，得到降噪信号。对于爆破地震波信号，真实信号成分的MPE阈值为0.6[20]。自适应VMD-MPE降噪流程如图1所示。

图1 自适应VMD-MPE降噪流程Fig.1 Denoising process with adaptive VMD-MPE

3 工程应用

3.1 信号降噪处理

邦山脉铁矿位于利比里亚中部，邦州西南地区，年开采量约为1 000万吨。矿区内岩层整体向北倾斜，倾角较陡65°～85°，主要岩性为含铁丰富片麻岩、片岩和石英岩，岩层完整，弱风化和中风化，裂隙不发育。矿山开采采用中深孔台阶爆破方式，爆破参数如表1所示。

表1 台阶爆破参数Tab.1 Bench blasting parameters

矿区范围内存在多个村庄，为评价矿山爆破对周围环境的影响，采用M20型测振仪进行爆破地震波监测。选取一条典型爆破地震波信号为研究对象(图2)。该实测信号的采样频率为4 000 sps，根据Nyquist采样定理，实测信号频率为2 000 Hz；采用时间为2 s，共采集8 000个采样点。

图2 爆破地震波信号(s1)Fig.2 Blasting seismic wave signal (s1)

自适应VMD算法的主要设定参数包括：模态数K、惩罚因子α、保真度τ和收敛条件，针对不同类型信号需要设置合适参数值。对于随机非平稳信号，保真度τ和收敛条件ε采用推荐值[21]；计算不同模态数K分解的能量差参数ξ和|ξK+1-ξK|如表2所示。根据自适应VMD算法的判断依据，该爆破地震波信号的最适合模态数K为6，其IMF分量如图3所示。然后通过计算不同惩罚因子α条件下信号VMD分解的信噪比SNR，得到适合本工程α=2 000，减少了分解过程的细节特征损失。

表2 不同模态数K分解的能量差参数ξTab.2 Energy difference parameter ξ with different modes number K

由表2和图3可知，采用VMD算法分解得到实测信号共6个IMF分量，IMF1～IMF6的中心频率逐渐增加。由IMF分量的波形和中心频率变化，可以推断IMF4～IMF6的幅值较小且频率较高可能是爆破地震波信号中的高频噪声分量，而IMF1～IMF3则为真实信号成分。

为准确区分真实信号成分和噪声，计算各IMF分量的MPE值。经过多次试算确定合适的嵌入维数m=6，时间延迟τ=1，尺度因子s=5。得到各IMF分量的MPE平均值如表3所示。

表3 IMF分量的MPE均值Tab.3 Mean MPE of IMF

由表3可知，IMF1～IMF6的MPE均值逐渐增加，IMF分量中所含噪声成分逐渐减少，说明自适应VMD算法能够很好地将真实信号与噪声信号成分分离。对于爆破地震波信号，真实信号成分的MPE阈值为0.6，IMF4～IMF6的MPE均值大于阈值为噪声信号，与上述波形分析结果一致。因此，将IMF4～IMF6从原信号中剔除，保留爆破地震波真实信号成分，达到降噪的目的。同时，自适应算法避免了VMD分解过程中的人工干预，降低了人为误差。

降噪后的爆破地震波信号如图4所示。与原信号波形(图2)相比，二者的波形一致，噪声得到明显消除，且对爆破地震波峰值振速影响较小。

图4 降噪后的爆破地震波信号Fig.4 Blasting seismic wave signal after denoising

3.2 时频能量分析

自适应最优核(adaptive optimal kernel, AOK)时频技术采用了短时模糊函数和时变自适应最优核函数，能够获得更多的细节信息，较好反映了爆破地震波信号的时频特性。为了进一步验证自适应VMD-MPE算法的降噪效果，对原信号和降噪信号进行AOK时频谱对比，如图5所示。图中X为峰值能量，Y为信号主频。

(a) 原信号

(b) 降噪信号图5 AOK时频谱Fig.5 AOK time-frequency spectrum

由图5可知，实测矿山爆破地震波信号的时频能量主要分布在频率0～250 Hz和时间0.2～0.9 s范围内，信号主频为10.71 Hz。原信号在频率500 Hz以上含有明显的噪声成分，是导致信号失真的主要原因，如图5(a)中虚线框中所示。利用自适应VMD-MPE算法获得的降噪信号，较好地剔除了原有噪声成分；降噪信号与原信号的主频一致，峰值能量和主频带能量(0～250 Hz)均没有明显降低。说明自适应VMD-MPE算法不仅能有效剔除实测矿山爆破地震波信号中的高频噪声，且对低频信号能量影响较小。

4 降噪效果对比分析

4.1 评价指标

实测信号在VMD分解过程中不可避免导致原有信号成分的丢失，为反映原信号与重构信号的差异性，定义了重构标准差ESD作为评价信号丢失和失真情况的评价指标。采用降噪信号与原信号的信噪比SNR、均方根误差ε作为降噪效果评价指标，分析实测矿山爆破地震波信号的降噪效果。重构标准差ESD和均方根误差ε定义如下

(1) 重构标准差ESD

(14)

(2) 均方根误差ε

(15)

此外，对比降噪前后信号波形特征，能够定性判断降噪效果，以确保特征波形的一致性和明显噪点去除干净。

4.2 降噪效果对比

为了进一步分析自适应VMD-MPE算法的降噪效果，对3组矿山爆破地震波信号(如图6)进行降噪处理，并与EEMD-MPE、CEEMDAN-MPE降噪算法进行对比。三种算法获得的降噪信号如图6所示，并分别计算重构标准差ESD、信噪比SNR、均方根误差ε，如表4所示。

图6 原信号与降噪信号对比Fig.6 Comparison between original signals and denoised signals

图6 原信号与降噪信号对比Fig.6 Comparison between original signals and denoised signals

表4 爆破地震波重构和降噪效果指标Tab.4 Reconstruct and denoised effect index of blasting seismic wave signals

对比3组实测矿山爆破地震波信号的降噪结果，由波形曲线(图6)可知，EEMD-MPE、CEEMDAN-MPE和自适应VMD-MPE三种算法均取得了一定程度降噪效果，去除了部分高频噪声，其中EEMD-MPE算法的降噪信号仍然具有较为明显的噪声，降噪效果相对较差。由表3可知，在信号分解与重构过程中，原信号与重构信号的标准差ESD在0.008 3～0.020 3，均属于较低范围，说明三种算法的保真性均较好。然而由于EEMD分解过程中，以增加白噪声方式消除模态混叠问题，导致了其重构标准差ESD均最大。对比两个降噪效果指标，自适应VMD-MPE算法的信噪比SNR均最大，CEEMDAN-MPE算法略小于EMD自适应VMD-MPE算法；自适应VMD-MPE算法的均方根误差ε最低。

因此，自适应VMD-MPE算法的降噪波形和降噪效果指标均优于其他两种算法，能够更好地去除高频噪声信号，同时能够更好地保留信号中的有效信息，在矿山爆破地震波信号降噪处理上具有明显的优越性。

5 结 论

受到矿山复杂环境的影响，矿山爆破地震波信号所含高频噪声较多，导致真实信号内容被噪声掩盖。针对这一问题提出了基于自适应VMD-MPE算法降噪方法，用于矿山爆破地震波信号降噪处理。

(1) 对爆破地震波信号进行VMD分解，得到不同频率的IMF，有效避免了模态混叠与端点效应问题；提出了基于能量差参数ξ确定模态数K的自适应方法，实现了最优分解，其结果比主观决策更为可靠。

(2) 对实测矿山爆破地震波信号进行处理表明，自适应VMD分解可以获得精细的IMF，MPE熵值识别出噪声IMF，达到了较好去除降噪的目的。对比分析降噪前后振动信号的AOK时频能量谱，自适应VMD-MPE降噪方法不仅能成功地去除500 Hz以上的高频噪声能量，而且对低频信号能量影响较小。

(3) 对比EEMD-MPE、CEEMDAN-MPE和自适应VMD-MPE三种降噪算法，三种方法都具有较好的保真性和一定的降噪效果。自适应VMD-MPE算法的降噪信号波形和降噪效果指标均优于EEMD-MPE、CEEMDAN-MPE算法，验证了基于自适应VMD-MPE算法的有效性，对矿山爆破地震波信号降噪和爆破有效效应分析具有指导意义。