SGWP算法在爆破振動信號分析中的應(yīng)用研究

謝全民，龍 源，郭 濤，張懷智，路 亮，張洋溢

(1.解放軍理工大學(xué) 工程兵工程學(xué)院，南京 210007;2.武漢軍械士官學(xué)校 彈藥修理與銷毀教研室，武漢 430075;

3.廣州軍區(qū)工程科研設(shè)計(jì)所，廣州 510515)

爆破振動信號分析是研究控制爆破振動危害的基礎(chǔ)，也是科學(xué)制定抗震措施的前提［1－3］。因爆破振動導(dǎo)致結(jié)構(gòu)破壞是動力響應(yīng)過程，為使爆破振動安全標(biāo)準(zhǔn)更加科學(xué)化、合理化，須建立綜合考慮振動強(qiáng)度、頻率和持續(xù)時間等多因素的復(fù)合判據(jù)，要求在爆破振動信號分析中對信號本身蘊(yùn)含的時頻、能量等重要特征信息進(jìn)行準(zhǔn)確提取?？紤]到爆破振動監(jiān)測過程中易受外部環(huán)境和測試系統(tǒng)的影響，實(shí)測數(shù)據(jù)中常含噪聲干擾［4］，因此為測試數(shù)據(jù)探尋有效的噪聲去除算法，可提高爆破振動測試信號特征提取的可靠性和準(zhǔn)確性。

工程爆破中為建立大型網(wǎng)絡(luò)化測試系統(tǒng)與控制平臺，實(shí)現(xiàn)爆破振動智能化監(jiān)測與預(yù)報(bào)的工程應(yīng)用目的，要求在傅里葉變換、一代小波變換等傳統(tǒng)信號處理技術(shù)基礎(chǔ)上提高信號分析處理的效率和質(zhì)量，滿足在線處理對算法效率、精度等方面的要求。

Sweldens［5］提出的提升算法被稱作二代小波變換，在繼承經(jīng)典小波多分辨分析特性的基礎(chǔ)上，可實(shí)現(xiàn)數(shù)據(jù)原地運(yùn)算，具有占用空間小，變換速度快等優(yōu)勢。相同數(shù)據(jù)長度下運(yùn)算速度提高近50%。采用二代小波變換技術(shù)可通過設(shè)計(jì)不同的預(yù)測算子和更新算子得到具有某些特殊性能的小波函數(shù)［5－6］，脫離傅里葉變換的局限性，更好地適應(yīng)待分析信號特征。

目前二代小波變換主要應(yīng)用于機(jī)械振動信號分析、故障診斷等方面［7－10］，但在爆破振動信號分析中的應(yīng)用尚未見相關(guān)研究報(bào)道?？紤]到小波包變換技術(shù)在信號處理方面能夠更加完整地保留爆破振動信號局部化特征，同時可提取指定頻率區(qū)間的小波包分量進(jìn)行主分量特征分析［2，11］，本文構(gòu)造了基于插值細(xì)分法的二代小波SGW(6，6)，對二代小波包(Second Generation Wavelet Packet，SGWP)變換在爆破振動信號噪聲濾除及頻譜、能量等重要特征提取中的應(yīng)用進(jìn)行了有益的探索，取得了較為滿意的效果。

1 SGWP算法及SGW構(gòu)造

1.1 SGWP算法

(1)剖分:對第(j，n)(n=1，2，…)，2j節(jié)點(diǎn)系數(shù)進(jìn)行奇偶分裂為

(2)預(yù)測:由式(1)得到第(j+1，2n)節(jié)點(diǎn)的系數(shù)。

式中:p［l］(l=1，2，… ，N)為預(yù)測器系數(shù)。

(3)更新:由式(2)得到第(j+1，2n+1)節(jié)點(diǎn)的系數(shù)。

式中:u［l］(l=1，2，… ，N～)為更新器系數(shù)。

1.1.2 重構(gòu)算法

(1)反更新:由(j+1，2n)，(j+1，2n+1)節(jié)點(diǎn)系數(shù)求 (j，n)節(jié)點(diǎn)的偶系數(shù)。

(2)反預(yù)測:由(j+1，2n+1)節(jié)點(diǎn)系數(shù)和(j，n)節(jié)點(diǎn)的偶系數(shù)求(j，n)節(jié)點(diǎn)的奇系數(shù)。

SGWP變換兩級分解與重構(gòu)過程如圖1所示。

圖1 SGWP兩級分解與重構(gòu)示意圖Fig.1 Two levels decomposition and reconstruction based on SGWP

表1 預(yù)測系數(shù)Tab.1 Predict coefficients

