隧道爆破振動信號時頻譜交叉項干擾抑制方法

付曉強， 俞 縉， 戴良玉， 秦雙雙

(1.三明學院 建筑工程學院，福建 三明 365004； 2.華僑大學 福建省隧道與城市地下空間工程技術研究中心，福建 廈門 361021； 3.三明科飛產(chǎn)氣新材料股份有限公司，福建 三明 365500)

鉆爆法是現(xiàn)階段土石方開挖中應用最為廣泛的施工方法之一[1-2]。爆破振動監(jiān)測作為爆破損傷效應最行之有效的評價手段，其時頻分析對于爆破方案優(yōu)化和參數(shù)調(diào)整具有至關重要的作用。受測試環(huán)境、儀器性能和地質(zhì)條件等因素的影響，爆破信號中普遍存在不相干分量，致使其時頻譜中均不同程度存在交叉項干擾，從而對信號分析精度產(chǎn)生不利影響。而傳統(tǒng)的以傅里葉變換為基礎的信號分析存在表示結果不稀疏，不具備自適應性且時頻分辨率差的缺陷，未能有效解決交叉項對信號辨識度的影響[3]。

近年來，在改變傳統(tǒng)爆破信號表示方面取得了很大的進展，一定程度上豐富和推動了爆破信號時頻分析理論的發(fā)展。如關曉磊等[4]采用極值點對稱延拓，消除了希爾伯特-黃變換方法中的端點效應缺陷，改善了算法中瞬時頻率出現(xiàn)負值問題，提高了爆破信號時頻分析精度。郭濤等[5]將頻率切片小波變換應用于爆破信號分析過程中，通過細化信號頻率區(qū)間，提取到更為真實的信號時頻特征，一定程度上抑制了交叉項干擾。楊仁樹等[6]對隧道爆破信號進行集合經(jīng)驗模態(tài)分解，利用各模態(tài)信號分形特征確定優(yōu)勢分量并重構，采用平滑偽魏格納分布(smooth pseudo Wigner-Ville distribution,SPWVD)對重構信號進行了時頻分析，削弱了低頻交叉項干擾的影響。同時應注意到，由于受爆破復雜環(huán)境和儀器參數(shù)設置等的影響，信號中含有的趨勢項和噪聲干擾成分會進一步污染其時頻平面上的能量分布。趨勢項和噪聲通常為時頻譜產(chǎn)生交叉項干擾的主要來源，趨勢項的剔除對信號時頻細節(jié)的表征具有明確的物理意義。如龍源等[7]分別采用最小二乘法、小波及經(jīng)驗模態(tài)分解(empirical mode decomposition,EMD)3種方法對爆破信號進行趨勢項消除，并對3種方法的分析效果進行了對比，初步探討了EMD算法在信號趨勢項處理方面的優(yōu)良特性，但EMD算法對信號端點處理不當卻會引起模態(tài)混疊現(xiàn)象。上述諸多算法對爆破信號分析均取了良好的效果，但參數(shù)設置及處理過程較為復雜，面對海量的監(jiān)測數(shù)據(jù)，算法缺乏自適應性嚴重制約了腦力時間的合理分配，甚至無法窮盡時間去精細刻畫信號的時頻特征。

針對上述難題，相關學者將匹配追蹤(matching pursuit，MP)算法引入到爆破信號分析領域[8-9]。該算法實現(xiàn)了信號靈活表達，其基函數(shù)采用被稱之為原子字典的過完備冗余函數(shù)(過完備原子庫或原子字典)集合取代，可將信號表示為原子字典中少數(shù)向量的線性組合。這樣，信號在過完備字典上的分解結果將是稀疏的，信號信息或能量集中在少數(shù)原子上，從而便于對信息的進一步處理。本文采用匹配追蹤-平滑偽魏格納分布(matching pursuit-smooth pseudo Wigner Ville distribution, MP-SPWVD)組合方法，對實測爆破振動信號時頻譜中交叉項干擾消除難題進行了分析與實踐，驗證了該算法用于多因素導致的爆破信號時頻交叉項消除的效果。

