吳 鋒，徐小明，鐘萬勰

（大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室，工程力學(xué)系，遼寧大連 116024）

廣義特征值問題的快速傅里葉變換法

吳 鋒，徐小明，鐘萬勰

（大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室，工程力學(xué)系，遼寧大連 116024）

針對廣義特征值問題提出離散傅里葉變換法。該方法把結(jié)構(gòu)的動(dòng)力響應(yīng)看作是一種信號，利用快速傅里葉變換進(jìn)行分析，從而得到結(jié)構(gòu)的振動(dòng)頻率。該方法避免對剛度矩陣求逆，可同時(shí)計(jì)算出所有的特征值，是一種直接方法。數(shù)值算例驗(yàn)證了該方法的正確性。

特征值；動(dòng)態(tài)結(jié)構(gòu)響應(yīng)；快速傅里葉變換；采樣

在工程中，結(jié)構(gòu)振動(dòng)頻率分析和穩(wěn)定分析都涉及到廣義特征值問題。關(guān)于廣義特征值問題已經(jīng)提出許多算法［1－9］，其中常用的算法有子空間迭代法、Ritz向量法、Lanczos算法等［4，6］。這些方法均為迭代方法，通常是利用一組基，把大規(guī)模的廣義特征值問題轉(zhuǎn)化成一個(gè)小規(guī)模的特征值問題，然后再用直接解法求解這一小規(guī)模特征值問題。對于小規(guī)模的特征值問題，通常采用Jacobi算法［2］。一般說來，迭代方法必須結(jié)合直接算法才可以進(jìn)行，而且需要對剛度矩陣求逆，往往只能求解部分特征值。

對于一個(gè)結(jié)構(gòu)，其結(jié)構(gòu)振動(dòng)響應(yīng)中包含著關(guān)于結(jié)構(gòu)振動(dòng)頻率的所有信息，如果把結(jié)構(gòu)的動(dòng)力響應(yīng)看作是一組信號，則可以利用信號處理中的快速傅里葉變換進(jìn)行分析［10，11］，從而找到結(jié)構(gòu)的振動(dòng)頻率。在實(shí)際結(jié)構(gòu)的模態(tài)分析中，常通過對測量得到的結(jié)構(gòu)振動(dòng)響應(yīng)進(jìn)行傅里葉變換，分析其結(jié)構(gòu)的振動(dòng)頻率［12］。實(shí)際上，這一思想也可用于計(jì)算廣義特征值問題。對于一組給定的剛度矩陣和質(zhì)量矩陣，可以人為構(gòu)造一種工況，采用計(jì)算機(jī)仿真計(jì)算其動(dòng)力響應(yīng)，然后利用快速傅里葉變換進(jìn)行分析，從而得到特征值。關(guān)于結(jié)構(gòu)的動(dòng)力響應(yīng)計(jì)算有許多成熟算法［13－16］?；谶@一認(rèn)識，本文針對廣義特征值問題提出基于結(jié)構(gòu)動(dòng)力響應(yīng)的離散傅里葉變換法。該方法把結(jié)構(gòu)的動(dòng)力響應(yīng)看作是一種信號，而結(jié)構(gòu)的動(dòng)力響應(yīng)可以利用成熟的算法如中心差分法、精細(xì)積分法等求解，也可以實(shí)際測量，從而避免對剛度矩陣求逆。得到結(jié)構(gòu)動(dòng)力響應(yīng)后，利用快速傅里葉變換進(jìn)行分析。本文方法可以同時(shí)計(jì)算出所有的特征值，是一種直接方法，也可以與迭代方法結(jié)合求解，處理迭代法中需要解決的小規(guī)模特征值問題。

1 基本思想

本文考慮結(jié)構(gòu)的廣義特征值問題

式中：K，M為結(jié)構(gòu)剛度、質(zhì)量矩陣；x為特征向量；ωi為振動(dòng)頻率。且有

