基于CEEMD與TQWT組合方法的爆破振動(dòng)信號(hào)精細(xì)化特征提取

楊仁樹， 付曉強(qiáng)， 楊國梁， 陳 駿

(中國礦業(yè)大學(xué)(北京)力學(xué)與建筑工程學(xué)院，北京 100083)

基于CEEMD與TQWT組合方法的爆破振動(dòng)信號(hào)精細(xì)化特征提取

楊仁樹， 付曉強(qiáng)， 楊國梁， 陳 駿

(中國礦業(yè)大學(xué)(北京)力學(xué)與建筑工程學(xué)院，北京 100083)

針對(duì)傳統(tǒng)小波在爆破振動(dòng)信號(hào)特征提取和分析方面的局限性，提出了基于CEEMD和TQWT組合的信號(hào)精細(xì)化特征提取方法。預(yù)先設(shè)定可調(diào)品質(zhì)因子小波TQWT高、低品質(zhì)因子參數(shù)對(duì)CEEMD分解優(yōu)勢分量重組信號(hào)進(jìn)行分解，并引入相對(duì)權(quán)重因子θ，優(yōu)化了分解過程，實(shí)現(xiàn)了爆破振動(dòng)信號(hào)特征的精細(xì)化提取。分析結(jié)果表明：組合方法對(duì)爆破振動(dòng)信號(hào)的分析不依賴于先驗(yàn)小波基的選擇，分解過程實(shí)現(xiàn)了信號(hào)的二次濾波。通過連續(xù)小波多尺度三維譜和時(shí)頻小波脊線對(duì)比，說明組合算法分解得到的最佳分析信號(hào)可真實(shí)反映振動(dòng)信號(hào)的細(xì)節(jié)信息，時(shí)頻分辨率更高。該組合方法抑制了雜波分量對(duì)信號(hào)特征的干擾，可精確地提取復(fù)雜環(huán)境下的爆破振動(dòng)信號(hào)特征信息。

爆破振動(dòng)；總體平均經(jīng)驗(yàn)?zāi)B(tài)分解；可調(diào)品質(zhì)因子小波變換；能量分布；時(shí)頻脊線

工程爆破振動(dòng)本質(zhì)上是一個(gè)非平穩(wěn)的隨機(jī)過程，爆破振動(dòng)信號(hào)經(jīng)過復(fù)雜的巖土介質(zhì)傳播后，往往摻雜著各種頻率成分的干擾波[1]。爆破振動(dòng)波形是爆破振動(dòng)特征最直觀的體現(xiàn)形式，目前針對(duì)爆破振動(dòng)波形信號(hào)，最主要的分析方法有：傅里葉變換(Fourier Transform，F(xiàn)FT)、小波變換(Wavelet Transform，WT)、希爾伯特黃變換(Hilbert Huang Transform，HHT)、奇異值分解(Singular Value Decomposition,SVD)等，但這些方法在實(shí)際工程應(yīng)用中都存在一定的局限性，使其在爆破振動(dòng)信號(hào)特征提取和分析方面存在局限[2]。

近幾年，國內(nèi)外學(xué)者對(duì)爆破振動(dòng)信號(hào)的分析開展了大量的研究。其中，趙明生等[3]利用RSPWVD(Reassigned of Smooth Pseudo-Wigner Ville Distrbution)二次型時(shí)頻分析與小波分析相結(jié)合的方法對(duì)實(shí)測的爆破振動(dòng)信號(hào)進(jìn)行了時(shí)頻分析，獲得了更加直觀、時(shí)頻分辨率更高的信號(hào)時(shí)頻特征；趙國彥等[4]采用頻率切片小波變換(Frequency Slice Wavelet Transform，F(xiàn)SWT)對(duì)比了巖體微震和爆破振動(dòng)兩類信號(hào)的不同頻帶的能量分布特性差異，得到了在高頻區(qū)域，爆破振動(dòng)比巖體微震信號(hào)所占能量比例更大的重要結(jié)論；龔敏等[5]以重慶渝中隧道為例，采用HHT方法研究了瞬時(shí)能量法和EMD(Empirical Mode Decomposition)法識(shí)別雷管延時(shí)時(shí)間效果，總結(jié)了兩種識(shí)別方法的適用條件和精度。本文中將以HHT分析方法中的CEEMD分解為基礎(chǔ)，結(jié)合可調(diào)品質(zhì)小波變換TQWT，實(shí)現(xiàn)了爆破振動(dòng)信號(hào)的精細(xì)化特征提取，為深入研究爆破振動(dòng)傳播特性和地震效應(yīng)下的振動(dòng)響應(yīng)提供了理論依據(jù)。

1 CEEMD與TQWT基本理論

1.1 CEEMD法[6-8]

