巖石聲發射信號能量分布特征的EMD分析

凌同華，張 勝，易志強，李品鈺

(長沙理工大學 土木與建筑學院，長沙 410004)

聲發射技術的研究對象已由金屬為主擴展到巖石、混凝土、復合材料等固體材料。作為一種無損檢測手段，聲發射技術已廣泛地應用于地下采礦、石油化工、材料試驗、民用工程、航天和航空等眾多領域。

對聲發射信號的分析處理是聲發射技術能夠得到更好地發展和推廣的重要前提，巖石聲發射信號具有瞬態性和多樣性的特點，屬于典型的非平穩信號。用傅里葉頻譜、小波(包)分析處理聲發射信號已得到廣泛應用。劉希靈［1］編制快速傅里葉變換頻譜分析程序對沖擊荷載作用下巖石聲發射信號進行頻譜分析，得出沖擊荷載作用下巖石聲發射信號采集的時間短，信號的能量、信號強度、絕對能量在某一個波擊數處有很大的值以及低頻較顯著、很集中，且低頻和其他諧波成分的振幅都很大。王余剛［2］提取了混凝土材料不同破壞階段的全波形聲發射信號并分析了其頻譜特性。康玉梅［3］基于小波變換，研究了聲發射信號時頻能量特征提取的原理和方法，利用時頻能量分析技術來實現聲發射信號的時延估計。趙奎［4］采用小波包頻帶分解方法，對巖石聲發射Kaiser點信號的能量分布特征進行了研究，分析了砂巖聲發射信號不同頻帶的能量分布規律。凌同華［5］利用小波包分析技術對沖擊荷載作用下巖石聲發射信號的能量分布特征進行了研究，重點討論了沖擊荷載作用下不同巖石對聲發射信號頻帶能量分布的影響。聶鵬［6］基于聲發射信號特點和小波包分解理論對不平穩信號特征提取的優勢，提出了一種利用聲發射信號的能量變化來監測刀具磨損狀態的方法。趙元喜［7］根據滾動軸承聲發射信號在各頻段的能量分布與軸承的故障類型相關性，利用諧波小波包將不同故障滾動軸承的聲發射信號分解到不同頻段，進而將各頻段的能量組成特征向量輸入BP神經網絡，通過神經網絡判別滾動軸承的故障類型。

HHT(Hilbert-Huang Transform)法是近年提出的適合處理穩態信號和非平穩信號的有效方法，它由EMD方法和Hilbert變換兩部分組成，其核心是EMD方法。相對小波(包)等分析方法而言，EMD(Empirical Mode Decomposition)方法是一種更具適應性的時頻局部化分析方法，它沒有固定的基函數，是自適應的［8－10］。本文針對沖擊荷載作用下巖石聲發射信號的特點，用EMD方法進行分析，探討了沖擊荷載作用下巖石聲發射信號的能量分布特征。

1 EMD方法的原理

EMD方法依據信號本身的局部特征信息進行自適應地分解，得到一系列具有不同特征時間尺度的IMF(Intrinsic Mode Function)分量。IMF分量必須滿足下面兩個條件:① 在整個時間序列內，極值點的個數和過零點的個數必須相等或相差最多不能超過一個;②在任意時刻，由局部極大值點、局部極小值點分別形成上、下包絡線的平均值為零。

采用EMD方法通過以下步驟對巖石聲發射信號x(t)進行分解［11－12］:

(1)確定信號所有的局部極值點，然后用三次樣條曲線將所有的局部極大值點連接起來形成上包絡線;

(2)用三次樣條曲線將所有的局部極小值點連接起來形成下包絡線，上、下包絡線應該包絡所有的數據點;

(3)上、下包絡線的平均值記為m1(t)，求出:

在理想狀況下，如果h1(t)滿足IMF分量的兩個基本條件，那么h1(t)就是x(t)的第1個IMF分量;

(4)如果h1(t)不滿足IMF分量的兩個基本條件，把h1(t)作為原始數據，重復步驟(1)～步驟(3)，得到上、下包絡線的平均值m11(t)，再判斷h11(t)=h1(t)－m11(t)是否滿足IMF分量的兩個基本條件，若不滿足，則重復循環k次，得到 h1k(t)=h1(k－1)(t)－m1k(t)，使得h1k(t)滿足IMF分量的兩個基本條件。記c1(t)=h1k(t)，則c1(t)為信號x(t)的第1個IMF分量;

