爆破振動信號的提升小波包分解及能量分布特征＊

路 亮，龍 源，謝全民，2，李興華，紀 沖，趙長嘯

（1.解放軍理工大學工程兵工程學院，江蘇 南京 210007；2.武漢軍械士官學校彈藥修理與銷毀教研室，湖北 武漢 430075）

隨著掘進爆破技術在城市隧道開挖中的大力開發(fā)和利用，爆破施工對既有相鄰隧道結(jié)構和設施的安全影響也逐漸成為人們?nèi)找骊P注的問題，其中爆破產(chǎn)生的地震效應是對隧道結(jié)構安全最具威脅的危害之一。為最大限度地降低爆破地震效應對結(jié)構設施等可能帶來的影響，必須解決好對爆破振動特征的認識、結(jié)構響應預測和結(jié)構破壞或損傷的判定等方面的問題［1－2］，其中，對爆破振動特征的認識是解決爆破地震危害效應的基礎，因此，對爆破振動特征的研究具有非常重要的理論和應用價值。

爆破振動本質(zhì)上是一個非平穩(wěn)的隨機過程，而振動信號經(jīng)過復雜的巖土介質(zhì)傳播后，往往摻雜著各種頻率成分的干擾波，因此，建立在平穩(wěn)隨機過程基礎上的傅里葉變換方法無法準確提取和分析信號的振動特征［3］。近年來，小波變換以其良好的時頻局部化特性，為爆破振動特征的提取和分析提供了一條有效的途徑［4－5］。然而，小波方法由于不能根據(jù)振動信號的特點選取適應的小波基來準確識別振動信號特征，致使經(jīng)典（一代）小波在爆破振動信號的特征提取和分析方面存在局限。

提升算法（lifting scheme）是一種二代小波構造方法，它在繼承了經(jīng)典小波變換多分辨率特性的基礎上，通過選擇不同的提升算子來改變小波濾波器的特性，得到不同性質(zhì)的雙正交小波，使小波的設計具有更大的靈活性。根據(jù)工程應用的需要，通過提升方法可獲得具有某種希望特性的小波，從而構造出與信號特性更加匹配的小波濾波器組［6－9］。本文中將以提升算法為基礎，依據(jù)小波包變換的思想，對信號中的細節(jié)分量進一步分解，實現(xiàn)對爆破振動信號的提升小波包 （second generation wavelet packet，SGWP）變換，進而能夠準確地提取爆破地震波中不同頻帶下的振動分量，研究各頻帶下爆破地震波能量的分布特征，更好地揭示爆炸應力波的傳播及衰減規(guī)律，為深入研究城市隧道結(jié)構在爆破地震效應下的振動響應提供依據(jù)。

1 提升算法

提升算法的分解過程可歸納為［8］：

（1）分裂。將信號X 分為偶序列Xe和奇序列Xo。

（2）預測。用相鄰的偶序列Xe和預測器P估計信號的奇序列Xo，并將估計誤差d作為信號的細節(jié)分量，反映了信號中的高頻成分

（3）更新。用細節(jié)分量d和更新器U，更新偶序列Xe以獲得信號的逼近分量c，代表了信號中的低頻成分

提升算法的重構過程為分解過程的逆運算，可直接由分解過程得到，重構過程的預測器P和更新器U與分解過程相同。

2 基于提升小波包分解算法的爆破振動信號頻帶能量表征

2.1 提升小波包分解算法關鍵技術

2.1.1 小波基的選擇

小波基的選擇是用小波包方法對信號進行分析時必須要考慮的問題，因為不同的小波基分析同一個信號會產(chǎn)生不同的結(jié)果。對爆破振動信號進行小波包分析時，小波基函數(shù)的選擇一般要遵循具有緊支撐性、對稱性和光滑性的原則。

在提升算法框架下，采用插值細分原理［10］設計預測器和更新器，可獲得雙正交小波的小波函數(shù)和尺度函數(shù)，一般記這種小波為二代小波（second generation wavelet），記為SGW（N，N～），其中 N、N～分別為小波的預測器和更新器個數(shù)。這種小波是對稱的、緊支撐的，并且具有沖擊衰減形狀［8，11］，能較好地與爆破振動信號相匹配，從而可有效抑制變換時振動信號的相位失真。

