付云驍，賈利民，楊 杰，魏秀琨，秦 勇

(1.中車工業(yè)研究院有限公司，北京 100070；2.北京交通大學軌道交通控制與安全國家重點實驗室，北京 100044；3.北京交通大學北京市城市交通信息智能感知與服務工程技術(shù)研究中心，北京 100044)

無論是既有線列車、城軌列車還是高速列車，作為列車轉(zhuǎn)向架的重要組成部分，軸箱軸承的故障診斷技術(shù)一直是被長期重點研究的列車安全保障關(guān)鍵技術(shù)之一。然而由于列車在途運行中輪對轉(zhuǎn)速不恒定，因而軸箱軸承的轉(zhuǎn)速狀態(tài)也隨之不斷改變，變速工況給軸承故障診斷的實現(xiàn)增加了困難。目前在旋轉(zhuǎn)機械的故障診斷技術(shù)中，振動信息的采集是最常采用的手段之一[1]，該方法在列車轉(zhuǎn)向架故障診斷測試中已被采用[2]。目前，滾動軸承故障診斷技術(shù)包括特征提取和智能診斷兩個重要環(huán)節(jié)。智能故障識別算法的診斷結(jié)果容易出現(xiàn)過分類和欠分類問題。若能夠通過低維特征對軸承的健康狀態(tài)進行快速、直觀的觀測，則會降低智能算法自身的誤分類概率，同時也能提升故障診斷技術(shù)的人機交互能力，降低健康檢查成本，但這要求特征參量對故障具備更強的敏感性。

非線性的滾動軸承振動信號的時域和頻域的細節(jié)信息可由自適應時頻域分析方法得到[3]。傳統(tǒng)的經(jīng)驗模態(tài)分解EMD(Empirical Mode Decomposition)[4]是通過插值迭代的思想從非線性信號提取本征模態(tài)函數(shù)IMF(Intrinsic Mode Function)，然而EMD分解所得IMF通常具有明顯的模態(tài)混疊效應。后續(xù)關(guān)于抑制模態(tài)混疊效應的一系列EMD改進方法，由于缺乏理論論證因此未能得到普適化的應用推廣[5]。聚合經(jīng)驗模態(tài)分解EEMD(Ensemble Empirical Mode Decomposition)[6]通過添加白噪聲以輔助EMD算法消除間歇現(xiàn)象，盡管能達到良好的改善效果，然而算法推演的時間復雜度及空間復雜度均明顯增加，使得計算的實時性明顯降低。文獻[7]在總結(jié)EMD和EEMD計算復雜度的基礎上，結(jié)合理論分析與實驗，總結(jié)出增強EEMD算法實時性的快速經(jīng)驗模態(tài)分解FEEMD。該方法可明顯提升軸箱軸承故障診斷的信號分解速率[8]。

滾動軸承的故障特征在時域、頻域以及時頻域3種情形都有體現(xiàn)[9]。由于旋轉(zhuǎn)軸承振動信號的非線性和非平穩(wěn)性，提取軸承時域和頻域的故障特征進行故障診斷相對困難[10]。而采用時頻分量結(jié)合信息理論推導軸承故障統(tǒng)計特征，如能量矩[11]、近似熵[12]等則比單純的時域或頻域特征更能全面反映故障狀態(tài)與正常狀態(tài)的區(qū)別。然而振動能量成分較高的IMF分量，其時域幅值分布更為相似，以IMF分量與原信號的時域分布相關(guān)測度為特征，可以避免非平穩(wěn)運行工況對故障診斷精度的干擾。相關(guān)熵[13]是對變量的一種泛化相關(guān)性度量，對隨機分布的變量具有魯棒性相關(guān)性度量性質(zhì)，變量間的局部微弱差異均可通過相關(guān)熵的高靈敏特性映射出來。另一方面，皮爾遜積矩相關(guān)系數(shù)(以下簡稱相關(guān)系數(shù))是對信號線性相關(guān)性的反映，是最常見的全局相關(guān)性的度量。鑒于此，這里將相關(guān)系數(shù)作為相關(guān)熵的幅度調(diào)制系數(shù)，計算振動信號與其IMF分量的相關(guān)熵函數(shù)，得到初步特征相關(guān)熵矩陣CM[14]。主元分析法PCA(Principal Component Analysis)的運算思想是將原始數(shù)據(jù)空間進行線性變換，進而得到各維度線性無關(guān)的數(shù)據(jù)空間，可在該數(shù)據(jù)空間進行特征變量的篩分，因而PCA在高維特征的優(yōu)化提取中較為常用[15]。若提取的特征維度小于三維，則可通過PCA實現(xiàn)滾動軸承故障特征的可視化。

