李志寧，孫宜權(quán)，2 ，張英堂，田 昊，范紅波

（1.軍械工程學(xué)院，河北 石家莊 050003；2.66267部隊(duì)，河北 石家莊 050081）

發(fā)動(dòng)機(jī)的軸承故障影響著發(fā)動(dòng)機(jī)的經(jīng)濟(jì)性、動(dòng)力性及行駛安全性。由于發(fā)動(dòng)機(jī)結(jié)構(gòu)復(fù)雜，軸承故障特征微弱，在其他信號(hào)干擾下較難被檢測(cè)到，如何發(fā)現(xiàn)軸承故障并作出診斷，對(duì)于減小經(jīng)濟(jì)損失和排除安全隱患具有重要的意義。目前，在發(fā)動(dòng)機(jī)軸承故障診斷中，已有較多的信號(hào)分析方法應(yīng)用于單一故障診斷，且效果較好［1］，而對(duì)于并發(fā)故障，由于故障信號(hào)之間的耦合與調(diào)制，導(dǎo)致診斷的難度增加、準(zhǔn)確率下降［2］。

盲源分離（BSS）算法為發(fā)動(dòng)機(jī)軸承并發(fā)故障信號(hào)分離提供了一種有效的方法，但是，目前在發(fā)動(dòng)機(jī)振動(dòng)信號(hào)盲源分離中，通常假設(shè)觀測(cè)信號(hào)的數(shù)目不小于源信號(hào)的數(shù)目，而在實(shí)際應(yīng)用中，這個(gè)假設(shè)往往是不成立的。對(duì)于欠定盲源分離，已有的研究基本上是基于源信號(hào)的稀疏表示，借助聚類(lèi)算法實(shí)現(xiàn)。由于柴油機(jī)本身的復(fù)雜性及其工作環(huán)境的多樣性，從柴油機(jī)上獲取的振動(dòng)信號(hào)必然是多個(gè)振源信號(hào)以未知的混合形式混疊而成的復(fù)雜信號(hào)，所以對(duì)于柴油機(jī)振動(dòng)信號(hào)，源信號(hào)的稀疏性假設(shè)并不一定成立［3］。

本研究提出了基于平行因子分析（PARAFAC）的發(fā)動(dòng)機(jī)軸承并發(fā)故障信號(hào)盲分離方法，利用求得的信號(hào)傳播路徑響應(yīng)函數(shù)構(gòu)造了虛擬觀測(cè)信號(hào)，將盲源分離從欠定轉(zhuǎn)變?yōu)檫m定或超定，然后采用自適應(yīng)PARAFAC方法解決柴油機(jī)軸承并發(fā)故障信號(hào)相互耦合的故障模式。

1 平行因子分析方法

平行因子分析概念始于心理實(shí)驗(yàn)學(xué)領(lǐng)域的數(shù)據(jù)分析，是一種多維低秩分解方法。Caroll等人發(fā)展了平行因子分析模型，用于無(wú)線通信和盲源分離等領(lǐng)域［4］。本研究以3維數(shù)據(jù)為例，數(shù)據(jù)X∈CI×J×K集可以分解為M 個(gè)秩為1的3維矩陣之和，如式（1）所示［5］：

式中：am∈CI×L；bm∈CJ×L；cm∈CK×L；表示張量積。X中的元素可表示為

若A＝［a1，…aM］∈CI×M，B＝［b1，…bM］∈CJ×M，C＝［c1，…cM］∈CK×M，稱矩陣A，B，C為X 的模式矩陣。為了討論X的分解是否唯一，引入定義1和定理1。

定義1［6］：若矩陣A的任意kA個(gè)列線性獨(dú)立，則最大的kA值稱為矩陣A的Kruskal秩，簡(jiǎn)稱k秩。

