蔡改貧，李洋波，楊麗榮，黃祥海

（1.江西理工大學 機電工程學院，江西贛州 341000；2.江西省礦冶機電工程技術研究中心，江西贛州 341000）

聲發(fā)射特征參數(shù)蘊含著豐富的巖石損傷漸進破壞的前兆信息，通過研究聲發(fā)射特征，對于理解巖石破裂機理、預防因巖石破裂失穩(wěn)而造成的事故具有十分重要的意義[1]。聲發(fā)射技術可以提供巖體失穩(wěn)所釋放的大量信號源信息，因此可進一步判斷巖體內部結構的動態(tài)變化、細觀損傷的出現(xiàn)以及宏觀裂紋的發(fā)展，因此國內外普遍運用聲發(fā)射技術對巖體失穩(wěn)現(xiàn)象進行識別與預測[2]。

VMD 是一種可變尺度的非平穩(wěn)信號處理方法[3]，而ITD 對于大數(shù)據(jù)信號的實時處理有突出表現(xiàn)[4]，兩種算法因其特點受到眾多學者研究。Abdoos 等[5]將S 變換、VMD 與SVM 相結合，用于提取電能質量信號特征及模式識別，結果表明該方法識別分類效果良好；An 等[6]利用VMD 算法對軸系振動信號降噪處理，驗證該算法對其降噪的有效性；張雪英等[7]通過VMD 與PNN 相結合，有效的提取特征信息，將齒輪故障診斷準確率提高至96.875%；Frei 等[8]提出ITD 算法，該方法在計算效率、分析精度、自適應劃分等方面有突出表現(xiàn)，更適合對非平穩(wěn)信號的分析。此外，郭力和鄧喻[9]應用美國聲學物理公司PAC 的PCI-2聲發(fā)射實驗儀器測量磨削聲發(fā)射信號，并采用遺傳算法優(yōu)化BP 神經(jīng)網(wǎng)絡，使相對誤差可控制在6.22%以下。李啟月等[10]通過EEMD 與小波閾值的方法有效的消除噪聲、趨勢項的干擾，保留波形的真實信息，提高頻譜的分析精度。

從之前學者的降噪算法研究中可以發(fā)現(xiàn)，不論是小波變換與EMD/EEMD，還是VMD 和ITD，雖然都能夠實現(xiàn)對噪聲的抑制，但在對某些特殊信號時其缺陷就會展現(xiàn)出來，如VMD 算法的模態(tài)數(shù)確定問題、ITD 算法的端點效應等。由于巖石聲發(fā)射信號具有非平穩(wěn)、大樣本等特點，單一的降噪算法難以較好地抑制其噪聲，因此利用多種算法進行聯(lián)合降噪是值得考慮與有效的方法。

1 VMD、ITD 基本原理及改進VMD-ITD聯(lián)合降噪算法

1.1 VMD 基本原理

VMD 分解算法通過對變分問題的構造及求解，在頻域內實現(xiàn)信號的自適應分離，能夠對信號的局部特征進行細致地刻畫。VMD 算法將原始信號X分解為K個本征模態(tài)函數(shù)（Instrinsic mode function，IMF）uk(t)，假設每個uk(t)具有中心頻率為 ωk的有限帶寬，則可將其引入到變分問題中并對其進行求解，求解過程如下：

1）通過Hilbert 變換得到uk(t)的解析信號Xk(t)，并利用預估中心頻率對Xk(t)進行調制得到信號Xωk(t)，即

式中Xk(t)為解析信號。

式中 δ(t)為單位脈沖函數(shù)。

2）通過解調信號的高斯平滑度，即梯度的L2范數(shù)的平方對帶寬進行估計，其約束變分模型為：

式中：{uk}={u1,u2,···,uk}為IMF分量的集合；{ωk}={ω1,ω2,···,ωk}為分量對應的中心頻率，Hz；u(t) 為輸入信號；?t為對t的偏導。

3）通過引入二次懲罰項 α在復雜噪聲環(huán)境下提高信號的重構精度，并引入lagrange 算子 λ(t)保持約束條件嚴格性，進而將變分問題轉化為非約束問題，得到的增廣拉格朗日表達式為

式中：F(?)為傅里葉變換；n為迭代次數(shù)。

綜上所述，VMD 算法具體步驟如下：

1）將、、λ1以及n初始化；

2）根據(jù)式（5）、式（6）更新uk、ωk；

3）對λ 進行更新，更新式為

4）不斷迭代，直到滿足

式中 σ為判別精度。