(5)從原始信號x(t)中分解出第1個IMF分量c1(t)后，從x(t)中減去c1(t)，得到剩余值序列r1(t):

將r1(t)作為原始數據重復步驟(1)～步驟(4)，得到x(t)的第2個IMF分量c2(t)，以此類推可得第3、第4個直至第n個IMF分量，即為c3(t)，c4(t)，…，cn(t)，最后剩下原始信號的余項rn(t)。

這樣，原始信號x(t)可分解為若干IMF分量和一個余項rn(t)的和，即:

步驟(5)中的停止條件被稱為分解過程的停止準則，它可以是如下兩種條件之一:① 當最后一個IMF分量cn(t)或剩余分量rn(t)比預期值小時便停止;②當剩余分量rn(t)變成單調函數，從中不能再分解出IMF分量為止。

2 巖石聲發射信號能量的表征

巖石聲發射信號經EMD分解后的各IMF分量分別代表了一組特征尺度下的平穩信號，各平穩信號能量的變化就表征了巖石聲發射信號的特征情況。

基于EMD方法的巖石聲發射信號能量特征提取步驟如下［13－15］:

(1)對巖石聲發射信號進行EMD分解，求得各IMF分量ci(t);

(2)求各IMF分量的總能量Ei:

(3)以能量為元素構造一個特征向量T:

當能量較大時，Ei(i=1，2，…，n)通常是一個較大的數值，在數據分析上會帶來一些不方便的地方。由此，可以對特征向量T進行改進，即對向量進行歸一化處理。令:

向量T'即為歸一化后的向量。

3 沖擊荷載作用下巖石聲發射信號的采集

本次試驗采用PCI－2聲發射儀，其中有中心響應頻率為250(500)kHz的PICO型諧振式窄頻帶傳感器、2/4/6型前置放大器及AE–win聲發射采集軟件。試驗選用長徑比為1:2的三種不同的巖石試件(花崗巖、石灰巖、矽卡巖)，在桿徑為50 mm的SHPB(Split Hopkinson Pressure Bar)試驗系統上分別對它們進行沖擊荷載作用下聲發射試驗。表1為試件的物理力學實驗指標，其中花崗巖有輕微破碎，節理不發育，抗壓強度高，表面硬度大，化學穩定性好，耐久性強;石灰巖有輕微破碎，節理裂隙不發育，具有良好的加工性、磨光性和很好的膠結性能;矽卡巖呈深灰色，節理不發育，塊狀構造，相對密度較大，且三種巖石試件完整性較好。每個試件在加工過程中進行標記，如“y1”對應花崗巖試件，“y2”對應石灰巖試件，“y3”對應矽卡巖試件。

表1 試件的物理力學實驗指標Tab.1 Specimen of physical and mechanical test indexes

將聲發射探頭(傳感器)的接觸面涂抹少量的耦合劑(如黃油等)，其目的是獲得高質量的信噪比及減少背景干擾信號，然后用膠帶綁在試件表面居中位置，最后將巖石試件夾在SHPB試驗裝置入射桿和透射桿之間。試驗中花崗巖、矽卡巖的采樣頻率為4×107Hz，而石灰巖采樣頻率為5×106Hz，為便于分析，被分析信號的采樣頻率必須一致，對石灰巖信號進行重新采樣(重新采樣并不改變信號的本質特征)，重新采樣的頻率為4×107Hz。試驗時采集了多組信號，現取每次試驗聲發射起始點的信號進行分析。圖1為沖擊荷載作用下巖石聲發射信號的幅值時程曲線。

圖1 沖擊荷載作用下巖石聲發射信號的幅值時程曲線Fig.1 The amplitude versus time curves of rock acoustic emission signals under impact loading

4 聲發射信號的EMD分解及能量分布特征分析

4.1 預處理

實測沖擊荷載作用下巖石聲發射信號除具有反映有關巖石本身的特征信息外，還包含大量的背景噪聲。因此，在提取特征向量前有必要對信號進行預處理，以突出特征信息，提高信號的信噪比。

4.2 EMD分解層