定理1［7］：A，B，C 的定義見(jiàn)式（2）。存在 K 個(gè)矩陣Xk＝ADk（C）BT，k＝1，2，3，…K，Dk（C）是對(duì)矩陣C的第k行進(jìn)行對(duì)角化處理；若kA＋kB＋kC≥2（M＋1），則分解得到的矩陣估計(jì)值A(chǔ)，BT，C各列排列次序與比例因子不確定外，可唯一確定。

2 欠定盲源分離算法

2.1 信號(hào)傳播路徑響應(yīng)函數(shù)預(yù)測(cè)

本研究利用PARAFAC方法，預(yù)測(cè)傳播路徑響應(yīng)函數(shù)，從而消除延遲濾波對(duì)觀測(cè)信號(hào)的影響。卷積混合信號(hào)的數(shù)學(xué)描述見(jiàn)式（3）［8］：

式中：S（t）＝［s1（t），s2（t），…sI（t）］T為I個(gè)源信號(hào)構(gòu)成的向量；X（t）＝［x1（t），x2（t），…xJ（t）］T為J個(gè)混合信號(hào)構(gòu)成的向量；混合矩陣H為J×L維矩陣；＊為卷積運(yùn)算；L為傳輸路徑響應(yīng)長(zhǎng)度。

為預(yù)測(cè)傳播路徑響應(yīng)函數(shù)，將觀測(cè)信號(hào)變換到頻域分析。式（3）經(jīng)離散傅里葉變換后為［8］

式中：頻率標(biāo)志f＝1，…F；p為數(shù)據(jù)段標(biāo)志；X（f，p）＝ ［x1（f，p），…xJ（f，p）］T；S（f，p）＝［s1（f，p），…sI（f，p）］T；H（f）的第i列代表第i 個(gè)源信號(hào)在頻率f處的傳播路徑響應(yīng)。

假設(shè)每一個(gè)觀測(cè)信號(hào)包含N個(gè)采樣點(diǎn)，將其分割成P個(gè)不交疊的數(shù)據(jù)段，每個(gè)數(shù)據(jù)段包含NP＝ent（N/P）個(gè)數(shù)據(jù)點(diǎn)，用p＝1，…P 標(biāo)志數(shù)據(jù)段，當(dāng)采樣頻率為Fs時(shí)，每個(gè)數(shù)據(jù)段持續(xù)的時(shí)間Tp＝（Np/Fs）。在此基礎(chǔ)上，J×J維自相關(guān)矩陣Rx（f，p）＝E［X（f，p）XH（f，p）］可寫(xiě)為

式中：Rs（f，p）＝E［S（f，p）SH（f，p）］為p 數(shù)據(jù)段，頻率為f的源信號(hào)自相關(guān)矩陣。式（5）展開(kāi)可寫(xiě)為

記Rx（f，p）為γ（x）j1，j2，p（f），其中，j1＝1，…J，j2＝1，…J，p＝1，…P；H（f）中的元素記為hj，i（f）；構(gòu)建一個(gè)P×I維矩陣C（f），其中的元素記為cp，i（f），代表Rs（f，p）矩陣對(duì)角線上第i個(gè)元素。根據(jù)式（6）γ（x）j1，j2，p（f）可寫(xiě)為

與式（2）對(duì)比可知，式（7）為PARAFAC分解，其中，I為Rs（f，p）的秩。根據(jù)PARAFAC方法，最終可得唯一解｛H（f），f＝1，…F｝，但各列解的排列次序與比例因子不確定。通過(guò)對(duì)排列次序與比例因子進(jìn)行自適應(yīng)調(diào)整，從而解決上述不確定問(wèn)題，稱為自適應(yīng)PARAFAC算法，具體過(guò)程參考文獻(xiàn)［5］。

2.2 通道虛擬方法

考慮到觀測(cè)信號(hào)受傳播路徑的影響，圖1示出了源信號(hào)S（f，p）至觀測(cè)信號(hào)X（f，p）的傳播過(guò)程，其中路徑響應(yīng)函數(shù)H（f）中的元素簡(jiǎn)記為hj，i。

