基于变分模态分解的延时识别法在短延时微差爆破中的应用

贾贝， 凌天龙， 侯仕军， 刘殿书

(中国矿业大学(北京)，力学与建筑工程学院，北京 100083)

微差爆破是降低爆破震动效应、减小炸药单耗的一种有效技术手段，近年来，该技术在工程爆破实践中得到了广泛的应用，合理的微差延迟时间是保证良好爆破效果的关键[1-4].由于雷管的延期时间偏差和爆破网路的延时性，使得实际微差延时与设计值具有一定差距，进而影响工程技术人员准确掌握微差时间与爆破效果之间的关系.工程中为了避免串段现象通常采用跳段起爆，其微差延时大于50 ms，但部分工程因降振需要采用多种50 ms延时以下的雷管，而目前对于短延时识别的研究较少.因此，研究短延时下的微差识别具有重要的现实意义.

21世纪初，小波变换因具有突出分析能量突变的特点，首次被引入到爆破信号的微差识别中[5-6].在实际应用中，研究者发现以小波变换为基础的微差延迟识别方法都存在最优小波基的选择问题，由于传统小波基的波形固定不能完全匹配爆破振动信号，会导致微差识别的精确性、分辨率和稳定性下降.因此，凌同华、张胜等[7-8]将模式自适应小波和双正交小波引入到爆破微差识别中，并证明两种小波基对微差识别精确度和分别率的提升能力，但微差识别的计算量和耗时也随之增加.

因为经验模态分解算法能突出数据局部特征且具有自适应性，张义平等[9]提出了基于EMD的微差延迟时间识别方法.龚敏等[10]经过现场试验证明基于HHT的EMD识别法在微差延迟时间识别方面具有比小波变换识别法更高的精确性与分辨率.但随着EMD识别法的广泛应用，发现EMD法仍存在一定的缺陷：其通过包络线分解的结束判断标准没有科学依据，可能使IMF分量出现模态混叠问题，从而导致识别结果不准确；而且需要经过多次迭代才可得到IMF分量，耗时长、计算量大.

2014年Dragomiretskiy等[11]提出了一种自适应信号处理方法:变分模态分解(variational mode decomposition，VMD)，该方法通过迭代搜寻变分模型最优解来确定每个分量的频率中心及带宽，从而能够自适应地实现信号的频域剖分及各分量的有效分离.与EMD相比，VMD将信号分解转化为变分问题求解过程，具有坚实的理论基础.本文结合VMD的分解特性将其引入爆破信号处理中，提出一种新的微差爆破延迟识别方法:VMD识别法.以八达岭长城站工程中的微差爆破振动信号为例，分别应用EMD法、HHT法和VMD法识别实际微差延时，对比分析VMD法在小延时微差识别的准确性和稳定性.

1 变分模态分解原理

VMD算法把信号分解为多个本征模态IMF(intrinsic mode functions)分量，且将IMF分量重新定义为

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

(1)

式中：t为时间；uk(t)为各IMF分量；Ak(t)为瞬时幅值，且Ak(t)≥0；φk(t)为瞬时相位，且φ′k(t)≥0.

EMD算法获取IMF分量时采用循环筛分剥离的方式，分解非平稳随机信号过程中时常会出现模态混叠等缺陷.与EMD算法不同，VMD算法将信号分解过程转化为变分求解过程，即把分解问题转移到变分框架内处理，通过寻找变分模型的最优解获取IMF分量，算法核心包括变分问题的构造和变分问题的求解.

1.1 变分问题的构造

假设每个模态分量都紧凑的围绕于一个中心频率分布，且具有有限带宽，中心频率会随着分解变化而变化.VMD算法中变分问题的核心为：以输入信号f(t)等于IMF分量之和为前提，寻找最小的IMF分量的预估带宽之和，构造过程如下.

① 对于每个IMF分量uk(t)，利用希尔伯特变换构造解析信号后，通过混合指数调谐各自估计中心频率的方法，将每个IMF分量的频谱调制到相应的基频带上

(2)

式中：uk={u1,u2,…,uk}为分解得到的k个IMF分量；wk={w1,w2,…,wk}为各IMF分量的中心频率；j表示虚数；δ(t)为狄拉克函数.