1.2 SGWP改進(jìn)算法

對信號進(jìn)行SGWP變換是逐層隔點(diǎn)采樣，每進(jìn)行下一個分解尺度變換時，采樣率會降低1/2，當(dāng)高頻部分不滿足采樣定理［11］時易發(fā)生頻率折疊，同時由于SGWP算法等效的濾波器為有限長度，不可能具備理想的頻率截止特性［10］，上述因素導(dǎo)致變換過程易出現(xiàn)混頻現(xiàn)象，使分析結(jié)果失真。文獻(xiàn)［12－13］中研究了小波包變換的頻率混疊問題，并通過移頻處理技術(shù)改進(jìn)二代小波包變換算法。

分解算法式(1)、式(2)變?yōu)?

重構(gòu)算法式(3)、式(4)變?yōu)?

1.3 SGW 構(gòu)造

若已知序列 x0，x1，…，xN對應(yīng)函數(shù)值為 y0，y1，…，yN，且滿足yi=f(xi)(i=0，…，N)，則存在唯一1個次數(shù)不大于 N的多項(xiàng)式 Ln(x)，使得 Ln(xi)=f(xi)，那么:

其中:

每次進(jìn)行插值細(xì)分時，取N個(N=2D，D∈Z+)等間隔樣本 yj，k－D+1，…，yj，k，…，yj，k+D，，其對應(yīng)的采樣時刻為xk+1，xk+2，…，xk+N，xk為任意起始時間。通過細(xì)分產(chǎn)生新的采樣值處于已知樣本中間，插值點(diǎn)為:x=xk+(N+1)/2，則預(yù)測系數(shù)可由式(11)確定，即:

當(dāng)N=6時，根據(jù)式(11)計(jì)算其對應(yīng)預(yù)測系數(shù)如表1所示。

設(shè)信號S為一個σ序列，即:

首先構(gòu)造尺度函數(shù)。由剖分原理及圖1中的小波包重構(gòu)部分對序列S進(jìn)行插值細(xì)分，即利用k－1，k，k+1，k+2 四點(diǎn)的 sj，k值預(yù)測 sj+1，2k+1，插值邊界采用補(bǔ)零的方法進(jìn)行處理，經(jīng)過18次迭代可得尺度函數(shù)。

下面構(gòu)造小波函數(shù)。首先需要設(shè)計(jì)更新算子U，因更新算子是為保證小波變換前后信號的低階消失矩不變，當(dāng)N=N～(N為預(yù)測器的個數(shù)，為更新器的個數(shù))時，可直接將預(yù)測系數(shù)除以2作為更新系數(shù)［14］。

其中:

式中，U和D分別為更新算子運(yùn)算系數(shù)和細(xì)節(jié)信號序列。與尺度函數(shù)的構(gòu)造類似，設(shè)細(xì)節(jié)信號為σ序列，近似信號S為零序列，進(jìn)行一次提升格式反變換，信號的偶序列為:

該序列經(jīng)過預(yù)測算子p(·)進(jìn)行恢復(fù)預(yù)測運(yùn)算，得到信號的奇序列為:

對s0進(jìn)行插值細(xì)分，經(jīng)過18次迭代可得到相應(yīng)的小波函數(shù)。

基于提升算法構(gòu)造的尺度函數(shù)和小波函數(shù)是雙正交、對稱、緊支撐的，且具有沖擊形狀［5，14－15］。當(dāng) N 和N～較小時，尺度函數(shù)和小波函數(shù)的支撐區(qū)間較小;反之，支撐區(qū)間較大，具有較好的連續(xù)性。實(shí)際應(yīng)用表明，支撐區(qū)間較小的小波函數(shù)適宜于處理非平穩(wěn)信號，而支撐區(qū)間大且連續(xù)較好的小波適合于描述穩(wěn)態(tài)信號［13］。爆破振動信號具有典型的短時非平穩(wěn)特性［1－4，11］，研究過程中發(fā)現(xiàn) SGW(6，6)適宜于爆破振動信號分析。

2 基于SGWP改進(jìn)算法的爆破振動信號去噪

圖2為爆破振動實(shí)測信號S1，信號采樣頻率fs=1 024 Hz，采樣點(diǎn)數(shù)n=2 048。從爆破振動時程曲線圖2可知，振動波形中混雜有毛刺以及由于測試系統(tǒng)本身帶來的方波干擾。圖3為測試信號S1降噪前的三維時頻譜，分析可知該測試信號中混雜有200～250 Hz頻段內(nèi)的高頻噪聲分量。

圖2 實(shí)測爆破振動含噪信號Fig.2 Measured vibration signal with noise

2.1 SGWP降噪算法

