蒲黔辉， 洪 彧， 王高新， 李晓斌

(1.西南交通大学 土木工程学院，成都 610031； 2.东南大学 土木工程学院，南京 210096)

试验模态参数识别仍是目前研究结构动力学特性的主要方法，识别的结构模态参数可以用于结构有限元模型修正、损伤识别和结构的实时安全监测等。然而土木结构往往都具有结构形式复杂、体型巨大、周围环境复杂多变等特点，例如高楼和桥梁，如果运用传统的锤击或者激振器方式来获得结构强迫振动往往效果不佳，可能仅有激励点附近的布设的传感器才能采集到结构的响应信号，不利于进行整体结构的模态分析。而如果采用大能量的激励会增加试验成本，甚至导致结构产生永久性的损伤，所以一般情况下运营中的土木结构不允许采用此种模态试验方法。此外，运营中的土木结构，特别是桥梁，一般不能长时间完全封锁路段，而车辆行驶造成结构振动会对模态测试造成误差。因此，基于环境激励的模态测试技术成为了土木结构中常用的模态试验方法。

基于环境激励的模态测试是指采集结构在地脉动、风荷载、移动车辆荷载等未知荷载作用下产生的振动响应，并提取响应中的自由衰减成分进行模态参数识别。这种方法不需要专门的激励设备，也不影响结构的正常运营，是一种操作简单、成本低、安全性高的测试方法。James等[1-2]提出了自然激励技术(Natural Excitation Technique, NExT)，证明了结构上任意两个自由度上的位移响应之间的互相关函数满足结构自由振动方程。Ibrahim等[3-4]提出了ITD(Ibrahim Time Domain)算法，该算法利用自由衰减响应信号进行三次不同延时采样，构造自由响应数据的增广矩阵，并由响应与特征值之间的复指数关系，建立特征矩阵，求解特征值和特征向量，再根据模型特征值与振动系统特征值的关系，求解出系统的模态参数。但该方法只适用于单自由度激励情况。后来，Juang等[5-6]提出了特征系统实现算法(Eigensystem Realization Algorithm, ERA)，该算法是在系统最小实现理论的基础上发展起来的时域方法，可适用于多输入多输出(Multiple Input Multiple Ouput, MIMO)结构系统。该方法需要利用脉冲响应构建Hankel矩阵并对其进行奇异值分解，再利用分解结果构建系统矩阵，从而实现模态参数识别。祁泉泉等[7]提出了扩展特征实现算法(Extended Eigensystem Realization Algorithm, EERA)，把传统的ERA算法使用范围扩展至了任何随机信号。崔定宇等[8-9]提出了一种改进的ERA算法，该方法利用奇异值截断后重建Hankel矩阵来对信号降噪，使得结构的高阶模态更容易识别。包兴先等[10]也改进了ERA方法，提出了利用低秩Hankel矩阵逼近方法对响应信号降噪，通过降噪信号后进行模态参数识别。目前采用ERA来进行模态参数识别方法已被国内外的学者广泛使用[11-12]。

由于土木结构复杂几何形式，往往需要较多的传感器采集结构振动信号才能得到完整的模态信息。并且对于环境激励而言，充足的采样时间才能保证结构模态分析的可靠性[13]。因此，采用ERA算法构建的Hankel矩阵往往尺寸巨大，对于大型矩阵进行奇异值分解需要消耗大量的运算和储存时间。针对这一问题，本文提出的快速特征系统实现算法，采用特征值分解来代替奇异值分解，仅利用非零的特征值和对应的特征向量来表达系统矩阵，从而实现模态参数识别。本文采用了一个四层框架仿真模拟来验证该方法的实用性，并且把该方法运用到了一座人行桥的模态参数识别上。

1 环境激励下的结构模态参数识别理论

1.1 自然激励技术

对于一个具有n个自由度的线性结构，在时刻t的振动微分方程为

(1)

假设结构的输入外力和响应都是平稳的随机过程，在式(1)左右两端同时乘上s时刻自由度i上的位移响应qi(s)，并对等式的左右两端取期望，整理可得

(2)

