基于匹配延拓的ITD邊界效應抑制方法研究

胥永剛, 陸 明, 付 勝, 張建宇

(北京工業大學 機電學院 北京市先進制造技術重點實驗室，北京 100124)

大型機電設備的軸承齒輪等零部件發生故障時，振動信號常呈現出非線性非平穩特征，且淹沒于強烈的背景噪聲中。如何處理復雜的非平穩振動信號，一直是故障診斷領域的研究熱點之一[1]。

Frei等[2]提出了一種新的非線性非平穩信號處理方法——固有時間尺度分解法(Intrinsic time-scale decomposition, ITD),并將該方法應用于生物醫學信號處理中, 取得了理想的效果。ITD 分解法可以自適應地將任意復雜信號分解為若干具有實際物理意義的固有旋轉分量(Proper Rotation Component, PRC) 和一個單調趨勢項，然后再進一步求出所有PRC 的瞬時幅值和瞬時頻率，即可得到原始信號完整的時頻分布。該方法一經提出，便迅速在通信[3]、空間信息[4]及故障診斷[5-7]領域得到應用，均取得了不錯的效果。

在ITD方法的應用過程中，分解所得的固有旋轉分量常在邊界處出現畸變，使得分解結果嚴重失真，稱之為ITD方法的邊界效應。本文對比分析了五種數據延拓方法(自適應波形匹配法、基于AR模型的延拓方法、鏡像延拓方法、多項式延拓方法和反對稱周期延拓法)在ITD邊界效應抑制中的應用，并以趨勢分量絕對值累加和為評價指標，判定基于自適應波形匹配方法的數據延拓更適合抑制ITD邊界效應，分解誤差最小。最后，將其運用到實際信號的分析，驗證了該方法的有效性。

1 ITD方法及其邊界效應

1.1 ITD方法的基本原理

設待分解信號Xt為一實值離散信號，{τk,k=1,2,…}代表Xt中所有局部極值點對應的時刻，方便起見，定義τ0=0。定義L為Xt的基線提取算子，一次ITD分解是將信號Xt分解為基線分量Lt和固有旋轉分量Ht，即[2]：

Xt=LXt+(1-L)Xt=Lt+Ht

(1)

式中，Lt=LXt為基線分量，代表信號中的局部相對低頻成分；Ht=(1-L)Xt為固有旋轉分量，表示信號中的局部相對高頻成分。

簡便起見，令Xk和Lk分別表示X(tk)和L(tk)，設Lt)和Ht在[0,τk]上有定義，Xt在[0,τk+2]上有定義，則在連續極值點間隔[τk,τk+1]上可定義該區間內Xt的分段線性基線提取因子L：

(2)

(3)

式中，α用于控制提取固有轉動分量幅度的線性縮放，α∈[0,1]，通常取α=0.5[2]。則，固有旋轉因子Ht為：

HXt=(1-L)Xt=Ht=Xt-Lt

(4)

這種基線構造方法形成的基線分量Lt保留了信號在各個極值點之間的單調性，提取了各個極值點之間疊加的局部高頻分量，即固有于信號某種尺度上的分量——固有旋轉分量。將基線分量當做新的待分解信號重復以上的分解過程，可以得到一系列的固有旋轉分量，直到獲得一個單調趨勢信號。這就將原始信號Xt分解成若干個從高到低不同頻率段的固有旋轉分量之和與一個單調趨勢分量。整個過程可表示為：

(5)

式中，HLkXt是第(k+1)層固有旋轉分量，LpXt是單調的趨勢或提取的最低頻率基線。

1.2 ITD方法的邊界效應

設有一仿真信號：

x(t)=(1+0.5cos(2π×5t))cos(2π×100t+1.5sin(2π×7t))+0.5sin(2π×30t)

(6)

該信號包括兩個主要分量：一是調幅-調頻分量，二是純正弦分量。該信號的時域波形如圖1所示。

圖1 仿真信號的時域波形

