SGWT在爆破振動信號信噪分離中的應(yīng)用研究

謝全民，龍 源，鐘明壽，陸凡東，毛益明

（1.解放軍理工大學(xué) 工程兵工程學(xué)院，南京 210007；2.中國人民解放軍63988部隊，山西晉城 048000）

工程爆破中，減小和控制爆破振動對工程安全和保護(hù)環(huán)境具有重要意義。爆破振動信號分析既是研究爆破振動危害控制的基礎(chǔ)，也是控制爆破振動危害的前提［1-3］。由于爆破振動效應(yīng)狀態(tài)監(jiān)測過程中易受外部環(huán)境和測試系統(tǒng)的影響，測試信號中不可避免地含有噪聲。因此，濾除爆破振動信號中噪聲的影響，提高爆破振動數(shù)據(jù)的可靠性和準(zhǔn)確性，是爆破振動效應(yīng)分析研究的基礎(chǔ)［4］。爆破振動信號是一種短時非平穩(wěn)信號，有效信號和噪聲經(jīng)過小波變換各自系數(shù)將呈現(xiàn)不同的傳播特性，據(jù)此使用小波閾值去噪具有良好的效果［5-6］。工程爆破中為建立爆破振動預(yù)報與控制平臺，要求在傳統(tǒng)小波去噪算法的基礎(chǔ)上提高信號分析處理的效率和質(zhì)量，實現(xiàn)在線處理對實時性、準(zhǔn)確性的要求。

Sweldens于1996年提出的小波提升或者提升算法Lifting Scheme，被稱作二代小波變換（Second Generation Wavelet Transform，SGWT），繼承了經(jīng)典小波多分辨率特性，可以實現(xiàn)原地運算，占用空間小，變換速度快［7-8］。與一代小波構(gòu)造方法的主要區(qū)別在于它不依賴于傅氏變換，小波基函數(shù)不再由某一個函數(shù)的平移和伸縮而產(chǎn)生。提升小波算法可以通過設(shè)計不同的預(yù)測算子和更新算子得到具有某些特殊功能的小波函數(shù)［9-10］，更好地解決實際問題。

Daubechies和Sweldens證明了任何離散小波都可以用提升方案來實現(xiàn)，增強(qiáng)了二代小波的適用性。目前提升算法主要應(yīng)用在機(jī)械振動信號分析，故障診斷等方面［11-14］，SGWT在爆破振動信號分析方面的應(yīng)用還沒有相關(guān)報道。

本文構(gòu)造了基于插值細(xì)分法的提升小波，同提升db小波一起用于實測爆破振動信號的分解-閾值去噪-重構(gòu)，對SGWT在爆破振動信號分析預(yù)處理中的應(yīng)用進(jìn)行了嘗試，取得了滿意效果。

1 爆破振動信號的提升小波變換

設(shè)爆破振動信號采樣序列{s（k），k∈Z}，基于提升算法的小波分解由三步組成：

（1）剖分（split）：將采樣序列{s（k），k∈Z}剖分為奇樣本序列S0={s0（k），k∈Z}和偶樣本序列Se={se（k），k∈Z}，其中：

（2）預(yù)測（predict）：也稱作對偶提升，設(shè)P（·）為預(yù)測器，用se（k）預(yù)測s0（k），定義預(yù)測偏差為細(xì)節(jié)信號d（k）：

