基于疊加和頻譜分析的電子雷管延期時間研究

龔 敏， 石發(fā)才， 吳曉東

(北京科技大學(xué) 土木與資源工程學(xué)院，北京 100083)

在小凈距隧道爆破開挖中，爆破振動會對已建隧道襯砌和中夾巖造成很大的影響，甚至?xí)鹨欢ǖ陌踩[患[1-2]。對此，使用常規(guī)的導(dǎo)爆管雷管進(jìn)行控制時，由于可用段數(shù)較少，有時效果不好。這種情況下，延期時間可以靈活設(shè)置的電子雷管為有效降低爆破振動提供了可能性[3]。使用電子雷管進(jìn)行微差爆破時，電子雷管延期時間的合理設(shè)置對有效降低爆破振動尤為重要。一般來說，延期時間過大會降低爆破作用的協(xié)同性，延期時間過小會使爆破振動大幅提高。如何識別雷管的實際延期時間和獲得電子雷管的最佳延期時間，從而有效控制爆破振動仍然是當(dāng)前的熱門課題。學(xué)者們對此進(jìn)行了大量研究，取得了一定的成果。

準(zhǔn)確識別雷管的實際延期時間，可以為爆破設(shè)計優(yōu)化提供依據(jù)。目前主要采用基于EMD-HHT和小波分析的奇異性檢測法識別雷管的實際延期時間。如龔敏等[4]利用基于EMD-HHT的瞬時能量法和EMD法對雷管起爆時刻識別進(jìn)行了識別，對比了兩種方法的識別效果。張勝等[5]采用模式自適應(yīng)小波時能密度法對微差爆破的實際延期時間進(jìn)行了確定。對于電子雷管的最佳延期時間，目前主要采用理論計算法、數(shù)值模擬法和現(xiàn)場實測法進(jìn)行確定。理論計算法，即根據(jù)各種理論推導(dǎo)出相關(guān)理論公式，從而根據(jù)炸藥參數(shù)和現(xiàn)場的圍巖參數(shù)等進(jìn)行最佳延期時間的確定。如李清等[6]基于充分形成新自由面的原則，提出了精確控制爆破延期時間的理論計算公式，在工程中得到了很好的應(yīng)用。李順波等[7]基于彈性波理論給出了毫秒延期時間計算公式。數(shù)值模擬法即采用動力分析軟件模擬爆破，對比不同延期時間下監(jiān)測點處的爆破振動，從而進(jìn)行最佳延期時間的確定。如樓曉明等[8]利用LS-DYNA建立了0 ms、17 ms、25 ms和42 ms四種不同延期時間下的臺階爆破模型，發(fā)現(xiàn)延期時間為25 ms時，降振效果明顯?，F(xiàn)場實測法即通過現(xiàn)場做單孔試驗，再基于線性疊加原理和半周期減振原理確定電子雷管的最佳延期時間。如崔正榮等[9]通過疊加現(xiàn)場實測的單孔波形，得到最佳延期時間為23 ms。田振農(nóng)等[10]采用電子雷管起爆，認(rèn)為延期時間設(shè)置為1/2倍主周期時，能有效降低爆破振動。由于現(xiàn)場條件多變，理論計算和數(shù)值模擬有時難以得到可靠的結(jié)果，而現(xiàn)場實測法獲得的結(jié)果更符合實際情況。利用實測法確定電子雷管的最佳延期時間時，如果能建立線性疊加原理和半周期減振原理間的聯(lián)系，便可以綜合利用兩種方法，獲得更加可靠的電子雷管延期時間，從而更有效地降低爆破振動。

本文對蓮塘超小凈距隧道混合起爆網(wǎng)絡(luò)中的振動波形進(jìn)行研究，采用EMD-HHT法和小波分析法進(jìn)行奇異性檢測，從而識別出各段雷管起爆時刻，分離出單孔波形。然后對單孔波形進(jìn)行線性疊加分析和頻譜分析，分析合成波形的特點，觀察各單孔波形的主頻特點。最后探究合成波形振速峰值最小、較小和最大時對應(yīng)的延期時間和主周期間的聯(lián)系，從而得到一種更加可靠的電子雷管最佳延期時間確定方法，為有效降低爆破振動提供依據(jù)。

1 研究方法

1.1 EMD-HHT方法

EMD-HHT方法是在EMD分解的基礎(chǔ)上，利用Hilbert變換得到各分量的瞬時頻譜，從而可同時分析信號時頻特征的一種方法。具體操作步驟如下[11]：

1.1.1 EMD分解

對原始信號X(t)，首先找出所有極值點。然后用三次樣條曲線對所有的極大值點進(jìn)行插值，從而擬合出原始信號X(t)的上包絡(luò)線。同理，擬合出X(t)的下包絡(luò)線。上下包絡(luò)線包含了所有的原始信號數(shù)據(jù)。然后，求出上下包絡(luò)線的均值，并按順序連接得到一條均值線m1(t)，再用X(t)減去m1(t)，得到h1(t)如下

h1(t)=X(t)-m1(t)

(1)

若果h1(t)同時滿足如下條件(1)和(2)，即認(rèn)為h1(t)是一個IMF分量。

(1) 整個數(shù)據(jù)序列中，過零點的數(shù)量與極值點的數(shù)量相等或至多相差1。

(2) 信號上任意一點，由局部極小值點確定的包絡(luò)線和由局部極大值點確定的包絡(luò)線的均值均為0，即信號關(guān)于時間軸對稱。

1.1.2 Hilbert變換

EMD分解完成后，對每個IMF分量c(t)作Hilbert變換，得到H(c(t))

(2)

式中:PV代表柯西主值。然后構(gòu)造信號z(t)

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

(3)

式中:a(t)為幅值函數(shù)，Φ(t)為相位函數(shù)。

定義Hilbert譜的表達(dá)式為

(4)

定義Hilbert瞬時能量譜的表達(dá)式為

(5)

1.2 小波分析

小波分析方法是一種時間窗和頻率窗都可以改變的時頻局部化分析方法。

設(shè)ψ(t)∈L2(R)，其傅里葉變換為F(w)。當(dāng)F(w)滿足

(6)

時，稱ψ(t)為一個母小波。將母小波平移和伸縮后，就可以得到一個小波序列。

連續(xù)小波為

(7)

對于任意函數(shù)f(t)∈L2(R)的連續(xù)小波變換為

(8)

實際中常用的小波為離散2進(jìn)小波，即a=2j，b=2jk(j，k∈Z)。

1.3 傅里葉變換

傅里葉變換是可以將信號在時域和頻域間相互轉(zhuǎn)化的工具。傅里葉變換的實質(zhì)是將f(t)分解成不同頻率正弦波的疊加，這樣就把對f(t)的研究轉(zhuǎn)化為對其加權(quán)系數(shù)F(w)的研究。通過將爆破振動時域信號轉(zhuǎn)換為頻域信號，可以很好地分析信號的頻域特征。

函數(shù)f(t)∈L2(R)的連續(xù)傅里葉變換定義為

(9)

實際應(yīng)用中常用的是離散傅里葉變換(DFT)，離散時間序列{fn}的DFT定義如下

(10)

式中:k=0，1，…，N-1。

1.4 疊加分析

為了進(jìn)行疊加分析，從原振動波形中截取得到單孔波形。單孔波形的持續(xù)時間為T0，采樣間隔為dt，每個采樣點的振速值為g′(t)。由于采集到的波形為離散數(shù)據(jù)點，為了程序計算方便，將單孔波形數(shù)據(jù)點擬合成傅里葉級數(shù)形式

(11)

式中:t為時間；n為級數(shù)；w0為基頻；a0、an，bn為擬合參數(shù)。表達(dá)式如下

w0=2π/T0

(12)

(13)

將g(t)擴(kuò)展到整個時間域，得到

(16)

根據(jù)Anderson、Hizen和謝烽等[12～14]的研究，可知利用線性疊加原理預(yù)測爆破振動是可行的。假設(shè)每個孔藥量相同，單孔波形為f(t)，微差爆破的延期時間為Ti，則m個炮孔微差爆破后的振動波形可表示為如下形式

(17)

通過上述疊加，可以得到不同延期時間下合成波形f′(t)的振速峰值及其在合成波形中出現(xiàn)的時刻，從而可以研究延期時間和合成波形振速峰值間的關(guān)系。當(dāng)合成波形的振速峰值最小時，對應(yīng)的延期時間即為最佳延期時間。

2 信號的奇異性檢測

2.1 小波基函數(shù)的選取

選擇合適的小波基函數(shù)對于提高信號奇異性檢測效率有很大的幫助。一般可根據(jù)小波基函數(shù)的性質(zhì)(緊支性、正則性、對稱性和消失距)進(jìn)行小波基函數(shù)的選取。根據(jù)研究[15-17]，小波基函數(shù)衰減越快，局部化特性就越好。正則性表現(xiàn)為小波基函數(shù)的可微性，一定的正則性和足夠高的消失距對有效檢測信號奇異性非常重要。而過高階的消失距會增加計算量，并使分析結(jié)果模糊。對稱性可保證信號重構(gòu)時的精確性，在進(jìn)行信號奇異性檢測時，可放松這方面的要求。

根據(jù)上述要求，可選擇db3小波基函數(shù)進(jìn)行信號的奇異性檢測。相對于其他db系列小波，db3小波基函數(shù)緊支集較小，時域緊支性好，近似對稱，有3階消失距，有利于檢測信號的奇異性，同時計算量較小。

2.2 奇異性檢測實例

某次爆破開挖的實際炮孔布置見圖1。為了防止風(fēng)鉆鉆孔角度問題帶來的欠挖，設(shè)置圖1中所示的補(bǔ)償炮孔(未填充的孔)。為了使掌子面底部爆破達(dá)到單孔單響的要求，增加一個550 ms的電子雷管。1 050 ms和1 170 ms為設(shè)計的單孔試驗。為了防止串段和獲得足夠的單孔波形，設(shè)計非電雷管跳段使用。在水平方向距離開挖輪廓線3.5 m處測得的X方向振動波形見圖2，從圖2可以大概看出各段振動波形有明顯的區(qū)分，其中振動信號的采樣頻率為8 000 Hz。為了精確確定各段的起爆時刻，對X方向的振動波形進(jìn)行EMD-HHT和db3小波分析(分解8層)，得到Hilbert瞬時能量譜圖3、主成分分量IMF2的Hilbert包絡(luò)曲線圖4、小波分析時能密度曲線圖5和小波變換模值(第3層)圖6。各方法的識別結(jié)果統(tǒng)計于表1。由于電子雷管的誤差為±1.5 ms，所以電子雷管起爆時刻550 ms和1 050 ms可以作為很好的識別性能評價標(biāo)準(zhǔn)。從表1中不難看出，EMD法(主成分包絡(luò)法)和小波分析模極大值法較為精確地識別出了非電雷管各段的起爆時刻，識別結(jié)果符合非電雷管段數(shù)及延期時間誤差圖7，但是四種方法對于延期時間為6 ms的電子雷管(1 250～1 424 ms)的識別均較為粗糙。從圖3和圖5可以看出，電子雷管(1 250～1 424 ms)爆破的瞬時能量較大，所以有必要充分利用電子雷管的延時精確性，錯相減振，降低爆破振動，而合理設(shè)置電子雷管延期時間是達(dá)到上述要求的關(guān)鍵。結(jié)合表1和圖2可以看出，第五段和第十一段有100 ms左右的衰減時間，基本衰減為0，而且前一段振動基本衰減為0。因此可以認(rèn)為第五段、第十一段以及1 050 ms處的波形為完整的單孔波形，且沒有疊加其他波形(除噪音)。從而得到第五段、第十一段和1 050 ms處X方向的單孔波形見圖8，由于篇幅限制，其他方向的單孔波形不再給出。

圖1 實際炮孔布置圖

圖2 X方向振動波形

圖3 Hilbert瞬時能量譜

圖4 主成分包絡(luò)曲線(IMF2)

圖5 時能密度曲線

圖6 小波變換模值(第三層)

段別瞬時能量法/msEMD/法ms時能密度法/ms模極大值法/ms16.81.823.52.9 352.142 55.041.85104.5104.3125.5 102.77199.8196.8205.9 202.5 9278.1 266.8282.4 268.311457.5460.3461.3445.9550 ms558.3 550.3554.5 552.613629.1641.1635.2629.515850.8、930.5 845.5、932.1 839.2、933.1 838.9、928.71 050 ms1 065.61 052.31 058.51 052.4

3 疊加分析和頻譜分析

3.1 疊加分析

現(xiàn)場實測波形一般是多個炮孔爆炸后各波形的疊加結(jié)果，但由于爆破振動信號非周期平穩(wěn)信號，即使假設(shè)每段產(chǎn)生的單孔波形相同，段數(shù)不一樣，相同延期時間下疊加波形的振速峰值也可能不一樣。因此，利用疊加法分析延期時間和合成波形振速峰值間的關(guān)系時，有必要先消除疊加的波列數(shù)對合成波形振速峰值的影響。

圖7 非電雷管延期時間誤差

圖8 各段X方向單孔波形

依據(jù)上述疊加原理，編寫MATLAB程序，其中擬合級數(shù)N取300，每兩列波間的延期時間為固定值T，T在區(qū)間[1,T0]上取整數(shù)值，疊加的波列數(shù)取區(qū)間[2,80]內(nèi)的整數(shù)。疊加第五段單孔波形并統(tǒng)計任意延期時間下合成波形的振速峰值不再變化的波列數(shù)，得到圖9。從圖9中可以看出，不同延期時間下，當(dāng)波列數(shù)達(dá)到2列后，合成波形的振速峰值不再變化的概率為0.793 81；當(dāng)波列數(shù)達(dá)到15列后，合成波形的振速峰值不再變化的概率為0.986 25；當(dāng)波列數(shù)為47列后，合成波形的振速峰值不再變化的概率為1。因此，當(dāng)疊加的波列數(shù)為15列時，就可以不考慮波列數(shù)對合成波形振速峰值的影響，僅研究延期時間和合成波形振速峰值間的關(guān)系。

圖9 大于某個波列數(shù)后合成波形振速峰值不再變化的概率

Fig.9 Probability that the peak vibration velocity of synthetic waveform will not change when the number of single waveforms is greater than a certain value

對第五段X方向單孔波形進(jìn)行疊加分析，疊加結(jié)果見圖10、圖11、圖12和圖13。從圖12可以看出當(dāng)延期時間為17 ms時，合成波形的振動峰值首次出現(xiàn)最小值，振速峰值出現(xiàn)在合成波形中的3.75 ms處，說明合成波形的振速峰值即為第一個單孔波形的振速峰值，后續(xù)波形有所抵消。從圖10和圖12可以看出，當(dāng)延期時間為7 ms時，合成波形振速峰值為所有峰值中的較小值，振速峰值出現(xiàn)在合成波形中的8.875 ms處，說明第一個波形和第二個波形有所疊加，而后期波形有所抵消。從圖11和圖12可以看出，當(dāng)延期時間為13 ms(除1 ms)時，合成波形振速峰值為所有峰值中的最大值，振速峰值出現(xiàn)在合成波形中的68.75 ms處，說明后期波形有疊加增強(qiáng)。當(dāng)延期時間大于30 ms時，合成波形的振速峰值基本為最小值。

圖10 第五段X方向合成波形(延期時間7 ms)

Fig.10 Synthetic waveform of fifth stage inXdirection (delay time is 7 ms)

圖11 第五段X方向合成波形(延期時間13 ms)

Fig.11 Synthetic waveform of fifth stage inXdirection (delay time is 13 ms)

圖12 第五段X方向延期時間與疊加峰值

Fig.12 Delay time and superposition peak value of fifth stage inXdirection

圖13 第五段X方向延期時間與峰值出現(xiàn)時刻

Fig.13 Delay time and peak occurrence time of fifth stage inXdirection

同樣的方法疊加第五段、第十一段和1 050 ms三個方向上的單孔波形，統(tǒng)計疊加結(jié)果于表2。從表2和圖12可以看出，當(dāng)15列波形疊加后，合成波形的最小峰值一般為單孔波形的峰值，出現(xiàn)在首個單孔波形的峰值處。而最大或較大峰值一般為多個單孔波形疊加后產(chǎn)生的，大于單孔波形的峰值，說明最大峰值對應(yīng)的延期時間下，在一定時間內(nèi)，波形有所疊加增強(qiáng)。從表2還可以看出，不同方向上根據(jù)疊加原理得到的最佳延期時間是不同的。所以，利用疊加原理確定雷管的最佳延期時間時，不能僅考慮一個方向單孔波形的疊加結(jié)果，而應(yīng)綜合考慮多個方向單孔波形的疊加結(jié)果。

表2 各段各方向單孔波形疊加結(jié)果

3.2 頻譜分析

利用快速傅里葉變換，可以將信號的時域信息和頻域信息聯(lián)系起來，從而可以從頻域信息來考察爆破振動信號的特征。利用MATLAB編寫程序，實現(xiàn)振動信號的快速傅里葉變換(FFT)，從而可以得到信號的頻譜圖。圖14為第五段X方向的頻譜圖，從圖中可以看出，第五段X方向單孔波形的主頻率為78.2 Hz，即主周期為12.8 ms。對第五段、第十一段和1 050 ms不同方向上的單孔波形作快速傅里葉變換，頻譜分析結(jié)果統(tǒng)計于表3。從表3可以看出，不同段單孔波形在同一方向上的主頻基本穩(wěn)定，不同方向上的主頻不同。結(jié)合表2和表3，可以得到振動峰值最小、較小和最大時對應(yīng)的延期時間與主周期的比值，結(jié)果見表4。從表4中可以看出，當(dāng)延期時間為主周期的0.4倍～0.7倍(0.5倍左右)、1.4倍～1.7倍(1.5倍左右)和2.4倍～2.7倍(2.5倍左右)時峰值振速最小或較小，這說明錯相降振是可行的，但是合成波形振速峰值最小值通常為單孔波形的振速峰值。除延期時間為1 ms的情況外，當(dāng)延期時間為主周期的0.9倍～1.2倍(1倍左右)時，合成波形的振速峰值達(dá)到最大，說明此時各單孔波形疊加增強(qiáng)。

圖14 第五段X方向頻譜圖

段別方向主頻/Hz主周期/ms0.4倍～0.7倍主周期/ms第五段X78.212.85.12～8.96Y101.79.83.92～6.86Z109.59.13.64～6.37第十一段X78.212.85.12～8.96Y101.79.83.92～6.86Z109.59.13.64～6.371 050 msX54.7418.37.32～12.81Y101.79.83.92～6.86Z132.97.53～5.25

3.3 最佳延期時間的確定

根據(jù)上述分析，當(dāng)延期時間取半個周期時并不能完全保證合成波形振速峰值最小，只能說延期時間為主周期的0.4倍～0.7倍、1.4倍～1.7倍或2.4倍～2.7倍時，合成波形振速峰值會有較小或最小值。但是當(dāng)延期時間為主周期的1倍左右時，合成波形的振速峰值必然會疊加增強(qiáng)。為了提高減振可靠性，使合成波形振速峰值最大可能地降低，提出了如下確定最佳延期時間的方法。首先通過線性疊加原理來確定最佳延期時間，根據(jù)多個單孔波形的疊加結(jié)果，考慮主要的控制要求，得到備選延期時間。然后，根據(jù)快速傅里葉變換求出單孔波形的主周期，判斷備選延期時間是否符合半周期減振原理，即備選延期時間是否為各單孔波形主周期的0.4倍～0.7倍、1.4倍～1.7倍或2.4倍～2.7倍。若符合，即得到最佳延期時間。若不符合，重新選擇較優(yōu)的延期時間，并重復(fù)上述步驟，直至找到最佳延期時間，這樣通過兩種原理相輔相成，得到更加可靠的最佳延期時間，有效地降低爆破振動。具體流程,如圖15所示。

表4 相應(yīng)延期時間與主周期的比值

Tab.4 The ratio of the corresponding delay time to the main period

段別方向峰值最小峰值較小峰值最大第五段X1.330.551.02Y1.330.511.02Z0.66、1.430.550.99第十一段X0.55、1.331.09Y2.651.431.02Z0.43～0.771.42～1.651.211 050 msX1.530.770.99Y2.651.631.12Z0.40.670.93

圖15 最佳延期時間確定方法

本實驗單孔波形主周期的0.4倍～0.7倍如表3所示。根據(jù)上述方法，為了保證爆破協(xié)同性，最終確定的最佳延期時間為6 ms，滿足主要控制要求，且與各方向單孔波形主周期的比值基本為0.47～0.66，符合干擾減振原理。本次試驗中以6 ms作為電子雷管延期時間，在2 m進(jìn)尺爆破中，測得距離開挖輪廓線垂直距離1.8 m且水平距離3.5 m處各方向的振速不超過3 cm/s。

4 結(jié) 論

通過上述分析，可以得出以下結(jié)論：

(1) 以電子雷管起爆時刻的識別結(jié)果作為標(biāo)準(zhǔn)，可以較好地評價四種不同的起爆時刻識別方法。其中EMD法和小波變換模極大值法有很好的識別效果，但是這兩種方法需要人為地選擇反應(yīng)信號特征的IMF主成分分量和最佳變換尺度。

(2) 合成波形振速峰值最小、較小或最大時對應(yīng)的延期時間和主周期有很緊密的關(guān)系。一般延期時間為主周期的0.4倍～0.7倍、1.4倍～1.7倍和2.4倍～2.7倍時，得到的合成波形振速峰值為最小或較小，這說明錯相減振是可行的。但是在疊加波形個數(shù)較多，延期時間都相同時，一般只能使合成波形的振速峰值降低為單孔波形的振速峰值。除延期時間為1 ms的情況外，當(dāng)延期時間為主周期的0.9倍～1.2倍時，合成波形的振速峰值達(dá)到最大，各單孔波形疊加增強(qiáng)。基于上述聯(lián)系，提出了一種更加可靠的電子雷管最佳延期時間確定方法。本實驗根據(jù)上述方法得到的最佳延期時間為6 ms。在此延期時間下，測得在2 m進(jìn)尺爆破中距離開挖輪廓線垂直距離1.8 m且水平距離3.5 m處各方向的振速不超過3 cm/s，效果較好。

(3) 同一爆破場地同一方向上單孔波形的主頻基本相同，而不同方向上單孔波形的主頻有所不同。