求解式（1）的廣義特征值問題常用方法為子空間迭代、Lanczos迭代等。本文提出與此完全不同之方法，即由結(jié)構(gòu)振動(dòng)響應(yīng)出發(fā)，對結(jié)構(gòu)振動(dòng)響應(yīng)進(jìn)行傅里葉變換，獲得結(jié)構(gòu)振動(dòng)頻率。

1.1 自由振動(dòng)

設(shè)某工況下結(jié)構(gòu)振動(dòng)響應(yīng)只含結(jié)構(gòu)振動(dòng)頻率，不含其它任何雜質(zhì)，即無阻尼系統(tǒng)的自由振動(dòng)。設(shè)結(jié)構(gòu)位移向量a，其動(dòng)力微分方程為

由式（9）可見結(jié)構(gòu)位移響應(yīng)中只含結(jié)構(gòu)振動(dòng)頻率，對結(jié)構(gòu)位移響應(yīng)進(jìn)行傅里葉變換便可獲得結(jié)構(gòu)的振動(dòng)頻率。

1.2 時(shí)間步長及采樣長度

以上分析均基于連續(xù)情況，實(shí)際計(jì)算結(jié)構(gòu)動(dòng)力響應(yīng)時(shí)只能得到各時(shí)間點(diǎn)位移值，故可用快速離散傅里葉變換（FFT）。計(jì)算結(jié)構(gòu)動(dòng)力響應(yīng)應(yīng)選精度較好的動(dòng)力算法。對規(guī)模較大動(dòng)力問題需避免求逆。獲得結(jié)構(gòu)在各時(shí)間點(diǎn)位移值后須用FFT進(jìn)行分析，其中時(shí)間步長取決于位移中所含最大頻率ωmax，而采樣長度則取決于位移中兩相鄰頻率之差Δω。設(shè)需分析信號為

根據(jù)式（15）知，如果頻率w為整數(shù)時(shí)，am＝δmw。因此根據(jù)am是否等于1，就可以找到頻率w。但是因?yàn)閙≤0．5N，如果w＞0．5N，則w≠m，此時(shí)對于所有的m有am＝0，無法識別出信號的頻率，而其截?cái)郚項(xiàng)的傅里葉級數(shù)為0，傅里葉信號失真。因此僅當(dāng)w≤0．5N時(shí)才可以識別出信號的頻率，而其離散傅里葉變換在各個(gè)時(shí)間點(diǎn)上的值才不會(huì)失真。把w≤0．5N代入式（14）可得：

am是關(guān)于w的一個(gè)復(fù)函數(shù)，其模為：

因此，｜am（w）｜是一個(gè)關(guān)于w和m的函數(shù)，當(dāng)固定w時(shí)，m越接近w，｜am｜越大，而遠(yuǎn)離m時(shí)，｜am｜只有N－1的量級。現(xiàn)在假設(shè)一個(gè)信號有兩個(gè)頻率組成，一個(gè)為ω，另一個(gè)為ω＋Δω，該信號離散傅里葉變換后的各個(gè)頻率的振幅即為：

只要函數(shù)｜bm（w）｜出現(xiàn)兩個(gè)峰值，一個(gè)在w處，一個(gè)在w＋Δw處即可以識別出這兩個(gè)峰值。當(dāng)Δw很小時(shí)，由于w與w＋Δw相處太近，導(dǎo)致兩個(gè)峰值合成一個(gè)峰值，識別不出這兩個(gè)頻率。這是頻率分辨率的問題，解決的方法是增加采樣長度N，使得在［w］（［w］表示距離w最近的整數(shù)）和［w＋Δw］之間出現(xiàn)明顯的波谷，從而把兩個(gè)峰值明顯的隔開?，F(xiàn)在假設(shè)和［w＋Δw］的平均數(shù)的整數(shù)，要求為一種識別指標(biāo)），因?yàn)?/p>