圖1 插值細分示意圖Fig.1 Sketch of interpolation subdivision

2.1.2 基于插值細分的二代小波的構造

插值細分的本質(zhì)是在原始樣本的基礎上采用多項式插值方法來獲得新的樣本值，即當位于2個相鄰樣本中間位置的插值點已知時，若插值多項式確定即可求出新的樣本值，因此，預測器系數(shù)可以通過插值多項式來確定。插值細分示意圖如圖1所示。

根據(jù)Lagrange插值定理［12］，已知N＋1個點x0，x1，…，xN的函數(shù)值為y0，y1，…，yN，且 yi＝f（xi），i＝0，…，N，則存在唯一一個次數(shù)不大于N的多項式 Ln（x），使得 Ln（xi）＝f（xi），那 么

對于一個已知的插值點x，Ln，i（x）是常數(shù)，它與函數(shù)值yi無關，插值點x對應的函數(shù)值是已知的N＋1個函數(shù)值y0，y1，…，yN的線性組合。令Pi＝Ln，i（x），x處對應的函數(shù)值為y，則該式說明新的函數(shù)值是由已知的一組函數(shù)值預測得到的，當x0，x1，…，xN和x確定后，｛Pi｝是唯一的，且∑Pi＝1。

每一次細分時，取 N（N＝2D，D 為正整數(shù)）個已知的樣本yj，k－D＋1，…，yj，k，…，yj，k＋D，假設這些樣本是等時間間隔采樣的，它們對應的采樣時刻分別為xk＋1，xk＋2，…，xk＋N，其中xk為任意的起始時間，細分產(chǎn)生的新的采樣值處于這些已知樣本的中間位置。插值點為：x＝xk＋（N＋1）／2，這樣預測系數(shù)可由式Ln，i（x）確定，即

根據(jù)式（3）即可求得幾種二代小波的預測器系數(shù)，如表1所示。

表1 預測器系數(shù)Table1 Predict coefficient of second generation wavelets

文獻［9，13］中指出更新器系數(shù)為預測器系數(shù)的一半，因此，當更新器系數(shù)的個數(shù)確定后，便可根據(jù)預測器的系數(shù)來確定更新器的系數(shù)。

預測器P和更新器U確定后，分別根據(jù)式（1）、（2）經(jīng)過迭代運算后即可生成尺度函數(shù)φ（x）與小波函數(shù)φ（x）。圖2為經(jīng)過10次迭代后生成的SGW（6，6）的尺度函數(shù)圖和小波函數(shù)圖。

圖2 SGW（6，6）的尺度函數(shù)與小波函數(shù)Fig.2 Scaling function and wavelet function of SGW（6，6）

從圖2中可以看出，SGW（6，6）很好地滿足了小波基選擇必須遵循的原則，具有良好的對稱性和緊支撐性，且具有沖擊衰減形狀。文獻［8，11，14］中指出，當 N和較小時，尺度函數(shù)和小波函數(shù)的支撐區(qū)間較?。环粗?，支撐區(qū)間較大。一般支撐區(qū)間小的小波函數(shù)適合處理非平穩(wěn)信號，小波系數(shù)能夠有效地描述信號的瞬態(tài)分量，而支撐區(qū)間大且連續(xù)較好的小波適合于描述穩(wěn)態(tài)信號。因此，本文中將采用SGW（6，6）作為提升小波包變換的小波基。

2.1.3 提升小波包變換的移頻算法

小波包分解逐層隔點采樣時，采樣率會降低1／2。當采樣頻率低于奈奎斯特（Nyquist）頻率時，高頻信號會發(fā)生頻率折疊，對其繼續(xù)分解會造成頻率混疊現(xiàn)象。文獻［15－16］中根據(jù)混頻的原因提出了一種移頻算法，其原理為：將高頻成分進行移頻處理，使其所含的最高頻率低于1／2奈奎斯特頻率，以避免頻率折疊。根據(jù)傅里葉變換的移頻特性，要使信號x（t）的傅里葉變換x＾（t）頻率降低f0，只需將x（t）乘上eiπf0t；設信號的采樣頻率為fs，則其最高分析頻率為fs／2，由于隔點采樣，小波包分解到第j層，采樣頻率將為2－jfs，最高分析頻率將為2－j－1fs，則高頻小波包分解系數(shù)應移頻2－j－1fs。

