，

(石家莊鐵道大學(xué) 機(jī)械工程學(xué)院，河北 石家莊 050043)

齒輪是旋轉(zhuǎn)機(jī)械設(shè)備的重要組成部分，其運行狀態(tài)是否正常直接影響機(jī)械設(shè)備的整體性能。據(jù)統(tǒng)計，傳動機(jī)械中80%的故障是由齒輪引起的，旋轉(zhuǎn)機(jī)械中齒輪故障占其故障的10%左右[1]。當(dāng)齒輪發(fā)生故障時，故障信號受到外部環(huán)境影響而被噪聲信號淹沒，信噪比低，早期難以發(fā)現(xiàn)，加之故障沖擊引起機(jī)械系統(tǒng)的非線性振動，導(dǎo)致測得的齒輪故障信號也是非線性、非平穩(wěn)的，因此應(yīng)用傳統(tǒng)的信號處理方法無法獲得滿意的效果。

為了提取齒輪故障信息，消噪、降噪是信號處理的關(guān)鍵，許多專家學(xué)者提出一些有效的方法：傅里葉濾波、小波降噪，Hilbert-Huang變換等，但都有各自的局限性[2-4]。EMD作為一種自適應(yīng)的信號處理算法，它可以分解得到頻率由高到低的IMF分量，通過不同分量的重構(gòu)可以構(gòu)成不同種類的濾波器，從而達(dá)到降噪的目的。但是IMF分量的重構(gòu)一直沒有一個統(tǒng)一的標(biāo)準(zhǔn)，且重構(gòu)方法因人而異，得到的效果不盡相同[5]。

數(shù)學(xué)形態(tài)學(xué)[6]是Matheron和Serra共同提出的一種新的算法，該方法可以很好地抑制噪聲，其算法簡單、實用有效[7]，被廣泛應(yīng)用于信號處理領(lǐng)域。如郝如江等[8-9]提出了多尺度的形態(tài)濾波器的降噪算法，并在軸承故障特征提取方法中取得良好效果。但形態(tài)學(xué)濾波器存在如何選擇最優(yōu)結(jié)構(gòu)元素的問題。

針對上述問題，本文提出了基于模態(tài)相關(guān)法的分選準(zhǔn)則和基于峭度準(zhǔn)則優(yōu)化結(jié)構(gòu)元素尺度的自適應(yīng)形態(tài)濾波算法，二者結(jié)合更好地達(dá)到了降噪、消噪的目的。以仿真數(shù)據(jù)和齒輪故障實測信號進(jìn)行了驗證，結(jié)果表明該方法可以更好地提取出齒輪故障特征信息。

1 自適應(yīng)形態(tài)學(xué)濾波

1.1 形態(tài)學(xué)濾波器

形態(tài)濾波的原理是以信號的時域特征為基礎(chǔ)，通過選擇一定形狀的結(jié)構(gòu)元素與被分析的信號進(jìn)行適當(dāng)?shù)男螒B(tài)運算，從而達(dá)到降噪、突出信號特征的目的。它有4個基本算子：腐蝕、膨脹、形態(tài)開和形態(tài)閉[10]。由于論文的研究對象均是一維離散信號，因此僅給出一維的形態(tài)學(xué)算子定義。

設(shè)采樣得到的一維多值信號f(n)和一維結(jié)構(gòu)元素序列g(shù)(m)分別為定義域F=(0,1,2，…,N-1)和G=(0,1,2,…,M-1)上的離散函數(shù)，且N≥M，則f(n)關(guān)于g(n)的腐蝕和膨脹分別定義為

(fΘg)(n)=min[f(n+m)-g(m)],m∈G

(1)

(f⊕g)(n)=max[f(n-m)+g(m)],m∈G

(2)

式中，Θ為腐蝕算子；⊕為膨脹算子。

