k值優(yōu)化VMD-小波包分析聯(lián)合降噪方法在隧道爆破信號(hào)中的應(yīng)用?

王海龍 李云赫 趙 巖

①河北建筑工程學(xué)院 土木工程學(xué)院(河北張家口，075000)

②河北省土木工程診斷、改造與抗災(zāi)重點(diǎn)試驗(yàn)室(河北張家口，075000)

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

引言

爆破施工是隧道掘進(jìn)最常用的破巖方式，它所帶來的結(jié)構(gòu)振動(dòng)極有可能對(duì)在建隧道及周圍建筑物產(chǎn)生損害。爆破實(shí)測(cè)現(xiàn)場(chǎng)環(huán)境復(fù)雜，爆破振動(dòng)信號(hào)不可避免地會(huì)受到外界環(huán)境的影響。因此，對(duì)爆破振動(dòng)實(shí)測(cè)信號(hào)進(jìn)行降噪處理對(duì)后續(xù)的信號(hào)分析具有重要的實(shí)際意義[1]。

傳統(tǒng)傅里葉變換用于處理平穩(wěn)信號(hào)具有良好的分析效果，但不適用于非平穩(wěn)、非線性的爆破振動(dòng)信號(hào)[2]；短時(shí)傅里葉變換在傳統(tǒng)傅里葉變換的理論基礎(chǔ)上引入窗函數(shù)獲取時(shí)域信息，但仍無法滿足非穩(wěn)態(tài)信號(hào)變化的頻率需求[3]；小波變換方法解決了傅里葉變換窗口大小不能隨頻率變化的問題，但在去噪過程中要進(jìn)行閾值和基函數(shù)的選擇，缺乏自適應(yīng)性[4]；經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)與局域均值分解(LMD)克服了用戶設(shè)置基函數(shù)盲目性問題，但無法避免模態(tài)混疊等問題[5-6]；總體平均經(jīng)驗(yàn)?zāi)B(tài)分解(EEMD)在分解初始加入一定量的白噪聲，一定程度上減弱了模態(tài)混疊現(xiàn)象，但其不穩(wěn)定性和計(jì)算程序復(fù)雜的缺點(diǎn)仍未得到解決[7]；互補(bǔ)集合經(jīng)驗(yàn)?zāi)B(tài)分解(CEEMD)使用正白噪聲和負(fù)白噪聲來增加額外的噪聲，但是其在分解時(shí)會(huì)出現(xiàn)較多的偽分量[8]。

本文中，針對(duì)隧道爆破振動(dòng)信號(hào)采集過程中存在的強(qiáng)噪聲干擾問題和利用變分模態(tài)分解(variational mode decomposition，簡(jiǎn)稱VMD)方法進(jìn)行信號(hào)分析時(shí)關(guān)鍵參數(shù)非最優(yōu)問題，提出一種基于優(yōu)化k值的VMD聯(lián)合小波包分析的降噪方法[9-10]。該方法兼具VMD和小波包分析的優(yōu)點(diǎn)。從能量守恒的角度確定最優(yōu)分解層數(shù)，在非遞歸的理論框架下，通過構(gòu)造并求解約束變分問題實(shí)現(xiàn)信號(hào)分解?？朔四B(tài)混疊和計(jì)算量大等問題，表現(xiàn)出良好的噪聲魯棒性。

1 VMD參數(shù)優(yōu)化原理

1.1 VMD原理

VMD算法一次性將復(fù)雜信號(hào)分解為k個(gè)調(diào)幅-調(diào)頻信號(hào)[11]，表達(dá)式為

式中:Ak(t)為uk(t)的瞬時(shí)幅值；φk(t)作為中間參數(shù)，對(duì)其求導(dǎo)可得各個(gè)調(diào)幅-調(diào)頻信號(hào)的中心頻率ωk(t)，即ωk(t)＝φk′(t)。

通過控制帶寬來避免模態(tài)混疊現(xiàn)象，構(gòu)造約束變分問題，表達(dá)式[12]如下:

為實(shí)現(xiàn)將約束變分問題轉(zhuǎn)變?yōu)榉羌s束變分問題，引入懲罰因子α和Lagrange乘法算子λ(t)，構(gòu)造增廣Lagrange表達(dá)式:

經(jīng)過傅里葉等距變換等過程可得到

