桥梁时变系统的张量子空间识别方法研究及试验验证

2021-12-20张二华单德山

振动与冲击 2021年23期

张二华，吴涤，刘昊，单德山

(1.四川省公路规划勘察设计研究院有限公司,成都 610031；2.西南交通大学土木工程学院，成都 610031)

近年来，随着结构健康诊断技术的蓬勃发展并广泛应用于大跨度桥梁，基于健康监测实测信号的桥梁系统识别研究方兴未艾[1]。基于实测信号的系统识别方法易于从桥梁各种激励环境下高精度获得桥梁的模态信息，已成为桥梁状态评估的关键[2]。其中子空间系统识别方法因人为干预少，识别结果鲁棒性好，在桥梁结构系统识别中日益受到重视。目前应用于桥梁系统识别的子空间方法主要有：随机子空间方法、递归子空间方法、确定子空间方法。国内外学者对3类方法均开展了广泛深入的研究。

Zarbaf等[3]用随机子空间方法与层次聚类算法，对某斜拉桥斜拉索力进行了估计；徐健[4]将随机子空间方法与经验模态分解方法相结合，对某大跨斜拉桥进行了系统识别；Bui等[5]分别采用能量谱密度、协方差驱动随机子空间及数据驱动随机子空间的方法，并结合峰值拾取法和稳定图，对某混凝土简支梁进行了系统研究。Khan等[6]采用递归随机子空间方法，研究了某大跨斜拉桥温度对模态参数影响规律；Huang等[7]讨论了递归随机子空间方法中遗忘因子、系统矩阵维度大小等参数对地震激励下结构模态参数识别的影响。秦世强等[8]将确定性随机子空间方法用于菜园坝长江大桥的模态参数识别；陈栋军[9]采用改进后的确定-随机子空间方法识别了某混凝土斜拉桥的模态参数。Huang等[10]采用确定子空间方法识别了地震激励下多输出系统的模态参数。

综上所述，随机子空间方法常用结构一次完整动力响应来识别模态参数，而忽略非线性结构模态参数的时变特性。递归随机子空间算法，能反映结构模态参数的时变规律，但识别结果精度较低，难以用于大跨复杂桥梁的时变模态参数识别。确定性随机子空间识别精度较高，但对实际桥梁激励和响应数据的完备性要求较高[11]。并且，子空间系统识别方法的现有改进多集中于算法的运算效率、数据前处理技术、模态自动识别方面，尚未看到有从子空间方法的基础理论——矩阵计算方面加以改进和突破的研究文献。

针对矩阵运算的子空间系统识别理论与方法的缺陷，本文基于短时滑窗子空间系统识别框架，在原有2维Hankel矩阵基础上，引入时间尺度，建立3维时变Hankel张量；结合张量平行因子分解理论，提出基于张量的平行系统矩阵估计方法；结合稳定图方法从模态参数识别结果中筛选真实模态信息，提出一种新的张量子空间系统识别理论与方法。基于曲线斜拉模型桥、斜拉桥振动台试验的实测数据，以计算耗时、稳定图稳定轴及虚假模态为评判指标，通过与已有短时滑窗子空间系统识别方法对比，验证了本文所提方法较现有方法的优越性与工程实用性。

1 张量子空间系统识别理论

1.1 张量子空间系统识别数学模型

1.1.1 滑窗子空间系统识别框架

根据空间状态控制理论[12]，连续确定性系统的状态空间模型可写为如式(1)和式(2)所示

(1)

y(t)=Ccx(t)+Dcu(t)

(2)

实际运用中，系统采样数据均为离散形式，将连续状态空间模型转变为离散状态空间模型，式(1)～式(2)改写为

xk+1=Axk+Buk+wk

(3)

yk=Cxk+Duk+vk

(4)

式中:xk=x(kΔt);A,B,C,D为连续系统矩阵Ac,Bc,Cc,Dc的离散形式;uk,yk为离散输入和输出矩阵;wk为过程噪声;vk为测量噪声。

当系统处于环境激励下，其激励源往往不可测，通常将环境激励假定为与随机噪声影响相似，则含有uk项与噪声项可合并,随机子空间系统状态空间模型可写为

xk+1=Axk+wk

(5)

yk=Cxk+vk

(6)

无论是确定子空间识别或随机子空间系统识别问题，均可通过将输入-输出数据或输出数据组成特定的Hankel矩阵，计算该矩阵的行空间投影，对投影进行奇异值分解(singular value decomposition,SVD)，从而得到可观测矩阵和状态序列的卡尔曼滤波估计，进而确定系统矩阵以及噪声协方差矩阵。然后经特征值分解，获得系统的模态参数，其求解过程可参考Van Overschee等的研究。

