LEMD端點延拓方法及其在微震信號降噪中的應用

梁 喆，彭蘇萍，鄭 晶

(1.中國礦業大學(北京) 煤炭資源與安全開采國家重點實驗室，北京 100083;2.安徽理工大學 電氣與信息工程學院，安徽 淮南 232001)

微震監測技術在礦業工程、石油開采和防災減災監測方面已經獲得了廣泛的應用。微震事件發生時，采集儀器得到的是混有背景噪聲的信號。帶有噪聲的微震信號對定位精度有很大的影響，同時也對后續提取特征參數帶來誤差。為了能精確定位微震事件和計算信號特征參量，必須對微震監測系統記錄到的信號進行濾波、去噪。

李夕兵等［1－2］采用EMD方法對震動信號進行濾波、消噪。雖然EMD變換非常適合處理非平穩信號，但運用EMD方法對非平穩信號進行分解時，存在端點效應問題，造成整個數據序列分解后的結果嚴重失真，由此得到的分析結果的精度較低。

針對EMD的端點效應問題和實際應用，許多專家學者提出了不同的延拓方法，例如極值點鏡像延拓、包絡線延拓、局部波形延拓、波形匹配度延拓、支持向量機延拓、神經網絡延拓、AR模型延拓、最大Lyapunov指數邊界延拓方法，基于時間尺度的LS-SVM端點延拓以及最大相關波形延拓等［3－9］。以上延拓方法，對短序列數據都有很好的延拓精度，但是使用時需要大量的計算時間。微震信號是一種突變的非平穩信號，其自身是一個不確定的信號，沒有規律可循。對微震數據處理時，數據量大，且要求實時處理，在處理、解釋的過程中，檢測出突變發生的位置非常重要。采用上述方法對微震信號EMD分解時的端點效應處理時導致邊界誤差大，造成分解后突變信號衰減或丟失。

自適應濾波具有在未知環境下良好運行并跟蹤輸入統計量隨時間變化的能力，在沒有任何關于信號和噪聲的先驗知識的條件下，自適應濾波器利用前一時刻已獲得的濾波器參數來自動調節現時刻的濾波器參數，以適應信號和噪聲未知或隨機變化的統計特性。

結合 LMS自適應濾波器理論，本文對提出了LEMD(Empirical Mode Decomposition Based on LMS)端點效應處理方法。所提出的LEMD方法是利用LMS自適應濾波器根據端點以內數據所包含的信息，在數據序列兩端分別延拓出一段波形，使其符合數據的自然趨勢，從而產生EMD分解時邊界所需要的極大值和極小值。自適應濾波器通過自動地調節現時刻濾波器參數，適應微震信號和隨機噪聲隨時間變化的統計特性，同時根據數據段內部的信息和端點處數據的信息，實時修正延拓的結果。在整個EMD分解的過程中，實時跟蹤邊界信息，采用邊延拓邊分解的方法，不斷修正EMD篩分時插值算法帶來的誤差，最終實現最優延拓。

將LEMD方法應用于微震信號的降噪，成功獲取了微震信號定位和特征參數等信息，結果表明該方法是合理、有效的。

1 EMD分解及端點效應分析

在Hilbert-Huang中，為了計算瞬時頻率，定義了本征模態函數(IMF)，它必須滿足兩個條件:①在整個數據段內，其極值點的個數和過零點的個數必須相等或相差一個;②在任意時刻，由局部極大值點和極小值點形成的上、下包絡線平均值為零。對于不滿足IMF條件的復雜信號，可以用 EMD方法對其進行分解［10－11］。

EMD分解的基本思想是:對一給定信號，先獲得信號極值點，然后通過插值獲得信號上下包絡，計算出包絡均值，信號每減去1次包絡均值，便可得到1個固有模態函數。如此重復，直到將信號分解成有限個基本模式分量IMF和殘余項rn(t)的組合。

通過經驗模態分解，信號x(t)可表示為:

式中，imfi(t)表示具有不同頻率分量的IMF分量;rn稱為殘余函數，代表信號的平均趨勢。

對EMD分解后的不同IMF分量的取舍可以分別實現對信號的低通、帶通及高通濾波的功能［12］，根據式(1)可知，忽略趨勢項rn(t)，有