1 基本算法

1.1 MP算法

MP算法是在信號處理和函數(shù)逼近領域獨立發(fā)展起來的，是信號稀疏表示最重要的算法之一[10]。它將信號按照字典原子逐步分解，通過多次迭代并選取原子庫中匹配度高的原子近似逼近，最終獲得最稀疏且具有明確物理意義的信號表示，其分解過程如圖1所示：① 選擇與信號Xn相似度高的原子Ψn并計算其投影值ɑn和差值信號Xn+1；② 對殘差信號Xn+1重新選擇原子進行匹配，獲取其在最相近原子Ψn+1上的投影，得到差值信號Xn+2；③ 重復上述過程，當殘差信號能量小于設定的閾值則計算終止。

圖1 信號投影過程示意圖

設D為有限維Hilbert空間H中的一個超完備字典，假定f為待分解信號(f∈H)，數(shù)據(jù)長度為N，則D滿足[11]

(1)

式中:gγ為子波算子;Γ為伽馬函數(shù)。

(2)

(3)

(4)

(5)

1.2 MP-SPWVD算法實現(xiàn)過程

Cohen類非線性時頻分布中的魏格納分布(Wigner-Ville distribution，WVD)算法會產(chǎn)生交叉項，偽魏格納分布(pseudo Wigner-Ville distribution，PWVD)算法可一定程度上削弱交叉項，而相對WVD和PWVD分布，SPWVD時頻解析精度更高，具體算法詳見文獻[14]。對于任意給定爆破振動信號X(t)，利用Hilbert變換可預先提取其瞬時頻率和相位信息，從而很大程度上減少信號原子庫的匹配數(shù)量，提高運算效率。具體步驟如下[15-17]。

步驟1選擇Gabor原子并對其擴展，形成匹配過完備原子數(shù)據(jù)庫Di(i=1,2,…,I)。

步驟2令初始差值信號r0等于待分解的初始信號X(t)，并將其代入過完備原子庫中，篩選并預先確定子波分解原子庫范圍。

步驟3經(jīng)過m次迭代后，從確定的原子庫內(nèi)找出與差值信號rm局部特征最為匹配的原子，也就是內(nèi)積最大的原子dmi，將匹配到的原子從差值信號中減去，獲取新的差值信號rm+1。

步驟4重復步驟2和步驟3，直至差值信號滿足規(guī)定的條件。最終，X(t)被分解為

(6)

步驟5剔除信號中干擾成分R，分別對通過反復迭代獲取的各分量信號cmidmi進行SPWVD運算并疊加，從而獲得在時頻平面的最優(yōu)化特征分布。信號MP-SPWVD算法流程，如圖2所示。

圖2 MP-SPWVD信號處理流程

2 仿真信號分析

為了說明算法的優(yōu)良特性，采用以子波瞬時頻率變化表示信號信息的頻率調(diào)制方式建立含4個高斯分量復雜調(diào)制信號，如圖3(a)所示。各分量歸一化頻率中心分別為0.15,0.25,0.35和0.45，峰值均為2 cm/s，調(diào)制信號由各分量信號實部采用線性疊加形成，盡可能模擬多段別雷管起爆多頻率多振型特點,如圖3(b)所示。其中x1,x4波形起始于32 ms，x2,x3初振時刻為96 ms，為了提高運算效率，設定分量波形振動時長均為128 ms。調(diào)制信號WVD時頻分布，如圖3(c)所示。

(a) 調(diào)制信號

