基于改進EMD算法的木材聲發射信號特征研究

方賽銀， 邱榮祖， 李 明

(1.福建農林大學 交通與土木工程學院，福州 350002；2.西南林業大學 機械與制造工程學院，云南 650224)

材料在受外力或內力作用產生變形或斷裂時，以彈性波的形式釋放應變能的現象稱為聲發射(Acoustic Emission, AE)。AE信號能夠動態反映材料內部應力應變的產生和發展狀況，從而為材料損傷提供了一種主動的無損監測方法。Ohuchi等[1]通過拉伸試驗研究木材早晚材的AE特征，Ritschel等[2-3]利用AE事件數研究木材和層積材內部損傷的發展規律。孫建平等[4]研究木材在動態載荷下的AE事件的變化規律，郭曉磊等[5]利用AE特征研究木基復合材料內部損傷過程，丁小康等[6]根據AE信號特征研究木材干燥開裂過程。上述研究的主要依據是AE事件的相關統計參數，然而由于木材是多孔的各向異性材料，無論是AE信號的傳播還是衰減規律都與金屬等材料存在本質差異，僅靠AE信號的幅值來確定AE事件是不準確的。為此，徐鋒等[7-9]提出根據AE信號的波形特征研究木材和膠合板損傷過程。

受信號采集噪聲和衰減影響，信號降噪和波形析取成為木材AE信號分析的關鍵。木材AE信號是一種非線性、非平穩過程，目前主要采用小波技術進行AE信號時頻域分析，然而小波核函數的選擇及分解層數的確定都直接影響信號分析的結果。Huang于1998年提出了經驗模態分解(Empirical Mode Decomposition, EMD)算法，為非線性和非平穩過程分析提供了完全自適應的分析方法，特別是分解后的固有模態函數(Intrinsic Mode Functions, IMF)具有明確的物理意義。包絡擬合是EMD算法的核心，包絡的不準確將導致信號分解的不完全。經典的EMD算法采用三次樣條插值(Cubic Spline Interpolation, CSI)生產包絡線，但三次樣條插值算法容易產生過沖現象，進而影響EMD分解的準確性。為此，許多改進的包絡擬合方法被用于EMD分解算法，其一是采用其它插值算法進行包絡擬合，主要包括B樣條插值[10]、分段冪函數插值[11]、分段Hermite插值[12]等。其二是通過增加插值點的方法提高包絡精度[13]，但此類算法需要估算新的插值點，從而增加算法的復雜性。其三是通過不同的方法產生多條包絡線[14]，然后在其中選擇最優的包絡擬合。

現有改進的包絡擬合方法在一定程度上改善了包絡的精度，但是往往付出了復雜性顯著增加的代價，而且從本質上看，都是基于插值點上的函數值信息，沒有引入新的關于節點信息。為此，本文根據EMD包絡生成原理，增加插值點上的一階導數信息，進而提出改進的分段三次Hermite插值(Improved Piecewise Cubic Hermite Interpolation, IPCHI)算法用于包絡擬合。IPCHI算法能夠在不增加算法復雜度的同時提高包絡精度，同時有效避免三次樣條插值過程中誤差向數據內部傳播而“污染”整個數據序列的問題。隨后，基于IPCHI包絡擬合方法提出改進的EMD算法。最后，針對木材受力彎曲過程的AE信號，提出基于小波降噪和EMD分解聯合的AE信號分析方法，進而提出依據瞬時頻率進行AE事件統計的方法。

1 經典EMD算法的基本原理

EMD算法的目的在于將信號分解為一系列表征信號特征時間尺度的IMF，再通過Hilbert變換獲得信號的時頻關系。理論上EMD算法可以處理任意非線性且非平穩信號，實際受IMF基本要求的限制，可進行EMD分解的信號通常需要滿足以下三個假設：①信號至少有兩個極值，即一個極大值和一個極小值；②信號特征時間尺度是由極值間的時間間隔確定的；③若數據沒有極值點，可以將拐點的一階或高階導數值作為極值點。

此外，IMF還必須滿足兩個條件：①在整個數據長度上，極值點和過零點的數目必須相等或最多相差一個；②在任意數據點，局部最大值的包絡和局部最小值的包絡的平均必須為零。

EMD算法的基本步驟如下：

步驟1 確定給定信號x(t)所有的局部極大值點和極小值點。然后利用插值方法擬合出由所有極大值點構成的上包絡線u0(t)和所有極小值點構成的下包絡線v0(t)。

步驟2 根據下式計算上下包絡線的均值：

(1)

將信號x(t)與m0(t)的差記為h0(t)=x(t)-m0(t)。

