基于HHT法的煤沖擊破壞低頻電磁信號去噪

付玉凱，楊 威，李成武

(1.天地科技股份有限公司開采設(shè)計(jì)事業(yè)部，北京100013;2.煤炭科學(xué)研究總院開采設(shè)計(jì)研究分院，北京100013;3.煤炭科學(xué)研究總院煤炭資源高效開采與潔凈利用國家重點(diǎn)實(shí)驗(yàn)室，北京100013;4.中國礦業(yè)大學(xué) (北京)資源與安全工程學(xué)院，北京100083)

隨著煤礦開采深度的增加，煤礦煤巖體動力災(zāi)害日益增多，嚴(yán)重影響了煤礦的安全、高效生產(chǎn)。針對煤巖體動力災(zāi)害，預(yù)測預(yù)報(bào)是關(guān)鍵。目前，電磁輻射預(yù)測方法受到了國內(nèi)外學(xué)者的關(guān)注［1－3］。

由于煤巖沖擊破壞過程中產(chǎn)生的電磁信號是一種非線性、陣發(fā)性的脈沖信號，屬于典型的非平穩(wěn)隨機(jī)信號［4］，所以對其進(jìn)行去噪濾波顯得非常困難。因此，去噪濾波技術(shù)嚴(yán)重制約了電磁信號在預(yù)測預(yù)報(bào)煤巖動力災(zāi)害中的應(yīng)用。目前，學(xué)者主要采用傅里葉變換 (FFT)、小波變換 (WT)等信號處理方法對電磁信號進(jìn)行分析，但是這些方法只能分析信號的總體頻率，不能有效分析信號的時(shí)頻特性［5－6］。小波變換分析方法與其他分析方法相比，其時(shí)頻分析精度依賴于小波基函數(shù)的選取，有時(shí)由于基函數(shù)選取問題，使其對信號的精細(xì)分解受到限制［7－8］。

希爾伯特黃變換 (HHT)是1998年由Huang等人［9－10］提出了一種處理非線性、非平穩(wěn)信號的時(shí)頻分析方法，HHT分析方法主要包括2個(gè)部分:一是多分辨經(jīng)驗(yàn)?zāi)B(tài)分解 (EMD)和瞬時(shí)頻率變換;二是對EMD分解的分量進(jìn)行時(shí)頻分析。HHT變換實(shí)際上是一種以傅里葉變換為基礎(chǔ)的改進(jìn)型信號處理方法，該方法在電磁信號去噪中的應(yīng)用較少。

本文采用HHT時(shí)頻分析方法對煤沖擊破壞的低頻電磁信號進(jìn)行經(jīng)驗(yàn)?zāi)Ｊ椒纸?(EMD)，并通過對各個(gè)IMF(某一頻率尺度上的模態(tài)信號)分量進(jìn)行重構(gòu)，最后對重構(gòu)信號進(jìn)行時(shí)頻譜分析。

1 HHT變換原理

1.1 經(jīng)驗(yàn)?zāi)B(tài)分解法 (EMD)原理及算法

1.1.1 EMD原理

EMD法［11］是 HHT的核心，其主要有兩個(gè)作用:過濾疊加波和對稱化波形。EMD可以將信號從時(shí)間尺度上進(jìn)行IMF分量分解，但是需要滿足下面兩個(gè)條件:

(1)整個(gè)信號數(shù)據(jù)的極值點(diǎn)和過零點(diǎn)個(gè)數(shù)相差不超過1。

(2)由局部最大值所繪制的包絡(luò)線和由局部最小值點(diǎn)所繪制的包絡(luò)線的平均值為0，就是信號必須對稱于時(shí)間軸。

1.1.2 EMD算法