② 通过解调信号的高斯平滑度，即计算式(2)表示的信号梯度的平方L2范数，估计出各IMF分量的带宽，构造的变分问题如下所示：

1.2 变分问题的求解

为求取式(3)中的约束变分问题，引入惩罚因子α和Lagrange乘法算子λ(t)，其中惩罚因子α为较大的正数且在高斯噪声存在的情况下可保证信号的重构精度，Lagrange算子λ(t)使得约束条件保持严格性，构造的扩展Lagrange表达式如下:

k({uk},{wk},λ}):=

(4)

VMD中采用了乘法算子交替方向法(ADMM)解决以上变分模型，通过交替更新ukn+1、wkn+1、λkn+1寻求扩展的拉格朗日函数的鞍点，此点即为变分模型的最优解.具体步骤如下：

② 执行循环：n=n+1；

③ 对所有w>0的分量,更新uk,wk;

ukn+1的更新求解过程如下.

首先在频域内计算式(5)得出ukn+1对应的频域函数，而后对式(5)进行傅里叶逆变换，即可得到时域内的IMF分量为

k∈{1,k}

(5)

(6)

④ 而后更新λ：

⑤ 对于给定判别精度e>0，若满足式(8)条件，则停止迭代，否则返回步骤(2).

(8)

1.3 模态分量的瞬时时能识别

VMD将信号分解为多个IMF分量，每个IMF分量均为时间函数，对IMF分量c(t)进行Hilbert变换：

(9)

式中：H[c(t)]为Hilbert变换函数；pv为柯西主值.

构造IMF信号的解析信号z(t)，即

z(t)=c(t)+jH[c(t)]=a(t)ejφ(t)

(10)

式中a(t)为幅值函数

(11)

φ(t)为相位函数

(12)

对IMF分量c(t)进行Hilbert变换后，可将原始IMF分量表示为：

(13)

式中：Re为取实部，且dφ(t)=2πf(t)dt.

由式(13)可知，IMF分量可表达为时间、瞬时频率f(t)的函数，即

(14)

式中H(f,t)为Hilbert谱.

对频率积分，可定义Hilbert瞬时能量为

(15)

2 基于变分模态分解的微差爆破延时识别

2.1 VMD分解的参数确定

根据变分模态原理，处理信号前需要事先设置输入参数：惩罚因子α与分解层数k.其中惩罚因子α的取值会影响分解精度，其取值越低各个IMF分量的带宽越大，而取值越高IMF分量的带宽越小，甚至会使得程序进入死循环，无法得到结果[11-12].多数情况下，研究者[13,14-18]将惩罚因子α设置为2倍的输入信号长度即可满足分解需求.

与惩罚因子α相比，分解层数k的取值更为重要，因为其取值直接影响分解结果是否正确.现阶段针对VMD分解中变量k的研究还处于起步阶段，且没有形成统一的求解方法.研究者[12]认为当分解层数k为最佳分解层数时，各个IMF分量之间的频谱分布均匀且连续.根据此特点，把各IMF分量瞬时频率均值特征曲线的斜率作为判断依据，当曲线的左端出现水平或弯曲现象时，认为信号被过分解，若此时k=n，那么k=n-1为信号的最佳分解层数.

2.2 VMD法的IMF分量选取

信号被分解为K个IMF分量后，需要从中选出部分IMF分量进行瞬时能量别.图2为典型的电子雷管爆破振动信号，爆破参数如表1所示.对信号依次进行VMD分解，绘制不同K值对应的IMF分量瞬时频率均值特征曲线，图1为K=15～20的特征曲线.观察图1发现，当进行K≥19的VMD分解时，特征曲线的左端出现了水平或下弯现象，此时信号被过度分解，因此确定K=18为该信号的最佳分解层数，且各阶IMF分量的中心频率如表2所示.

表1 典型电子雷管爆破振动信号参数

图1 K=15～20的模态分量瞬时频率均值的特征曲线

表2 各阶IMF分量中心频率

图2 典型电子雷管爆破振动信号

