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

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

在上述的PCA應(yīng)用當(dāng)中，關(guān)鍵問題是如何確定有效主成分的個(gè)數(shù)。以往的研究策略習(xí)慣采用累積貢獻(xiàn)率[9]取定某個(gè)閾值的方法，閾值越大主成分的個(gè)數(shù)就越多，保留的信息量就越大，從而可能保留的噪聲成分也越多，但這種方法并不能任意提取感興趣的單個(gè)或者多個(gè)頻率成分。文獻(xiàn)[10]中通過仿真試驗(yàn)研究了特征值數(shù)量與信號(hào)基本參數(shù)(頻率、幅值和相位)的內(nèi)在關(guān)系，表明一個(gè)頻率最多產(chǎn)生兩個(gè)特征值，但并沒有嚴(yán)格的理論推導(dǎo)。

針對(duì)上述不足，文中將從理論上推導(dǎo)出有效特征值的數(shù)量規(guī)律及其排序規(guī)律，并將這兩種性質(zhì)統(tǒng)稱為PCA的幅值濾波特性?；诖颂匦蕴岢隽艘环N稱為主成分分析幅值濾波器的算法，并通過仿真信號(hào)及實(shí)際的轉(zhuǎn)子振動(dòng)信號(hào)的分析驗(yàn)證算法的有效性。最后將該算法應(yīng)用于大型滑動(dòng)軸承轉(zhuǎn)子的軸心軌跡提純，成功識(shí)別轉(zhuǎn)子的不對(duì)中故障。

1 PCA基本理論

1.1 PCA信號(hào)分解原理

假設(shè)xi∈R1×n、yi∈R1×n都是n維隨機(jī)向量，則主成分分析[11]是指用k個(gè)線性無關(guān)的新向量yi來線性表示m個(gè)n維初始向量xi，使得新變量的方差最大或降維損失最小。即有：

(1)

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

(2)

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

Cαi=λiαi

(3)

式中，矩陣C的計(jì)算式為

(4)

其中，cov(xi，xj)=E[(xi-E(xi))(xj-E(xj))T].

1.2 Hankel矩陣構(gòu)造

式(4)中矩陣X采用Hankel矩陣。對(duì)時(shí)域連續(xù)信號(hào)a(t)以采樣間隔Ts進(jìn)行離散化處理，采集N個(gè)數(shù)據(jù)，得到離散序列a=[a(1)，a(2)，…，a(N)]。利用此離散序列，可按下式構(gòu)造Hankel矩陣X：

(5)

令X=[x1，x2，…，xm]T，xi=[ai1，ai2，…，ain]，則矩陣X稱為Hankel矩陣。文獻(xiàn)[12]中證明了矩陣行數(shù)與列數(shù)越近似相等時(shí)，消噪效果越良好：當(dāng)信號(hào)長(zhǎng)度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矩陣。

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

求解式(3)可得協(xié)方差矩陣的特征向量αi及其特征值λi，將所有特征值降序排列有：λ1≥λ2≥…≥λr，r為Hankel矩陣X的秩且r= min(m，n)。用αi同時(shí)左乘式(1)等號(hào)兩邊，并只選擇感興趣的k個(gè)主成分相加可得近似矩陣：

(6)

(7)

至此，仍然存在一個(gè)問題需要解決：即式(6)中的k值(有效特征值的個(gè)數(shù))應(yīng)該如何確定，該問題一直是PCA應(yīng)用領(lǐng)域的研究熱點(diǎn)。文獻(xiàn)[10]中研究表明，特征值的數(shù)量與信號(hào)中的頻率個(gè)數(shù)存在兩倍關(guān)系，但僅僅是通過仿真試驗(yàn)的手段進(jìn)行研究，并沒有從理論上證明這一結(jié)論。因此，文中將從理論上對(duì)此進(jìn)行證明，并結(jié)合特征值大小與信號(hào)幅值大小的關(guān)系，進(jìn)一步研究主成分分析的本質(zhì)，提出一種主成分分析幅值濾波器算法。

2 有效特征值數(shù)量規(guī)律的理論推導(dǎo)

根據(jù)上一節(jié)可知，PCA算法是在特征值分解協(xié)方差矩陣的基礎(chǔ)上實(shí)現(xiàn)的，將其應(yīng)用于信號(hào)分析的關(guān)鍵是合理確定有效特征值的個(gè)數(shù)。文中將分3種情況對(duì)有效特征值的數(shù)量規(guī)律進(jìn)行討論，即分別討論信號(hào)包含單個(gè)頻率、多個(gè)頻率及隨機(jī)噪聲3種情況。

(1)單個(gè)頻率成分

若信號(hào)中僅含一個(gè)頻率成分，其表達(dá)式為：