對一個(gè)原始信號X(t)，首先找出X()t上全部的最大值和最小值，然后采用插值方法對曲線進(jìn)行極值擬合，從而繪制出曲線的包絡(luò)線Xmax()t。同理得出包絡(luò)線Xmin()t，2條包絡(luò)線包含了全部信號數(shù)據(jù)。對2條包絡(luò)線取其平均值得平均線m1()t，再用X()t減掉m1()t得到h1()t。如果信號不同，h1()t有可能產(chǎn)生一個(gè)IMF分量，也可能得到2個(gè)IMF分量。若分量不滿足限定條件，此時(shí)將h1()t當(dāng)做原信號，重復(fù)以上的程序，即得h11()t =h1()t －m11()t，m11()t是h1()t的2條包絡(luò)線平均值;若h11()t沒有變換成IMF分量，則接著計(jì)算，進(jìn)行k次計(jì)算，得到第k次計(jì)算的數(shù)據(jù)h1k(t)=h1(k－1)(t)－m1k(t)。判斷h1k(t)是不是可以滿足IMF分量，這樣就需要一個(gè)計(jì)算過程終止的法則，一般可以用2個(gè)連續(xù)結(jié)果的標(biāo)準(zhǔn)差SD作為判斷依據(jù)［12－13］:

在實(shí)際運(yùn)算中，可以通過對信號反復(fù)篩選來確定IMF分量，篩選結(jié)果通過SD值來確定。經(jīng)驗(yàn)表明，當(dāng)取SD=0.2～0.3時(shí)比較合適，既可確保IMF的線性和穩(wěn)定性，又可使IMF具有相應(yīng)的物理意義［11］。

1.2 HHT變換與時(shí)頻譜

1.2.1 HHT 變換

對信號進(jìn)行分解，得到IMF分量，對分量進(jìn)行HHT變換，計(jì)算出其瞬時(shí)頻率。對全部IMF分量進(jìn)行上述變換，即HHT譜［10］。HHT變換很好地刻畫了信號的局部性質(zhì)，可以很好地得到信號的瞬時(shí)頻率，避免了傅里葉變換中產(chǎn)生的不真實(shí)高低頻成分，其具有直觀的物理意義［10］。

1.2.2 HHT譜

HHT譜變換［14］可以把信號幅值變換為時(shí)間、頻率平面上的等高線圖，該變換稱之為HHT時(shí)頻譜。時(shí)頻譜主要有3種表達(dá)形式:灰度圖、等高線及三維空間圖形，其表達(dá)式如下:

式中，H為信號幅值;Re為累計(jì)相加;ai()t為每一階IMF的幅值;ωi(t)為每一階IMF的瞬時(shí)頻率。

如果對Hω，()t進(jìn)行時(shí)間上積分，可以得到信號的HHT邊際譜:

2 實(shí)驗(yàn)簡介

2.1 實(shí)驗(yàn)系統(tǒng)及裝置

試驗(yàn)系統(tǒng)包括兩部分，即霍普金森壓桿(SHPB)和電磁輻射測試系統(tǒng)，試驗(yàn)系統(tǒng)見圖1，圖2。

圖1 實(shí)驗(yàn)裝置

圖2 實(shí)驗(yàn)裝置實(shí)物

霍普金森壓桿系統(tǒng)由子彈、輸入桿和輸出桿組成，壓桿為鋼質(zhì)壓桿，直徑為50mm，子彈為φ50mm×400mm的圓柱體。被測試樣夾在輸入桿和輸出桿之間。

選用的電磁輻射接收裝置為ZDKT－1型瞬變磁振測試系統(tǒng) (圖2)，該系統(tǒng)包括磁場天線、信號采集系統(tǒng) (3000s－1)及計(jì)算機(jī)。實(shí)驗(yàn)時(shí)，天線正對煤試件，距其30～40mm。

2.2 煤樣制作

煤試件來源于某礦掘進(jìn)頭，煤樣采用圓柱體，尺寸為 φ50mm ×50mm，平行度0.02mm［15］，試樣兩端涂抹石墨，以減少其摩擦效應(yīng)［16－17］。

2.3 煤沖擊破壞低頻電磁信號去噪前分析

