主分量分析在激勵源識別中的應用研究

董建超，楊鐵軍，李新輝，代 路

機械系統(tǒng)的激勵源識別對于系統(tǒng)優(yōu)化設計、聲振控制及故障診斷等方面有著重要的意義。實際機械設備趨于復雜化和聯(lián)合化，傳感器采集的信號往往是多個振動源耦合的結果，這使得振動信號的分析過程尤為困難。這時，基于統(tǒng)計理論的盲信號處理方法便顯現(xiàn)出其優(yōu)越性。盲信號處理方法在激勵源特性未知以及傳遞通道先驗知識不足的情況下，以信號的統(tǒng)計特征為評價標準，僅僅利用采集的混合信號來估計源信號并辨識傳遞通道特性。然而我們所采集的信號往往數據量龐大，如何進行數據壓縮并保留有效信息，成為盲信號處理中需要首先解決的問題。

主分量分析(Principal Component Analysis，PCA)是盲信號處理中非常有效的預處理方法［1－5］，由于算法簡單實用，被廣泛應用于數據壓縮、特征提取和去除噪聲等方面。Sutton［6］應用PCA對汽車懸架振動進行了分析，討論車內噪聲來源并進行主動控制。Cho等［7］結合PCA與傳遞路徑分析進行了多源環(huán)境下的輸入識別。Serviere等［8～9］應用改進的 PCA方法進行數據的降噪與白化預處理。相關文獻中應用PCA方法時，往往根據經驗確定主分量貢獻率閾值(例如90%)來進行不相關源數的上界估計［1－2］。如何定量的截取主分量，分析PCA的影響因素，具有重要的實際意義。

針對上述問題，本文引入鄰階分量信噪比(Signal Noise Ratio，SNR)作為主分量截斷的依據，通過算例仿真研究了混合矩陣條件數對PCA識別結果的影響。當PCA應用于機械系統(tǒng)激勵源識別時，混合矩陣對應為系統(tǒng)的頻率響應函數矩陣，激勵點位置、測點位置與系統(tǒng)固有頻率都會對頻響矩陣的條件數產生影響。針對工程中的激勵源，分析其線譜特征對PCA識別結果的影響。最后進行了簡支梁結構激勵源數目識別實驗，驗證了PCA不相關激勵源數目識別的效果。

1 PCA幾何闡釋及原理

1.1 PCA幾何闡釋

PCA基于信號的二階統(tǒng)計特性，將原始變量投影到坐標互相正交的新坐標系，使得具有最大方差的變量集中到新坐標系的少數幾個坐標軸上［1］，從而在數據信息丟失最少的原則下，用少數新的變量替代原來的多維變量。圖1以二維數據為例給出了主分量分析的幾何解釋。圖1(a)為原始數據散點圖，可以看出數據沿原始坐標系的x軸和y軸方向并不具有最大方差，數據存在冗余。圖1(b)給出方差隨坐標軸旋轉角度變化的曲線，當坐標軸旋轉至方差對應的峰值，即x'y'坐標系時，數據沿x'軸的分量方差達到峰值，包含了原始數據的大量信息。截取沿x'軸方向的分量即實現(xiàn)了原始數據的壓縮處理。

圖1 PCA在幾何上的解釋Fig.1 Geometrically illustration of PCA

1.2 主分量求解

假設存在噪聲n的條件下，混合信號x由信源信號 s經矩陣 A 混合得到［2－3］，即:

假設信源數目為m且相互之間不相關，則信源的協(xié)方差矩陣 Rs=ssT=E=Diag［e1，e2…em］，為 m 維對角方陣，對角線元素為各個信源的方差，表征信源信號的能量。δn表示噪聲對信源造成的干擾。假設噪聲為不相關的隨機噪聲，則噪聲的協(xié)方差矩陣δRn=δnδnT=N=Diag［μ1，μ2…μm］，亦為 m 維對角方陣，對角線元素為各個干擾成分的方差，表征干擾成分的能量。假定測點數目為n且n≥m，即混合系統(tǒng)為恰定或者超定系統(tǒng)。對混合矩陣A進行奇異值分解［1］:

式中:上標T代表轉置，U和V是奇異值分解中的左右奇異陣，且分別是p維和q維正交歸一方陣，即UTU=UUT=I，VTV=VVT=I;Σ為n×m維準對角陣，其準對角線元素σi(i=1，…m)是矩陣A的各階奇異值，且按照降序排列。

x的協(xié)方差矩陣為:

式中:ΣΣT(E+N)為m維對角陣，對角線元素為Rx的特征值，且前m個對角線元素按降序排列，其余n－m個元素為零。左奇異陣U中的列向量為Rx的特征向量。

將UT左乘x，得:

計算y的協(xié)方差矩陣:

可見Ry為對角陣，說明上述變換去除了x的相關性。y的各行向量稱為主分量，協(xié)方差矩陣的各階特征值反映了各階主分量的能量。對主分量進行截斷，舍棄能量很低的分量，可實現(xiàn)數據的降維壓縮和信源數目的判定。

2 算例仿真及分析

根據式(5)可知，信源信號能量的差異、噪聲干擾的大小以及不同的混合矩陣，都會對PCA的結果產生影響。本節(jié)主要進行信源能量均等且存在弱噪聲的情況下的PCA算例仿真，并分析混合矩陣條件數的大小對分量截斷產生的影響。

2.1 PCA分量截斷

本節(jié)針對式(1)所示模型進行PCA仿真。其中:信源信號s為兩個獨立的均勻分布且具有單位方差的隨機噪聲信號;A=［1.8，－0.3;1.1，－0.9;－0.2，1.3］為3×2維的混合矩陣，矩陣條件數為 1.92;噪聲n為三個與信源不相關且方差為0.01的隨機噪聲。

對混合信號x進行主分量分析，得到主分量y=［y1，y2，y3］T?；旌闲盘柡透麟A主分量的方差及相關系數絕對值如表1所示。由表1可以看出，各混合信號之間存在一定的相關性。進行主分量分析后，各階分量間相關系數變?yōu)?，去除了相關性并且各階分量按照方差大小進行了排序。

表1 方差及相關系數Tab.1 Variances and correlation coefficients

結合主分量貢獻率g，引入鄰階分量信噪比來綜合分析，進行分量的截斷。前r階主分量的貢獻率、第i階與第i+1階主分量的鄰階分量信噪比分別由式(6)及式(7)計算。按照相關文獻的經驗［1］，本文取 g≥0.9，即當貢獻率超過90%時，判定前r階分量包含數據樣本的能量。當鄰階分量信噪比很高時，說明較小的分量可作為噪聲被截斷。

上述算例中，一階分量的貢獻率為78.35%，前兩階分量的貢獻率達到99.86%;一、二階信噪比為5.61 dB，二、三階信噪比為21.75 dB。從而得出結論:不相關信號源的數目為2個，截取前兩階分量可以達到降低數據維數和減少冗余的目的。

2.2 混合矩陣條件數的影響

對于矩陣A，條件數定義如下［1］:

其中:A+1=［AHA］－1AH，為A的偽逆陣?！ぁ硎揪仃嚨?－范數，等于矩陣的最大奇異值。

令σmax，σmin分別為矩陣A的最大與最小奇異值，則‖A‖ =σmax，‖A+1‖ =1/σmin，可得

由式(5)可知，各階分量的能量與混合矩陣A的奇異值的平方呈正比，又由式(9)可得，第一階分量的能量與最后一階分量的能量之比與條件數的平方呈正比。當混合矩陣條件數很大時，主分量貢獻率會變大，一階分量分量與最后一階分量的信噪比也會變高，這樣會導致不相關源數的錯誤識別，主分量截斷隨即失效。