本文結(jié)合FEEMD以及相關(guān)熵推理，提出一種新的滾動軸承故障特征——融合相關(guān)熵矩陣ICM，并利用PCA對相關(guān)熵矩陣進行特征變換及優(yōu)選。通過計算與模擬實驗，驗證了基于PCA的ICM對混合工況下滾動軸承故障辨識的魯棒能力，以及多維故障特征的降維可視化效果。

1 特征提取方法

融合相關(guān)熵矩陣的計算包括如下3個關(guān)鍵流程：

(1)IMF分量獲取。通過信號的非線性分解進行計算。

(2)獲取CM。分別計算各原始信號樣本與自身IMF分量的線性相關(guān)系數(shù)及相關(guān)熵，再將兩者相乘并整合成多維特征的樣本矩陣CM。

(3)獲取ICM。利用PCA對CM進行空間變換，并以方差貢獻率對變換后的特征矩陣進行篩分，得到ICM。

1.1 快速經(jīng)驗模態(tài)分解

EMD計算過程見文獻[4]。不做處理的EMD方法在對非線性信號進行信號分解時會出現(xiàn)間歇現(xiàn)象，導致模態(tài)混疊。為了使信號分解的結(jié)果更加準確，EEMD在EMD的基礎上采取添加高斯白噪聲的方法以消除模態(tài)混疊效應。在每次進行信號分解前，給原始信號加上不同幅值的高斯白噪聲。每次添加噪聲后的信號均需進行EMD分解，其多個EMD集成的過程，明顯增加了EEMD的計算復雜度。

EEMD過程如下：設初始信號為S(t)，對混有白噪聲Ni(t)的信號Si(t)進行EMD分解，得到若干IMF分量Cij(t)與一個余項ri(t)，其中Cij(t)表示第i次加入高斯白噪聲后，分解得到的第j個IMF分量。每次分解前加入均方根相等的不同白噪聲，并重復T次上述過程。計算IMF分量的總體均值，以消除高斯白噪聲的干擾效應，EEMD分解過程可描述為

( 1 )

式中：Cj(t)為整個EEMD分解后所得到的第j個IMF分量。

而FEEMD是在EEMD方法基礎上，通過簡化篩選準則，使計算效率有較大提升。首先，F(xiàn)EEMD給出固定常數(shù)作為硬閾值來約束IMF篩分過程，然后FEEMD以log2n作為一次EMD過程的停止篩分約束。此外，通過優(yōu)化計算EMD中的極值計算及三次樣條插值過程，也能明顯提升EEMD的計算效率。FEEMD的具體運算分析見文獻[7，16]。

令S(t)的長度為m，共k個故障類型，每個故障類型下采集l個S(t)，得到三維k×l×m的樣本矩陣SM。對SM進行FEEMD后，可得四維矩陣IMFM：Ck×l×n×m。

1.2 相關(guān)熵矩陣

隨機過程的相似性度量只能包含統(tǒng)計分布或時域結(jié)構(gòu)的信息，而不能包含兩者全部信息。相關(guān)熵(Correntropy)是通過Mercer核函數(shù)將數(shù)據(jù)從非線性空間投影到高維再生核希爾伯特空間，并計算高維空間的點的距離范數(shù)，得到高維空間數(shù)據(jù)線性相關(guān)的統(tǒng)計分布信息。Gauss核是典型的Mercer核，采用以核寬度為σ的Gauss核能夠很好地衡量隨機信號的相似度[14]。設m維向量X與Y的聯(lián)合概率為PXY，可得二維離散變量的相關(guān)熵估計為

( 2 )

