FFT譜連續(xù)細(xì)化分析在頻域模態(tài)參數(shù)識(shí)別中的應(yīng)用

楊 陳，代洪慧

（中國(guó)水利水電科學(xué)研究院 工程抗震研究中心，北京 100048）

1 研究背景

經(jīng)典的傳遞函數(shù)求解方法，一般都是對(duì)輸入輸出數(shù)據(jù)分段（可以有部分重疊），然后對(duì)每段數(shù)據(jù)分別進(jìn)行FFT變換，進(jìn)而得到自功率譜和互功率譜。將分段得到的自功率譜和互功率譜進(jìn)行平均，利用平均了的互譜與自譜的比值得到傳遞函數(shù)。在計(jì)算中，為了充分地利用測(cè)試數(shù)據(jù)、消除混雜在測(cè)試信號(hào)里的噪聲，一般在取段時(shí)都要使部分?jǐn)?shù)據(jù)重疊（一般取50%），且消除每段的趨勢(shì)項(xiàng)，加上合適的窗函數(shù)[1-2]。

由于譜域的分辨率為：df=fs/nFFT（fs為時(shí)域信號(hào)的采樣頻率，nFFT為做FFT變換時(shí)的規(guī)模），即頻域的分辨率取決于時(shí)域分段時(shí)的長(zhǎng)度，為了比較準(zhǔn)確地在頻域里估計(jì)出模態(tài)參數(shù)，就要求df不能過大，也就要求取段的數(shù)據(jù)長(zhǎng)度nFFT不能太小。在總數(shù)據(jù)量一定的情況下，這會(huì)造成分段數(shù)減小，最后得到的頻響函數(shù)由于平均次數(shù)不夠，會(huì)有很多毛刺，峰值處和模態(tài)耦合嚴(yán)重的地方不夠光滑。為了獲得比較平滑的頻響函數(shù)曲線，需數(shù)據(jù)分段較多，每段的數(shù)據(jù)量減少（nFFT比較?。＿@樣處理會(huì)使頻響曲線比較光滑，較大程度地消除了隨機(jī)誤差，使模態(tài)定階問題等容易解決，但卻帶來一個(gè)問題，在頻響函數(shù)峰值處和模態(tài)耦合比較嚴(yán)重的地方，頻響函數(shù)的數(shù)據(jù)點(diǎn)數(shù)較少，而在頻響函數(shù)幅值比較小的地方相對(duì)較密。本文提出的方法能對(duì)頻響函數(shù)峰值處和模態(tài)耦合密集處的傳遞函數(shù)進(jìn)行加密，使得該區(qū)段的數(shù)據(jù)量增加，而且在該區(qū)段同時(shí)滿足平滑性和高分辨率的特性，從而提高模態(tài)參數(shù)識(shí)別的精度。

2 利用FFT譜連續(xù)細(xì)化加密傳遞函數(shù)

傳遞函數(shù)反映的是結(jié)構(gòu)在頻域的特性：

式中：H(f)為傳遞函數(shù)；F(f)為輸入的復(fù)數(shù)幅值譜；Y(f)為響應(yīng)的復(fù)數(shù)幅值譜。

結(jié)構(gòu)在頻域里的響應(yīng)可以用輸入與傳遞函數(shù)的相乘來表示，而將實(shí)際采樣得到的時(shí)域信號(hào)變換到頻域，要利用到快速傅利葉變換FFT。但實(shí)際采樣得到的信號(hào)會(huì)混入各種噪聲干擾，為了最大程度地排除噪聲的影響，采用輸入輸出的互譜比輸入的自譜得到傳遞函數(shù)的一種估計(jì)H1(f)，以及用輸出的自譜比輸入輸出的互譜得到傳遞函數(shù)的另一種估計(jì)H2(f )等[1]。