可以看出：信號時頻面上除了4個真實高斯分量的能量分布外，亦包含任意兩個分量時頻中心連線處出現(xiàn)的交叉項，信號4個有效分量共產(chǎn)生6個交叉項，分別為交叉項1(由分量x1,x4產(chǎn)生)、交叉項2(由分量x3,x4產(chǎn)生)、交叉項3(由分量x2,x3產(chǎn)生)、交叉項4(由分量x1,x2產(chǎn)生)、交叉項5(由分量x1,x3產(chǎn)生)和交叉項6(由分量x2,x4產(chǎn)生)，且對角交叉項5、交叉項6互相疊加產(chǎn)生了重疊。信號PWVD時頻分布將交叉項個數(shù)減少至兩個，分別為x1,x4(交叉項1)和x2,x3(交叉項2)產(chǎn)生,如圖4(a)所示。PWVD一定程度上減少了交叉項的個數(shù)，但仍未能從根本上抑制交叉項的產(chǎn)生。對調(diào)制信號進行MP迭代分解后，可準確獲得4個高斯分量x1～x4的波形曲線，如圖4(b)所示。對各分量分別求取SPWVD后疊加，便得到信號真實的時-頻-能分布，如圖4(c)所示。

(a) PWVD時頻譜

調(diào)制信號WVD和PWVD在時頻域本不應有能量存在的局部出現(xiàn)交叉項。而MP-SPWVD算法將信號分解為相互獨立且包含局部結構特征的有限個原子的線性組合，由于原子間相互獨立，則其SPWVD分布亦是相互獨立且不存在交叉項，通過相互獨立原子SPWVD分布疊加得到信號時頻表達。因此，MP-SPWVD可有效避免交叉項的影響并改善多分量信號時頻分辨率，進而獲得更為真實的時頻譜分布。

3 實例分析

3.1 隧道爆破信號

圖5(a)為典型的隧道多段別爆破信號，該信號波峰值為3.36 cm/s,波谷值為2.77 cm/s,主頻為180.22 Hz。隧道掘進斷面41.2 m2，爆破采用MS1～MS15段別雷管跳段使用，其中，掏槽孔采用MS1、MS3段，單孔裝藥量0.8 kg；輔助孔采用MS5、MS7、MS9、MS11及MS13段，單孔裝藥量0.6 kg；周邊孔采用MS15段，單孔裝藥量0.4 kg；底孔為MS15，單孔裝藥量0.6 kg。采用為采用“短進尺+弱爆破”方案，掏槽眼深度1.5 m，其余炮眼1.2 m，單循環(huán)總裝藥量為56.2 kg。

圖5(b)功率譜表明：信號在低頻段存在明顯的趨勢項，在起始頻率0～3 Hz附近產(chǎn)生0.02 dB/Hz的奇異值，而在大于500 Hz的頻帶出現(xiàn)顯著的低幅高斯白噪聲，信號頻譜中存在明顯的偽信息，其在信號時頻譜上產(chǎn)生無法解釋的交叉項干擾，嚴重影響了對信號特征的解讀和判別。圖5(c)中MP重構信號功率譜有效克服了低頻趨勢項對信號真實頻譜的影響，獲得了具有明確物理意義的信號頻譜特征。

(a) 隧道爆破信號

3.2 信號MP分解與重構

在信號頻譜分析中，Gaussian函數(shù)具有良好的時頻聚集性[18]。因此，文中建立的Gabor原子庫中的每個原子都是由Gaussian函數(shù)經(jīng)伸縮、平移和調(diào)制得到的，如圖6(a)所示。

通過引入本征時間尺度分解得到信號一系列不同頻率段的PR(proper rotation)分量[19-20]。對各PR分量進行Hilbert變換，獲取信號的瞬時優(yōu)勢頻率和相位，隨后將其代入MP算法中，可大大減少程序循環(huán)步數(shù)。采用1.1節(jié)的MP算法對原始信號進行處理，建立與分析信號相似度高，長度為10 ms的過完備稀疏Gabor原子字典庫。實踐證明，迭代次數(shù)滿足10n(n≤3且為正整數(shù))指數(shù)關系時能取得較為理想的分析效果，n較小時(n=1)存在分解不徹底缺陷，而n較大時(n=3)易產(chǎn)生過分解，同時運算時長激增。根據(jù)爆破信號特點及子波分量頻帶限制，此處設置迭代次數(shù)為100(n=2)，從而將原信號分解為7個子波分量，如圖6(b)所示。由于MP算法是以信號時頻能量高低對其進行分解的，因此將迭代獲得的前4個能量占比高的分量進行線性疊加，得到重構后的真實信號如圖6(c)所示。