YEH等提出了一種補(bǔ)充的總體平均經(jīng)驗(yàn)?zāi)B(tài)分解方法CEEMD,其算法核心是通過向原始信號(hào)中添加兩對(duì)相反的白噪聲信號(hào)分別進(jìn)行EMD分解，將分解的結(jié)果進(jìn)行組合得到最終的IMF(Intrinsic Mode Function)。CEEMD在保證分解效果與總體平均經(jīng)驗(yàn)?zāi)B(tài)分解(Ensemble Empirical Mode Decomposition,EEMD)相當(dāng)?shù)那闆r下，抑制了由白噪聲引起的重構(gòu)誤差。

以EMD分解為基礎(chǔ)，由以下幾個(gè)步驟組成：

步驟1 向原始信號(hào)中加入n組正、負(fù)成對(duì)的輔助白噪聲，從而生成兩套集合IMF：

(1)

式中:S為原信號(hào)；N為輔助噪聲；M1，M2分別為加入正負(fù)成對(duì)噪聲后的信號(hào)，這樣得到集合信號(hào)的個(gè)數(shù)為2n。

步驟2 對(duì)集合中的每一個(gè)信號(hào)做EMD分解，每個(gè)信號(hào)得到一組IMF分量，其中第i個(gè)信號(hào)的第j個(gè)IMF分量表示為cij。

步驟3 通過多組分量組合的方式得到分解結(jié)果：

(2)

式中：cj為CEEMD分解最終得到的第j個(gè)IMF分量。

1.2 可調(diào)品質(zhì)因子小波變換(TQWT)[9-10]

2011年，SELESNICK從頻域?yàn)V波角度出發(fā)設(shè)計(jì)出了可調(diào)品質(zhì)因子小波變換(Tunable Q-factor Wavelet Transform,TQWT)，其利用帶通濾波器組以迭代方式實(shí)現(xiàn)信號(hào)的分解重構(gòu)。以3層分解重構(gòu)為例，對(duì)應(yīng)的雙通道濾波器組,如圖1所示[11]。

圖1 雙通道濾波器組Fig.1 Dual channel filter banks

與傳統(tǒng)恒定品質(zhì)因子(品質(zhì)因子定義為中心頻率與帶寬的比值)小波變換相比，TQWT的最大優(yōu)勢就是其品質(zhì)因子可以進(jìn)行自由調(diào)節(jié)，更具靈活性。設(shè)置TQWT的品質(zhì)因子Q和冗余因子r，利用可調(diào)品質(zhì)因子小波對(duì)原信號(hào)進(jìn)行處理，分解層數(shù)j為理論允許的最大值，計(jì)算公式如下：

(3)

式中：Ns為信號(hào)長度；[?]表示向下取整。

設(shè)信號(hào)x由高共振屬性信號(hào)x1及低共振屬性信號(hào)x2組成：

x=x1+x2

(4)

信號(hào)x必定為非線性信號(hào)，故信號(hào)x1與x2的分離不能基于頻率濾波方法，應(yīng)分別構(gòu)建高、低兩種品質(zhì)因子小波基(分別用TQWT1、TQWT2表示)對(duì)信號(hào)進(jìn)行分解，此方法為基于分形主成分分析方法。對(duì)信號(hào)x分離可轉(zhuǎn)化為約束最優(yōu)化問題：

(5)

式中：w1,j為TQWTi(i=1,2)子帶。計(jì)算獲得w1，w2，分離出的信號(hào)x1，x2可表示為

(6)

振動(dòng)信號(hào)中，波形除了介質(zhì)響應(yīng)頻率、瞬態(tài)沖擊信號(hào)外，亦包含有雜頻環(huán)境噪聲，記為N，分析信號(hào)y可表示為：

y=x1+x2+N

(7)

對(duì)信號(hào)y的分離可將式(6)轉(zhuǎn)化為：

(8)

式中：Φ1，Φ2為高、低品質(zhì)因子逆小波變換；λ1，λ2為正則化參數(shù)，據(jù)噪聲能量高低選取。

從前述信號(hào)的分解可知，信號(hào)正則化分解也可能存在高、低品質(zhì)成分之間，其中某個(gè)品質(zhì)成分分配的能量太多，而剩余的另一品質(zhì)成分分配的能量不足的問題。為了優(yōu)化分配過程，在高、低品質(zhì)兩成分中引入相對(duì)權(quán)重因子θ(0<θ<1)，于是，信號(hào)y分離優(yōu)化問題轉(zhuǎn)化為：

(9)

通過tqwt-mca-eq函數(shù)對(duì)x1，x1進(jìn)行優(yōu)化，得到優(yōu)化后的w1，w2,記為w1s，w2s，w1,js為TQWTi(i=1,2)優(yōu)化子帶，運(yùn)用分裂增廣拉格朗日收縮算法(Split Augmented Lagrangian Shrinkage Algorithm,SALSA)進(jìn)行迭代計(jì)算，分離出優(yōu)化后高共振分量和低共振分量分別為：

