郭航伸，郭連軍,2，楊 巍，潘 博，徐振洋

(1.遼寧科技大學(xué)礦業(yè)工程學(xué)院，遼寧 鞍山 114051；2.沈陽工業(yè)大學(xué)建筑與土木工程學(xué)院，沈陽 110870；3.成遠(yuǎn)礦業(yè)開發(fā)股份有限公司，遼寧 遼陽 111000；4.北京科技大學(xué)土木與資源工程學(xué)院，北京 100083)

在爆破作業(yè)環(huán)境采集的振動信號通常會含有大量噪聲信號，且噪聲來源復(fù)雜多樣。目前，常見的振動信號去噪方法有濾波法[1]、形態(tài)學(xué)法[2]、小波法[3-4]、經(jīng)驗?zāi)B(tài)分解法[5]等。但是這些方法的去噪效果一般且不易精確重構(gòu)原始爆破振動信號。小波技術(shù)具有多分辨分析特性和良好的時頻局部性，原始信號分解后，真實信號和噪聲可依據(jù)小波系數(shù)的不同特性實現(xiàn)分離，被廣泛應(yīng)用于爆破振動信號去噪領(lǐng)域[6-7]。局部均值分解法[8](LMD)是2005年由Jonathan Smith提出的一種自適應(yīng)的信號時頻分析方法，可以將一個復(fù)雜的非平穩(wěn)信號分解成若干個PF(Product Function)分量和殘余量之和的形式，每一個PF分量由其對應(yīng)層的包絡(luò)函數(shù)和純調(diào)頻函數(shù)相乘得到，所有PF分量的瞬時頻率和瞬時幅值疊加便得到了原始信號完整的時頻分布。LMD方法由于其良好的分析非線性、非平穩(wěn)信號的能力，已在腦電信號分析[9-10]和機械故障診斷領(lǐng)域[11-12]得到應(yīng)用。LMD法雖然比經(jīng)驗?zāi)B(tài)分解法更為溫和，但同樣存在模態(tài)混疊和端點效應(yīng)問題[13]。

通過引入邊界局部特征尺度匹配延拓算法[14]、自適應(yīng)互補集合算法與小波閾值法對局部均值分解算法進(jìn)行改進(jìn)，以實現(xiàn)對爆破振動信號的有效降噪。

1 自適應(yīng)互補集合局部均值分解算法

1.1 總體平均局部均值分解方法

總體平均局部均值分解方法[15]是對局部均值分解方法的改進(jìn)，通過添加輔助噪聲進(jìn)行數(shù)據(jù)分析，向原始信號中添加n次等幅值隨機白噪聲后進(jìn)行LMD分解，加入的白噪聲使得數(shù)據(jù)極值點的分布更加均勻，進(jìn)而改善了LMD算法中存在的模態(tài)混疊問題，同時利用多次隨機白噪聲的均值接近于0的特性消除引入白噪聲的影響。

1.2 自適應(yīng)互補集合局部均值分解算法

自適應(yīng)互補集合局部均值分解算法(CELMDAN)是對ELMD算法的改進(jìn)。ELMD算法中加入了多次白噪聲，最終數(shù)據(jù)重構(gòu)誤差取決于加入白噪聲的次數(shù)，次數(shù)越高則重構(gòu)誤差越小，但是程序的運行時間更長。CELMDAN算法通過添加有限次數(shù)的自適應(yīng)白噪聲，使信號能在加入較少次的噪聲下達(dá)到重構(gòu)誤差為0，所以CELMDAN不但可以降低LMD算法中存在的模態(tài)混疊狀態(tài)，還具有重構(gòu)誤差小，計算效率高的優(yōu)點。其算法的具體流程如下：

在原始信號x(t)中加入N個不同的白噪聲，進(jìn)行LMD分解，對獲得的PF分量取平均值，得到第1個模態(tài)分量，記為

(1)

第1個唯一的剩余分量記做r1(t)，記為

r1(t)=x(t)-PF1(t)

(2)

定義Ek(·)為信號經(jīng)過LMD分解產(chǎn)生的第k個PF分量，則第2個模態(tài)分量如下

(3)

第k個剩余分量記為

rk(t)=rk-1(t)-PFk(t)

(4)

第k+1個模態(tài)分量記為

(5)

執(zhí)行公式(4)和公式(5)，直到余量信號不能被分解時停止。分解終止時，信號可以表示為

(6)

式中：R(t)為最后的剩余分量。

從上述的算法中可以看出，CELMDAN的算法分解過程是完整的，能夠?qū)π盘栠M(jìn)行精確重構(gòu)，且在分解的過程中能通過系數(shù)εk選擇合適的噪聲。

1.3 端點效應(yīng)的處理方法

LMD算法也存在端點效應(yīng)，算法中由于振動信號兩側(cè)的極值是未知的，在插值擬合的過程中會出現(xiàn)較大的誤差，進(jìn)而出現(xiàn)端點效應(yīng)。針對LMD算法中的端點效應(yīng)問題，相關(guān)學(xué)者提出了一些改進(jìn)方法[16-17]。