x(t)=a1sin(ω1t+φ1)

(8)

式中，a1為信號(hào)的幅值，ω1為頻率，φ1為相位。

以采樣間隔Ts對(duì)信號(hào)x(t)進(jìn)行離散化處理，采集N個(gè)數(shù)據(jù)，得到離散序列X=[x(1)，x(2)，…，x(N)]。利用此離散序列，根據(jù)式(5)可構(gòu)造Hankel矩陣H，即：

(9)

根據(jù)歐拉公式可得：

(10)

將式(10)代入(9)并化解可得：

(11)

由式(11)可知，矩陣H1等于兩個(gè)子矩陣H11與H12之和，由于這兩個(gè)矩陣的每一行都線性相關(guān)，可以化簡(jiǎn)為秩為1的矩陣，即r(H11)=r(H12)=1。而根據(jù)文獻(xiàn)[12]的研究可知，矩陣H1與分量矩陣H11和H12的秩之間滿足r(H1)≤r(H11)+r(H12)=2。為了進(jìn)一步確定矩陣H1的秩的大小，根據(jù)矩陣秩的定義來進(jìn)行研究。

對(duì)矩陣H1存在2階行列式：

(12)

可見2階行列式不等于0，因此有r(H1)≥2，結(jié)合前述分析結(jié)果(r(H1)≤ 2)，可知r(H1)=2。而協(xié)方差矩陣C的秩r(C)=r(H1)=2，故對(duì)僅含有一個(gè)頻率成分的信號(hào)進(jìn)行PCA處理之后，協(xié)方差矩陣的特征值恒為2，即經(jīng)PCA處理之后一個(gè)頻率成分僅產(chǎn)生兩個(gè)非零特征值。

(2)k個(gè)頻率成分

若原始信號(hào)x(t)含有k個(gè)頻率成分，即：

(13)

以采樣間隔Ts對(duì)信號(hào)x(t)進(jìn)行離散化處理，采集N個(gè)數(shù)據(jù)，得到離散序列X=[x(1)，x(2)，…，x(N)]。利用此離散序列，根據(jù)式(5)可構(gòu)造Hankel矩陣H：

(14)

式中，Hi是利用第i個(gè)頻率分量對(duì)應(yīng)的Hankel矩陣，其表達(dá)式如下：

(15)

式中，N為信號(hào)的長(zhǎng)度，m=N-n+1。

同理，利用歐拉公式可將Hi改寫為如下形式：

(16)

將式(16)代入式(14)中，可得含有k個(gè)頻率成分信號(hào)的Hankel矩陣如下：

(17)

由式(17)可知，矩陣H由k個(gè)Hankel矩陣Hi求和得到，且其秩滿足r(Hi)=2，與頻率ωi、幅值ai及相角φi都無關(guān)。由于ω1≠ω2≠…≠ωk，可知各Hi之間是互相線性無關(guān)的，因此可知當(dāng)原始信號(hào)含有k個(gè)頻率成分時(shí)，H的秩為r(H)=2k。而協(xié)方差矩陣C的秩r(C)=r(H)=2k，故對(duì)含有k個(gè)頻率成分的信號(hào)進(jìn)行PCA處理之后，協(xié)方差矩陣的特征值恒為2k。

(3)噪聲對(duì)奇異值的影響

文獻(xiàn)[12]研究表明，隨機(jī)噪聲僅對(duì)奇異值大小產(chǎn)生影響，而不會(huì)改變奇異值的數(shù)量。而文獻(xiàn)[14]證明了特征值與奇異值之間存在平方關(guān)系(λi=σi2)。因此，隨機(jī)噪聲也只會(huì)影響特征值的大小，而不會(huì)改變其數(shù)量。

綜上所述，可得到如下重要結(jié)論：非零特征值的數(shù)量是信號(hào)x(t)所含的頻率個(gè)數(shù)的兩倍，若假設(shè)x(t)中含有k個(gè)頻率成分，對(duì)其進(jìn)行PCA處理之后可得到2k個(gè)非零特征值。簡(jiǎn)言之，一個(gè)頻率成分對(duì)應(yīng)兩個(gè)特征值，此即特征值的數(shù)量規(guī)律。

3 特征值排序規(guī)律的證明

文中第2節(jié)已經(jīng)證明了每個(gè)頻率成分對(duì)應(yīng)兩個(gè)特征值，那么對(duì)于一個(gè)確定的頻率應(yīng)該選擇哪兩個(gè)特征值，下面將會(huì)證明特征值的排序規(guī)律：特征值的大小取決于對(duì)應(yīng)頻率分量信號(hào)的幅值大小，信號(hào)的幅值越大，則其對(duì)應(yīng)的兩個(gè)特征值也越大。假設(shè)某個(gè)特征頻率在幅值譜中的幅值排序?yàn)閕，則其對(duì)應(yīng)第2i-1與第2i個(gè)特征值，選擇這兩個(gè)特征值對(duì)應(yīng)的特征向量重構(gòu)分量矩陣，即：