(10)

2 爆破振動(dòng)信號(hào)精細(xì)化特征提取

2.1 信號(hào)CEEMD法分解重構(gòu)

圖2(a)為某一實(shí)測典型微差爆破振動(dòng)信號(hào)速度-時(shí)程曲線。其主振頻率為64.697 Hz，最大振速為1.919 cm/s,采樣頻率10 kHz，采樣時(shí)間長度為0.700 3 s，具有爆破振動(dòng)信號(hào)隨機(jī)和非線性的典型特征。

圖2 爆破振動(dòng)信號(hào)時(shí)程波形曲線及頻譜Fig. 2 Time history waveform curve and spectrum of blasting vibration sub signal

從信號(hào)頻譜圖2(b)中可以看出微差爆破振動(dòng)信號(hào)的頻譜極為豐富，對(duì)于多段別微差爆破，爆破振動(dòng)信號(hào)在頻譜上體現(xiàn)出多振型多峰值的振動(dòng)特性。該爆破振動(dòng)信號(hào)頻率主要集中在250 Hz以內(nèi)，同時(shí)在300～400 Hz范圍內(nèi)包含明顯的強(qiáng)噪聲，強(qiáng)噪聲的存在對(duì)信號(hào)有效信息的精確提取會(huì)造成很大干擾。

將圖2(a)中爆破振動(dòng)信號(hào)按照1.1節(jié)所述步驟進(jìn)行分解，得到12個(gè)IMF分量和一個(gè)周期趨于無窮的趨勢項(xiàng)r。圖3為圖2(a)原始信號(hào)及其CEEMD分解得到的IMF分量?？梢钥闯觯篊EEMD分解和EMD分解一樣具有二進(jìn)濾波特征，在依次排列的IMF分量中，前面幾個(gè)也是高頻分量，通常隨機(jī)噪聲包含其中。

圖3 原信號(hào)CEEMD分解結(jié)果Fig.3 The CEEMD decomposition results of original signal

在CEEMD分解中IMF=ceemdan(s,Nstd,NR,MaxIter),其中s為分析信號(hào)；Nstd指的是噪聲標(biāo)準(zhǔn)偏差，一般取值為：0.2～0.3，此處取為0.2；NR為分解實(shí)現(xiàn)次數(shù)，此處取為5，MaxIter是篩分迭代計(jì)算允許的最大次數(shù)，取值100。圖4為5次篩分各模態(tài)迭代次數(shù)分布。

圖4 各模態(tài)篩分迭代次數(shù)分布Fig.4 The iteration number distribution of each mode sifting

圖4表明：分析信號(hào)的IMF分量越復(fù)雜，篩分次數(shù)也越多，包含的信息也更為豐富，與原信號(hào)的相似度也就越高。為了確定信號(hào)的優(yōu)勢分量，采用互相關(guān)性指標(biāo)來客觀評(píng)價(jià)IMF分量與原信號(hào)的相關(guān)度。表1為圖2(a)中爆破振動(dòng)信號(hào)CEEMD分解得到的12個(gè)IMF分量(r分量除外)與原信號(hào)的互相關(guān)系數(shù)。

表1 各IMF分量與原信號(hào)的互相關(guān)性

從表1可知，IMF3～I(xiàn)MF7分量與原信號(hào)的互相關(guān)系數(shù)遠(yuǎn)大于其余高階和低階IMF分量的互相關(guān)系數(shù)，可確定為信號(hào)的優(yōu)勢分量。因此，這里選擇這幾階分量作為優(yōu)勢IMF分量，將其合成得到爆破振動(dòng)特征信號(hào)。圖5為信號(hào)IMF3～I(xiàn)MF7階分量重構(gòu)信號(hào)及分離出的自相似性高頻噪聲?？梢钥闯觯簝?yōu)勢IMF分量重構(gòu)的信號(hào)幅值顯著減小，但多段別雷管微差爆破振動(dòng)波形態(tài)清晰出現(xiàn)。去噪后的爆破振動(dòng)信號(hào)相對(duì)光滑，且局部奇異性更易辨識(shí)。

圖5 IMF3～I(xiàn)MF7分量合成的特征信號(hào)與分離噪聲Fig.5 Characteristic of IMF3～I(xiàn)MF7 component synthesis signal and separation noise

信號(hào)經(jīng)過CEEMD法優(yōu)勢分量重構(gòu)后，高頻噪聲分量得到了有效的抑制，信號(hào)的辨識(shí)度進(jìn)一步得到提高，有利于對(duì)信號(hào)的后續(xù)分析。

2.2 重構(gòu)信號(hào)TQWT分解與二次合成

2.2.1 重構(gòu)信號(hào)高、低品質(zhì)TQWT分解

