基于改进确定-随机子空间算法的桥梁结构模态参数识别

陈永高，钟振宇,2，何 杰

(1. 浙江工业职业技术学院 建筑工程学院，浙江 绍兴 312000； 2. 浙江大学 建筑工程学院，杭州 310058；3. 中铁二院工程集团有限责任公司，成都 610031)

随机子空间算法[1](Stochastic Subspace Identification,SSI)作为现阶段最为广泛使用的模态参数识别算法之一，能有效识别出环境激励[2]下桥梁结构的固有频率、阻尼比及模态振型。随着该算法的不断推广，其缺陷也逐渐凸显；虽然已有不少学者[3]对SSI存在的缺陷提出了相应的改进算法，但由于该识别系统的输入信号为环境激励(信号未知)，以致在模态参数的整个求解过程中均无法建立明确的激励-系统-响应模型，导致系统特征值解空间的扩大，使得最终的识别结果中存在模态密集和虚假模态的现象，降低了识别结果的精确度。

确定-随机子空间系统识别(Combined Determine-Stochastic Subspace Identification，CDSI)算法[4]是在SSI算法理论基础上建立的，CDSI是利用结构的输入信号(激励)和输出信号(响应)作为系统的“输入”进行参数识别，能够有效地弥补SSI算法的不足之处，提高参数识别的精确度。同时由于CDSI算法的特点在于能够考虑结构的输入、输出及噪声，所以得到了国外不少学者的青睐。Belleri等[5]将CDSI运用于混凝土框架结构，实现了地震激励下不同损伤状态下模态参数的准确识别；Bakir[6]将CDSI运用于识别地震激励下的框架结构，识别结果表明该算法不仅能有效识别出结构的固有频率，识别的模态振型同样具有可靠性；Huang等[7]将CDSI运用于多层钢框架振动台试验中，并通过对模态参数变化情况的辨识实现了对结构损伤状态的判定。

虽然SSI算法被广泛运用于桥梁结构的模态参数识别中，但CDSI算法却几乎没被运用于识别桥梁结构的模态参数。主要原因归于：CDSI的核心算法是线性算法[8]，并不适用于像桥梁结构这种非线性时变结构。基于此，提出将滑窗技术[9]嵌套于CDSI算法中，以实现对桥梁结构进行时变模态参数跟踪识别；同时为了避免人为参与真假模态的区分，提出了将模糊C均值聚类算法(Fuzzy C-Means，FCM)[10]运用于模态参数结果的智能化筛选，以完成模态参数的在线跟踪智能化识别。

1 CDSI算法基本原理

图1是CDSI算法的基本流程图，其中：uk为系统的输入信号；B为输入矩阵；wk为环境因素带来的各种过程噪声；xk为系统状态向量；A为状态矩阵；C为输出矩阵；vk为因传感器带来的测量噪声；D为中途传递矩阵；yk为系统的输出数据。

图1 确定-随机系统框图Fig.1 Combined determine-stochastic system flowchart

得到确定-随机状态方程如式(1)所示

(1)

步骤1根据输入信号和输出信号建立Hankel矩阵，其中：Up作为过去的输入数据；Uf作为将来的输出数据。

(2)

步骤2确定卡尔曼滤波状态序列。

Kk-1=(G-APk-1CT)(Λ0-CPk-1CT)-1

Pk=APk-1AT+(G-APk-1CT)×

(Λ0-CPk-1CT)-1(G-APk-1CT)T

(3)

式中：K为Kalman滤波增益矩阵；P为向前状态矩阵。

(4)

得到卡尔曼滤波状态序列为

(5)

步骤3建立含有状态矩阵(A)、输入矩阵(B)、输出矩阵(C)及中途传递矩阵(D)的线性方程式。

(6)

步骤4对矩阵A进行特征值分解。

(7)

步骤5基于式(8)求解结构的固有频率值(f)、阻尼比(w)及模态振型(Φ)。

(8)

基于步骤1～步骤5便能实现桥梁结构的模态参数识别，其中就如何确定状态序列、矩阵A的特征值分解以及各步骤中各参数的具体含义可见文献[11]。

2 滑窗技术的运用

2.1 加窗原理分析

对响应信号进行模态参数识别时需控制信号的长短以提高参数识别的效率，实际运用中可通过对信号作加窗截断处理以控制信号的长短。对信号作加窗处理时，因窗函数会直接响应信号的频谱，且不同的窗口形式产生的频谱泄漏情况不同，以致最终的频率分辨效果不同。基于此，如何选择窗函数的种类、窗口大小及滑窗步长显得尤为重要，以下将依次详细论述如何确定这三项参数。

2.2 窗函数选择