(18)

(19)

將式(2)代入式(19)可得：

(20)

由文獻(xiàn)[15]研究結(jié)果可知，Hankel矩陣X的能量滿足：

(21)

將特征值λi=σi2及式(21)代入式(20)可得：

(22)

4 PCA幅值濾波特性研究

4.1 PCA幅值濾波器概念的提出

從前述分析可知，PCA具有以下兩個(gè)重要的性質(zhì)：(1)特征值的數(shù)量規(guī)律，即一個(gè)頻率成分對(duì)應(yīng)兩個(gè)非零特征值；(2)特征值排序規(guī)律，即頻率分量的幅值越大，其對(duì)應(yīng)的兩個(gè)特征值也越大。文中將PCA的以上兩個(gè)性質(zhì)統(tǒng)稱為PCA的幅值濾波特性。

與其他的濾波方法相比，PCA不像陷波器那樣在濾除干擾信號(hào)的同時(shí)會(huì)造成有用信號(hào)的損失，也不像自適應(yīng)濾波器那樣需要參考輸入，更不像小波變換等帶通濾波器那樣濾波效果的好壞依賴于基函數(shù)的特性[12]。由PCA的濾波特性可知，其本質(zhì)上是一種幅值濾波器，根據(jù)各頻率分量信號(hào)的幅值大小對(duì)其進(jìn)行濾波，在滿足采樣定理的條件下，理論上可以提取出頻譜圖中的任意頻率。

4.2 基于PCA幅值濾波特性的信號(hào)分離算法

根據(jù)PCA的幅值濾波特性，提出一種基于PCA幅值濾波特性的信號(hào)分離方法，簡(jiǎn)稱PCA幅值濾波器(PCA-AF)。其具體步驟如下：

1)對(duì)給定的信號(hào)a(t)以采樣間隔Ts進(jìn)行離散化采樣并進(jìn)行零均值化處理，得到離散序列a=[a(1)，a(2)，…，a(N)]，利用該序列按式(5)構(gòu)造Hankel矩陣X。

2)計(jì)算協(xié)方差對(duì)矩陣C并對(duì)其進(jìn)行特征值分解，將得到的特征值按降序排列λ1≥λ2≥…≥λr，其對(duì)應(yīng)的特征向量為α1，α2，…，αr。

4.3 實(shí)例分析

(1)仿真信號(hào)分析

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

(19)

以1 024 Hz的采樣頻率對(duì)信號(hào)x(t)進(jìn)行離散化，并疊加信噪比為1.5 dB的高斯白噪聲，結(jié)果如圖1所示，可見噪聲干擾較嚴(yán)重。下面采用文中提出的PCA幅值濾波器算法進(jìn)行處理。

(a)波形圖

(b)頻譜圖

圖2 特征值曲線

(a)波形圖

(b)頻譜圖

(a)第1、2個(gè)特征分量

(b)第3、4個(gè)特征分量

(c)第5、6個(gè)特征分量

(d)第7、8個(gè)特征分量

圖4 不同分量的重構(gòu)結(jié)果

Fig.4 Reconstructed results of different components

綜上所述，上述分析結(jié)果驗(yàn)證了PCA確實(shí)能夠利用頻率分量的幅值大小進(jìn)行幅值濾波：依據(jù)特征值的數(shù)量規(guī)律和排序規(guī)律能夠準(zhǔn)確提取出感興趣的單個(gè)或者多個(gè)頻率成分，且幅值精度高、無相位偏移。

(2)轉(zhuǎn)子振動(dòng)特征提取

試驗(yàn)數(shù)據(jù)來源于課題組自主研發(fā)的滑動(dòng)軸承試驗(yàn)臺(tái)[16]，電渦流傳感器(D1、D2)安裝在轉(zhuǎn)子同一軸截面兩側(cè)相互垂直的斜45°方向(如圖5所示)。以采樣頻率2 048 Hz采集4 096個(gè)點(diǎn)，結(jié)果如圖6所示。由圖6(c)、(d)可知，位移信號(hào)的主要成分為1X和2X，同時(shí)存在轉(zhuǎn)頻的高次諧波成分、工頻干擾及隨機(jī)噪聲。

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

(a)D1的波形圖

(b)D2的波形圖

(c)D1的頻譜圖

(d)D2的頻譜圖