當 σ滿足設定的終止條件時，迭代不再繼續(xù)，得到K個窄帶模態(tài)分量；否則返回繼續(xù)更新uk、ωk，直至滿足式（8）。

1.2 ITD 基本原理

固有時間尺度分解（ITD）是一種適宜于大數(shù)據(jù)的實時處理以及非平穩(wěn)信號分析的方法，它能夠精確地提取信號的瞬時特征[11-12]。它通過構建基線信號將含噪聲信號分解為多個固有旋轉分量（PRC）以及一個殘差分量，每個PRC 分量都包含信號的時頻信息。ITD 分解原理如下：

1）定義Λ 為基線提取算子，則對于含噪待分解信號Xt的表達式為

式中：Lt=ΛXt為基線信號；Ht=(1-Λ)Xt，Ht為PRC分量。

2）求取信號Ht的所有局部極值點并將其所對應的時刻記為 τk,k=1,2,···N（N為極值點個數(shù)），且令τ0=0。為便于觀察，分別用Xk和Lk表示X(τk)和L(τk)。若在 [0,τk]上Lt、Ht有定義，則信號Xt在時間[0,τk+2]范圍內有意義，在連續(xù)極值點的區(qū)間(τk,τk+1]上定義分段線性提取算子 Λ，基線信號表達式為：

式中：α為線性縮放因子，用于控制PRC 分量的幅值，α ∈(0,1)，在工業(yè)應用中，α的大小與噪聲水平有關，通常α 取0.5[13]。

ITD 分解全過程表達式為

式中：p為分解層數(shù)；為殘余分量。

1.3 互信息的加權分量信號重構方法

互信息法（MI）能夠反映某個分量包含原始信號的信息量，不同于相關系數(shù)，它不限于實值變量即具備一般性，且能較好的衡量非線性關系[14-15]。因此，本文采用互信息法進行信號重構，將各IMF 分量與原信號的互信息作為加權系數(shù)進行IMF 分量重構以達到降噪的目的?；バ畔⒌幕驹砣缦拢?/p>

1）對于兩個隨機變量X、Y，它們的互信息I(X,Y)利用概率與熵的概念進行定義，定義式的表達式為

式中：P(x,y)為X、Y的聯(lián)合概率分布函數(shù)；P(x)和P(y)分別為X、Y的邊緣概率分布函數(shù)。

2）由于P(x,y)在一般情況下難以得到，通常采用核密度估計進行替代，本文中采用正態(tài)多元核密度估計，因此式（13）變換為

式中：n為樣本數(shù)；h為核函數(shù)寬度；d為變量X的維數(shù)；T與d有關，當d=1時是方差，當d>1時是協(xié)方差矩陣，|T|為矩陣行列式。

互信息能夠描述兩個變量之間相關性，將其作為加權系數(shù)進行信號重構不會因舍棄某些分量造成信息丟失，并且能夠抑制相關性小的含噪聲分量，重構后的信號Xt表示為

式中：K為模態(tài)數(shù)；uk(t)為第K個分量；Ik(X、Y)為uk(t)與原信號的互信息。

1.4 改進VMD-ITD 的聲發(fā)射信號聯(lián)合降噪方法

VMD 分解過程中常用中心頻率法[16-17]來確定模態(tài)數(shù)K的值，但這一方法容易受到人為主觀因素的影響。針對上述問題，改進模態(tài)數(shù)K確定過程，并將其應用到改進VMD-ITD 的聲發(fā)射信號聯(lián)合降噪主要過程中，過程如下：

1）對原始信號X(t)進行VMD 分解，且將模態(tài)數(shù)K初始化為1，則得到一個IMF 分量u1(t)，計算u1(t)的能量E1的表達式為

其余參數(shù)都設為默認值。

2）將得到的前K個IMF 分量逐次累加，利用2-范數(shù)計算信號X(t)的能量E0，即

式中||X(t)為信號X(t)的2-范數(shù)。

分量能量比η 的計算式為

3）判斷 η的值，由于所設定的能量閾值 θ=0.95。因此，若η<95%，則更新模態(tài)數(shù)即K=K+1并重復上述步驟直至達到終止條件；若 η ≥95%，則停止VMD 分解并輸出K，得到K個分量uk(t)；

4）計算各分量與原信號的互信息Ik(X、Y)，并將其作為加權系數(shù)對各分量進行加權重構得到信號；

綜上所述，改進VMD-ITD 聯(lián)合降噪算法流程如圖1所示。

圖1 聯(lián)合降噪算法流程Fig.1 Joint denoising algorithm flowchart

2 仿真信號分析

2.1 仿真信號的構建