图3为信号的频谱图，信号的振动频率主要分布于90～350 Hz之间，表明此信号的振动集中发生于90～350 Hz之间，因此选取中心频率处于此频率范围的IMF分量进行瞬时能量识别，此方法在还原爆破信号携带信息的同时，排除了噪声等不利因素的影响，使得识别结果更加精确.选取1阶和2阶IMF分量求和，并进行VMD法识别，结果如图4(a)所示，可以从中分辨出23个明显的波峰，与实际爆破雷管数相同.图4(b)和图4(c)分别为1阶和4阶IMF分量的识别结果，都无法精确识别所有雷管信息.因此，选择中心频率处于信号主频范围的IMF分量作为识别分量.

图3 信号的频谱图

图4 VMD法识别结果

3 工程实例分析

八达岭长城站位于八达岭滚天沟停车场下方新八达岭隧道内，是国内首座采用矿山法施工的深埋高铁地下车站.车站主体长度450 m，中心处线路埋深约102.55 m，分三连拱区段和三洞分离标准段，其中三洞分离段398 m，三洞分离段通过保留岩体分割，岩体最小厚度2.23 m，最大厚度6 m，需严格控制爆破振动.因此采用高精度电子雷管控制振动，根据中夹岩层厚度改变炮孔间的延期时间以达到控制振动的目的.

2017年11月21日的爆破施工中设置的电子雷管延迟时间分别为8，15，20 ms，详细的测点布置和爆破参数分别如图5、表3所示.以此次爆破为例，分别利用EMD，HHT和VMD识别方法处理测点的数据，得到小延时下三种方法的结果，其中171测点识别结果如图6所示.

图5 测点布置图

表3 2017年11月21日三洞分离爆破参数

图6～8包含三种方法对171测点数据的识别结果，并将识别结果按照炮孔编号分为1#～6#.由图6可知，HHT法在1#(延期15 ms)、2#(20 ms)、3#(20 ms)和4#(20 ms)识别结果的波峰杂乱，且明显多于实际炮孔数，无法精确找到起爆时间；在5#(8 ms)可识别出15个波峰多于实际起爆的11个雷管；在6#(8 ms)共起爆51个电子雷管，识别出40个.

图7中，EMD分别识别出14、14、15、8、2和9个电子雷管，相对于实际起爆数24、18、20、15、11和50而言，识别出的雷管信息较少.图8中，VMD分别识别出20、17、20、13、3和23个电子雷管，虽然仍未识别出全部雷管信息，但与实际起爆数相差较小.

图6 HHT法对171#测点电子雷管微差识别结果

图7 EMD法对171#测点电子雷管微差识别结果

图8 VMD法对171#测点电子雷管微差识别结果

图9为23日EMD法和VMD法的电子雷管识别率对比图.

由图9(a)可发现，EMD法对延期时间为8，15和20 ms的电子雷管识别率分别为6%～24%、33%～58%以及47%～70%，其识别能力随着延期时间的增大而增强，但总体识别率高并不能满足实际需求；图9(b)中，VMD法对延期时间为8，15和20 ms的电子雷管识别率分别为20%～60%、84%～94%以及82%～95%，其识别能力在8 ms时较差，但在15 ms和20 ms较强.对比可看出，在小延时雷管的识别中，VMD法的识别能力更强.

为验证结论的普适性，随机选取某日三洞分离段爆破测量数据进行实验，处理结果如图10所示.

对比图10(a)和图10(b)可得，在8 ms和20 ms的电子雷管识别中，VMD法的识别能力强于EMD法；在25 ms的电子雷管识别中，VMD法与EMD法识别能力基本相同，但总体而言，VMD法对小延时雷管的识别能力强于EMD法.

图9 23日电子雷管识别率对比图

图10 某日电子雷管识别率对比图

4 结 论

针对变分模态分解在理论与处理低频信号的优势，结合八达岭长城站现场监测数据，得出以下结论：

① 基于变分模态分解原理，采用IMF分量的瞬时频率均值特征曲线作为判断标准，求解最佳分解层数K值，选取所有中心频率处于主频分布范围的IMF分量进行识别，其识别结果为对应的雷管起爆时间.

② 在短延时爆破识别中，HHT法识别结果的波峰杂乱且明显多于实际起爆炮孔数，无法精确找到起爆时间；EMD法和VMD法的识别率波动较大处于10%～90%之间，且随着雷管延时的增加而增加.

③ 当延时低于20 ms时，EMD法的识别率约为VMD法的3/4；当延时高于20 ms时，两种方法识别率都在80%以上.总体而言，VMD法对短延时雷管的识别能力强于EMD法.