進(jìn)一步，上式也可以寫成ΔωΔt≥8h／N，這一不等式就是測不準(zhǔn)原理。當(dāng)采樣點(diǎn)數(shù)N固定后，ΔωΔt大于一定值，因此，時(shí)域的分辨率Δt越小，則頻域的分辨率Δω就大，反之亦然，要想兩者都小，則必須增加采樣點(diǎn)數(shù)。

1.3 最大頻率

確定時(shí)間步長時(shí)需要用到結(jié)構(gòu)的最大振動(dòng)頻率，這一最大振動(dòng)頻率可以根據(jù)結(jié)構(gòu)的最小單元的最大特征值確定。當(dāng)結(jié)構(gòu)的質(zhì)量矩陣是集中質(zhì)量陣時(shí)，這里基于廣義蓋爾圓定理［17］，給出一個(gè)確定最大頻率的方法。在標(biāo)準(zhǔn)特征值問題中有個(gè)蓋爾圓定理，那么對于如式（1）所示的廣義特征值問題，我們可以很容易的建立一個(gè)廣義蓋爾圓定理，現(xiàn)在設(shè)：

剛度矩陣K中第i行第j列的元素為kij，質(zhì)量矩陣M為對角矩陣，第i個(gè)對角元素為mi，x中第i個(gè)元素為xi，則矩陣形式的特征方程可以展開成如下所示的N個(gè)方程：

2 算法設(shè)計(jì)

根據(jù)式（38）和式（31）確定時(shí)間步長Δt和采樣長度N。然后根據(jù)式（3）和式（4）計(jì)算結(jié)構(gòu)的響應(yīng)，把響應(yīng)進(jìn)行離散傅里葉變換，既可以找到結(jié)構(gòu)的各階振動(dòng)頻率。在式（38）中，h是一種識別指標(biāo)，描述了兩個(gè)相鄰頻率w和w＋Δw的分辨率。當(dāng)h取得大時(shí)，｜bm（w）｜在這兩個(gè)頻率之間會(huì)出現(xiàn)一個(gè)明顯的波谷，從而可以識別出這兩個(gè)頻率。但是h取值很大時(shí)，采樣長度就會(huì)變大。經(jīng)過計(jì)算發(fā)現(xiàn)，取h＝5較為合理。當(dāng)通過傅里葉變換找到各個(gè)峰值時(shí)，此峰值所對應(yīng)的頻率w為無量綱頻率，還需要通過式（14）化成圓頻率，也就是結(jié)構(gòu)的振動(dòng)頻率。下面講述具體計(jì)算過程。

（1）根據(jù)式（38）和式（31）確定時(shí)間步長和采樣長度N；

（2）選定初始位移M－1a0，根據(jù)式（3）計(jì)算結(jié)構(gòu)的自由振動(dòng)響應(yīng)；

（3）計(jì)算出結(jié)構(gòu)的動(dòng)力響應(yīng)之后，任意提取某節(jié)點(diǎn)的位移響應(yīng)作為信號，對該信號進(jìn)行快速傅里葉變換，得到頻域信號；

（4）對頻域信號取絕對值，找出頻域信號峰值所在的無量綱頻率w；

（5）根據(jù)式（14），把w化成圓頻率，即得到結(jié)構(gòu)的振動(dòng)頻率。

3 算例

計(jì)算圖1所示20層平面鋼框架結(jié)構(gòu)的振動(dòng)頻率。已知框架柱為450×450×25 mm箱型截面，彈性模量E＝2．06× 105MPa。按剛性樓板假設(shè)，采用層剪切模型，集中質(zhì)量，每層的質(zhì)量m＝168 750 kg。分別采用傳統(tǒng)求解特征值問題的QR法、Jacobi方法和本文方法計(jì)算。QR法和Jacobi方法屬于迭代法，把QR法相對誤差取10－13時(shí)的解作為標(biāo)準(zhǔn)解。本文方法計(jì)算時(shí)，初始位移中a0的各項(xiàng)元素均取為1，識別指標(biāo)取h＝5，計(jì)算所用的時(shí)間步長按式（38）選取，而頻率分辨率取為Δω＝0．5。結(jié)構(gòu)的位移響應(yīng)分別采用精細(xì)積分法［11］、中心差分法和文獻(xiàn)［12］中所提出的的高精度攝動(dòng)法計(jì)算。經(jīng)過計(jì)算發(fā)現(xiàn)，把任一層的位移響應(yīng)作為時(shí)域信號進(jìn)行快速傅里葉變換，計(jì)算結(jié)果完全一樣，這是因?yàn)槠鋾r(shí)域信號中所包含的頻率分量完全相同，這里只給出采用第一層位移響應(yīng)作為時(shí)域信號分析的結(jié)果。所有的計(jì)算結(jié)果列于表1。