從圖1可知，當(dāng)路徑響應(yīng)函數(shù)各列解的排列次序與比例因子正確時(shí)，觀測(cè)信號(hào)為各通道傳感器拾取的信號(hào)；當(dāng)傳播路徑?jīng)_擊響應(yīng)函數(shù)各列解隨機(jī)排列次序時(shí)，觀測(cè)信號(hào)可看作虛擬通道傳感器拾取的信號(hào)，稱之為虛擬觀測(cè)信號(hào)。

2.3 源數(shù)估計(jì)

為實(shí)現(xiàn)信號(hào)欠定盲分離，首先要估計(jì)系統(tǒng)的源信號(hào)數(shù)。提出了虛擬觀測(cè)信號(hào)與主分量分析（PCA）結(jié)合的源數(shù)估計(jì)方法［9］。

利用通道虛擬方法得到虛擬觀測(cè)信號(hào)，將觀測(cè)信號(hào)與虛擬觀測(cè)信號(hào)組成新的多維觀測(cè)信號(hào)XN（t）＝［x1（t），…xJ（t），x1V（t），x2V（t），…］T，其中，［x1V（t），x2V（t），…］T為虛擬觀測(cè)信號(hào)，解決了源信號(hào)數(shù)大于觀測(cè)信號(hào)數(shù)的難題。

PCA可以確定信號(hào)和噪聲子空間的維數(shù)，假設(shè)將XN（t）∈Rm建模成

式中：H∈Rm×n是一個(gè)列滿秩混合矩陣，且m＞n；S（t）∈Rn是一個(gè)零均值高斯源，其協(xié)方差矩陣Rss＝E｛S（t）·SH（t）｝；V（t）∈Rm是高斯零均值噪聲，其協(xié)方差矩陣Rw＝σ2vIm。而且，S（t）與V（t）不相關(guān)。由上述假設(shè)可知，XN（t）的協(xié)方差矩陣可表示為

式中：H·Rss·HH＝Vs·Λs·VHs，Vs包含與n個(gè)主特征值Λs＝diag｛λ1≥λ2≥λ3≥…≥λn｝對(duì)應(yīng)的特征向量；VN包含（m－n）個(gè)對(duì)應(yīng)于噪聲特征值ΛN＝diag｛λn＋1，λn＋2，…λm｝＝σ2vIm－n的特 征 向 量。理 論上Rxx的（m－n）個(gè)最小特征值等于σ2v，因此，從最小特征值的重復(fù)個(gè)數(shù)得到信號(hào)子空間的維數(shù)。但是在實(shí)際應(yīng)用中，樣本數(shù)是有限的，得到的最小特征值一般是不同的。為能得到信號(hào)自空間的維數(shù)，引入Bayesian信息準(zhǔn)則［9］：

2.4 算法仿真研究

發(fā)動(dòng)機(jī)振動(dòng)特性比較復(fù)雜，無(wú)法準(zhǔn)確描述。通常情況下，可簡(jiǎn)單用多個(gè)單自由度彈簧－質(zhì)量－阻尼系統(tǒng)模型描述［8］。對(duì)上述系統(tǒng)物理模型進(jìn)行力學(xué)分析得到：

式中：M為系統(tǒng)質(zhì)量塊的質(zhì)量；K為系統(tǒng)的彈性系數(shù)；C為系統(tǒng)的阻尼系數(shù)；X（t）為系統(tǒng)產(chǎn)生的位移；F（t）為與發(fā)動(dòng)機(jī)轉(zhuǎn)速相關(guān)的激勵(lì)力。F（t）可表示為

式中：mu為轉(zhuǎn)軸不平衡質(zhì)量；ru為轉(zhuǎn)軸不平衡質(zhì)量距離轉(zhuǎn)軸軸心的距離；ω為轉(zhuǎn)軸的角速度。