因确定子空间与随机子空间系统识别方法均基于线性系统假设，近年来，Moaven等[13]将滑窗技术引入到上述方法中，以解决线性子空间系统识别方法在追踪复杂非线性时变系统模态参数方面的不足。其核心思想为：基于短时平稳性假设，对测试数据进行分段，各窗口内将结构视为时不变系统；然后对每个窗口内动力数据重复使用子空间识别方法求解结构模态参数，达到对非线性系统时变模态参数的追踪。算法流程图,如图1所示。

图1 滑窗子空间系统识别算法流程图Fig.1 The sliding window based subspace system identification flow chart

由图1可知：上述方法未从根本上对基于矩阵运算的子空间系统识别理论和方法进行改进，而是不同时段内对方法的重复使用，因此其耗时为N次子空间求解过程花费时间的总和。

1.1.2 张量子空间数学模型

为保持和张量表达形式一致，滑窗子空间中t时刻的输入输出Hankel矩阵重新写为U2i,j,t、Y2i,j,t，其中第一下标2i表示t时刻Hankel矩阵的行数，第二下标j表示Hankel矩阵的列数，第三下标t表示当前的时刻。假设t+1时刻新的数据列为ut+1和yt+1，形成如式(7)的行向量空间

(7)

将式(7)行向量空间附于t时刻Hankel矩阵的最右一列，同时将原U2i,j,t、Y2i,j,t最左一列剔除，保持t+1时刻的Hankel矩阵与t时刻Hankel矩阵维度相同，形成t+1时刻的输入输出Hankel矩阵，如式(8)所示

(8)

假设信号时间窗口为L，可将信号时程为T的时变振动信号划分为N个时间窗，则t=1,2,…,N，按信号采集时间先后，构建时变3维Hankel张量，3维张量沿时间维度的切片为

(9)

基于上述时变3维Hankel张量，对系统矩阵进行估计，首先需要解决的问题是剖析其数学本质，即建立对应的数学模型进行求解。为此，以下将首先通过理论推导分析非线性系统的系统矩阵的时变性，在此基础上，建立系统矩阵估计的数学模型。

(10)

式中:矩阵U∈RM×M,V∈RN×N为正交矩阵;UsU(1:M,1:r);VsV(1:N,1:r);Σs=diag(σ1,…,σr),σ1≥σ2≥…≥σr≥σr+1=…σM=0。

由式(10)经推导，可得到

(11)

在X中加入白噪声项Ψ，则可得到受噪声污染的实测信号为

Y=X+Ψ

(12)

式中，Ψ为测量噪声分量。

当数据长度足够大时，存在如式(13)关系

(13)

进一步可得到式(14)

(14)

(15)

式中,S=diag(s1,…,sM),si>0。si表达式为

(16)

由式(15)和式(16)可知，含噪数据组成矩阵Y经SVD后左奇异矩阵U和对角阵S包含结构时变动力指纹信息，对于时变的非线性系统，左奇异矩阵U和对角阵S集中反映了系统时变的特性

由张量展开和高阶张量SVD，时变Hankel张量沿时间维度展开的切片可类似看作某时段如式(12)含噪信号Y，则由式(9)～式(15)可得，张量子空间系统矩阵求解的数学模型为

H(k)=UkSkVT+Rk

(17)

式中:Rk为拟合的残差平方和;Uk为酋矩阵;Sk为非负对角阵;k为沿时间轴展开的3维Hankel张量H“切片”下角标。

式(17)在形式上与张量平行因子分解理论中的PARAFAC模型改进形式PARAFAC2模型保持一致。因此，以下将引入改进的张量平行因子分解理论，求解上述数学问题。

1.2 基于张量分解的系统矩阵估计

1.2.1 PARAFAC2模型简介

Harshman[14]在1972年提出了PARAFAC2模型，其数学表达为：定义张量的展开式Xk,k=1,…,K为矩阵集合，其中，Xk是Ik×J的矩阵，Ik的值随着k的改变而改变，假设PARAFAC2分解的维度为R，此时，PARAFAC2具有

Xk≈UkSkVT,k=1,…,K

(18)

式中:Uk为Ik×R的矩阵;Sk为R×R的对角矩阵;k=1,…,K;V为J×R的因子矩阵。

PARAFAC2分解示意图,如图2所示。

图2 PARAFAC2分解示意图Fig.2 Illustration of PARAFAC2

Xk≈QkHSkVT,k=1,…,K

(19)

(20)

1.2.2 系统矩阵估计

基于约束的PARAFAC2模型，式(17)可视为广义的主成分分析问题，其问题可表达为