從圖6中可知，重構信號很好地繼承了原始信號的波動特征，同時對信號波形局部進行了平滑擬合。為了便于分析，分別標識原始信號、MP處理信號為X1,X2，采用文獻[21]中的交叉小波變換方法，求取兩者的相關性凝聚譜,如圖7所示。

(a) Gobor原子庫

圖7 信號相關性譜

相關性凝聚譜X2和X1在時域和主振頻域上均為顯著正相關。X2與X1在低頻部分(<32 Hz)和高頻部分(>512 Hz)為負相關，有效濾除了信號中的低頻趨勢項和高頻噪聲且相關性在時間軸上連續(xù)統(tǒng)一，顯著性檢驗基本貫穿小波影響椎范圍內(nèi)的整個時頻域。定義互相關系數(shù)為

λs1,s2=max[Rs1(τ)]、max[Rs2(τ)]，

λs1,s2∈(0,1)

(7)

式中:Rs1(τ),Rs2(τ)為分別為X1,X2的自相關值;λs1,s2值越大，則兩信號間互相關性越強。計算X2與X1之間的互相關系數(shù)為0.968，交叉小波變換清晰描述X2與X1在時域和頻域不同層次的波動性和共振相位相關度。MP算法在信號處理過程中對低頻趨勢項進行了有效校正，同時對高頻噪聲分量起到了良好的抑制效果，體現(xiàn)了算法信號處理優(yōu)勢。

3.3 時頻分析

為了分析爆破信號交叉項干擾的抑制效果，分別采用WVD，PWVD和MP-SPWVD算法求取信號的時頻譜，如圖8所示。隧道爆破信號WVD時頻譜交叉項干擾嚴重，時頻譜交叉項出現(xiàn)在原本不存在能量的位置，時頻分布的可辨性差。同時，段間延期時間越短，交叉項越明顯。PWVD譜削弱并剔除了信號時頻面上的高頻交叉項，但對于低頻交叉項分布卻有一定程度放大，表現(xiàn)在時頻面上的分布密度增強且分布形態(tài)趨于一致。MP-SPWVD時頻譜能量分布清晰可辨，實現(xiàn)了多段別微差爆破起爆能量的有效分割。爆破產(chǎn)生的能量主要聚集在掏槽段起爆0～35 ms時程內(nèi)，能量峰值位于歸一化頻率0.05～0.07處(采樣頻率為3 000 Hz，約150～210 Hz內(nèi))，低頻交叉項位于0～35 ms及歸一化頻率0～0.02處(約0～60 Hz)，而高頻交叉項主要位于波形時程尾部邊緣，且集中在歸一化頻率0.15～0.25處(約450～750 Hz)。MP-SPWVD算法時頻分布能夠隨著被分析信號局部特征自適應地調(diào)節(jié)原子基函數(shù)的特征參數(shù)，如時寬和帶寬，避免了傳統(tǒng)方法中窗效應的影響，還可很好地匹配信號的局部特征，得到的時頻譜對于低頻漂零和高頻邊緣交叉項均有優(yōu)良的抑制效果。

(a) WVD時頻譜

4 可靠性驗證

