基于變分模態(tài)分解的延時識別法在短延時微差爆破中的應(yīng)用

賈貝， 凌天龍， 侯仕軍， 劉殿書

(中國礦業(yè)大學(xué)(北京)，力學(xué)與建筑工程學(xué)院，北京 100083)

微差爆破是降低爆破震動效應(yīng)、減小炸藥單耗的一種有效技術(shù)手段，近年來，該技術(shù)在工程爆破實踐中得到了廣泛的應(yīng)用，合理的微差延遲時間是保證良好爆破效果的關(guān)鍵[1-4].由于雷管的延期時間偏差和爆破網(wǎng)路的延時性，使得實際微差延時與設(shè)計值具有一定差距，進而影響工程技術(shù)人員準確掌握微差時間與爆破效果之間的關(guān)系.工程中為了避免串段現(xiàn)象通常采用跳段起爆，其微差延時大于50 ms，但部分工程因降振需要采用多種50 ms延時以下的雷管，而目前對于短延時識別的研究較少.因此，研究短延時下的微差識別具有重要的現(xiàn)實意義.

21世紀初，小波變換因具有突出分析能量突變的特點，首次被引入到爆破信號的微差識別中[5-6].在實際應(yīng)用中，研究者發(fā)現(xiàn)以小波變換為基礎(chǔ)的微差延遲識別方法都存在最優(yōu)小波基的選擇問題，由于傳統(tǒng)小波基的波形固定不能完全匹配爆破振動信號，會導(dǎo)致微差識別的精確性、分辨率和穩(wěn)定性下降.因此，凌同華、張勝等[7-8]將模式自適應(yīng)小波和雙正交小波引入到爆破微差識別中，并證明兩種小波基對微差識別精確度和分別率的提升能力，但微差識別的計算量和耗時也隨之增加.

因為經(jīng)驗?zāi)B(tài)分解算法能突出數(shù)據(jù)局部特征且具有自適應(yīng)性，張義平等[9]提出了基于EMD的微差延遲時間識別方法.龔敏等[10]經(jīng)過現(xiàn)場試驗證明基于HHT的EMD識別法在微差延遲時間識別方面具有比小波變換識別法更高的精確性與分辨率.但隨著EMD識別法的廣泛應(yīng)用，發(fā)現(xiàn)EMD法仍存在一定的缺陷：其通過包絡(luò)線分解的結(jié)束判斷標準沒有科學(xué)依據(jù)，可能使IMF分量出現(xiàn)模態(tài)混疊問題，從而導(dǎo)致識別結(jié)果不準確；而且需要經(jīng)過多次迭代才可得到IMF分量，耗時長、計算量大.

2014年Dragomiretskiy等[11]提出了一種自適應(yīng)信號處理方法:變分模態(tài)分解(variational mode decomposition，VMD)，該方法通過迭代搜尋變分模型最優(yōu)解來確定每個分量的頻率中心及帶寬，從而能夠自適應(yīng)地實現(xiàn)信號的頻域剖分及各分量的有效分離.與EMD相比，VMD將信號分解轉(zhuǎn)化為變分問題求解過程，具有堅實的理論基礎(chǔ).本文結(jié)合VMD的分解特性將其引入爆破信號處理中，提出一種新的微差爆破延遲識別方法:VMD識別法.以八達嶺長城站工程中的微差爆破振動信號為例，分別應(yīng)用EMD法、HHT法和VMD法識別實際微差延時，對比分析VMD法在小延時微差識別的準確性和穩(wěn)定性.

1 變分模態(tài)分解原理

VMD算法把信號分解為多個本征模態(tài)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分量時采用循環(huán)篩分剝離的方式，分解非平穩(wěn)隨機信號過程中時常會出現(xiàn)模態(tài)混疊等缺陷.與EMD算法不同，VMD算法將信號分解過程轉(zhuǎn)化為變分求解過程，即把分解問題轉(zhuǎn)移到變分框架內(nèi)處理，通過尋找變分模型的最優(yōu)解獲取IMF分量，算法核心包括變分問題的構(gòu)造和變分問題的求解.

1.1 變分問題的構(gòu)造

假設(shè)每個模態(tài)分量都緊湊的圍繞于一個中心頻率分布，且具有有限帶寬，中心頻率會隨著分解變化而變化.VMD算法中變分問題的核心為：以輸入信號f(t)等于IMF分量之和為前提，尋找最小的IMF分量的預(yù)估帶寬之和，構(gòu)造過程如下.

① 對于每個IMF分量uk(t)，利用希爾伯特變換構(gòu)造解析信號后，通過混合指數(shù)調(diào)諧各自估計中心頻率的方法，將每個IMF分量的頻譜調(diào)制到相應(yīng)的基頻帶上