當存在噪聲等干擾因素時，條件數很大的混合矩陣會放大干擾對信源產生的不良影響。由式(5)可知，ΣΣTN是主分量中干擾因素所對應的部分，當條件數很大時，第一階主分量中的噪聲成分會隨之放大，進行分量截斷時會造成信源成分的丟失。

針對2.1節(jié)中的兩個具有單位方差隨機信源信號，在弱噪聲條件下(噪聲方差為0.01)產生100個隨機混合矩陣A進行混合，使前50個矩陣的條件數在1.0～1.5之間，后50個矩陣條件數在500～1 000之間。對每次混合后的數據分別進行PCA壓縮降維及源數識別。100次PCA的主分量貢獻率和鄰階分量信噪比結果如圖2所示。

如圖2所示，當混合矩陣條件數較小時(前50次)，一階分量貢獻率在65%左右，前兩階貢獻率基本達到100%，從而可以判定信源數目為2個。其中一、二階分量信噪比在3 dB左右，二、三階分量信噪比在22 dB左右，前兩階分量量級接近且均遠高于第三階分量。截斷前兩階分量，即可在不損失數據信息的前提下實現(xiàn)數據的壓縮。

圖2 PCA仿真結果Fig.2 Simulation results of PCA

當混合矩陣條件數很大時(后50次)，一階分量貢獻率已達到100%，而且一、二階分量信噪比就高達27 dB左右，截斷前一階分量即包含所有混合數據的能量，判定信源數為1個，會得出錯誤結論。此時，截斷前一階分量勢必會造成數據有效信息的丟失。

可見，當混合矩陣的條件數很大時，數據壓縮降維處理及源數識別會失效，所以在進行主分量分析時應避免混合矩陣條件數過大的情況。當PCA應用于實際的機械系統(tǒng)時，混合矩陣對應為系統(tǒng)的傳遞函數矩陣，情況變得更為復雜。下面將對機械系統(tǒng)中PCA激勵源識別問題進行分析。

3 頻域下振動響應主分量分析

機械振動在傳播過程中往往存在著耦合以及延遲等諸多效應，振動響應往往由多個激勵源卷積混合而成。對于穩(wěn)態(tài)過程的振動響應，可以在時域或者頻域內進行PCA激勵源識別［5］。時域法利用信號整體時間歷程的統(tǒng)計特性，處理過程相對簡便。比較而言，頻域法能夠在各個頻段去除相關性，更有利于定量分析。因此，本文在頻域下應用PCA識別機械系統(tǒng)激勵源數目。

3.1 頻域下機械系統(tǒng)PCA模型

對于多輸入多輸出的機械振動系統(tǒng)，假定激勵數為m，響應數為n，頻域下的響應為:

式中:ω為頻域內各個頻率值，xi(ω)表示第i個響應的傅里葉變換，hij(ω)表示第i個響應點與第j個激勵點之間的傳遞函數，頻域下對應為系統(tǒng)的頻率響應函數，fj(ω)表示第j個激勵的傅里葉變換。采用矩陣形式表示如下:

式中:

分別為系統(tǒng)的響應頻譜、激勵力頻譜和頻響函數矩陣。機械系統(tǒng)通?？山茷檎承宰枘嵯到y(tǒng)，由模態(tài)分析理論，頻率響應函數hij(ω)可表示為:

式中:r為系統(tǒng)模態(tài)階數，ωr、ξr、φr(L)分別為系統(tǒng)的 r階固有頻率、阻尼系數和振型函數，Lxi和Lfj分別為第i個響應點與第j個激勵點對應的位置。響應的功率譜密度矩陣為:

對SX進行特征值分解，特征值λi(ω)反映頻域內各階分量的能量。截取特征值的階數，即為PCA識別的不相關源數目。

為了保證PCA結果的有效性，希望頻響函數矩陣的條件數不要過大。由hij(ω)的表達式(13)可知，頻響矩陣與激勵點位置、測點位置和系統(tǒng)的固有頻率有關。下面分別分析這些因素對頻響矩陣條件數的影響。