基于SGWP變換的爆破振動測試信號降噪過程如下:

(1)采用移頻算法改進(jìn)的SGWP技術(shù)對測試信號進(jìn)行滿二叉樹的多分辨分解;

(2)對于給定的信息代價函數(shù)，選擇最優(yōu)的二代小波包基;

(3)對最優(yōu)基分解所得的各節(jié)點(diǎn)系數(shù)進(jìn)行閾值量化;

(4)逐層對閾值處理后的節(jié)點(diǎn)系數(shù)進(jìn)行重構(gòu);(5)去噪效果分析。

2.2 基于最優(yōu)基的SGWP分解

采用1.3節(jié)中構(gòu)造的插值小波SGW(6，6)，根據(jù)式(5)、式(6)對S1進(jìn)行分解層數(shù)為3層的滿二叉樹分解。選用對數(shù)能量熵［13］作為信息代價函數(shù)，通過對各小波包變換系數(shù)的信息熵進(jìn)行最優(yōu)基搜索，完成二叉樹的修剪，確定最優(yōu)分解結(jié)構(gòu)如圖4所示。

2.3 節(jié)點(diǎn)系數(shù)閾值量化

為充分保留爆破振動測試信號時頻特征的局部特性，選擇Donoho等［16］提出的軟閾值法進(jìn)行節(jié)點(diǎn)系數(shù)的閾值量化。

其中，N為小波包系數(shù)長度，j為分解尺度。

式(18)中爆破振動信號的噪聲方差σ未知，按下式進(jìn)行估計(jì):

其中，median［·］為中值函數(shù)。

根據(jù)式(17)～式(19)，可對各節(jié)點(diǎn)小波包系數(shù)進(jìn)行軟閾值處理。

2.4 去噪效果分析

圖 5 為圖 4 中確定的最優(yōu)基 d10，d11，d20，d21，d2，2，d23，d30，d31，d32，d33，d36，d37對應(yīng)的二代小波包變換系數(shù)。對上述對應(yīng)的小波包系數(shù)閾值處理后再進(jìn)行逐層重構(gòu)，即得到降噪后的爆破振動信號，波形曲線如圖6所示。

圖5 基于SGW(6，6)最優(yōu)基分解小波包變換系數(shù)Fig.5 Optimum wavelet packets coefficients based on SGW(6，6)

從圖6中濾波后的振動波形可以看出，采用基于SGWP變換的軟閾值去噪算法已基本消除了爆破振動監(jiān)測信號中由于測試系統(tǒng)和環(huán)境噪聲帶來的干擾。降噪后的時程曲線圖6相對圖2中降噪前的爆破振動波形曲線而言更為光滑，且峰值、波形上升和衰減等振動信號局部特征在降噪后的振動信號中表現(xiàn)得更加清晰，所得的有用信息圖6中包含了更多能夠客觀準(zhǔn)確地反映實(shí)測信號的主要成分，為進(jìn)一步提取爆破振動信號特征提供了更加有效的信息。

圖6 降噪后的爆破振動信號Fig.6 Blasting vibration signal after de－noising

為定量評價SGWP算法在爆破振動信號去噪分析中的應(yīng)用效果，引入均方根誤差(RMSE)、信噪比(SNR)、峰值誤差(PE)等三項(xiàng)評價指標(biāo)。

根據(jù)式(20)～式(22)計(jì)算得到圖6中降噪后的爆破振動信號相對于圖2中的實(shí)測爆破振動信號:RMSE= 0. 091 9， SNR =25.239 3，PE=0.028 15?？梢姡?SGWP算法的爆破振動信號去噪獲得了較高的信噪比，且降噪前后誤差較小。

圖7 降噪后時頻譜Fig.7 Time－frequency energy spectrum after de－nosing

3 基于SGWP的爆破振動信號特征提取

3.1 能量特征提取

爆破振動輸入到附近建(構(gòu))筑物中的振動能量是導(dǎo)致結(jié)構(gòu)失效、破壞的重要因素，采用SGWP改進(jìn)算法研究爆破振動信號小波包頻帶相對能量分布特征的思路如下:

設(shè)爆破振動采樣信號{s［i］}(i=1，2，3，4，…，2N)，N∈Z+。經(jīng)分解層數(shù)j=3的SGWP改進(jìn)算法分解后，第(j，n)節(jié)點(diǎn)的系數(shù)記為djn［k］(k=1，2，3，4，…，2N－k)。

定義歸一化能量:

式中，E為信號總能量，即利用E對各頻帶內(nèi)能量進(jìn)行歸一化，相應(yīng)的特征向量稱為歸一化特征向量:

