基于固有時間尺度分解與最大相關峭度解卷積的軸承故障特征提取

王銀玲 尹顯明 李華聰

(1. 西南科技大學工程技術中心 四川綿陽 621010； 2. 西北工業大學動力與能源學院 西安 710072)

滾動軸承廣泛應用在旋轉機械中，受負荷、摩擦、沖擊等非線性因素影響會出現損傷。通常由于滾動軸承工作環境復雜，傳感器采集到的信號含有大量的背景噪聲，造成損傷信號的非線性、非平穩、沖擊性特征被淹沒。因此，如何通過對損傷過程中產生的非線性、非平穩隨機信號進行分析，準確檢測出設備的故障信號是滾動軸承故障檢測研究的熱點及難點[1-3]。

在軸承故障診斷中，損傷信號通常被看作是周期性沖擊信號，其與系統響應的卷積即為傳感器測量的輸出信號。若不考慮噪聲的影響，可以尋找一個逆濾波器，通過傳感器輸出的觀測信號，利用解卷積恢復出原始的故障信號。Wiggin[4]結合熵值最小卷積濾波思想，提出最小熵解卷積方法(Minimum Entropy Deconvolution，MED)，在解卷積過程中使得峭度最大化，該方法又稱為最大峭度解卷積。但MED只能處理單一沖擊信號，并不適合旋轉機械的周期性沖擊信號的處理。

故障信號通過EMD分解[5]過程中存在端點效應以及模態混疊現象。為了解決這些問題，學者們提出了不同的解決方案。陳哲等[6]提出一種集合經驗模態分解(Ensemble Empirical Mode Decomposition，EEMD)[7]與多尺度排列熵(Multiscale Permutation Entropy，MPE)的艦船輻射噪聲復雜度特征提取方法，解決復雜海洋環境中艦船輻射噪聲的特征提取問題。耿讀艷等[8]為了有效識別心沖擊(Ballistocardiogram)信號，利用噪聲自適應完備總體平均經驗模態分解(Complete EEMD with Adaptive Noise，CEEMDAN)的自適應降噪優勢，提出一種基于自適應噪聲的CEEMDAN聯合排列熵(PE)的BCG降噪方法，降噪后信號的幅頻特性得到明顯改善且信噪比較傳統方法有明顯提高。趙榮珍[9]提出一種基于經驗小波變換、多尺度排列熵、GG聚類算法相結合的故障診斷方法。劉興教[10]通過EEMD和最大相關峭度解卷積(Maximum Correlated Kurtosis Deconvolution，MCKD)提取柔性薄壁軸承故障特征。

由于以上方法不能產生適當的旋轉分量的余量，且在信號分解過程中為了降低噪聲及模態混疊，在模態分解過程引入了隨機噪聲，但是引入隨機噪聲后不能從根本上消除模態混疊，并可能會產生虛假分量。為了解決噪聲及模態混疊對故障信號的不利影響，本文提出一種改進的固有時間尺度分解(Intrinsic Time-scale Decomposition, ITD)與最大相關峭度解卷積(MCKD)[11]相結合的故障信號特征提取方法。首先通過ITD將原始信號分解為一系列固有旋轉分量和一個剩余的單調趨勢分量，然后選擇與原信號相似度高的旋轉分量作為濾波信號的主要成分，利用粒子群算法優化MCKD的參數，使得峭度值最大，最后對信號進行包絡解調，提取滾動軸承信號的故障特征。

1 改進的固有時間尺度分解

固有時間尺度分解(ITD)是由Frei等[12]提出的一種非平穩信號的時頻分析方法。ITD可以實時提取信號的瞬時特征，將信號分解為一系列適當的旋轉分量和一個趨勢分量。對于t>0的待分解信號Xt，ITD的分解過程是定義一個運算算子L，從Xt提取出基線信號Lt，再通過Xt-Lt得到旋轉分量Ht，用公式表示為：

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

(1)

