孫鑫暉， 郝木明， 王淮維

(中國(guó)石油大學(xué)(華東)化學(xué)工程學(xué)院，山東 東營(yíng) 266555)

實(shí)驗(yàn)?zāi)B(tài)分析經(jīng)歷了幾十年的發(fā)展歷程，已經(jīng)從單自由度發(fā)展為多自由度，由單輸入單輸出發(fā)展為多輸入多輸出，由局部估計(jì)發(fā)展為整體估計(jì)［1]，并且新的識(shí)別方法不斷出現(xiàn)。其中PloyMAX方法［2，3]是近幾年發(fā)展起來(lái)的頻域識(shí)別算法，該方法也稱作多參考點(diǎn)最小二乘復(fù)頻域算法(poly－reference least square complex frequencydomain)［4]。由于比較好的解決了頻域算法中數(shù)值病態(tài)問(wèn)題，該方法可以處理很寬的頻率范圍和很高的模態(tài)階次。并且能夠產(chǎn)生非常清晰的穩(wěn)定圖，適用于大阻尼與高噪聲數(shù)據(jù)，目前已經(jīng)成為一種工業(yè)標(biāo)準(zhǔn)。

本文通過(guò)分析該算法實(shí)現(xiàn)過(guò)程中的矩陣結(jié)構(gòu)，給出了縮減正則方程的近似計(jì)算方法，從而避免了矩陣求逆以及乘法運(yùn)算，使得計(jì)算量大為降低，提高了模態(tài)參數(shù)識(shí)別速度。

1 PolyMAX算法描述

1.1 頻響函數(shù)的右矩陣分式模型

對(duì)于一個(gè)具有No個(gè)輸出，Ni個(gè)輸入的多輸入多數(shù)出系統(tǒng)(MIMO)，系統(tǒng)的頻響函數(shù) H(ω)∈CNo×Ni可以通過(guò)右矩陣分式模型進(jìn)行描述:

其中矩陣多項(xiàng)式 B(ω)∈CNo×Ni、A(ω)∈CNi×Ni定義為:

其中 Ωj(ω)=exp(－iωTs·j)為基函數(shù)(Ts為采樣周期)，n為模型階次，Aj、Bj(j=0…n)矩陣多項(xiàng)式系數(shù)。

1.2 正則方程的縮減

通過(guò)對(duì)式(1)右乘A(ω)進(jìn)行線性化處理，然后利用最小二乘法求解系數(shù)矩陣。但是直接求解存在以下缺點(diǎn)是:

(1)當(dāng)頻響函數(shù)規(guī)模較大、模型階次較高時(shí)導(dǎo)致正則方程的維數(shù)太高，直接求解難以實(shí)現(xiàn)。

(2)分母系數(shù)矩陣Aj不能獨(dú)立于分子系數(shù)矩陣Bj求解，而模態(tài)參數(shù)的求解只需要分母系數(shù)矩陣Aj。

文獻(xiàn)［4] 給出了一種通過(guò)縮減正則方程求解分母系數(shù)矩陣的方法，該方法能夠大幅減小正則方程的維數(shù)且無(wú)需求解分子系數(shù)矩陣Bj。這里給出縮減正則方程的表達(dá)式，具體推導(dǎo)參見(jiàn)文獻(xiàn)［4] 。

式(4)中Xo、Yo的表達(dá)式如下:

其中?為 kronecker積，Nf為識(shí)別頻段內(nèi)包含的譜線數(shù)。

1.3 通過(guò)系數(shù)矩陣得到模態(tài)參數(shù)

當(dāng)系數(shù)矩陣Aj確定之后，模態(tài)參數(shù)可以通過(guò)式(6)中友矩陣進(jìn)行特征值分解得到:

其中極點(diǎn)e－λiTs位于Λ的對(duì)角線上，頻率與阻尼可以通過(guò)式(7)求得，模態(tài)參數(shù)與向量位于V最后Ni行對(duì)應(yīng)的列向量。

2 PolyMAX算法的快速實(shí)現(xiàn)

PolyMAX算法的識(shí)別可以按照以下步驟進(jìn)行:

① 確定出最高計(jì)算階次nmax、最低計(jì)算階次nmin。

② 通過(guò)式(4)得到階次n對(duì)應(yīng)的Ro、So、To。

③ 計(jì)算縮減的正則方程M，然后通過(guò)最小二乘求解系數(shù)矩陣Aj。根據(jù)式(6)、式(7)得到模態(tài)參數(shù)。

④ 改變系統(tǒng)階次建立穩(wěn)定圖。

在計(jì)算 Ro、So、To的過(guò)程中，當(dāng) Xo、Yo維數(shù)很高時(shí)，一般不直接通過(guò)式(4)中矩陣乘法計(jì)算，而是采用FFT 運(yùn)算代替矩陣乘法［4，5]。當(dāng) Ro、So、To計(jì)算完成之后，計(jì)算量主要集中在正則方程中矩陣M上，因?yàn)樵诮⒎€(wěn)定圖的過(guò)程中，每一個(gè)階次的改變，都需要重新計(jì)算矩陣M，通過(guò)式(3)可以看出，這需要大量的矩陣求逆以及矩陣乘法運(yùn)算。對(duì)于復(fù)雜結(jié)構(gòu)，輸出數(shù)目N0很多、計(jì)算階次n很高，在這種情況下計(jì)算量尤其突出。