共加工12個(gè)試樣，分為4組進(jìn)行實(shí)驗(yàn)，分別記為 A1，A2，A3;B1，B2，B3;C1，C2，C3;D1，D2，D3。對每一組進(jìn)行相同速率下的實(shí)驗(yàn)，由于煤體強(qiáng)度較低，經(jīng)多次沖擊測試發(fā)現(xiàn)，沖擊速率大于3m/s即可破壞煤試樣，且沖擊速率過大會導(dǎo)致應(yīng)力－應(yīng)變曲線失真。當(dāng)沖擊速率在3～10m/s之間時(shí)，測試結(jié)果可靠性較高。由于子彈沖擊速率是由動力系統(tǒng)中的高壓氮?dú)馑┘拥模渌俾士刂朴幸欢ǖ恼`差，所以沖擊加載速率以平行光源測試結(jié)果為準(zhǔn)。根據(jù)平行光源測定結(jié)果，實(shí)驗(yàn)的沖擊速 率 分 別 為 3.287m/s，6.251m/s，6.950m/s，8.714m/s。實(shí)驗(yàn)結(jié)果發(fā)現(xiàn)，同一組實(shí)驗(yàn)結(jié)果的重復(fù)性較好。鑒于文章篇幅有限，以D1為典型信號進(jìn)行分析，沖擊速率為8.714 m/s時(shí)，最大應(yīng)變率為166.35s－1，采集的低頻磁場原始信號見圖3。

圖3 原始信號

由圖3可以看出，共采集了10s的低頻磁場信號，但突變的低頻磁場信號介于5～7s之間，持續(xù)時(shí)間較短 (小于2s)，并且信號中伴隨著大量的背景噪聲信號，這些噪聲信號主要來自外界環(huán)境和采集系統(tǒng)自身［18］。存在噪聲的低頻電磁信號對于預(yù)測煤巖沖擊破壞十分不利。為了能清楚地認(rèn)識低頻磁場信號的特征，需要對原始低頻電磁信號進(jìn)行去噪分析。為了驗(yàn)證HHT分析煤沖擊破壞的低頻磁場信號的有效性和突顯信號非線性的能力，首先將原始低頻磁場信號進(jìn)行FFT頻譜分析和Morlet時(shí)頻分析［19－20］。Morlet時(shí)頻分析是小波變換的一種形式，分析結(jié)果見圖4和圖5。

由圖4和圖5可以看出，信號的能量主要集中在600Hz以內(nèi)，原始低頻信號中存在明顯的噪聲成分，并且Morlet的時(shí)頻分析譜的有效信號也被背景噪聲信號掩蓋，從上面2個(gè)圖很難清楚認(rèn)識有效信號隨時(shí)域和頻域的動態(tài)變化特征，這就需要進(jìn)一步采用HHT法對信號進(jìn)行去噪分析。

圖4 原始信號的FFT譜

圖5 原始信號的Morlet時(shí)頻譜

3 低頻磁場信號的HHT分析

3.1 磁場信號的EMD分解及重構(gòu)

低頻電磁信號屬于非平穩(wěn)脈沖信號，在煤巖沖擊破壞過程中有時(shí)某一時(shí)間段的信號強(qiáng)度大，那么信噪比就較高;而另一時(shí)間段信號強(qiáng)度較弱，這時(shí)信噪比就很低。如果采集到的信號強(qiáng)度一直較低，那么會嚴(yán)重影響信號的采集質(zhì)量［21］。如果在信號分析時(shí)，我們能選擇能量強(qiáng)的信號進(jìn)行分析，而舍棄那些能量弱的信號 (信號強(qiáng)弱跟煤巖體破裂程度相關(guān))，這樣才能得到原始信號中的有效信號。其他信號處理方法難以完成上述問題，這也是HHT方法能完成這一問題的關(guān)鍵。

首先運(yùn)用HHT中的經(jīng)驗(yàn)?zāi)B(tài)分解方法，即EMD分解，可以得到原始低頻磁場信號的有限數(shù)目的IMF分量。由于信號的特征尺度參數(shù)都是基于所采集信號，所以，EMD所篩分出來的IMF分量都具有實(shí)際物理意義，每個(gè)IMF分量代表了某一頻率尺度上的模態(tài)。圖6是低頻電磁信號的EMD分解圖，該信號共分為11階IMF分量，在時(shí)間域上表示從小尺度到大尺度的層層篩選濾波。

從圖6可以看出，IMF_h1～I(xiàn)MF_h5分量包含大量的噪聲信號 (外界環(huán)境噪聲和采集系統(tǒng)噪聲);IMF_h6～I(xiàn)MF_h10分量噪聲信號較低，其屬于低頻磁場信號的優(yōu)勢分量，為信號優(yōu)勢頻率分量;IMF_residual分量表示磁場信號的變化趨勢，通常稱之為殘余分量。

將IMF_h6～I(xiàn)MF_h10分量進(jìn)行重構(gòu)，重構(gòu)信號見圖7所示。由圖7可以看出，重構(gòu)的信號能很清晰地刻畫出低頻磁場信號的非線性、非平穩(wěn)性，很好地表示出了信號的時(shí)頻特征——磁場信號初期線性上升，然后指數(shù)衰減，最后尾部小幅震蕩［22－23］。

圖6 原始信號EMD分解

圖7 去噪重構(gòu)信號

將重構(gòu)的信號進(jìn)行FFT變換，得到信號的幅值－頻率圖 (圖8)。由圖8可以看出，煤沖擊破壞產(chǎn)生的低頻磁場信號的優(yōu)勢頻率較低，范圍在0～40Hz之間，對另外幾組煤樣進(jìn)行實(shí)驗(yàn)。采用上述分析方法后，分析結(jié)果基本相同。

圖8 重構(gòu)信號的頻域

3.2 重構(gòu)信號的時(shí)頻譜分析

圖9 (a)是對重構(gòu)信號進(jìn)行的HHT時(shí)頻譜，從圖中可以看出，顏色越亮表示其能量越高，反之越低;HHT時(shí)頻譜直觀地表現(xiàn)出了信號的聚集性，與原始信號的時(shí)頻譜 (圖5)相比，可以很清晰地看出能量隨時(shí)間和頻率的變化。顯然，在0～100Hz頻段內(nèi)，整個(gè)時(shí)域上的能量都較強(qiáng);在時(shí)間域上，第5s至第6s時(shí)間域上的能量較強(qiáng)，該時(shí)間域?qū)?yīng)于煤的沖擊破壞時(shí)間。

圖9(b)表示重構(gòu)信號的邊際時(shí)間，表示每個(gè)時(shí)間在全局上的幅度之和，是表示時(shí)間點(diǎn)上信號能量強(qiáng)弱的一個(gè)指標(biāo)，從該圖中也可以看出，第5s至第6s時(shí)間域上的能量較強(qiáng)，很好地驗(yàn)證了HHT時(shí)頻譜分析的正確性。

圖9(c)表示重構(gòu)信號的邊際譜，表示在每個(gè)頻域上全部幅值的總和，由圖可以看出，0～40Hz范圍的邊際譜比較大，這也說明了信號能量的主要頻域集中在0～40Hz，分析結(jié)果與圖8相吻合。

圖9 重構(gòu)信號的HHT時(shí)頻譜、邊際時(shí)間和邊際譜

總之，HHT方法可以有效地對煤沖擊破壞的低頻磁場信號進(jìn)行濾波去噪，可以很好地表示磁場信號的時(shí)頻動態(tài)非線性、非平穩(wěn)性及脈沖特性，可以表現(xiàn)出信號時(shí)域上的頻率和能量差異，與其他信號處理方法相比，有其獨(dú)特的優(yōu)勢。由上述分析可以看出，HHT時(shí)頻分析方法為煤巖體沖擊破壞磁場信號的濾波處理提供了一種新的方法。

4 結(jié)論

(1)采集到的煤沖擊破壞低頻電磁信號持續(xù)時(shí)間較短 (1～2s)，并且信號中存在大量的背景噪聲。通過對原始信號進(jìn)行FFT頻譜分析和Morlet時(shí)頻分析，得出原始信號的能量主要集中在600Hz以內(nèi);并且Morlet時(shí)頻分析譜的有效信號被背景噪聲信號完全掩蓋。

(2)HHT法中的EMD可以很好地分解原始磁場信號，然后得出其IMF分量，每個(gè)IMF分量都有其特定的物理意義，通過對有效IMF分量進(jìn)行重構(gòu)，重構(gòu)信號能很好地刻畫出低頻電磁信號的非線性、非平穩(wěn)性及脈沖特性。

(3)通過對重構(gòu)信號進(jìn)行時(shí)頻譜分析，可以很清晰地看出能量隨時(shí)間和頻率的變化，并且得出低頻磁場信號的優(yōu)勢頻率在0～40Hz之間，能量最強(qiáng)點(diǎn)位于信號的第5～6s之間。

(4)與其他分析方法相比，HHT分析方法具有完全的局部時(shí)頻特征，可以準(zhǔn)確地刻畫低頻磁場信號的動態(tài)變化特征，且可以很好地刻畫信號的突變點(diǎn)信息，這為低頻電磁信號的檢測、信號濾波、數(shù)據(jù)篩選等提供了有利條件。

［1］何學(xué)秋，劉明舉.含瓦斯煤巖破壞電磁動力學(xué)［M］.徐州:中國礦業(yè)大學(xué)出版社，1995.

［2］王恩元，何學(xué)秋.煤巖變形破裂電磁輻射的實(shí)驗(yàn)研究 ［J］.地球物理學(xué)報(bào)，2000，43(1):131－137.

［3］聶百勝，何學(xué)秋，王恩元，等.煤體剪切破壞過程電磁輻射與聲發(fā)射研究 ［J］.中國礦業(yè)大學(xué)學(xué)報(bào)，2002，31(2):609－611.

［4］朱郴韋.煤體破裂電磁輻射信號波形特征及降噪方法研究［D］.北京:中國礦業(yè)大學(xué) (北京)，2009.

［5］Cohen L.Time－frequency analysis［M］.New Jersey:Prentice Hall，1995.

［6］謝 中，程迎軍，徐清燕.電解氣泡析出時(shí)電位波動的頻譜分析［J］.中國有色金屬學(xué)報(bào)，2003，13(4):1011－1016.

［7］劉希靈，李夕兵，洪 亮，等.基于離散小波變換的巖石SHPB測試信號去噪［J］.爆炸與沖擊，2009，29(1):67－72.

［8］周子龍，李夕兵，龍八軍.巖石SHPB試驗(yàn)信號的小波包去噪［J］.巖石力學(xué)與工程學(xué)報(bào)，2005，24(S1):4780－4783.

［9］Huang N E，Shen Z，Long S R，et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non－station time series analysis［C］.Proceeding of the Royal Society of London，1998.

［10］Yue H Y，Guo H D.A SAR.Interferogram filter based on the empirical mode decomposition method［C］.Proceedings of Geosoience and Remote Sensing Symposium.IGARSS O1，2001.

［11］Semion Kizhner，Thomas P Flatley，et al.On the Hilbert－ Huang transform data processing system development［A］.2004 IEEE Aero pace Conference Procedings［C］.2004:1961－1979.

［12］李夕兵，凌同華，張義平.爆破震動信號分析理論與技術(shù)［M］.北京:科學(xué)出版社，2009.

［13］Norden E Huang，Zheng Shen，Steven R Long，et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non－stationary time series analysis［J］.The Royal Society，1998，454:903－995.

［14］武安緒，吳培稚，蘭從欣，等.Hilbert－Huang變換與地震信號的時(shí)頻分析［J］.中國地震，2005，21(2):207－216.

［15］李勝林，劉殿書，李祥龍，等.75mm分離式霍普金森壓桿試件長度效應(yīng)的實(shí)驗(yàn)研究［J］.中國礦業(yè)大學(xué)學(xué)報(bào)，2010，39(1):93－97.

［16］ZENCKER U，CLOS R.Limiting conditions for compression testing of flat specimens in the Hopkinson pressure bar［J］.Experimental Mechaincs，1998，39(4):343－348.

［17］WEINONG W，CHEN B S.Split Hopkinson(Kolsky)bar design，testing and applications［M］.Springer New York Dordrecht Heidelberg London，2011:7－17.

［18］李成武，解北京，楊 威.基于HHT法的煤沖擊破壞SHPB測試信號去噪 ［J］.煤炭學(xué)報(bào)，2012，37(11):1796－1801.

［19］Daubechies I.小波十講［M］.李建平，楊萬年，譯.北京:國防工業(yè)出版社，2004:38－40.

［20］付玉凱，李成武，段昌瑞，等.煤體失穩(wěn)破壞過程中的低頻磁場變化特征研究 ［J］.煤礦開采，2014，19(4):13－17，76.

［21］林 君，項(xiàng)葵葵，朱寶龍，等.MT信號現(xiàn)場處理的實(shí)現(xiàn)技術(shù)研究［J］.數(shù)據(jù)采集與處理，1997，12(1):52－55.

［22］李成武，解北京，楊 威.煤沖擊破壞過程中的近距離瞬變磁場變化特征研究 ［J］.巖石力學(xué)與工程學(xué)報(bào)，2012，31(5):973－981.

［23］程 磊，瞿偉廉.基于Hilbert－Huang變換理論的非平穩(wěn)數(shù)據(jù)處理［J］.建筑科學(xué)與工程學(xué)報(bào)，2007，24(1):26－30.