傳統(tǒng)的頻域分析方法是利用離散的傅利葉變換，該方法有分辨率的限制且運(yùn)算速度較慢。目前國(guó)內(nèi)外的信號(hào)處理設(shè)備一般都采用快速傅利葉變換（FFT）進(jìn)行頻域分析。傳遞函數(shù)的對(duì)應(yīng)的頻率值依賴于FFT變換。為了解決上述計(jì)算傳遞函數(shù)時(shí)的問題，引入FFT譜連續(xù)細(xì)化分析算法。

根據(jù)文獻(xiàn)[4-6]提出的算法，對(duì)時(shí)間序列x(t)，任意頻率f處的離散傅利葉變換（DFT）系數(shù)：

其中，fs為x(tk)的采樣率；f為指定的頻率值（應(yīng)不超過頻譜分析上限0.5fs）。

利用傅氏級(jí)數(shù)的復(fù)指數(shù)形式可以寫成：

式中：N為x(tk)的數(shù)據(jù)量；c(f)為傅利葉系數(shù)的復(fù)數(shù)表達(dá)。

將上面c(f)的表達(dá)式代到H1(f)的表達(dá)式里，可以得到加密處的傳遞函數(shù)估計(jì)值：

上述3個(gè)式子中：N為響應(yīng)和輸入的分段數(shù)（分段間一般會(huì)有重疊）；fi為指定處的頻率值；yk為響應(yīng)的第k段數(shù)據(jù)；yk(tj)為響應(yīng)的第k段數(shù)據(jù)的第j+1個(gè)數(shù)據(jù)點(diǎn)；fk為輸入的第k段數(shù)據(jù)；fk(tj)為輸入的第k段數(shù)據(jù)的第j+1個(gè)數(shù)據(jù)點(diǎn)；m為分段后，每段的數(shù)據(jù)點(diǎn)數(shù)；Yk為yk的傅利葉變換；Fk為fk的傅利葉變換。

以上以H1(f)估計(jì)為例推導(dǎo)了加密算法，H2(f)估計(jì)的推導(dǎo)過程類似。在算法實(shí)施過程中，由于fi是任意給定的，并不是基頻(1 / T)的整數(shù)倍，m也不一定是2的整次冪，F(xiàn)FT已經(jīng)不能使用了。因此在求H(fi)時(shí)，只能應(yīng)用離散傅利葉變換DFT，其計(jì)算量為N2，在計(jì)算點(diǎn)數(shù)比較多時(shí)，計(jì)算量增加很快。為了加快計(jì)算傳遞函數(shù)時(shí)，可以只在峰值處和頻率耦合比較嚴(yán)重的地方加密。

3 頻域模態(tài)參數(shù)識(shí)別方法

頻域模態(tài)參數(shù)識(shí)別有單模態(tài)識(shí)別法，如峰值拾取法、擬合導(dǎo)納圓法。多模態(tài)識(shí)別法，如最小二乘法、Klosterman迭代識(shí)別法、Levy法和正交多項(xiàng)式擬合法[1]等。本文采用的Levy法對(duì)模態(tài)耦合比較嚴(yán)重的情況有良好的識(shí)別精度，其計(jì)算原理是：在給定的傳遞函數(shù)離散數(shù)據(jù)點(diǎn)H(w)處，擬合有理多項(xiàng)式,確定多項(xiàng)式的系數(shù)，最后由分母多項(xiàng)式的系數(shù)求得模態(tài)頻率和阻尼比，分式的極點(diǎn)值得出留數(shù)，最后由不同測(cè)點(diǎn)的留數(shù)值估計(jì)出振型[1-3]。傳遞函數(shù)的表達(dá)式為：

式中：N為模態(tài)階次；當(dāng)傳遞函數(shù)為加速度傳遞函數(shù)時(shí)，m為2N；當(dāng)傳遞函數(shù)為速度傳遞函數(shù)時(shí)，m為2N-1；當(dāng)傳遞函數(shù)為位移傳遞函數(shù)時(shí)，m為2N-2。由歸一化可以令a2N為1。