圖1 平面框架Fig．1 Plane frame

表1 計(jì)算的頻率及相對誤差Tab.1 Com puted frequencies and relative errors

由表1可見，不管采用哪種方法計(jì)算結(jié)構(gòu)的自由振動(dòng)響應(yīng)，計(jì)算得到的頻率完全一樣，這說明當(dāng)時(shí)間步長和采用長度選定后，結(jié)構(gòu)的自由振動(dòng)響應(yīng)計(jì)算手段對計(jì)算結(jié)果的影響不大。與標(biāo)準(zhǔn)解相比，各階頻率的誤差均很小。誤差最大的為第1個(gè)頻率，也只有0．42%，這表明本文方法的可靠性。為進(jìn)一步比較本文方法和經(jīng)典的QR法和Jacobi方法的性能，分別取QR法和Jacobi方法的迭代相對誤差為0．4%，比較計(jì)算時(shí)間，結(jié)果列于表2。表2中，本文方法計(jì)算結(jié)構(gòu)動(dòng)力響應(yīng)時(shí)采用的是精細(xì)積分方法。根據(jù)表2可見，對于小規(guī)模的廣義特征值問題，當(dāng)相對誤差均為0．4%時(shí)，本文方法僅需要很少的計(jì)算時(shí)間即可以給出十分滿意的結(jié)果。這說明本文方法比較適用于小規(guī)模廣義特征值問題的求解。根據(jù)表1的計(jì)算結(jié)果還可見，本文方法計(jì)算的高階頻率的相對誤差要小于低階頻率的相對誤差，這一特點(diǎn)不同于傳統(tǒng)的迭代方法，在傳統(tǒng)迭代方法中，往往高階頻率的相對誤差大于低階頻率的相對誤差。

表2 計(jì)算時(shí)間Tab.1 Computation times

4 結(jié) 論

本文針對廣義特征值問題提出基于結(jié)構(gòu)動(dòng)力響應(yīng)的離散傅里葉變換法。該方法把結(jié)構(gòu)的動(dòng)力響應(yīng)看作是一種信號，利用快速傅里葉變換進(jìn)行分析，從而得到結(jié)構(gòu)的振動(dòng)頻率。在本文方法中，結(jié)構(gòu)的動(dòng)力響應(yīng)可以利用成熟的算法如中心差分法、精細(xì)積分法等求解，也可以實(shí)際測量，從而避免對剛度矩陣求逆。本文方法可以同時(shí)計(jì)算出所有的特征值，是一種直接方法，也可以與迭代方法結(jié)合求解，處理迭代法中需要解決的小規(guī)模特征值問題。必須指出的是，本文方法是把工程中模態(tài)分析的做法進(jìn)行數(shù)值化，從而用于計(jì)算廣義特征值問題。對于一個(gè)復(fù)雜的真實(shí)結(jié)構(gòu)，其動(dòng)力響應(yīng)計(jì)算往往要比動(dòng)力特性計(jì)算費(fèi)時(shí)，從這一點(diǎn)上講，本文方法比較適用于小規(guī)模廣義特征值問題。但是從兩點(diǎn)上講本文方法仍然是有價(jià)值的。

