基于經(jīng)驗?zāi)B(tài)分解的爆破延期識別優(yōu)化方法

易文華 ， 劉連生， 閆 雷， 董斌斌， 劉 偉， 楊 硯

(江西理工大學(xué) 資源與環(huán)境工程學(xué)院，江西 贛州 341000)

隨著爆破行業(yè)的發(fā)展，延期爆破作為降低爆破地震效應(yīng)的一種控制技術(shù)越來越多地應(yīng)用到工程實踐，對爆破振動信號進行延期分析是爆破參數(shù)設(shè)計優(yōu)化和爆破盲炮識別中常用的分析方法[1-2]；但在爆破工程作業(yè)中，由于受復(fù)雜地質(zhì)、地形條件、施工質(zhì)量等因素的影響，爆破信號中蘊含著豐富的振動狀態(tài)信息且彼此之間互相干擾，對信號延期識別精度造成了影響。

目前延期爆破延期識別的主要方法有八線示波器和金屬拉線法[3]、小波時-能密度[4-5]和經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition, EMD)識別等方法。但八線示波器和金屬拉線法操作復(fù)雜，成本較高，且精度較低。隨著振動信號分析理論在爆破工程中的發(fā)展，凌同華等利用小波變換的時-能密度法分析信號能量突變的優(yōu)勢，有效地識別出延期爆破延期時間，但小波變換分解的精度依賴小波基的選擇，選擇不同的小波基會產(chǎn)生不同精度的誤差。隨后張義平等[6]利用EMD法將爆破振動信號分解成若干個本征模函數(shù)(intrinsic mode function, IMF)分量,繼而通過Hilbert變換提取主IMF分量的包絡(luò)幅值對實際爆破延期時間進行識別。

由于EMD能自適應(yīng)地將信號按不同時間尺度進行分解，可以很好地提取非平穩(wěn)信號變化的特征[7]；與小波時-能密度法相比，EMD識別法不需要選擇基函數(shù)且自適應(yīng)性強。但EMD識別法在分解信號過程中，分解出的IMF分量之間會出現(xiàn)模態(tài)混疊現(xiàn)象[8-11]。針對這個問題，Wu等[8]提出了集總經(jīng)驗?zāi)B(tài)分解(ensemble empirical mode decomposition，EEMD)方法來抑制模態(tài)混疊現(xiàn)象，得到了較好的效果，但對低頻率比例的混合信號抑制效果不佳。司莉等[9]鏡像延拓方法以及高頻諧波法對EMD分解不完全現(xiàn)象進行了抑制，但選擇的頻率需要根據(jù)異常事件及基礎(chǔ)信號的頻率反復(fù)測試得出，具有較大的不穩(wěn)定性。詹瀛魚等[10]提出了基于解相關(guān)多頻率經(jīng)驗?zāi)B(tài)分解方法，通過添加掩蔽信號和嵌入相關(guān)系數(shù)處理， 抑制了模態(tài)混疊現(xiàn)象。李曉斌[11]采用了Ort判別法研究出對IMF分量進行完全正交化處理能夠有效地消除模態(tài)混疊現(xiàn)象。因此分解出的IMF分量之間是否具有混疊現(xiàn)象可由正交性來判斷。

針對上述問題的優(yōu)缺點，本文提出了基于主成分分析(principal component analysis，PCA)的完全正交經(jīng)驗?zāi)B(tài)分解(principal empirical mode decomposition, PEMD) 方法。由于模態(tài)混疊是由于EMD分解出的IMF分量之間的信息相互耦合，本質(zhì)為 IMF 分量之間不完全正交[12]，故可利用主成分分析[13-14]能將具有相關(guān)性的數(shù)組轉(zhuǎn)化為正交數(shù)組的特性，對IMF分量進行主成分分析，實現(xiàn)各IMF分量之間完全正交化，從而抑制模態(tài)混疊現(xiàn)象[15]，最終提升EMD方法對爆破振動信號延期識別的精度；并通過相似物理模擬試驗和德興露天臺階延期爆破試驗，與EMD方法進行對比，檢驗和評價PEMD方法的爆破延期識別效果。

1 基本原理簡介

經(jīng)驗?zāi)B(tài)分解是由Huang等[16]提出的一種信號時頻分析方法，該方法能夠自適應(yīng)地將原始信號x(t)分解成一系列IMF分量和一個殘余量。其中IMF分量突出了信號的局部特征，殘余分量則代表信號中的緩慢變化量。分解過程[17]可表達(dá)為

(1)

式中：ci(t)為第i個IMF分量；rn(t)為殘差，也稱信號的趨勢項。