數(shù)學(xué)形態(tài)學(xué)的這4種基本運算都可以提取信號的輪廓信息，但其作用不同。腐蝕運算可以抑制正沖擊，平滑負(fù)沖擊；膨脹運算則相反；開運算用于濾除信號上方的峰值噪聲，去除信號邊緣的毛刺；閉運算則用于平滑或抑制信號下方的波谷噪聲，填補信號的漏洞和裂紋。也就是說，腐蝕和開運算的效果類似，膨脹和閉運算的效果則與之相反。實際應(yīng)用中，常用的形態(tài)濾波器是這幾種算子的組合。

(1)平均(AVG)濾波器

(AVG)(f)=(f·g+f°g)/2

(3)

(2)差值(DIF)濾波器

(DIF)(f)=f·g-f°g

(4)

開運算抑制正沖擊、閉運算抑制負(fù)沖擊，兩種濾波結(jié)果進(jìn)行平均可以起到降噪的目的。差值濾波器可以提取信號中的正、負(fù)沖擊。

1.2 基于峭度準(zhǔn)則優(yōu)化形態(tài)學(xué)結(jié)構(gòu)元素

結(jié)構(gòu)元素是形態(tài)學(xué)的基本算子，形態(tài)學(xué)濾波效果的好壞很大程度上取決于結(jié)構(gòu)元素，而結(jié)構(gòu)元素的選擇在于其形狀和尺寸(應(yīng)根據(jù)待分析信號的特點而定)的確定，結(jié)構(gòu)元素與待提取的特征形狀越相似，提取效果越好。結(jié)構(gòu)元素可以選擇為圓盤形、三角形、余弦、正方形、扁平形等。文獻(xiàn)[11]的研究結(jié)果表明，結(jié)構(gòu)元素形狀的選取對濾波效果影響不大，但整體來說余弦、半圓、三角優(yōu)于其它形狀的結(jié)構(gòu)元素。因此為了計算簡單，選取幅值為1的余弦結(jié)構(gòu)元素。當(dāng)其形狀和高度確定之后，結(jié)構(gòu)元素長度的最優(yōu)解就成了形態(tài)濾波中的重要工作。

對于一個離散變量x，歸一化的4階中心矩稱為峭度(kurtosis)，定義為

(5)

式中，E(x-μ)4是4階數(shù)學(xué)期望，μ是均值；σ是標(biāo)準(zhǔn)差。在一定范圍內(nèi)，故障信號的峭度值隨著故障的嚴(yán)重程度而增加，峭度的這一特性適用于齒輪故障特征提取。

本文利用峭度指標(biāo)(峭度越大，濾波誤差越小，則提取效果越好)作為衡量標(biāo)準(zhǔn)來自適應(yīng)地尋優(yōu)最佳尺度值。結(jié)合實驗工況和理論計算的齒輪故障頻率，將結(jié)構(gòu)元素尺度設(shè)定在[2,100]，即從2循環(huán)到100，找到峭度值最大時對應(yīng)的尺度值，即為最優(yōu)解。

2 EMD濾波理論

EMD分解方法是將一個待分析信號分解為若干個頻率由高到低的內(nèi)稟模態(tài)函數(shù)，分解的整個過程就是一個“篩分”的過程，表現(xiàn)出EMD本身的一個自適應(yīng)濾波特性。根據(jù)它的這個特點，我們可以將所需要的IMF分量重構(gòu)，從而得到某一個頻段的信息，這便相當(dāng)于一個濾波器。

在實際應(yīng)用中，EMD算法的濾波效果具有一定的局限性，具體表現(xiàn)為：

(1)EMD將信號分解為N個頻率由高到低的IMF分量，其中階數(shù)低的對應(yīng)高頻成分，一般認(rèn)為這些高頻成分是由噪聲主導(dǎo)的，重構(gòu)信號時直接去除，但是具體去除哪幾個，沒有相應(yīng)的準(zhǔn)則來判斷。

(2)即使人為地判斷出哪些是噪聲主導(dǎo)的模態(tài)分量，但對實際測得的信號而言，這些分量中也包含一些故障信號的高頻成分，直接去除這些分量會導(dǎo)致故障信息的缺失，影響后續(xù)齒輪故障的判斷。