TQWT可以不依賴基函數(shù)而根據(jù)實(shí)際需要設(shè)定其品質(zhì)因子，通過確定品質(zhì)因子Q和冗余因子r，使生成的小波基函數(shù)庫中的所有小波都具有相同的品質(zhì)因子，因此TQWT本質(zhì)上是一種具有一定冗余度的恒定品質(zhì)因子小波變換[12]。

可調(diào)品質(zhì)因子小波變換分析中品質(zhì)因子Q可預(yù)先設(shè)定，冗余因子r代表TQWT變換信號(hào)的過采樣率(信號(hào)分解小波系數(shù)總和nt/信號(hào)長度l)，通常取r=3(r>1)。圖2(a)中爆破振動(dòng)信號(hào)采樣頻率fs(每秒采樣點(diǎn)數(shù))為10 kHz，因此其Nyquist頻率為fs/2=5 000。信號(hào)經(jīng)j層可調(diào)品質(zhì)因子小波變換處理后，最終得到j(luò)+1個(gè)子帶信號(hào)，隨分解層數(shù)j的增加，中心頻率隨之降低，帶寬亦隨之減小。分別表示為：

(11)

(12)

式中：fs為采樣頻率；j為分解層數(shù)。

采用高品質(zhì)因子Q=4，冗余因子r=3，爆破信號(hào)數(shù)據(jù)長度為7 004，分解層數(shù)j由式(3)確定，計(jì)算得j=40，最終得到41個(gè)子帶信號(hào)。因重構(gòu)信號(hào)高品質(zhì)因子分解時(shí)，前13個(gè)子帶能量趨于0，這里僅輸出14～41子帶。信號(hào)TQWT高品質(zhì)因子分解頻率響應(yīng)見圖6(a)，分解子帶小波時(shí)域波形見圖6(b)，頻率能量百分比分布見圖6(c)；采用低品質(zhì)因子Q=1，冗余因子r=3，由式(3)確定，計(jì)算得j=16時(shí)，最終得到17個(gè)子帶信號(hào)。低品質(zhì)因子小波變換頻率響應(yīng)見圖7(a)，分解子帶小波時(shí)域波形見圖7(b)，頻率能量百分比分布見圖7(c)。

圖6(a)和圖7(a)中頻率響應(yīng)為一組非恒定帶寬濾波器組，且相鄰頻帶不正交；圖6(b)和圖7(b)中分解子帶小波時(shí)域波形表明：隨著爆破振動(dòng)信號(hào)高、低品質(zhì)因子小波分解層數(shù)的增加，子帶波形振動(dòng)時(shí)間隨之變長。對(duì)比圖6(c)和圖7(c)中爆破振動(dòng)信號(hào)高、低品質(zhì)因子分解頻率能量百分比分布可以發(fā)現(xiàn)：高品質(zhì)因子分解得到的響應(yīng)頻帶比低品質(zhì)因子分解得到的響應(yīng)頻帶更集中，頻率響應(yīng)范圍更緊湊。

2.2.2 重構(gòu)信號(hào)優(yōu)化TQWT分解

圖6 高品質(zhì)因子TQWT分解(Q=4，r=3，j=40)Fig.6 High Q factor wavelet transform (Q=4，r=3，j=40)

圖7 低品質(zhì)因子TQWT分解(Q=1,r=3,j=16)Fig.7 Low Q factor wavelet transform (Q=1,r=3,j=16)

優(yōu)化后的高、低品質(zhì)因子成分各自的子帶能量百分比分布狀態(tài)見圖8(b)和8(c)，通過SALA收縮算法迭代計(jì)算后的爆破振動(dòng)信號(hào)高、低品質(zhì)因子成分能量集中在有限的子帶區(qū)域內(nèi)，能量的表征更為清晰。優(yōu)化后TQWT分解得到的高品質(zhì)因子成分信號(hào)各子帶時(shí)域波形中：1～4子帶小波系數(shù)為4 096；5～9子帶為2 048；10～13子帶為1 024；14～18子帶為512；19～23子帶為256；24～28子帶為128；29～33子帶為64；34～37子帶為32；38～40子帶為16；41子帶為32，所有子帶小波系數(shù)總和nt為35 728。信號(hào)數(shù)據(jù)長度l為7 004，計(jì)算數(shù)據(jù)冗余參數(shù)r=35 616/7 004=5.10，大于設(shè)定的r初值為3，證明了信號(hào)分離的有效性。

圖8 重構(gòu)信號(hào)優(yōu)化分解及相應(yīng)的高、低品質(zhì)因子成分各子帶能量百分比分布Fig.8 Optimal decomposition of the reconstructed signal and energy percentage distribution of each sub band of the high and low quality factor components

2.3 最佳分析信號(hào)獲取與雷管延時(shí)識(shí)別

圖9 兩組子帶與合成最佳分析信號(hào)(Q=4，r=3，j=40)Fig.9 Two groups of sub bands and synthetic optimal analysis signals (Q=4，r=3，j=40)

