唐貴基，丁 傲,，王曉龍，張 曄

(1.華北電力大學(xué)機(jī)械工程系，河北 保定 071003；2.國(guó)網(wǎng)吉林省電力有限公司長(zhǎng)春供電公司，吉林 長(zhǎng)春 130021)

0 引言

滾動(dòng)軸承作為旋轉(zhuǎn)機(jī)械中的重要組成部分，支撐著整個(gè)設(shè)備的安全、可靠運(yùn)行[1]。因其常常處于高溫重載的惡劣工作環(huán)境之中，易出現(xiàn)表面損傷。軸承損傷會(huì)引發(fā)整個(gè)機(jī)械系統(tǒng)異動(dòng)，造成設(shè)備停機(jī)，威脅人身安全。軸承的早期故障信號(hào)中沖擊成分微弱，而實(shí)際運(yùn)行中背景噪聲廣泛存在，且信號(hào)在傳輸中會(huì)被衰減，所以采集到的信號(hào)信噪比低，故障特征提取困難。在準(zhǔn)確識(shí)別軸承故障類型方面，從原始振動(dòng)信號(hào)中提取微弱故障特征信息一直是軸承故障診斷領(lǐng)域的熱點(diǎn)[2]。

目前，經(jīng)驗(yàn)?zāi)B(tài)分解、局部均值分解、固有時(shí)間尺度分解等頻帶劃分信號(hào)處理算法，難以精細(xì)分析故障信號(hào)高頻區(qū)域，且存在虛假分量、曲線失真等問題；共振稀疏分解、經(jīng)驗(yàn)小波變換、變分模態(tài)分解等非遞歸式信號(hào)處理算法，分解效果受復(fù)雜的參數(shù)設(shè)置制約。2014年，Bonizzi等[3]提出奇異譜分解(singular spectrum decomposition，SSD)。該方法通過重新定義軌跡矩陣來增強(qiáng)信號(hào)的振蕩分量：根據(jù)信號(hào)特征自適應(yīng)選取軌道矩陣嵌入維數(shù)，對(duì)軌跡矩陣主成分進(jìn)行奇異值分解及特征重組，將信號(hào)分割為由高頻至低頻的若干奇異譜分量(singular spectrum component，SSC)。SSD算法為非平穩(wěn)非線性機(jī)械故障信號(hào)處理提供了新思路，目前已有學(xué)者對(duì)其進(jìn)行研究。文獻(xiàn)[4]將SSD算法與一維卷積神經(jīng)網(wǎng)絡(luò)(1 dimension convolutional neural network,1DCNN)算法相結(jié)合，診斷恒速工況下滾動(dòng)軸承故障。文獻(xiàn)[5]將SSD算法應(yīng)用于齒輪箱復(fù)合故障診斷中。文獻(xiàn)[6]將SSD算法與奇異值分解相結(jié)合，用于齒輪故障特征提取。文獻(xiàn)[7]將SSD算法用于風(fēng)電機(jī)組驅(qū)動(dòng)系統(tǒng)多通道故障診斷。

然而，SSD算法需人為設(shè)置奇異譜分量個(gè)數(shù)K。若K值過小，信號(hào)分解不充分；若K值過大，會(huì)出現(xiàn)虛假分量。針對(duì)該問題，本文提出了自適應(yīng)奇異譜分解(adaptive singular spectrum decomposition，ASSD)。ASSD采用合成峭度與斯皮爾曼等級(jí)相關(guān)系數(shù)(Spearman rank correlation coefficient,SRCC)作為聯(lián)合判據(jù)，實(shí)現(xiàn)奇異譜分量個(gè)數(shù)K的自適應(yīng)確定。