根據(jù)轉(zhuǎn)子動(dòng)力學(xué)理論，信號(hào)中二倍頻幅值較大，同時(shí)存在高次諧波成分，這是轉(zhuǎn)子不對(duì)中的典型特征。下面采用文中提出的算法，提取出這兩個(gè)頻率成分(1X和2X)。

利用圖6中的信號(hào)D1、D2分別構(gòu)造Hankel矩陣X1、X2∈R1 025×1 024，對(duì)其進(jìn)行特征值分解，得到對(duì)應(yīng)的1 024個(gè)特征分量，其對(duì)應(yīng)的特征值曲線如圖7所示(圖中只列出前50個(gè)值)。

(a)D1

(b)D2

(a)D1的波形圖

(b)D2的波形圖

(c)D1的頻譜圖

(d)D2的頻譜圖

狀況f/HzD1/mmD2/mm原始180.017150.01502360.019410.02567提純180.017110.01505360.019420.02569

5 基于軸心軌跡提純的轉(zhuǎn)子故障診斷

研究表明，軸心軌跡的形狀與故障類型存在密切關(guān)系，如外“8”字形或香蕉形對(duì)應(yīng)不對(duì)中故障，內(nèi)“8”字形對(duì)應(yīng)油膜渦動(dòng)故障等[10]。直接利用圖6中的兩個(gè)信號(hào)D1和D2合成的軸心軌跡如圖9所示，由圖可知，由于受到轉(zhuǎn)頻的高次諧波成分、電源工頻成分(50 Hz及150 Hz)及隨機(jī)噪聲的干擾，軸心軌跡雜亂無章，無法據(jù)此辨別轉(zhuǎn)子的故障類型。

圖9 轉(zhuǎn)子的原始軸心軌跡

因此，需要對(duì)圖9中的軸心軌跡進(jìn)行提純，其本質(zhì)是分別對(duì)信號(hào)D1和D2進(jìn)行降噪及濾除多余的頻率成分之后，再利用信號(hào)的提純結(jié)果合成軸心軌跡。本文中，信號(hào)的提純工作在4.3小節(jié)中已經(jīng)完成(如圖8所示)，因此，直接利用圖8中的提純結(jié)果合成軸心軌跡，結(jié)果如圖10所示，由圖可知，軸心軌跡為“8”字形，說明轉(zhuǎn)子存在不對(duì)中故障。這是由于在升速過程中，轉(zhuǎn)子兩端的油膜力的大小不一致引起的轉(zhuǎn)子中心與軸承中心發(fā)生偏離所致。

圖10 轉(zhuǎn)子的提純軸心軌跡

此外，由PCA幅值濾波器算法流程可知，該算法是一種依據(jù)實(shí)際信號(hào)的頻譜特征及各頻率的幅值大小進(jìn)行特征提取的，文中雖然只分析了轉(zhuǎn)子的不對(duì)中故障的情況，但是對(duì)于轉(zhuǎn)子的不平衡、動(dòng)靜件、油膜渦動(dòng)等其他故障類型是同樣適用的。顯然，文中算法也可以推廣到其他類型的旋轉(zhuǎn)機(jī)械設(shè)備的故障診斷。

6 結(jié)論

(1)從理論上推導(dǎo)了非零特征值的數(shù)量規(guī)律，研究表明：對(duì)含有k(k=1，2，…)個(gè)頻率的一維信號(hào)進(jìn)行PCA處理，恒產(chǎn)生2k個(gè)非零特征值，且特征值的數(shù)量不受信號(hào)的幅值、相位及隨機(jī)噪聲的影響。

(2)研究了特征值的排序規(guī)律，即特征值大小與信號(hào)幅值的對(duì)應(yīng)關(guān)系，指出某個(gè)頻率分量對(duì)應(yīng)的幅值越大，該頻率對(duì)應(yīng)的特征值也越大。利用該性質(zhì)可以確定特征值出現(xiàn)的順序。

(3)根據(jù)特征值的數(shù)量規(guī)律及排序規(guī)律可以看出，特征值分解本質(zhì)上是一種幅值濾波器，文中將這兩個(gè)性質(zhì)統(tǒng)稱為PCA的幅值濾波特性?；诖颂匦?，提出了一種稱為PCA幅值濾波器(PCA-AF)的信號(hào)分離算法，通過仿真信號(hào)及實(shí)際的轉(zhuǎn)子信號(hào)驗(yàn)證了算法的有效性。

(4)將PCA-AF算法應(yīng)用于大型滑動(dòng)軸承試驗(yàn)臺(tái)的軸心軌跡提純，成功識(shí)別了轉(zhuǎn)子的不對(duì)中故障。文中提出的算法同樣適用于其他旋轉(zhuǎn)機(jī)械的故障診斷及特征提取。