文獻(xiàn)[13]通過對(duì)實(shí)測爆破振動(dòng)信號(hào)EMD分解的主分量進(jìn)行Hilbert變換，提取其包絡(luò)曲線準(zhǔn)確識(shí)別了隧道微差爆破雷管實(shí)際延時(shí)時(shí)間。同樣，按照該方法對(duì)最佳分析信號(hào)進(jìn)行希爾伯特包絡(luò)解調(diào)分析，得到的包絡(luò)譜,如圖10所示。

圖10 最佳分析信號(hào)Hilbert包絡(luò)譜Fig.10 The optimal analysis signal Hilbert envelope spectrum

根據(jù)現(xiàn)場爆破炮孔布置形式并結(jié)合雷管實(shí)驗(yàn)數(shù)據(jù)，取包絡(luò)譜中各段最大波峰所處時(shí)刻作為該段別雷管實(shí)際延遲時(shí)間，則各段雷管起爆時(shí)間分別為：0.008 3 s，0.046 1 s，0.115 6 s，0.223 8 s，0.328 5 s，0.494 8 s。說明圖2(a)所示微差爆破振動(dòng)信號(hào)是由6段爆破振動(dòng)波形疊加后形成的，通過與雷管廠家提供的實(shí)驗(yàn)數(shù)據(jù)比較，起爆時(shí)間點(diǎn)分別對(duì)應(yīng)雷管1,3,5,7,9,11段，這與實(shí)際施工設(shè)計(jì)的爆破網(wǎng)絡(luò)中規(guī)定使用的雷管段別是一致的。微差爆破延遲時(shí)間間隔可定義為前后相鄰兩段雷管起爆時(shí)刻間的時(shí)差間隔。這里，為了識(shí)別的方便，將包絡(luò)譜中第一個(gè)峰值點(diǎn)出現(xiàn)的時(shí)間位置定為最低段別雷管(MS1)的起爆時(shí)刻，則該識(shí)別方法得到的各段別雷管實(shí)際起爆時(shí)刻分別為：37.8 ms，69.5 ms，108.2 ms，104.7 ms，166.3 ms，延時(shí)間隔位于雷管實(shí)驗(yàn)起爆延時(shí)誤差區(qū)間范圍內(nèi)。同時(shí)，選用db8小波基對(duì)圖2(a)中原爆破信號(hào)在尺度a為16時(shí)進(jìn)行連續(xù)小波變換取模極大值結(jié)果,如圖11所示。通過對(duì)比可以發(fā)現(xiàn)：原信號(hào)小波變換模極大值方法對(duì)于高段別雷管微差識(shí)別出現(xiàn)“多峰值”現(xiàn)象，識(shí)別精度差，并且過度依賴小波基和尺度的選擇，不具有自適應(yīng)性。

圖11 原爆破信號(hào)連續(xù)小波變換模值(a=16)Fig.11 Continuous wavelet transform modulus maxima of origin Blasting vibration signal(a=16)

上述分析也說明在該次微差爆破中，雷管未發(fā)生跳段、串段等安全隱患，微差延時(shí)識(shí)別一定程度上反映了雷管網(wǎng)絡(luò)的可靠度，為后續(xù)爆破參數(shù)優(yōu)化提供理論支撐。

3 信號(hào)精細(xì)化提取結(jié)果對(duì)比分析

3.1 小波多尺度對(duì)比分析

連續(xù)小波變換是信號(hào)處理中常用的一類積分變換，通過小波基函數(shù)的平移和伸縮操作，可以實(shí)現(xiàn)信號(hào)的多尺度分析。連續(xù)小波變換是將傅里葉積分變換中的積分核變成了小波基函數(shù)，這個(gè)小波基函數(shù)是在標(biāo)準(zhǔn)基函數(shù)的基礎(chǔ)上平移和伸縮得到的，連續(xù)小波變換的定義為[14]：

(13)

式中：ψ*(t;a,b)為小波基函數(shù)；a稱為伸縮因子；b稱為平移因子。信號(hào)連續(xù)小波變換多尺度分析常用的小波基函數(shù)包括三類：Morlet小波、Gauss小波、MexHat小波[15]。這里選用Morlet小波基函數(shù)進(jìn)行分析，并指定波數(shù)k=5對(duì)圖2(a)中原始信號(hào)和圖9中通過組合算法得到的最佳分析信號(hào)同時(shí)進(jìn)行連續(xù)小波變換，獲得兩個(gè)信號(hào)的小波系數(shù)矩陣，該矩陣的各行分別對(duì)應(yīng)一個(gè)伸縮尺度，各列為不同位置的頻率信息，將該矩陣用三維形式,如圖12所示。

圖12 原信號(hào)與最佳分析信號(hào)連續(xù)小波多尺度三維譜 Fig.12 The original signal and the optimal analysis signal wavelet continuous multi-scale 3D spectrum