式中：τ=t-s>0，Rqqi(τ)为结构所有自由度的位移响应与i自由度位移响应的互相关函数。由式(2)可以看出结构位移响应之间的互相关函数是与结构自由振动有类似的表达公式，James等也详细地证明了利用采集的结构位移响应之间的互相关函数可以进行模态参数识别。同样的，按照类似的推导方式也可以得出：加速度响应之间的互相关函数也可以进行模态参数识别。

1.2 传统特征系统实现算法

为了避免处理二阶求导问题，式(1)可以在状态空间里表示为一阶微分方程。而试验模态分析是要利用结构振动响应信号反推其模态参数，但是结构响应是传感器以一定的采样频率间断采集得到的。所以式(1)需要在状态空间中表示为离散的线性时不变系统，如下所示

x[k+1]=Adx[k]+Bdp[k]

(3a)

y[k]=Cdx[k]+Ddp[k]

(3b)

式中：Ad∈R2n×2n,Bd∈R2n×np,Cd∈Rny×2n,Dd∈Rny×np分别为离散系统的系统矩阵、控制矩阵、观测矩阵和输入观测矩阵；x∈R2n×1为状态空间变量；y∈Rny×1为系统输出，ny为结构中传感器的个数，系统输出可以为加速度、速度和位移，取决于观测矩阵和输入观测矩阵的改变；p∈Rnp×1为系统输入，np为输入的个数。

假设结构没有初始位移和速度，即x[0]=0，把式(3a)代入式(3b)，经过多次迭代可得

(4)

定义Hankel矩阵

(5)

其中，

(6)

H[k-1],H[k]∈Rnyr×nps且s,r≥2n；式(6)是将频响函数矩阵表达式代入式(5)整理后得到。刘宇飞等提出精确的模态分析要求Hankel矩阵尺寸需要足够大，至少要包括信号振动的几个周期。

当k=0时，有H[0]=ZW，对H[0]进行奇异值分解

H[0]=U∑VT

(7)

定义H[0]的广义逆矩阵为H[0]+

H[0]+=V∑-1UT

(8)

为了找到Hankel矩阵与频响函数(包含结构系统矩阵)之间的关系，定义Ey=[I0 … 0]T∈Rnyr×ny，Ep=[I0 … 0]T∈Rnps×np，它们有以下关系存在

(9)

把式(5)～式(8)代入式(9)的左边进行计算分解，可以得到Ad，Bd，Cd的表达式

(10)

1.3 快速特征系统实现算法

通常情况下，由于H[0]矩阵的维数较大，对其进行奇异值分解需要占用计算机相当部分内存，因此本文提出了一种快速特征系统实现算法。利用特征值分解(一种更经济的计算方法)来代替奇异值分解，仅采用非零的特征值与相应的特征向量来构建系统矩阵，相比传统的ERA算法，既提高了运算速度，又减少了数据储存。

定义一个新的矩阵S=H[0]TH[0]，把式(7)代入S后整理可得

S=V∑T∑VT

(11)

在式(11)的左部分右乘V，又由左右奇异向量矩阵的正交性可知，VTV=I，UTU=I，经过整理后可得

SV=V( ∑T∑)=VΛ

(12)

式(12)为S矩阵的特征方程，其中Λ为对角矩阵，对角线上的元素为S矩阵的特征值。由式(12)可以看出，S矩阵的特征值为Hankel矩阵H[0]特征值的平方，并且S矩阵的特征向量矩阵等同于H[0]的右奇异向量矩阵。

矩阵H[0]与矩阵H[1]存在以下关系，定义传递矩阵T

H[1]=H[0]T

(13)

把式(7)、式(12)和式(13)代入式(10)便可利用S矩阵特征分解的结果重新表达系统矩阵Ad，控制矩阵Bd，观测矩阵Cd，

(14)

值得注意：式(13)中的传递矩阵T不能直接按照广义求逆的方法得到(H[0]和H[1]都有可能不是满秩矩阵)。把T矩阵分块表示

(15)

(16)

T2=VΛ-1VTH[0]T

(17)

得到了离散系统的最小实现矩阵Ad后，对其进行特征值分解

AdΦd=ΦdΔd

(18)

