曹艳君 王 皓

（复旦大学航空航天系，上海 200433）

基于快速Fourier变换法的广义特征值问题重根辨识方法∗

曹艳君 王 皓†

（复旦大学航空航天系，上海 200433）

对具有重根的广义特征值问题，采用基于快速Fourier变换的方法进行求解，实现重根辨识.文章中采用多次单点初始激励的方式，仿真计算测点上的自由振动响应，对响应进行快速Fourier变换后得到频域数据.而后对频域数据分析，得到固有频率和多组测点振型数据.根据单频和重频处的振型特性，引入振型的余弦相似度为判别参数，辨识重根.数值算例表明，该方法可有效实现重根辨识，同时特征值的计算能达到较高精度.

广义特征值问题， 重根辨识， 快速Fourier变换法， 固有频率， 动力学响应

引言

广义特征值问题的求解是结构动力分析和稳定性分析的关键.求解广义特征值问题的传统算法有Lanczos 法、子空间迭代法、Ritz 向量法［1，2］和 Jacobi算法等.近年来也涌现出一些新兴算法，如混合人工鱼群算法［3］、神经网络算法［4，5］、瀑布型多重网格法［6］等.

2014年，吴锋、徐小明和钟万勰提出将快速Fourier变换法用于求解广义特征值问题［7］.该方法的主要思想是，采用计算机仿真计算结构的自由振动响应，对响应进行快速Fourier变换，从频域数据可直观地得到特征值信息.但对于具有重频的广义特征值问题，作者并未叙及.

重频是多自由度系统中固有频率重叠的现象，在动力学分析中较为常见.在求解广义特征值问题过程中若没有考虑重频存在的可能性，求解时会遗漏重频处的振型.这会导致分析结果误差增大，甚至得到完全错误的结果［8］.因此，在广义特征值问题的求解中，有必要对重频进行辨识.

本文主要对具有重频的广义特征值问题进行研究，采用多次单点初始激励，仿真计算测点上的动力学响应.在对频域数据进行分析时，除提取固有频率信息以外，还同时得到其对应的振型.得到固有频率和多组振型数据后，根据单频和重频处振型的不同特性，以振型向量的余弦相似度为判别参数辨识重频.

1 重频结构与仿真计算

1.1 重频结构的自由振动问题

无阻尼系统的自由振动方程为:

其中x为结构的位移向量，M为质量矩阵，K为刚度矩阵，n为系统的总自由度数．

本文仅讨论K为正定矩阵的情形，此时ωi大于零．

若系统固有频率中存在重频，假设ωm为k重频（k＞1），即:

对应的 k 个主振型记为 φm，φm+1，…，φm+k-1．系统的自由振动的解写成以下形式:x（t） ＝

对于非零位移、零速度初始条件，系统自由振动的响应式为:

对于第j个自由度，有响应:

式（6）中，单频 ωs前系数为而重频ωm前系数为在仿真计算中，这些系数可通过模态参数辨识方法得到．

第二次计算 ωs处的振型为＝，而ωm处的振型为:

由此可知，对于具有二重频的问题，单频处两次计算得到的振型归一化后是一致的．而重频处，只要:

得到的两个振型就线性独立．于是，根据此特性可以对重频进行判别．

自由振动响应算法可采用精细积分法［9］、Runge-Kutta法［10］等．实际计算中，无需在所有的自由度上得到完整的振型．只需在结构中选取若干个测点，得到这些测点上的振型，便可识别出重频．

1.2 快速Fourier变换及振型提取

在利用快速Fourier变换进行仿真计算时，要求数据在时域和频率均是离散的，且为有限长［11］．若采样间隔为Δt，采样点数为Nt，采样得到的响应数据序列为 xj（r）（r＝ 0，1，2，…，Nt-1），xj（r）做快速Fourier变换后得到的频域数据为Xj（l），有:

其中 Aj（l）为 Xj（l）的幅值，φj（l）为相位．

对Xj（l）提取峰值便可得到结构的固有频率．同时，分析频域数据可得到振型．此处的振型指在选取的测点上的振型，而非完整振型．

1.3 采样时间间隔Δt和采样点数Nt

仿真计算时，采样时间间隔Δt和采样点数Nt的选取需要满足一定要求．

选取采样间隔Δt之前，要对结构的最大频率作估计．由采样定理可知，采样频率应至少是最大频率的两倍．基于广义盖尔圆定理［12］，有最大特征值:

仿真计算时同样需要考虑频率分辨率．频率分辨率指算法能将信号中两个相近的谱峰保持分开的能力．若信号由两个相近频率组成，分别为ω0和ω0+Δω0．当 Δω0很小时，可能导致两个峰合在一起．此时需要增加采样点数Nt，使得频域数据两峰之间出现明显的波谷．根据文献［7］，采样点数Nt可通过下式确定:

1.4 振型向量相似度度量

两次计算在不同的点上给定初始条件，得到两组频率和测点上振型的数据．满足一定条件时，可望重频处得到的两个振型向量线性无关．此处引入振型向量的余弦相似度刻画振型的相似度．若有两个 n 维向量 μ＝｛μ1，μ2，…，μn｝T和 η＝｛η1，η2，…，ηn｝T，则有μ和η夹角的余弦值为:

cosθ的值越接近于1，两向量的相似度越高．| cosθ＝1时， μ 和 η 可以互相线性表出；cosθ＝0时，μ和η相互正交．

理论上，单频处两个振型向量的夹角余弦值|cosθ|＝1，而重频处，夹角余弦值|cosθ|＜1．由此，以振型向量的余弦相似度为判别参数，可实现重频的判别．

2 算例验证与分析

2.1 算法

以具有二重根的广义特征值问题为例，具体的算法如下:

（1）选定s个测点，确定识别指标h和频率分辨率Δω0的值，由式（13）和式（14）确定时间间隔Δt和采样点数Nt．

（2）取一个初始激励点，给定初始条件．

（3）计算给定初始条件下测点的位移响应，对响应数据进行快速Fourier变换，得到频域数据．

（4）从频域数据中提取固有频率和振型信息，振型做归一化处理．

（5）另取一个初始激励点，给定初始条件，执行步骤（3）和（4），得到另一组频率和振型数据．

（6）计算每一阶频率两振型向量的夹角余弦值，判别重频．

2.2 算例验证

如图1所示，有一四边固支的对称钢板结构，尺寸为 1m×1m×0.002m，弹性模量 E＝2×1011Pa，泊松比 μ＝0．3，密度 ρ＝ 7800kg／m3．采用矩形板单元进行网格划分（板单元的每个节点上有3个自由度，分别为挠度和x方向、y方向的转角）．如图，划分网格后共有8×8＝64个单元，81个节点．节点编号方式为:左下角节点为1号节点，沿左侧边向上编号，至左上角点为9号节点，而后以底边左面第二个节点为10号节点，依次向上编号，以此类推，直至编至右上角点为81号节点．

图1 钢板结构Fig.1 Steel plate

选择5个节点计算响应，节点编号为20、29、51、62和71号．第一次计算时，在22号节点挠度方向给定大小为0．05m初位移，第二次在65号节点挠度方向给定同样大小的初位移．两次计算初速度均为零．

取 Δω0＝0．5，h＝5，相应的采样时间间隔 Δt＝2×10-4s，计算响应的数值方法为四阶Runge-Kutta法．

提取频域数据的峰值得到固有频率，综合多个测点的数据，得到如表1第二列和第四列所示的固有频率结果．

根据式（10）和式（11）从频域数据中提取振型，对于每个测点，其挠度方向计入振型．

其振型向量夹角余弦值为 cosθ1＝1．000．

夹角余弦值cosθ2＝-0．235．由此，可以判定处存在重频．

由于篇幅原因，其余各阶频率处的振型不再一一列出，振型向量的夹角余弦值如表1第三列和第六列所示．

表1 固有频率和振型夹角余弦值Table 1 Natural frequency and cosine values

表2 计算结果与精确解Table 2 Calculated and exact results

由表 1 可知，频率 35.138Hz、76.550Hz、101.663Hz、133.225Hz和160.050Hz是重频.用Matlab的eig函数计算精确解，将计算结果和精确解进行对比.表2中列出了前20阶固有频率的结果（带∗的为重频）.由此可见，本方法能较好地实现重频辨识，同时固有频率计算达到较高的精度.

2.3 分析与讨论

对于具有二重频的问题，以上算法需要进行两次响应计算．理论上，只计算一次响应也可实现重频的辨识．计算响应的初始条件为:在一个点上给定初位移，同时另一个点上给定初速度．此时重频处实部分量得到的振型˜φreal和虚部分量得到的振型˜φimag线性无关．理论上可据此实现重频的判别，但实际操作时，重频辨识结果受初始条件的影响较大，效果不理想．因此，对于具有二重频的广义特征值问题，本文仍选择求两次响应的方式．

若研究具有k重频（k≥2）的问题，首先，可按照二重频的方法，计算固有频率并辨识重频．在每阶重频处，可得两个线性无关的振型．而后继续选用新的初始激励点，计算动力学响应，并得到新的一组振型数据．若重频处的振型不可用之前的振型线性表出，则表示该阶频率处又找到一个新的振型．继续选用新的初始激励点，计算响应，直到重频处的振型均可用之前的振型线性表出．此时，重频处得到的线性无关的振型的数目即是该阶频率的重数．

3 结论