＾uk(ω)經(jīng)過兩次傅里葉逆變換，由頻域轉(zhuǎn)換到時(shí)域，進(jìn)而求解實(shí)部uk(t)，實(shí)現(xiàn)將信號(hào)分解為k個(gè)模態(tài)分量。

1.2 模態(tài)數(shù)k優(yōu)化準(zhǔn)則

VMD過程遵循能量守恒原則，故可從各模態(tài)分量與原始信號(hào)之間的能量關(guān)系入手，對(duì)VMD算法參數(shù)進(jìn)行優(yōu)化。當(dāng)VMD分解程度過大，即k值過大時(shí)，由于存在過分解的分量，分量能量線性之和大于原始信號(hào)的能量，并不合理。因此，從能量守恒的角度優(yōu)化VMD模態(tài)數(shù)k顯得尤為重要。

原信號(hào)或分量信號(hào)的能量計(jì)算公式為

式中:E為信號(hào)的能量；x(i)為信號(hào)序列；n為采樣點(diǎn)數(shù)。

同一信號(hào)在VMD算法中設(shè)定不同模態(tài)數(shù)k分解后，模態(tài)分量能量線性總和大小不同，為評(píng)價(jià)信號(hào)在不同模態(tài)數(shù)k值分解后分解程度水平，引入分解能量差值參數(shù)λ，公式為:

式中:Ex為原信號(hào)的總能量；Ej為原信號(hào)經(jīng)VMD處理后對(duì)應(yīng)的第j個(gè)模態(tài)分量的能量。

由式(7)可知:λ的大小與分解程度呈正相關(guān)，λ越大，過分解現(xiàn)象越嚴(yán)重；而λ越接近或等于零，則說明分解程度趨于合適。

繪制k-λ曲線?；诘饶芰糠纸庠?，VMD模態(tài)數(shù)k優(yōu)化步驟如下:在k逐次增加的條件下依次對(duì)信號(hào)進(jìn)行VMD預(yù)處理，分別計(jì)算分解能量差值參數(shù)λ；當(dāng)λ出現(xiàn)明顯突變時(shí)停止計(jì)算，取突變點(diǎn)處對(duì)應(yīng)的k為最優(yōu)模態(tài)數(shù)。

2 k值優(yōu)化VMD-小波包分析聯(lián)合降噪原理

2.1 小波包降噪原理

實(shí)現(xiàn)小波降噪過程時(shí)，首先自定義小波基函數(shù)對(duì)含噪信號(hào)進(jìn)行分解，再對(duì)高頻分量進(jìn)行閾值濾波，最后重構(gòu)信號(hào)得到降噪信號(hào)[13]。小波包分析以小波變換為基礎(chǔ)，對(duì)小波變換中未處理的高頻分量再次細(xì)化分解。相比于小波降噪，小波包分析具有更高的時(shí)頻分辨率，可以進(jìn)一步消除高頻部分存在的噪聲余量，提高去噪精度[14]。

2.2 k值優(yōu)化VMD-小波包分析聯(lián)合降噪

圖1 為基于k值優(yōu)化的VMD-小波包分析聯(lián)合降噪流程圖。首先，將原始爆破振動(dòng)信號(hào)利用基于k值優(yōu)化的VMD法進(jìn)行分解，得到k個(gè)模態(tài)分量；通過比較相關(guān)系數(shù)、方差貢獻(xiàn)率兩個(gè)指標(biāo)，篩選出包含噪聲的模態(tài)分量，并對(duì)篩選出的含噪分量利用小波包分析做閾值降噪處理；最后，將降噪處理后的含噪分量與優(yōu)勢(shì)分量重構(gòu)，得到純凈信號(hào)。

圖1 降噪流程圖Fig.1 Flow diagram of denoising

3 仿真分析

構(gòu)造多頻信號(hào)疊加仿真信號(hào)，表達(dá)式為

式中:γ(t)為人工噪聲。

VMD算法中分解完備性平衡參數(shù)α取默認(rèn)值2000[15]，利用Matlab軟件做仿真信號(hào)分析，信號(hào)時(shí)域波形如圖2所示。在不同k值條件下，依次對(duì)仿真信號(hào)進(jìn)行VMD預(yù)處理，并分別計(jì)算原信號(hào)能量和各模態(tài)分量能量之和，結(jié)果如表1所示。

圖2 多頻疊加信號(hào)時(shí)域圖Fig.2 Time-domain diagram of multi-frequency superposition signals