σk(Q1,…,Qk,H,V,S1,…,Sk)=

(21)

首先求解Qk，式(21)的最小化问题，等价于式(22)最大化问题

(22)

(23)

则列正交的Qk为

(24)

通过H,S1,…,Sk,V，最小化式(21)，转变为下述问题

σk(H,V,S1,…,Sk|Q1,…,Qk)=

(25)

综上所述，利用PARAFAC2模型平行求解张量系统矩阵过程可归纳为：

步骤4由式(25)判断求解过程是否达到收敛准则，若σold,k-σnew,k>εσold,k，ε为收敛准则(经试算，本文中ε取值为0.1)，则回到步骤3，否则，进入下一步。

至此，基于PARAFAC2模型，包含结构时变系统参数的系统矩阵Uk和Sk已求解得到。

则每个时间窗内的扩展观测矩阵可按式(26)进行求解

(26)

式中:Wk为权重影响因子;Γk为扩展可观测矩阵。由扩展观测矩阵的移位不变性，可计算获得系统矩阵A,C，经特征值分解运算，即可获得系统的频率、阻尼比和振型，具体计算过程可参考Van Overschee等的研究，由于篇幅所限，不在赘述。

1.3 求解过程对比分析

本文所提出的张量子空间系统识别方法与两种滑窗子空间系统(确定系统或随机系统)识别方法的求解过程流程图,如图3～图4所示。

图3 滑窗子空间系统识别方法流程示意图Fig.3 Schematic diagram of the sliding window-based subspace system identification method

图4 张量子空间系统识别方法求解过程示意图Fig.4 Schematic diagram of the tensor subspace system identification method

由图3～图4可知：假设测试信号划分时间窗口数为K，稳定图中最大系统阶次为N，对于滑窗子空间系统识别方法，每一窗口内数据均需进行一次完整的基于矩阵运算的子空间系统识别过程，则需进行的SVD次数总共为K·N/2；而本文中张量子空间系统识别方法，通过K个时间窗口构建时变3维Hankel张量，在系统识别过程中，仅在稳定图分析过程中，每一阶次进行两次SVD，系统识别全过程SVD次数总共为2·N/2=N次，随着K增大(K为大于2的整数)，SVD次数将成倍减少。而已有研究表明[16]，子空间系统识别过程中SVD是耗时严重的步骤之一，且SVD过程极易受到噪声影响导致奇异值发生变化，致使最后模态参数识别结果不准确，在稳定图上表现为稳定点散乱，稳定轴扭曲。

综上所述，本文所提出的张量子空间系统识别方法较滑窗子空间系统识别方法过程，计算效率更高，计算结果受噪声影响更少，稳定轴将更加清晰可辨，更易识别到系统更多阶次的模态参数信息。以下通过两座斜拉模型桥进一步验证上述理论分析的结论。

2 试验验证

2.1 曲线斜拉模型桥试验

采用一座几何缩尺比例为1∶20的曲线斜拉模型桥，验证本文所提方法用于桥梁随机子空间系统识别问题的可行性。

模型斜拉桥由5跨构成，跨度布置为(2.45+4.05+14.24+4.05+2.45)m，桥梁总长27.25 m，以中跨跨中为界分为曲线梁段和直线梁段，曲线梁段曲率半径为27.5 m。主梁为Π型钢梁，宽1.1 m，高0.081 m，由板厚6 mm的Q235钢板焊接制作。斜拉索采用双幅索面布置，每幅60对，由10 mm钢丝绳制作。辅助墩与过渡墩由C30混凝土与Φ6 mm光面钢筋制作完成。桥塔为钻石型桥塔，由板厚10 mm的Q235钢板焊接为矩形空心截面，桥塔每侧钢板中轴线位置由6 mm的Q235钢板焊接成加紧肋。桥梁模型示意图,如图5所示。

图5 曲线斜拉模型桥示意图(m)Fig.5 A model cable-stayed bridge (m)

为获得主梁损伤后的动力测试信号，在主梁截面中线布置941B型加速度传感器，传感器布置方式为：分别在两侧边跨跨中各布置1个测点，中跨八分点各布置1个测点，共计11个测点。信号采样频率为256 Hz，在中跨跨中位置处施加脉冲激励后，采集主梁自由振动信号，采集时长为40 s，如图6所示。

