郭明軍 李偉光? 楊期江 趙學(xué)智

(1.華南理工大學(xué) 機(jī)械與汽車工程學(xué)院，廣東 廣州 510640；2.廣州航海學(xué)院 輪機(jī)工程學(xué)院，廣東 廣州 510725)

主成分分析[1](Principal Component Analysis，PCA)是一種基于二階統(tǒng)計的數(shù)據(jù)分析方法，通過分析各變量之間的相關(guān)關(guān)系，利用一組較少的、互不相關(guān)的新變量(即主元)來線性表示初始變量。通過PCA方法可以消除隨機(jī)變量之間的相關(guān)性、剔除冗余信息，從而達(dá)到突出原始數(shù)據(jù)中隱含特性的目的[2]。目前，PCA在數(shù)字水印[2]、圖像處理[3]、模式識別[4]、特征提取[5- 6]及故障診斷[7- 8]等領(lǐng)域都得以廣泛應(yīng)用。Li等[5]將PCA方法應(yīng)用于工業(yè)齒輪箱的狀態(tài)監(jiān)測和故障診斷，獲得很好的降維和分類效果。Wang等[6]采用PCA提取局部放電故障(PD)信號的多尺度分散熵(MDE)中的有效主成分作為輸入特征，然后采用超球面多類支持向量機(jī)(HMSVM)進(jìn)行PD模式識別，實現(xiàn)電氣設(shè)備絕緣狀態(tài)診斷。Gharavian等[7]將PCA方法應(yīng)用于汽車變速器故障診斷。Gao等[8]提出一種融合了EEMD、PCA和PNN 3種方法的改進(jìn)滾動軸承的故障診斷方法，其中PCA用于故障特征的降維和去除相關(guān)性。

確定有效主元的個數(shù)是主成分分析的關(guān)鍵問題，常規(guī)方法是根據(jù)累計貢獻(xiàn)率[5- 8]取定某個閾值的方法選擇有效主元的個數(shù)，其表達(dá)式為

(1)

式中，Li表示累計貢獻(xiàn)率，λi為協(xié)方差矩陣的特征值，δ為某個選定的閾值。由此可見，閾值越大主成分的個數(shù)就越多，保留的信息量就越大，從而可能保留的噪聲成分也越多。聶振國等[9]提出一種改進(jìn)的PCA算法，并應(yīng)用于滾動軸承信號的提純，取得了較好的效果。李振等[10]研究表明，信號中的每個頻率對應(yīng)兩個有效特征值，其幅值決定協(xié)方差矩陣的特征值在其分布圖中出現(xiàn)的先后順序，并提出一種基于PCA的算法用于軸心軌跡的提純。

然而，上述應(yīng)用當(dāng)中，都是基于特征值分解協(xié)方差矩陣實現(xiàn)PCA算法(姑且稱為傳統(tǒng)PCA算法)，這種方法計算量較大且存在一定的舍入誤差。文獻(xiàn)[9]研究結(jié)果表明，奇異值分解(Singular value decomposition，SVD)因無需計算協(xié)方差矩陣，可以避免舍入誤差，且信號的重構(gòu)誤差較傳統(tǒng)PCA算法小。鑒于此，文中提出一種全新的思路：基于SVD分解矩陣的原理實現(xiàn)PCA算法。首先，通過理論推導(dǎo)得到協(xié)方差矩陣特征值與矩陣奇異值之間的代數(shù)關(guān)系及協(xié)方差矩陣特征向量和矩陣奇異向量之間的對應(yīng)關(guān)系；然后，利用有效特征值與頻率分量個數(shù)的關(guān)系，提出一種基于SVD原理的PCA特征頻率提取算法，并通過仿真信號驗證算法的有效性；最后，將該算法應(yīng)用于大型滑動軸承試驗臺主軸的軸心軌跡提純中。

1 Hankel矩陣方式下PCA信號處理原理

1.1 PCA信號處理原理

主成分分析[11]是指通過正交變換將一組相關(guān)變量xi(i=1，2，…，m)轉(zhuǎn)換為一組線性不相關(guān)的變量yj(j=1，2，…，k)，即有：

(2)

式中，系數(shù)αi=(αi1，αi2，…，αim)T(i=1，2，…，k)為協(xié)方差矩陣C中降序排列的第i個特征值λi對應(yīng)的特征向量，且αi滿足：

(3)

協(xié)方差矩陣C的計算式為

(4)

式中，cov(xi，xj)=E[(xi-E(xi))(xj-E(xj))T]。

協(xié)方差矩陣C的特征方程為

Cαi=λiαi

(5)