现阶段常用的窗函数[12]包括矩形窗、汉宁窗、海明窗、平顶窗、高斯窗等。对比各窗函数的优缺点，最终选择采用矩形窗作为本文响应信号的加窗函数。原因如下：矩形窗作为时间变量的零次幂窗，具有主瓣集中的优点；同时相比其它窗函数，矩形窗能精确识别出频率结果，但识别的幅值结果精确度不高。在实际运用中，如果仅需精确地识别出主瓣的频率，而无需在乎幅值的识别精度，比如识别结构的自振频率，则可以优先选择矩形函数对信号进行加窗处理。矩形窗窗函数的表达式为

w(t)=1

(9)

2.3 窗口大小确定

每次输入系统的初始运算数据点个数-窗口大小会在一定程度上影响识别结果的精度。当窗口太小时，会因输入信号不足导致无法识别出准确的稳定图，且容易造成模态遗漏现象。当窗口太大时，会因输入信号过多导致稳定图中存在过多的虚假模态，且会降低参数识别的效率。为了确定窗口的大小，叶飞等[13]认为应该综合考虑结构自身的特点、信号的采样频率及信号含噪程度等因素。但笔者认为根据这些因素来综合考虑窗口的大小不具可行性，首先是因为涉及的因素太多，其次是因为各结构都具有自身的独有特性。基于此，可以借鉴稳定图定阶的原理来定义窗口的大小，具体分析步骤如下：

步骤1基于随机子空间法中的稳定图定阶[14]确定结构的真实阶次为N。

步骤2假定初始窗口数据长度为L1，对应窗口时间为t1s，采样频率为cf，三者之间满足L1=t1×cf。

步骤3假定前后窗口时间差均为Ts，即第i个窗口对应的数据长度为Li=[t1+(i-1)×T]×cf。

步骤5分析BFi数值大小随信号长度增加的变化情况，可知该百分比呈现先逐渐增加最后保持平稳的变化趋势，如图2所示。认为当百分比超过90%时对应的信号长度为合适的窗口大小。

图2 稳定点百分比趋势图Fig.2 Percentage trend chart of stable points

2.4 窗口步长

窗口大小确定后，则需进一步确定窗口的滑动步长，即每次窗口向前推进的数据点数，实现对输入信号的实时更新。假定信号每次更新的数据量为V，则各窗口间数据的更替、系统的输入及输出间的关系图如图3所示。

图3 滑窗步长示意Fig.3 Schematic diagram of step length

根据Hankel矩阵的定义[15]可知，过去的输入信号长度等于将来的输出信号长度，基于此，将定义为1/2窗口对应的的数据量，即窗口间的重叠部分为1/2窗口。

3 FCM聚类算法的运用

无论是SSI算法还是CDSI算法均涉及人为筛选稳定图中的有效模态[16]，为了避免因人为参与有效模态辨识而带来的主观性，提出了基于FCM聚类算法的有效模态智能筛选法。以下将首先介绍FCM聚类算法的基本原理及步骤，其次介绍如何将FCM算法和CDSI算法进行结合以实现有效模态的智能化辨识。

3.1 FCM聚类算法的基本原理

FCM聚类算法[17]相比K-means聚类和谱系聚类等而言，其优点在于：该算法的聚类原理是依据数据点的隶属度来辨识其分类情况，是一种模糊性聚类算法[18]；且能实现不同维度数据间的聚类。

利用FCM算法将数据集X分为c类，假定X由N组样本构成，即X={x1,x2,…，xN}；其中每个样本由n维向量构成，即xj={xj1,xj2,…,xjn}。定义数据集X的隶属度矩阵为U=[uij]，U为C×N的矩阵，其中C为分类个数，uij为X中的第j个样本归属于第i个类的隶属度。uij需满足条件

(10)

其中,uij的计算公式为

(11)

式中：m为模糊指数，m∈[1,∞]；dij为数据对象xj和第i个聚类中心(vi)间的欧氏距离,dij∈‖xj-vi‖；{v1,v2,…，vi}为聚类中心集合，其计算公式为

(12)

FCM算法的目标函数值为J(u,v)，其迭代过程的目的即是找寻到使J(u,v)函数为最小值的隶属度矩阵和聚类中心。

(13)

FCM聚类算法的具体步骤如下：

步骤1确定分类数目c和最大迭代次数m(模糊指数)。

步骤2从数据集X中随机选取一组数据为初始聚类的中心点。

步骤3根据式(11)求解隶属度矩阵U。

步骤4根据式(12)求解出聚类中心v。

步骤5根据式(13)找寻使J(u,v)函数为最小值的隶属度矩阵和聚类中心，并完成聚类。

3.2 FCM与CDSI的结合

为了实现稳定图中有效模态的智能化辨识，提出将FCM算法与CDSI算法进行结合以得到新的参数识别算法(C-CDSI)，该算法的具体步骤如下所示：

