改進PCA算法及其在轉子特征提取中的應用*

李偉光, 郭明軍, 楊期江, 趙學智, 李國臣

(1.華南理工大學機械與汽車工程學院 廣州，510640) (2.廣州航海學院輪機工程學院 廣州，510725) (3.東莞職業技術學院實訓中心 東莞，523808)

引 言

開展大型滑動軸承試驗臺的設計和特性研究，對各種類型、結構和參數的大型滑動軸承的研究具有重要意義[1]。 筆者所在團隊自主研發了一種大型變支點滑動軸承試驗臺，采用一種基于特征值差分譜理論的PCA信號處理方法來提純其軸心軌跡，進而識別其工作狀態。

滑動軸承的失效形式多種多樣，常用的診斷方法包括油樣分析、聲發射檢測和振動分析等，其中通過監測轉子位移信號的振動分析方法最為簡單有效[2]。采集轉子同一軸截面相互垂直布置的兩個位移傳感器的振動位移信號，將其合成軸心軌跡，軸心軌跡的辨識主要研究軸心軌跡圖像的特征提取和識別問題，其前提也需要有清晰的軸心軌跡圖作為依據[3]。然而實際采集的振動信號通常會受到干擾因素的影響，所以，需要對原始信號合成的軸心軌跡進行提純，其本質是對原始振動信號進行降噪處理后合成軸心軌跡[4]。常用的提純方法包括：數字或模擬低通濾波法[3]、小波變換和小波包變換[5]、粒子群算法、EMD降噪[6]及形態濾波[7]等。上述降噪方法一般會存在相位偏移、降噪畸變大及頻帶選擇不明確等問題[4]，因此，改進已有方法或尋求新的方法都是一種可行的方案。

近年來，主成分分析技術被廣泛應用于消噪、故障診斷、特征提取及模式識別等領域[8-11]。劉永斌等[8]根據累積貢獻率選取主元個數，對雙層圓柱殼體機械噪聲數據進行降維。尚前明等[9]將PCA直接應用到船舶柴油機的故障監測中，通過PCA實現柴油機熱工參數的降維，從而準確識別出柴油機的異常狀態。Li等[10]采用PCA方法對風力機滾動軸承故障信號的特征矩陣進行降維，通過濾除特征中的冗余信息使得支持向量機(support vector machine,簡稱SVM)的分類取到了更加有效的信息，精度得以提高并極大減少了計算量。Seghouane等[11]針對主載荷矢量的密集結構使得主載荷矢量的降維難以解釋的問題，提出了一種自適應塊稀疏PCA方法。

主成分分析的關鍵是有效主元個數的確定，目前多數研究都是根據累計貢獻率[8-11]取定某個閾值的方法選擇有效主元的個數，閾值越大主成分的個數就越多，保留的信息量就越大，從而可能保留的噪聲成分也越多[12]。趙學智等[13]提出用奇異值差分譜來描述有用信號和噪聲的奇異值性質差異性，對于去除直流分量的信號，根據其差分譜首個峰值的位置可以自動選擇有效奇異值的個數，并能有效抑制噪聲成分的影響。將差分譜方法與 PCA技術融合, 用于描述主成分和次要成分的協方差矩陣特征值的差異性，以及探討奇異值與特征值之間的關聯等兩個問題目前還鮮有文獻報道。筆者結合裝配了新型結構滑動軸承的試驗對這兩個問題開展研究, 并將具體的 PCA 算法應用到大型滑動軸承試驗臺的轉子軸心軌跡的提純上，提純效果優于傳統PCA算法。

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

1.1 PCA信號分解原理

主成分分析[14]是指用k個n維的新變量y1,y2,…,yk來線性表示n維初始變量x1,x2,…,xm(m≥k),使得新變量的方差最大或降維損失最小,即有

(1)

其中：系數αi=(αi1,αi2,…,αim)T(i=1,2,…,k)為協方差矩陣C中降序排列的第i特征值λi對應的特征向量，且αi滿足

(2)

將式(1)改寫為分量形式有

(3)

其中：yi∈R1×n,i=1,2,…,k。

協方差矩陣C的特征方程為

Cαi=λiαi

(4)

協方差矩陣C的計算式為

C=

(5)

其中：cov(xi,xj)=E[(xi-E(xi))(xj-E(xj))T]。

1.2 Hankel矩陣構造

在式(5)中矩陣X采用Hankel矩陣，可利用一維零均值的離散信號a=[a(1),a(2),…,a(N)] ，按照以下方法構造Hankel矩陣

(6)

矩陣X稱為重構吸引子軌跡矩陣，也稱為Hankel矩陣。

1.3 奇異值與特征值的關系推導

對任意實矩陣X,其奇異值分解表示[12]為

X=UDVT

(7)

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

(8)

根據式(7)可得

(9)

根據U，V的正交性有

VTV=E

(10)

其中：E為單位矩陣。

又因為

(11)

把式(9～11)代入式(8)得

(12)

對比式(8)與式(12)可得

(13)

式(13)表明, 奇異值和特征值之間是一種平方關系。

1.4 有效主元分量選取