以一段模擬信號為例，信號的波形及時頻圖如圖3所示。選取SGW（6，6），按照提升小波包分解算法，分別采用移頻和不移頻算法對模擬信號進行2層分解，分解后第2層的節(jié)點系數(shù)如圖4所示。

圖3 模擬信號波形曲線及三維時頻圖Fig.3 Wave curve and time－frequency map of simulated signal

由圖4可知，兩種分解算法得到不同的分解結(jié)果，為說明移頻算法的合理性，選?。?，0）節(jié)點的小波包分解系數(shù)進行單支重構，即將其余節(jié)點的系數(shù)置零，按照提升算法重構過程進行小波包重構。單支重構后的信號如圖5（a）所示。

對利用節(jié)點（2，0）單支重構后的信號進行短時傅里葉變換［17］，得到如圖5（b）所示的時頻圖。對比后可知，采用移頻算法的提升小波包變換可以將原信號按頻帶精確劃分，重構后的信號僅包含0～128Hz的成分；而采用不移頻算法重構的信號，由于混頻現(xiàn)象致使不能對原信號按頻帶準確劃分，分解結(jié)果存在頻率失真。

圖4 移頻與不移頻分解結(jié)果對比Fig.4 Decomposition results with frequency－shift and no frequency－shift

圖5 節(jié)點（2，0）單支重構結(jié)果Fig.5 The signal－branch reconstruction result of node（2，0）

2.2 爆破振動信號各頻帶能量表征

信號s經(jīng)提升小波包i層分解后，可以得到2i個頻帶上的子空間信號，則s可由這些子空間的正交和表示，即

式中：fi，j為爆破振動信號小波包分解到第i層第j個節(jié)點上的重構信號。若s的最低頻率為1，最高頻率為ω，則第i層每個頻帶的頻率寬度為ω／2i。

根據(jù)Parseval定理［18］，由式（4）可得爆破振動信號的能量譜為

式中：m 為爆破振動信號的采樣點數(shù)，xj，k（j＝0，1，2，…，2i－1，k＝1，2，…，m）為重構信號fi，j的幅值，Ei，j為爆破振動信號分解到第i層第j個節(jié)點的頻帶能量。

由式（5）可知，信號s的總能量為

信號分解到第i層時，各頻帶能量占總能量的百分比為

由式（4）～（5）可知，爆破振動信號由提升小波包分解成不同頻帶的振動分量，從而能反映頻率在爆破振動時的影響，且頻帶能量又能同時反映該頻帶爆破振動分量的強度和作用時間。因此，在此基礎上獲得的頻帶能量能同時反映爆破振動信號強度、頻率及作用時間的共同作用影響［19］。

3 基于SGWP算法的爆破振動信號能量分布特征

3.1 振動測試

圖6 實測振動信號的時程曲線Fig.6 Time－h(huán)istory curve of measured vibration signal

圖6為結(jié)合某城市隧道開挖工程采集的一段爆破振動波形。本次實驗采用EXP－4850型爆破振動測試儀，其最小工作頻率為1Hz。因此，儀器的頻率響應性能完全滿足信號的測量要求。

3.2 爆破振動信號的提升小波包分解

選用2.1.2節(jié)構造的SGW（6，6）作為小波基，按照提升算法，對爆破振動實測信號進行5層小波包分解，分解結(jié)果見圖7。根據(jù)式（5）編寫計算程序，可將前16節(jié)點對應的不同頻帶內(nèi)的能量值列于表2中，圖8為按照式（6）～（7）求得各頻帶內(nèi)的能量百分比分布，通過該圖可直觀比較振動信號中的能量分布情況。