稳定图是子空间方法中用于筛选结构真实模态信息的常用方法，本文选取10 s长度信号作为初始数据，滑窗步长为0.4 s，由滑窗随机子空间识别方法与张量子空间系统识别方法分别获得系统稳定图,如图7所示，前3阶识别结果的对比,如图8所示。图7中，“◇”稳定极轴代表结构频率、阻尼比与振型结果均稳定，“*”极轴代表仅频率稳定，“○”极轴代表仅频率与阻尼比稳定。稳定图方法以频率、阻尼比与振型结果均稳定的极轴代表结构的真实模态结果。稳定图稳定极轴条件及判别方法可参考Khan等的方法。本节中曲线斜拉模型桥的理论模态参数识别结果可参考文献[17]。

(a)频率

由图7～图8可知：① 由张量子空间方法所得稳定图红色稳定轴的数量较滑窗随机子空间方法更多，如2.12 Hz处，滑窗随机子空间方法并未得到“◇”稳定轴；张量子空间方法所得虚假模态散点数量更少，如在3～6 Hz内，滑窗子空间方法所得结果出现多个“*”零星散点，张量子空间方法结果零星散点较少，鲁棒性更好；② 张量子空间方法识别到的频率与振型结果与传统子空间方法结果均极为接近，两种方法识别出的阻尼比值差异较大，这主要是由于基于子空间的时域类系统识别方法在识别阻尼比过程中对噪声极为敏感造成的。③ 滑窗随机子空间方法得到上述结果用时337.49 s，而张量子空间识别方法仅用时280.31 s(计算机配置CPU均为I7-6700K，内存均为16 G)。以上结果表明，针对随机系统识别问题，本文所建立的张量子空间方法较滑窗随机子空间方法在大幅提高计算效率的同时，能够从信号中挖掘更多系统模态信息，且识别精度与滑窗随机子空间方法相同。

2.2 斜拉模型桥振动台试验

采用一座缩尺比例为1∶20的振动台斜拉模型桥，验证本文所提方法用于确定子空间系统识别问题的可行性。

斜拉模型桥由5跨构成，跨度布置为(2.9+3.6+19.0+3.6+2.9)m，桥梁总长32 m。斜拉模型桥索塔、过渡墩及辅助墩均选用M15微粒混凝土模拟原桥混凝土材料，钢筋采用φ6圆钢作为纵筋、10#钢丝作为箍筋。模型桥主梁设计为矩形空心断面，采用厚度为5 mm钢板制作。桥塔为门式框架结构，高度为4.55 m；边墩、辅助墩高度为1.9 m。拉索采用直径10 mm钢丝绳模拟。斜拉模型桥示意图,如图9所示。

图9 振动台斜拉模型桥示意图(cm)Fig.9 A model cable-stayed bridge tested on a shaking table (cm)

主梁与桥塔分别布置941B型加速度传感器，测点布置方式为：沿桥塔两侧塔柱间隔1 m布置纵桥向和横桥向加速度计，每侧塔柱布置5个测点；中跨主梁在1/8测点各布置1个竖向和横向加速度计，两边跨跨中及各支点处均布置1个竖向和横向加速度计；各振动台分别布置1个纵横向加速度计；辅助墩墩顶各布置1个纵向和横向加速度计。信号采样频率256 Hz，采样时长为47 s，桥塔地震加速度响应,如图10所示。

本节同样采用稳定图作为结构真实模态的筛选方法，滑窗步长为5 s，由滑窗确定子空间识别方法与张量子空间系统识别方法分别获得系统稳定图，如图11所示，前3阶识别结果的对比如图12所示。图11中，“◇”极轴、“*”极轴、“○”极轴含义与2.1节相同，稳定图方法以频率、阻尼比与振型结果均稳定极轴代表结构的真实模态结果。稳定图稳定极轴条件及判别方法可参考Huang等的研究。本节中斜拉模型桥的理论模态参数识别结果可参考文献[18]。

(a)频率

由图11～图12可知：① 针对确定性系统而言，张量子空间方法所得稳定图中“◇”稳定极轴上的极点数量较滑窗确定子空间方法更多，如0～10 Hz内，滑窗确定子空间方法虽然也有“◇”稳定轴，但稳定极点数量较张量子空间方法少，这表明张量子空间方法求解过程鲁棒性更好，结果更加稳定；② 张量子空间方法识别到的频率与振型结果与传统确定子空间方法结果均极为接近，两种方法识别出的阻尼比值差异较大，其原因与2.1节相同，不再赘述；③ 滑窗确定子空间方法得到全部结果用时179.28 s，而张量子空间方法仅用时61.79 s，耗时减少约2/3(计算机配置CPU均为I7-6700K，内存均为16 G)。以上结果表明，针对确定性系统识别问题，本文所提出的张量子空间方法较滑窗确定子空间方法能大幅提高计算效率，并能从信号中挖掘更多系统模态信息，且识别精度相同。