將式(2)改寫為分量形式有：

(6)

其中：yj∈R；j=1，2，…，k。

1.2 Hankel矩陣的構(gòu)造

在式(4)中，矩陣A采用Hankel矩陣。對于一維零均值的離散信號a=[a(1)a(2) …a(N)]。利用此信號按照以下方法構(gòu)造Hankel矩陣：

(7)

矩陣A稱為重構(gòu)吸引子軌跡矩陣，也稱為Hankel矩陣。針對Hankel矩陣的構(gòu)造問題，趙學(xué)智等[12]證明了噪聲去除量與列數(shù)呈拋物線的對稱關(guān)系，并得出結(jié)論：如果信號長度N為偶數(shù)，應(yīng)取列數(shù)n=N/2、行數(shù)m=N/2+1來構(gòu)造Hankel矩陣；如果N為奇數(shù)，應(yīng)取列數(shù)n=(N+1)/2、行數(shù)m=(N+1)/2來構(gòu)造Hankel矩陣。令A(yù)=[x1x2…xm]T，其中xi=[ai1ai2…ain]。

1.3 PCA與SVD的關(guān)系

由式(5)可見，PCA計算需要用到特征值計算，而且需要計算矩陣A的協(xié)方差矩陣。當(dāng)矩陣A的維數(shù)較大時，協(xié)方差矩陣的計算量很大，且協(xié)方差計算時存在舍入誤差[9]。鑒于此，文中提出一種基于SVD的PCA的計算方法，避免了計算協(xié)方差矩陣，減少了計算量，同時也避免了協(xié)方差矩陣計算時的舍入誤差。

對任意實矩陣，其奇異值分解表示為[13]

A=UΣVT

(8)

式中：U=(u1，u2，…，um)∈Rm×m、V=(v1，v2，…，vn)∈Rn×n分別稱為左奇異矩陣和右奇異矩陣，且都為正交矩陣，其中ui∈Rm×1、vi∈Rn×1分別稱為左奇異向量和右奇異向量；Σ=diag(σ1，σ2，…，σr)為對角矩陣，其元素為按降序排列的奇異值，即σ1≥σ2≥…≥σr≥0(r=min(m，n)是矩陣A的秩)。

構(gòu)造對稱矩陣AAT∈Rm×m，其特征值分解形式為

(9)

根據(jù)式(7)可得：

AAT=UΣVT(UΣVT)T=UΣVTVΣTUT

(10)

根據(jù)U、V的正交性有：

VTV=E

(11)

式(11)中E為單位矩陣。

又有：

(12)

把式(10)-(12)代入式(9)，得：

(13)

對比式(9)與式(13)可得：

(14)

其中，i=1，2，…，m。

同理可證，對于ATA∈Rn×n存在：

(15)

其中，i=1，2，…，n。

通過上面的推導(dǎo)可知：矩陣AAT或ATA的特征值是矩陣A或AT奇異值的平方，A的左奇異向量是AAT的特征向量，而A的右奇異向量是ATA的特征向量。

1.4 有效特征值數(shù)量規(guī)律

李振等[10]研究了協(xié)方差矩陣的有效特征值(即：非零特征值)與信號參數(shù)之間的關(guān)系，發(fā)現(xiàn)協(xié)方差矩陣的有效特征值數(shù)量取決于頻率個數(shù)，而與頻率大小及其幅值和相位無關(guān)，每一個頻率成分總是最多只產(chǎn)生兩個非零特征值，頻率的幅值越大，則其對應(yīng)的兩個特征值也越大，且這兩個特征值總是一前一后緊密排列。據(jù)此，根據(jù)信號中各頻率在幅值譜中的幅值大小排序，就可確定每一個頻率對應(yīng)的兩個特征值。

1.5 分量信號重構(gòu)

根據(jù)式(8)和式(14)求出與奇異值對應(yīng)的降序排列的特征值：λ1≥λ2≥…≥λr≥0。式(6)兩邊同時左乘αi并求和，可得：

(16)

假設(shè)某個特征頻率在幅值譜中的幅值大小排序為i，則選擇特征值分布圖中的第2i-1與第2i個特征值對應(yīng)的特征向量重構(gòu)分量信號，將待提取的k個特征頻率成分疊加得到近似矩陣Ak：

(17)

(18)

2 基于SVD原理的PCA算法

2.1 算法流程

根據(jù)協(xié)方差矩陣C的特征值與矩陣A的奇異值及信號頻率與有效特征值之間的關(guān)系，提出一種基于SVD原理的PCA特征頻率提取算法，其具體步驟如下。