（1）傳統(tǒng)的特征值計(jì)算方法如QR法、Jacobi方法、子空間迭代法等均為迭代方法，其收斂的理論依據(jù)均為反冪法，而本文方法則是基于傅里葉變換理論，因此本文方法是對現(xiàn)有的廣義特征值問題的數(shù)值算法的一種補(bǔ)充，從數(shù)值算法角度上講，本文方法是有價(jià)值的。

（2）本文算法把廣義特征值問題看成是結(jié)構(gòu)動(dòng)力響應(yīng)問題和信號分析問題來處理，具有鮮明的工程背景。數(shù)值算例表明，對于小規(guī)模的廣義特征值問題，本文方法僅需要很少的計(jì)算時(shí)間即可以給出十分滿意的結(jié)果。在本文方法的求解過程中，結(jié)構(gòu)動(dòng)力響應(yīng)的計(jì)算手段對特征值的計(jì)算影響不大。

［1］薛長峰，周樹荃．非對稱廣義特征值問題的并行連續(xù)同倫算法［J］．計(jì)算物理，1997，14（4／5）：619－621．

XUE Chang-feng，ZHOU Shu-quan．Parallel homotopy continuation algorithm for nonsymmetric generalized eigenvalue problem［J］．Chinese Journal of Computational Physics，1997，14（4／5）：619－621．

［2］王順緒，戴華．廣義特征值問題的并行塊Jacobi-Davidson方法及應(yīng)用［J］．計(jì)算力學(xué)學(xué)報(bào)，2008，25（4）：428－433．

WANG Shun-xu，DAI Hua．Parallel block Jacobi-Davidson method for solving large generalized eigenvalue problems and it's application［J］．Chinese Journal of Computational Mechanics，2008，25（4）：428－433．

［3］劉寒冰，龔國慶，劉建設(shè)．一種有效的廣義特征值分析方法［J］．固體力學(xué)學(xué)報(bào)，2003，24（4）：419－428．

LIU Han-bing，GONG Guo-qing，LIU Jian-she．An effective algorithm for generalized eigenvalue problems［J］．ACTA Mechanica Solida Sinica，2003，24（4）：419－428．

［4］李秀梅，吳鋒．廣義特征值問題求解的改進(jìn)Ritz向量法［J］．振動(dòng)與沖擊，2012，31（7）：19－23．

LI Xiu-mei，WU Feng．Improved Ritz vector method for generalized eigenvalue problems［J］．Journal of Vibration and Shock，2012，31（7）：19－23．

［5］李海濱，黃洪鐘，趙明揚(yáng)．有限元結(jié)構(gòu)動(dòng)力分析的廣義特征值的神經(jīng)計(jì)算［J］．哈爾濱工業(yè)大學(xué)學(xué)報(bào)，2006，38（9）：1523－1527．

LI Hai-bin，HUANG Hong-zhong，ZHAO Ming-yang．Finite element approach for structural dynamic analysis using generalized eigenvalue-based neurocomputing［J］．Journal of Harbin Iistitute of Technology，2006，38（9）：1523－1527．

［6］黃吉鋒．求解廣義特征值問題的多重Ritz向量法［J］．力學(xué)學(xué)報(bào)，1999，31（5）：585－595．

HUANG Ji-feng．The multi-ritz-vector method in generalized eigenvalue problems［J］．Chinese Journal of Theoretical and Applied Mechanics，1999，31（5）：585－595．

［7］黃華娟，周永權(quán)，韋杏瓊，等．求解矩陣特征值的混合人工魚群算法［J］．計(jì)算機(jī)工程與應(yīng)用，2010，46（6）：56－59．

HUANG Hua-juan，ZHOU Yong-quan，WEIXing-qiong，et al．Hybrid artificial fish school algorithm for solving matrix eigenvalues［J］．Computer Engineering and Applications，2010，46（6）：56－59．

［8］Nandy S，Sharma R，Bhattacharyya SP．Solving symmetric eigenvalue problem via genetic algorithms：serial versus parallel implementation［J］．Applied Soft Computing，2011，11（5）：3946－3961．