針對上述問題，提出基于模態(tài)相關(guān)分選準(zhǔn)則和形態(tài)學(xué)濾波相結(jié)合的降噪方法，該方法的具體步驟以及用于齒輪故障信號的分析如下。

3 基于模態(tài)相關(guān)分選準(zhǔn)則和形態(tài)學(xué)濾波相結(jié)合的降噪方法

相關(guān)系數(shù)主要反映兩個函數(shù)的相關(guān)程度，相關(guān)系數(shù)越大，表明兩個函數(shù)之間的相關(guān)程度越高，相反，相關(guān)程度越小。經(jīng)過大量實驗發(fā)現(xiàn)，各階模態(tài)與原始信號的相關(guān)系數(shù)具有較好的數(shù)據(jù)穩(wěn)定性。本文將選用各階模態(tài)與原始信號的相關(guān)系數(shù)作為噪聲模態(tài)與信號模態(tài)的判斷準(zhǔn)則，稱之為模態(tài)相關(guān)分選準(zhǔn)則。其相關(guān)系數(shù)定義如下[12]

(6)

式中， cov()代表協(xié)方差。

文獻(xiàn)[12]用具體的數(shù)學(xué)推導(dǎo)過程證明了該準(zhǔn)則的合理性。另外，本文采用各個IMF分量的自相關(guān)函數(shù)對其進(jìn)行了進(jìn)一步地驗證。現(xiàn)就自相關(guān)函數(shù)的定義和相關(guān)特性介紹如下：

自相關(guān)函數(shù)作為信號時域中的一種統(tǒng)計特性，主要反映信號在兩個不同時刻時的相關(guān)程度，其定義為

Rx(t1,t2)=E(x(t1)x(t2))

(7)

實際應(yīng)用中一般對其歸一化處理，即

(8)

式中，x(t)為隨機(jī)信號；τ為時間差，即τ=t1-t2。

按式(7)和式(8)分別得到高斯白噪聲和一般信號的自相關(guān)函數(shù)，如圖1所示。其中，一般信號是由頻率均為10 Hz的正弦和余弦成分疊加而成，并將其歸一化處理。

圖1 自相關(guān)函數(shù)圖

由圖1(a)可知，高斯白噪聲的自相關(guān)函數(shù)在零點處為1，其它點處迅速衰減到零；由圖1(b)可知，一般信號的自相關(guān)函數(shù)在零點處取得最大值，其它點處不是迅速衰減到零，而是緩慢變化的，這與噪聲信號存在顯著的區(qū)別，因此可以根據(jù)各個IMF分量的自相關(guān)函數(shù)來判斷模態(tài)相關(guān)法分選準(zhǔn)則的正確性。

對于判斷出的噪聲起主導(dǎo)作用的IMF分量，傳統(tǒng)的EMD濾波算法是將其直接去除，事實上這些分量中還包含少量有用信號的高頻成分，為了盡可能多地保留有用信號的特征信息，本文利用自適應(yīng)形態(tài)學(xué)濾波器對其進(jìn)行降噪處理，然后將處理后的IMF分量與未處理的IMF分量進(jìn)行重構(gòu)得到濾波重構(gòu)信號。

綜上，所述方法流程如圖2所示：

(1)給定被噪聲污染的信號x(t)，將其進(jìn)行EMD分解，得到N個IMF分量；

(2)分別計算各個分量與原始信號的相關(guān)系數(shù)，做出相關(guān)系數(shù)分布圖，以及各個分量的自相關(guān)函數(shù)圖；

(3)從相關(guān)系數(shù)圖中找到噪聲主導(dǎo)分量與信號主導(dǎo)分量的分界點k，并根據(jù)自相關(guān)函數(shù)的特點，驗證所得分界點的正確性；

(4)對判斷出的噪聲主導(dǎo)的分量IMF1～I(xiàn)MFk，進(jìn)行自適應(yīng)形態(tài)學(xué)濾波處理，得到處理后的分量IMF1′～I(xiàn)MFk′；

圖2 基于EMD模態(tài)相關(guān)和形態(tài)學(xué)降噪方法流程圖