(14)

其中:bi形成的序列B=(b1,b2,…,br-1)為特征值差分譜。

差分譜反映了相鄰兩個特征值之間的變化趨勢，如果特征值在某個位置變化較大，在差分譜中將出現一個峰值，而在所有的峰值當中必然存在一個最大值bk。此時，k處發生一個最大突變，這種突變反映了有效主成分和次要成分的相關性的差異，代表主成分和次要成分的分界，據此確定的k值即為主成分的個數。

為了研究奇異值差分譜與特征值差分譜的聯系，將式(13)代入式(14)可得

(15)

從式(15)可以看出，特征值差分譜可以通過奇異值的平方差計算出來 ，這種平方運算可以將奇異值差分譜放大，從而使主成分峰值特征更加突出。從這一點而言，特征值差分譜是奇異值差分譜理論的升華。

1.5 分量信號重構

求解式(4)可得協方差矩陣的特征向量αi和λi，將所有特征值降序排列有λ1≥λ2≥…≥λm。用特征向量αi同時左乘式(3)等號兩邊，并只選擇前k個主成分相加可得

(16)

(17)

2 基于差分譜理論的PCA算法

2.1 算法流程

根據特征值差分譜的性質，提出一種基于特征值差分譜理論的PCA算法(簡稱差分PCA算法)，具體步驟如下：

1) 對經零均值化處理的一維離散信號a，按式(6)構造Hankel矩陣X；

2.2 信號分析實例

為了驗證文中算法，任意構造一個信號

x(t)=1.5sin(10πt+0.5)+2.0sin(20πt+0.8)

(18)

以1 024 Hz的采樣頻率對信號x(t)進行離散化,并疊加信噪比為0.15 dB的高斯白噪聲,結果如圖1所示。

圖1 原始含噪信號Fig.1 Original noisy signal

圖2 使用差分譜提取不同分量信號的過程Fig.2 Process to extract signals using difference spectrum

為了弄清楚差分譜不同峰值位置代表的含義，筆者逐個進行分析。從圖2(a)可知，差分譜的第1個峰值出現在第2個序號，因此利用前2個分量進行重構，結果如圖2(c)所示，圖中虛線為理想信號中的2.0sin(20πt+0.8)成分，該提純信號代表原始信號x(t)中的2.0sin(20πt+0.8)信號分量。另外已知，最大峰值序號對應的前4個分量提純的信號包含了x(t)中的所有頻率成分，那么利用最大峰值與其前面的第1峰值之間的第3，4個分量重構的信號必然為x(t)中的1.5sin(10πt+0.5)；利用第3，4個分量重構的信號如圖2(d)所示，其中虛線為理想的1.5sin(10πt+0.5)，說明該提純信號為原信號x(t)中的1.5sin(10πt+0.5)，提純結果與前述分析一致。

綜上所述，根據特征值差分譜最大峰值位置能夠可靠地提取純信號中的主成分，通過不同峰值之間的分量信號的組合可以提取出不同的頻率成分，且頻率對應的幅值越大其對應的峰值位置越靠前。

3 試驗分析

3.1 試驗裝置簡介

試驗采用自主研發的滑動軸承試驗臺裝置[1]，如圖 3 所示，主要包括潤滑系統、驅動系統、機械結

圖3 試驗臺實物圖Fig.3 Rotor test-bed

構等部分組成。

3.2 位移傳感器測點布置

在轉子兩端垂線斜45°位置各布置一個電渦流傳感器，分別標記為di(i=1,2,3,4)，對應的信號為Di(i=1,2,3,4)，其中靠驅動裝置端安裝的傳感器如圖4所示。試驗所用的便攜式LMS SCADASD多功能數據采集系統如圖5所示。

圖4 驅動裝置端位移傳感器測點布置Fig.4 Location of displacement sensors near the driver

圖5 LMS數據采集系統Fig.5 LMS data acquisition system

3.3 振動信號采集

筆者對某次試驗中(采樣頻率為 1 024 Hz)試驗臺主軸轉速分別為 450 r/min 和 1 080 r/min 的驅動端主軸振動位移信號(D1和D2)進行分析，其時域波形及頻譜圖分別如圖 6、圖7 所示。在兩種轉速工況下，信號的頻譜成分在整個頻帶內存在隨機噪聲和系統工作時激發出的高頻諧波分量。

圖6 信號去除直流分量的結果(450 r/min)Fig.6 Signals with D C components removed at 450 r/min

圖7 信號去除直流分量的結果(1 080 r/min)Fig.7 Signals with D C components removed at 1 080 r/min

4 軸心軌跡提純

4.1 采用差分PCA算法提純軸心軌跡

利用處理過的D1，D2兩組信號合成的轉子軸心軌跡如圖8所示，但不能清晰表明軸心軌跡。

圖8 原始軸心軌跡圖Fig.8 Original orbits of the rotor