步骤1利用CDSI算法识别得到各窗口对应的参数结果Xi={Fi,wi,Φi}，其中F为频率值，w为阻尼比，Φ为模态振型，各参数结果均为m×n的矩阵，其中m数值大小等于1/2系统真实阶次。

Fi=fi(m,n)
wi=ξi(m,n)
Φi=ψi(m,n)

(14)

步骤2频率聚类。

对Xi={Fi,wi,Φi}和Xj={Fj,wj,Φj}分别在频率、阻尼比及振型三项参数上进行FCM聚类处理——①确定聚类数据集，假定对f1(m,n)与f2(m,n)中的两两频率模态进行聚类；②确定聚类分类数，因每次聚类的数据集只有2阶频率模态，则设定分类数c=2；③确定最大迭代次数，因参数结果矩阵的行数均为1/2系统真实阶次(N)，则设定最大迭代次数m=N/2；④首先对f1中第1阶频率结果f1(m,1)和f2中的第1阶频率结果f2(m,1)进行聚类，当两者为同一类时，则取平均值作为新的第1阶频率结果，即f1-2(m,1)=(f1(m,1)+f2(m,1))/2；依次类推完成两频率矩阵中所有模态间的聚类，并将所有聚类后的新模态与未聚类的模态构成新的频率矩阵f1-2。

步骤3阻尼比聚类。

为了验证频率矩阵f1-2中聚类的2阶模态具有真实性，则对ξ1(m,n)和ξ2(m,n)中相应的2阶模态进行阻尼比聚类辨识——①假定f1-2中f1(m,1)和f2(m,1)是同类项，则基于步骤2的算法流程完成对ξ1(m,1)和ξ2(m,1)的聚类；②当两者为同类项时，取平均值作为新的1阶阻尼比结果,ξ1-2(m,1)=(ξ1(m,1)+ξ2(m,1))/2；③当两者不为同类项时，则将f1-2(m,1)拆分为最初的f1(m,1)和f2(m,1)，并重新构建新的频率矩阵ff1-2；④将所有聚类后的阻尼比与未聚类的阻尼比构成新的阻尼比矩阵ξ1-2。

步骤4模态振型聚类。

为验证ff1-2中频率聚类和ξ1-2中阻尼比聚类的结果具有真实性，则可对ψ1(m,n)和ψ2(m,n)中相应2阶模态进行模态振型聚类——①假定ff1-2中f1(m,1)和f2(m,1)是同类项，则基于步骤2中频率聚类的步骤对ψ1(m,1)和ψ2(m,1)进行振型聚类辨识；②当两者为同类项时，取平均值作为新的1阶模态振型结果：ψ1-2(m,1)=(ψ1(m,1)+ψ2(m,1))/2；③当两者不为同类项时，则将f1-2(m,1)分为最初的f1(m,1)和f2(m,1)，并重新构建新的频率矩阵fff1-2；同时将ξ1-2(m,1)分为最初的ξ1(m,1)和ξ2(m,1)，并重新构建新的阻尼比矩阵ξξ1-2；④最后将所有聚类后的振型与未聚类的振型构成新的振型矩阵ψ1-2。

步骤5将fff1-2,ξξ1-2和ψ1-2组建为新的一组参数结果X1-2={fff1-2,ξ1-2,ψ1-2}，并基于步骤2～步骤4完成X1-2与X3间的同类模态聚类，得到新的聚类结果X2-3；以此类推，依次完成所有窗口间的参数结果聚类，假定窗口总数为N，则最后的聚类参数结果为XN-1-N。

步骤6为筛选出XN-1-N中有效模态(频率、阻尼比及振型均稳定的模态)，可统计XN-1-N中各聚类模态是由多少窗口中相应模态聚类而成的；并将聚类窗口数大于0.8N的模态作为有效模态绘制于稳定图中。

综合上述可知，可将滑窗技术、FCM聚类算法及CDSI识别算法融汇到一起构建新的模态参数在线智能化识别算法(SC-CDSI)，该算法的具体流程图如图4所示。

图4 SC-CDSI识别算法流程图Fig.4 Flowchart of SC-CDSI

4 振动台试验桥

为了验证前文所提改进确定-随机子空间算法(SC-CDSI)的可靠性，以某大型斜拉桥振动台试验为背景，分析其在地震激励下的模态参数变化情况。

4.1 工程概况

该试验桥为斜拉桥，按照1∶20的比例对实际桥梁进行缩小，实验室中桥梁跨度为(6.5+19+6.5) m。图5为该桥梁结构的桥型布置图、主梁的平面图，其中索塔和桥墩采用M15微粒混凝土模拟；普通钢筋采用直径6 mm的圆钢模拟；抗剪钢筋采用10#钢丝模拟；主梁和钢箱梁均采用厚度为5 mm的钢板模拟；拉索采用直径为10 mm的钢丝绳模拟。