式中：Δd=diag(σd1σd2…σd2n)∈C2n×2n为特征值矩阵；Φd=[Φd1Φd2…Φd2n]∈C2n×2n为特征向量矩阵。

根据离散系统与连续时间系统之间的关系，可以求得结构的复模态参数。

(19)

(20)

式中：“abs”为绝对值；“real”为实部。

1.4 快速特征系统实现算法用于环境激励下的结构模态参数识别的步骤

步骤1 对采集信号进行滤波，降噪等预处理；

步骤2 选定参考点，并计算其它响应与参考点的互相关函数；

步骤3 利用互相关函数构建Hankel矩阵H[0]与对称矩阵S；

步骤4 对S矩阵进行特征值分解求得特征值对角阵Λ和特征向量矩阵V；

步骤5 求解T2矩阵，从而求得T矩阵；

步骤6 组建离散系统的系统矩阵Ad，控制矩阵Bd，观测矩阵Cd，并对Ad进行特征值分解；

步骤7 把离散系统的模态参数转换为连续系统(结构)的模态参数。

2 数值仿真

采用文献[14]中的四层框架模型，仿真其在缩小比例的El Centro地震加速度下的各种结构振动，框架模型见图1。仿真模型每层楼板质量均为4.99 kg；层间刚度k1=1 401.1 N/M、k2=1 576.1 N/M、k3=1 225.9 N/m、k4=1 050.8 N/m，假设每层有黏滞阻尼：c1=1.40 Ns/m、c2=1.31 Ns/m、c3=1.49 Ns/m、c4=1.23 Ns/m，其它结构信息参可参见该文。

图1 用于仿真的框架模型Fig.1 Simulated frame model

假设在框架每层楼板中心位置均安装可测得沿图中地震方向的位移计、速度计和加速度计，采样时间为40 s，采样频率为100 Hz。地震加速度与仿真得到的#4测点的位移、速度、加速度响应如图2所示。为了避开模态节点，选择#1测点为参考点，可以得到#1测点和其它点间的互相关函数，现以#1测点与#4测点的互相关函数为例，见图3。

对于自然激励试验而言，信号时长的选用可直接根据互相关函数图像考虑，选用时长应尽量包含振动衰减的完整过程，在本次仿真算例中选择时长0～8 s的数据。利用互相关函数组建Hankel矩阵，再按照1.4节中步骤进行模态参数识别，结果见表1。表1中的理论值是按照质量、刚度、阻尼矩阵直接求解得出的模态参数精确值。从表1可以看出，对于该四层结构，FERA算法能采用位移、速度、和加速度进行模态参数识别(利用相同的MATLAB程序，只改变响应的类型)，且均能得到较好结果，可计算出频率的最大误差仅为0.5%；阻尼比误差相对较大，最大误差约为25%；模态振型误差几乎为0，MAC值全为1，MAC值算法见文献[15]。与传统的ERA算法相比，FERA算法的精度相当，并均显示出利用加速度响应进行模态参数识别结果与理论值最为接近。仅以此小结构为例，FERA算法就较ERA算法速度提高了约一半。

图2 地震激励和#4测点结构响应Fig.2 Earthquake excitation and structural responses on #4

图3 #1测点与#4测点的互相关函数Fig.3 Cross correlation functions of responses on #1 and #4

模态1模态2模态3模态4时间/s理论值a0．92262．47793．84404．9461b0．00290．00870．01330．0152c1．00001．00001．00001．0000快速特征实现算法位移a0．92652．47393．83914．9509b0．00230．00860．01020．01260．28601c1．00001．00001．00001．0000速度a0．92792．47363．83964．9509b0．00240．00890．00910．01150．28292c1．00001．00001．00001．0000加速度a0．92332．47413．84094．9508b0．00330．00910．01230．01210．27863c1．00001．00001．00001．0000特征实现算法位移a0．92652．47393．83914．9511b0．00230．00860．01020．01200．38067c1．00001．00001．00001．0000速度a0．92812．47343．84294．9545b0．00260．0090．00970．01160．41137c1．00001．00001．00001．0000加速度a0．92392．47243．84064．9551b0．00310．00980．01110．01200．41140c1．00001．00001．00001．0000注：a为频率/Hz；b为阻尼比；c为MAC值

