0 引言

模态分析技术可分为3大类^[1]：其一是基于计算机仿真的有限元分析法(Finite Element Analysis,FEA)；其二为基于输入(激励)/输出(响应)模态试验模态分析法(Experimental Modal Analysis,EMA)；其三是基于仅有输出(响应)模态试验的运行模态分析法(Operational Modal Analysis,OMA)。第一类属于结构动力学正问题，第二类与第三类属于结构动力学反问题。基于仅有输出(响应)的模态参数识别技术(OMA)，因其无需增加特殊的激励设备，基本不影响结构的正常运行，对结构损伤小，测试费用低等优点，被广泛应用于土木工程领域^[2]。

模态参数(固有频率、阻尼比、振型)识别是结构动力学特性分析、结构健康监测、损伤识别等方面的关键内容，特别是对于服役期桥梁健康监测，能否利用在线监测数据实时识别模态参数，对识别方法及算法提出来更高的挑战。随机子空间识别法(Stochastic Subspace Identification，SSI)以其输入参数少，可程序化程度高等优点受到国内外学者的热捧。SSI假定输入信号为随机白噪声，利用自然环境激励条件下进行结构模态参数识别，基于离散时间状态空间方程，直接处理时间序列的时域方法。比利时鲁汶大学的Bart Peeter等人提出基于数据的子空间识别的概念及思想，并把系统控制理论与矩阵理论联系起来，利用稳定图确定系统的阶次^[3-4]；福州大学任伟新教授团队提出基于稳定图的平均正则化稳定图算法，使得模态参数识别效率和识别精度得到提高，提炼出经验模式分解(EMD)的SSI，并应用于桥梁和建筑物等工程实例^[5-6]；同济大学的孙利民、张启伟、常军等人提出两阶段稳定图方法剔除虚假模态，消除模态遗漏现象，对Hankel矩阵中的“过去”输出数据缩减，以提高计算效率^[7-8]；重庆大学的章国稳、汤宝平课题组利用谱系聚类法实现模态参数的自动识别，依据模态能量水平剔除虚假模态，提出基于数据缩减的分频段小波变换模态参数快速自动识别^[9]；刘兴汉、王跃宇等人提出基于Cholesky分解的改进随机子空间识别法，通过对维数大大减少的Cholesky分解来得到投影矩阵，提高计算效率^[10]；此外，香港理工大学、东南大学、清华大学、西北工业大学等多所科研院校对SSI进行了深入的研究和探讨。

经过40余年的发展，国内外研究学者在SSI方法上取得很大成就，但是仍面临许多问题，特别是SSI方法的计算效率有待提高。基于数据驱动SSI方法，需要构造单无限大的Hankel矩阵，需要处理的数据非常庞大，并且涉及到高维QR分解和SVD分解，需要的时间和内存开销很大程度上影响了其在实际工程中的应用。针对以上问题，本文对基于数据驱动SSI方法进行优化，首先，以部分测点信息作为“过去”输出信号代替全部测点信号，精简Hankel矩阵中的“过去”输出数据Y_p，减少计算量，针对文献^[7]可能出现模态遗漏的问题，本研究通过基于模态置信度(Modal Assurance Criterion，MAC)的数据优化，有效避免模态遗漏。其次，通过公式推导，避开求解投影矩阵简化计算步骤，避免高维矩阵的分解和存储，改善计算机使用内存，提高计算效率。最后，以西宁北川河桥为工程实例，验证该优化算法的实用性和有效性。

1 结构振动的状态空间模型的计算法

结构振动微分方程为:

(1)

式中，M∈R^n×n，C∈R^n×n，K∈R^n×n分别为质量矩阵、阻尼矩阵、刚度矩阵;u(t)∈Rⁿ为位移向量;F(t)∈R^n×l外部激励向量。

令

(2)

则式(1)可转化为：

(3)

引入观测方程：

(4)

式中，y(t)∈R^l×1为观测向量(l表示测点数)；A_c∈R^2n×2n为连续状态矩阵；B_c∈R^2n×2n为连续时间输入矩阵；C_c∈R^l×2n为连续时间输出矩阵；D_c∈R^l×m为直馈矩阵。