采用筆者提出的差分PCA算法,分別利用處理過的振動位移信號D1和D2構造513×512維的Hankel矩陣Xi(i=1,2,3,4)；然后，按式(5)分別計算協方差矩陣Ci(i=1,2,3,4)，并對其進行特征值分解，經換算得到原始信號的特征值及其差分譜曲線如圖9所示(E1，E2分別代表信號D1，D2的特征值)，圖中只給出了前30個值的結果。由圖9可知：前4個特征值較大，第5個及其以后的特征值相對較小; 差分譜中前兩峰值的幅值較大。這與轉子旋轉時的1倍頻與2倍頻幅值最為突出、且2倍頻與1倍頻幅值接近[16]的結論吻合。

圖9 特征值及其差分譜曲線Fig.9 Eigenvalue curves and their difference spectrum

若選取差分譜中第2個峰值對應的前4個分量，按照上述步驟對信號進行同樣的處理，結果如圖11所示。由圖11可知，此時提取的是原始信號的1X(7.5 Hz)及3X(22.5 Hz)成分，而與圖10的結果比較可知，前2個峰值之間的第3，4分量對應的是3X(22.5 Hz)，從而驗證了2.2節中的“通過差分譜不同峰值之間的分量信號組合可以提取出不同的頻率成分”的結論。

圖11 前4個分量信號的重構結果Fig.11 Reconstruction results of the first four components

利用圖10，11中的信號合成的軸心軌跡分別如圖12(a，b)所示。由圖12(a)可知軸心軌跡為長短軸相差不大的橢圓形，說明轉子存在輕微不平衡現象。由圖12(b)可知軸心軌跡為花瓣形，其特征是各小瓣到瓣心的距離r較大且互相分離。

圖12 工況1的提純軸心軌跡Fig.12 The purified axis trajectories at condition 1

圖13 前4個分量信號的重構結果Fig.13 Reconstruction results of the leading four components

由圖13中的信號合成的軸心軌跡如圖14(a)所示，該圖為各花瓣相交的花瓣形，但瓣心到花瓣的距離r是減小的，由此表明相對于工況1，工況2轉子的軸心軌跡波動量更大，但是軸心振動幅值是減小的。本試驗臺裝配的是一種不同于常規圓軸承的特殊結構滑動軸承，由于本研究主要側重于該軸承結構轉子的振動軸心軌跡特征的提取，對低轉速下產生花瓣型軸心軌跡的機理需要進一步深入研究。若選取差分譜中前2個峰值之間對應的第3，4個分量，按照上述步驟對信號進行同樣的處理，然后合成的軸心軌跡如圖14(b)所示，軸心軌跡的橢圓形長短軸與工況1相近，表明與工況1的動不平衡量基本一致。兩組試驗未調整過軸承與轉子系統的任何參數，其動不平衡量不會發生改變，但不同組試驗測試采集的數據存在誤差，會導致微小的差異。

圖14 工況2的提純軸心軌跡Fig.14 The purified axis trajectories at condition 2

4.2 與傳統PCA算法提純效果比較

傳統PCA算法具有依賴于研究者的信號分析經驗的缺陷，其根據累積貢獻率[8-11]來確定主成分的個數，累積貢獻率的表達式為

(19)

其中：Li為累計貢獻率；k為主特征值個數；m為有效特征值個數；δ為某個選定的閾值。

現以工況1的信號為例，利用傳統PCA方法進行處理。得到的特征值分布圖與圖9(a，b)中的一致。由圖可知，當i≤14時有λi>1，故取m=14。為了便于分析，將k=1,2,…，6時，相應的特征值與累積貢獻率列于表1中，此時分別利用提純得到前k個分量合成的軸心軌跡如圖15所示。

表1 前6個特征值及其貢獻率

圖15 使用傳統PCA算法提純的軸心軌跡Fig.15 Purified axis orbits with traditional PCA algorithm

由圖15可知，貢獻率不同，提純得到的軸心軌跡形狀也不同。對比圖12與15可知，圖15(b，d)分別與圖12(a，b)基本一致，說明只有當k=2，4時，根據貢獻率確定的分量才是原始信號的主成分，而當k取其他值時，提純的結果并不準確。

在實際應用當中，累積貢獻率的閾值并不容易確定，嚴重依賴于研究者的信號分析經驗，這使得其應用具有一定的局限性。而筆者提出的特征值差分譜方法，一方面可以根據差分譜最大峰值位置自動篩選主成分的個數；另一方面，通過不同譜峰之間的分量信號的組合可以提取出不同的頻率成分。這些是傳統PCA方法所不具備的。

5 結束語

筆者提出利用協方差矩陣特征值差分譜的概念來描述有效主成分與次要成分的特征值差異性，根據差分譜的最大峰值位置可自動選擇有效主成分的個數，從而解決了主成分分析的關鍵問題。研究了Hankel矩陣方式下PCA方法處理信號的原理，并提出一種基于特征值差分譜理論的PCA算法，通過仿真信號驗證了該算法的有效性。研究表明：通過差分譜不同譜峰之間的分量信號的組合可以提取出不同的頻率成分。將筆者提出的差分PCA算法用于裝配了特殊結構滑動軸承的轉子軸心軌跡提純，效果優于傳統PCA算法，可直觀表明轉子的動不平衡特征狀態。