本研究以振動(dòng)加速度的形式采集信號(hào)，通過(guò)改變式（12）中的ω（分別為100r/min，160r/min，120r/min）產(chǎn)生3個(gè)不同的激勵(lì)力，進(jìn)而得到3個(gè)源信號(hào)（圖2前3個(gè)），最后1個(gè)源信號(hào)（圖2最后1個(gè)）通過(guò)下式得到：

通過(guò)設(shè)置源信號(hào)與觀測(cè)信號(hào)間的傳播路徑?jīng)_擊響應(yīng)［5］，并利用卷積運(yùn)算得到3個(gè)觀測(cè)信號(hào)（圖3的前3個(gè)觀測(cè)信號(hào)），通過(guò)虛擬通道的方法，得到4個(gè)虛擬觀測(cè)信號(hào)（圖3的后4個(gè)觀測(cè)信號(hào)）。各觀測(cè)信號(hào)的特征值貢獻(xiàn)率見(jiàn)圖4。從中可以看出，前4個(gè)特征值貢獻(xiàn)率之和接近100%，后3個(gè)特征值貢獻(xiàn)率很小，進(jìn)一步由Bayesian信息準(zhǔn)則可以判斷，仿真信號(hào)的源數(shù)估計(jì)為4，與實(shí)際的源信號(hào)數(shù)一致。

以實(shí)際觀測(cè)信號(hào)與虛擬觀測(cè)信號(hào)組成新的觀測(cè)信號(hào)進(jìn)行盲源分離，使盲源分離從欠定變?yōu)檫m定或超定。圖5示出了分離信號(hào)，對(duì)比圖2可知，觀測(cè)信號(hào)得到較好的分離。

3 發(fā)動(dòng)機(jī)曲軸連桿軸承并發(fā)故障振動(dòng)信號(hào)盲分離

3.1 試驗(yàn)系統(tǒng)相關(guān)設(shè)置

試驗(yàn)在東風(fēng)EQ6BT 6缸四沖程發(fā)動(dòng)機(jī)上進(jìn)行，為得到反映曲軸連桿軸承并發(fā)故障振動(dòng)信號(hào)，減少其他缸的振動(dòng)干擾，根據(jù)文獻(xiàn)［10］可知，傳感器最佳安裝部位為發(fā)動(dòng)機(jī)缸體和油底殼結(jié)合處左右側(cè)及油底殼處。在這3個(gè)部位安裝加速度振動(dòng)傳感器（見(jiàn)圖6）。振動(dòng)信號(hào)在加速狀態(tài)下同步采集，采樣頻率為12.8kHz，采樣點(diǎn)數(shù)為16 384，采集時(shí)利用設(shè)定的轉(zhuǎn)速（1 800r/min）觸發(fā)采集。試驗(yàn)設(shè)置的曲軸軸承與連桿軸承技術(shù)狀態(tài)見(jiàn)表1［2］。本研究是在參考文獻(xiàn)［2］的基礎(chǔ)上對(duì)軸承并發(fā)故障診斷方法進(jìn)行更深入的研究。

表1 曲軸連桿軸承并發(fā)故障參數(shù)設(shè)置

3.2 軸承并發(fā)故障振動(dòng)信號(hào)盲分離

試驗(yàn)安裝的3個(gè)傳感器離軸承并發(fā)故障振源最近，測(cè)量得到的信號(hào)對(duì)軸承配合間隙變化敏感［11］。由于發(fā)動(dòng)機(jī)燃爆振動(dòng)信號(hào)的主要頻率為1kHz，3kHz，5kHz附近［12］，軸承故障信號(hào)特征頻率主要分布在3kHz以下低頻段［2］，為減弱燃爆振動(dòng)信號(hào)的干擾，對(duì)觀測(cè)信號(hào)進(jìn)行低通濾波，正常狀態(tài)下濾波后的觀測(cè)信號(hào)見(jiàn)圖7（前3個(gè)為實(shí)際觀測(cè)信號(hào)，后2個(gè)為虛擬觀測(cè)信號(hào)），源數(shù)估計(jì)為5。