很明顯，式（7）為非線性擬合。經(jīng)典的Levy法中誤差向量寫為：

式中：Dk=a0+a1(jwk)+a2(jwk)2+a3(jwk)3…+a2N(jwk)2N；Nk=b0+b1(j wk)+b2(jwk)2+b3(jwk)3…+bm(jwk)m；Hk為理論值；為實(shí)測(cè)的傳遞函數(shù)值。

化非線性為線性誤差為：

則誤差的平方和是：

式中：C為常數(shù)；{θ}為未知系數(shù){a}構(gòu)成的列向量：

根據(jù)E最小的要求，應(yīng)使

由此可以得到：

最后由式（13）求出系數(shù){a}、，由分母的復(fù)根求出頻率和阻尼比，由傳遞函數(shù)在極點(diǎn)處的值得到留數(shù)，最終求得模態(tài)參數(shù)。

4 FFT譜連續(xù)細(xì)化和Levy法聯(lián)合運(yùn)用

聯(lián)合運(yùn)用譜連續(xù)細(xì)化算法和Levy法進(jìn)行模態(tài)識(shí)別，主要步驟如下：（1）將實(shí)測(cè)的輸入輸出數(shù)據(jù)用經(jīng)典的傳遞函數(shù)算法給出估計(jì)。這一步中，可以使分段數(shù)比較大，頻域分辨率比較低，但可以得到比較光滑的傳遞函數(shù)曲線。從而容易地分辨出模態(tài)階次和固有頻率的大致范圍，以及模態(tài)耦合比較嚴(yán)重的區(qū)段；（2）選取需要加密的主模態(tài)區(qū)間以及模態(tài)耦合嚴(yán)重的區(qū)段。在區(qū)段里，設(shè)定出需要達(dá)到的加密后的頻率分辨率，以及加密的點(diǎn)數(shù)。在每個(gè)加密點(diǎn)對(duì)應(yīng)的頻率值處，根據(jù)式（4）所給出的內(nèi)插算法，計(jì)算出該點(diǎn)的傳遞函數(shù)值。則加密后的傳遞函數(shù)值與之前給出的～H 1構(gòu)成了最終的傳遞函數(shù)估計(jì)值～H ；（3）用Levy法算出實(shí)測(cè)頻響曲線～H 所對(duì)應(yīng)的數(shù)學(xué)模型式（8）的系數(shù)值：X0=[b0，b1，b2，…bm；a0，a1，a2，…a2N]；（4）根據(jù)經(jīng)典的Levy法算出模態(tài)參數(shù)；（5）如果要求更精確地求出模態(tài)參數(shù)，可以在步驟（2）中調(diào)整加密區(qū)段和加密后的分辨率，然后運(yùn)用Levy法求出模態(tài)參數(shù)，考察模態(tài)參數(shù)的穩(wěn)定性，確定比較穩(wěn)定可靠的模態(tài)參數(shù)。

從上述算法過程可知加密后的傳遞函數(shù)，其主模態(tài)區(qū)間以及模態(tài)耦合嚴(yán)重的區(qū)段數(shù)據(jù)量大大增加，而且在該區(qū)段同時(shí)滿足了平滑性和高分辨率的特性。由于運(yùn)用Levy法求解模態(tài)時(shí)，著重強(qiáng)調(diào)了主模態(tài)區(qū)間以及模態(tài)耦合嚴(yán)重區(qū)段的擬合誤差中，因此模態(tài)識(shí)別精度得以提高。而傳統(tǒng)Levy法的擬合誤差僅僅強(qiáng)調(diào)了未加密前的傳遞函數(shù)實(shí)測(cè)值處的誤差，這是與本文算法的最大的不同點(diǎn)。

5 算例

本文給出兩個(gè)算例，一個(gè)數(shù)值仿真和一個(gè)實(shí)測(cè)算例來驗(yàn)證本算法的有效性。