式中，x(t)為原始信號，X(t)為去噪后的信號。由式(2)可知，從原始信號中除去噪聲信號，就完成了噪聲信號的分離。

EMD是基于數據極值點進行分解的，對數據進行分析的時候，由于時間窗的選取，分析的數據往往是被截斷后的數據，這就造成信號在端點處不一定都是極值點，更不可能同時為極大值點和極小值點，致使利用樣條曲線插值求包絡運算的時候，在端點處形成較大的擺動，出現誤差，隨著分解過程的進行，重復進行插值求包絡線，導致誤差不斷累積并向內傳播污染數據，從而影響分解結果，這就是 EMD方法的端點效應問題。

微震信號的主頻率一般在幾赫茲到一百赫茲左右，信號在傳播的過程中，經過不同介質的吸收，其頻率總是在低頻段分布。對于處于低頻段微震信號，EMD分解時其極值點之間的距離較大，即時間跨度大，端部的邊緣效應更容易影響到信號內部，所以在用EMD對微震信號進行處理時必須抑制端點效應。

2 改善端點效應新方法

2.1 自適應濾波器(預測算法)基本原理［13］

在目前的研究及應用領域中，自適應濾波器可分為四類:系統識別、系統求逆、信號預測和干擾抵消。

假設輸入為某一特定矢量

自適應算法即根據第k個采樣周期所獲得的誤差e(kT)，計算下一個采樣周期計算式中的各個權系數:

式中，μ為步長因子。式(5)與式(6)的算法稱為LMS算法。這種算法對于每一個輸入樣本，只需對其進行式(7)中的兩個乘法與兩個加法運算，因此該算法易于用實時系統實現。

2.2 利用自適應預測算法對信號進行延拓

