謝全民,龍 源,鐘明壽,陸凡東,毛益明
(1.解放軍理工大學 工程兵工程學院,南京 210007;2.中國人民解放軍63988部隊,山西晉城 048000)
工程爆破中,減小和控制爆破振動對工程安全和保護環(huán)境具有重要意義。爆破振動信號分析既是研究爆破振動危害控制的基礎,也是控制爆破振動危害的前提[1-3]。由于爆破振動效應狀態(tài)監(jiān)測過程中易受外部環(huán)境和測試系統(tǒng)的影響,測試信號中不可避免地含有噪聲。因此,濾除爆破振動信號中噪聲的影響,提高爆破振動數(shù)據(jù)的可靠性和準確性,是爆破振動效應分析研究的基礎[4]。爆破振動信號是一種短時非平穩(wěn)信號,有效信號和噪聲經(jīng)過小波變換各自系數(shù)將呈現(xiàn)不同的傳播特性,據(jù)此使用小波閾值去噪具有良好的效果[5-6]。工程爆破中為建立爆破振動預報與控制平臺,要求在傳統(tǒng)小波去噪算法的基礎上提高信號分析處理的效率和質(zhì)量,實現(xiàn)在線處理對實時性、準確性的要求。
Sweldens于1996年提出的小波提升或者提升算法Lifting Scheme,被稱作二代小波變換(Second Generation Wavelet Transform,SGWT),繼承了經(jīng)典小波多分辨率特性,可以實現(xiàn)原地運算,占用空間小,變換速度快[7-8]。與一代小波構造方法的主要區(qū)別在于它不依賴于傅氏變換,小波基函數(shù)不再由某一個函數(shù)的平移和伸縮而產(chǎn)生。提升小波算法可以通過設計不同的預測算子和更新算子得到具有某些特殊功能的小波函數(shù)[9-10],更好地解決實際問題。
Daubechies和Sweldens證明了任何離散小波都可以用提升方案來實現(xiàn),增強了二代小波的適用性。目前提升算法主要應用在機械振動信號分析,故障診斷等方面[11-14],SGWT在爆破振動信號分析方面的應用還沒有相關報道。
本文構造了基于插值細分法的提升小波,同提升db小波一起用于實測爆破振動信號的分解-閾值去噪-重構,對SGWT在爆破振動信號分析預處理中的應用進行了嘗試,取得了滿意效果。
設爆破振動信號采樣序列{s(k),k∈Z},基于提升算法的小波分解由三步組成:
(1)剖分(split):將采樣序列{s(k),k∈Z}剖分為奇樣本序列S0={s0(k),k∈Z}和偶樣本序列Se={se(k),k∈Z},其中:
(2)預測(predict):也稱作對偶提升,設P(·)為預測器,用se(k)預測s0(k),定義預測偏差為細節(jié)信號d(k):
(3)更新(update):也稱作提升,設U(·)為更新器,在細節(jié)信號序列D={d(k),k∈Z)的基礎上更新se(k),其結果定義為逼近信號c(k):
小波變換的重構過程由恢復更新、恢復預測、奇偶采樣序列合并三步構成,重構公式可由分解式(1)~式(4)變換得到。
圖1 提升小波分解重構示意圖Fig.1 Lifting wavelet decomposition and reconstruction
Sweldens指出在等間隔采樣和無加權函數(shù)的情況下,任何一種提升格式都對應著一個廣義小波[15],研究者可以根據(jù)需要,構造新的小波。SGWT的關鍵在于確定預測系數(shù)和更新系數(shù),插值細分法是提升模式構造廣義小波的常用方法,由它構造的尺度函數(shù)和小波函數(shù)具有良好的緊支撐性,小波變換具有線性相位,可以有效避免變換過程中的相位失真[16]。
已知N+1個點x0,x1,…,xN的函數(shù)值為y0,y1,…,yN且yi=f(xi),i=0,…,N,則存在唯一一個次數(shù)不大于N的多項式Ln(x),使Ln(xi)=f(xi),那么:
式中:
每次插值細分,取N個(N=2D,D∈Z+)等間隔樣本yj,k-D+1,…,yj,k,…,yj,k+D,其對應的采樣時刻為xk+1,xk+2,…,xk+N,xk為任意的起始時間,通過細分產(chǎn)生新的采樣值處于已知樣本的中間,插值點為:x=xk+(N+1)/2,則預測系數(shù)可由式(7)確定,即:
N=6.8時,根據(jù)式(7)求得預測系數(shù)如表1所示。
設信號S為一個σ序列,即:
由圖1小波重構部分,重構關系簡化為:
根據(jù)上式對S進行插值細分,即利用k-1,k,k+1,k+2 四點的sj,k值預測sj+1,2k+1,插值邊界采用補零的方法進行處理,經(jīng)過20次迭代后生成尺度函數(shù)φ(x)波形見圖2、圖3。
構造小波函數(shù),需要先設計更新算子U,更新算子的作用是保證小波變換前后信號的低階消失矩不變。當N=(N為預測器的個數(shù),為更新器的個數(shù))時,可以直接將預測系數(shù)除以2作為提升系數(shù)[17]。
其中:
式中,U和D分別為更新算子運算系數(shù)和細節(jié)信號序列。與尺度函數(shù)的構造類似,設細節(jié)信號為σ序列,近似信號S為零序列,進行一次提升格式反變換,信號的偶序列為:
該序列經(jīng)過預測算子P(·)進行恢復預測運算,得到信號的奇序列為:
對s0按照式(9)進行插值細分,經(jīng)過20次迭代可得到相應的小波函數(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]。當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)進行爆破振動信號分析。
Daubechies系列小波具有較好的緊支性、光滑性及近似對稱性,并已成功應用于包括爆破振動信號在內(nèi)的非平穩(wěn)信號問題[6,19-21]。該小波系列按正整數(shù)N具有不同的序列(dbN),本研究采用Matlab中提供db小波的提升方案(Daubechies lifting schemes)liftingdb6、liftingdb8,同基于插值細分的二代小波 SGW(6,6)、SGW(8,8)進行爆破振動信號去噪。
基于插值細分的二代小波變換,當(N,)給定后,運算量與log2(N-1)/(Nmax-1)成正比,Mallat小波分解運算量與Nlog2N成正比,在同樣數(shù)據(jù)長度下,提升db小波比傳統(tǒng)db小波變換速度提高很多,但和相同支集的插值細分小波相比,運算量較大,運算速度較慢[5,7]?;诓逯导毞值亩〔ㄗ儞Q具有更高的信號處理效率。
表1 預測系數(shù)Tab.1 Predict coefficients
傳統(tǒng)小波去噪中,小波分解是采用信號與小波基進行卷積運算,因此小波分解的結果與所采用小波基的形狀密切相關,一旦選用不適當?shù)男〔ɑ瘮?shù)會沖淡振動信號的局部特征信息,造成原始信號的細節(jié)信息丟失。為克服上述缺陷,依據(jù)基偶空間的相關性[22],利用上述構造的插值細分小波SGW和提升db小波進行爆破振動信號的二代小波閾值去噪。
圖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
對爆破振動信號小波降噪過程如下:
(1)采用提升算法對信號進行多尺度分解;
(2)對分解所得高頻細節(jié)信號進行閾值處理;
(3)將逼近信號和經(jīng)閾值處理后的細節(jié)信號進行重構得到降噪后的爆破振動信號;
(4)去噪效果分析。
采用Matiab編程,根據(jù)提升小波變換原理及插值細分小波構造算法,分別使用上述構造的基于插值細分的 SGW(6,6)、SGW(8,8)和提升 db6、db8 小波進行提升小波分解,分解層數(shù)為3層,如圖6~圖9所示。a3為第3層逼近信號,體現(xiàn)爆破振動信號的低頻分量,d3~d1為細節(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
對含噪信號進行二代小波變換時,噪聲也會產(chǎn)生高頻系數(shù),故高頻細節(jié)分量是有用信號和噪聲高頻系數(shù)的疊加,需要選擇合適的閾值,使噪聲能與有效信號相對應的小波系數(shù)合理地區(qū)分,實現(xiàn)爆破振動信號的信噪分離。
Donoho和Johnstone提出軟閾值處理法[23]:
圖9 基于liftingdb8的爆破振動小波分解Fig.9 Wavelet decomposition of blasting induced vibration based on liftingdb8
其中各個尺度的閾值按下式確定:
Tj為各分解尺度對應的閾值,j為分解尺度。
爆破振動信號的噪聲方差 σ未知,由式(16)進行估計,median(·)為中值函數(shù)。
根據(jù)式(14)~式(16),可在提升小波分解的基礎上對高頻系數(shù)實現(xiàn)閾值去噪。
圖10表示分別采用基于插值細分的二代小波SGW(6,6)、SGW(8,8)和基于提升變換的 db小波 liftingdb6、liftingdb8經(jīng)過3層提升小波變換之后,由尺度j=3的逼近信號a3和各個尺度閾值處理后的細節(jié)信號3、2、1重構獲得去噪處理后的爆破振動波形。由圖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ī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個二代小波進行爆破振動信號去噪前后的相對誤差分析,可以看出提升db小波相對插值細分小波SGW去噪前后的相對誤差較小。為定量研究SGWT在爆破振動信號去噪分析中的應用效果,引入均方根誤差(RMSE)、信噪比(SNR)、峰值誤差(PE)作為評價標準。
式(17)~式(19)中{si}(i=1,2,…,N)為實測爆破振動信號采樣值,{i}(i=1,2,…,N)為經(jīng)過去噪后的重構信號,N為采樣點數(shù)。
圖13 降噪效果對比Fig.13 Comparison of de-noising effect
分析圖10中4個提升小波去噪后的振動信號重構圖和圖13中降噪效果對比,可以得出:
(1)所采用的基于4個提升小波基的二代小波閾值去噪均能有效濾除爆破振動測試信號中包含的高頻噪聲。對同一爆破振動信號進行去噪分析,在選定同樣閾值條件下,liftingdb6小波可以得到最小的RMSE,最高的SNR和最小的PE,去噪效果較其它提升小波改善明顯。對于插值細分小波系列,SGW(6,6)優(yōu)于SGW(8,8);對于提升db小波系列,liftingdb6優(yōu)于liftingdb8。
(2)基于提升算法的db小波相對基于插值細分法構造的SGW在去噪過程中獲得相近的RMSE、PE,但提升db相對插值SGW提高了SNR。主要原因在于提升db的小波函數(shù)波形[16]比基于插值細分法構造SGW的小波函數(shù)波形與爆破振動波形具有更強的相似性。
(1)根據(jù)二代小波變換(SGWT)原理,通過插值細分法構造的二代小波SGW(6,6)、SGW(8,8)同提升db小波liftingdb6、liftingdb8用于爆破振動信號信噪分離,研究表明,插值細分小波SGW和提升db小波均可有效濾除實測爆破振動信號中包含的噪聲,達到了預期的研究目的。
(2)小波基的選擇會影響去噪效果,根據(jù)爆破振動信號短時非平穩(wěn)特性的要求,在SGWT閾值去噪過程中,選取提升小波基的支集不宜過長。
(3)提升db小波函數(shù)與爆破振動波形具有更強的相似性,在二代小波閾值去噪過程中能提高信噪比。插值細分小波相對提升db小波有更高的信號處理效率,因此在滿足信號處理效率和質(zhì)量要求的范圍內(nèi),采用SGWT對爆破振動信號進行分析時,應當合理利用二者的優(yōu)勢。
[1]張耀平,曹 平,高賽紅.爆破振動信號的小波分解及各頻段的能量分布特征[J].金屬礦山,2007,377(11):42-43.
[2]史秀志,薛劍光,陳壽如.爆破振動特征參量的粗糙集模糊神經(jīng)網(wǎng)絡預測[J].振動與沖擊,2009,28(7):73-76.
[3]晏俊偉,龍 源,方 向.基于小波變換的爆破振動信號能量分布特征分析[J].爆炸與沖擊,2007,27(5):405-410.
[4]中國生,徐國元,趙建平.基于小波變換的爆破地震信號閾值去噪的應用研究[J].巖土工程學報,2005,27(9):1055-1059.
[5]孔國杰,張培林,曹建軍.基于提升小波變換的信號降噪及其工程應用[J].計算機工程與應用,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)基分解算法及在振動信號降噪中的應用[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]耿艷峰,馮叔初.小波構造綜述[J].石油大學學報(自然科學版),2004,28(1):127-131.
[11]粟 鳴,郭東敏,權建峰.基于提升小波的改進半軟閾值降噪方法[J].探測與控制學報,2009,31(4):54-57.
[12]段晨東,何正嘉.基于提升模式的特征小波構造及其應用[J].振動工程學報,2007,20(1):85-90.
[13]段晨東,李凌均,何正嘉.第二代小波變換在旋轉(zhuǎn)機械故障診斷中的應用[J].機械科學與技術,2004,23(2):224-229.
[14]荊雙喜,李勝華,華 偉.第二代小波降噪法在故障診斷中的應用[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]何正嘉.小波有限元理論及其工程應用[M].北京:科學出版社,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]婁建武,龍 源.爆破震動信號的特征提取及識別技術研究[J].振動與沖擊,2003,22(3):80-82.
[19]龍 源,婁建武,徐全軍.小波分析在結構物對爆破振動響應的能量分析法中的應用[J].爆破器材,2001,30(3):1-5.
[20]何 軍,于亞倫,梁文基.爆破振動信號的小波分析[J].巖土工程學報,1998,20(1):47-50.
[21]黃文華,徐全軍.小波變換在判斷爆破地震危害中的應用[J].工程爆破,2001,7(1):24-27.
[22]劉樹春,潘紫微,宋 淼.第二代小波在振動信號去噪中新方法的研究[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.