5.1 算例1選取一個(gè)四自由度的彈簧質(zhì)點(diǎn)模型（見圖1），用白噪聲激勵(lì)，根據(jù)加速度響應(yīng)進(jìn)行算法驗(yàn)證。在質(zhì)量單元2處施加160s的白噪聲激勵(lì)，并在第一質(zhì)量單元m1的加速度響應(yīng)中混有噪聲干擾。利用質(zhì)量單元1的加速度響應(yīng)求傳遞函數(shù)H12(f)。激振力的采樣頻率和響應(yīng)的采樣頻率均為250Hz。

采用H1估計(jì)傳遞函數(shù)，分段長(zhǎng)度為2 048點(diǎn)，頻率分辨率為：df=250/2048=0.1221Hz，每段重疊了1 024點(diǎn)，每段里施加“hanning”窗。得到輸入點(diǎn)2和響應(yīng)點(diǎn)1間的傳遞函數(shù)H12(f)，如圖2所示。將傳遞函數(shù)H12(f)的6.5～32Hz間的部分的加密，共加密了1 275個(gè)點(diǎn)，加密后的傳遞函數(shù)的分辨率接近0.02Hz，如圖3所示。

應(yīng)用Levy法擬合加密前、后的傳遞函數(shù)H12(f)，如圖4、圖5所示。實(shí)線是根據(jù)輸入、輸出響應(yīng)計(jì)算得到的傳遞函數(shù)H12(f)，點(diǎn)劃線是Levy法擬合H12(f)的值。采用Levy法擬合傳遞函數(shù)H12的結(jié)果如表1所示。

由圖4、圖5及表1可以看出，F(xiàn)FT譜連續(xù)細(xì)化分析和Levy法聯(lián)合運(yùn)用的優(yōu)勢(shì)，其4階模態(tài)都能夠比較準(zhǔn)確地識(shí)別出來，未加密的算法的誤差比較大。上述的加密算法不僅能用在傳遞函的計(jì)算方面，還可以用在自功率譜、互功率譜估計(jì)等方面，能夠給出更加平滑、飽滿的曲線（峰值處半功率帶寬里可以達(dá)到8個(gè)點(diǎn)以上）。

圖2 傳遞函數(shù)H12（f）幅頻

圖3 加密了的傳遞函數(shù)H12（f）幅頻（圓點(diǎn)為內(nèi)插值）

圖4 Levy法擬合未加密的傳遞函數(shù)H12幅頻

圖5 Levy法擬合加密的傳遞函數(shù)H12幅頻

表1 Levy法擬合傳遞函數(shù)H12的結(jié)果

5.2 算例2下面給出了本文的算法在某大型升船機(jī)振動(dòng)臺(tái)模型試驗(yàn)研究[7]中的應(yīng)用。該升船機(jī)承重塔柱長(zhǎng)119m，寬57.8m，高145m，由對(duì)稱布置4個(gè)塔柱構(gòu)成。模型的長(zhǎng)度比尺為1∶25，在塔柱模型上布置了51個(gè)加速度傳感器，模型底板上表面分別沿X、Y、Z方向各布置1個(gè)，模型塔柱四角沿豎向隔層布置1～2個(gè)。

在第一階段小震情工況時(shí)，沿橫河向輸入了50gal的白噪聲激勵(lì)，得到了試驗(yàn)前模型的動(dòng)力特性。選取臺(tái)面順河向的輸入時(shí)程和塔柱上橫河向第23#測(cè)點(diǎn)的響應(yīng)驗(yàn)證本文的算法。

采用H1估計(jì)傳遞函數(shù)，并將傳遞函數(shù)H(f)的1～60Hz間的部分加密，共加密了3 000個(gè)點(diǎn)。加密前的頻率分辨率為0.25Hz，加密后為0.019 7Hz。如圖6、圖7所示。