4.2 加速度传感器位置

该振动台试验，在基础底部共安装4处加速度传感器以监测结构的输入信号；在主梁上共安装11处加速度传感器以收集结构的响应信号。图6为传感器的布置图。

4.3 地震输入

输入激励的各项参数见表1所示，可知峰值加速度均为0.1g，地震波的方向均为横桥向，波形分别为Chichi波和场地波。图7为C1工况下Chichi地震波的加速度时程曲线，传感器的采样频率均为256 Hz。

图5 桥梁模型立面图及平面图(cm)Fig.5 Plans and elevations of bridge model(cm)

图6 加速度传感器布置图(cm)Fig.6 Arrangement of acceleration sensors(cm)

表1 试验工况表Tab.1 Working conditions of test table

图7 加速度时程曲线(C1)Fig.7 Time-curve of acceleration(C1)

4.4 MIDAS有限元结果

采用MIDAS 2017建立该桥梁结构的仿真分析模型。桥墩、索塔以及主梁均采用非线性梁柱单元模拟，拉索采用桁架单元模拟；支座均采用盆式橡胶支座，且采用零长度单元模拟。全桥模型如图8所示。对该桥作特征值分析，获得其前5阶频率值，如表2所示，图9为前4阶模态振型图。

图8 全桥模型(MIDAS)Fig.8 Bridge model(MIDAS)

表2 自振频率Tab.2 Natural frequency Hz

图9 前4阶模态振型图(MIDAS)Fig.9 Modal pattern of the first four orders (MIDAS)

4.5 各窗口参数识别

根据滑窗技术和CDSI算法识别该斜拉桥在两种地震激励下的模态参数，设定初始窗口时间为1 s，试算出合适窗口大小为4 s；同时对应的滑窗步长为2 s，即每隔2 s更新一次窗口数据。对各窗口的信号进行参数识别，得到两种工况下各窗口对应的稳定图结果，图10为C1工况下前两幅稳定图结果。

图10 稳定图(C1工况)Fig.10 Stable diagram (working condition: C1)

识别出C1和S1两种工况下频率值在各窗口中的具体数值并绘制前5阶频率值随时间的趋势图，结果如图11所示。分析图中各阶频率的变化情况，可知两种工况下桥梁结构在前5阶的频率值均处于稳定的状态，即桥梁结构在两种地震波作用下，其自身的固有频率具有稳定性，也表明该斜拉桥在两种工况下未产生损伤。

图11 频率走势图Fig.11 Frequency chart

4.6 有效模态智能辨识

基于3.2节所提FCM聚类算法对两种工况下各窗口参数结果进行聚类分析，实现稳定图中有效模态的智能化辨识，图12为C1和S1工况下对应的有效稳定图。

图12 有效稳定图Fig.12 Effective stable diagram

将MIDAS有限元结果、C1和S1工况下前5阶频率结果作对比分析，结果如表3所示。根据表中两两结果间的差值百分比，可知SC-CDSI算法所得有效模态频率值与有限元结果间的误差在7%以内，即所提SC-CDSI算法的频率识别结果具有可靠性。

表3 频率对比分析结果Tab.3 Comparative analysis results of frequency

4.7 前3阶模态振型结果

为验证SC-CDSI算法不仅能精确地筛选出稳定图中的真实模态，还能识别出结构的模态振型。得到该斜拉桥的前3阶振型图，如图13所示，将其与图9中的前3阶振型图进行对比，可知，SC-CDSI算法所得模态振型与有限元结果具有90%以上的相似度，即SC-CDSI算法所得振型图具有可靠性。

图13 前3阶振型图Fig.13 Modal pattern of the first three orders

5 结 论

考虑到确定-随机子空间算法并不适用于桥梁结构这种时变结构，提出了将滑窗技术运用于CDSI算法中实现对桥梁结构模态参数的实时跟踪识别；同时考虑到需人为参与稳定图中有效模态的辨识，提出了将FCM运用于辨识稳定图中的真假模态，以实现模态参数结果的智能化筛选。

最后以某大型斜拉桥振动台试验为背景，将所提SC-CDSI算法运用于识别该桥的模态参数，并将所得结果与MIDAS有限元结果进行对比分析，结果表明，所提SC-CDSI算法不仅能精确地识别出桥梁结构的固有频率值，还能精确地识别出结构的模态振型结果，能够实现试验桥梁结构模态参数的在线智能化跟踪识别。因缺乏实际地震的输入信号和输出信号，所以该算法在实桥中的运用还需进一步分析论证。