表2 主要頻帶能量分布表Table2 Energy distribution of main frequency band

圖7 提升小波包分解結(jié)果Fig.7 Decomposition result of lifting wavelet packet

圖8 爆破振動信號不同頻帶能量百分比分布Fig.8 The percentage of energy for blasting vibration signal in different frequency bands

3.3 爆破振動信號不同頻帶能量分布規(guī)律

從表2、圖8中可以看出，信號在0～256Hz間的能量占到總能量的99.4%，表明所選爆破振動信號的能量在頻域范圍內(nèi)雖然分布較廣，但大多數(shù)集中在0～256Hz之間。

信號中頻率在256Hz以上的能量在總能量中的所占比例只有0.6%，說明隨著雷管段數(shù)的增加和振動作用時間的延長，信號中高頻成分所占比例減少，能量更加集中在主振頻帶附近。

（3）通過表2可知，信號在其頻域中除含有一個主振頻帶外，還存在多個子頻帶，子頻帶中又有各自主頻且能量分布很不均勻，直接反映在爆破振動信號中出現(xiàn)多個峰值。由此可見，在爆破振動安全判據(jù)中，只參考單個主振頻率有欠準確，而應將多主頻的影響考慮進去。

圖9 重構信號與實測信號的相對誤差Fig.9 Relative error of reconstruction and measurement signal

3.4 SGWP算法有效性檢驗

依據(jù)SGWP分解算法逆運算對圖8中的各節(jié)點系數(shù)進行信號重構，得到重構信號與實測信號的相對誤差ε如圖9所示。從中可以看出，重構信號與實測信號的誤差ε量級均在10－6以上，對工程分析、計算的影響可忽略不計。因此，本文中用SGWP算法對爆破振動信號進行信號分解的過程中，信號的能量損失可忽略不計，因此本文中所用的分析方法是有效的，可以保證真實反應爆破振動信號中的能量分布情況。

4 結(jié) 論

（1）根據(jù)小波包變換具有多分辨率分析的優(yōu)點，通過設計基于插值細分的二代提升小波基，實現(xiàn)對爆破振動信號的提升小波包分解，較好地解決了變換時振動信號的相位失真；

（2）針對信號經(jīng)小波包分解采樣率會降低的特點，通過引入移頻算法，有效避免了采樣頻率低于Nyquist頻率時，高頻信號頻率發(fā)生折疊造成的頻率混疊現(xiàn)象，通過模擬信號的單支重構證明了移頻算法的合理性；

（3）依據(jù)提升小波包分解算法，通過將爆破振動信號分解到不同頻帶上，可以對不同頻率范圍內(nèi)振動分量的時間變化規(guī)律加以分析，以給出爆破振動信號的時－頻分布特征，為爆破振動研究提供了新的分析手段；

（4）實測信號的能量特征研究表明，基于提升小波包算法的能量特征分析方法可以很好地與爆破振動等非平穩(wěn)信號相匹配，不僅能準確地了解爆破振動信號的能量－頻率分布，還可以給出不同頻帶上的振動分量分布和時間衰減規(guī)律，為指導結(jié)構抗震設計和工程爆破監(jiān)測提供分析依據(jù)。

［1］Dowding C H.Blasting vibration monitoring and control［M］.Englewood Cliffs：Prentice－Hall，1985：1－5.

［2］晏俊偉，龍源，方向，等.基于小波變換的爆破振動信號能量分布特征分析［J］.爆炸與沖擊，2007，27（5）：405－409.Yan Jun－wei，Long Yuan，F(xiàn)ang Xiang，et al.Analysis on features of energy distribution for blasting seismic wave based on wavelet transform［J］.Explosion and Shock Waves，2007，27（5）：405－409.

［3］何軍，于亞倫，梁文基.爆破振動信號的小波分析［J］.巖土工程學報，1998，20（1）：47－50.He Jun，Yu Ya－lun，Liang Wen－ji.Wavelet analysis for blasting seismic signals［J］.Chinese Journal of Geotechnical Engineering，1998，20（1）：47－50.