3.2 不良激勵位置與測點位置的影響

由式(13)中的分子項可知，系統(tǒng)的頻響矩陣與激勵點位置和測點位置有關。如圖3(a)所示，當激勵源fi與fj靠近時，源的位置與各個測點位置之間的頻響函數非常接近，頻響矩陣H的第i列與第j列線性相關，矩陣條件數非常大，在激勵源識別過程中會將兩個源視作一個“組合源”來進行處理。如圖3(b)所示，當響應xi與xj靠近時，各個激勵源位置與兩個響應測點位置之間的頻響函數非常接近，頻響矩陣H的第i行與第j行線性相關，條件數也會過大，使得采集的數據冗余。實際測試分析中，為了保證激勵源識別結果有效并且降低測試的復雜程度，測點位置應該保持一定間隔。

圖3 兩種情況下的頻響函數矩陣Fig.3 Frequency response function matrices of two situations

本文針對簡支梁結構進行了不同激勵位置與測點位置的仿真研究，結果表明當激勵源位置靠近或者測點位置靠近時，PCA識別過程失效，由于篇幅所限，并未在此列出結果。

3.3 激勵源特性的影響

PCA方法基于信號的不相關性，所以激勵源的特性也會對識別結果造成一定的影響。工程中的激勵源具有明顯的線譜特征，如式(15)所示，激勵源f1(t)由p個不同頻率ω1i不同幅值a1i的正弦波構成，f2(t)由q個不同頻率ω2j不同幅值a2j的正弦波構成:

在此分析PCA能夠得到正確識別結果時，激勵源需要滿足的條件。

對于不同頻率的正弦信號，當采樣時間足夠長時，其不相關性與相位無關;同頻率的正弦信號，只有當相位差為π/2時，才滿足不相關性［10］，如式(16)所示。

因此，當 f1(t)和 f2(t)滿足 ω1i≠ω2j(?i，j)時，兩激勵源互不相關。由式(13)中的分母項可知，當激勵源的頻率接近系統(tǒng)固有頻率時，頻響函數矩陣的條件數非常大，系統(tǒng)的模態(tài)響應被激發(fā)出來并占據主導地位，固有頻率附近的分離結果可信度不高，激勵源識別分析時應當避開。

綜上所述，當激勵源f1(t)和f2(t)中頻率互不相等并且均不接近系統(tǒng)的固有頻率時，PCA能夠有效地識別激勵源數目。

4 PCA源數識別實驗

為了驗證寬頻范圍PCA的識別效果，本文分別采用了相關的白噪聲信號與不相關的白噪聲信號作為激勵源，進行簡支梁結構的激勵源數目識別實驗。兩種激勵源激勵的情況下，分別布置兩個測點與三個測點，進行兩點響應與三點響應的PCA激勵源數目識別研究。通過測點數目與不相關源數目相等、測點數目高于不相關源數目的情況下的識別結果，來驗證PCA的激勵源數目識別能力。

圖4 簡支梁實驗臺Fig.4 Simply supported beam rig

簡支梁實驗臺及測試儀器如圖4所示，梁結構的材料為45 號鋼，尺寸為0.64 m ×0.056 m ×0.008 m，邊界條件為簡支。為避免頻響函數矩陣條件數過大，激勵源位置相隔一定距離，前兩個測點布置在兩個激勵點位置，第三個測點遠離前兩個測點。定義圖中簡支梁左端為坐標原點，激勵點位置為:0.18 m，0.32 m;三個測點位置確定為:0.18 m，0.32 m，0.46 m。通過仿真計算理論頻響函數矩陣的條件數，除結構的固有頻率外，其余頻域范圍的矩陣條件數均較小，表明所選定的激勵位置與測點位置不會對識別過程產生不良影響。加速度響應的采樣頻率為4 096 Hz，采樣時間為4 s。采用Welch法［10］計算響應的自功率譜和互功率譜，分段數據的點數為2 048點，重疊率為50%，窗函數使用漢寧窗。