步驟1 對一維離散信號a進(jìn)行零均值處理，去除直流分量，按式(7)構(gòu)造Hankel矩陣A；

步驟2 對矩陣A進(jìn)行SVD處理得到降序排列的奇異值σi及左奇異向量ui；

步驟3 根據(jù)協(xié)方差矩陣C的特征值(向量)與矩陣A的奇異值(向量)之間的關(guān)系，由式(14)求得降序排列的特征值λ1≥λ2≥…≥λr≥0及對應(yīng)的特征向量αi；

步驟4 根據(jù)特征值λi分布圖及特征頻率的幅值排序，選取待提取的k個頻率分量對應(yīng)的特征向量，按式(17)疊加得到近似矩陣Ak；

2.2 信號分析實例

構(gòu)造一個信號如下：

(19)

從圖1(c)、1(d)可知，采用文中的PCA算法提取的單個頻率成分與理想信號沒有相位偏移，幾乎與理想信號完全重合。由此說明，文中提出的PCA算法提取頻率的效果很好，且具有零相位特性。

圖1 提取含噪信號中的單個頻率

3 試驗分析

3.1 試驗裝置

采用課題組自主研發(fā)的能夠模擬汽輪機(jī)實際工況的大型滑動軸承試驗臺[14]進(jìn)行滑動軸承的減振性能試驗，該試驗臺主要包括驅(qū)動系統(tǒng)、潤滑系統(tǒng)、控制系統(tǒng)、采集系統(tǒng)及機(jī)械結(jié)構(gòu)等部分;轉(zhuǎn)子兩端由試驗滑動軸承支撐。為了有效抑制環(huán)境干擾采取以下措施：(1)采用空氣彈簧有效隔離地基傳遞到試驗臺的振動;(2)在小鑄鐵臺面和大鑄鐵臺面之間以及水泥地基與大鑄鐵臺面之間都增加了減振阻尼材料，以降低外界干擾對軸承-轉(zhuǎn)子系統(tǒng)的影響；(3)驅(qū)動電機(jī)與大鑄鐵臺面之間加丁晴橡膠，削弱電機(jī)振動傳遞的能量。試驗臺示意圖及實物圖如圖2、圖3所示。

圖2 試驗臺示意圖

圖3 試驗臺實物圖

3.2 位移傳感器測點布置

分別在試驗臺主軸兩端的A、B軸截面兩邊斜45°位置各布置一個電渦流傳感器，所測得的位移信號分別標(biāo)記為D1、D2、D3、D4，其中A端傳感器安裝如圖4所示。試驗所用的便攜式LMS SCADASD多功能數(shù)據(jù)采集系統(tǒng)如圖5所示。

圖4 位移傳感器測點布置(A端)

圖5 LMS數(shù)據(jù)采集系統(tǒng)

3.3 試驗方案

正式開展試驗前，必須進(jìn)行預(yù)備性試驗，目的是確保試驗裝置能夠正常運(yùn)行，滿足后續(xù)試驗要求。安裝好軸承試驗件且試驗臺調(diào)試完成后，設(shè)定液壓站潤滑油壓力為0.2 MPa、流量為20 L/min，調(diào)整轉(zhuǎn)子軸的轉(zhuǎn)速至200 r/min，穩(wěn)定運(yùn)行5～10 min，觀察監(jiān)測各參數(shù)及數(shù)據(jù)采集系統(tǒng)。若無異常，則由低轉(zhuǎn)速逐步提升到各工作轉(zhuǎn)速，由LMS數(shù)據(jù)采集系統(tǒng)采集升速過程中各試驗工況下的電渦流位移振動信號。文中僅對主軸A端兩個測點位移信號D1和D2進(jìn)行分析，采樣頻率為1 024 Hz、采集1 024個數(shù)據(jù)點。取其中任意3組試驗進(jìn)行分析：(1) 在試驗臺主軸兩端安裝軸承試驗件1，調(diào)試完成后，將轉(zhuǎn)速升至450 r/min (7.5 Hz)；(2) 在試驗臺主軸兩端安裝軸承試驗件2，調(diào)試完成后，將轉(zhuǎn)速升至540 r/min (9 Hz);(3)在試驗臺主軸兩端安裝軸承試驗件3，調(diào)試完成后，將轉(zhuǎn)速升至1 380 r/min (23 Hz)。

試驗中得到的信號如圖6所示，由圖6可知信號的波形難以辨別。各信號的頻譜圖如圖7所示，由圖7可知，頻譜中包含了隨機(jī)噪聲、轉(zhuǎn)頻及其高次諧波等眾多頻率成分。

圖6 各試驗中位移信號的時域波形