表1 仿真信號(hào)VMD各能量參數(shù)Tab.1 Energy parameters of VMD of simulation signals

原信號(hào)中由于存在人工噪聲的干擾，所以在不同的k值條件下總能量Ex不嚴(yán)格相等，存在波動(dòng)現(xiàn)象。直觀分析表1可得:當(dāng)k≤4時(shí)，λ恒等于0，即信號(hào)經(jīng)VMD處理后各模態(tài)分量能量之和嚴(yán)格等于原信號(hào)能量；從k＝5開始，由于存在虛構(gòu)的分量，λ逐漸增大。

根據(jù)上述模態(tài)數(shù)和選擇準(zhǔn)則，將突變點(diǎn)k＝4取為最佳分解層數(shù)。探討信號(hào)在此參數(shù)下的分解水平，進(jìn)一步分析VMD信號(hào)的頻譜特征，如圖3、圖4所示。結(jié)合圖4中的頻譜圖，仿真信號(hào)由4組多頻信號(hào)和人工噪聲組成，在k＝4的條件下，經(jīng)過VMD處理的信號(hào)u1、u2、u3、u4分別為主頻為45、81、95、110 Hz的單信號(hào)，各分量信號(hào)的特征頻率與原始仿真信號(hào)有1 Hz的差別，均在允許誤差范圍內(nèi)。由此可以驗(yàn)證，基于k值優(yōu)化的VMD法分解效果較好。

圖3 仿真信號(hào)VMD各分量的時(shí)域Fig.3 Time-domain of VMD of simulation signals

圖4 仿真信號(hào)VMD各分量的頻域Fig.4 Frequency-domain of VMD of simulation signals

針對(duì)于上述模擬信號(hào)，利用基于k值優(yōu)化的VMD將其分解為一系列本征模態(tài)分量，利用Matlab中互相關(guān)函數(shù)的概念計(jì)算各個(gè)模態(tài)分量與原始信號(hào)的相關(guān)性，并依據(jù)相關(guān)性絕對(duì)值的大小對(duì)篩選出需要進(jìn)行處理的含噪分量進(jìn)行降噪處理，將處理后的結(jié)果與未經(jīng)處理的模態(tài)分量重構(gòu)得到最終的純凈信號(hào)[16]。降噪處理后的純凈信號(hào)與原始仿真信號(hào)及初始信號(hào)的對(duì)比關(guān)系如圖5所示。

由圖5(a)和圖5(c)可看出，處理后的純凈信號(hào)在保留仿真信號(hào)特征信息的同時(shí)，基本剔除了隱藏于其中的噪聲信息，且降噪后的信號(hào)較原始波形曲線更為光滑，降噪效果良好。從圖5(b)和圖5(c)可以看出，去噪后的純凈信號(hào)與初始信號(hào)在信號(hào)峰值及局部波形特征上的相似吻合度較好。

圖5 經(jīng)降噪處理后的純凈信號(hào)與原始仿真信號(hào)及初始信號(hào)的對(duì)比Fig.5 Comparison of the pure signal after denoising with the original simulation signal and the initial signal

為了進(jìn)一步研究兩種信號(hào)的相似程度，利用Matlab中互相關(guān)函數(shù)corrcoef進(jìn)行計(jì)算，得到純凈信號(hào)與仿真信號(hào)的相關(guān)系數(shù)為0.962 3，去噪后的信號(hào)與初始信號(hào)表現(xiàn)出較好的相關(guān)性。綜上所述，基于k值優(yōu)化的VMD聯(lián)合小波包閾值去噪方法在仿真信號(hào)領(lǐng)域應(yīng)用效果良好，為后續(xù)對(duì)實(shí)測(cè)爆破振動(dòng)信號(hào)的降噪處理提供了一定的理論基礎(chǔ)。

4工程實(shí)例分析

4.1 工程背景

試驗(yàn)以新建河北省張家口市太子城至內(nèi)蒙古自治區(qū)錫林浩特段崇禮隧道爆破施工工程為背景。全斷面爆破施工，使用1～13段非電毫秒雷管和2＃巖石乳化炸藥，最大單段藥量43.2 kg，循環(huán)爆破總藥量204 kg。采用TC-4850N爆破振動(dòng)監(jiān)測(cè)儀采集到的某次爆破施工中爆破振速原始信號(hào)，如圖6所示。從圖6可知，受爆破施工現(xiàn)場(chǎng)復(fù)雜環(huán)境的影響，爆破振速時(shí)程曲線中夾雜著各種毛刺噪聲信號(hào)。

