一种自动识别结构模态参数的随机子空间方法张小宁 段忠东
张小宁 段忠东
摘要： 对结构进行实时在线监测的要求提出了对结构模态参数进行自动识别的需求。目前发展的结构模态参数识别方法均需要人工干预，为实现无人值守的结构实时监测目标，试图发展一种结构模态参数自动识别方法。随机子空间法唯一需要确定的参数是“系统的阶次”，因此，提出了一种基于频率稳定性和振型稳定性自动判别系统阶次的方法，基于此，建立了基于随机子空间法的模态参数自动识别方法；通过两个桥梁算例，对该方法的适用性和鲁棒性进行了验证。关键词： 模态参数识别； 随机子空间法； 斜拉桥； 健康监测； 系统阶次
中图分类号：TU311.3； U441+.3文献标志码： A文章编号： 10044523（2017）04054207
DOI：10.16385/j.cnki.issn.10044523.2017.04.003
引言
随着超高层建筑和大跨桥梁的兴建，健康监测系统越来越多的应用于此类大型土木工程结构中。这些结构的损伤诊断和受损结构评定、实时监测和安全预警、有限元模型修正等成为了健康检测系统的重点，而结构的模态参数识别是需要首先解决的问题[13]。为了能够及时掌握结构的运行状况，及时发现结构损伤，提前做出安全预警，因此需要对结构进行实时在线监测，而这也对结构模态参数自动识别提出了要求。
随机子空间法（SSI）是近年来发展起来的模态参数识别的方法[45]。传统的模态参数识别方法需耗费大量工作来确定具有最小量参数的模型即规范模型，而随机子空间法是基于随机状态空间模型，只需确定“系统的阶次”这一个参数；传统方法都需要进行迭代运算，往往产生发散或者收敛缓慢等问题，而随机子空间法由于算法中的QR分解，矩阵正交投影，奇异值分解等运算不需进行迭代运算，故随机子空间法不存在收敛问题；传统方法往往需要计算高阶次的模型，而随机子空间法通过协方差运算或矩阵正交投影使系统阶次明顯降低，从而运算更加快捷[67]。
1基本原理〖2〗1.1传统SSI定阶方法随机子空间法经典的确定系统阶次的方法主要有奇异值跳跃法和稳定图法。奇异值跳跃法是根据奇异值的跳跃性来确定系统的阶次，奇异值的跳跃点就是真实非零奇异值与零奇异值的分界点，系统的阶次就是奇异值跳跃点之前所有奇异值个数的一半。稳定图法假定系统有多个阶次，计算每一阶次系统的模态参数，对应于某一阶模态，将相邻阶次的模态参数作比较，如果频率、阻尼比和振型的差异小于提前预设的限值，则该点就称为稳定点，稳定点组成稳定轴，从而组成了稳定图。稳定轴对应的阶次为系统的阶次。
这两种方法都需要人工参与，且在人工识别过程中容易受到噪声信号的干扰，很容易确定出错误的系统阶次[8]。
1.2SSI方法系统阶次自动确定原理
系统阶次自动判别方法首先通过求取系统从1到nmax每个阶次的模态参数，用求得的模态参数构造下三角矩阵，然后判定频率的稳定性——模态置信因子（MAF）和振型的稳定性——模态保证准则（MAC），当MAF·∑MAC最大时，该阶次就是系统的最佳阶次，从而实现了自动定阶。
1.2.1求取下三角矩阵
假定系统的阶次依次取1，2，3，4，…，nmax，nmax的取值至少要大于系统的真实阶次。求得每个阶次对应的系统状态空间模型，此时系统矩阵的阶数为2，4，6，8，…，2nmax。
随机子空间法是基于随机系统状态空间模型的方法[10]，随机系统状态空间方程如下所示xk+1=Axk+wk
yk=Cxk+vk （1）式中A为系统的状态矩阵；C为系统输出矩阵；wk为k时刻的过程噪声；vk为k时刻的测量噪声；x（t）为系统的状态向量；y（t）为系统的输出向量。
本文采用的是基于数据驱动的随机子空间方法（SSIDATA），首先将输出数据直接组成Hankel矩阵；然后对Hankel矩阵进行QR分解，得到投影矩阵；对投影矩阵进行奇异值分解（SVD），得到扩展的可观测矩阵和系统状态的卡尔曼滤波；最后由卡尔曼滤波序列和系统输出采用最小二乘法求得系统的状态矩阵A和输出矩阵C。对系统状态矩阵A进行特征值分解A=ΦΛΦ-1（2）式中Λ为包含复特征值的对角矩阵； Φ为特征向量组成的矩阵。
上式中，Λ=diag[λi]，其特征值两两共轭。令λi=ai±jbi （3）可算得系统的振动频率fi和阻尼比εi：fi = a2i + b2i /（2π）（4）
εi = -ai a2i + b2i （5）将所求模态参数构造下三角矩阵，故最终求得了3个nmax阶的下三角方阵Fnmax=f1，100…0
f2，1f2，20…0
f3，1f3，2f3，3…0
……………
fnmax，1fnmax，2fnmax，3…fnmax，nmax（6）式中Fnmax为 振动频率的下三角矩阵。fa，b为 第a阶次的第b阶振动频率。
同时也可以得到模态保证准则（MAC）和阻尼比的下三角矩阵。
1.2.2计算模态置信因子（MAF）
第4期张小宁，等：一种自动识别结构模态参数的随机子空间方法振 动 工 程 学 报第30卷欧氏距离也被称为欧几里得度量，是比较常用的对距离的定义，它表示n维空间中两点之间的距离，对于n维空间，欧氏距离可以表示为d（x1，x2，…，xn）=[（x1-x、1）2+
（x2-x、2）2+…+（xn-x、n）2]12（7）式中d（x1，x2，…，xn）为n维空间的欧氏距离；xn为所求阶次的第n阶模态参数（频率、阻尼比、MAC）。
在这里，定义一种专门针对于振动频率稳定性的因子——模态置信因子（MAF），模态置信因子在数值上等于欧氏距离的倒数（d（x1，x2，…，xn））-1：MAFa，b=[（xa，1-xb，1）2+（xa，2-xb，2）2+
…+（xa，n-xb，n）2]-12（8）式中MAFa，b为a阶次模态（频率）和b阶次模态之间的模态置信因子；xa，b为a阶次的第b阶振动频率或阻尼比；
因为随机子空间法所求出的各阶次的模态参数之间的差异先是随着系统阶次的升高而减小，差异达到最小后又会随着系统阶次的升高而增大，欧氏距离表示了相邻阶次模态参数之间的差异，欧氏距离和系统阶次之间的关系如图1所示。所以，MAFa，b越小（欧氏距离越大）表示该阶模态的置信度越低，MAFa，b越大（欧氏距离越小）表示该阶模态的置信度越高，越接近系统真实模态。图1欧氏距离随系统阶次的变化
Fig.1The variation of Euclidean distance with the order of system
实际应用中，往往只取系统的前几阶模态参数，式（8）中n的取值为所取模态参数的最高阶数。计算出nmax-1个阶次的欧氏距离d1，d2，…，dnmax-1，求出对应的nmax-1个模态置信因子MAF2，1，MAF3，2，…，MAFnmax，nmax-1，然后选取最大的m个模态置信因子MAFi，i-1，MAFj，j-1，MAFk，k-1，…，其所对应的阶次就是系统的近似阶次，建议m取至少大于3。
1.2.3模态保证准则（MAC）
在上一步骤中，选取了模态置信度最高的m个阶次，这m个阶次可以看作是频率稳定的阶次，但它们都是系统的近似阶次，可能含有噪声模态的阶次。为了避免噪声模态，还应该对这m个频率稳定的阶次进行振型的判定，看其振型是否稳定。最后选取这m个阶次中频率稳定性和振型稳定性最高的阶次作为系统的最佳阶次。振型稳定性用模态保证准则（Modal Assurance Criterion）来定义MACa，b=φTa，jφb，j2（φTa，jφa，j）（φTb，jφb，j）（9）式中MACa，b为 两组振型（理论和试验、完好和损伤）的相关性振型；φa，j为 a阶次第j阶模态的振型向量；φb，j為 b阶次第j阶模态的振型向量。
MAC的值处在0和1之间。当MAC=0时表示这两个模态之间完全无关；当MAC=1表示这两个模态之间完全相关。MAF表示某阶次前n阶的频率稳定性，∑MAC表示该阶次前n阶振型模态保证准则的和。∑MAC越高则该阶次的振型稳定性越高。
1.3自动定阶方法具体步骤
系统阶次自动判别的具体步骤如下：
（1） 假定系统的阶次依次取1，2，3，4，…，nmax，nmax的取值至少要大于系统的真实阶次。求得每个阶次对应的系统状态空间模型，此时系统矩阵的阶数为2，4，6，8，…，2nmax，将所求模态参数依次写入下三角矩阵。
（2） 计算相邻阶次（频率矩阵相邻行空间）之间的欧式距离di，从而求得每一阶模态参数的模态置信因子（MAF），选取模态置信因子（MAF）最大的m个MAFi，i-1，MAFj，j-1，MAFk，k-1，…对应的阶次作为系统近似阶次（第一次筛选）。
（3） 计算第二步中的m个系统近似阶次的模态保证准则MACi，MACj， MACk，…。
（4） 最后计算MAF·∑MAC，该值最大时对应的系统阶次就是系统的最佳阶次。求得系统最佳阶次后，频率下三角矩阵和阻尼比下三角矩阵对应行的频率值和阻尼比值以及对应的振型就是系统的真实模态参数（第二次筛选）。
通过MAF及MAF·∑MAC对系统阶次进行两次筛选之后，保证了系统最佳阶次的唯一性。本文将通过工程实验对系统最佳阶次的正确性进行验证。
2算例〖2/3〗2.1三跨连续梁桥模型模态参数识别该三跨连续梁桥模型是在实验室制作而成，图2是实物图。该模型尺寸为3800 mm×200 mm×30 mm，桥总长为3800 mm，桥面宽为200 mm，桥面厚为30 mm。三跨的跨度分别为1000，1700，1000 mm。图2三跨连续梁桥实物图
Fig.2A threespan continuous beam bridge
本次试验数据采集系统使用北京东方振动与噪声技术研究所的整套采集仪器，加速度器的布置如图3所示。试验用加速度传感器可用频率范围为0.5 Hz～8 kHz，采样频率为512 Hz，采样时长为30 s。用小铁锤对模型沿着A，B，C，D，E共5个点依次敲击。图3加速度计布置位置（单位：mm）
Fig.3The position of acceleration sensorJ（Unit：mm）
2.1.1与传统SSI方法识别结果比较
应用本文提出的自动识别方法与传统方法的识别结果比较，这三种方法所求得的前4阶频率如表1所示。
表 1随机子空间法三种方法识别结果
Tab.1The results of three SSI methods
阶次SSI自动识别方法SSI稳定图法SSI奇异值跳跃法126.36426.36326.362257.54157.76957.770364.87964.83264.832490.54490.41390.437
由表1知，SSI模态参数自动识别方法的识别结果与传统的SSI方法识别的结果非常接近，它们之间的偏差不到1%。
2.1.2与其他模态参数识别方法识别结果比较
将SSI模态参数自动识别方法、频域分解法（FDD）、希尔伯特黄变换（HHT）、ANSYS有限元模态分析的结果汇于表2。表2各方法频率识别结果
Tab.2The results of various related methods
阶次ANSYSSSI自动识别FDDHHT频率偏差/%频率偏差/%频率偏差/%125.00026.3635.45226.3805.52026.4115.644254.68457.7795.66058.1306.30158.2416.505364.40464.8570.70365.0000.92064.3210.129492.82390.4812.52390.7502.23392.9810.170
由表2可知，SSI模態参数自动识别方法、FDD和HHT这三种方法识别的振动频率非常接近，第1和2阶频率与ANSYS有限元结果相差5.5%左右，第3和4阶频率与ANSYS有限元结果最大仅相差2.5%。SSI模态参数自动识别方法识别结果与FDD和HHT识别结果仅相差1%以内。
2.2滨州黄河大桥缩尺模型
滨州黄河大桥实验室缩尺模型根据刚度相似准则确定构件尺寸，选用1/40 的缩尺比例。缩尺模型桥面长15.2 m，桥面宽0.82 m，中塔高3.1 m，边塔高1.9 m，桥面板最大板厚20 mm，斜腹板厚6 mm。主梁材料采用铝合金，拉索采用拔丝钢筋。主梁、桥塔和斜拉索的EA和EI刚度相似比误差都控制在5%以内。
用ANSYS对滨州黄河大桥的缩尺模型进行了建模，主梁采用BEAM4单元，弹性模量为6.86×104 N/mm2，泊松比取0.34，密度为3.2×10-3 kg/mm3；主塔、边塔和横隔梁采用BEAM188单元，弹性模量为6.86×107 N/mm2，泊松比取0.34，密度为2.7×10-3 kg/mm3；斜拉索采用LINK8单元，弹性模量为20×104 N/mm2，泊松比取0.3，密度为7.85×10-3 kg/mm3。ANSYS有限元模型如图4所示。图4缩尺模型的ANSYS有限元模型
Fig.4The ANSYS finite element model2.2.1振动试验过程描述
本次试验的加速度传感器可用频率范围为0.5 Hz～8 kHz，采样频率为256 Hz，采样时长为50 s。由于该桥沿主塔严格对称，其振型也是沿主塔的主轴线对称或反对称。所以将加速度计布置在桥的左半侧或右半侧求得半桥的振型即可。本实验将加速度计布置在桥的左半侧，第1个加速度计布置在桥左侧边缘0.25 m处，加速度计之间间隔0.45 m，第16个加速度计布置在中塔左侧0.5 m处，如图5所示。图5加速度计布置图（单位：m）
Fig.5The position of acceleration sensor（Unit：m）
通过人在桥上走动作为桥的外部激励方式。本次试验共采集到215组数据。
2.2.2模态参数识别
将采集到的数据用本文所提出的SSI模态参数自动识别法进行识别，求得的系统阶次为20，如图6所示。图 6系统的最佳阶次
Fig.6The optimal order of the system
用SSI模态参数自动识别法对每一组数据进行识别，表3是一组数据所求系统的前6阶模态与ANSYS有限元结果的对比。
SSI模态参数自动识别法识别了该桥的左半桥的振型，与ANSYS识别的前6阶振型对比如图7所示。图7前6阶振型
Fig.7The first six mode shapes
表3SSI自动识别法与ANSYS有限元频率识别结果
Tab.3The results of SSI automatic identification method and ANSYS finite element for frequency
阶次SSI自动
识别法/HzANSYS
理论值/Hz振型描述频率值
偏差/%14.0714.110一阶反对称0.97528.6299.540一阶对称10.557310.68410.995二阶反对称2.911411.89612.822二阶对称7.784515.26915.494三阶反对称1.452616.32217.286三阶对称5.906
由表3可知，SSI模态参数自动识别方法识别的滨州黄河大桥缩尺模型前6阶频率与ANSYS模态分析结果相比较，相差均在10%以内。由图7可知，SSI模态参数自动识别方法识别的滨州黄河大桥缩尺模型前6阶振型与ANSYS模态分析结果在走势和形状上比较一致。可见本文提出的SSI模态参数自动识别方法可以很好地识别出该模型的模态参数。
3结论
本文提出了基于随机子空间法的模态参数识别方法SSI模态参数自动识别法， 将SSI模态参数自动识别法识别结果与ANSYS模态分析的结果以及其他模态参数的识别结果进行了对比：
（1） 本文提出了一种系统阶次自动判别方法，并对该方法从理论上进行了推导，给出了自动定阶方法计算流程，论证了系统最佳阶次识别结果的唯一性，通过实验验证了识别结果的正确性。
（2）相比较于奇异值跳跃法和稳定图法，自动判别方法克服了它们的缺点，由于不需要人工参与，该方法实现了系统阶次的自动识别，其识别结果是不受人的主观判断的影响。
（3） 提出了SSI模态参数自动识别方法，对3个不同类型结构的模态参数进行识别，该方法识别出的振动频率与ANSYS有限元结果、传统SSI法、FDD和HHT识别结果很接近，表明SSI模态参数自动识别法可以很好地识别出这3种结构的振动频率和振型，且具有较高的鲁棒性。因此SSI模态参数自动识别法可以作为一种模态参数自动识别方法运用于结构实时在线参数识别。
总之，通过这两个算例的计算，本文提出的SSI模态参数自动识别法能很好地识别出结构的频率和振型，表现出了良好地适用性和鲁棒性。
参考文献：
[1]常军， 张启伟， 孙利民. 随机子空间产生虚假模态及模态遗漏的原因分析[J].工程力学， 2007， 24（11）： 57—62.
CHANG Jun， ZHANG Qiwei， SUN Limin.Analysis how stochastic subspace identification brings false modes and mode absence[J]. Engineering Mechanics， 2007， 24（11）： 57—62.
[2]Lennart L . System Identification： Theory for the User[M]. 北京： 清华大学出版社， 2002： 219—221.
[3]Herman V A， Bart P. International research projects on structural health monitoring： an overview[J]. Structural Health Monitoring， 2003， 2（4）： 341—358.
[4]Bart P. System identification and damage detection in civil engineering[D]. Belgium： Department of Civil Engineering， K.U. Leuven， 2000： 123—127.
[5]Bart P， Roeck D. Referencebased stochastic subspace identification for outputonly modal analysis [J]. Mechanical Systems and Signal Processing， 1999， 13（6）： 855—878.
[6]單德山， 徐敏. 数据驱动随机子空间算法的桥梁运营模态分析[J]. 桥梁建设， 2011， 6： 16—21.
SHAN Deshan， XU Min. Operational modal analysis of bridge structure based on datadriven stochastic subspace identification algorithm[J].Bridge Construction， 2011， 6： 16—21.
[7]Bart P， Roeck D. Stochastic system identification for operational modal analysis： a review[J]. Journal of Dynamic Systems， Measurement and Control， 2001， 123（4）： 659—667.
[8]Filipe M， Cunha A， Caetano E. Online automatic identification of the modal parameters of a long span arch bridge[J]. Mechanical Systems and Signal Processing， 2009， 23 （5）： 316—329.
[9]Andersen P， Brincker R， Goursat M， et al. Automated modal parameter estimation for operational modal analysis of large systems[C]. Proceedings of the 2nd International Operational Modal Analysis Conference， Copenhagen： 2007： 299—308.
[10]肖祥， 任伟新. 实时工作模态参数数据驱动随机子空间识别[J]. 振动与冲击， 2009， 28（8）： 148—153.
XIAO Xiang，REN Weixin.Improved datadriven stochastic subspace identification of online operational modal parameters[J]. Journal of Vibration and Shock， 2009， 28（8）： 148—153.
A stochastic subspace identification method for automotive identification
of structural modal parameters
ZHANG Xiaoning1， DUAN Zhongdong2
（1.The Ningxia Hui Autonomous Region Electric Power Design Institute，Yinchuan 750000， China；
2. Harbin Institute of Technology（Shenzhen）， Shenzhen 518055， China）
Abstract： Realtime online monitoring of the structure put forward the demand for the automatic identification of the structural modal parameters. The current methods of the structural modal parameter identification require manual intervention， in order to realize realtime monitoring of the structure， this article attempts to develop an automatic structural modal parameters identification method. Stochastic subspace method only needs to determine one parameter， which is "system order". This paper proposes a method for automatic identification system order based on the frequency stability and vibration mode stability， based on those， we established an automatic identification method for structural modal parameters. Through identification two bridge examples， the reliability and applicability of the method was verified.Key words： modal parameters identification； stochastic subspace method； cablestayed bridge； health monitoring； system order作者简介：张小宁（1988—），男，硕士研究生。电话：18295390955。Email： zxn880625@163.com