步驟3 判斷h0(t)是否滿足IMF的兩個條件。若滿足，則h0(t)為篩選出來的第一個IMF分量；如果不滿足，則將h0(t)作為原始數據重復進行上述兩個步驟，直到h0(t)滿足IMF的兩個條件。第一次篩選得到的IMF分量記為c1(t)。

步驟4 將c1(t)從原始信號x(t)中分離出來得到剩余信號r1(t)=x(t)-c1(t)。然后將r1(t)作為新的待分析信號，重復步驟1～步驟3，通過篩選得到第二個IMF分量c2(t)，再計算出余項r2(t)=r1(t)-c2(t)。重復上述步驟，直至得到的余項rn(t)是一個單調信號或rn(t)的值小于預先給定的閾值，分解結束。

經過上述4個步驟的計算，最終得到n個IMF分量，c1(t),…,cn(t)，余項為rn(t)，這樣，原始信號x(t)可以表示為

(2)

由于篩選過程會平滑信號的幅值，從而“抹殺”信號的物理特性，所以實際應用時可以如下定義的閾值SD作為判斷篩選結果是否為IMF分量的依據

(3)

其中T為信號序列的長度。SD也稱為篩分閾值，一般在0.2～0.3之間取值。

2 改進的EMD算法

EMD算法主要通過插值擬合的方法生成信號的上下包絡線。根據插值擬合理論，受龍格現象的影響，主要可以通過兩個途徑提高插值精度:①采用分段插值方法，即對插值區間進行細分，在每個小區間上采用低階插值函數；②提供更多關于插值節點的有效信息。由于EMD算法的上下包絡線是通過對極大和極小值點插值擬合生成的，所以插值區間取決于相鄰兩個極大或極小值點的距離，人為細分插值區間存在較大困難。因此，充分挖掘更多關于插值節點的信息成為改善擬合精度的有效途徑。

根據EMD算法基本原理可知，包絡線的所有插值節點都是信號函數的局部極大和極小值點，既然是函數的極值點，那么就可以合理地認為在插值節點上函數的一階導數值為零，從而將關于插值節點的可用信息量提高一倍。由于三次樣條插值是在整個信號區間上估算插值節點上的二階導數值，所有在求解的過程中，各個點上計算值相互影響，必然會導致擬合誤差的累積和傳播。與三次樣條插值不同，三階Hermite插值僅涉及兩個插值點，所以，在增加插值點一階導數信息的基礎上，采用分段三階Hermite插值方法，不僅可以保證擬合的光滑性及精度，同時可以有效避免擬合誤差的傳播。

2.1 改進PCHI擬合算法

對于n個插值節點x=(x1,x2,…,xn)的CSI擬合問題，主要是利用n個插值節點上的函數值f(x)=(f(x1),f(x2),…,f(xn))估算所有內點上的二階導數信息，具體通過求解一個關于插值節點二階導數值的三對角線性方程組來實現，而且為了保證線性方程組求解的唯一性，還需要補充兩個邊界條件，顯然在求解過程中，邊界點上的誤差將不斷向信號內部傳播，進而影響擬合精度。

PCHI擬合算法的基本原理是在每個小區間[xj-1,xj]，j=2,…,n上，依據兩個插值節點xj-1,xj上的函數值f(xj-1),f(xj)和一階導數值f′(xj-1),f′(xj)，采用三次多項式函數進行插值擬合，每個區間上的四個插值基函數如下所示

(4)

式中:x∈[xj-1,xj],hj=xj-xj-1,j=2,…,n。

因此，每個區間上的PCHI函數的表示

(5)

因為EMD算法中的信號包絡線的插值節點都是局部極大和極小值點，根據極值點的定義，可以認為插值節點上信號的一階導數值等于零，這樣就增加了以下n個插值條件

f′(xi)=0,i=1,…,n

(6)

將上述插值條件代入式(5)，可以得到IPCHI擬合函數為

Hj(x)=f(xj-1)l0(x)+f(xj)l1(x)

(7)

根據上述推導，在每個小區間上，IPCHI算法僅需要構造兩個插值基函數，從而有效降低了算法的復雜程度。此外，由于所有內點上的一階導數值都相等，所以IPCHI至少能夠保證與CSI擬合同樣的光滑性，而且每個區間上的插值函數僅與該區間的兩個端點信息有關，各區間上的插值函數保持相對獨立，從而避免了CSI中擬合誤差向信號內部傳播的問題。

另外，根據插值理論，PCHI在整個信號區間上的擬合誤差是

(8)

同樣情況下的CSI誤差是

(9)

通過上述兩式的比較發現，CSI的擬合誤差上限是PCHI的5倍，這也進一步說明了PCHI能夠有效提高擬合精度。

2.2 數值試驗