滾動(dòng)軸承故障信號(hào)周期性沖擊提取可視為源信號(hào)與信道解卷積的過程。Wiggins提出最小熵解卷積(minimum entropy deconvolution，MED)算法。但MED算法傾向于單點(diǎn)脈沖而非周期性脈沖的提取。此后，McDonald提出了最大相關(guān)峭度解卷積(maximum correlate kurtosis deconvolution，MCKD)算法。相較MED算法而言，MCKD算法增加了局部有限脈沖數(shù)。2017年，McDonald進(jìn)一步提出了多點(diǎn)最優(yōu)最小熵反褶積(multipoint optimal minimum entropy deconvolution adjusted，MOMEDA)算法，以D-范數(shù)最大為目標(biāo)，根據(jù)預(yù)先估計(jì)的故障頻率構(gòu)建目標(biāo)向量。與MED算法和MCKD算法相比，MOMEDA算法在對(duì)周期性沖擊信號(hào)特征提取方面的優(yōu)勢(shì)更為明顯[8]。

2014年，O’Tooley提出頻率加權(quán)能量算子(frequency weighted energy operator，F(xiàn)WEO)算法，利用信號(hào)的諧波特性，通過計(jì)算信號(hào)導(dǎo)數(shù)的包絡(luò)，得到信號(hào)的瞬時(shí)能量估計(jì)。與希爾伯特變換和傳統(tǒng)能量算子相比，F(xiàn)WEO算法的計(jì)算復(fù)雜性低，具有更好的抗干擾特性[9]。

鑒于上述分析，針對(duì)強(qiáng)噪聲背景下的滾動(dòng)軸承故障診斷問題，本文提出1種ASSD-MOMEDA-FWEO相結(jié)合的方法：利用ASSD算法對(duì)原始信號(hào)進(jìn)行預(yù)處理，從中篩選出信噪比較高的最佳奇異譜分量；利用MOMEDA算法和FWEO算法作為后處理策略，進(jìn)一步放大故障特征，實(shí)現(xiàn)滾動(dòng)軸承微弱故障的準(zhǔn)確判定。

1 基本原理介紹

1.1 SSD算法

SSD作為一種新型信號(hào)處理方法，能將非平穩(wěn)非線性信號(hào)從高頻至低頻依次分解為一系列奇異譜分量和殘余分量的線性疊加。具體分解過程如下。

為更好地估計(jì)主峰帶寬Δf，利用基于3個(gè)高斯函數(shù)混合的譜模型擬合功率譜密度輪廓，即：

(1)

式中：Ai為第i個(gè)高斯函數(shù)的幅值；ui為位置；σi為寬度；θ=[Aσ]T為參數(shù)向量且滿足A=[A1A2A3]和σ=[σ1σ2σ3]。

(2)

當(dāng)比重小于給定閾值時(shí)，則終止迭代。否則，將殘余分量作為原始信號(hào)重復(fù)迭代過程，直至滿足停止條件。最終獲得K個(gè)SSC分量。原始信號(hào)可表示為如下形式：

(3)

1.2 MOMEDA算法

MOMEDA算法是在MED算法的基礎(chǔ)上開發(fā)的。濾波后的實(shí)際信號(hào)yn為：

(4)

(5)

MOMEDA算法利用目標(biāo)向量和D-范數(shù)重新定義了多點(diǎn)D范數(shù)(multi D-norm，MDN)指標(biāo)：

(6)

將最小熵解卷積問題變?yōu)榍蠼舛帱c(diǎn)D范數(shù)最大化問題：

(7)

式中：t為確定沖擊成分的位置和權(quán)重的目標(biāo)向量。

對(duì)式(7)求解，可得：

(8)

取其特解作為一組最優(yōu)濾波器：

(9)

將式(9)代入式(4)，得到原始沖擊信號(hào)為：

(10)

1.3 FWEO算法

FWEO算法的基本原理如下：

(11)

式中：S()為信號(hào)的復(fù)包絡(luò)；H()為希爾伯特變換；φ[x(t)]為頻率加權(quán)能量算子。

振動(dòng)信號(hào)的調(diào)幅-調(diào)頻信號(hào)的表達(dá)式為：

x(t)=Acos(ωt+φ)

(12)

根據(jù)式(12)可知：

(13)

(14)