3 人行桥环境激励试验

为了验证FERA算法可以用于实际结构的模态参数识别，本文对美国佐治亚理工校园内的一座人行桥进行了模态参数识别，如图4所示。该桥为简支梁桥，长30.18 m、宽2.13 m、高2.74 m，其它结构信息、传感器信息等可参见文献[16-17]。该桥测试数据来自佐治亚理工智能结构实验室，为了获得较为完整的结构模态振型，在桥的几何节点处都布置了加速度测点(共44个)，采集每个测点的竖向和横向加速度，传感器布置参见图5。为了解决测点数量多而传感器数量受限制的问题，对该桥进行了测试分段。首先在测点1～14安装双向无线加速度传感器，采集自然激励下的加速度响应，然后把加速度传感器移动至测点11～26，再采集加速度响应，然后按照相同方法依次采集测点23～34和31～44的加速度响应。每次环境激励测试的采样时间为150 s，采样频率为100 Hz。以#8测点为例，加速度响应如图6所示。

图4 实验人行桥Fig.4 Experimental pedestrian bridge

图5 传感器布置Fig.5 The allocation of accelerometers

图6 #8测点横向加速度Fig.6 Transverse accelerations on #8

测试结束后需要对各个分段结构单独进行模态参数识别，然后利用平均的方法得出整体桥梁的频率和阻尼比，并且利用共节点拼接分段模态振型来得到结构的完整模态振型。以第一分段为例，选择#8测点横向自由度为参考点，计算其他所有自由度上加速度响应与该自由度加速度响应的互相关函数，示例如图7所示。利用互相关函数组建Hankel矩阵，再按照1.4节中步骤对人行桥进行模态参数识别。最终识别出前4阶振型，结果如图8所示，第1阶振型为横向模态、第2阶振型为竖向模态、第3阶振型为横向模态、第4阶振型为扭转模态。为了与模态参数识别结果对比，利用ANSYS对人行桥进行了仿真模态分析，结果显示试验模态与理论模态的频率值有较小差距，但是二者的

模态振型能很好对应，说明FERA算法同样适用于真实结构的模态参数识别(见图9)。本研究同样使用NExT+ERA方法对该桥进行了模态参数识别，与NExT+FERA方法的识别结果很接近，所以不再赘述。但总体上看，本文提出的FERA算法较ERA方法快2～3倍。

图7 #8测点横向与竖向加速度的互相关函数Fig.7 Cross correlation functions of the transverseand vertical accelerations on #8

图8 模态参数识别结果Fig.8 Modal identification results

图9 ANSYS模态参数识别结果Fig.9 Modal identification results from ANSYS analysis

4 结 论

本文提出的FERA模态参数识别算法是对传统ERA算法的一种改进，利用特征值分解代替了复杂的奇异值分解。本文通过四层框架的数值仿真试验和人行桥的自然激励试验验证了FERA算法的准确性和优势。有如下结论：

(1)四层框架数值模拟表明FERA算法较ERA算法有着相同的精度，适用于各种动力响应。但较其他两种响应而言，利用加速度响应进行模态参数分析能得到更精确的结果，所以推荐使用加速度响应进行模态参数识别。

(2)人行桥的自然激励试验表明了FERA算法同样适用于真实结构，并且速度统计显示，结构越大，FERA算法较ERA算法的速度优势更明显。

(3)FERA算法可以提高结构安全监测、结构模型修正、损伤识别的效率。

[ 1 ] JAMES G H, CARNE T G, LAUFFER J P. The natural excitation technique (NExT) for modal parameter extraction from operating wind turbines: SAND-92-1666[R]. 1993.

[ 2 ] JAMES G H, CARNE T G, LAUFFER J P, et al. Modal testing using natural excitation[C]∥Proceedings of the International Modal Analysis Conference. Albuquerque: SEM Society for Experimental Mechanics Inc., 1992: 1208.

[ 3 ] IBRAHIM S R, MIKULCIK E C. A method for the direct identification of vibration parameters from the free response[J]. The Shock and Vibration Bulletin, 1977, 47(4): 183-

198.