邊界局部特征尺度匹配延拓法可以有效地解決信號兩端出現(xiàn)的端點效應(yīng)。設(shè)Mi、Ni分別為信號的極大值和極小值，對應(yīng)的時間點為tmi、tni，S1為信號左邊界處的點。

首先以S1-M1-N1為三角形，找到3個頂點對應(yīng)的時間坐標(biāo)，找到S1對應(yīng)的點Si，Si對應(yīng)的時間點tsi由下式可得：

(7)

若求得的tsi不在原始信號的波形上，采用插值的方法求出Si。然后計算匹配誤差Ei：

Ei=|Si-S1|+|Ni-N1|+|Mi-M1|

+|Mi+1-M2|

(8)

當(dāng)min(Ei)<βl時，認(rèn)為信號的規(guī)律性較強。其中β為常數(shù)，β取值既不能太大，也不能太??；l為S1-M1-N1構(gòu)成的三角形數(shù)據(jù)長度。若計算出的匹配誤差在所規(guī)定的閾值內(nèi)，取所有計算出的匹配誤差中的最小值；當(dāng)計算出的匹配誤差有兩個或兩個以上時，取與起始點距離最遠(yuǎn)的波形。當(dāng)匹配誤差不滿足條件范圍時，則認(rèn)為信號規(guī)律性不強。此時根據(jù)信號的端點極值點的特征進(jìn)行延拓，對信號的另一端采用同樣的方法進(jìn)行計算，兩側(cè)的計算完成后，就可以采用完成后的數(shù)據(jù)進(jìn)行分解。以上過程可以看出，經(jīng)過延拓方法處理的數(shù)據(jù)與原始信號具有很好的一致性。

2 小波閾值降噪方法

在爆破振動信號去噪中，小波閾值去噪有較為廣泛的應(yīng)用。小波閾值去噪方法是將爆破振動的時域信號轉(zhuǎn)換到小波域，在小波域進(jìn)行閾值處理。由于隨機噪聲具有較小的小波系數(shù)，可以設(shè)置相應(yīng)的閾值將其去除，再把處理后的剩余信號重構(gòu)，可以獲得去噪后的信號。小波閾值去噪效果的優(yōu)劣與所選的小波基、分解層數(shù)、閾值估計方法和閾值函數(shù)有關(guān)。閾值選取法主要包括默認(rèn)法、硬閾值法和軟閾值法，默認(rèn)法適應(yīng)性差，可信度較低；硬閾值保留較大的小波系數(shù)，能夠保留信號有效部分的特征；軟閾值將大于閾值的系數(shù)向零收縮，去噪效果更為平滑。軟閾值的具體處理過程如下：

(9)

式中：xs(i,j)為去噪后的小波系數(shù)；x(i,j)為去噪前的小波系數(shù)；thr(i)是由各類閾值準(zhǔn)則確定的閾值。

3 基于CELMDAN和小波閾值的降噪方法

單獨使用小波閾值降噪方法雖然能夠消除爆破振動信號的噪聲，但是同時也消除了一部分有用信息。單獨使用CELMDAN對爆破振動信號進(jìn)行強制降噪時，通常的做法是直接去除前幾個PF分量，這樣雖然去掉了高頻噪聲，但是前幾個PF分量也包含了爆破振動信號內(nèi)在的高頻信息。為了最大程度去除信號噪聲同時盡可能保留信號中的有用信息，提出基于CELMDAN和小波閾值相結(jié)合的降噪方法。

首先采用CELMDAN算法將爆破振動信號分解為一系列的PF分量，然后分別計算各PF分量的自相關(guān)系數(shù)，根據(jù)自相關(guān)系數(shù)判斷包含隨機噪聲的PF分量，對由噪聲主導(dǎo)的PF分量利用小波閾值降噪方法進(jìn)行降噪處理，最后將小波閾值降噪后的PF分量與剩余PF分量重構(gòu)得到降噪后的爆破振動信號。

4 仿真分析

對提出的基于CELMDAN和小波閾值去噪方法進(jìn)行仿真測試。測試中，采用兩組正余弦波相加來模擬爆破振動信號，兩組正余弦波分別為2sin (100πt)和cos (50πt)，頻率分別為50 Hz和25 Hz，采樣頻率為2 048 Hz，采樣時間為2 s(見圖1)。在仿真信號數(shù)據(jù)中加入信噪比(SNR)為2的白噪聲，得到假設(shè)的含噪爆破振動信號(見圖2)。

圖1 仿真信號Fig.1 Original signal

圖2 含噪仿真信號Fig.2 Vibration signal containing noise

將含噪仿真爆破振動信號進(jìn)行CELMDAN分解，共得到11個PF分量以及1個趨勢項，PF分量如圖3所示。分別計算各個PF分量的自相關(guān)系數(shù)，當(dāng)某個PF的自相關(guān)系數(shù)較小時，則認(rèn)為該分量的主要成分為噪聲，將判斷為噪聲PF分量再進(jìn)行小波閾值降噪處理。降噪后的結(jié)果如圖4所示。

圖3 CELMDAN分解結(jié)果Fig.3 CELMDAN decomposition results