該信號的ITD分解結果如圖2所示。由圖可知，原始信號x(t)被分解為2個固有旋轉分量(記為PRC1～PRC2)和一個趨勢分量(記為L)，其中PRC1和PRC2能量較大，分別對應原始信號中的調幅-調頻分量和正弦分量，L能量很小，為ITD分解過程中最終剩余的趨勢分量。可見，ITD方法能夠將復雜的非平穩信號分解為若干個簡單的固有旋轉分量和一個趨勢分量。

仔細觀察所有的PRC分量，還可發現ITD分解過程中幾乎所有的PRC分量在邊界處均出現了不同程度的畸變，與理論的分量信號存在差別，同時趨勢分量L理論上應為常數0，但實際所得L分量在波形中部出現細微波動，兩端出現較大的抖動，與理論分析嚴重不符，此即ITD分解過程中的邊界效應，應設法予以抑制或消除。

圖2 仿真信號的ITD分解結果

2 基于數據延拓的ITD邊界效應的抑制方法

2.1 自適應波形匹配延拓法

抑制邊界效應最有效的方法是對端點處的數據進行延拓，延拓的波形須符合原始數據在端點處的變化趨勢，否則易產生附加的極值點。在機械設備故障診斷中，原始振動信號在端點處的變化趨勢通常在信號的內部也會存在，故可用信號內部某一段與信號端點處變化趨勢相同的子波對端點以外的信號進行延拓，從而最大限度維持信號在端點處的變化規律。

(7)

自適應波形匹配方法的具體步驟如下[8]：

(1) 設x(t)為原始信號，確定x(t)最左端的兩個相鄰極值點，不妨設其分別為極大值點和極小值點，分別記為P0和P1，從起始點到P1的這段波形記為w0，設其長度為l；

(2) 設Emax為x(t)的極大值點集合，以Emax-{P0}中的每一個極大值點Pi作為參考點，計算該段相同長度的波形wi和w0的匹配度M(w0,w1,Pi)；記

M(w0,wi0,Pi0)=min{M(w0,wi,Pi),i=1,2,…}

(8)

(3) 若M(w0,wi0,Pi0)<α·l,其中α為一常數，則取wi0左側包含了一個極大值和極小值的子波, 作為原始x(t)左端的延拓, 延拓完畢，否則轉(4)；

(4) 直接指定端點處的極大和極小值：取原始信號最左端的兩個相鄰極大值點的均值作為左端點的極大值，取信號最左端的兩個相鄰極小值點的均值作為左端點的極小值，延拓完畢。

仍采用1.2節中的仿真信號，首先對該信號進行自適應波形匹配延拓，然后再進行ITD分解，最后對分解結果進行截斷以保留與原始信號相對應的部分，結果如圖3所示。

圖2中，趨勢分量的幅值約在-0.6～0之間，而圖3中趨勢分量的幅值約在-0.025～0.025之間，分解誤差顯著減小，故對原始信號進行自適應波形匹配延拓后再做ITD分解，可以有效地抑制ITD分解過程中的邊界效應，減少因為邊界效應產生的誤差。

圖3 仿真信號的ITD分解結果

2.2 鏡像延拓法

鏡像延拓相當于放一面鏡子在信號的端點，這樣在鏡子里就會出現一個與原始信號相反的信號。若原始信號為x(n)，其中n∈(1,2,3,…,N)，則鏡中的信號xj(m)為，

xj(m)=x(N+1-m)m∈(1,2,3,…,N)

(9)

將xj(m)與x(n)尾首相連，就得到了鏡像延拓后的信號。

仍采用前述仿真信號，用ITD分解經鏡像延拓后的信號，再截取與原始信號相對應的部分，得到結果如圖4所示。使用鏡像延拓法，ITD分解的邊界效應在一定程度上得到抑制，趨勢分量畸變的波動范圍較圖2有所改善，但其波動幅值遠大于圖3，邊界效應抑制不太理想。

圖4 仿真信號的ITD分解結果

2.3 多項式擬合延拓法

多項式擬合法是通過信號極大值或者極小值本身特征擬合多項式，來確定延拓的極值點[9]，具體方法如下：

(1) 在數據的極大值或者極小值的集合中，提取靠近兩端的幾個極值點。

(2) 用提取出來的極大值或是極小值擬合多項式。

(3) 選取在數據兩端之外的兩個時刻，用多項式計算出這兩個時刻所對應的函數值。將該函數值作為延拓的極值點，對應的時刻作為出現的時間。

