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

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

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

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

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

針對(duì)上述問題的優(yōu)缺點(diǎn)，本文提出了基于主成分分析(principal component analysis，PCA)的完全正交經(jīng)驗(yàn)?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ù)組的特性，對(duì)IMF分量進(jìn)行主成分分析，實(shí)現(xiàn)各IMF分量之間完全正交化，從而抑制模態(tài)混疊現(xiàn)象[15]，最終提升EMD方法對(duì)爆破振動(dòng)信號(hào)延期識(shí)別的精度；并通過相似物理模擬試驗(yàn)和德興露天臺(tái)階延期爆破試驗(yàn)，與EMD方法進(jìn)行對(duì)比，檢驗(yàn)和評(píng)價(jià)PEMD方法的爆破延期識(shí)別效果。

1 基本原理簡介

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

(1)

式中：ci(t)為第i個(gè)IMF分量；rn(t)為殘差，也稱信號(hào)的趨勢(shì)項(xiàng)。

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

(2)

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

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

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

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

具體步驟如下:

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

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

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

(3)

R=(rij)m×n

(4)

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

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

(5)

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

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

(6)

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

(7)

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

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

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

3 相似物理模擬試驗(yàn)

3.1 相似物理模型

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

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

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

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

3.2 爆破延期識(shí)別

3.2.1 EMD識(shí)別

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

圖3 速度時(shí)程曲線Fig.3 Velocity time-history curve

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

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

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

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

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

3.2.2 改進(jìn)算法PEMD識(shí)別

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

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

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

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

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

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

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

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

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

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

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

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

圖8 速度時(shí)程曲線Fig.8 Velocity time-history curve

首先對(duì)同一高程不同延期時(shí)間爆破振動(dòng)信號(hào)進(jìn)行EMD，PEMD識(shí)別，如圖9所示。

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

由圖9可知，EMD識(shí)別法對(duì)10 ms，16 ms，21 ms延期信號(hào)分別識(shí)別出2個(gè)、4個(gè)、4個(gè)峰值點(diǎn)，表明EMD對(duì)3段爆破振動(dòng)波的識(shí)別發(fā)生誤差；而PEMD在2號(hào)、4號(hào)、6號(hào)臺(tái)階處均識(shí)別出3個(gè)峰值點(diǎn)，精確識(shí)別出了3段爆破振動(dòng)信號(hào)。

繼而對(duì)同一延期時(shí)間不同高程爆破振動(dòng)信號(hào)進(jìn)行EMD和PEMD識(shí)別，如圖10所示。

由圖10可知，EMD識(shí)別法在2號(hào)、4號(hào)、6號(hào)臺(tái)階處分別識(shí)別出2個(gè)、4個(gè)、2個(gè)峰值點(diǎn)，而PEMD在2號(hào)、4號(hào)、6號(hào)臺(tái)階處均識(shí)別出3個(gè)峰值點(diǎn)，因此PEMD與EMD相比，識(shí)別精度不受延期時(shí)間和高程因素的影響。為了進(jìn)一步量化對(duì)比兩者的識(shí)別誤差，計(jì)算出各峰值點(diǎn)的延期時(shí)間如表3所示。

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

表3 峰值點(diǎn)延期時(shí)間Tab.3 The delay time of the peak point

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

4 露天邊坡延期爆破試驗(yàn)

4.1 試驗(yàn)背景

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

圖11 測(cè)點(diǎn)布置圖(m)Fig.11 Layout of measuring points(m)

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

表4 測(cè)點(diǎn)高程差與爆心距Tab.4 The difference between the burst distance and elevation cm

4.2 爆破延期識(shí)別

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

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

取2號(hào)、4號(hào)、6號(hào)、8號(hào)水平測(cè)點(diǎn)爆破振動(dòng)信號(hào)進(jìn)行EMD、PEMD識(shí)別，并與爆破延期設(shè)計(jì)值進(jìn)行對(duì)比，結(jié)果如圖13所示。

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

圖13 識(shí)別精度圖Fig.13 Identification precision diagram

(8)

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

計(jì)算得到各測(cè)點(diǎn)的延期時(shí)間識(shí)別率，如表5所示。

表5 測(cè)點(diǎn)延期時(shí)間識(shí)別率Tab.5 Recognition rate of delay time of measuring point %

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

5 結(jié) 論

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

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