(2)

式中：uk={u1,u2,…,uk}為分解得到的k個IMF分量；wk={w1,w2,…,wk}為各IMF分量的中心頻率；j表示虛數(shù)；δ(t)為狄拉克函數(shù).

② 通過解調(diào)信號的高斯平滑度，即計算式(2)表示的信號梯度的平方L2范數(shù)，估計出各IMF分量的帶寬，構(gòu)造的變分問題如下所示：

1.2 變分問題的求解

為求取式(3)中的約束變分問題，引入懲罰因子α和Lagrange乘法算子λ(t)，其中懲罰因子α為較大的正數(shù)且在高斯噪聲存在的情況下可保證信號的重構(gòu)精度，Lagrange算子λ(t)使得約束條件保持嚴格性，構(gòu)造的擴展Lagrange表達式如下:

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

(4)

VMD中采用了乘法算子交替方向法(ADMM)解決以上變分模型，通過交替更新ukn+1、wkn+1、λkn+1尋求擴展的拉格朗日函數(shù)的鞍點，此點即為變分模型的最優(yōu)解.具體步驟如下：

② 執(zhí)行循環(huán)：n=n+1；

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

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

首先在頻域內(nèi)計算式(5)得出ukn+1對應(yīng)的頻域函數(shù)，而后對式(5)進行傅里葉逆變換，即可得到時域內(nèi)的IMF分量為

k∈{1,k}

(5)

(6)

④ 而后更新λ：

⑤ 對于給定判別精度e>0，若滿足式(8)條件，則停止迭代，否則返回步驟(2).

(8)

1.3 模態(tài)分量的瞬時時能識別

VMD將信號分解為多個IMF分量，每個IMF分量均為時間函數(shù)，對IMF分量c(t)進行Hilbert變換：

(9)

式中：H[c(t)]為Hilbert變換函數(shù)；pv為柯西主值.

構(gòu)造IMF信號的解析信號z(t)，即

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

(10)

式中a(t)為幅值函數(shù)

(11)

φ(t)為相位函數(shù)

(12)

對IMF分量c(t)進行Hilbert變換后，可將原始IMF分量表示為：

(13)

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

由式(13)可知，IMF分量可表達為時間、瞬時頻率f(t)的函數(shù)，即

(14)

式中H(f,t)為Hilbert譜.

對頻率積分，可定義Hilbert瞬時能量為

(15)

2 基于變分模態(tài)分解的微差爆破延時識別

2.1 VMD分解的參數(shù)確定

根據(jù)變分模態(tài)原理，處理信號前需要事先設(shè)置輸入?yún)?shù)：懲罰因子α與分解層數(shù)k.其中懲罰因子α的取值會影響分解精度，其取值越低各個IMF分量的帶寬越大，而取值越高IMF分量的帶寬越小，甚至?xí)沟贸绦蜻M入死循環(huán)，無法得到結(jié)果[11-12].多數(shù)情況下，研究者[13,14-18]將懲罰因子α設(shè)置為2倍的輸入信號長度即可滿足分解需求.

與懲罰因子α相比，分解層數(shù)k的取值更為重要，因為其取值直接影響分解結(jié)果是否正確.現(xiàn)階段針對VMD分解中變量k的研究還處于起步階段，且沒有形成統(tǒng)一的求解方法.研究者[12]認為當分解層數(shù)k為最佳分解層數(shù)時，各個IMF分量之間的頻譜分布均勻且連續(xù).根據(jù)此特點，把各IMF分量瞬時頻率均值特征曲線的斜率作為判斷依據(jù)，當曲線的左端出現(xiàn)水平或彎曲現(xiàn)象時，認為信號被過分解，若此時k=n，那么k=n-1為信號的最佳分解層數(shù).

2.2 VMD法的IMF分量選取

信號被分解為K個IMF分量后，需要從中選出部分IMF分量進行瞬時能量別.圖2為典型的電子雷管爆破振動信號，爆破參數(shù)如表1所示.對信號依次進行VMD分解，繪制不同K值對應(yīng)的IMF分量瞬時頻率均值特征曲線，圖1為K=15～20的特征曲線.觀察圖1發(fā)現(xiàn)，當進行K≥19的VMD分解時，特征曲線的左端出現(xiàn)了水平或下彎現(xiàn)象，此時信號被過度分解，因此確定K=18為該信號的最佳分解層數(shù)，且各階IMF分量的中心頻率如表2所示.

表1 典型電子雷管爆破振動信號參數(shù)

圖1 K=15～20的模態(tài)分量瞬時頻率均值的特征曲線

表2 各階IMF分量中心頻率

圖2 典型電子雷管爆破振動信號