同樣采用上述的仿真信號，得到的結果如圖5所示。ITD分解的邊界效應得到了有效的抑制，并且分解得到的誤差項非常小，說明這種方法對于模擬信號是有效的。

圖5 仿真信號的ITD分解結果

2.4 基于AR模型預測延拓法

AR模型可以根據數據點之間的相關性，預測數據的未來值[10]。對于一個時間序列x，若其第n個點x(n)的取值與其前m個點(x(n-1)…x(n-m))都有關，就可以用這m個點建立AR模型，來預測x(n)。其AR模型為

(10)

基于AR模型信號的邊界延拓的具體步驟為：

(1) 選取數據點建立AR模型。

(2) 根據模型預測數據，進行延拓。

采用前述同樣的仿真信號，首先對其進行AR模型延拓，然后再進行ITD分解，得到的結果如圖6所示。信號分解后的分解誤差較小，而且PRC分量只有略微的畸變，說明該方法抑制模擬信號邊界問題的效果較好。

圖6 仿真信號的ITD分解結果

2.5 反對稱周期延拓法

如果信號x(n)滿足下列情況，可以用反對稱周期延拓法[11]。

x(1)=0
x(2)-x(1)≠0

(11)

反對稱周期延拓法的具體步驟如下：

(1) 對于原始信號取負處理，即

y(n)=-x(n)

(12)

(2) 將得到的信號y(n)以n/2為中心對調，使信號y(n)的首尾對調，得到新的信號z(n)。

(3) 將z(n)與x(n)首尾相連得到延拓后的信號。

對前述仿真信號進行反對稱周期延拓之后再進行ITD分解，也得到了兩個PRC分量和一個趨勢分量，結果如圖7所示。雖然第二個PRC分量中還含有一定的邊界效應，該方法仍對于原始的邊界效應有抑制效果，使得ITD得到的分解誤差也較小。

圖7 仿真信號的ITD分解結果

2.6 延拓方法的數學模擬比較

同一個仿真信號，經過上述五種波形延拓方法之后再進行ITD分解，得到了五組不同的ITD分解結果。這五組結果都是由兩個PRC分量和一個趨勢分量組成。造成趨勢分量畸變的主要原因在于邊界效應，在該仿真信號中，趨勢分量的理論值應為常數0，故可用分解所得趨勢分量與理論值(常數0)之間的差異來定量表征不同延拓方法對ITD分解邊界效應的抑制效果。設系數ξ為

(13)

其中：L為ITD分解的趨勢分量，N為趨勢分量的數據長度。ξ用來定量表征趨勢分量L與理論值之間的差距。ξ值越大，說明趨勢分量L較大，邊界效應抑制效果差。ξ值越小，說明趨勢分量L較小，邊界效應抑制效果好。前述5種不同數據延拓方法中的ξ值如表1所示。

通過比較可以看出，ξ值最小為13.2，其對應方法為自適應波形匹配法，所以抑制ITD分解邊界效應最好的數據延拓方法為自適應波形匹配法。

表1 不同方法ξ值統計

3 工程應用

實驗系統由軸承實驗臺、壓電式加速度傳感器、數據采集儀、筆記本電腦組成，如圖8所示。將帶有人工模擬故障的軸承安裝在實驗臺上，進行實驗數據的采集和分析。

圖8 軸承故障模擬實驗臺

該實驗的滾動軸承型號為NU205EM，上有人工模擬的內圈點蝕故障，電機轉速為310 r/min，皮帶輪傳動比為1∶1。振動加速度信號的采樣頻率為5 000 Hz，采樣點數為8 192，經過計算故障的特征頻率為39.7 Hz。圖9(a)和(b)分別為滾動軸承內圈點蝕故障的時域波形及其幅值譜。

圖9 原始信號的時域波形及其幅值譜

由圖9可以看出，時域波形中出現等間隔的沖擊成分，幅值譜中亦含有明顯的邊頻帶，說明原始信號中存在著典型的調制現象。利用自適應波形匹配延拓法對該數據進行延拓，然后再進行ITD分解并截取，得到8個PRC分量和一個趨勢分量，其中前三個PRC分量包含有明顯的調制信息，如圖10所示。