主成分分析，是Hotelling[18]在1993年提出的一種多變量統(tǒng)計方法，通過正交變換將一組不完全正交的原始數(shù)據(jù)矩陣轉(zhuǎn)X換為完全正交矩陣V，在不損失原信號信息的前提下，達(dá)到正交化的目的。即

(2)

式中：Λ為非零對角元素組成的特征值矩陣；V特征向量組成的正交矩陣；n為采樣點的個數(shù)。

2 完全正交經(jīng)驗?zāi)B(tài)分解

模態(tài)混疊現(xiàn)象的產(chǎn)生與 EMD 不完全正交分解有關(guān)，分解出的各IMF分量之間信息相互耦合。因此需要采用一種處理方法使各IMF分量之間滿足完全正交性，就能夠?qū)δB(tài)混疊現(xiàn)象進行有效地抑制。本文提出一種基于PCA對EMD進行輔助分解的算法PEMD，通過對混疊的IMF分量進行主成分分析，實現(xiàn)各IMF分量之間完全正交化，從而提高EMD延期識別精度。該算法的實現(xiàn)流程如圖1所示。

圖1 PEMD算法流程圖Fig.1 PEMD algorithm flow chart

具體步驟如下:

步驟1將原始信號x(t)通過EMD分解成m個IMF分量，每個分量都取n個采樣點。

步驟2假設(shè)各個IMF分量分別為x1,x2,…,xm，第i個采樣點的第j個分量的取值為aij

步驟3計算IMF分量協(xié)方差矩陣R

(3)

R=(rij)m×n

(4)

式中，R為相關(guān)系數(shù)矩陣。

步驟4計算相關(guān)系數(shù)矩陣R的特征值λ和特征向量u，由特征向量組成m個正交IMF分量yi,i=1,2,…,m為

(5)

步驟5根據(jù)正交IMF分量的幅值和波形衰減特征選擇主分量[19-20]。

步驟6對主分量進行Hilbert變換，提取包絡(luò)線

(6)

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

(7)

式中：pv為柯西主值(cauchy principal value)；a(t)為解析信號z(t)的幅值,也稱為信號的包絡(luò)。

步驟7對主分量包絡(luò)線的峰值點進行識別，得到實際的延期時間。

通過將EMD的自適應(yīng)性和PCA的完全正交性特點融合起來，對EMD分解出的IMF分量進行完全正交化處理，可以有效地降低各IMF分量之間的相關(guān)性，從而達(dá)到抑制信號的模態(tài)混疊現(xiàn)象，最終提高振動信號EMD爆破延期識別精度。

3 相似物理模擬試驗

3.1 相似物理模型

由于露天邊坡現(xiàn)場地質(zhì)地形條件復(fù)雜，外界干擾因素對試驗研究影響較大，為了盡量控制其他外界的干擾因素，可根據(jù)爆破相似理論[21]，在滿足與露天邊坡現(xiàn)場的幾何、材料、動力以及邊界條件相似的情況下得到相似的結(jié)果，對新方法的初步研究具有重要意義。因此，為了構(gòu)建露天邊坡相似物理模型，結(jié)合露天邊坡幾何形狀將模型設(shè)計成四周高邊坡臺階、中間凹陷的形態(tài)，選用的幾何縮比為1 ∶100，同時為盡量滿足露天邊坡巖體材料特征要求，選用了等級為P·O 42.5的硅酸鹽水泥、細(xì)沙及水三種材料，選用的配比(水 ∶水泥 ∶砂子)為0.44 ∶1 ∶1.5，最后為使模型的動力以及邊界條件與露天邊坡相似，設(shè)計過程中盡可能加大了模型的尺寸長度以改善邊界效應(yīng)，設(shè)計模型的尺寸為270 cm×270 cm×110 cm，其剖面圖如圖2所示。

圖2 相似物理模型(mm)Fig.2 Similar physical model(mm)

其中，設(shè)計模型底座高38 cm，在底座之上共設(shè)計4個臺階面，各臺階面之間高程差均為24 cm，臺階總高72 cm，模型總高度110 cm，臺階坡面角均為67°，最終邊坡角42°。模型底座中心區(qū)域設(shè)計為炮孔布置區(qū)，以12.5 cm的孔間距及排間距共布置25個炮孔，炮孔深度13.5 cm。平臺包括最底層一共4個，由最底層到最上層臺階分別編號為0號、2號、4號和6號，以0號平臺為基準(zhǔn)，各臺階爆心距和高程差見表1。

表1 爆心距與高程差Tab.1 The difference between the burst distance and elevation cm

3.2 爆破延期識別

3.2.1 EMD識別