圖4 降噪后信號Fig.4 Signal after noise reduction

為驗證CELMDAN-WT方法對比其他方法的優(yōu)勢，分別采用小波閾值直接降噪方法、EMD分解、ELMD分解去除第1個分量降噪法對仿真信號降噪處理。采用均方根誤差(RMSE)和信噪比(SNR)兩項指標(biāo)對算法的性能進(jìn)行檢驗，均方根誤差和信噪比定義如下[18]：

(10)

(11)

表1 不同去噪方法的性能

對比4種方法的性能指標(biāo)可以得出：單獨使用小波閾值去噪方法在去除噪聲時也帶來嚴(yán)重的信號失真；EMD分解后單獨舍棄IMF1分量，由于噪聲信號去除不完全，導(dǎo)致信噪比不高，信號局部不光滑；總體平均局部均值分解和經(jīng)驗?zāi)B(tài)分解相比較而言，一定程度上降低了模態(tài)混疊現(xiàn)象，去除了更多的噪聲，相比于前兩種方法信噪比提升很大；而采用CELMDAN-WT去噪方法處理信號時可以得到更高的信噪比和更小的均方根誤差，信號獲得了良好的降噪效果。

5 應(yīng)用實例

以某工程現(xiàn)場實測爆破振動信號為例，爆破振動信號來自鞍鋼關(guān)寶山露天鐵礦爆破現(xiàn)場，該區(qū)域地質(zhì)條件較差，測點周邊臨近省道、村莊、民房，由于測試環(huán)境復(fù)雜，獲取的爆破振動信號受噪聲污染嚴(yán)重。監(jiān)測儀器為成都中科測控有限公司生產(chǎn)的TC-4850測振儀，測量范圍：0～35 cm/s；A/D分辨率：16 Bit；頻響范圍：0～10 kHz；測量精度：0.01 cm/s。隨機選取了在不同時間內(nèi)測得的兩組數(shù)據(jù)進(jìn)行處理，原始數(shù)據(jù)如圖5所示。原始爆破振動信號波形復(fù)雜、受噪聲污染嚴(yán)重，需要去噪處理獲得真實信號。采用CELMDAN-WT方法對所測信號進(jìn)行去噪處理，去噪后信號如圖6所示。

圖5 現(xiàn)場監(jiān)測原始信號Fig.5 Signals measured in situ

圖6 去噪后信號Fig.6 Denoised signal

由圖6可以明顯看出，經(jīng)過CELMDAN-WT方法去噪后的爆破振動信號波形較去噪前清晰，計算含噪爆破振動信號去噪后的信噪比、均方根誤差，計算結(jié)果如表2所示。

表2 爆破振動信號去噪效果

從表2中可以看出，去噪后爆破振動信號的信噪比都有不同程度的提升，受噪聲干擾最嚴(yán)重的是第1組信號，信噪比提升了9.2 dB。對去噪前后的信號做FFT時頻譜分析，去噪前后頻域如圖7所示。

圖7 去噪前后信號頻域Fig.7 Frequency domain of signal before and after denoising

從圖7可以看出，原始含噪信號的頻譜相對雜亂。100 Hz以下的低頻部分是爆破振動的優(yōu)勢頻段，但其在頻譜圖中優(yōu)勢不明顯，說明信號受到了嚴(yán)重的噪聲干擾，經(jīng)過去噪處理后的信號頻譜優(yōu)勢頻段更加明顯。經(jīng)過CELMDAN-WT方法處理后的爆破振動信號，信噪比有了較大的提升，信號的原始信息得到了更好的保留，為下一步的爆破振動信號的分析提供了可信的基礎(chǔ)。

由原始爆破振動信號三維能量(見圖8)可以看出，爆破振動信號能量主要集中在0～100 Hz，在100～200 Hz范圍內(nèi)分布著大量的低能量信號，這些低能量信號主要是采樣環(huán)境中的其他噪聲信號。通過與采用CELMDAN-WT方法去噪后爆破振動信號的三維能量(見圖9)，對比可以看出，原始爆破振動信號經(jīng)過去噪后，100～200 Hz內(nèi)噪聲信號基本完全除去。由于噪聲信號的存在，信號優(yōu)勢頻段在三維能量圖中表現(xiàn)不明顯，而去噪后則可以明顯看出信號有用成分能量在時域和頻域的分布。

圖8 原始爆破振動信號三維能量Fig.8 Three dimensional energy of original blasting vibration signal

圖9 去噪后爆破振動信號三維能量Fig.9 Three dimensional energy of blasting vibration signal after denoising

6 結(jié)論

1)對含有噪聲的爆破振動信號直接使用小波閾值法降噪效果不佳，信號嚴(yán)重失真。

2)EMD和ELMD方法分解后直接剔除IMF1分量雖然能夠去除一部分噪聲，但因為存在模態(tài)混疊的問題會減少信號高頻部分的有用成分。

3)CELMDAN-WT方法能夠有效地消除爆破振動信號中的噪聲，實現(xiàn)了信號與噪聲的分離，信噪比最高提升9.2 dB。