自然环境激励通常为风、车辆或行人等，假定是理想的白噪声激励，并与输入噪声和输出噪声合并，将式(3)与式(4)离散化可得离散时间状态空间模型：

(5)

(6)

式中，y_k∈R^l×1为第k(k∈N)个采样间隔(Δt)的输出向量；x_k∈R^2n×1为状态向量；A∈R^2n×2n为离散状态矩阵；C∈R^l×2n为离散输出矩阵；E为数学期望算子，δ_pq为Kronecker Delta函数；w_k∈R^l×1，v_k∈R^2n×1为测量噪声、过程噪声；Q∈R^2n×2n，S∈R^2n×l和R∈R^l×l为噪声序列w_k，v_k的协方差矩阵。

2 随机子空间优化算法

SSI可分为：基于协方差驱动随机子空间识别法(Cov-SSI)与基于数据驱动随机子空间识别法(Data-SSI)，其基本原理都是以线性系统的最小实现理论为基础，利用Hankel 矩阵与系统可观测性、可控制性矩阵的特殊关系求得系统矩阵A和输出矩阵C，在状态空间辨识系统的模态参数。SSI利用矩阵分析QR分解使得数据缩减，采用奇异值分解(SVD)剔除噪声干扰，运用最小二乘法来识别系统的状态矩阵，最后通过特征值分解直接确定结构模态参数：自振频率、阻尼比和振型。

2.1 构造Hankel矩阵

构造Hankel矩阵时，数据量很大(单无限)，虽然现在计算机硬件得到很好的改善，但Hankel矩阵的QR分解会消耗计算机很大内存，计算耗时长。

2.1.1 基于MAC准则的数据优化

根据结构响应的特点，各个测点的输出信号都包含相同的频率信号-固有频率，由投影原理可知，所有的输出信号都可以单独或者部分组合作为被投影信号，即精简“过去”输出数据Y_p，但是牺牲参与信息以减少计算量伴随有模态遗漏等问题。为避免各向量之间的交角偏小让系统遗漏关键模态的情况发生，MAC算法可以在选择测点数据时，使量测的模态向量保持较大的空间交角，增大模态向量的识别效果，尽可能地保留原模型的特性，模态置信度MAC矩阵(简称M阵)表达式如下：

(7)

式中,Φ_i和Φ_j分别为第i阶和第j阶模态向量。通过检查各模态在量测自由度上形成的向量M阵的非对角元，即可判断出相应两模态向量的交角情况。当M阵的某一元素M_ij (i≠j)等于1时表明第i向量与第j向量交角为零，两向量不可分辨；而当M_ij(i≠j)等于零时，则表明第i向量与第j向量相互正交，两向量可以轻易识别。故选择数据的测点布置应力求使M阵非对角元向最小化发展，一般建议非对角元最大取值为0.25^[11]。

当从全部测点数据中删除k行，式(17)变为：

(8)

从已有测点数据里删除某点数据是一个反复迭代，不断寻优的过程，其根本目的是删除使M阵非对角元最大的测点数据。

2.1.2 构造数据精简的Hankel矩阵

从l个测量数据中选择r个作为“过去”输出数据，即：

(9)

(10)

式中，L为选择矩阵;I_r为r阶单位对角矩阵；y_i^r为i时刻“过去”输出信号的块元素；y_i为使得M阵非对角元减小最快测点至最慢测点重新排列后的i时刻所有输出信号的块元素。

构造Hankel矩阵：

(11)

式中，“过去”输出数据Y_p∈R^ri×j；“将来”输出数据Y_f∈R^li×j。

2.2 Hankel矩阵QR分解

SSI通过QR分解对Hankel矩阵进行数据缩减：

(12)

式中，Q为正交矩阵：Q^TQ=QQ^T=I；R为下三角矩阵。

Y_f向Y_p的投影矩阵^[7]：

(13)

定义可观测矩阵：

(14)

则投影矩阵可表示为：

