改進(jìn)PCA算法及其在轉(zhuǎn)子特征提取中的應(yīng)用*

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

(1.華南理工大學(xué)機(jī)械與汽車工程學(xué)院 廣州，510640) (2.廣州航海學(xué)院輪機(jī)工程學(xué)院 廣州，510725) (3.東莞職業(yè)技術(shù)學(xué)院實(shí)訓(xùn)中心 東莞，523808)

引 言

開展大型滑動軸承試驗(yàn)臺的設(shè)計(jì)和特性研究，對各種類型、結(jié)構(gòu)和參數(shù)的大型滑動軸承的研究具有重要意義[1]。 筆者所在團(tuán)隊(duì)自主研發(fā)了一種大型變支點(diǎn)滑動軸承試驗(yàn)臺，采用一種基于特征值差分譜理論的PCA信號處理方法來提純其軸心軌跡，進(jìn)而識別其工作狀態(tài)。

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

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

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

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

1.1 PCA信號分解原理

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

(1)

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

(2)

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

(3)

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

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

Cαi=λiαi

(4)

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

C=

(5)

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

1.2 Hankel矩陣構(gòu)造

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

(6)

矩陣X稱為重構(gòu)吸引子軌跡矩陣，也稱為Hankel矩陣。

1.3 奇異值與特征值的關(guān)系推導(dǎo)

對任意實(shí)矩陣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)

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

(9)

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

VTV=E

(10)

其中：E為單位矩陣。

又因?yàn)?/p>

(11)

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

(12)

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

(13)

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

1.4 有效主元分量選取

(14)

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

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

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

(15)

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

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

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

(16)

(17)

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

2.1 算法流程

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

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

2.2 信號分析實(shí)例

為了驗(yàn)證文中算法，任意構(gòu)造一個(gè)信號

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

(18)

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

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

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

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

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

3 試驗(yàn)分析

3.1 試驗(yàn)裝置簡介

試驗(yàn)采用自主研發(fā)的滑動軸承試驗(yàn)臺裝置[1]，如圖 3 所示，主要包括潤滑系統(tǒng)、驅(qū)動系統(tǒng)、機(jī)械結(jié)

圖3 試驗(yàn)臺實(shí)物圖Fig.3 Rotor test-bed

構(gòu)等部分組成。

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

在轉(zhuǎn)子兩端垂線斜45°位置各布置一個(gè)電渦流傳感器，分別標(biāo)記為di(i=1,2,3,4)，對應(yīng)的信號為Di(i=1,2,3,4)，其中靠驅(qū)動裝置端安裝的傳感器如圖4所示。試驗(yàn)所用的便攜式LMS SCADASD多功能數(shù)據(jù)采集系統(tǒng)如圖5所示。

圖4 驅(qū)動裝置端位移傳感器測點(diǎn)布置Fig.4 Location of displacement sensors near the driver

圖5 LMS數(shù)據(jù)采集系統(tǒng)Fig.5 LMS data acquisition system

3.3 振動信號采集

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

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

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

4 軸心軌跡提純

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

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

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

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

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

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

圖11 前4個(gè)分量信號的重構(gòu)結(jié)果Fig.11 Reconstruction results of the first four components

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

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

圖13 前4個(gè)分量信號的重構(gòu)結(jié)果Fig.13 Reconstruction results of the leading four components

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

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

4.2 與傳統(tǒng)PCA算法提純效果比較

傳統(tǒng)PCA算法具有依賴于研究者的信號分析經(jīng)驗(yàn)的缺陷，其根據(jù)累積貢獻(xiàn)率[8-11]來確定主成分的個(gè)數(shù)，累積貢獻(xiàn)率的表達(dá)式為

(19)

其中：Li為累計(jì)貢獻(xiàn)率；k為主特征值個(gè)數(shù)；m為有效特征值個(gè)數(shù)；δ為某個(gè)選定的閾值。

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

表1 前6個(gè)特征值及其貢獻(xiàn)率

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

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

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

5 結(jié)束語

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