下面基于相似物理模型，利用本文提出的PEMD方法與EMD方法對一組典型的爆破振動信號進行延期識別對比。本次試驗采用隆芯一號系列數(shù)碼電子雷管四孔12 ms延期爆破起爆，在6號臺階處布置振動監(jiān)測點，試驗期間所有爆破均設(shè)置為單孔單發(fā)電雷管起爆。選擇一組典型爆破振動信號進行分析，其速度時程曲線如圖3所示。

圖3 速度時程曲線Fig.3 Velocity time-history curve

對其進行EMD分解，獲取了13個頻率從高到低依次排列的IMF分量，(篇幅有限，僅列前6個IMF分量)，對其進行EMD分解，獲取了13個頻率從高到低依次排列的IMF分量，并通過Hilbert變換做出各分量的包絡(luò)線，如圖4所示。

圖4 IMF分量及包絡(luò)線Fig.4 IMF components and envelope

由圖4可以看出，IMF1幅值較大，波形衰減明顯，攜帶了爆破的大部分信息，且包絡(luò)線波峰最接近四孔爆破效果，因此考慮選IMF1分量作為主分量，對其包絡(luò)線進行進一步分析，如圖5所示。

由圖5可知，EMD識別法識別出了6個峰值點，分別是(0.216，0.277),(0.256，0.844),(0.267, 1.257),(0.278, 1.837),(0.290, 2.260)和(0.302,0.361)。延期時間分別為40 ms,11 ms,11 ms,12 ms和12 ms。由于本次爆破試驗設(shè)計為12 ms延期的四孔爆破，在此EMD識別法識別出了6個峰值點，很明顯EMD分解出的主IMF分量受到了外界因素的干擾，出現(xiàn)模態(tài)混疊現(xiàn)象，影響了EMD識別法的精度。

圖5 主分量包絡(luò)線Fig.5 The envelope of the principal component

3.2.2 改進算法PEMD識別

由于EMD不完全正交分解，導(dǎo)致分解出的各IMF分量出現(xiàn)模態(tài)混疊現(xiàn)象，造成延期識別精度的降低，而PCA能夠?qū)⒉煌耆坏亩鄠€向量轉(zhuǎn)變成完全正交向量，故在此引入PCA方法對EMD進行改進。

將原始信號進行EMD分解，得到13個IMF分量xj，j=1,2,…,13，計算各分量的相關(guān)系數(shù)矩陣R如表2所示。

表2 各IMF分量相關(guān)系數(shù)Tab.2 The correlation coefficients of IMF components

繼而計算相關(guān)系數(shù)矩陣的特征向量μ

則由特征向量構(gòu)建正交IMF分量yj

對各正交IMF分量進行Hilbert變換，提取包絡(luò)線如圖6所示。

圖6 正交IMF分量及包絡(luò)線Fig.6 Orthogonal IMF component and envelope

由圖6可以看出，IMF1 幅值較大，波形衰減明顯，故選擇為主分量，將其與EMD主分量包絡(luò)線進行對比，如圖7所示。

圖7 主分量包絡(luò)線對比Fig.7 Comparison of principal component envelope

由圖7可知，經(jīng)過PEMD處理后的幅值包絡(luò)圖識別出了四個峰值點，分別是(0.256，0.844)，(0.267, 1.257)，(0.278, 1.837)，(0.290, 2.260)，延期時間分別為12 ms，12 ms，12 ms，誤差為0。與EMD識別結(jié)果相比，有效地去除了模態(tài)混疊現(xiàn)象，提高了識別精度。

3.3 PEMD算法穩(wěn)定性檢驗

試驗控制高程和延期時間兩種變量進行穩(wěn)定性檢驗，選擇同一高程(4號臺階)不同延期時間(10 ms,16 ms,21 ms)起爆和不同臺階(2號，4號，6號臺階)同一延期時間(16 ms)起爆，雷管個數(shù)均設(shè)為3發(fā)，理論上爆破振動信號將由3段爆破振動波疊加而成。其爆破振動速度時程曲線如圖8所示。

圖8 速度時程曲線Fig.8 Velocity time-history curve

首先對同一高程不同延期時間爆破振動信號進行EMD，PEMD識別，如圖9所示。

圖9 不同延期時間包絡(luò)線對比Fig.9 Comparison of envelops with different delay time

由圖9可知，EMD識別法對10 ms，16 ms，21 ms延期信號分別識別出2個、4個、4個峰值點，表明EMD對3段爆破振動波的識別發(fā)生誤差；而PEMD在2號、4號、6號臺階處均識別出3個峰值點，精確識別出了3段爆破振動信號。

繼而對同一延期時間不同高程爆破振動信號進行EMD和PEMD識別，如圖10所示。