［9］Li H，Yang Y．The adaptive finite elementmethod based on multi-scale discretizations for eigenvalue problems［J］．Computers＆Mathematics with Applications，2013，65（7）：1086－1102．

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

LIU Jin-ming，YING Huai-qiao．Zoom FFT spectrum by Fourier transform［J］．Journal of Vibration Engineering，1995，8（2）：162－166．

［11］易立強(qiáng)，鄺繼順．一種基于FFT的實(shí)時(shí)諧波分析算法［J］．電力系統(tǒng)及其自動(dòng)化學(xué)報(bào)，2007，19（2）：98－102．

YILi-qiang，KUANG Ji-shun．FFT-based algorithm for realtime harmonic analysis［J］．Proceedings of The CSU－EPSA，2007，19（2）：98－102．

［12］Schwarz B，Richardson M．Experimentalmodal analysis［J］．CSIReliability Week，1999，36（6）：30－32．

［13］周杰梁，鄭百林．基于中心差分法的纖維結(jié)構(gòu)碰撞動(dòng)力學(xué)分析［J］．力學(xué)季刊，2011，32（3）：466－472．

ZHOU Jie-liang，ZHENG Bai-lin．Impact dynamics analysisof fiber-composed structure based on central difference method［J］．Chinese Quarterly of Mechanics，2011，32（3）：466－472．

［14］鐘萬勰．子域精細(xì)積分及偏微分方程數(shù)值解［J］．計(jì)算結(jié)構(gòu)力學(xué)及其應(yīng)用，1995，12（3）：253－260．

ZHONG Wan-xie．Subdomain precise integration method and numerical solution of partial differential equations［J］．Computational Structural Mechanics and Applications，1995，12（3）：253－260．［15］李秀梅，吳鋒，張克實(shí)．結(jié)構(gòu)動(dòng)力微分方程的一種高精度攝動(dòng)解［J］．工程力學(xué)，2013，30（5）：8－12．

LI Xiu-mei，WU Feng，ZHANG Ke-shi．A high-precision perturbation solution for structural dynam ic differential equations［J］．Engineering Mechanics，2013，30（5）：8－12．

［16］李秀梅，吳鋒，黃哲華．結(jié)構(gòu)動(dòng)力方程一種新的級數(shù)形式的解析解［J］．振動(dòng)與沖擊，2010，29（6）：219－222．

LIXiu-mei，WU Feng，HUANG Zhe-hua．New series form of analytical solution for the structural dynamic equations［J］．Journal of Vibration and Shock，2010，29（6）：219－222．

［17］李秀梅，吳鋒，黃哲華．PCG法的理論解釋及在結(jié)構(gòu)分析中的應(yīng)用［J］．廣西大學(xué)學(xué)報(bào)（自然科學(xué)版），2009（4）：463－468．

LI Xiu-mei，WU Feng，HUANG Zhe-hua．Theoretical explanation for PCG method and its application in structural analysis［J］．Journal of Guangxi University（Natural Science Edition），2009（4）：463－468．

Fast Fourier transform method for generalized eigenvalue problems

WU Feng，XU Xiao-ming，ZHONGWan-xie
（State Key Laboratory of Structural Analysis of Industrial Equipment，Department of Engineering Mechanics，Dalian University of Technology，Dalian 116024，China）

For generalized eigenvalue problems，a method based on the fast Fourier transform（FFT）was developed．In the proposed method，the dynamical structural response was viewed as a signal containing all information about the vibrational frequencies．Using FFT to the signal，the vibrational frequencies can be obtained．The method is a kind of direct solution method which can compute all eigenvalues without the matrix inversion．A numerical example manifests the correctness of the proposed method．

eigenvalue；dynamical structural response；fast Fourier transform；sampling

O327；O214．6

：A

10．13465／j．cnki．jvs．2014．22．012

國家基礎(chǔ)研究發(fā)展計(jì)劃973（2009CB918501）

2013－08－26 修改稿收到日期：2013－11－15

吳鋒男，博士生，1985年生