(15)

2.3 求解观测矩阵Γ_i

Data-SSI传统的计算方法，投影矩阵O_i贯穿整个计算过程，但是投影矩阵的列数与Hankel矩阵列数相同，依然是单无限，并且求解过程很复杂，如果能避开求解投影矩阵，不仅减少计算步骤，还能改善计算机内存，从而提高计算效率。

可观测矩阵Γ_i与投影矩阵O_i相关联，投影矩阵O_i又可用矩阵R₂₁表示，如果能找到可观测矩阵Γ_i与矩阵R₂₁的关系式，就可以避免求解投影矩阵，推导如下：

对矩阵R₂₁进行奇异值(SVD)分解：

(16)

式中，U，V属于酉矩阵；S=diag(σ₁,σ₂,…,σ_n)，σ_i为从大到小的全部非零奇异值。

则式(13)可表示为：

(17)

根据矩阵运算的性质^[12]，Q₁^T属于正交矩阵，必属于酉矩阵；两个酉矩阵的乘积Q₁V₁必属于酉矩阵，由此推知式(17)为投影矩阵O_i奇异值分解的一种形式。

则可观测矩阵可表示为^[7]：

(18)

至此，可观测矩阵Γ_i与矩阵R₂₁的SVD分解建立关系，优化算法避免了求解投影矩阵O_i，Hankel矩阵QR分解时只存储下三角矩阵R，利用子矩阵R₂₁的SVD分解求得可观测矩阵Γ_i。计算过程中避开了大矩阵的存储和分解(R₂₁要比O_i小很多)，并且不用存储正交矩阵Q，改善了计算机的使用内存，提高了计算效率。

2.4 求解系统矩阵A和C

定义观测矩阵Γ_i的两个子阵Γ₁和Γ₂：

(19)

系统矩阵A：

(20)

式中，(•)⁺表示矩阵的伪逆。

通过观测矩阵Γ_i求得系统矩阵A与C，实际工程应用中，在确定系统阶次时结合稳定图的方法计算不同的A，C从而识别系统参数。

3 工程算例

本文采用与文献^{[7, 13-14]}相同的工程实例，中承式钢管混凝土系杆拱桥-西宁北川河桥为工程实例，便于比较计算方法的有效性。

3.1 工程概况

西宁北川河桥，净跨90 m，矢跨比1/5，拱轴线为悬链线，拱轴系数m=1.167；桥面净宽21.6m，共16对吊杆；主拱圈采用4根φ650 mm×10 mm 钢管，内灌注C50混凝土，横梁采用C40钢筋混凝土，吊杆采用127×φ5高强碳素钢丝，水平系杆选用φ5高强钢丝束。

环境激励是利用自然风载、随机交通荷载或其他环境激励及其组合等，测量桥梁的动力特性，具有简单易行、不中断交通、对桥梁结构的损伤小、测试费用低等优点，应用越来越广泛。全桥上下游每侧布置竖向测点16个，共32个，位于每根吊杆下侧，一个固定参考点，位于桥面起始处，传感器横向布置如图 2所示。全桥共分4组进行测试，每组8个测点，信号采样频率80 Hz。以20 Hz重采样、组成36行18176列的原始数据。

3.2 数据优选

首先确定竖向弯曲振动特征作为有效振型，提取对桥梁的竖向振动有贡献的前4阶振型数据来考虑作为数据优选的原始数据，并将这4阶振型组成模态向量矩阵(模态向量选用文献^[14]方法识别结果)。然后以全部测点为初始状态进行筛选，利用MAC法，依次删除使得MAC阵非对角元减少最小的测点数据，剩余测点数据为候选，反复迭代，逐步删除更多测点数据，直到满足预设值要求(即MAC阵非对角元小于0.05)，并依次记录删除测点信息。最后，按照测点删除顺序反向排列测点数据，即最先删除测点排列到最后，保留测点排列前几行。最终优化测点布置方案如图 3所示。

3.3 检验与对比分析