將式(13)、式(14)代入式(11)，可得：

φ[x(t)]=[Aωsin(ωt+φ)]2+

[Aωcos(ωt+φ)]2=A2ω2

(15)

2 自適應(yīng)奇異譜分解方法

2.1 合成峭度

軸承發(fā)生故障時(shí)，故障表面與其他元件發(fā)生周期性接觸。此時(shí)，振動(dòng)信號(hào)中故障成分具有沖擊性和循環(huán)平穩(wěn)性雙重特性。時(shí)域峭度對(duì)沖擊性敏感，忽略了循環(huán)平穩(wěn)性，而包絡(luò)譜峭度可評(píng)價(jià)信號(hào)的循環(huán)平穩(wěn)性。本文采用包絡(luò)譜峭度和時(shí)域峭度構(gòu)成得到的合成峭度指標(biāo)[10]。其定義如下：

KE=KES·Kku

(16)

(17)

式中：KE為合成峭度；Kku為峭度；KES為包絡(luò)幅值峭度；SE(j)為采樣點(diǎn)j的信號(hào)包絡(luò)譜；p為包絡(luò)譜的采樣點(diǎn)數(shù)。

2.2 斯皮爾曼等級(jí)相關(guān)系數(shù)

SRCC是用于評(píng)價(jià)2個(gè)變量之間相關(guān)性的1種非參數(shù)、無(wú)分布的秩統(tǒng)計(jì)度量方法[11]，將樣本轉(zhuǎn)換成等級(jí)，無(wú)需考察樣本規(guī)模及總體分布特性，具有快捷、穩(wěn)健的特點(diǎn)。首先，把數(shù)列x={x1，x2，…，xn}按升序或降序排列，得到排序數(shù)列a={a1，a2，…，an}，將x內(nèi)各元素xi在數(shù)列a中的位置記為ri，稱其為元素xi的秩次，從而得到秩次數(shù)列r。將另一數(shù)列y={y1，y2，…，yn}按同樣方式排列，得到排序數(shù)列b={b1，b2，…，bn}，繼而得到y(tǒng)對(duì)應(yīng)的秩次數(shù)列s。將數(shù)列r和s內(nèi)每個(gè)元素對(duì)應(yīng)地相減，得到秩次差數(shù)列d={d1，d2，…，dn}。則SRCC計(jì)算表達(dá)式如下：

(18)

式中：ρ為SRCC值；n為數(shù)列點(diǎn)數(shù)，對(duì)應(yīng)采樣點(diǎn)數(shù)。

式(18)中，分子為2個(gè)序列的誤差和，反映2個(gè)變量之間的差異；分母為與序列長(zhǎng)度相關(guān)的1個(gè)常數(shù)。

2.3 具體實(shí)現(xiàn)流程

通過研究發(fā)現(xiàn)，SSD算法需人為設(shè)置奇異譜分量個(gè)數(shù)K。若K設(shè)置得過小，信號(hào)未能分解完全，蘊(yùn)含敏感特征的分量無(wú)法被分離，故障信息將被掩蓋；若K設(shè)置得過大，則易出現(xiàn)虛假分量，同時(shí)運(yùn)算量加大。針對(duì)該問題，本文提出1種自適應(yīng)奇異譜分解方法，采用合成峭度與SRCC作為聯(lián)合判據(jù)，可實(shí)現(xiàn)奇異譜分量個(gè)數(shù)K的自適應(yīng)確定。具體過程如下。

首先，初始化分解數(shù)量K=2，應(yīng)用SSD將原始信號(hào)分解為若干個(gè)SSC分量。計(jì)算各SSC分量與原始信號(hào)的SRCC，提取最小值，記為ρm1。將ρm1作為停止迭代條件，即判斷ρm1是否小于0.1。若小于0.1，即認(rèn)為通過SSD算法分解產(chǎn)生了與原信號(hào)不相關(guān)的分量，視為原始信號(hào)被過度分解，則輸出(K-1)為最佳奇異譜分量個(gè)數(shù)。