從圖12可知:通過連續(xù)小波變換多尺度分析，可以得到爆破振動(dòng)信號(hào)在任意尺度任意位置的頻率特征，具有自動(dòng)適應(yīng)信號(hào)時(shí)頻特性的優(yōu)點(diǎn)[16]?？梢跃劢沟饺我獬叨?，從而獲取信號(hào)的細(xì)節(jié)，這對(duì)于信號(hào)分析，尤其是爆破振動(dòng)信號(hào)這類非穩(wěn)態(tài)信號(hào)是非常強(qiáng)大的一個(gè)功能。通過對(duì)比發(fā)現(xiàn)，原信號(hào)含有的強(qiáng)噪聲不規(guī)則地分布在信號(hào)不同尺度上，其小波多尺度譜受Heisenberg測不準(zhǔn)原理的極限制約，在信號(hào)小尺度(高頻段)的頻域集中性差，在信號(hào)的大尺度(低頻段)存在大量的交叉干擾項(xiàng)，這對(duì)信號(hào)多尺度譜的分析帶來了一定的影響。而組合算法二次合成得到的最佳分析信號(hào)可以有效濾除高頻噪聲污染(尺度2～4區(qū)間)和低頻交叉項(xiàng)干擾(尺度8～10區(qū)間)，其尺度譜將信號(hào)能量集中到所關(guān)心的有限的尺度上，更加清晰地判別爆破振動(dòng)信號(hào)的突變信息，使信號(hào)的奇異點(diǎn)位置得以凸現(xiàn)，并且消除了交叉項(xiàng)影響。最佳分析信號(hào)很好地表現(xiàn)了爆破振動(dòng)信號(hào)中一些能量相對(duì)較小的分量，提高了信號(hào)多尺度譜的集中性，更好地反映了爆破振動(dòng)信號(hào)的內(nèi)在特征。

3.2 小波尺度譜時(shí)頻脊提取對(duì)比分析

小波尺度譜時(shí)頻脊線是小波尺度譜時(shí)頻面上不同時(shí)刻信號(hào)小波系數(shù)的模取極大值點(diǎn)的分布，這些點(diǎn)稱為小波脊點(diǎn)[16]。小波脊點(diǎn)分布反映了被分析信號(hào)所包含的瞬時(shí)頻率和瞬時(shí)幅值信息[17]。關(guān)于小波尺度譜時(shí)頻脊線的基礎(chǔ)理論詳見文獻(xiàn)[17]，對(duì)圖2(a)中原信號(hào)和圖9中最佳分析信號(hào)的小波尺度譜時(shí)頻脊提取,如圖13所示?？v坐標(biāo)為歸一化頻率幅值。

圖13 原信號(hào)與最佳分析信號(hào)小波尺度譜時(shí)頻脊線Fig.13 The original signal and the optimal analysis signal wavelet scale spectrum time-frequency ridge line

圖13(a)較為準(zhǔn)確地體現(xiàn)了爆破振動(dòng)信號(hào)的時(shí)頻特征，但受到嚴(yán)重的噪聲干擾影響，在高頻區(qū)域的頻率分辨率很差。而圖13(b)中最佳分析信號(hào)的小波尺度譜時(shí)頻脊分布離散度較小，信號(hào)能量主要集中在0.05 Hz歸一化頻率(250 Hz)范圍內(nèi)，進(jìn)一步提高了其頻域聚集性，很好地保持了時(shí)頻分辨率，且抑制了強(qiáng)噪聲對(duì)爆破信號(hào)時(shí)頻特征的污染和弱噪聲對(duì)時(shí)頻特征的干擾，能更好地提取低頻區(qū)域的調(diào)頻分量信息。組合方法得到的最佳分析信號(hào)的小波尺度譜時(shí)頻脊分布驗(yàn)證了組合分析方法在信號(hào)時(shí)頻特征表征方面的優(yōu)越性，更有助于提高爆破振動(dòng)波數(shù)據(jù)解釋的準(zhǔn)確性。

4 結(jié) 論

(1)可調(diào)品質(zhì)因子小波變換(TQWT)克服了傳統(tǒng)小波依賴先驗(yàn)小波基選擇的局限性，通過確定高、低品質(zhì)因子Q值對(duì)爆破振動(dòng)信號(hào)自適應(yīng)分解，并引入相對(duì)權(quán)重因子θ可實(shí)現(xiàn)分解過程的優(yōu)化，對(duì)優(yōu)化后的高品質(zhì)因子成分(Q=4，r=3，j=40)關(guān)鍵子帶進(jìn)行時(shí)域疊加可獲取原信號(hào)的最佳分析信號(hào)。