4.1 相關白噪聲激勵下源數識別

采用相關白噪聲信號對簡支梁結構進行激勵，加速度響應的自功率譜密度如圖5所示(其中由于第三個測點離激勵位置稍遠，所以其響應量級較小，但不影響識別結果)。

由圖5可以看出，響應在大部分頻段內相互重疊，僅僅根據響應的頻譜無法得到有效的信息來識別激勵源的數目。

圖5 相關激勵源情況下的加速度響應自功率譜Fig.5 Power spectral density of acceleration responses under correlated excitation sources

圖6 相關激勵源情況下的功率譜矩陣特征值Fig.6 Eigenvalues of power spectral density matrices under corrlated excitation sources

分別對兩點響應、三點響應構成的功率譜矩陣進行特征值分解，各階特征值如圖6所示。觀察兩圖，第一階特征值均比其余特征值高出近60 dB，圖6(b)中二、三階特征值比較接近，相差5 dB左右，并且數值都很小。由此說明，在測點數目與源數目相等或者高于源數目兩種情況下，一階分量量級都遠高于其余階分量，初步判定不相關激勵源數目為1個。

截取特征值得到的主分量貢獻率如圖7所示?？梢钥闯鰞煞N情況下截取前一階特征值，貢獻率在頻域范圍內接近100%，說明當測點數等于或高于源數時，一階分量基本包含了響應信號的能量。

鄰階分量信噪比如圖8所示。在固有頻率處，由于結構模態(tài)響應占據主導地位，所以結果稍差。除了固有頻率外的其余頻域范圍，一、二階分量信噪比為30 dB左右，圖8(b)中二、三階信噪比為3 dB左右。結果表明一階分量量級很高，可判定為激勵源，二、三階分量量級都很低，可以作為噪聲進行截斷。

圖7 相關激勵源情況下的的主分量貢獻率Fig.7 Contribution coefficients of principal components under correlated excitation sources

綜合圖6～8可得，對于相關白噪聲信號進行兩點激勵的情況下，當測點數目等于或高于源數目時，應用PCA進行激勵源識別，結合主分量貢獻率及鄰階信噪比等評價參數可以確定:不相關激勵源數目為1個。

4.2 不相關白噪聲激勵下源數識別

圖8 相關激勵源情況下的鄰階分量信噪比Fig.8 SNR between neighboring components under correlated excitation sources

采用兩個不相關白噪聲信號對簡支梁結構進行激勵，加速度響應的自功率譜密度如圖9所示。與4.1節(jié)圖5類似，響應在大部分頻段內相互重疊，僅僅根據響應的頻譜無法準確識別激勵源數。

分別對兩點響應、三點響應構成的功率譜矩陣進行特征值分解，進行不相關源數目識別分析，各階特征值如圖10所示。固有頻率處模態(tài)響應占主導地位，其余頻段內，一、二階特征值均相差5 dB左右，而圖10(b)中三階特征值與前兩階相比低近60 dB，并且數值很小。初步判定不相關激勵源數目為2個。

主分量貢獻率如圖11所示。由圖可看出截取前一階特征值，兩種情況下貢獻率在除固有頻率的頻段內約為65～75%左右，右圖中截取前兩階特征值的貢獻率基本達到了100%。由此說明，測點數目與源數目相等和高于源數目兩種情況下，一階分量都未能包含響應的所有能量，而前兩階分量已經基本包含了所有能量。

鄰階分量信噪比如圖12所示。圖中除固有頻率處較高外，一、二階分量信噪比為3～6 dB左右;右圖中二、三階信噪比達到30 dB，說明一、二階分量量級都很高，可判定為激勵源，三階分量量級很低，可以作為噪聲進行截斷。

圖9 不相關激勵源情況下的加速度響應自譜Fig.9 Power spectral density of acceleration responses under uncorrelated excitation sources