EMD方法能對非平穩、非線性信號進行平穩化、線性化處理，并在分解的過程中保留原始信號的固有特性。因此，在Matlab語言平臺上采用EMD方法對采集的聲發射信號進行分解得到n個IMF分量c1(t)，c2(t)，…，cn(t)，每一個IMF分量都包含了不同的特征尺度信息，這樣通過EMD分解，原始信號的特征完全可以由這n個IMF分量c1(t)，c2(t)，…，cn(t)來表征。考慮到原始信號的特征信號主要集中在高頻帶，因此可以只選取前8個IMF分量作為進一步的研究對象。巖石聲發射信號y1，y2，y3的EMD分解結果見圖2～圖4。

圖2 巖石聲發射信號y1的EMD分解結果Fig.2 The EMD results of rock acoustic emission signal y1

圖3 巖石聲發射信號y2的EMD分解結果Fig.3 The EMD results of rock acoustic emission signal y2

圖4 巖石聲發射信號y3的EMD分解結果Fig.4 The EMD results of rock acoustic emission signal y3

從圖中可以看出，第1個IMF分量是原始信號中分解出的時間尺度最短、頻率最高的分量，代表信號中的高頻成分，且其振幅大，表明其所占能量大，然后依次分解出余下各IMF分量。隨著分解的進行，所得各IMF分量時間尺度越來越長、頻率越來越低，直到分解出時間尺度最長、頻率最低的最后一個IMF分量。

EMD方法把沖擊荷載作用下巖石聲發射信號分解成8個IMF分量，不同的IMF分量包含不同的特征時間尺度，可以使信號的特征信息在不同的分辨率下顯示出來，表明EMD分解中分辨率是自適應的，相對小波的多分辨率分析，EMD方法在分解信號時更加簡單。EMD分解中沒有固定的基函數，每個IMF分量的提取都是由信號本性所決定的，分解出的各IMF分量與原始信號相似。在數字信號處理中常用相關系數來反映兩個信號的相似程度［16］，各IMF分量與原始信號相關系數計算結果見表2。

表2 IMF分量與原始信號的相關系數Tab.2 Correlation coefficients between the IMF components and the original signals

4.3 巖石聲發射信號的能量分布特征

沖擊荷載作用下巖石聲發射信號經EMD方法分解后得到8個IMF分量，按式(4)～式(7)分別提取每個IMF分量的特征向量T'。計算結果見表3、圖5。

表3 巖石聲發射信號的能量分布Tab.3 Energy distribution of rock acoustic emission signals

(1)從表3和圖5可以看出，3種巖石聲發射信號經EMD分解出的前4個IMF分量的能量占該信號總能量的比例分別為 93.96%，85.43%，85.50%，表明沖擊荷載作用下巖石聲發射信號的能量在頻域上分布比較廣泛，但絕大部分能量集中在前4個IMF分量內，且分布不均勻。

(2)為了顯示巖石聲發射信號各IMF分量相關系數與能量分布的關系，以各IMF分量為橫坐標，將表2和表3中的數值繪制成曲線，如圖6所示。從圖6中可以看出各IMF分量相關系數與能量分布存在很大相關性，兩者的相關系數分別為0.94，0.84，0.96，表明各 IMF分量相關系數與能量分布的相關程度高，反映了EMD方法是依據信號本身的固有特性進行自適應地分解。

圖5 巖石聲發射信號的能量分布Fig.5 The energy distribution of rock acoustic emission signals

圖6 各IMF分量的相關系數與能量分布Fig.6 Correlation coefficients and energy distribution of IMF components

圖7 巖石聲發射信號的功率譜Fig.7 The power spectrum of rock acoustic emission signals

圖8 巖石聲發射信號y1的IMF分量(c1～c3)的功率譜Fig.8 The power spectrum of IMF components(c1～ c3)of rock acoustic emission signal y1

(3)沖擊荷載作用下巖石聲發射信號的頻譜豐富，且絕大部分分布在500 kHz以下。從圖7中可知3種巖石聲發射信號在0～500 kHz的能量占該信號總能量的比例分別為 99.56%，97.42%，85.32%，并可分成多個子震頻帶，如圖8中分量c1，c2及c3各自的主頻率所示，這一點與圖7中聲發射信號y1的頻譜基本一致。與此同時，圖7中花崗巖聲發射信號的優勢頻率較石灰巖、矽卡巖更集中，這一特征同參考文獻［5］是一致的，但EMD方法沒有固定的基函數，具有自適應性，使信號分析更加靈活多變;而小波包分析中不同的小波基具有不同的特性，對信號的分析效果也不相同，對同一個信號采用不同的小波基得到的結果基本上也沒有可比性。