(2)對(duì)最佳分析信號(hào)進(jìn)行Hilbert包絡(luò)解調(diào)分析，得到的包絡(luò)譜峰值點(diǎn)對(duì)應(yīng)于多段微差爆破雷管實(shí)際起爆時(shí)刻，據(jù)此可得到各段別雷管實(shí)際延時(shí)時(shí)間，為爆破效果評(píng)價(jià)和參數(shù)優(yōu)化提供了新的思路。

(3)連續(xù)小波變換多尺度分析，可以聚焦到任意尺度，從而獲取信號(hào)的細(xì)節(jié)，更加清晰判別爆破振動(dòng)信號(hào)的突變信息，使信號(hào)的奇異點(diǎn)位置得以凸現(xiàn)，并且消除了交叉項(xiàng)影響。而小波尺度時(shí)頻脊線提取清晰界定了爆破振動(dòng)頻率和能量的分布范圍，驗(yàn)證了組合方法在信號(hào)精細(xì)化分析方面的優(yōu)越性。

[ 1 ] 路亮,龍?jiān)?郭濤,等. 爆破地震波作用下城市隧道結(jié)構(gòu)動(dòng)力響應(yīng)的敏感性[J]. 爆炸與沖擊,2014,34(6):701-708. LU Liang,LONG Yuan,GUO Tao,et al.Dynamic response sensitivity of urban tunnel structures under blasting seismic waves to parameters [J].Explosion and Shock Waves,2014,34(6):701-708.

[ 2 ] 李舜酩,郭海東,李殿榮. 振動(dòng)信號(hào)處理方法綜述[J]. 儀器儀表學(xué)報(bào),2013,34(8):1907-1915. LI Shunming,GUO Haidong,LI Dianrong.Review of vibration signal processing methods [J]. Chinese Journal of Scientific Instrument,2013,34(8):1907-1915.

[ 3 ] 趙明生,張建華,易長平. 基于小波分解的爆破振動(dòng)信號(hào)RSPWVD二次型時(shí)頻分析[J]. 振動(dòng)與沖擊,2011,30(2):44-47. ZHAO Mingsheng,ZHANG Jianhua,YI Changping.Blasting vibration signal RSPWVD quadratic time-frequency analysis based on wavelet decomposition[J]. Journal of Vibration and Shock, 2011, 30(2): 44-47.

[ 4 ] 趙國彥,鄧青林,馬舉. 基于FSWT時(shí)頻分析的礦山微震信號(hào)分析與識(shí)別 [J]. 巖土工程學(xué)報(bào),2015,37(2):306-312. ZHAO Guoyan,DENG Qinglin,MA Ju.Recognition of mine microseismic signals based on FSWT time-frequency analysis [J].Chinese Journal of Geotechnical Engineering，2015,37(2)：306-312.

[ 5 ] 龔敏,邱燚可可,孟祥棟,等. 基于HHT的雷管實(shí)際延時(shí)識(shí)別法在城市環(huán)境微差爆破中的應(yīng)用[J]. 振動(dòng)與沖擊,2015,34(10):206-212. GONG Min,QIU Yikeke,MENG Xiangdong,et al.Identification method of detonator’s actual firing time delay based on HHT and its application in millisecond blasting under urban environment [J].Journal of Vibration and Shock, 2015, 34(10): 206-212.

[ 6 ] YEH J R,SHIEH J S.Complementary ensemble empirical mode decomposition: a novel noise enhanced data analysis method[J].Advances in Adaptive Data Analysis，2010,2 (2):135-156．

[ 7 ] 周濤濤,朱顯明,彭偉才,等. 基于CEEMD和排列熵的故障數(shù)據(jù)小波閾值降噪方法[J]. 振動(dòng)與沖擊,2015,34(23):207-211. ZHOU Taotao,ZHU Xianming,PENG Weicai,et al.A wavelet threshold denoising method for fault data based on CEEMD and permutation entropy [J]. Journal of Vibration and Shock, 2015, 34(23): 207-211.

[ 8 ] LEI Y G，HE Z J，ZI Y Y.Application of the EEMD method to rotor fault diagnosis of rotating machinery [J]. Mechanical Systems and Signal Processing,2009，23(4) : 1327-1338．

[ 9 ] PACHORI R B, NISHAD A.Cross-terms reduction in the Wigner-Ville distribution using tunable-Q wavelet transform [J].Signal Processing, 2016,120:288-304.

[10] PATIDAR S, PACHORI R B, GARG N.Automated diagnosis of coronary artery disease using tunable-Q wavelet transform applied on heart rate signals [J]. Expert Systems with Applications, 2015, 42:3315-3326.

[11] 王宏超,陳進(jìn),董廣明,等. 可調(diào)品質(zhì)因子小波變換在轉(zhuǎn)子早期碰摩故障診斷中應(yīng)用[J]. 振動(dòng)與沖擊,2014,33(10):77-80. WANG Hongchao,CHEN Jin,DONG Guangming,et al.Early rub-impact diagnosis of rotors based on tunable Q-factor wavelet transformation [J].Journal of Vibration and Shock, 2014, 33(10): 77-80.