（3）更新（update）：也稱作提升，設(shè)U（·）為更新器，在細(xì)節(jié)信號序列D={d（k），k∈Z）的基礎(chǔ)上更新se（k），其結(jié)果定義為逼近信號c（k）：

小波變換的重構(gòu)過程由恢復(fù)更新、恢復(fù)預(yù)測、奇偶采樣序列合并三步構(gòu)成，重構(gòu)公式可由分解式（1）～式（4）變換得到。

圖1 提升小波分解重構(gòu)示意圖Fig.1 Lifting wavelet decomposition and reconstruction

2 廣義小波與提升小波構(gòu)造

Sweldens指出在等間隔采樣和無加權(quán)函數(shù)的情況下，任何一種提升格式都對應(yīng)著一個廣義小波［15］，研究者可以根據(jù)需要，構(gòu)造新的小波。SGWT的關(guān)鍵在于確定預(yù)測系數(shù)和更新系數(shù)，插值細(xì)分法是提升模式構(gòu)造廣義小波的常用方法，由它構(gòu)造的尺度函數(shù)和小波函數(shù)具有良好的緊支撐性，小波變換具有線性相位，可以有效避免變換過程中的相位失真［16］。

2.1 基于插值細(xì)分小波

已知N+1個點x0，x1，…，xN的函數(shù)值為y0，y1，…，yN且yi=f（xi），i=0，…，N，則存在唯一一個次數(shù)不大于N的多項式Ln（x），使Ln（xi）=f（xi），那么：

式中：

每次插值細(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)生新的采樣值處于已知樣本的中間，插值點為：x=xk+（N+1）/2，則預(yù)測系數(shù)可由式（7）確定，即：

N=6.8時，根據(jù)式（7）求得預(yù)測系數(shù)如表1所示。

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

由圖1小波重構(gòu)部分，重構(gòu)關(guān)系簡化為：

根據(jù)上式對S進(jìn)行插值細(xì)分，即利用k-1，k，k+1，k+2 四點的sj，k值預(yù)測sj+1，2k+1，插值邊界采用補(bǔ)零的方法進(jìn)行處理，經(jīng)過20次迭代后生成尺度函數(shù)φ（x）波形見圖2、圖3。

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

其中：

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

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

對s0按照式（9）進(jìn)行插值細(xì)分，經(jīng)過20次迭代可得到相應(yīng)的小波函數(shù)，如圖2、圖3所示。

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

圖3 SGW（8，8）尺度函數(shù)φ（x）與小波函數(shù)ψ（x）Fig.3 Scaling function and wavelet function of SGW（8，8）

基于提升模式的尺度函數(shù)和小波函數(shù)是雙正交的、對稱的、緊支撐的，并且具有沖擊形狀［9，15］。當(dāng)N和較小時，尺度函數(shù)和小波函數(shù)的支撐區(qū)間較?。环粗?，支撐區(qū)間較大，具有較好的連續(xù)性。一般而言，支撐區(qū)間小的小波函數(shù)適合于處理非平穩(wěn)信號，小波系數(shù)能夠有效地描述信號的瞬態(tài)分量，而支撐區(qū)間大且連續(xù)較好的小波適合于描述穩(wěn)態(tài)信號。爆破振動信號具有典型的短時非平穩(wěn)特性［18，19］，本文選用 SGW（6，6）和SGW（8，8）進(jìn)行爆破振動信號分析。

2.2 提升db小波

Daubechies系列小波具有較好的緊支性、光滑性及近似對稱性，并已成功應(yīng)用于包括爆破振動信號在內(nèi)的非平穩(wěn)信號問題［6，19-21］。該小波系列按正整數(shù)N具有不同的序列（dbN），本研究采用Matlab中提供db小波的提升方案（Daubechies lifting schemes）liftingdb6、liftingdb8，同基于插值細(xì)分的二代小波 SGW（6，6）、SGW（8，8）進(jìn)行爆破振動信號去噪。

基于插值細(xì)分的二代小波變換，當(dāng)（N，）給定后，運算量與log2（N-1）/（Nmax-1）成正比，Mallat小波分解運算量與Nlog2N成正比，在同樣數(shù)據(jù)長度下，提升db小波比傳統(tǒng)db小波變換速度提高很多，但和相同支集的插值細(xì)分小波相比，運算量較大，運算速度較慢［5，7］?；诓逯导?xì)分的二代小波變換具有更高的信號處理效率。

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

3 基于提升算法的爆破振動信號去噪

傳統(tǒng)小波去噪中，小波分解是采用信號與小波基進(jìn)行卷積運算，因此小波分解的結(jié)果與所采用小波基的形狀密切相關(guān)，一旦選用不適當(dāng)?shù)男〔ɑ瘮?shù)會沖淡振動信號的局部特征信息，造成原始信號的細(xì)節(jié)信息丟失。為克服上述缺陷，依據(jù)基偶空間的相關(guān)性［22］，利用上述構(gòu)造的插值細(xì)分小波SGW和提升db小波進(jìn)行爆破振動信號的二代小波閾值去噪。

圖4為典型的爆破振動實測信號，信號采樣頻率為1 024 Hz，采樣時間長度取2 000個采樣點，具體的爆破振動信號如圖4所示。振動波形中混雜有毛刺及由于測試系統(tǒng)原因而帶來的方波干擾信號。圖5為采用SGW（6，6）分解三層后繪制的三維時頻能量譜，從圖中可以看出該信號中含有200 Hz～250 Hz頻段內(nèi)的高頻噪聲分量。

圖4 實測爆破振動含噪信號Fig.4 Measured blasting induced vibration signal

圖5 爆破振動含噪信號三維時頻譜Fig.5 3D time-frequency spectrum of blasting induced vibration signal

3.1 提升小波降噪算法

對爆破振動信號小波降噪過程如下：

（1）采用提升算法對信號進(jìn)行多尺度分解；

（2）對分解所得高頻細(xì)節(jié)信號進(jìn)行閾值處理；

（3）將逼近信號和經(jīng)閾值處理后的細(xì)節(jié)信號進(jìn)行重構(gòu)得到降噪后的爆破振動信號；