4 實驗分析研究

4.1 仿真數(shù)據(jù)分析

為了驗證所述方法的有效性，下面以仿真信號為例來說明，其中采樣頻率為1 024 Hz，采樣點數(shù)為1 000，假設(shè)原始信號為

圖3 仿真信號時域圖

x1(t)=sin(2π×10t)+cos(2π×10t)

(9)

該信號由頻率均為10 Hz的正弦信號和余弦信號疊加而成，如圖3所示。其加入噪聲之后的信號為

x(t)=x1(t)+x2(t)

(10)

圖4 仿真信號加噪時域圖

式中，x2(t)為高斯白噪聲；x(t)為加噪之后的信號(信噪比為2)，如圖4所示。其降噪的目的就是要得到一個最接近信號x(t)的估計信號x′(t)，使得x(t)-x′(t)最小。其具體仿真步驟及結(jié)果分析如下：

(1)EMD將信號分解為9個IMF分量，計算每個IMF分量與原始信號的相關(guān)系數(shù)，如圖5所示。最后，做出各個IMF分量的自相關(guān)函數(shù)，如圖6所示。

(2)由圖5可知，噪聲模態(tài)和信號模態(tài)之間的分界點為4，即前4個IMF分量是噪聲占主導(dǎo)的分量，后5個分量是信號占主導(dǎo)的分量，并且由圖6可知，前4個分量的自相關(guān)函數(shù)圖滿足噪聲自相關(guān)函數(shù)的衰減特性，從第5個分量開始滿足一般信號自相關(guān)函數(shù)的特性，進(jìn)一步證明所提方法的準(zhǔn)確性。

圖5 各IMF分量與原始信號的相關(guān)系數(shù)

圖6 前5個IMF分量的自相關(guān)函數(shù)

(3)由于前4個分量中也包含一些有用信號的高頻成分，所以對其進(jìn)行自適應(yīng)形態(tài)學(xué)濾波處理，通過峭度準(zhǔn)則優(yōu)化形態(tài)學(xué)結(jié)構(gòu)元素尺度，選擇峭度最大時對應(yīng)的尺度值，如圖7所示，選取尺度為13的余弦結(jié)構(gòu)元素。最后，將濾波之后的分量與后面的5個分量進(jìn)行重構(gòu)，得到濾波重構(gòu)信號，如圖8所示，由圖可以看出，重構(gòu)信號與原始信號之間相差并不大，幾乎重合，充分證明了所提方法在降噪方面的可行性。

圖7 不同尺度形態(tài)濾波的峭度值

圖8 濾波重構(gòu)信號與原始信號對比圖

4.2 齒輪故障信號的特征提取分析

為了進(jìn)一步證明所述方法的可行性，采用實測齒輪故障信號進(jìn)行分析研究，實驗采用江蘇千鵬診斷系統(tǒng)-QPZZ-II，其中輸入軸轉(zhuǎn)速為870 r/min，采樣頻率為5 kHz，采樣點數(shù)為65 536，故障齒輪(小齒輪)為3.2 mm的裂紋，齒輪參數(shù)和故障頻率如表1所示，采集到的齒輪故障信號的時域波形和頻譜圖如圖9所示。

表1 齒輪參數(shù)及故障頻率

圖9 實測信號時域波形和頻譜圖

由圖9可以看出，信號的頻率成分被淹沒在噪聲之中，無法找到信號的故障頻率，對其直接進(jìn)行自適應(yīng)形態(tài)學(xué)濾波處理，得到的時域波形和經(jīng)形態(tài)學(xué)差值濾波器解調(diào)之后的結(jié)果如圖10所示。從圖10(b)中可以找到信號的故障特征頻率(14.5 Hz)和倍頻及其邊頻信息，但幅值較小，二倍頻和三倍頻不是特別明顯。