圖10 前三個PRC分量的時域波形

對這三個PRC分量進行Hilbert包絡解調，得到其包絡譜，如圖11所示。在三個包絡譜中能明顯地看到軸承內圈的故障特征頻率39.7Hz，在前兩個PRC分量的包絡譜中還能清晰地看到39.7 Hz的倍頻成分，因此，可以肯定地判斷該軸承存在內圈故障，與預設的故障類型完全吻合。

圖11 前三個PRC分量的包絡譜

4 結 論

(1) ITD法可將一個復雜的非平穩信號分解為若干個固有旋轉分量和一個趨勢分量，便于融合其它現代信號處理方法進行故障特征的深層次挖掘。

(2) 自適應波形匹配法、基于AR模型的延拓方法、鏡像延拓方法、多項式延拓方法、反對稱周期延拓法，均可有效地抑制邊界效應，使分解所得的PRC分量端點處畸變小。通過仿真實驗比較，自適應波形匹配的抑制效果最佳。

(3) 將自適應波形匹配應用于滾動軸承內圈點蝕故障振動信號特征提取，實驗表明該方法能精確提取信號典型故障特征，在故障診斷領域具有很好的應用前景。

參 考 文 獻

[1]Yang Y, Peng Z K, Meng G. Characterize highly oscillating frequency modulation using generalized warblet transform[J]. Mechanical Systems and Signal Processing, 2012, 26: 128-140.

[2]Frei M G, Osorio I. Intrinsic time-scale decomposition：time-frequency-energy analysis and real-time filtering of non-stationary signals[J]．Proceedings of the Royal Society A, 2007, 463: 321-342.

[3]安金坤, 田 斌, 孫永軍,等. 一種基于ITD算法的直擴信號檢測算法[J]. 電子與信息學報，2010,32(5): 1178-1182.

AN Jin-kun, TIAN Bin, SUN Yong-jun, et al. An alagorithm for direct sequence spread spectrum signal detection based on intrinsic time-scale decomposition[J]. Journal of Electronics & information Technology, 2010, 32(5): 1178 -1182.

[4]顧小昕. 基于固有時間尺度分解的信號分析與干擾抑制技術研究[D]. 西安：西安電子科技大學, 2010.

[5]AN Xue-li, JIANG Dong-xiang, CHEN Jie, et al. Application of the intrinsic time-scale decomposition method to fault diagnosis of wind turbine bearing[J]. Journal of Vibration and Control, 2012,18(2): 240-245.

[6]LIN Jin-shan. Improved intrinsic time-scale decomposition method and its simulation[J].Frontiers of Manufacturing and Design Science II, 2011,121: 2045-2048.

[7]HAN Jing-tao, JIAO Si-hai, JIANG Zheng-yi. The de-noising algorithm based on intrinsic time-scale decomposition[J]. Advanced Materials Research, 2011,422: 347-352.

[8]邵晨曦, 王劍, 范金鋒,等. 一種自適應的EMD端點延拓方法[J]. 電子學報, 2007，35(10): 1944-1948.

SHAO Chen-xi, WANG Jian, FAN Jin-feng, et al. A self-adaptive method Dealing with the end Issue of EMD[J]. Acta Electronica Sinica, 2007，35(10): 1944-1948.

[9]許寶杰, 張建民, 徐小力,等.抑制EMD端點效應方法的研究[J]. 北京理工大學學報,2006,26(3): 196-200.

XU Bao-jie, ZHANG Jian-min, XU Xiao-li, et al. A study on the method of restraining the ending effect of empirical mode decomposition[J]. Transactions of Beijing Institute of Technology,2006,26(3): 196-200.

[10]胡勁松, 楊世錫. EMD方法基于AR模型預測的數據延拓與應用[J]. 振動、測試與診斷,2007,27(2): 116-120.

HU Jin-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-120.

[11]趙娜. HHT經驗模式分解的周期延拓方法[J]. 計算機仿真,2008,25(12): 346-350.

ZHAO Na. A periodic extension approach for HHT empirical mode decomposition[J]. Computer Simulation, 2008, 25(12): 346-350.