(4)從表1、圖6和圖7還可以看出，隨著巖石單軸壓縮強度增加，沖擊荷載作用下巖石聲發射信號分解出的各IMF分量相關系數與能量分布的相關程度越高;隨著巖石的密度、縱波波速、彈性模量的降低，沖擊荷載作用下巖石聲發射信號的優勢頻率越來越集中，且其優勢頻率有往低頻發展的趨勢。

5 結論

(1)EMD方法把沖擊荷載作用下巖石聲發射信號分解為有限個不同特征時間尺度的IMF分量，其信號能量主要分布在前4個IMF分量內，且分布不均勻;

(2)EMD方法分解出的IMF分量大都有清晰的物理意義，且各IMF分量的頻譜與原始信號的頻譜基本一致，表明EMD方法分解出的各IMF分量集中了原始信號中最顯著的特征信息;

(3)隨著巖石的密度、縱波波速、彈性模量的降低，沖擊荷載作用下巖石聲發射信號的優勢頻率越來越集中，且其優勢頻率有往低頻發展的趨勢。相比小波法分析［5］，EMD分析不需預設基函數，具有自適應和高效的特點。

［1］劉希靈.基于激光三維探測的空區穩定性分析及安全預警的研究［D］.長沙:中南大學，2008.

［2］王余剛，駱 英，柳祖亭.全波形聲發射技術用于混凝土材料損傷監測研究［J］.巖石力學與工程學報，2005，24(5):803－807.

［3］康玉梅，朱萬成，白 泉，等.基于小波變換時頻能量分析技術的巖石聲發射信號時延估計［J］.巖石力學與工程學報，2010，29(5):1010 －1016.

［4］趙 奎，王更峰，王曉軍，等.巖石聲發射Kaiser點信號頻帶能量分布和分形特征研究［J］.巖土力學，2008，29(11):3082－3088.

［5］凌同華，廖艷程，張 勝.沖擊荷載下巖石聲發射信號能量特征的小波包分析［J］.振動與沖擊，2010，29(10):127－130.

［6］聶 鵬，王東磊，徐 濤，等.頻帶能量特征法在聲發射刀具磨損監測系統中的應用［J］.工具技術，2009，43(2):24－27.

［7］趙元喜，胥永剛，高立新，等.基于諧波小波包和BP神經網絡的滾動軸承聲發射故障模式識別技術［J］.振動與沖擊，2010，29(10):162 －165.

［8］ Huang N E，Shen Z，Long S R，et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis［J］.Proc.R.Soc.Lond.A，1998，454:903 －995.

［9］Huang N E，Shen Z，Long S R.A new view of nonlinear water waves:the hilbert spectrum［J］.Annual Review of Fluid Mechanics，1999，31:417 －457.

［10］ Ma S，Zhang R.Empirical mode decomposition of the 1994 northridge earthquake and its interpretation for seismic source mechanism［A］.The 10th International Conference on Soil Dynamics and Earthquake Engineering.USA，Philadelphia:［s.n.］，2001:7 －10.

［11］于德介，程軍圣，楊 宇.機械故障診斷的Hilbert－Huang變換方法［M］.北京:科學出版社，2006.

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

［13］楊 宇，于德介，程軍圣.基于EMD與神經網絡的滾動軸承故障診斷方法［J］.振動與沖擊，2005，24(1):85－88.

［14］胡昌華，李國華，周 濤.基于MATLAB7.X的系統分析與設計－小波分析(第三版)［M］.西安:西安電子科技大學出版社，2008.

［15］凌同華，李夕兵.單段爆破振動信號頻帶能量分布特征的小波包分析［J］.振動與沖擊，2007，26(5):41－43.

［16］樓順天，劉小東，李博涵.基于MATLAB7.X的系統分析與設計－信號處理(第二版)［M］.西安:西安電子科技大學出版社，2005.