（4）去噪效果分析。

3.2 提升小波變換

采用Matiab編程，根據(jù)提升小波變換原理及插值細(xì)分小波構(gòu)造算法，分別使用上述構(gòu)造的基于插值細(xì)分的 SGW（6，6）、SGW（8，8）和提升 db6、db8 小波進(jìn)行提升小波分解，分解層數(shù)為3層，如圖6～圖9所示。a3為第3層逼近信號，體現(xiàn)爆破振動信號的低頻分量，d3～d1為細(xì)節(jié)信號，體現(xiàn)爆破振動信號的高頻分量。

圖6 基于SGW（6，6）的爆破振動小波分解Fig.6 Wavelet decomposition of blasting induced vibration based on SGW（6，6）

圖7 基于SGW（8，8）的爆破振動小波分解Fig.7 Wavelet decomposition of blasting induced vibration based on SGW（8，8）

圖8 基于liftingdb6的爆破振動小波分解Fig.8 Wavelet decomposition of blasting induced vibration based on liftingdb6

3.3 閾值去噪

對含噪信號進(jìn)行二代小波變換時，噪聲也會產(chǎn)生高頻系數(shù)，故高頻細(xì)節(jié)分量是有用信號和噪聲高頻系數(shù)的疊加，需要選擇合適的閾值，使噪聲能與有效信號相對應(yīng)的小波系數(shù)合理地區(qū)分，實現(xiàn)爆破振動信號的信噪分離。

Donoho和Johnstone提出軟閾值處理法［23］：

圖9 基于liftingdb8的爆破振動小波分解Fig.9 Wavelet decomposition of blasting induced vibration based on liftingdb8

其中各個尺度的閾值按下式確定：

Tj為各分解尺度對應(yīng)的閾值，j為分解尺度。

爆破振動信號的噪聲方差 σ未知，由式（16）進(jìn)行估計，median（·）為中值函數(shù)。

根據(jù)式（14）～式（16），可在提升小波分解的基礎(chǔ)上對高頻系數(shù)實現(xiàn)閾值去噪。

圖10表示分別采用基于插值細(xì)分的二代小波SGW（6，6）、SGW（8，8）和基于提升變換的 db小波 liftingdb6、liftingdb8經(jīng)過3層提升小波變換之后，由尺度j=3的逼近信號a3和各個尺度閾值處理后的細(xì)節(jié)信號3、2、1重構(gòu)獲得去噪處理后的爆破振動波形。由圖10中二代小波濾波后的振動波形可以看出，已經(jīng)基本消除了爆破振動監(jiān)測信號中由于測試系統(tǒng)和環(huán)境噪聲帶來的干擾。圖10中降噪后的波形曲線相對圖4中含噪的爆破振動波形曲線更加光滑，峰值上升和衰減等振動特征在去噪后的振動信號中表現(xiàn)得更加清晰，得到的有用信息圖10客觀反映了實測爆破振動信號的主要成分。

圖10 降噪后的爆破振動信號Fig.10 De-noised blasting induced vibration signals

圖11為采用SGW（6，6）小波去噪后的三維時頻譜，與圖5中原爆破振動測試信號的三維時頻譜進(jìn)行比較，可以看出，經(jīng)過二代小波去噪后200 Hz～250 Hz頻段內(nèi)的高頻噪聲被基本濾除，爆破振動信號主要能量集中在0 Hz～100 Hz頻段內(nèi)，爆破振動信號能量隨時間、頻率變化而體現(xiàn)的起伏衰減特征可以明顯地得到識別。

圖11 降噪后的爆破振動三維時頻圖Fig.11 3D time-frequency spectrum of de-noised blasting induced vibration signal

圖12 去噪前后相對誤差Fig.12 Relative error after and before de-noising

圖12為采用上述4個二代小波進(jìn)行爆破振動信號去噪前后的相對誤差分析，可以看出提升db小波相對插值細(xì)分小波SGW去噪前后的相對誤差較小。為定量研究SGWT在爆破振動信號去噪分析中的應(yīng)用效果，引入均方根誤差（RMSE）、信噪比（SNR）、峰值誤差（PE）作為評價標(biāo)準(zhǔn)。

式（17）～式（19）中{si}（i=1，2，…，N）為實測爆破振動信號采樣值，{i}（i=1，2，…，N）為經(jīng)過去噪后的重構(gòu)信號，N為采樣點數(shù)。

圖13 降噪效果對比Fig.13 Comparison of de-noising effect

分析圖10中4個提升小波去噪后的振動信號重構(gòu)圖和圖13中降噪效果對比，可以得出：