為了說明IPCHI的擬合效果，首先產生長度為100的隨機信號，然后分別采用CSI和IPCHI生成信號的上下包絡線,如圖1所示。為了確保上下包絡線能夠覆蓋全部信號，本文將信號兩側的端點視為極值點，具體處理過程是：若端點處信號為正值，則視為極大值點，同時對該點再賦以零值作為極小值點；反之，若端點處信號為負值，則視為極小值點，同時對該點再賦以零值作為極大值點。

圖1 CSI和IPCHI包絡擬合曲線

從圖1可知,CSI擬合出現了明顯的過沖現象，譬如在圖1中的A和B處，CSI擬合后出現了新的極大和極小值點，類似的情況還在其它地方多次出現。

為了進一步說明包絡的不準確對EMD算法的影響，本文分別采用CSI和IPCHI方法產生信號上下包絡線，并進行一次信號分解，IMF分量的判斷標準設定為SD=0.2，EMD分解后得到的第一個IMF分量,如圖2所示。

從圖2可知，兩種插值算法下得到的分量都能滿足IMF的兩個基本條件，而且從幅值上可以明顯看出EMD算法對信號的平滑作用。然而，通過比較發現，基于CSI擬合得到的IMF分量已經不能反映原始信號特征，特別是信號的后半部分已完全喪失了原始信號特征，顯然這種不完全分解必然會影響其后的信號分析。與之相反，基于IPCHI方法得到的IMF分量能夠很好地保留原始信號特征。

圖2 基于CSI和IPCHI的EMD分解

3 基于改進EMD算法的木材AE信號分析

常規的參數分析法主要依據AE事件的屬性，而由于缺乏有效的AE信號時頻特性，通常只能根據AE信號幅值判斷是否發生AE事件。然而，木材作為一種多孔的各向異性材料，受力過程釋放的AE信號不僅類型豐富，而且存在更為明顯的衰減特性，甚至會“淹沒”在測量噪聲中，實際使用時AE事件的確定和統計非常困難。為此，本文將采用EMD算法研究木材AE信號特征，首先利用小波技術對采集的原始AE信號進行濾波和波形重構，然后采用本文提出的改進EMD算法對對重構后的AE信號進行IMF分解，在信號相關性分析的基礎上，確定能夠反映木材AE信號特征的IMF主分量。最后通過對IMF主分量進行Hilbert變換，獲得木材AE信號的特征頻率，并進一步根據瞬時頻率確定AE事件，從而賦予AE事件更為明確的物理意義。

3.1 木材AE信號采集平臺

為了獲得可靠的原始AE信號，本文采用專門的AE傳感器和高速采集設備自行搭建木材AE信號采集系統，同時設計相應的軟件平臺。

AE信號采集系統主要由AE傳感器、前置放大器、高速數據采集設備等組成,如圖3所示。其中傳感器選用聲華SR150N單端諧振AE傳感器，采集頻率范圍為22～220 kHz，且配備了40 dB的PAI前置放大器。數據采集設備選用NI的8通道的高速數據采集設備USB-6366，該設備最大采樣頻率可達2 MHz。現有研究表明木材AE信號特征頻率通常在150 kHz以內，所以根據香農采樣定理并保留一定的裕度，本論文在試驗中設定的采樣頻率為500 kHz，采集的多通道AE原始信號在軟件平臺中分離后以數據文件形式保存。

3.2 試驗設計

圖3 木材AE信號采集系統

本論文主要通過三點彎曲力學試驗獲取木材AE信號。首先制作尺寸為400 mm×40 mm×13 mm的云南松(Pinus yunnanensis)試件，并確保試件尺寸均勻且無干燥缺陷。然后依據GB/T17657標準在力學試驗機上進行三點彎曲力學試驗,如圖4所示。設定試件跨度為300 mm，載荷從5 N開始，以2 mm/min的進給速率進行加載，直至試件斷裂為止。在試件表面施力點的兩側分別布置一個AE傳感器。

圖4 木材三點彎曲試驗

3.3 木材AE信號分析

由于試驗過程設定的采樣頻率較高(500 kHz)，所以完整的試驗數據非常龐大，而本文主要研究EMD算法在木材AE信號分析中的應用，所以這里僅截取了試驗后期一個長度為2 ms的數據片段,如圖5所示。該片段共包含10 000個數據，所有的AE信號分析均在MATLAB中完成。

圖5 AE原始信號和IMF波形

文獻[15]研究表明，對AE信號進行小波預處理后，再進行EMD分解能夠取得更好的效果。為此本文首先對原始木材AE信號進行小波降噪和波形重構，在小波分析時選擇正交db小波作為核函數，并設定小波分解層數為8，由于小波分析不是本論文的重點，所以這里不再給出小波分析的詳細過程。