由式( 2 )可知，計算每個IMF分量與原始信號的相關(guān)熵，可以獲得兩者在希爾伯特空間的細節(jié)相似度量。物理隨機過程的度量可以分析信號的整體趨勢特征，最直接的方法是線性相關(guān)系數(shù)法LCC(Linear Correlation Coefficient)，LCC不會隨變量尺度及順序的變化而改變。用LCC作為調(diào)幅系數(shù)可以穩(wěn)定相關(guān)熵的局部敏感變化，更平穩(wěn)地展現(xiàn)統(tǒng)計分布信息。若IMF分量Ci(t)與原始振動信號S(t)的LCC用γ表示，將γ帶入式( 1 )，長度為m的C(t)和S(t)的相關(guān)熵向量可以表示為

Φ·(Cn×m-Sl×m)

( 3 )

為便于計算管理，需對數(shù)值較大的CM進行歸一化處理。

( 4 )

將式( 3 )和式( 4 )推廣到IMFM，令Δ=k×l，Ω=Δ×n，則Ω維的相關(guān)熵矩陣(CM)為

( 5 )

式中：IMFΔ、SΔ表示二維k×l的樣本矩陣。為了便于進一步計算的表述，需要對φΩ部分轉(zhuǎn)置，令Γ=(k×l)′×n，則φΓ為所求CM。以矩陣表達的CM可以簡化計算過程，然而相關(guān)熵向量維度較高，仍有待進一步簡化處理。

1.3 融合特征

故障特征中冗余分量會降低故障識別率。為精簡故障特征，提高故障診斷精度，這里對φΓ做進一步的特征降維。因為k與l兩個參量可保持不變，因此特征降維需通過減小n來實現(xiàn)。這里采用主成分分析方法PCA對n進行處理。PCA計算過程保持樣本總方差不變，依統(tǒng)計方差降序排列統(tǒng)計特征，且樣本間始終保持相對獨立性。PCA處理后的數(shù)據(jù)空間，以特征向量方差貢獻率大小為依據(jù)，對特征值進行降序排列。若降維到可視空間，特征空間的組成應為3個以內(nèi)的主成分分量，由此得到希爾伯特空間的特征分布。

令p∈{1，2，3}，相關(guān)熵矩陣可通過PCA處理后，得到融合相關(guān)熵矩陣ICM。

( 6 )

這里的φl×k×p是滾動軸承故障辨識的特征全集。

經(jīng)PCA處理后ICM的類內(nèi)散度及類間散度均為非奇異矩陣[14]。為通過特征空間的距離測度實現(xiàn)故障可分性提供了理論支撐[17]。類內(nèi)及類間散度為

( 7 )

可由Sw、Sb構(gòu)造如下可分性判據(jù)[18]

( 8 )

式中：tr(·)為矩陣的跡。J1、J2和J3在任何非奇異變換下保持不變，J1，J2，J3越大，特征的類內(nèi)聚類性及類間可分性越好，反之聚類性及可分性越差。

2 城軌列車輪對軸承在途模擬試驗

為了檢驗ICM的可視化故障辨識能力，同時檢驗其對軸承轉(zhuǎn)速變化的抗擾性，滾動軸承故障辨識試驗模擬了城軌列車輪對軸承的不同轉(zhuǎn)速(勻速運行、加速啟動、減速制動)工況。并通過城軌列車輪對轉(zhuǎn)速模擬，驗證該方法對在途軸承故障診斷技術(shù)實現(xiàn)的有效性。

2.1 準備工作

試驗使用軸承型號為SKFN205M，其形狀為外圈固定的圓柱滾子軸承。具體參數(shù)見表1。

表1 軸承尺寸參數(shù)

試驗所用軸承采用外圈固定內(nèi)圈滾動的方式，利用放電加工技術(shù)分別在正常軸承各自的內(nèi)圈、外圈及滾子處理出0.3 mm的故障深度。外圈故障位于外圈垂向直徑的最下方。圖1為3種軸承故障模式。

圖1 內(nèi)圈故障、外圈故障及滾子故障示意圖

圖2是軸承振動試驗臺，可用于模擬城軌列車輪對運行過程。試驗平臺的電機額定功率為0.55 kW；最大轉(zhuǎn)速為1 450 r/min。試驗臺采用的速度控制方法為交流變頻器單相輸入三相輸出，其變速范圍為75～1 450 r/min，速度波動量在±2%以內(nèi)，可通過調(diào)節(jié)變頻器調(diào)整電機轉(zhuǎn)速。