ITD分解[13]首先找出Xt上所有的局部極值點所對應的時間集合{τk,k=1,2,3…}，時間集合包含信號起始時間點。為了表示方便，用Xk和Lk代替X(τk)和L(τk)。假設Ht和Lt在時間區間[0,τk]上有定義，Xt中時間t∈[0,τk+2]，這樣就可以在時間(τk,τk+1]上定義一個分段線性的基提取算子L：

t∈(τk,τk+1]

(2)

其中：

(1-α)Xk+1

(3)

式中，α∈(0,1)，通常取值為0.5。以這種方式構造的基線信號Lt，可以使其在極值之間保持單調性，并在極值點的包絡內。由式(2)和式(3)得到基線點信號后，通過對其進行線性插值，得到基信號Lt。用待分解信號減去基信號Lt，就得到一個適當的旋轉分量Ht。將提取出的基信號Lt作為原始信號，繼續進行ITD分解，直到所分解的基信號小于某設定值時，ITD分解結束。最終ITD分解可以表示為：

H1+H2+…HP+LP

(4)

式中，Hi表示第i(i=1,2,…p)個適當的旋轉分量，而Lp表示P次后的殘余基線分量。

由于ITD分解過程中采用線性插值，線性插值只顧及其附近兩點的影響，函數上兩點之間的近似隨著所近似函數的二階導數增大而逐漸變差，為了解決這個問題，本文采用Akima插值方法。Akima插值方法是一種5點求導分段三次多項式插值算法，其5點求導法利用所求點與前后連續兩點組成5個點，然后利用點之間延長線得到交點，最后得到所求點的導數[14]。

為了驗證基于Akima插值ITD分解的有效性，給出一個合成的數字信號，其數學表達式為：

y(t)=y1(t)+y2(t)+0.07·random

(5)

其中random為隨機信號：

y1(t)=[1+0.5cos(8πt)]·cos[200πt+

2cos(10πt)]

y2(t)=0.8sin(πt)sin(30πt)

圖1為原始合成信號及其分量波形，圖2為通過EMD分解得到前3個主要IMF分量的時域波形，圖3為通過ITD分解后得到的特征波形。可以看到ITD分解方法要優于EMD分解方法，驗證了該方法的有效性。

圖1 原始合成信號及其分量Fig.1 Original composite signal and its components

圖2 合成信號EMD分解后前3個IMF分量Fig.2 The first three IMF components ofcomposite signal after EMD decomposition

圖3 合成信號ITD分解后前3個分量Fig.3 The first three components of the composite signal after ITD decomposition

2 最大相關峭度解卷積

傳感器采集到的滾動軸承故障信號x(n)可以表示為：

x(n)=h(n)·y(n)+e(n)

(6)

式中，h(n)為系統的沖擊響應，y(n)為系統故障信號，e(n)為噪聲。在不考慮噪聲的情況下，尋找最優濾波器f=[f1,f2,…fL]T，對系統采集的信號x(n)進行解卷積處理，可還原出故障信號y(n)，即

y(n)=f(n)·x(n)

(7)

熵表示系統內信息的混亂程度，熵值越大則系統越難預測。MED方法將熵值最小與卷積濾波思想進行了結合，對傳感器采集到的原始信號進行解卷積，以反卷積信號x(n)峭度作為目標函數求解最優濾波器，以突出信號的沖擊成分。但是濾波器選擇不當將會提取出虛假成分，因此出現了以相關峭度(Correlation Kurtosis，CK)為最優指標的MCKD，相關峭度的計算式為：

(8)

(9)

式中L為濾波器長度，而f=[f1,f2,…fL]T。對目標函數求導，使得：

(10)

最終得到濾波器為：

(11)

其中：

(12)

r=[0,T,2T,…,mT]

(13)

(14)

從以上分析可知，通過MCKD將傳感器采集信號x(n)進行解卷積濾波時，需要設定周期T、偏移量M、濾波器長度L以及迭代次數N。