為了驗證1.2節(jié)組合方法的有效性，選取煤礦凍結立井某次爆破井壁振動信號，如圖9(a)所示。信號采用文獻[22]中井壁傳感器預埋法采集，測點位于爆破掌子面距離為6 m的井壁處。由于測點位于爆破近區(qū)，該信號含有明顯的偏離基線的趨勢項，導致其在時間軸上0.1 s后出現(xiàn)了“甩尾”(基線漂零)現(xiàn)象。爆破采用雷管為MS1～MS5 5個段別，其中，周邊眼采用的MS5段雷管的起爆及誤差區(qū)間為(110±15) ms，與趨勢項在時間軸上的出現(xiàn)區(qū)間重疊。因此，趨勢項引起的時頻譜交叉項干擾對其頻譜特征的解讀產(chǎn)生誤判，特別是0.1 s后的時頻平面上的能量分布尤其值得關注。

對圖9(a)中井壁振動信號分別采用小波方法和MP算法進行信號恢復，其結果見圖9(b)和圖9(c)。圖9對比表明：小波方法信號重構過程中，閾值選取不當易使信號產(chǎn)生了微幅的高頻振蕩噪聲，同時信號幅值大幅降低，損失了信號部分有用信息。MP算法校正并剔除了信號中的趨勢項分量，信號特征辨識度高。

(a) 立井爆破信號

求取原始信號WVD,PWVD和MP-SPWVD時頻譜分別如圖10(a)～圖(c)所示。

(a) WVD時頻譜

從圖10可知：WVD時頻譜窗函數(shù)不可調(diào)，致使其時頻分辨率是固定的，存在嚴重的交叉項干擾。PWVD時頻譜雖然改善了交叉項干擾的影響程度，但對于時頻譜在頻率邊界處的交叉項干擾有一定程度的放大。MP-SPWVD時頻譜是通過相互獨立的原子重構子信號SPWVD相互疊加構成的，可有效抑制交叉項的產(chǎn)生，其更加符合微差爆破信號的局部特征，能夠更清晰揭示立井信號的時頻分布特性。

同時應注意到，MP-SPWVD算法也有一定不足之處，即當待提取振動信號數(shù)據(jù)較長時，原子數(shù)量將急劇增加，這種情形下會嚴重降低運算速度。選取文獻[23]中高斯調(diào)制信號、三正弦和典型爆破信號在時間軸上分別進行人為延拓至1 s,3 s,5 s，來比較不同信號長度處理所耗機時，具體如表1所示。

由表1可知，隨著信號長度和復雜性的增大，分析所耗機時也會增加，尤其對于采樣點密集或持續(xù)時長較大(大于5 s)的爆破信號計算量過大，嚴重時甚至導致內(nèi)存溢出。因此求解過程需耗費很長機時，對求解器配置和性能有特殊要求。

表1 不同信號CPU運行機時

5 結 論

(1) 爆破信號時頻譜提供了信號在時域和頻域的聯(lián)合分布，清晰描述了能量在時頻域上的變化關系。受測試環(huán)境及儀器自身原因的影響，爆破信號時頻譜中普遍存在交叉項干擾，交叉項的存在嚴重影響了對爆破信號信息特征的判別，致使時頻域能量分布在局部物理意義不明確，影響后續(xù)爆破參數(shù)調(diào)整優(yōu)化和效果評價。

(2) MP算法實現(xiàn)了爆破信號的稀疏化分解和降噪重構，MP-SPWVD分布可有效抑制了時頻譜交叉項，降低不確定因素對信號譜特征的干擾，精細刻畫爆破信號能量分布的非平穩(wěn)時變特征。交叉項抑制效果與雷管段別及延期時間密切相關。

(3) MP-SPWVD算法雖是一種高精度的信號恢復和時頻表征方法，但受信號時長的限制較為嚴重，尤其對于采樣密集或持續(xù)時長大于5 s以上的信號數(shù)據(jù)，求解速度受到一定的制約，所耗機時激增，對求解器配置和性能有特殊要求。后續(xù)研究過程中，應重點優(yōu)化長時信號計算速度，提高算法的高效性。