應(yīng)用Levy法擬合加密前、后的傳遞函數(shù)H(f)，得到前3階模態(tài)參數(shù)，如圖8、圖9和表2所示。

由表2可以很明顯地看出本文提出的細(xì)化算法的優(yōu)勢(shì)。由圖8、圖9可以看出，加密后的傳遞函數(shù)的擬合值與實(shí)測(cè)傳遞函數(shù)值相當(dāng)吻合，而且擬合的傳遞函數(shù)曲線比較平滑，阻尼比的識(shí)別精度有很大的提高，且第四階模態(tài)參數(shù)的識(shí)別結(jié)果也優(yōu)于未加密前的結(jié)果。本文的算法在實(shí)際應(yīng)用中得到了驗(yàn)證。

圖6 未加密的傳遞函數(shù)幅頻

圖7 加密后的傳遞函數(shù)幅頻（圓點(diǎn)內(nèi)插點(diǎn)）

圖8 Levy法擬合未加密的傳遞函數(shù)幅頻

圖9 Levy法擬合加密后的傳遞函數(shù)幅頻

表2 某大型升船機(jī)塔柱模型橫河向前3階模態(tài)參數(shù)

6 結(jié)語(yǔ)

本文提出了一種利用FFT譜連續(xù)細(xì)化算法對(duì)傳遞函數(shù)在峰值處和模態(tài)耦合密集處進(jìn)行加密，然后利用Levy法進(jìn)行模態(tài)參數(shù)識(shí)別的算法。數(shù)值仿真結(jié)果表明，此算法在區(qū)分密集模態(tài)和提高模態(tài)參數(shù)識(shí)別結(jié)果方面有較大的改進(jìn)，而且在抗噪聲方面有很大的優(yōu)勢(shì)。但改進(jìn)算法在加密傳遞函數(shù)時(shí)，由于不能采用快速FFT變換，導(dǎo)致計(jì)算量比較大，為了加快計(jì)算傳遞函數(shù)，可以選擇在峰值處和模態(tài)耦合比較嚴(yán)重的地方加密。應(yīng)用本文算法在某大型升船機(jī)振動(dòng)臺(tái)模型試驗(yàn)中取得了良好的結(jié)果，驗(yàn)證了其實(shí)際應(yīng)用價(jià)值。數(shù)值實(shí)驗(yàn)和實(shí)測(cè)實(shí)驗(yàn)都很好地驗(yàn)證了本文算法的有效性。

[1]李德葆，陸秋海.實(shí)驗(yàn)?zāi)B(tài)分析及其應(yīng)用[M].北京：科學(xué)出版社，2001.

[2]王濟(jì)，胡曉.MATLAB在振動(dòng)信號(hào)處理中的應(yīng)用[M].北京：中國(guó)水利水電出版社，2005.

[3]傅志方.振動(dòng)模態(tài)分析與參數(shù)辯識(shí)用[M].北京：機(jī)械工業(yè)出版社，1990.

[4]劉進(jìn)明，應(yīng)懷樵.FFT譜連續(xù)細(xì)化分析的富里葉變換法[J].振動(dòng)工程學(xué)報(bào)，1995，8（2）：162-166.

[5]孫鶴泉.實(shí)用Fourier變換及C++實(shí)現(xiàn)[M].北京：科學(xué)出版社，2009.

[6]Stefan Hollos，Richard Hollos.C Program Magnifies Spectrum When An FFT Can’t Hack It Highlights[OL].http：//electronicdesign.com/article/embedded/c-program-magnifies-spectrum-when-an-fft-can-t-hac.aspx.August 18，2003.

[7]中國(guó)水利水電科學(xué)研究院.某升船機(jī)塔柱大型振動(dòng)臺(tái)抗震試驗(yàn)研究報(bào)告[R].北京：中國(guó)水利水電科學(xué)研究院，2009.