圖3為信號的頻譜圖，信號的振動頻率主要分布于90～350 Hz之間，表明此信號的振動集中發(fā)生于90～350 Hz之間，因此選取中心頻率處于此頻率范圍的IMF分量進行瞬時能量識別，此方法在還原爆破信號攜帶信息的同時，排除了噪聲等不利因素的影響，使得識別結(jié)果更加精確.選取1階和2階IMF分量求和，并進行VMD法識別，結(jié)果如圖4(a)所示，可以從中分辨出23個明顯的波峰，與實際爆破雷管數(shù)相同.圖4(b)和圖4(c)分別為1階和4階IMF分量的識別結(jié)果，都無法精確識別所有雷管信息.因此，選擇中心頻率處于信號主頻范圍的IMF分量作為識別分量.

圖3 信號的頻譜圖

圖4 VMD法識別結(jié)果

3 工程實例分析

八達嶺長城站位于八達嶺滾天溝停車場下方新八達嶺隧道內(nèi)，是國內(nèi)首座采用礦山法施工的深埋高鐵地下車站.車站主體長度450 m，中心處線路埋深約102.55 m，分三連拱區(qū)段和三洞分離標準段，其中三洞分離段398 m，三洞分離段通過保留巖體分割，巖體最小厚度2.23 m，最大厚度6 m，需嚴格控制爆破振動.因此采用高精度電子雷管控制振動，根據(jù)中夾巖層厚度改變炮孔間的延期時間以達到控制振動的目的.

2017年11月21日的爆破施工中設(shè)置的電子雷管延遲時間分別為8，15，20 ms，詳細的測點布置和爆破參數(shù)分別如圖5、表3所示.以此次爆破為例，分別利用EMD，HHT和VMD識別方法處理測點的數(shù)據(jù)，得到小延時下三種方法的結(jié)果，其中171測點識別結(jié)果如圖6所示.

圖5 測點布置圖

表3 2017年11月21日三洞分離爆破參數(shù)

圖6～8包含三種方法對171測點數(shù)據(jù)的識別結(jié)果，并將識別結(jié)果按照炮孔編號分為1#～6#.由圖6可知，HHT法在1#(延期15 ms)、2#(20 ms)、3#(20 ms)和4#(20 ms)識別結(jié)果的波峰雜亂，且明顯多于實際炮孔數(shù)，無法精確找到起爆時間；在5#(8 ms)可識別出15個波峰多于實際起爆的11個雷管；在6#(8 ms)共起爆51個電子雷管，識別出40個.

圖7中，EMD分別識別出14、14、15、8、2和9個電子雷管，相對于實際起爆數(shù)24、18、20、15、11和50而言，識別出的雷管信息較少.圖8中，VMD分別識別出20、17、20、13、3和23個電子雷管，雖然仍未識別出全部雷管信息，但與實際起爆數(shù)相差較小.

圖6 HHT法對171#測點電子雷管微差識別結(jié)果

圖7 EMD法對171#測點電子雷管微差識別結(jié)果

圖8 VMD法對171#測點電子雷管微差識別結(jié)果

圖9為23日EMD法和VMD法的電子雷管識別率對比圖.

由圖9(a)可發(fā)現(xiàn)，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法的識別能力更強.

為驗證結(jié)論的普適性，隨機選取某日三洞分離段爆破測量數(shù)據(jù)進行實驗，處理結(jié)果如圖10所示.

對比圖10(a)和圖10(b)可得，在8 ms和20 ms的電子雷管識別中，VMD法的識別能力強于EMD法；在25 ms的電子雷管識別中，VMD法與EMD法識別能力基本相同，但總體而言，VMD法對小延時雷管的識別能力強于EMD法.

圖9 23日電子雷管識別率對比圖

圖10 某日電子雷管識別率對比圖

4 結(jié) 論

針對變分模態(tài)分解在理論與處理低頻信號的優(yōu)勢，結(jié)合八達嶺長城站現(xiàn)場監(jiān)測數(shù)據(jù)，得出以下結(jié)論：

① 基于變分模態(tài)分解原理，采用IMF分量的瞬時頻率均值特征曲線作為判斷標準，求解最佳分解層數(shù)K值，選取所有中心頻率處于主頻分布范圍的IMF分量進行識別，其識別結(jié)果為對應(yīng)的雷管起爆時間.

② 在短延時爆破識別中，HHT法識別結(jié)果的波峰雜亂且明顯多于實際起爆炮孔數(shù)，無法精確找到起爆時間；EMD法和VMD法的識別率波動較大處于10%～90%之間，且隨著雷管延時的增加而增加.

③ 當延時低于20 ms時，EMD法的識別率約為VMD法的3/4；當延時高于20 ms時，兩種方法識別率都在80%以上.總體而言，VMD法對短延時雷管的識別能力強于EMD法.