圖6 原始爆破振動(dòng)信號(hào)時(shí)域波形圖Fig.6 Time-domain waveform of original blasting vibration signals

4.2 基于k值優(yōu)化VMD-小波包閾值聯(lián)合降噪實(shí)現(xiàn)

采用上述分析理論，對(duì)采集到的隧道爆破振動(dòng)信號(hào)在不同的k值條件下進(jìn)行VMD預(yù)處理，并分別計(jì)算原信號(hào)能量Ex、各分量的能量Ej總和∑Ej、分解能量差值參數(shù)λ，直到λ出現(xiàn)較明顯的變化，停止計(jì)算。

計(jì)算結(jié)果如表2所示。

表2 爆破振動(dòng)信號(hào)VMD處理后的各能量參數(shù)Tab.2 Energy parameters of blasting vibration signals after VMD

由表2可看出:當(dāng)k＝2時(shí)，λ嚴(yán)格等于0；k＝3～6時(shí)，λ出現(xiàn)0.003 3～0.011 8之間的微小波動(dòng)，但波動(dòng)范圍較小，可認(rèn)為仍處于平穩(wěn)無突變狀態(tài)；而在k＝7時(shí)，λ陡增至0.095 6，隨后不斷地增加。故判定k＝6是λ的變化轉(zhuǎn)折點(diǎn)。按照模態(tài)數(shù)k選取準(zhǔn)則，將k＝6取為VMD最優(yōu)分解層數(shù)，并由此對(duì)爆破振動(dòng)信號(hào)進(jìn)行VMD處理。各個(gè)模態(tài)分量的時(shí)程曲線見圖7。

圖7 爆破振動(dòng)信號(hào)各模態(tài)分量時(shí)程曲線Fig.7 Time history curve of each modal component of blasting vibration signals

根據(jù)式(9)計(jì)算每一個(gè)模態(tài)分量與原始爆破信號(hào)的相關(guān)系數(shù)ri，計(jì)算結(jié)果見表3。

3.2.1 基本信息管理模塊。完成校企雙方指導(dǎo)老師、頂崗實(shí)習(xí)學(xué)生個(gè)人基本信息的維護(hù)和查詢功能。校內(nèi)指導(dǎo)老師的基本信息應(yīng)包括：系部、姓名、郵箱、QQ、聯(lián)系電話等。企業(yè)指導(dǎo)老師的基本信息應(yīng)包括：姓名、聯(lián)系電話、郵箱、QQ、企業(yè)名稱、企業(yè)地址、工作崗位等。學(xué)生的基本信息應(yīng)包括：學(xué)號(hào)、姓名、系別、班級(jí)、聯(lián)系電話、郵箱、QQ、實(shí)習(xí)單位名稱、實(shí)習(xí)單位地址及實(shí)習(xí)崗位等。

表3 模態(tài)分量的相關(guān)系數(shù)Tab.3 Correlation coefficient of intrinsic mode function

式中:xi表示各模態(tài)分量；y表示原始爆破信號(hào)；ˉxi、ˉy分別為xi、y的平均值。

分析表3可得，經(jīng)VMD處理后的各模態(tài)分量與原始信號(hào)的相關(guān)系數(shù)差異較大。其中，u4、u5、u6相關(guān)系數(shù)均小于0.15，初步判定為含有較多噪聲的含噪分量；其余分量可以看作是優(yōu)勢(shì)模態(tài)分量。

定義模態(tài)分量平方的算術(shù)平均值與其均值的平方之差為模態(tài)分量的方差s，通過各個(gè)模態(tài)分量的方差貢獻(xiàn)率ε校核上述選擇的可行性[12，17]。模態(tài)分量的方差及方差貢獻(xiàn)率計(jì)算如下:

計(jì)算得到各模態(tài)分量的方差貢獻(xiàn)率，見表4。

表4 模態(tài)分量的方差貢獻(xiàn)率Tab.4 Variance contribution rate of intrinsic mode function