采用本文所提方法進(jìn)行處理，首先將采集得到的信號進(jìn)行EMD分解，得到19個IMF分量，并利用公式(6)計算每一個分量與原信號的相關(guān)系數(shù)，做出相關(guān)系數(shù)圖，如圖11所示。由圖可知，噪聲模態(tài)和信號模態(tài)的分界點是第9個分量，并對每一個IMF分量做了自相關(guān)分析驗證了本文所提方法的準(zhǔn)確性。由于前9個分量中也可能包含齒輪故障信號的高頻成分，直接去除將導(dǎo)致故障頻譜信息不完整，所以引入自適應(yīng)形態(tài)濾波算法，將其進(jìn)行濾波處理，通過峭度準(zhǔn)則優(yōu)化形態(tài)學(xué)結(jié)構(gòu)元素尺度，如圖12所示，選取尺度為37的余弦結(jié)構(gòu)元素。最后，將濾波之后的分量與未處理的分量進(jìn)行重構(gòu)，得到濾波重構(gòu)信號，然后對其進(jìn)行Hilbert包絡(luò)解調(diào)，結(jié)果如圖13所示。由圖可知，本文所述方法有效地提取出了齒輪的故障頻率、倍頻以及邊頻信息。邊頻信息數(shù)量多、幅值不高、均勻分布，倍頻成分比較明顯，幅值較高、數(shù)量多，與所設(shè)的齒輪裂紋故障相吻合。并與圖10(b)對比發(fā)現(xiàn)，本文方法明顯比直接進(jìn)行形態(tài)學(xué)濾波的效果更佳，得到的頻譜圖中也更能凸出齒輪的故障特征頻率，其前三階故障頻率幅值分別提高107%、243%和190%(見表2)。充分證明，本文所提基于模態(tài)相關(guān)分選準(zhǔn)則和自適應(yīng)形態(tài)學(xué)濾波的降噪新方法可以更好地對信號進(jìn)行降噪處理，可以更有效地應(yīng)用于齒輪故障診斷研究。

圖10 形態(tài)濾波處理后的結(jié)果

圖11 各階分量與原始齒輪信號的相關(guān)系數(shù)

圖12 不同尺度形態(tài)濾波的峭度值

圖13 濾波重構(gòu)解調(diào)后的頻譜圖

表2 故障特征頻率幅值比較

5 結(jié)論

將EMD模態(tài)相關(guān)法和形態(tài)學(xué)濾波相結(jié)合，首先對信號進(jìn)行EMD分解，計算分解得到的IMF分量與原始信號的相關(guān)系數(shù)，找到噪聲主導(dǎo)模態(tài)與信號主導(dǎo)模態(tài)的分界點，并利用各個IMF分量的自相關(guān)函數(shù)驗證該分界點的準(zhǔn)確性。其次將選取的噪聲模態(tài)利用自適應(yīng)形態(tài)學(xué)濾波進(jìn)行降噪處理。最后，將降噪后的模態(tài)與剩余模態(tài)重構(gòu)，并對其進(jìn)行Hilbert包絡(luò)解調(diào)得到故障特征頻率。仿真數(shù)據(jù)和齒輪裂紋故障實驗分析均表明該方法可以更加有效地提取齒輪故障特征頻率。主要結(jié)論有：

(1)針對IMF重構(gòu)中噪聲模態(tài)和信號模態(tài)分界點的選取問題，提出了基于相關(guān)系數(shù)的模態(tài)分選準(zhǔn)則，并利用各個IMF分量的自相關(guān)函數(shù)進(jìn)行驗證，得出該準(zhǔn)則的準(zhǔn)確性。

(2)針對形態(tài)學(xué)濾波器中結(jié)構(gòu)元素的選取問題，提出了基于峭度準(zhǔn)則優(yōu)化結(jié)構(gòu)元素的自適應(yīng)算法，實驗結(jié)果表明該方法可以提高形態(tài)學(xué)的濾波能力。

(3)將EMD模態(tài)相關(guān)分選準(zhǔn)則和自適應(yīng)形態(tài)學(xué)濾波相結(jié)合能夠明顯地從齒輪振動信號中提取出故障特征頻率，比直接進(jìn)行形態(tài)學(xué)濾波效果更佳。