基于SVD原理的PCA特征頻率提取算法及其應(yīng)用

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

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

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

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

(1)

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

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

1 Hankel矩陣方式下PCA信號(hào)處理原理

1.1 PCA信號(hào)處理原理

主成分分析[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個(gè)特征值λi對(duì)應(yīng)的特征向量，且αi滿足：

(3)

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

(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矩陣。對(duì)于一維零均值的離散信號(hào)a=[a(1)a(2) …a(N)]。利用此信號(hào)按照以下方法構(gòu)造Hankel矩陣：

(7)

矩陣A稱為重構(gòu)吸引子軌跡矩陣，也稱為Hankel矩陣。針對(duì)Hankel矩陣的構(gòu)造問題，趙學(xué)智等[12]證明了噪聲去除量與列數(shù)呈拋物線的對(duì)稱關(guān)系，并得出結(jié)論：如果信號(hào)長度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計(jì)算需要用到特征值計(jì)算，而且需要計(jì)算矩陣A的協(xié)方差矩陣。當(dāng)矩陣A的維數(shù)較大時(shí)，協(xié)方差矩陣的計(jì)算量很大，且協(xié)方差計(jì)算時(shí)存在舍入誤差[9]。鑒于此，文中提出一種基于SVD的PCA的計(jì)算方法，避免了計(jì)算協(xié)方差矩陣，減少了計(jì)算量，同時(shí)也避免了協(xié)方差矩陣計(jì)算時(shí)的舍入誤差。

對(duì)任意實(shí)矩陣，其奇異值分解表示為[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)為對(duì)角矩陣，其元素為按降序排列的奇異值，即σ1≥σ2≥…≥σr≥0(r=min(m，n)是矩陣A的秩)。

構(gòu)造對(duì)稱矩陣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)

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

(14)

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

同理可證，對(duì)于ATA∈Rn×n存在：

(15)

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

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

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

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

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

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

(16)

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

(17)

(18)

2 基于SVD原理的PCA算法

2.1 算法流程

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

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

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

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

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

2.2 信號(hào)分析實(shí)例

構(gòu)造一個(gè)信號(hào)如下：

(19)

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

圖1 提取含噪信號(hào)中的單個(gè)頻率

3 試驗(yàn)分析

3.1 試驗(yàn)裝置

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

圖2 試驗(yàn)臺(tái)示意圖

圖3 試驗(yàn)臺(tái)實(shí)物圖

3.2 位移傳感器測(cè)點(diǎn)布置

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

圖4 位移傳感器測(cè)點(diǎn)布置(A端)

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

3.3 試驗(yàn)方案

正式開展試驗(yàn)前，必須進(jìn)行預(yù)備性試驗(yàn)，目的是確保試驗(yàn)裝置能夠正常運(yùn)行，滿足后續(xù)試驗(yàn)要求。安裝好軸承試驗(yàn)件且試驗(yàn)臺(tái)調(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)測(cè)各參數(shù)及數(shù)據(jù)采集系統(tǒng)。若無異常，則由低轉(zhuǎn)速逐步提升到各工作轉(zhuǎn)速，由LMS數(shù)據(jù)采集系統(tǒng)采集升速過程中各試驗(yàn)工況下的電渦流位移振動(dòng)信號(hào)。文中僅對(duì)主軸A端兩個(gè)測(cè)點(diǎn)位移信號(hào)D1和D2進(jìn)行分析，采樣頻率為1 024 Hz、采集1 024個(gè)數(shù)據(jù)點(diǎn)。取其中任意3組試驗(yàn)進(jìn)行分析：(1) 在試驗(yàn)臺(tái)主軸兩端安裝軸承試驗(yàn)件1，調(diào)試完成后，將轉(zhuǎn)速升至450 r/min (7.5 Hz)；(2) 在試驗(yàn)臺(tái)主軸兩端安裝軸承試驗(yàn)件2，調(diào)試完成后，將轉(zhuǎn)速升至540 r/min (9 Hz);(3)在試驗(yàn)臺(tái)主軸兩端安裝軸承試驗(yàn)件3，調(diào)試完成后，將轉(zhuǎn)速升至1 380 r/min (23 Hz)。

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

圖6 各試驗(yàn)中位移信號(hào)的時(shí)域波形

圖7 各試驗(yàn)中位移信號(hào)的頻譜圖

4 軸心軌跡提純

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

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

圖8 原始軸心軌跡圖

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

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

圖9 特征值分布圖

圖10 從試驗(yàn)1和2中提取的特征信號(hào)

圖11 試驗(yàn)1和2中提純的軸心軌跡

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

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

圖12 從試驗(yàn)3中提取的特征信號(hào)

圖13 試驗(yàn)3中提純的軸心軌跡

5 結(jié)語

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