由表4可知，u1、u2、u3模態(tài)分量的方差貢獻(xiàn)率較大。將相關(guān)性較差、方差貢獻(xiàn)率較小的u4、u5、u6重構(gòu)得到新的分量u′，以u(píng)′作為新信號(hào)進(jìn)行小波包分析。工程爆破中的振動(dòng)頻率一般在100 Hz左右，既要保證采集到的信號(hào)完整，又要避免引入高頻信號(hào)，故需要將采樣頻率設(shè)為信號(hào)頻率的10～100倍。結(jié)合前期爆破振動(dòng)監(jiān)測(cè)試驗(yàn)，本次監(jiān)測(cè)將爆破振動(dòng)監(jiān)測(cè)儀TC-4850N的信號(hào)采樣頻率設(shè)置為5 000 Hz。

綜合考慮支撐長度、消失矩、對(duì)稱性、正則性以及相似性等因素，結(jié)合采樣定理[16]，選擇具有良好緊支撐系性、光滑性及近似對(duì)稱性的db8小波基函數(shù)。選取分解層數(shù)為3、5、7和9進(jìn)行仿真計(jì)算。結(jié)果顯示，在分解層數(shù)為7時(shí)，去噪效果最佳，然后利用默認(rèn)的軟閾值函數(shù)進(jìn)行降噪處理[14]。

最后，將降噪處理過的含噪分量u′與優(yōu)勢(shì)分量u1、u2、u3重構(gòu)，得到最終的純凈信號(hào)，如圖8所示。由圖8可知，與原始信號(hào)相比，經(jīng)VMD-小波包閾值聯(lián)合降噪后的純凈信號(hào)在保證局部波形特征及峰值不變的基礎(chǔ)上，基本消除了原始信號(hào)中存在的噪聲分量。

圖8 基于k值優(yōu)化VMD聯(lián)合小波包閾值降噪前、后信號(hào)的對(duì)比Fig.8 Comparison of signals before and after threshold denoising by k-value optimization VMD combined with wavelet packet

4.3 降噪效果評(píng)價(jià)

通常利用信噪比η和均方根差μ兩種指標(biāo)衡量爆破振動(dòng)信號(hào)降噪效果。信噪比越大，均方根差越小，則降噪效果越好[18]。信噪比、均方根差的計(jì)算公式為

式中:zi(t)為原始信號(hào)為經(jīng)降噪處理后的信號(hào)；n為信號(hào)長度。

評(píng)價(jià)VMD聯(lián)合小波包降噪效果的同時(shí)，再次驗(yàn)證VMD法中模態(tài)數(shù)選擇是否為最佳。分別選取不同k值的VMD聯(lián)合小波包對(duì)信號(hào)進(jìn)行上述降噪操作，并計(jì)算信噪比、均方根誤差，計(jì)算結(jié)果見表5。從表5中可知，k＝6時(shí)，降噪效果最好。

表5 不同k值的VMD聯(lián)合小波包閾值降噪效果對(duì)比Tab.5 Comparison of threshold denoising effects of VMD with different k-values combined with wavelet packet

為進(jìn)一步驗(yàn)證此方法的降噪效果，對(duì)實(shí)測(cè)爆破振動(dòng)信號(hào)分別采用小波包閾值降噪法、EMD-小波包聯(lián)合降噪法、CEEMD-小波包聯(lián)合降噪法進(jìn)行分析，與VMD-小波包聯(lián)合降噪法的對(duì)比結(jié)果見表6。

表6 降噪效果對(duì)比Tab.6 Denoising effect comparison

對(duì)比分析可得，經(jīng)k值優(yōu)化的VMD-小波包聯(lián)合降噪法的信噪比(95.230 1)最大，均方根差最小(2.503 12×10-5)，降噪效果較好。

5 結(jié)論

1)VMD方法基于等能量理論模態(tài)數(shù)選取準(zhǔn)則確定參數(shù)k，能有效避免信號(hào)分解不足或過分解。

2)實(shí)際工程中采集到的原始含噪信號(hào)經(jīng)過k值優(yōu)化VMD-小波包閾值法聯(lián)合降噪后，可以在消除噪聲的同時(shí)，有效地保留原始信號(hào)中的細(xì)節(jié)特征，為后續(xù)準(zhǔn)確分析爆破信號(hào)奠定了基礎(chǔ)。

3)與小波包閾值降噪法、EMD-小波包聯(lián)合降噪法、CEEMD-小波包聯(lián)合降噪法相比，k值優(yōu)化的VMD-小波包聯(lián)合降噪方法信噪比大，均方根差小，降噪效果較好。