圖8 典型的爆破振動時程曲線Fig.8 Time－curve of blasting vibration

圖8為場地實(shí)驗(yàn)中采用單孔單段起爆方式時采集到的爆破振動信號S2，采樣率為1 024 Hz，采采樣點(diǎn)數(shù)n=1 024。S2已按第2節(jié)中的去噪算法濾除了噪聲分量，表現(xiàn)出典型的沖擊加載瞬態(tài)響應(yīng)和運(yùn)動特征。

圖9為該信號的功率譜曲線。由功率譜圖分析可知 S2的主振頻率處于38～80 Hz、80～122 Hz及150～180 Hz三個頻帶區(qū)間。200 Hz以上的頻率部分包含的功率譜密度較小，譜圖中的分辨率較低，因此圖中只給出了0～256 Hz范圍內(nèi)的功率譜密度分布情況。

根據(jù)式(5)、式(6)，采用 SGW(6，6)對 S2進(jìn)行三層小波包分解，得到8個小波包d30～d37，相應(yīng)的小波包系數(shù)如圖10所示。

圖9 功率譜Fig.9 Power spectrum

圖10 d30～d37層提升小波包系數(shù)圖Fig.10 Lifting wavelet packets coefficients of d30～ d37

圖11是按式(23)、式(24)計(jì)算得到8個小波包對應(yīng)的歸一化能量分布。根據(jù)圖11相對能量分布易知，爆破振動信號的優(yōu)勢能量主要分布在d30，d31，d32三個小波包對應(yīng)的頻帶內(nèi)，d33～d37對應(yīng)頻帶包含的能量較小，與圖9中功率譜密度分布情況一致。由于建(構(gòu))筑物的固有頻率一般較低，當(dāng)爆破振動能量分布趨于低頻帶時容易引發(fā)建(構(gòu))筑物共振而加劇破壞。因此，在爆破工程中可采用SGWP技術(shù)快速準(zhǔn)確地獲取原始信號中不同頻率成份對爆破振動作用的影響，并用以指導(dǎo)爆破安全設(shè)計(jì)。

3.2 時頻特征提取

圖11 相對能量分布Fig.11 Relative energy distribution

為研究爆破振動信號中不同頻率成份分量的振動響應(yīng)機(jī)理，需掌握不同頻率分量隨時程變化產(chǎn)生的幅值衰減、頻譜變化等規(guī)律，可通過提取指定頻帶范圍內(nèi)的振動分量進(jìn)行研究。

實(shí)際應(yīng)用中選取能量占優(yōu)、且頻段處于原爆破振動信號主振頻帶范圍的二代小波包作為爆破振動信號的主分析小波包。根據(jù)圖11相對能量分布，選取d30，d31，d32三個主分析小波包作為該信號的特征包。

對上述三個主分析小波包分量，根據(jù)式(7)、式(8)分別進(jìn)行單支重構(gòu)，即將其余節(jié)點(diǎn)系數(shù)置0，按照移頻算法改進(jìn)的二代小波包重構(gòu)步驟進(jìn)行小波包重構(gòu)，分別獲得對應(yīng)的振動分量如圖12(a)所示。該圖體現(xiàn)了爆破振動測試信號S2主振頻帶分量的峰值出現(xiàn)時刻及時程衰減規(guī)律。d30代表的低頻分量衰減較慢，振動持續(xù)時間較長;d32代表的高頻分量衰減迅速，振動持續(xù)時間較短。

圖12 三個主分析小波包單支重構(gòu)及頻譜分析Fig.12 Reconstruction of main wavelet packets and spectral analysis

分別對d30～d32重構(gòu)所得的振動分量再進(jìn)行FFT變換作頻譜分析如圖12(b)所示。由圖12(b)中頻譜特征可知三個小波包代表的振動分量均具有一定的頻帶寬度，且采用改進(jìn)后的SGWP算法將原信號按獨(dú)立的頻帶區(qū)間進(jìn)行了比較精確地劃分，基本避免了傳統(tǒng)小波包變換過程中容易出現(xiàn)的混頻現(xiàn)象，提高了頻域分析精度。

綜合分析圖9、圖12發(fā)現(xiàn):三個主分析特征包分別對應(yīng)頻段d30(0～64 Hz)、d31(64～128 Hz)、d32(128～192 Hz)，三個特征包重構(gòu)分量的頻譜分別具有峰值f0=40 Hz，f1=96 Hz，f2=161 Hz。由此可見，爆破振動信號的主分析頻帶內(nèi)存在兩個或多個明顯的優(yōu)勢頻帶區(qū)間，共同構(gòu)成爆破振動信號的主振頻帶;同時根據(jù)能量的集中程度，可定義爆破振動信號中存在第一、第二、…主頻率。盡管三個主分析小波包所表征的振動分量幅值和能量存在差異，但上述小波包代表的振動分量處于原信號中能量占優(yōu)的頻帶區(qū)間，且體現(xiàn)的振動衰減特性和頻譜特征均符合爆炸沖擊響應(yīng)規(guī)律，因此可采用其作為原信號的特征包描述爆破振動測試信號的時頻特征。