3 試驗驗證

實驗數據來源于某變速器軸承的外圈故障，變速器轉速為n=923 r/min，軸承節徑D=57.50 mm，滾動體直徑d=14.22 mm，壓力角α=0°，滾動體個數Z=7。外圈故障的特征頻率滿足以下關系：

(15)

通過上式可以計算出軸承外圈故障頻率fr=17.43 Hz。

采集變速器軸承外圈故障信號，采樣頻率為40 kHz，采集數據長度為84×103。首先對信號做去均值處理，處理后信號的時域與頻域特征如圖4所示。

圖4 原始故障信號時域及頻域特征Fig.4 Time domain and frequency domain characteristics of the original fault signal

圖5為原始信號EMD分解后一階IMF分量的包絡譜，由圖5可以看出，故障信號的特征頻率在35.4 Hz，105 Hz處都有較大的峰值，但是由于噪聲的影響，不能判斷是否含有沖擊成分，也就是故障信號的特征包絡譜。所以通過對原始信號的包絡譜分析，不能有效提取故障特征頻率。

圖5 信號EMD分解后一階IMF信號包絡譜Fig.5 First-order IMF signal envelope spectrum after signal EMD decomposition

接下來需要去除系統的噪聲，突出損傷沖擊信號特征。對信號進行基于Akima插值ITD分解，分解后前4個分量的時域波形和其頻譜如圖6所示。

圖6 ITD分解后各分量波形及頻域特征Fig.6 Waveform and frequency domain characteristics of each component after ITD decomposition

ITD分解后各分量的能量比、峭度、峭度因子和相似性的統計特征如圖7所示。由圖7可知，第一階次的能量比最大，是信號的主要成分；第一、二階次的峭度最大，證明其變化最為劇烈；通過相似性分析可知第一階次與原始信號最接近，故將第一階分量作為提取信號。

圖7 ITD分解后各分量的統計特征Fig.7 Statistical characteristics of each component after ITD decomposition

通過ITD分解得到信號的主要成分后，再對信號進行MCKD處理以得到信號的沖擊成分。通過MCKD對信號進行解卷積時，需要設定周期T、偏移量M、濾波器長度L及迭代次數N。其中：周期T計算值并不一定最優，與實際還有偏差；而濾波器長度L直接影響著整個計算存儲數據大小；偏移量M一般取1～7之間；迭代次數太小有可能不是最優解，太大影響計算效率，一般選取30。以上參數的選取直接影響著峭度，為了獲得最優參數，本文以最大峭度為目標函數，采用粒子群算法(PSO)對上述參數進行優化，選取最優組合。最終選擇濾波器尺寸L=100、最終迭代次數N=30、解卷積周期為T=30、偏移量M=5時取得最優解。圖8比較了MCKD處理前后的信號譜峭度，圖9為MCKD處理后信號的包絡譜。圖9清晰顯示出故障頻率為17.09，外圈故障頻率及其倍頻上出現峰值，能夠清晰地診斷出設備的故障特征。

圖8 MCKD處理前后信號譜峭度Fig.8 Signal spectral kurtosis before and after MCKD processing

圖9 MCKD處理后信號的包絡譜Fig.9 Envelope spectrum of signal processed by MCKD

4 結束語

對于非平穩性、非線性的軸承故障信號，本文提出了基于固有時間尺度分解和相關峭度解卷積相結合的特征提取方法。借助固有尺度分解，提取滾動軸承故障信號的主要特征，為了使得到的波形更加光滑平穩，在尺度分解中采用Akima插值方法代替線性插值。借助MCKD提取信號的沖擊成分過程中，為了選擇合適的參數，通過粒子群算法對參數進行優化，以最大峭度為目標函數，得到最優參數，最后對MCKD后信號進行包絡譜分解，可以清晰地提取出故障特征頻率。實驗結果與理論一致，有利于對噪聲環境下軸承的狀態進行特征提取。