圖10 不相關激勵源情況下的功率譜矩陣特征值Fig.10 Eigenvalues of power spectral density matrices under uncorrelated excitation sources

圖11 不相關激勵源情況下的主分量貢獻率Fig.11 Contribution coefficients of principal components under uncorrelated excitation sources

綜合圖10～12可知，對于不相關白噪聲信號進行兩點激勵的情況下，當測點數等于或高于源數時，應用PCA進行激勵源識別，結合主分量貢獻率及鄰階信噪比，可以確定:不相關激勵源數目為2個。

當測點數目增加時，相關白噪聲激勵、不相關白噪聲激勵情況下的多點響應PCA識別結果類似，前一階或者前兩階主分量的特征值及貢獻率基本不變，后幾階分量特征值和貢獻率都很小，激勵源數目識別結果不變。由于篇幅所限，本文未將多點響應的識別結果列出。

圖12 不相關激勵源情況的鄰階分量信噪比Fig.12 SNR between neighboring components under uncorrelated excitation sources

綜合對比相關白噪聲與不相關白噪聲兩種激勵情況下的源數目識別結果可得:當測點數目等于或高于激勵源數目時，應用PCA進行激勵源識別，可以有效地分離出能量量級高的主分量作為激勵源，小量級分量由于其貢獻率低且鄰階分量信噪比高，可作為噪聲進行截斷。

5 結論

分析了混合矩陣條件數對PCA識別結果的影響，主分量貢獻率、一階分量與最后一階分量的信噪比以及噪聲干擾對信源產生的不良影響均隨條件數的增大而增大，從而導致主分量截斷失效。

當PCA應用于機械系統(tǒng)時，不良的激勵點位置與測點位置情況下的分離結果以及系統(tǒng)固有頻率附近的分離結果可信度都不高。研究了工程中激勵源的線譜特征對識別結果的影響，發(fā)現(xiàn)各個激勵源的頻率互不相等并且均不接近系統(tǒng)的固有頻率時，激勵源滿足不相關性條件，能夠保證識別結果有效。

進行了簡支梁結構激勵源數目識別實驗，引入鄰階分量信噪比作為主分量截取的依據，研究了測點數目等于或高于激勵源數目情況下的識別能力。結果表明:在多源激勵環(huán)境下，當測點數目等于或高于源數目時，對響應進行主分量分析能夠準確識別不相關激勵源的數目。以此為基礎的預處理，能夠確保激勵源盲識別過程更加可靠。

［1］Jolliffe I T.Principal component analysis［M］.Springer-Verlag，New York，2002.

［2］ Comon P，Jutten C.Handbook of blind source separation:independent component analysis and applications［M］.Elsevier，USA，2010.

［3］Hyvarinen A，Karhunen J，Oja E.獨立成分分析［D］.北京:電子工業(yè)出版社，2007.

［4］ Cichocki A，Amari S.Adaptive blind signal and image processing［M］.John Wiley＆ Sons，Ltd，2002.

［5］Liu X.Blind source separation methods and their mechanical application［D］.Australia:The University of New South Wales，2006.

［6］ Sutton T J.Active control of road noise inside vehicles［J］.Noise Control Eng，1994，42(4):137－147.

［7］ Cho M S，Kim K J.Indirect input identification in multisource environments by principal component analysis［J］.Mechanical Systems and Signal Processing，2002，16(5):873－883.

［8］ Serviere C，F(xiàn)abry P.Blind source separation of noisy harmonic signals for rotating machine diagnosis［J］.Journal of Sound and Vibration，2004，272:317－339.

［9］ Serviere C，F(xiàn)abry P.Principal component analysis and blind source separation of modulated sources for electro-mechanical systems diagnostic［J］.Mechanical Systems and Signal Processing，2005，19:1293－1311.

［10］ 胡廣書.數字信號處理［M］.北京:清華大學出版社，2003.