基于HHT法的煤沖擊破壞低頻電磁信號(hào)去噪

付玉凱，楊 威，李成武

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

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

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

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

本文采用HHT時(shí)頻分析方法對(duì)煤沖擊破壞的低頻電磁信號(hào)進(jìn)行經(jīng)驗(yàn)?zāi)Ｊ椒纸?(EMD)，并通過對(duì)各個(gè)IMF(某一頻率尺度上的模態(tài)信號(hào))分量進(jìn)行重構(gòu)，最后對(duì)重構(gòu)信號(hào)進(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è)作用:過濾疊加波和對(duì)稱化波形。EMD可以將信號(hào)從時(shí)間尺度上進(jìn)行IMF分量分解，但是需要滿足下面兩個(gè)條件:

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

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

1.1.2 EMD算法

對(duì)一個(gè)原始信號(hào)X(t)，首先找出X()t上全部的最大值和最小值，然后采用插值方法對(duì)曲線進(jìn)行極值擬合，從而繪制出曲線的包絡(luò)線Xmax()t。同理得出包絡(luò)線Xmin()t，2條包絡(luò)線包含了全部信號(hào)數(shù)據(jù)。對(duì)2條包絡(luò)線取其平均值得平均線m1()t，再用X()t減掉m1()t得到h1()t。如果信號(hào)不同，h1()t有可能產(chǎn)生一個(gè)IMF分量，也可能得到2個(gè)IMF分量。若分量不滿足限定條件，此時(shí)將h1()t當(dāng)做原信號(hào)，重復(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)算中，可以通過對(duì)信號(hào)反復(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 變換

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

1.2.2 HHT譜

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

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

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

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

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

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

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

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

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

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

2.2 煤樣制作

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

2.3 煤沖擊破壞低頻電磁信號(hào)去噪前分析

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

圖3 原始信號(hào)

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

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

圖4 原始信號(hào)的FFT譜

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

3 低頻磁場(chǎng)信號(hào)的HHT分析

3.1 磁場(chǎng)信號(hào)的EMD分解及重構(gòu)

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

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

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

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

圖6 原始信號(hào)EMD分解

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

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

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

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

圖9 (a)是對(duì)重構(gòu)信號(hào)進(jìn)行的HHT時(shí)頻譜，從圖中可以看出，顏色越亮表示其能量越高，反之越低;HHT時(shí)頻譜直觀地表現(xiàn)出了信號(hào)的聚集性，與原始信號(hào)的時(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)信號(hào)的邊際時(shí)間，表示每個(gè)時(shí)間在全局上的幅度之和，是表示時(shí)間點(diǎn)上信號(hào)能量強(qiáng)弱的一個(gè)指標(biāo)，從該圖中也可以看出，第5s至第6s時(shí)間域上的能量較強(qiáng)，很好地驗(yàn)證了HHT時(shí)頻譜分析的正確性。

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

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

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

4 結(jié)論

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

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

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

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

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

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

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

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

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

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

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

［8］周子龍，李夕兵，龍八軍.巖石SHPB試驗(yàn)信號(hào)的小波包去噪［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］李夕兵，凌同華，張義平.爆破震動(dòng)信號(hào)分析理論與技術(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變換與地震信號(hào)的時(shí)頻分析［J］.中國(guó)地震，2005，21(2):207－216.

［15］李勝林，劉殿書，李祥龍，等.75mm分離式霍普金森壓桿試件長(zhǎng)度效應(yīng)的實(shí)驗(yàn)研究［J］.中國(guó)礦業(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測(cè)試信號(hào)去噪 ［J］.煤炭學(xué)報(bào)，2012，37(11):1796－1801.

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

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

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

［22］李成武，解北京，楊 威.煤沖擊破壞過程中的近距離瞬變磁場(chǎng)變化特征研究 ［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.