自適應濾波器的權系數向量迭代更新的過程，實質上是一種自學習和自訓練的過程。當使用LMS自適應濾波器進行端點延拓時，如果直接進行延拓，權系數向量會影響延拓的結果，造成結果不準確。在延拓之前可以首先對濾波器進行訓練，對權系數向量反復地優化，增強其在動態條件下的適應性。訓練結束后將訓練得到的權系數作為延拓的初始值，從端點處開始進行延拓。具體方法為:設一個數據序列內有n個數據{x(1)，x(2)，…，x(n)］，首先根據待延拓的數據點數m產生學習樣本，根據需要的目標輸出利用式(5)～式(7)對自適應濾波器進行訓練，得到對應的權系數Wj。然后用訓練好的權系數Wj作為延拓時候的初始值，將端點x(n)前的n－1個數據代入式(5)計算出端點處的第1個延拓值x^(n+1)。并以x^(n+1)為原數據序列新的端點，利用x^(n+1)為端點構成新的數據n個數據組合{x(2)，x(3)，…，x(m)，x^(n+1)］ 可以得到第2個數據序列延拓值x^(n+2)。以此類推，根據所需延拓的數據個數m可得到全部延拓序列x^(n+j)(j=1，2，…，m)。最后將新的延拓數據與原數據序列連接起來，就可以得到延拓后的全部數據序列。實際使用時數據點延拓個數可以根據需要的極值數確定，由延拓算法自適應決定。

綜上分析，LEMD算法的具體步驟如下:

(1)輸入數據序列;

(2)根據數據序列的長度，選取一部分數據按式(5)～式(7)對自適應濾波器訓練，訓練結束的自適應濾波器在數據序列的兩端延拓一段波形，波形的長度自適應決定，停止條件為兩邊各延拓出2個極大值和2個極小值;

(3)EMD分解，得到第一層分量，保留剩余信號中間對應延拓前的部分，采用自適應濾波器再次進行延拓，停止條件和第一層相同;

(4)每層分解結束后均對剩余信號再次延拓，直到整個分解結束。

LEMD算法采用延拓的數據完成EMD分解，信號包絡插值計算采用延拓后的全部數據，而在EMD分解完成后，舍棄兩端的延拓數據，只保留中間對應延拓前的部分，作為最終分解結果。

為了便于分析比較，采用文獻［14］中的仿真信號作為實例來分析，該仿真信號為:

原始信號和上下包絡線如圖1所示。為了和文獻［14］中的結果對比，采樣頻率設為200 Hz，數據點數為N=181。圖1中可見，該數據序列僅存在3個極大值點和4個極小值點，包絡線不僅沒有包含所有的數據，且上包絡線在數據的兩端出現了嚴重的失真。因此在采用EMD方法對信號進行分解時必須對端點進行處理。采用本文的方法，在信號的兩邊各延拓40個數據。首先對自適應濾波器進行訓練。以右邊的延拓為例:用信號的第1～160點的值按順序分為4組，每組40個數據，將1、3組作為學習樣本，2、4組作為期望輸出。也就是1組作為訓練輸入，2組作為期望輸出，通過構造的學習樣本和期望輸出對自適應濾波器訓練，訓練結束后，利用得到的權系數，將信號的第82～121點的值代入式(4)進行預測運算，得到第122點的值，再利用第83～122點的值預測得到第123點的值，以此類推預測得到第161點的值，即得到40個數據點。同樣，在左邊也可以得到40個數據點。延拓的波形和包絡線如圖2所示。

圖1 原始信號和包絡線Fig.1 original signal and the envelope

圖2 延拓后波形和包絡線Fig.2 Waveform and envelope after extension

圖2 (a)是對圖1所示的信號采用LEMD方法在端點處延拓后得到的信號，其中實線是對數據延拓后的信號波形，虛線為真實信號波形。雖然延拓的波形有稍許畸變，但和真實值相比，其極值的大小非常接近，誤差較小。圖2(b)為延拓后的包絡線，可以看出本文方法有效的抑制了包絡線在端點處的發散，所有信息都被包含在包絡線內。

為了比較LEMD方法的特點，將文獻［14］中的的幾種延拓方法性能比較的結果一并列出。從表1中可見，LEMD端點延拓方法與支持向量機、BP神經網絡相比，LEMD方法的學習時間僅需要0.025 8 s，誤差的均方值為0.143，LEMD方法的學習時間明顯比支持向量機和BP神經網絡都短很多。LEMD方法延拓的信號與原始信號誤差的均方值要比上述2種方法的稍大，但是本文采用的是邊分解邊延拓的方法，誤差不會累積，因此對數據精度沒有太高的影響。同時由于支持向量機不能預測太多的數據，BP網絡法訓練時間過大無法實時處理，而本文的LEMD方法則沒有這方面的限制，在保證系統精度的條件下數據延拓的實時跟蹤性強，更適合處理非平穩的突變信號。

表1 各種方法性能對比(仿真信號)Tab.1 Performance comparison of various methods(simulation signal)

采用LEMD方法對信號外延后進行分解得到的結果如圖3所示，其中虛線為分解后的IMF分量，實線為真實信號中的各分量，延拓后的IFM分量和真實信號中的各分量之間的相關系數分別為0.991 1、0.995 4、0.940 3，兩者非常接近。說明LEMD方法能夠有效地對數據序列進行延拓。

圖3 EMD分解圖Fig.3 Diagrams of EMD decomposition

3 實驗結果及其分析

選取某煤礦的一段現場監測微震信號數據來驗證本文方法的有效性，信號采樣頻率為1 kHz，采用ESG的采集儀，用地震檢波器拾取微震信號。現場監測是24 h監測，這里選用1個通道10 000個點的監測數據，截取中間1.5～5.596 s的一段，時域波形如圖4(a)所示，圖4(b)是其頻譜圖。

圖4 原始信號和頻譜Fig.4 Original signal and frequency spectrum

從圖4(a)中可見，由于背景噪聲的存在，從時域波形中無法獲取明顯地特征參數或定位信息，從圖4(b)中可觀測，微震信號特征頻率不明顯，整個頻譜圖中分布各種頻率的信號。

對微震信號進行分析，就是從數據序列中找出或分離出能夠反映定位或提取特征參數的信號。由于微震信號非平穩性強，且信號微弱，為了突出信號中的沖擊成分的特征，應該先消除噪聲，可以采用EMD方法對信號降噪，然后進行相應特征參數的提取。

圖5 延拓效果圖Fig.5 Extension effect diagram

為了對比本文方法的延拓效果，對截取的這一段數據向兩端各延拓100個數據點，以右邊延拓為例，在原數據從1.5～5.596 s這一段數據向右繼續取100點數據，即5.596～5.696 s時間段數據，用于和延拓后的結果進行對比。效果對比如圖5所示，其中虛線為原始信號，實線為基于 LEMD方法端點延拓后信號。LEMD方法在使用之前選取內部3 000點的數據對自適應濾波器進行訓練，訓練完成后進行數據延拓。

從圖5可以看出，5.671～5.696 s時間段數據屬于非平穩的突變信號，雖然延拓后的波形在幅值上和原數據有較大的誤差，但延拓后的數據與原數據極值點跳變規律基本保持一致。

為驗證LEMD方法抑制端點效應的優越性，以相同的分解終止條件，分別采用鏡像延拓法、LEMD方法對圖4的原始信號延拓并分解。圖5是該信號利用鏡像延拓抑制端點效應以后分解得到的各階IMF分量;圖6為采用LEMD抑制端點效應后分解得到的各階IMF分量。

圖6 采用鏡像延拓后的EMD分解圖Fig.6 EMD decomposition diagrams after mirror extension

圖7 采用LEMD端點延拓后的分解圖Fig.7 Decomposition diagrams after LEMD endpoint extension

通過圖6和圖7對比可以看出，利用鏡像延拓的方法和LEM方法來抑制經驗模態分解的端點效應后分解同樣的信號，得到相同的IMF分量，都是13個IMF分量和一個殘余分量。由圖6可以看出，imf8～imf13分量和殘余分量的右端點產生了較大的偏離，端點擺動現象很明顯。而圖7中經過LEMD處理端點效應以后，信號的端點效應得到了有效抑制，得到的IMF信號在端點處的偏離程度較小，尤其是信號的右端點更為明顯。imf8～imf11分量、imf13分量和殘余分量的右端點近似為0，分解得到的殘差更為合理。綜合分析比較，本文提出的LEMD算法分解精度更高，分解后的信號更加真實。

EMD分解后，對信號進行自適應低通濾波。利用EMD分解的結果重構信號，自適應的去掉信號高頻成分。圖8是采用鏡像延拓后的EMD分解重構后信號及頻譜圖，圖9為 LEMD方法分解重構后信號及頻譜圖。對比圖8(b)和圖9(b)可看出2種延拓方法分析結果的差異所在。盡管在圖8(b)和圖9(b)頻譜圖中均能明顯地看出微震信號，但圖9(b)噪聲抑制的效果高于圖8(b)。即基于本文方法對微震信號降噪的效果更好，對于定位信息的提取更加精確。

圖8 鏡像延拓后的EMD分解重構后信號及頻譜Fig.8 Reconstruction signal and frequency spectrum after Mirror extension

圖9 LEMD分解重構后信號及頻譜Fig.9 Reconstruction signal and frequency spectrum after LEMD

4 結論

采用EMD方法對微震信號降噪，由于數據選取時的截斷，造成EMD分解時引起端點效應問題，提出一種基于自適應濾波器延拓的LEMD端點延拓方法，該方法既考慮微震信號端點處的信號走向特點，又兼顧信號整體的發展趨勢，仿真結果表明該方法適應性強，有效延拓距離長，延拓實時性好。將該方法應用于微震信號降噪、分析上，在抑制端點效應的同時，較好的提取出反映微震特性的信號，為微震監控與定位提供數據支持，具有較高的實用價值。

［1］李夕兵，張義平，劉志祥，等.爆破震動信號的小波分析與HHT變換［J］.爆炸與沖擊，2005，25(6):528－535.LI Xi-bing，ZHANG Yi-ping，LIU Zhi-xiang，et al.Wavelet analysis and Hilbert-Huang Transform of blasting vibration signal［J］.Explosion and Shock Waves，2005，25(6):528－535.

［2］李夕兵，張義平，左宇軍，等.巖石爆破振動信號的EMD濾波與消噪［J］.中南大學學報(自然科學版)，2006，37(1):150－154.LI Xi-Bing，ZHANG Yi-ping，ZUO Yu-Jun，et al.Filtering and denoising of rock blasting vibration signal with EMD［J］.Journal of Central South University (Science and Technology)，2006，37(1):150 －154.

［3］胡勁松，楊世錫.EMD方法基于AR模型預測的數據延拓與應用［J］.振動、測試與診斷，2007，27(2):116－119.HU Jing-song，YANG Shi-xi.AR model prediction-based EMD method and its application to data extension［J］.Journal of Vibration，Measurement ＆ Diagnosis，2007，27(2):116－119.

［4］鄧擁軍，王偉，錢成春，等.EMD方法及Hilbert變換中邊界問題的處理［J］.科學通報，2001，46(3):257－263.DENG Yong-jun，WANG Wei，QIAN Cheng-chun，et al.Boundary processing technique in EMD method and Hilbert transfor m［J］.Chinese Science Bulletin，2001，46(3):257－263.

［5］張郁山，梁建文，胡聿賢.應用自回歸模型處理EMD方法中的邊界問題［J］.自然科學進展，2003，13(10):1054－1059.ZHANG Yu-shan，LIANG Jian-wen，HU Yu-xian.The processing of end effects in EMDmethod by autoregressive model［J］.Progresses in Nature Science，2003，13(10):1054－1059.

［6］程軍圣，于德介，楊宇.基于支持矢量回歸機的Hilbert-Huang變換端點效應問題的處理方法［J］.機械工程學報，2006，42(4):23 －31.CHENG Jun-sheng，YU De-jie，YANG Yu.Process method for end-effects of Hilbert-Huang transform based on support vector regression machines［J］. Chinese Journal of Mechanical Engineering，2006，42(4):23 －31.

［7］蔡艷平，李艾華，張瑋，等.HHT端點效應的最大Lyapunov指數邊界延拓方法［J］.儀器儀表學報，2011，32(6):1330－1336.CAI Yan-ping，LI Ai-hua，ZHANG Wei，et al.HHT end effect processing method based on maximum Lyapunov index boundary extension model［J］.Chinese Journal of Scientific Instrument，2011，32(6):1330 －1336.

［8］郭明威，倪世宏，朱家海，等.振動信號中 HHT/EMD端點延拓方法研究［J］.振動與沖擊，2012，31(8):62－66.GUO Ming-wei，NI Shi-hong，ZHU Jia-hai，et al.HHT/EMD end extension method in vibration signal analysis［J］.Journal of Vibration and Shock，2012，31(8):62 －66.

［9］高強，段晨東，趙艷青，等.基于最大相關波形延拓的經驗模式分解端點效應抑制方法［J］.振動與沖擊，2013，32(2):62－66.GAO Qiang，DUAN Chen-dong，ZHAO Yan-qing，et al.A maximal correlation waveform extension method for end effects reduction of empirical mode decomposition［J］.Journal of Vibration and Shock，2013，32(2):62 －66.

［10］ Huang N E，Shen Z，Steven R L，et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonsta-tionary time series analysis［J］.Proceedings of the.Royal Society A Mathematical，Physical and Engineering Sciences，1998，454:899 －995.

［11］ Huang N E，Wu M L C，Long SR，et al.A confidence limit for the empirical mode decomposition and Hilbert spectral analysis［J］.Proceedings of the Royal Society.Series A:Mathematical，Physical and Engineering Sciences，2003，459:2317－ 2345.

［12］譚善文，秦樹人，湯寶平.Hilbert-Huang變換的濾波特性及其應用［J］.重慶大學學報，2004，27(2):9－12.TAN Shan-wen，QIN Shu-ren，TANG Bao-ping.The filtering character of Hilbert-Huang transform and its application［J］.Journal of Chongqing University，2004，27(2):9 －12.

［13］高晉占.微弱信號檢測(第二版)［M］.北京:清華大學出版社.

［14］蘇玉香，劉志剛，李科亮，等.一種改善EMD端點效應的新方法及其在諧波分析中的應用［J］.電工電能新技術，2008，27(2):33－37.SU Yu-xiang，LIU Zhi-gang，LI Ke-liang，et al.A new method for end effect of EMD and its application to harmonic analysis［J］.Advanced Technology of Electrical Engineering and Energy，2008，27(2):33－37.