本文通過(guò)分析正則方程中矩陣R0的結(jié)構(gòu)，給出其逆矩陣的近似計(jì)算方法。從而避免了矩陣的求逆以及乘法運(yùn)算。通過(guò)將基函數(shù)表達(dá)式 Ω.j(ωk)=exp·(－iωkTs·j)=exp(－i(k－1)/(Nf－1)·j)帶入到X0，然后由式(4)得到Ro的表達(dá)式如下:

由式(8)可以得出:當(dāng) m=n時(shí)，［Ro]mn=Nf

當(dāng)m－n為偶數(shù)時(shí):

當(dāng)m－n為奇數(shù)時(shí):

因此Ro具有如下形式:

由于PolyMax為寬頻帶模態(tài)參數(shù)識(shí)別算法，識(shí)別范圍內(nèi)的譜線數(shù)Nf?1，因此:

其中In+1為n+1維的單位矩陣。則M可以寫成:

3 仿真算例

圖1 七自由度系統(tǒng)Fig.1 Seven DOF system

圖2 穩(wěn)定圖Fig.2 Stabilization diagram

圖3 渦輪盤測(cè)點(diǎn)布置(箭頭為參考點(diǎn))Fig.3 Measurement setup of gas turbine disk

表1 常規(guī)實(shí)現(xiàn)與快速實(shí)現(xiàn)的識(shí)別結(jié)果比較Tab.1 Identification results between normal and fast implementation

4 實(shí)測(cè)算例

實(shí)測(cè)算例來(lái)自于某航空發(fā)動(dòng)機(jī)渦輪盤模態(tài)實(shí)驗(yàn)，目的是驗(yàn)證當(dāng)數(shù)據(jù)規(guī)模比較大時(shí)，快速實(shí)現(xiàn)方法所能帶來(lái)的效率提高。實(shí)驗(yàn)通過(guò)錘擊法進(jìn)行測(cè)試，測(cè)點(diǎn)布置如圖3所示。其中FRF中包含168個(gè)輸出，3個(gè)輸入。分析頻率范圍0－8 000 Hz，譜線數(shù)Nf=1 600。在分析范圍內(nèi)大約有50多個(gè)結(jié)構(gòu)模態(tài)，根據(jù)式(14)中模型階次n與結(jié)構(gòu)模態(tài)數(shù)目Nm之間的關(guān)系，可以得出所需的模型階次n大約為100階。

圖4為在不同階次下的計(jì)算時(shí)間，每一次計(jì)算中階次范圍都為20(運(yùn)行環(huán)境為Pentium M 、主頻1.6 GHz、1 G內(nèi)存、Matlab6.5編程)。從圖中可以看出在計(jì)算階次不高的情況下，二者的計(jì)算時(shí)間相差不大。隨著計(jì)算階次的提高，快速實(shí)現(xiàn)方法所用時(shí)間明顯少于常規(guī)實(shí)現(xiàn)方法，并且階次越高，差距越明顯。這說(shuō)明快速實(shí)現(xiàn)方法適合于復(fù)雜結(jié)構(gòu)的大規(guī)模運(yùn)算。圖5為在最高計(jì)算階次100下得到的穩(wěn)定圖。在此情況下，常規(guī)實(shí)現(xiàn)所用時(shí)間為123.6 s，而快速實(shí)現(xiàn)方法所用時(shí)間為42.7s。圖6為識(shí)別的部分振型(前9階)。

圖4 計(jì)算時(shí)間的比較Fig.4 Comparison of computation time

圖5 穩(wěn)定圖Fig.5 Stabilization diagram

圖6 識(shí)別的振型Fig.6 Identified modeshapes

5 結(jié)論

本文給出了PolyMAX算法的一種快速實(shí)現(xiàn)方法，該方法通過(guò)避免縮減正則方程中的求逆運(yùn)算與矩陣乘法，減少了運(yùn)算量。實(shí)測(cè)算例表明:在數(shù)據(jù)量比較大的情況下，與常規(guī)實(shí)現(xiàn)相比較，采用快速實(shí)現(xiàn)后參數(shù)識(shí)別時(shí)間顯著減少，證明了該方法的可行性。并且輸出數(shù)目越多、系統(tǒng)階次越高的情況效果越明顯，因此該方法適合于復(fù)雜結(jié)構(gòu)的模態(tài)參數(shù)識(shí)別。

［1] 傅志方，華宏星.模態(tài)分析理論與應(yīng)用［M] .上海:上海交通大學(xué)出版社，2000.

［2] Peeters B，Van Der Auweraer H，Guillaume P，et al.The PolyMAX frequency-domain method:a new standard for modal parameter estimation ［J] .Shock and Vibration，2004，11:395－409.

［3] Peeters B，Guillaume P.Automotive and aerospace applictions of the LMS PolyMAX modal parameter estimation method［C] //Proc of the 22th International Modal Analysis Conference，Dearborn，USA，January，2004.

［4] Guillaume P，Verboven P，Vanlanduit S，et al.A polyreference implementation of the least-squares complex frequency domain estimator［C] // Proc of the 21th International Modal Analysis Conference.Kissimmee，USA，F(xiàn)ebruary，2003.

［5] Schoukens J，Rolain Y，Gustafsson F，et al.Fast calculation of least-squares for system identification［C] //Proc on the 37th IEEE Conference on Decision ＆ Control，F(xiàn)lorida，USA，December，1998，3408 －3410.

［6] Li Z F，Hua H X，Song H W，et al.Extracting modal parameters from structures undergoing ambient excitation excitation［J] .Journal of Shanghai Jiaotong University，2001，2:117－122.

［7] Shih C Y，Tsuei Y G，Allemang R J.Complex mode indication function and its applications to spatial domain parameter estimation［J] .Mechanical System and Signal Processing，1988，2(4):367 －377.