其次，計(jì)算各SSC分量的合成峭度，按合成峭度最大原則從各奇異譜分量中剝離出包含敏感特征信息的最佳奇異譜分量SSC-M。計(jì)算SSC-M分量與其余各分解分量的SRCC值，提取最大值，記為ρm2。將ρm2作為第二層停止迭代條件，即判斷ρm2是否大于0.1。若大于0.1，則認(rèn)為SSD分解后的信號(hào)中已產(chǎn)生與SSC-M相關(guān)的分量，視為原始信號(hào)被過度分解，則輸出(K-1)為最佳奇異譜分量個(gè)數(shù)。

若ρm1大于0.1且ρm2小于0.1，則K=K+1。重復(fù)上述過程，最終可自適應(yīng)地得到SSD的最佳奇異譜分量個(gè)數(shù)。奇異譜分量個(gè)數(shù)自適應(yīng)選取流程如圖1所示。

圖1 奇異譜分量個(gè)數(shù)自適應(yīng)選取流程

3 故障特征提取

基于ASSD-MOMEDA-FWEO的軸承微弱故障特征提取流程如圖2所示。

圖2 故障特征提取流程

故障特征提取流程具體步驟如下。

①采用ASSD分解原始信號(hào)，得到一系列SSC分量。

②按合成峭度最大原則，從ASSD處理結(jié)果中篩選出蘊(yùn)含豐富特征信息的最佳奇異譜分量。

③利用MOMEDA算法對(duì)最佳奇異譜分量進(jìn)行解卷積處理，實(shí)現(xiàn)微弱故障特征信息的進(jìn)一步強(qiáng)化。

④執(zhí)行FWEO運(yùn)算，計(jì)算所得解卷積信號(hào)的瞬時(shí)能量并進(jìn)行傅里葉變換，得到其瞬時(shí)能量譜。

⑤從瞬時(shí)能量譜中識(shí)別峰值明顯的頻率成分，與理論故障特征頻率進(jìn)行對(duì)比，最終判定軸承運(yùn)行狀態(tài)。

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

為評(píng)估本文所提ASSD-MOMEDA-FWEO方法的故障特征辨識(shí)能力，采用以下模型[12]對(duì)強(qiáng)噪聲背景下的滾動(dòng)軸承內(nèi)圈故障進(jìn)行模擬。其中，采樣頻率為12 kHz。

(19)

式中：s(t)為軸承故障周期性沖擊信號(hào)；T為故障周期，取0.006 3 s，則故障頻率為159.7 Hz;A0為沖擊信號(hào)幅值，取1；B為衰減系數(shù)，取800；fr為轉(zhuǎn)頻，取29.5 Hz；fn為局部缺陷沖擊所激發(fā)的諧振頻率，取3 500 Hz；u(t)為單位階躍函數(shù)；n(t)為高斯白噪聲，添加噪聲后的仿真信號(hào)信噪比為-14.8 dB。

仿真信號(hào)時(shí)域波形及包絡(luò)譜如圖3所示。由圖3可知，內(nèi)圈故障特征頻率及其倍頻均淹沒在強(qiáng)噪聲中，表明傳統(tǒng)包絡(luò)方法對(duì)該故障仿真信號(hào)特征提取失效。

圖3 仿真信號(hào)時(shí)域波形及包絡(luò)譜

根據(jù)ASSD算法實(shí)現(xiàn)流程，當(dāng)K為2～3時(shí)，ρm1均大于0.1，且ρm2均小于0.1，因此繼續(xù)分解；當(dāng)K=4時(shí)，ρm1為0.075 1，小于0.1，因此停止分解，輸出最佳奇異譜分量個(gè)數(shù)為3。

本文算法仿真信號(hào)處理結(jié)果如圖4所示。

圖4 本文算法仿真信號(hào)處理結(jié)果