圖7 各試驗中位移信號的頻譜圖

4 軸心軌跡提純

根據(jù)軸心軌跡的形狀可以判斷旋轉(zhuǎn)機(jī)械的故障類型。研究表明：外“8”字形或香蕉形對應(yīng)不對中故障，內(nèi)“8”字形對應(yīng)油膜渦動故障，規(guī)則或不規(guī)則的花瓣形對應(yīng)轉(zhuǎn)子的碰磨故障等[15]。

利用圖6和圖7的各組試驗中兩個互相垂直的信號合成的軸心軌跡如圖8所示。由圖8可知，原始軸心軌跡雜亂無章，這是由于現(xiàn)場測試的振動信號受到了嚴(yán)重的噪聲干擾造成的。

圖8 原始軸心軌跡圖

對于旋轉(zhuǎn)機(jī)械的軸心軌跡提純，工程中普遍采用的方法是提取轉(zhuǎn)頻及其倍頻成分中幅值靠前的幾個(如 1×、2×、3×等)來合成軸心軌跡，才能得到清晰的圖形，進(jìn)而根據(jù)圖形進(jìn)行轉(zhuǎn)子的故障識別。而具體選取哪幾個特征頻率要結(jié)合信號的頻譜特征進(jìn)行分析。

下面采用文中提出的PCA算法對各信號進(jìn)行分析。首先，利用圖6中的各信號按式(7)構(gòu)造Hankel矩陣Aij(i=1，2，3;j=1，2)，Aij表示第i次試驗中的第j個信號所構(gòu)造的Hankel矩陣；其次，對矩陣Aij進(jìn)行SVD處理得到降序排列的奇異值及左奇異向量；然后，根據(jù)式(14)中的代數(shù)關(guān)系可求得對應(yīng)的特征值(如圖9所示)及特征向量。由圖7(a)、7(b)可知試驗1和試驗2中的信號1×和2×幅值較大，而其他的頻率幅值較??；故對這兩組試驗中的信號，選取前4個特征值及對應(yīng)的特征向量，采用平均法[12]重構(gòu)得到的信號如圖10所示，由它們合成的軸心軌跡如圖11所示。

圖9 特征值分布圖

圖10 從試驗1和2中提取的特征信號

圖11 試驗1和2中提純的軸心軌跡

圖11(a)中的軸心軌跡為橢圓形，由圖11(a)可知，軸心軌跡為橢圓形(長短軸尺寸相近)，說明轉(zhuǎn)子存在輕微不平衡故障，這是因轉(zhuǎn)子的加工工藝、制造誤差等因素造成且無法避免的，通常視為正常狀態(tài)。而11(b)中的軸心軌跡為“香蕉”形，說明轉(zhuǎn)子存在不對中故障。由圖10中的特征信號可以發(fā)現(xiàn)，正常信號(圖10(a))的特點是1×信號的幅值比2×的大得多，而不對中信號(圖10(b))中的2×幅值大于1×。

接下來分析試驗3的振動情況。由圖7(c)可知，試驗中信號的1×和3×幅值較大，而其他的頻率成分幅值較小。故對試驗3中的信號，選取前4個特征值及對應(yīng)的特征向量，采用平均法[12]重構(gòu)得到的信號如圖12所示。由圖12中的信號合成的軸心軌跡如圖13所示，由圖13可知，提純的軸心軌跡為花瓣形，說明轉(zhuǎn)子存在碰磨故障。

圖12 從試驗3中提取的特征信號

圖13 試驗3中提純的軸心軌跡

5 結(jié)語

本研究首先從理論上推導(dǎo)了矩陣奇異值(奇異向量)與協(xié)方差矩陣特征值(特征向量)之間的對應(yīng)關(guān)系，經(jīng)推導(dǎo)得出，協(xié)方差矩陣的特征值在數(shù)值上等于矩陣的奇異值的平方，左奇異向量對應(yīng)特征向量；然后，研究了Hankel矩陣方式下，主成分分析方法處理信號的分解及重構(gòu)原理，根據(jù)上述結(jié)論及有效特征值與頻率參數(shù)之間的關(guān)系，提出一種基于SVD原理的PCA特征頻率提取算法，并通過仿真信號驗證了算法的有效性；最后，將本研究提出的PCA算法應(yīng)用于大型滑動軸承試驗臺主軸的軸心軌跡提純，得到清晰的軸心軌跡，成功識別了轉(zhuǎn)子的不對中故障和碰磨故障。該算法同樣適用于其他旋轉(zhuǎn)機(jī)械，具有廣泛的應(yīng)用前景。