本文采用多次单点初始激励的方式，从响应的频域数据中得到多组测点上的振型，引入振型的余弦相似度为判别参数辨识重根.以具有二重频的问题为例进行特征值计算和重频辨识.计算结果表明，振型向量的余弦值在单频和重频处有明显差异，能够根据余弦值有效辨识重频，同时特征值计算可达到较高精度.

1 Li R C，Zhang L H.Convergence of the block Lanczos method for eigenvalue clusters.Numerische Mathematik，2015，131（1）:83～113

2 宫玉才，周洪伟，陈璞等.快速子空间迭代法、迭代Ritz向量法与迭代Lanczos法的比较.振动工程学报，2005，18（2）:227～232 （Gong Y C， Zhou H W， Chen P， et al.Comparative study of fast subspace iteration method，Ritz method and Lanczos method.Journal of Vibration Engineering， 2005，18（2）:227～232 （in Chinese））

3 黄华娟，周永权，韦杏琼等.求解矩阵特征值的混合人工鱼群算法.计算机工程与应用，2010，46（6）:56～59（Huang H J， Zhou Y Q， Wei X Q， et al.Artificial fish school method for eigenvalue problems.Computer Engineering and Applications， 2010，46（6）:56～59 （in Chinese））

4 Nandy S，Sharma R，Bhattacharyya S P.Solving symmetric eigenvalue problem via genetic algorithms:Serial versus parallel implementation.Applied Soft Computing， 2011，11（5）:3946～3961

5 Tang Y，Li J P.Another neural network based approach for computing eigenvalues and eigenvectors of real skew-symmetric matrices.Computers＆Mathematics with Applications， 2010，60（5）:1385～1392

6 Chen H T，He Y H，Li Y，et al.A multigrid method for eigenvalue problems based on shifted-inverse power technique.European Journal of Mathematics， 2015，1（1）:207～228

7 吴锋，徐小明，钟万勰.广义特征值问题的快速傅里叶变换法.振动与冲击， 2014，33（22）:67～71 （Wu F， Xu X M，Zhong W X.Fast Fourier transform method for generalized eigenvalue problems.Journal of Vibration and Shock，2014，33（22）:67～71 （in Chinese））

8 刘福林，闫维明，何浩祥.重频结构的振型识别及处理方法.工业建筑， 2011，41（S1）:368～372 （Liu F L， Yan W M，He H X.Mode shape identification and disposal in the structure of overlap frequencies.Industrial Construction， 2011，41（S1）:368～372 （in Chinese））

9 张继锋，邓子辰，徐方暖等.一种新的改进精细直接积分法.动力学与控制学报， 2015，13（4）:241～245 （Zhang J F，Deng Z C，Xu F N，et al.A new improved precise direct integration method.Journal of Dynamics and Control，2015，13（4）:241～245 （in Chinese））

10 吴志桥，高普云，任钧国.Runge-Kutta方法求解结构动力学方程.系统仿真学报，2010，22（9）:2085～2090（Wu Z Q，Gao P Y，Ren J G.Runge-Kutta methods for time integration in structural dynamics.Journal of System Simulation， 2010，22（9）:2085～2090 （in Chinese））

11 陈恒亮，蒋勇.基于DSP的实数FFT算法研究与实现.动力学与控制学报， 2005，3（2）:52～55 （Chen H L，Jiang Y.Design and realization of real FFT based on DSP.Journal of Dynamics and Control， 2005，3（2）:52～55 （in Chinese））

12 Nakatsukasa Y.Gerschgorin′s theorem for generalized ei-genvalue problems in the Euclidean metric.Mathematics of Computation， 2011，80（276）:2127～2142

A METHOD OF MULTIPLE-FREQUENCY IDENTIFICATION FOR GENERALIZED EIGENVALUE PROBLEMS BASED ON FAST FOURIER TRANSFORM

Cao Yanjun Wang Hao†
（Department of Aeronautics and Astronautics，Fudan University，Shanghai 200433，China）

A method based on the fast Fourier transform was proposed to solve the generalized eigenvalue problems with multiple roots.This paper studied the dynamic responses of the measure nodes with nonzero initial condition on one point.Both natural frequencies and mode shapes were extracted from the data in the frequency domain.Responses under different initial conditions were calculated to get several sets of mode shapes.Taking the cosine similarity of mode shapes as the discriminant parameter，multiple roots were then identified.The numerical example shows that this method can identify the multiple roots efficiently，and the result reaches a high accuracy.

generalized eigenvalue problem， multi-frequency identification， fast Fourier transform， natural frequency， dynamic response

13 September 2016，revised 21 December 2016.

10.6052／1672-6553-2017-015

2016-09-13收到第1稿，2016-12-21收到修改稿.

∗国家自然科学基金资助项目（11572089）

†通讯作者 E-mail:wanghao＠fudan.edu.cn

∗The project supported by the National Natural Science Foundation of China（11572089）

† Corresponding author E-mail:wanghao＠fudan.edu.cn