由圖10可知，EMD識別法在2號、4號、6號臺階處分別識別出2個、4個、2個峰值點，而PEMD在2號、4號、6號臺階處均識別出3個峰值點，因此PEMD與EMD相比，識別精度不受延期時間和高程因素的影響。為了進一步量化對比兩者的識別誤差，計算出各峰值點的延期時間如表3所示。

圖10 不同高程包絡(luò)線對比Fig.10 Comparison of envelope at different elevations

表3 峰值點延期時間Tab.3 The delay time of the peak point

由表3可知，從不同延期時間分析，EMD識別法在4號臺階處10 ms，16 ms，21 ms的平均誤差為4 ms，PEMD為0；從不同高程分析，在16 ms延期時間的2號、4號和6號臺階處，EMD平均誤差為6 ms，PEMD為0。故無論是從延期時間還是高程進行分析，PEMD識別精度均優(yōu)于EMD。

4 露天邊坡延期爆破試驗

4.1 試驗背景

德興銅礦位于江西省德興市泗洲鎮(zhèn)，黃牛前露天采礦邊坡位于銅廠采區(qū)的東部，邊坡設(shè)計總體坡角為46°，邊坡走向144°～152°，傾向234°～242°，傾角41°～42°。本文采用德興黃牛前露天采礦邊坡開挖爆破數(shù)據(jù)，爆破振動測試所采用的是Blastmate III型號爆破振動監(jiān)測儀，采樣頻率為4 096 Hz。爆破區(qū)域位于銅廠采區(qū)215 m邊坡臺階，爆破測振過程中在黃牛前西側(cè)邊坡170 m，230 m，260 m，290 m水平上分別布置了2號、4號、6號、8號測點，如圖11所示。

圖11 測點布置圖(m)Fig.11 Layout of measuring points(m)

以爆區(qū)坐標(biāo)為基準(zhǔn)，各水平測點的爆心距和高程差，如表4所示。

表4 測點高程差與爆心距Tab.4 The difference between the burst distance and elevation cm

4.2 爆破延期識別

本次爆破為逐孔起爆，延期時間在42～694 ms，炮孔數(shù)為38個，其相應(yīng)的爆破網(wǎng)路及延期設(shè)計時間,如圖12所示。

圖12 爆破網(wǎng)路及延期時間(ms)Fig.12 Blasting network diagram and delay time(ms)

取2號、4號、6號、8號水平測點爆破振動信號進行EMD、PEMD識別，并與爆破延期設(shè)計值進行對比，結(jié)果如圖13所示。

由圖13可知，PEMD與EMD相比，對各水平測點爆破振動信號延期的識別結(jié)果與設(shè)計值更為接近，為了進一步量化對比兩者的精度，定義延期時間的識別率為

圖13 識別精度圖Fig.13 Identification precision diagram

(8)

式中:δ為識別值與理論值的相對誤差；n為炮孔個數(shù)。

計算得到各測點的延期時間識別率，如表5所示。

表5 測點延期時間識別率Tab.5 Recognition rate of delay time of measuring point %

從表5可知，PEMD在各測點處的識別率均高于EMD，且EMD識別率在74%～91%內(nèi)波動，是因為EMD識別信號時會產(chǎn)生模態(tài)混疊現(xiàn)象，隨著混疊現(xiàn)象程度的不同，識別率也隨著發(fā)生波動；而PEMD由于解決了EMD模態(tài)混疊的問題，識別率大幅提升，且穩(wěn)定保持在90%以上。由于PEMD對炮孔的延期識別率大幅提升，將各炮孔的延期識別值與設(shè)計值進行對比，可識別爆破過程中是否出現(xiàn)盲炮，同時結(jié)合炮孔布置圖，能夠進一步確定盲炮的具體位置，對解決爆破工程中的盲炮識別問題具有重要意義。

5 結(jié) 論

為了提高EMD爆破延期識別精度，本文提出了一種基于PCA的PEMD方法，利用PCA能將具有相關(guān)性的數(shù)組轉(zhuǎn)化為正交數(shù)組的特性，抑制了振動信號EMD分解過程中的模態(tài)混疊現(xiàn)象，提高了爆破延期識別的精度。

通過相似物理模型試驗分析結(jié)果可知，PEMD 能夠有效地抑制振動信號EMD分解時出現(xiàn)模態(tài)混疊現(xiàn)象，爆破振動波峰值點的識別精度顯著提高，并通過控制高程和延期時間變量對PEMD方法的穩(wěn)定性進行了檢驗；同時以露天邊坡延期爆破試驗為例，結(jié)果表明PEMD能夠?qū)Ω魉綔y點爆破振動信號實現(xiàn)90%以上的識別率，對后續(xù)爆破工程中爆破參數(shù)設(shè)計優(yōu)化和盲炮識別具有重要的意義。