[ 4 ] 杨佑发，李帅，李海龙. 环境激励下结构模态参数识别的改进ITD法[J]. 振动与冲击, 2014, 33(1): 194-199.

YANG Youfa, LI Shuai, LI Hailong. Improved ITD method for structural modal parameter identification under ambient excitation [J]. Journal of Vibration and Shock, 2014, 33(1): 194-199.

[ 5 ] JUANG J N, PAPPA R S. An eigensystem realization algorithm for modal parameter identification and model reduction [J]. Journal of Guidance, Control, and Dynamics, 1985, 8(5): 620-627.

[ 6 ] JUANG J N, PAPPA R S. Effects of noise on modal parameters identified by the eigensystem realization algorithm [J]. Journal of Guidance Control and Dynamics, 1986, 9(3): 294-303.

[ 7 ] 祁泉泉，辛克贵，崔定宇.扩展特征系统实现算法在结构模态参数识别中的应用[J]. 工程力学, 2011, 28(3): 29-

34.

QI Quanquan, XIN Kegui, CUI Dingyu. The application of the extended eigensystem realization algorithm for structural modal parameter identification [J]. Engineering Mechanics, 2011, 28(3): 29-34.

[ 8 ] 崔定宇，辛克贵，祁泉泉.扩展特征系统实现算法的模态参数识别特性研究 [J]. 工程力学, 2013, 30(8): 49-53.

CUI Dingyu,XIN Kegui,QI Quanquan. Modal parameter identification research for the extended eigensystem realization algorithm [J]. Engineering Mechanics, 2013, 30(8): 49-

53.

[ 9 ] 辛峻峰, 张永波, 盛进路. 一种改进的特征系统实现算法及其在海洋平台上的应用[J]. 振动与冲击, 2015, 34(24): 197-201.

XIN Junfeng, ZHANG Yongbo, SHENG Jinlu. Improved eigen system realization algorithm and its application in the analysis of ocean platform [J]. Journal of Vibration and Shock, 2015, 34(24): 197-201.

[10] 包兴先，李昌良，刘志慧. 基于低秩Hankel矩阵逼近的模态参数识别方法 [J]. 振动与冲击,2014, 33(20): 57-

62.

BAO Xingxian, LI Changliang, LIU Zhihui. Modal parameters identification based on low rank approximation of a Hankel matrix [J]. Journal of Vibration and Shock, 2014, 33(20): 57-62.

[11] REYNDERS E. System identification methods for (operational) modal analysis: review and comparison[J]. Archives of Computational Methods in Engineering, 2012, 19(1): 51-124.

[12] 刘宇飞, 辛克贵, 樊健生, 等. 环境激励下结构模态参数识别方法综述[J]. 工程力学, 2014, 31(4):46-53.

LIU Yufei , XIN Kegui, FAN Jiansheng, et al. A review of structure modal identification methods through ambient excitation [J]. Engineering Mechanics, 2014, 31(4): 46-

53.

[13] CAICEDO J M. Practical guidelines for the natural excitation technique (NExT) and the eigensystem realization algorithm (ERA) for modal identification using ambient vibration[J]. Experimental Techniques, 2011, 35(4): 52-58.

[14] HONG Y, LIU X, DONG X, et al. Experimental model updating using frequency response functions[C]∥Proceedings of SPIE: Sensors and Smart Structures Technologies for Civil, Mechanical, and Aerospace Systems. Las Vegas: SPIE， 2016.

[15] ALLEMANG R J, BROWN D L. A correlation coefficient for modal vector analysis[C]∥ Proceedings of the 1st International Modal Analysis Conference. Orlando: IMAC， 1982： 110-116.

[16] ZHU D, GUO J, WANG Y, et al. Field validation of flexure-based mobile sensing nodes on a space frame bridge[C]∥ Proceedings of the 8th International Workshop on Structural Health Monitoring.Stanford: IWSHM, 2011: 13-15.

[17] ZHU D, CHO C, WANG Y. Finite element model updating of a space frame bridge with mobile sensing data[C]∥ Proceedings of SPIE, Sensors and Smart Structures Technologies for Civil, Mechanical, and Aerospace Systems. San Diego: SPIE， 2012.