[12] 唐貴基,王曉龍. 可調(diào)品質(zhì)因子小波變換在滾動(dòng)軸承微弱故障特征提取中的應(yīng)用[J]. 中國電機(jī)工程學(xué)報(bào),2016,36(3):746-754. TANG Guiji,WANG Xiaolong.Application of tunable q-factor wavelet transform to feature extraction of weak fault for rolling bearing [J].Proceedings of the CSEE,2016,36(3):746-754.

[13] 龔敏,邱燚可可,李永強(qiáng),等. 定制雷管微差時(shí)間實(shí)測與識(shí)別法在城市隧道爆破設(shè)計(jì)中的應(yīng)用[J]. 巖石力學(xué)與工程學(xué)報(bào),2015,34(6):1179-1187. GONG Min,QIU Yikeke,LI Yongqiang,et al.Application of measuring and identifying detonation time of customized detonator in blasting design of urban tunnels [J].Chinese Journal of Rock Mechanics and Engineering,2015,34(6):1179-1187.

[14] 趙學(xué)智,陳統(tǒng)堅(jiān),葉邦彥. 奇異值分解對(duì)連續(xù)Morlet小波變換的壓縮和提純[J]. 機(jī)械工程學(xué)報(bào),2015,51(16):57-70. ZHAO Xuezhi,CHEN Tongjian,YE Bangyan.Purification and compression of continuous morlet wavelet transform based on singular value decomposition [J].Journal of Mechanical Engineering, 2015, 51(16):57-70.

[15] 龔海健,黃偉國,趙凱,等. 基于Wigner-Ville分布與小波尺度譜融合的時(shí)頻特征提取方法[J]. 振動(dòng)與沖擊,2011,30(12):35-38. GONG Haijian,HUANG Weiguo,ZHAO Kai,et al.Time-frequency feature extraction based on fusion of Wigner-Ville distribution and wavelet scalogram[J].Journal of Vibration and Shock, 2011, 30(12):35-38.

[16] 路亮,龍?jiān)?謝全民,等. 爆破振動(dòng)信號(hào)的提升小波包分解及能量分布特征[J]. 爆炸與沖擊,2013,33(2):140-147. LU Liang,LONG Yuan,XIE Quanmin,et al.Decomposition and energy distribution of blasting vibration signal based on second generation wavelet packet (SGWP) [J].Explosion and Shock Waves,2013,33(2):140-147.

[17] 李宏坤,徐福健,高巧紅,等. 利用小波尺度譜同步平均的時(shí)頻脊故障特征提取[J]. 振動(dòng)工程學(xué)報(bào),2015,28(3):487-494. LI Hongkun,XU Fujian,GAO Qiaohong,et al . Fault feature extraction for synchronous averaging wavelet scalogram based on time-frequency ridge [J].Journal of Vibration Engineering,2015,28(3):487-494.

Precise feature extraction of blasting vibration signals based on combined method of CEEMD and TQWT

YANG Renshu, FU Xiaoqiang, YANG Guoliang, CHEN Jun

(School of Mechanics and Civil Engineering, China University of Mining and Technology(Beijing), Beijing 100083, China)

Aiming at the limitations of traditional wavelet transformation in feature extraction and analysis of blasting vibration signals, a combination method of Tunable Q-factor wavelet Transform(TQWT) and Complementary Ensemble Empirical Mode Decomposition(CEEMD) signal refine feature extraction was proposed. Presetting tunable wavelet TQWT high and low quality factors’ parameters, the reconstructed signals obtained with dominant components of CEEMD were decomposed. The relative weighted factorθwas introduced to optimize the decomposition process. The fine feature extraction of blasting vibration signals was realized. Analysis results showed that the combination method is not dependent on choices of priori wavelet basis to analyze blasting vibration signals, two filterings of signals are realized in the decomposition process; through comparing continuous wavelet multi-scale 3D spectrum and time frequency wavelet ridge lines, the optimal analysis signals obtained with the decomposition of the combined algorithm can truly reflect the details of vibration signals, the time-frequency resolution rate is higher; the combined method suppresses the interference of clutter wave components to signals’ characteristics, it can extract accurately the characteristics information of blasting vibration signals under complex environment.

blasting signal; CEEMD; TQWT; energy distribution; time-frequency ridge line

國家自然科學(xué)基金面上資助(51274203);高等學(xué)校學(xué)科創(chuàng)新引智計(jì)劃(B14006)

2016-05-09 修改稿收到日期：2016-06-21

楊仁樹 男，教授，博士生導(dǎo)師，1963年生

付曉強(qiáng) 男，博士生，1984年生 E-mail：fuxiaoqiang1984@163.com

TD235.1

A

10.13465/j.cnki.jvs.2017.03.007