原始仿真信號(hào)經(jīng)ASSD處理后，得到3個(gè)分量，計(jì)算各SSC分量的合成峭度值如圖4(a)所示。其中，SSC-3分量合成峭度值最大，將其作為最優(yōu)分量。對(duì)SSC-3分量進(jìn)行MOMEDA處理后的解卷積信號(hào)時(shí)域波形如圖4(b)所示，周期性脈沖成分得以加強(qiáng)，通過FWEO進(jìn)一步追蹤解卷積信號(hào)的瞬時(shí)能量，并計(jì)算得到解卷積信號(hào)瞬時(shí)能量譜如圖4(c)所示。其中，故障特征頻率基頻～6倍頻成分均清晰可辨。

EMD算法仿真信號(hào)處理結(jié)果如圖5所示。

圖5 EMD算法仿真信號(hào)處理結(jié)果

為驗(yàn)證本文方法在強(qiáng)噪聲環(huán)境下故障信息剝離及周期性沖擊特征增強(qiáng)的優(yōu)勢(shì)，采用經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition,EMD)方法進(jìn)行對(duì)比研究。仿真信號(hào)經(jīng)過EMD處理后得到8個(gè)IMF分量，計(jì)算各IMF分量的合成峭度值，如圖5(a)所示。其中，IMF-3分量的合成峭度值最大，因此將IMF-3作為最優(yōu)分量。圖5(b)為IMF-3的包絡(luò)譜，無(wú)法提取出故障特征頻率。

5 試驗(yàn)數(shù)據(jù)分析

為進(jìn)一步探究本文方法在實(shí)際工況下的有效性，采用美國(guó)凱斯西儲(chǔ)大學(xué)軸承數(shù)據(jù)中心的試驗(yàn)數(shù)據(jù)進(jìn)行分析。滾動(dòng)軸承故障模擬試驗(yàn)臺(tái)包括電機(jī)、扭矩傳感器、功率計(jì)和控制元件4大部分。驅(qū)動(dòng)端軸承型號(hào)為SKF6205，軸承節(jié)徑為39.04 mm，滾動(dòng)體直徑為7.94 mm，滾動(dòng)體為9個(gè)，接觸角為0°。試驗(yàn)軸承故障通過電火花在內(nèi)、外圈植入，故障特征頻率計(jì)算表達(dá)式為：

(20)

(21)

5.1 外圈故障信號(hào)分析

外圈故障試驗(yàn)信號(hào)時(shí)域波形與包絡(luò)譜如圖6所示。

圖6 外圈故障試驗(yàn)信號(hào)時(shí)域波形與包絡(luò)譜

選取電機(jī)驅(qū)動(dòng)端處軸承外圈點(diǎn)蝕直徑為0.533 4 mm的故障數(shù)據(jù)進(jìn)行分析：試驗(yàn)采樣頻率為12 kHz，電機(jī)轉(zhuǎn)速為1 718 r/min，轉(zhuǎn)頻fr=28.63 Hz。根據(jù)式(20)計(jì)算出外圈理論故障特征頻率f0=102.6 Hz。由圖6可知，外圈故障特征頻率基頻可被識(shí)別，但譜圖中無(wú)關(guān)干擾頻率成分過多，并且特征頻率的倍頻均無(wú)法識(shí)別。整體來看，傳統(tǒng)包絡(luò)譜分析方法診斷效果不佳。

根據(jù)ASSD算法實(shí)現(xiàn)流程，當(dāng)K為2、3時(shí)，ρm1均大于0.1，且ρm2均小于0.1，繼續(xù)分解；當(dāng)K=4時(shí)，ρm1為0.245 3，大于0.1，但ρm2為0.162 4，大于0.1，停止分解。因此，可確定最佳奇異譜分量個(gè)數(shù)為3。