圖8示出了該觀測(cè)信號(hào)的頻譜，圖9示出了正常狀態(tài)分離得到的時(shí)域信號(hào)，圖10示出了其分離信號(hào)的頻譜，由于篇幅有限，在此只列出其他技術(shù)狀態(tài)的頻譜圖，圖11至圖14分別對(duì)應(yīng)曲軸軸承輕微故障＋連桿軸承輕微故障頻譜、曲軸軸承輕微故障＋連桿軸承嚴(yán)重故障頻譜、曲軸軸承嚴(yán)重故障＋連桿軸承輕微故障頻譜、曲軸軸承嚴(yán)重故障＋連桿軸承嚴(yán)重故障頻譜。

從圖8可以看出，對(duì)于軸承并發(fā)故障振動(dòng)信號(hào)，即使觀測(cè)信號(hào)離振源很近，也無(wú)法從頻譜中獲取有效的特征信號(hào)，主要由于軸承并發(fā)故障振動(dòng)是微弱信號(hào)，淹沒(méi)在其他信號(hào)中。通過(guò)自適應(yīng)PARAFAC算法，得到分離信號(hào)，并得到它的頻譜分布。已有的研究表明，發(fā)動(dòng)機(jī)加速運(yùn)轉(zhuǎn)時(shí)，當(dāng)轉(zhuǎn)速為1 800r/min時(shí)觸發(fā)采集，曲軸軸承故障的特征頻率段在1 500Hz附近，連桿軸承故障的特征頻率段在1 000Hz附近［13］。由圖11可知，分離信號(hào)1與信號(hào)5的能量分別集中在1 414Hz和2 516Hz附近，其中分離信號(hào)1反映了曲軸軸承故障的特征頻率分布；由于2 516Hz為1 000Hz附近頻率和1 500Hz附近頻率的和頻，因此，分離信號(hào)5反映了連桿軸承和曲軸軸承故障的特征頻率分布。由圖12可知，分離信號(hào)2與信號(hào)4的能量分別集中在1 629Hz和2 274Hz附近，其中分離信號(hào)2反映了曲軸軸承故障的特征頻率分布；由于2 274Hz為1 000Hz附近頻率的倍頻，因此，分離信號(hào)4反映了連桿軸承故障的特征頻率分布。由圖13可知，分離信號(hào)2與信號(hào)5的能量分別集中在2 168Hz和625Hz附近，由于2 168Hz為1 000Hz附近頻率的倍頻，因此，分離信號(hào)2反映了連桿軸承故障的特征頻率分布；由于625Hz為1 500Hz附近頻率和1 000Hz附近頻率的差頻，分離信號(hào)5反映了連桿軸承和曲軸軸承故障的特征頻率分布。同理易知，圖14中，分離信號(hào)1與信號(hào)4反映了曲軸軸承和連桿軸承故障的特征頻率分布。曲軸軸承和連桿軸承的故障特征頻率來(lái)自經(jīng)驗(yàn)值，并沒(méi)有計(jì)算公式可參考。

由上述分析可知，雖然自適應(yīng)PARAFAC方法可消除并發(fā)故障信號(hào)間的耦合，但對(duì)于振源較多、故障模式復(fù)雜、耦合嚴(yán)重的信號(hào)盲分離，各工況下分離信號(hào)的頻譜分布和文獻(xiàn)［13］的特征頻率分布并非完全一致，分離信號(hào)中包含了故障特征頻率的倍頻、和頻、差頻，其余各諧波為其他頻率的非線性組合，根據(jù)以上征兆，可診斷發(fā)動(dòng)機(jī)軸承出現(xiàn)了并發(fā)故障。

4 結(jié)論

a）基于平行因子分析法，實(shí)現(xiàn)了信號(hào)傳播路徑響應(yīng)函數(shù)的預(yù)測(cè)，并以此構(gòu)造了虛擬觀測(cè)信號(hào)，使盲源分離從欠定變?yōu)檫m定或超定；