［4］謝全民，龍源，鐘明壽，等.基于小波、小波包兩種方法的爆破振動信號對比分析［J］.工程爆破，2009，15（1）：5－9.Xie Quan－min，Long Yuan，Zhong Ming－shou，et al.Comparative analysis of blasting vibration signal based on wavelet and wavelet packets transform［J］.Engineering Blasting，2009，15（1）：5－9.

［5］張耀平，曹平，高賽紅.爆破振動信號的小波包分解及各頻段的能量分布特征［J］.金屬礦山，2007（11）：42－43.Zhang Yao－ping，Cao Ping，Gao Sai－h(huán)ong.Wavelet packet decomposition of blasting vibration signals and energy distribution characteristics of frequency bands［J］.Metal Mine，2007（11）：42－43.

［6］張德豐.MATLAB小波分析［M］.北京：機械工業(yè)出版社，2009.

［7］姜洪開，何正嘉，段晨東，等.基于提升方法的小波構造及早期故障特征提?。跩］.西安交通大學學報，2005，39（5）：494－498.Jiang Hong－kai，He Zheng－jia，Duan Chen－dong，et al.Wavelet construction based on lifting scheme and incipient fault feature extraction［J］.Journal of Xi’an Jiaotong University，2005，39（5）：494－498.

［8］Sweldens W.The lifting scheme：A construction of second generation wavelet［J］.SIAM Journal on Mathematics Analysis，1997，29（2）：511－546.

［9］Claypoole R L，Davis G M，Sweldens W，et al.Nonlinear wavelet transforms for image coding via lifting［J］.IEEE Transactions on Image Processing，2003，12（12）：1449－1458.

［10］Daubechies I，Sweldens W.Factoring wavelet transforms into lifting steps［J］.Journal of Fourier Analysis and Application，1998，4（3）：247－269.

［11］段晨東，張建丁.基于第二代小波變換的轉(zhuǎn)子碰摩故障特征提取方法研究［J］.機械科學與技術，2006，25（10）：1229－1232.Duan Chen－dong，Zhang Jian－ding.Study of the fault feature extraction method for rotor rub－impact based on the second generation wavelet transform［J］.Mechanical Science and Technology，2006，25（10）：1229－1232.

［12］鄧建中，劉之行.計算方法［M］.2版.西安：西安交通大學出版社，2001.

［13］Sweldens W，Schr?der P.Building your own wavelets at home［DB／OL］.http：／／cm.bell－labs.com／who／wim／papes.html／at home，1998－01－05.

［14］Sweldens W.The lifting scheme：A custom－design construction of biorthogonal wavelets［J］.Applied and Computational Harmonic Analysis，1996，15（3）：186－200.

［15］耿中行，屈梁生.小波包的移頻算法與振動信號處理［J］.振動工程學報，1996，19（2）：145－152.Geng Zhong－xing，Qu Liang－sheng.Frequency shift algorithms of wavelet packet and vibration signal analysis［J］.Journal of Vibration Engineering，1996，19（2）：145－152.

［16］劉世元，杜潤生，楊叔子.小波包改進算法及其在柴油機振動診斷中的應用［J］.內(nèi)燃機學報，2000，18（1）：11－16.Liu Shi－yuan，Du Run－sheng，Yang Shu－zi.An improved algorithm for wavelet packets and its applications to vibration diagnosis in diesel engines［J］.Transaction of CSICE，2000，18（1）：11－16.

［17］Alan V O，Ronald W S，John R B.Discrete－time signal processing［M］.劉樹堂，黃建國，譯.西安：西安交通大學出版社，2001：557－582.

［18］周德廉，邵國友.現(xiàn)代測試技術與信號分析［M］.徐州：中國礦業(yè)大學出版社，2005.

［19］中國生，熊正明.基于小波包能量譜的建（構）筑物爆破地震安全評估［J］.巖土力學，2010，31（5）：1522－1528.Zhong Guo－sheng，Xiong Zheng－ming.Safety assessment of structure by blasting seism based on wavelet packet energy spectra［J］.Rock and Soil Mechanics，2010，31（5）：1522－1528.