圖2 滾動軸承振動試驗平臺

通過模擬城軌列車輪對運行工況，驗證模擬工況下的輪對軸承狀態(tài)辨識效果?？赏ㄟ^式( 9 )實現(xiàn)地鐵運行速度與軸承轉(zhuǎn)速的轉(zhuǎn)換。

( 9 )

式中：n為軸承轉(zhuǎn)速；v為列車運行速度；d為輪對直徑。

通過文獻[19]廣州地鐵列車的軸箱轉(zhuǎn)速工況，來設計試驗臺的軸承轉(zhuǎn)速，見表2。

2.2 試驗流程

測試勻加速、勻速及勻減速三組不同工況下的軸承辨識效果。除表2所示轉(zhuǎn)速差別外，其他運行參數(shù)均保持一致(為保證模擬實驗的合理性，選擇電機滿載運行，載荷恒定為1 kW)。試驗過程如圖3所示。

步驟1樣本矩陣生成。振動加速度傳感器采樣頻率為12 kHz，每段截取的樣本長度均為1 024。每組工況中，每個故障狀態(tài)下均提取160個樣本段，各工況的全樣本矩陣SM維度均為k×l×m=4×640×1 024。

表2 滾動軸承試驗參數(shù)

圖3 故障特征提取流程圖

步驟2對SM1、SM2和SM3分別進行FEEMD，得到k×l×n×m=4×640×10×1 024的IMFM1、IMFM2和IMFM3。各工況下分別對各故障類別取一段樣本，圖4(a)所示為各樣本的故障信號時域波形圖，圖4(b)～圖4(d)所示為內(nèi)圈故障在加速、勻速和減速時的IMF分量示意圖。

步驟3構(gòu)建融合特征矩陣。IMFM1、IMFM2和IMFM3通過相關(guān)熵運算分別得到矩陣CM1、CM2和CM3，三者具有同樣的維度k×l×n=4×640×6。根據(jù)經(jīng)驗CM的核尺度參數(shù)σ設為0.5。再通過PCA得到維度為k×l×n=4×640×p(p∈{1，2，3})的融合相關(guān)熵矩陣ICM1、ICM2和ICM3。

圖4 原信號及內(nèi)圈故障樣本IMF分量

以上為計算ICM的具體過程。為了驗證ICM的有效性，依據(jù)圖3所示流程，將相關(guān)熵算子改為能量矩算子、譜峭度算子及樣本熵算子，分別對傳統(tǒng)優(yōu)勢特征能量矩[11]、譜峭度[9]及樣本熵[12]進行運算，各自得到融合能量矩IEM(Integrated energy moment)、融合譜峭度ISK(Integrated spectral kurtosis)和融合樣本熵矩陣ISM(Integrated sample-entropy matrix)。分別提取勻加速、勻速及勻減速3個運行狀態(tài)下的IEM、ISK及ISM。這里分別提取二維融合的ICM、IEM、ISK及ISM特征，在二維特征下，PCA的累積方差貢獻率已經(jīng)達到90%。

2.3 結(jié)果分析

從圖5可以得出結(jié)論：ICM的特征分布中，滾動軸承的特征在正常與故障狀態(tài)下分別呈現(xiàn)出明顯的可分性。不同故障狀態(tài)(內(nèi)圈故障、外圈故障、滾子故障)下的樣本分布也各自形成一定的分布區(qū)域。而作為對比的三組特征分布混亂，難以辨識相應的故障狀態(tài)。尤其IEM、ISK兩組特征分布不具備類間可分度，因此這三組故障特征明顯不可辨識正常軸承與故障軸承。故測試結(jié)果表明，在勻加速工況下，二維ICM可用于故障特征辨識軸承的正常與故障狀態(tài)，二維IEM、ISK和ISM不具備故障可辨識性，因而不可作為故障特征。

圖5 模擬列車加速階段的融合特征分布