本文算法外圈故障信號(hào)處理結(jié)果如圖7所示。外圈故障試驗(yàn)信號(hào)經(jīng)過ASSD處理得到3個(gè)分量。經(jīng)計(jì)算，得到各SSC分量的合成峭度值如圖7(a)所示，選擇合成峭度最大的SSC-2為最優(yōu)分量。對(duì)SSC-2分量進(jìn)行MOMEDA解卷積處理，解卷積信號(hào)時(shí)域波形如圖7(b)所示。其周期性沖擊成分得以明顯加強(qiáng)。進(jìn)一步執(zhí)行FWEO獲得解卷積信號(hào)的瞬時(shí)能量，通過頻域變換得到的解卷積信號(hào)瞬時(shí)能量譜，如圖7(c)所示。外圈故障特征頻率的基頻～9倍頻成分均清晰呈現(xiàn)，分析效果與圖6(b)相比顯著提升。

圖7 本文算法外圈故障信號(hào)處理結(jié)果

5.2 內(nèi)圈故障信號(hào)分析

內(nèi)圈故障試驗(yàn)信號(hào)時(shí)域波形與包絡(luò)譜如圖8所示。

圖8 內(nèi)圈故障試驗(yàn)信號(hào)時(shí)域波形與包絡(luò)譜

選取電機(jī)驅(qū)動(dòng)端處軸承點(diǎn)蝕直徑為0.533 4 mm的內(nèi)圈故障數(shù)據(jù)，電機(jī)轉(zhuǎn)速為1 772 r/min，轉(zhuǎn)頻fr=29.53 Hz，根據(jù)式(21)計(jì)算得到內(nèi)圈理論故障特征頻率fi=159.9 Hz。圖8所示的譜圖中僅可識(shí)別出內(nèi)圈故障特征頻率的基頻成分，其倍頻成分均無(wú)法辨識(shí)，分析效果存在很大提升空間。

利用ASSD算法對(duì)該信號(hào)進(jìn)行處理。當(dāng)K為2～4時(shí)，ρm1均大于0.1，且ρm2均小于0.1，因此繼續(xù)分解；當(dāng)K=5時(shí)，ρm1為0.109 6，大于0.1，但ρm2為0.269 0，大于0.1，停止分解。因此，可確定最佳奇異譜分量個(gè)數(shù)為4。

本文算法內(nèi)圈故障信號(hào)處理結(jié)果如圖9所示。

圖9 本文算法內(nèi)圈故障信號(hào)處理結(jié)果

內(nèi)圈故障試驗(yàn)信號(hào)經(jīng)過分解處理得到4個(gè)SSC分量，各SSC分量的合成峭度值如圖9(a)所示。其中，SSC-3合成峭度值最大。因此，將SSC-3作為最優(yōu)分量作進(jìn)一步解卷積處理，所得解卷積信號(hào)時(shí)域波形如圖9(b)所示，最終所得解卷積信號(hào)瞬時(shí)能量譜如圖9(c)所示。其中，內(nèi)圈故障特征頻率基頻及其倍頻成分均被有效提取出來，據(jù)此可實(shí)現(xiàn)軸承狀態(tài)準(zhǔn)確判定。

6 結(jié)論

針對(duì)軸承微弱損傷難以準(zhǔn)確識(shí)別的問題，本文提出了基于ASSD-MOMEDA-FWEO的滾動(dòng)軸承故障診斷方法，并通過仿真及試驗(yàn)信號(hào)進(jìn)行驗(yàn)證，總結(jié)如下。

①采用合成峭度及斯皮爾曼等級(jí)相關(guān)系數(shù)，可自適應(yīng)確定SSD的奇異譜分量個(gè)數(shù)，實(shí)現(xiàn)原始信號(hào)的合理分解。通過合成峭度指標(biāo)可直接從ASSD分解結(jié)果中篩選出包含豐富故障信息的最優(yōu)分量，有利于后續(xù)分析處理進(jìn)程。

②MOMEDA算法能夠有效增強(qiáng)被噪聲淹沒的周期性沖擊成分。將其與FWEO算法進(jìn)行融合，作為后處理策略，可在強(qiáng)背景噪聲干擾下進(jìn)一步強(qiáng)化、放大微弱故障特征，有助于實(shí)現(xiàn)軸承故障的精確診斷。