可見，SGWP算法利用其在信號分析領(lǐng)域中特有的局部化分析能力，快速準(zhǔn)確地提取到了爆破振動信號中主振頻帶內(nèi)的振動分量，可為爆破振動作用下結(jié)構(gòu)的動態(tài)響應(yīng)分析及抗震對策制定等方面提供數(shù)據(jù)參考。

4 結(jié)論

(1)基于改進(jìn)的SGWP技術(shù)能夠快速有效地濾除爆破振動測試信號中包含的噪聲分量，更多地保留測試信號的有效成分，為爆破振動信號特征的精確提取奠定了基礎(chǔ)。

(2)采用SGWP算法成功獲取了爆破振動信號的時頻特征、能量特征，其良好的局部化處理能力適宜于具有短時、非平穩(wěn)特性的爆破振動信號特征提取。

(3)SGWP算法具有速度快，精度高，易于實(shí)現(xiàn)等優(yōu)點(diǎn)，在爆破振動效應(yīng)分析領(lǐng)域具有較好的應(yīng)用前景，為構(gòu)建爆破振動大型網(wǎng)絡(luò)化測試系統(tǒng)與控制平臺提供了適宜的數(shù)據(jù)分析技術(shù)。

［1］ 謝全民，龍 源，田作威.爆破振動信號時頻特征的三維分形特性研究［J］.振動與沖擊，2010，29(12):122－126.

［2］ 謝全民，龍 源，鐘明壽.小波包與分形組合技術(shù)在爆破振動信號分析中的應(yīng)用研究［J］.振動與沖擊，2011，30(1):11－15.

［3］ 婁建武，龍 源.爆破震動信號的特征提取及識別技術(shù)研究［J］.振動與沖擊，2003，22(3):80－83.

［4］ 中國生，徐國元，江文武.基于小波變換的爆破地震信號去噪的應(yīng)用［J］.中南大學(xué)學(xué)報(bào)(自然科學(xué)版)，2006，37(1):155－159.

［5］ Sweldens W. The lifting scheme:a custom－design construction of biorthogonal wavelets［J］.Applied and Computational Harmonic Analysis，1996，15(3):186－200.

［6］ 耿艷峰，馮叔初.小波構(gòu)造綜述［J］.石油大學(xué)學(xué)報(bào)(自然科學(xué)版)，2004，28(1):127－131.

［7］ 粟 鳴，郭東敏，權(quán)建峰.基于提升小波的改進(jìn)半軟閾值降噪方法［J］.探測與控制學(xué)報(bào)，2009，31(4):54－57.

［8］ 段晨東，何正嘉.基于提升模式的特征小波構(gòu)造及其應(yīng)用［J］.振動工程學(xué)報(bào)，2007，20(1):85－90.

［9］ 段晨東，李凌均，何正嘉.第二代小波變換在旋轉(zhuǎn)機(jī)械故障診斷中的應(yīng)用［J］.機(jī)械科學(xué)與技術(shù)，2004，23(2):224－229.

［10］ 周 瑞.基于第二代小波的機(jī)械故障信號處理方法研究［D］.哈爾濱:哈爾濱工業(yè)大學(xué)，2009.

［11］ Xie Q M，Long Y.Blasting vibration signal comparative analysis based on wavelet and wavelet packet technology［C］.The 2th Aisa－Pacific Symposium on Blasting Technology，2009:513－519.

［12］ 劉世元，杜潤生，楊叔子.小波包改進(jìn)算法及其在柴油機(jī)振動診斷中的應(yīng)用［J］.內(nèi)燃機(jī)學(xué)報(bào)，2000，18(1):11－16.

［13］ 曹建軍.基于提升小波包和改進(jìn)蟻群算法的自行火炮在線診斷研究［D］.石家莊:軍械工程學(xué)院，2007.

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

［15］ Sweldens W.The lifting scheme:a construction of second generation wavelets［J］.SIMA J Math Anal，1996，29(2):511－546.

［16］ Donoho D L，Johnstone I M.Adapting to unkown smoothless via wavelet shrinkage［J］.J.Amer.Statist.Assoc.，1995，90(432):1200－1224.