觀察圖6可知，勻速運行階段的ICM具備正常狀態(tài)和故障狀態(tài)的明顯可分性。在故障狀態(tài)下，內(nèi)圈故障相對容易辨識，而外圈故障和滾子故障下的樣本可分性相對較弱。反觀另外三組對照特征的樣本分布，無論是IEM、ISK還是ISM，均不能從樣本狀態(tài)分布中辨識出正常樣本區(qū)域。相比較四類融合特征，ICM的樣本特征同樣具有顯著聚類性和故障可分性。因此可知，二維ICM可用于勻速工況下的滾動軸承故障檢測，而二維IEM、ISK和ISM不具備故障可辨識度，不可用于勻速工況下的軸承故障檢測。

觀察圖7可知，模擬減速階段的ICM特征分布同樣具備故障可辨識性。減速階段的正常狀態(tài)、內(nèi)圈故障、外圈故障以及滾子故障4種狀態(tài)均具備較好的聚類性，且類間間距明顯。減速階段的IEM在特征分布中可以看出一定的故障可辨識性，而正常狀態(tài)和內(nèi)圈故障狀態(tài)之間的可辨識性較差，相較于安全與故障的高識別率要求，IEM并不能很好地體現(xiàn)這一點。而另外兩組對照特征ISK和ISM則明顯不具備故障可分能力。因此在減速過程，ICM和IEM可作為故障特征對滾動軸承進行故障診斷；而IEM的安全辨識能力弱，ICM相較而言仍然是辨識效果最好的融合特征。

圖6 模擬列車勻速階段的融合特征分布

圖7 模擬列車減速階段的融合特征分布

以上是觀察特征空間分布圖得出的定性結(jié)論。以二維融合特征的空間分布參量作為定量指標，分析融合特征是否具備故障可分性。

ICM、IEM、ISK和ISM故障辨識能力的定量判定J1、J2和J3由式( 8 )計算得到，具體計算結(jié)果見表3。J1、J2、J3越大，特征的故障辨識能力越強。通過分析表3所示的數(shù)值，可得出如下結(jié)論：除勻加速時ICM的J1小于IEM的J1外，在勻加速階段、勻速階段、勻減速階段中，除勻加速時ICM的J1小于IEM的J1外，其他情況下ICM的J1、J2和J3均最高。這表明ICM的故障辨識能力優(yōu)于另外三組對照特征。而IEM在勻加速和勻減速階段的定量判據(jù)均大于ISK和ISM；在勻速階段ISM的定量判據(jù)大于IEM和ISK，定量分析的結(jié)論與定性分析的結(jié)論相吻合。

上述從定性和定量兩個方面綜合判斷了ICM具備精確的故障可分能力，是在勻變速及勻速環(huán)節(jié)均可用于軸承故障辨識的特征依據(jù)。

表3 故障特征有效性判據(jù)

3 結(jié)束語

本文提出基于FEEMD-ICM-PCA的滾動軸承故障魯棒辨識方法，驗證了該方法有效隔離工況變化對軸承故障辨識的影響。為了使計算及試驗的分析更加簡便，采用矩陣對軸承采集的樣本分塊處理。提取滾動軸承振動信號的SM進行FEEMD分解，得到IMFM，進一步通過相關(guān)熵算子計算，得到CM，最后用PCA對CM進行空間變換，提取前兩維特征，得到融合相關(guān)熵分量，即ICM。ICM的樣本空間分布，可以直觀地辨識軸承的故障類型。

通過對比勻加速、勻速及勻減速三組模擬運行階段的滾動軸承振動融合特征，從定性和定量的角度綜合論證了二維ICM適合于滾動軸承故障診斷的結(jié)論，且ICM對軸承轉(zhuǎn)速的勻速變化具有明顯的魯棒性。

為驗證現(xiàn)階段基于FEEMD-ICM-PCA的列車軸箱軸承可視化故障辨識方法的可行性，實現(xiàn)了單一故障辨識，同時驗證了在變轉(zhuǎn)速工況下進行故障辨識的可行性。這對處理滾動軸承的并發(fā)故障，以及安全預警具有重要的指導意義，但仍需要進一步試驗進行驗證。此外，基于試驗臺模擬的列車軸箱軸承的故障辨識試驗，充分表明ICM作為軸箱軸承故障特征的可視化及魯棒性優(yōu)勢，為下一步的實車試驗提供了可行的參考依據(jù)。