為了說明小波預處理的作用，本文分別對小波處理前后兩組AE信號進行EMD分解處理，經過本文提出的改進EMD算法自適應分解，分別獲得13個IMF分量，為簡單起見，這里僅列出了小波處理后AE信號EMD分解后的所有IMF分量及余項,如圖6所示。

因為EMD分解具有非常明確的物理意義，IMF分量反映了信號在每個時刻的固有特性，因此，本文根據IMF分量與原信號的相關性來確定主分量，表1和表2分別為小波處理前后兩組AE信號經EMD分解后，前4個IMF分量與對應原信號的相關系數。

表1小波處理前IMF分量與原始信號的相關系數

Tab.1CorrelationcoefficientsbetweenIMFsandoriginalAEsignals

IMFIMF1IMF2IMF3IMF4相關系數0.553 80.590 20.472 60.106 9

表2小波處理后IMF分量與原始信號的相關系數

Tab.2CorrelationcoefficientsbetweenIMFsandAEsignalstreatedbywavelet

IMFIMF1IMF2IMF3IMF4相關系數0.507 20.806 60.197 10.067 5

依據相關性比較，兩組AE信號都可以將第2個IMF分量作為主分量,如圖5所示。表1中數據顯示所有IMF分量與原信號的相關性都不超過0.6，甚至前兩個IMF分量相差無幾，這主要是因為受噪聲影響，EMD算法只能對AE信號進行很不完全的分解，介于這種原因，以相關性作為選擇IMF主分量的依據也不再充分。與之相反，表2中數據說明經過小波預處理后，再進行EMD分解的效果明顯提高，而且可以根據相關性條件容易選定第2個分量作為IMF主分量。

為了進一步研究AE信號的頻率特性，本文分別對原始AE信號及其EMD分解后的IMF主分量、小波處理后AE信號的IMF主分量進行了頻域分析，相應的頻譜圖,如圖7所示。

圖7顯示未經小波預處理的信號頻率分散，很難確定木材AE信號的特征頻率范圍，然而，AE信號經小波預處理后，IMF主分量的頻域特性明顯集中。根據IMF主分量的頻率分布可以確定木材AE信號的特征頻率中心位于46.5 kHz，為了通過瞬時頻率判斷AE事件的發生，本文設定木材AE信號的特征頻率范圍為45～47 kHz，即瞬時頻率在此范圍內，均視為發生了AE事件。雖然特征頻率范圍的設定將直接影響AE事件的數量，但是只要在同一個試驗中保持同樣的標準，對后繼分析沒有本質的影響，本文以頻譜圖中最大幅值的80%作為設定依據。

圖6 EMD分解后的IMF分量及余項

在確定木材AE信號的特征頻率范圍后，就能夠依據瞬時頻率統計AE事件，圖8是小波預處理后AE信號的IMF主分量的瞬時頻率變化曲線，依據以上設定的特征頻率范圍，經統計在0.02 s的信號片段中，共發生了373次AE事件，顯然，依據瞬時頻率判斷AE事件具有明確的物理意義。

(a) 頻譜圖(原始信號)

(b) 頻譜圖(原始AE信號IMF主分量)

(c) 頻譜圖(小波處理后AE信號的IMF主分量)

圖8 IMF主分量的瞬時頻率曲線

4 結 論

AE技術為木材損傷提供了有效的主動無損檢測方法，受復雜的物理結構影響，木材AE信號呈現較強的非線性和非平穩性特征，為此，本文結合小波預處理和EMD分解相結合的方法，在降噪的同時提取木材AE信號的頻率特征，并且從瞬時頻率的角度判定AE事件的發生，使得AE事件的定義更具物理意義。

本文根據生成信號包絡線的插值節點的極值性質，將插值節點處的一階導數值為零增加為新的插值條件，從而形成改進的分段三階Hermite插值擬合算法，不僅有效解決了三次樣條插值容易出現過沖的現象，還有效避免了擬合誤差的傳播，從而提高包絡擬合的精度，進而在IPCHI包絡擬合的基礎上提出改進的EMD算法。

其后以云南松試件為對象，通過三點彎曲試驗獲取木材損傷過程的AE信號，并截取0.02 s的信號片段進行分析，通過比較分析，AE信號經過小波預處理后，再進行EMD分解可以得到較為理想的IMF主分量，并且根據IMF主分量的頻率分布容易確定木材AE信號的特征頻率范圍，最后利用瞬時頻率判定并統計AE事件，從而為木材AE信號特征研究提供了一種有效的分析方法。

未來可以將本文提出的方法進一步推廣到整個木材損傷過程，研究木材損傷程度與AE事件累積量之間的本構關系，從而建立基于AE特征的木材損傷預測本構模型。