b）針對(duì)信號(hào)傳輸延時(shí)濾波導(dǎo)致盲源分離信號(hào)的難題，利用自適應(yīng)PARAFAC算法得到的分離信號(hào)準(zhǔn)確診斷出曲軸軸承和連桿軸承耦合故障，為診斷發(fā)動(dòng)機(jī)軸承并發(fā)耦合故障提供了一種新方法；

c）為診斷發(fā)動(dòng)機(jī)每缸所對(duì)應(yīng)的曲軸和連桿軸承的并發(fā)故障，需安裝較多的振動(dòng)加速度傳感器，并且提出的故障診斷方法需人工判斷，如何實(shí)現(xiàn)故障診斷的自動(dòng)化仍需進(jìn)一步研究。

［1］ 田 昊，唐力偉，田 廣.基于盲源分離的齒輪箱復(fù)合故障診斷研究［J］.兵工學(xué)報(bào)，2010，31（5）：646－649.

［2］ 梅檢民，張玲玲，肖云魁，等.基于高階累積量的軸承并發(fā)故障振動(dòng)信號(hào)分析［J］.內(nèi)燃機(jī)學(xué)報(bào)，2011，29（4）：327－331.

［3］ 毋文峰，陳小虎，蘇勛家.基于經(jīng)驗(yàn)?zāi)Ｊ椒纸獾膯瓮ǖ罊C(jī)械信號(hào)盲分離［J］.機(jī)械工程學(xué)報(bào)，2011，47（4）：13－15.

［4］ 梁軍利，冀邦杰，趙 峰，等.一種基于平行因子分析的近場(chǎng) 源 定 位 新 方 法 ［J］.電 子 學(xué) 報(bào)，2007，35（10）：1909－1914.

［5］ Nion D，Mokios K N，Sidiropoulos N D，et al.Batch and Adaptive PARAFAC－based Blind Separation of Convolutive Speech Mixtures［J］.IEEE Trans.on Audio，Speech and Language Processing，2010，18 （6）：1193－1207.

［6］ Sidropoulos N D，Giannakis G B.Parallel factor analysis in sensor array processing［J］.IEEE Trans Signal Processing，2000，48（8）：2377－2388.

［7］ Sidropoulos N D，Giannakis G B，BIOR.Blind PARAFAC receivers for DS－CDMA systems［J］.IEEE Trans Signal Processing，2000，48（3）：810－823.

［8］ LI W，Gu，BALL A D，et al.A study of the noise from diesel engines using the independent component analysis［J］.Mechanical System and Signal Processing，2001，15（6）：1165－1184.

［9］ 贠亞男，許明圣.未知信源個(gè)數(shù)的自適應(yīng)盲信號(hào)分離算法［J］.科學(xué)技術(shù)與工程，2011，32（11）：7931－7933.

［10］ 肖云魁，李世義，王建新，等.基于小波包——AR譜技術(shù)提取柴油發(fā)動(dòng)機(jī)曲軸軸承故障特征［J］.北京理工大學(xué)學(xué)報(bào)，2004，24（6）：508－511.

［11］ 肖云魁，曹亞娟，吳 曉，等.用循環(huán)譜理論分析發(fā)動(dòng)機(jī)曲軸軸承加速振動(dòng)信號(hào)［J］.振動(dòng)、測(cè)試與診斷，2008，28（2）：117－121.

［12］ 史玉鵬，劉建敏，許世永，等.基于缸蓋振動(dòng)的柴油機(jī)噴油、燃燒信息檢測(cè)研究［J］.車(chē)用發(fā)動(dòng)機(jī)，2011（3）：74－76.

［13］ 趙慧敏.發(fā)動(dòng)機(jī)非穩(wěn)態(tài)振動(dòng)信號(hào)分析與智能故障診斷研究［D］.天津：天津大學(xué)，2010.