改進VMD-ITD 聯(lián)合降噪算法從能量角度確定其模態(tài)數(shù)，并根據(jù)各分量與原始信號的互信息量進行加權重構，最后對重構信號進行ITD 二次降噪。為說明該方法對非平穩(wěn)信號降噪的可行性與有效性，通過構造仿真信號進行降噪分析評價。實測信號中往往是多種類型信號進行疊加，因此本文結合Mitrakovic 等[18]與白瑞[19]的聲發(fā)射信號數(shù)學模型，構建的仿真信號y(t)由3個分信號及一個強度為0.1的隨機噪聲 ζ(t)構成，其表達式為

式中：Ai為第i個信號的幅值，V；αi為第i個信號的衰減系數(shù)；fi為第i個信號的頻率，Hz。

根據(jù)構造的聲發(fā)射數(shù)學模型，可以得到其加噪仿真的信號，并計算其頻譜，加噪仿真信號及其頻譜如圖2所示。

圖2 加噪仿真信號及其頻譜Fig.2 A simulated noisy signal and its spectrum

2.2 仿真信號的降噪結果分析

首先利用本文提出的基于分量能量比的改進VMD 對仿真信號進行分解，當達到終止條件（η ≥95%）時，得到一系列的IMF 分量如圖3所示。

圖3 仿真信號的IMFFig.3 IMF of the simulated signal

由圖3可知，通過分量能量比來確定模態(tài)數(shù)k=10，分解結果表明該方法能夠有效將各分信號分離；并且部分分量信號與原信號差異較大，如IMF1、IMF7、IMF8、IMF9和IMF10，這表明這些分量信號是與原信號相關性較低的噪聲分量。計算各分量的頻譜，結果如圖4所示。

圖4 各IMF 的頻譜Fig.4 Spectrum of each IMF

由圖4可知，分量IMF2～IMF7基本都在仿真信號的頻率上（200、400、600 Hz），其余分量頻率分布在更低的100 Hz 附近或更高的900 Hz 以上頻段內，說明IMF2～IMF7含有原始仿真信號的最主要信息，而其余分量包含的大多是噪聲信息。為抑制這些噪聲且盡可能保留完整信息，將所有分量與原信號的互信息量作為加權系數(shù)來重構信號，得到的各IMF 與加噪仿真信號的互信息如表1所示。

表1 各IMF 與加噪仿真信號的互信息量Tab.1 Mutual Information between each IMF and the simulated noisy signal

由表1可知，IMF1、IMF7、IMF8、IMF9和IMF10的互信息量明顯小于其他分量，說明這些分量包含原信號有效信息更少，這與實際波形也相符；在10個IMF 中互信息量最大的IMF5其值為0.591，說明IMF5包含較多的有效信息。根據(jù)式（17）對各分量進行加權得到重構信號，再通過ITD 分解并二次重構以達到進一步降噪的目的，得到的二次降噪信號及其頻譜如圖5所示。

圖5 二次降噪信號及其頻譜Fig.5 The second-order denoised signal and its spectrum

對比圖2和圖5可知，經(jīng)過ITD 二次降噪的信號其頻譜在高頻和低頻部分更加光滑，毛刺明顯減少，說明在改進VMD-ITD 能夠將信號的高頻和低頻噪聲濾除。為評價改進VMD-ITD 降噪性能，通過計算其降噪后的信噪比和均方根誤差并與VMD、ITD 算法進行對比，得到3種算法對同一信號降噪后的信噪比與均方根誤差對比結果如表2所示。

表2 3種算法降噪后信號的信噪比與均方根誤差Tab.2 Signal-to-noise ratio (SNR) and root mean square error (RMSE) of the denoised signals using three algorithms

由表2可知，改進VMD-ITD（（IP）VMD-ITD）算法的均方根誤差為0.158 6，在3種方法中最小，且其信噪比為14.287 7，明顯高于另外兩種方法，說明相比于單一的VMD、ITD 降噪方法本文提出的方法降噪性能更好，能夠在保留更多有效信息的前提下提高信噪比。

3 改進VMD-ITD 降噪算法在巖石聲發(fā)射信號處理中的應用

本次試驗由單軸壓縮試驗平臺和聲發(fā)射信號采集平臺兩部分組成，其中單軸壓縮試驗平臺采用無錫歐凱電子的萬能材料試驗機搭建的，該試驗機由控制臺和壓力機兩部分組成，控制系統(tǒng)由液壓動力源提供穩(wěn)定輸出，控制系統(tǒng)通過液壓源控制壓力機的壓頭進行試驗。單軸試驗平臺如圖6所示。