（1）所采用的基于4個提升小波基的二代小波閾值去噪均能有效濾除爆破振動測試信號中包含的高頻噪聲。對同一爆破振動信號進(jìn)行去噪分析，在選定同樣閾值條件下，liftingdb6小波可以得到最小的RMSE，最高的SNR和最小的PE，去噪效果較其它提升小波改善明顯。對于插值細(xì)分小波系列，SGW（6，6）優(yōu)于SGW（8，8）；對于提升db小波系列，liftingdb6優(yōu)于liftingdb8。

（2）基于提升算法的db小波相對基于插值細(xì)分法構(gòu)造的SGW在去噪過程中獲得相近的RMSE、PE，但提升db相對插值SGW提高了SNR。主要原因在于提升db的小波函數(shù)波形［16］比基于插值細(xì)分法構(gòu)造SGW的小波函數(shù)波形與爆破振動波形具有更強(qiáng)的相似性。

4 結(jié)論

（1）根據(jù)二代小波變換（SGWT）原理，通過插值細(xì)分法構(gòu)造的二代小波SGW（6，6）、SGW（8，8）同提升db小波liftingdb6、liftingdb8用于爆破振動信號信噪分離，研究表明，插值細(xì)分小波SGW和提升db小波均可有效濾除實測爆破振動信號中包含的噪聲，達(dá)到了預(yù)期的研究目的。

（2）小波基的選擇會影響去噪效果，根據(jù)爆破振動信號短時非平穩(wěn)特性的要求，在SGWT閾值去噪過程中，選取提升小波基的支集不宜過長。

（3）提升db小波函數(shù)與爆破振動波形具有更強(qiáng)的相似性，在二代小波閾值去噪過程中能提高信噪比。插值細(xì)分小波相對提升db小波有更高的信號處理效率，因此在滿足信號處理效率和質(zhì)量要求的范圍內(nèi)，采用SGWT對爆破振動信號進(jìn)行分析時，應(yīng)當(dāng)合理利用二者的優(yōu)勢。

［1］張耀平，曹 平，高賽紅.爆破振動信號的小波分解及各頻段的能量分布特征［J］.金屬礦山，2007，377（11）：42-43.

［2］史秀志，薛劍光，陳壽如.爆破振動特征參量的粗糙集模糊神經(jīng)網(wǎng)絡(luò)預(yù)測［J］.振動與沖擊，2009，28（7）：73-76.

［3］晏俊偉，龍 源，方 向.基于小波變換的爆破振動信號能量分布特征分析［J］.爆炸與沖擊，2007，27（5）：405-410.

［4］中國生，徐國元，趙建平.基于小波變換的爆破地震信號閾值去噪的應(yīng)用研究［J］.巖土工程學(xué)報，2005，27（9）：1055-1059.

［5］孔國杰，張培林，曹建軍.基于提升小波變換的信號降噪及其工程應(yīng)用［J］.計算機(jī)工程與應(yīng)用，2008，44（10）：234-237.

［6］謝全民，龍 源，鐘明壽.基于小波、小波包兩種方法的爆破振動信號對比分析［J］.工程爆破，2009，15（1）：5-9.

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

［8］曹建軍，張培林，任國全.提升小波包最優(yōu)基分解算法及在振動信號降噪中的應(yīng)用［J］.振動與沖擊，2008，27（8）：115-116.

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

［10］耿艷峰，馮叔初.小波構(gòu)造綜述［J］.石油大學(xué)學(xué)報（自然科學(xué)版），2004，28（1）：127-131.

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

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

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

［14］荊雙喜，李勝華，華 偉.第二代小波降噪法在故障診斷中的應(yīng)用［J］.煤礦機(jī)電，2007，2：3-5.

［15］ Sweldens W.The lifting scheme：a construction of second generation wavelets［J］.SIMA J Math Anal，1996，29（2）：511-546.

［16］何正嘉.小波有限元理論及其工程應(yīng)用［M］.北京：科學(xué)出版社，2006.

［17］ 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.

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

［19］龍 源，婁建武，徐全軍.小波分析在結(jié)構(gòu)物對爆破振動響應(yīng)的能量分析法中的應(yīng)用［J］.爆破器材，2001，30（3）：1-5.

［20］何 軍，于亞倫，梁文基.爆破振動信號的小波分析［J］.巖土工程學(xué)報，1998，20（1）：47-50.

［21］黃文華，徐全軍.小波變換在判斷爆破地震危害中的應(yīng)用［J］.工程爆破，2001，7（1）：24-27.

［22］劉樹春，潘紫微，宋 淼.第二代小波在振動信號去噪中新方法的研究［J］.機(jī)械傳動，2008，32（3）：64-69.

［23］ Donoho D L，Johnstone I M，Adapting to unkown smoothless via wavelet shrinkage［J］.J.Amer.Statist.Assoc.，1995，90：1200-1224.