基于Matlab平台，使用C语言编写程序，识别桥梁结构的频率，以相同的标准：最小阶次取10，最大阶次取50；稳定点个数取15；判定标准：频率0.01、阻尼0.4、振型0.05。图 4中连续曲线为平均正则化功率谱，峰值法提取的频率与SSI识别频率形成直观对比。

当r=4，i=50，无数据优选，即取前4个测点数据作为“过去”输出数据：Y_p∈R^{200×18 176},Y_f∈R^{800×18 176}，形成的稳定图见图 4(a)。

当r=8，i=50时：Y_p∈R^{400×18 176},Y_f∈R^{800×18 176}，本文优化算法形成的稳定图见图 4(b)。

文献^[14]选用桥梁一侧全部的16个测点数据，无数据优选，i=50时：Y_p∈R^{800×18 176},Y_f∈R^{800×18 176}，形成的稳定图见图 4(c)。

观察图 4(a)第4阶频率(4.589)完全稳定次数只有9次，无法提取，形成遗漏。对比图 4(a)与(b)，前4阶频率完全稳定次数都有所增加，特别是第4阶频率全部完全稳定，将成为优势频率。观察图 4(c)可知，文献^[14]算法的第4阶频率(4.589)完全稳定次数12次，第一阶频率(2.014)完全稳定次数14次，如果要提取前4阶频率，稳定点个数标准应降低到12次。观察图 4(b)，本文优化算法前4阶频率完全稳定次数依次是：17，22，22，22，3个稳定图中都是最多。通过以上分析，r取值受到振型稳定的影响，不宜太低并且与选择的位置相关，本文优化算法通过数据优选有效地解决了此问题，避免了模态遗漏；在识别高阶次模态时，为避免形成更多的虚假模态，完全稳定次数的提高为判定标准提供更广的选择空间。

各文献及本文优化算法识别的频率如表 1所示，对比各组数据，本文优化算法结构与其他方法差别不大，误差范围在5%以内。优化算法识别结果比较理想。优化算法的前4阶竖向振型如图 5所示。

耗时方面，以下所有耗时均采用同一计算机：windous XP系统、3G内存。

如果Hankel矩阵j取5 000列，i取50块行，文献^[14]算法需要14.641s，本文优化算法仅需3.045 s，识别时间仅为前者的20.8%。如果Hankel矩阵的列取5 000固定不变，二者耗时对比如图 6所示，随着块行数的增加，本文优化算法优势越来越明显，当块行数i大于120时，文献^[14]的算法由于内存超界(out of memory)而无法计算。如果Hankel矩阵i取50不变，文献^[14]算法的识别时间约是本文算法的5倍，如图 7所示，并且当列数j大于12 500时，前者算法因内存超界而无法计算。观察图 6、图 7，块行i增加对识别时间的影响远远大于列j，识别时间随i呈指数增长，随j呈线性趋势增长；Hankel矩阵的行数为2li或者(r+l)i行，Hankel矩阵的大小对i的变化更敏感所致。

4 结论

本文推导了SSI优化算法，并通过工程算例分析及对比验证了其实用性和有效性。在稳定图中包含了平均正则化功率谱，相互校核，确保识别结果的可靠性，本文优化算法具有如下优点：

(1) 通过MAC准则优化选取数据，缩减了Hankel矩阵的“过去”输出数据Y_p，有效避免了由于数据减少带来模态遗漏；

(2) 避开求解投影矩阵O_i，利用Hankel矩阵R₂₁求得可观测矩阵Γ_i；

(3) 计算过程中避开了高维矩阵的存储、QR分解和SVD分解，减少了计算量，并且不用存储正交矩阵Q，很大程度上改善了计算机的使用内存。

通过图 6和图 7可以明显看出，本文优化算法耗时更短，在相同计算机配置下，该算法处理的数据能力更强大，特别是对于梁拱组合桥、斜拉桥、悬索桥等需要处理海量数据的柔性结构，该算法更具优势。同时，由于识别速度快、方法稳定、精度高等优点，成为桥梁的健康监测、桥梁在线损伤识别和评估等领域更有力的工具。