圖6 萬能材料試驗系統(tǒng)Fig.6 Universal material testing system

聲發(fā)射信號采集平臺由聲華科技的雙通道聲發(fā)射儀搭建而成，主要有聲發(fā)射傳感器、前置放大器、信號采集卡以及采集系統(tǒng)等幾個部分。其中聲發(fā)射傳感器是感知聲發(fā)射現(xiàn)象的必備器件，它的性能參數(shù)影響著信號的采集，本次試驗選用的傳感器參數(shù)如表3所示。

表3 傳感器SR150M 的基本參數(shù)Tab.3 Basic parameters of the sensors

聲發(fā)射傳感器安裝在試件中部軸線位置，為保證傳感器與試件緊密貼合，兩者的接觸面采用真空脂作為耦合劑并用透明膠布進行固定，傳感器布局如圖7所示。所準備的巖樣為鎢巖，如圖8所示。

圖7 傳感器布局Fig.7 Sensor layout

圖8 鎢巖巖樣Fig.8 Tungsten rock sample

為驗證本文提出方法在巖樣破裂聲發(fā)射信號噪聲抑制的有效性，從鎢巖破裂試驗中隨機選取一組聲發(fā)射信號進行降噪分析，選取的聲發(fā)射信號波形及其頻譜如圖9所示。

圖9 含噪原始信號及其頻譜Fig.9 Noise-containing original signal and its spectrum

運用改進VMD 對原始信號進行分解，根據(jù)分量能量比得到各IMF 分量如圖10所示。

圖10 含噪原始信號的IMFFig.10 IMFs of the noise-containing original signal

由圖10可知，分量能量比大于95%時原信號分解為8個IMF 分量，即分解模態(tài)數(shù)k=8；在所有分量中，IMF1與IMF2能較好地反映原信號波形，其余分量存在部分重疊。計算所有分量與原信號的對應互信息，結果如表4所示。

表4 各IMF 與含噪原始信號的互信息量Tab.4 Mutual information between each IMF and the noise-containing original signal

由表4可知，IMF1的互信息量為0.620 0，在所有分量中最大，說明IMF 包含原信號有效信息最多，這與實際波形也相符，而其他分量的互信息量相差不大，說明包含的原始信號有效信息較少。將各分量進行加權得到重構信號，再對重構信號進行ITD 分解并重構得到二次降噪信號，兩次重構信號的波形及頻譜如圖11和圖12所示。

圖11 基于加權分量一次重構信號及其頻譜Fig.11 First-reconstructed signal using weighted components and the signal’s spectrum

圖12 二次降噪重構信號及其頻譜Fig.12 Second denoised reconstructed signal and its spectrum

對比圖9、圖11及圖12可知，經(jīng)過ITD 二次降噪的信號其頻譜在高頻部分更加光滑，毛刺明顯減少，說明在改進VMD-ITD 能夠將信號的高頻噪聲濾除。為評價改進VMD-ITD 的降噪性能，通過計算其降噪后的信噪比和均方根誤差并與VMD、ITD 算法進行對比，得到3種算法對同一實測含噪信號降噪后的信噪比與均方根誤差結果對比如表5所示。

表5 3種算法降噪后信號的信噪比與均方根誤差Tab.5 Signal-to-noise ratio (SNR) and root mean square Error (RMSE) of the denoised signals using three algorithms

由表5可知，改進VMD-ITD（（IP）VMD-ITD）算法的均方根誤差為4.552 1，在3種方法中最小，且其信噪比為10.012 8，高于另外兩種方法，說明本文提出的改進VMD-ITD 方法在巖樣破裂聲發(fā)射信號降噪中的可行性與有效性，且相比于VMD、ITD 算法對噪聲抑制的能力更強。

4 結論

1）采用分量能量比從能量角度來確定VMD 的分解模態(tài)數(shù)k，并利用互信息法對分量進行加權重構，相比中心頻率法與相關系數(shù)法能更好地突出有效分量同時抑制噪聲分量，減少因剔除某些分量而造成的誤差。

2）運用改進VMD-ITD 聯(lián)合降噪算法抑制非平穩(wěn)信號中的噪聲相比單一的VMD、ITD 算法，其均方根誤差更小且信噪比更大，說明利用改進VMD-ITD對信號進行降噪能夠有效抑制噪聲，進一步提高信號的可靠性。

3）改進VMD-ITD 算法對巖樣聲發(fā)射破裂信號降噪處理能夠有效地抑制巖樣在破裂過程中的噪聲信號，提取出原始聲發(fā)射信號中的大部分有效信息，提高后續(xù)信號分析的有效性。