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

2022-06-17 03:04单豪良陈永高孙泽阳

振动与冲击 2022年11期

单豪良，陈永高，孙泽阳

(1.浙江工业职业技术学院，浙江绍兴 312000； 2.东南大学土木工程学院，南京 211189)

随着桥梁健康检测系统[1]的不断发展和完善，现阶段，几乎所有大型桥梁结构上均布置了一套完整的检测系统。该检测系统的基本原理是，首先通过布置于桥梁结构上不同位置的传感器采集该结构的各种响应信号[2]；其次将采集的响应信号作为某种数学算法的输入以实现对该结构的模态参数识别[3]，获取对应的固有频率值、阻尼比及模态振型；最后对所得模态参数结果进行对比分析，以判别该结构的健康状态。可见桥梁结构健康检测系统的核心计算算法在于模态参数识别算法，不同的识别算法对应着不同的识别进度，以致直接影响结构健康状态的评估结果。

目前，已被广泛运用于桥梁结构的参数识别算法有两种，分别是随机子空间算法(stochastic subspace identification,SSI)[4]和确定随机子空间系统识别方法(combined determine-stochastic subspace identification,CDSI)[5]。通过对比两种算法的基本原理，可知CDSI算法是基于SSI算法所得，实际应用中因该算法不仅能够考虑结构的输入信号和输出信号，并且还能考虑系统的过程噪声[6]和测试噪声，在一定程度上弥补了SSI算法在实际运用中的不足之处。

随着CDSI算法被广泛运用于桥梁结构模态参数识别中，人们发现该算法和SSI算法一样存在缺陷，即需要识别者参与到稳定图[7]中模态的辨识，这一环节大大的降低了参数识别的效率，同时还会因为识别者的主观意识导致最终辨识出的模态带有一定的主观性。为了避免模态参数辨识过程中的“主观性”，汤宝平等[8]提出将固有频率、阻尼比、模态振型以及模态能量作为谱系聚类算法的聚类因子对真实模态进行智能化筛选；周思达等[9]提出对模态参数结果进行模糊聚类分析以实现有效模态的自动化筛选；Ubertini等[10]直接对稳定图进行谱系聚类分析以筛选真实模态。虽然这些自动化识别算法均能在一定程度上完成有效模态参数的智能化识别，但均是基于SSI算法建立的识别算法。

鉴于此，本文通过将基于密度的全局领域算法(density-based spatial clustering of application with noise,DBSCAN)[11]运用于CDSI算法所得稳定图中模态参数的辨识，避免人为参与带来的弊端，提高参数识别效率，进一步推进CDSI算法在桥梁结构模态参数识别中的实际运用。文中将首先分别简单介绍CDSI算法和DBSCAN聚类算法的基本原理，其次介绍DBSCAN聚类算法的大致实现步骤及如何将该算法与CDSI进行融合，最后将改进的CDSI算法运用于识别某大型桥梁结构的模态参数结果，并通过将所得结果与理论值进行对比分析，以验证所提改进算法的可靠性。

1 CDSI算法基本原理

确定-随机子空间方法[12]能够将系统的输入信号(u1,u2,…,uk)和输出信号(y1,y2,…,yk)用于构建特定的Hankel矩阵，并对该矩阵进行行空间投影，同时对投影矩阵作奇异值分解，进而获得可观测矩阵和状态序列的卡尔曼滤波估计[13]，确定最终的系统矩阵A，B，C，D以求得系统的固有频率值、阻尼比及振型。具体实现步骤如下：

步骤1图1为确定-随机子空间算法的系统图，通过使用状态空间模型来描述该系统图，可得式(1)如下

(1)

图1 确定-随机子空间算法的系统图Fig.1 Combined determine-stochastic system flowchart

式中:uk∈Rm和yk∈Rl分别代表系统的输入和输出数据;wk∈Rn×l和vk∈Rl×l分别表示环境激励和传感器自身带来的噪声数据;k表示该系统所处的具体时刻;xk∈Rn为系统状态向量;A∈Rn×n为状态矩阵；B∈Rn×m为输入矩阵;C∈Rl×n为输出矩阵;D∈Rl×m为中途传递矩阵。

步骤2定义系统的输入矩阵和输出矩阵，式(2)为该系统的输入矩阵，同理可得该系统的输出矩阵，其中i为马可科夫更新过程步数，j代表数据的具体长度。

(2)

步骤3根据以下式子定义系统的状态序列。

Xi=(xi,xi+1,…,xi+j-1)

(3)

步骤4基于输入和输出信号及状态序列，建立含系统矩阵的线性方程

(4)

步骤5基于矩阵A求解得到该系统的固有频率值和阻尼比，并结合矩阵A和矩阵C求得该系统的模态振型结果。

2 DBSCAN聚类算法

在利用CDSI算法对桥梁结构进行模态参数识别时，需要人工参与到最后稳定图中模态的辨识，为了避免这一环节带来的弊端，引入DBSCAN聚类算法，该算法是一种基于密度的聚类算法，通过将紧密相连的样本划为一类，进而得到一个聚类类别；同时也通过将所有各组紧密相连的样本划为各个不同的类别，进而得到最终的所有聚类类别结果。

选择该聚类算法的主要原因如下：

(1) 相比K-Means算法[14]而言，该算法无需输入类别数；同时该算法能够实现任意形状的聚类簇，而K-Means仅能实现凸的样本集聚类；

(2) 相比谱系聚类和模糊聚类算法[15]而言，该算法无需定义各项参数间的相对权重值；

(3) 该算法的主要优点在于：可对任意形状的稠密数据集进行聚类；聚类的同时可找出异常点，且对异常点不敏感；不受初始值的影响，即聚类结果不存在偏倚。

以下将首先分析DBSCAN聚类算法中的相关定义[16]，其次将详细介绍如何将该算法运用于CDSI识别算法中。

2.1 DBSCAN的相关定义

定义1：领域阀值(Eps)，可将其看作为半径数据，即以点p为中心对象，以Eps为半径的区域：

NEps(p)={q∈D|Dist(p,q)≤Eps}

(5)

式中:D代表需要聚类的数据集合;Dist(p,q)表示对象q距离中心对象p间的距离值。

定义2：距离Dist(p,q)，假定对象p为n维向量(p1,p2,…,pn)，q为n维向量(q1,q2,…,qn)，则两对象间的距离公式采用欧几里德函数求解，即：

(6)

定义3：核心对象，对于聚类数据集D，假定密度阀值MinPts，当存在对象p∈D，且满足式(7)时，则可定义对象p为核心对象。

|NEps(p)≥MinPts|

(7)

式中，|NEps(p)|为p在Eps领域内的对象数量。

定义4：密度直达，对于聚类数据集D，若存在q∈NEps(p)且满足式(7)，则称定义q是从p密度直达的。

定义5：密度可达，对于聚类数据集D，若存在链p1,p2,…,pn∈D，其中pi(0

定义6：密度相连，对于聚类数据集D，若对象o∈D，使得p从q是从o密度可达的，则定义p和q是密度相连的。

定义7：簇，对于聚类数据集D，当数据集合C是属于D的一个非空子集，且满足如下两个条件时，则可定义集合C为一个簇：

(1) 对任意q，若核心p∈C且q是从核心p密度可达的，则定义对象q∈C。

(2) 对任意p,q∈C，p和q是密度相连的。

定义8：噪声点，对于聚类数据集D，若p游离于任何一个簇，则将其定义为噪声点。

2.2 参数配置

在利用DBSCAN聚类算法对数据集进行聚类时，仅需要确定两个参数即可，分别是领域阀值(Eps)和密度阀值(MinPts)。现阶段，传统的DBSCAN聚类算法中的这两项参数，均需根据经验进行设置。为了避免根据经验定参数带来的弊端，以下将结合CDSI算法分析如何根据稳定图来智能化确定这两项参数，实现聚类分析和模态参数识别的一体化智能操作。

2.3 稳定图定阶

稳定图定阶的基本原理[17]是通过分析相邻两阶模态对应参数结果是否满足事先假定的容差来辨识系统阶次，计算式如下，式中MAC代表模态置信因子。

(1-MAC(j,j+1))×100%<Δψ(振型容差)

(8)

式中:fj代表第j阶的频率值；ξj代表第j阶的阻尼比。

根据式(8)可知，稳定图选择的辨识因子包括频率值、阻尼比以及振型。

2.4 DBSACAN与CDSI的结合

本节将详细介绍如何将DBSCAN聚类算法与CDSI算法进行有效的嵌套使用，并合理地确定DBSCAN聚类算法中两大参数值-领域阀值(Eps)和密度阀值(MinPts)，进而实现模态参数的高效识别。

步骤1通过事先布置在桥梁结构上的加速度传感器采集对应的输入信号和输出信号，假定单位时间t内的输入信号为u1,u2,…,un，输出信号为y1,y2,…,yn，其中下标n代表单位时间t内采集的数据长度。

步骤2采集多组单位时间(t1,t2,…,tN)对应的输入和输出信号数据，且各组时间长度均相同，其中N代表共采集的时间组数；

步骤3将各组输入信号(u1,u2,…,un)和输出信号(y1,y2,…,yn)作为CDSI算法的系统输入，并基于第1章所示算法流程识别得到共N组的模态参数结果，即得到N幅稳定图。假定每组参数结果为Xi={Fi,wi,Φi}(i=1,2,…,N)，其中F为频率值，w为阻尼比，Φ为模态振型，各参数结果均为m×n的矩阵，其中m数值大小等于1/2系统真实阶次。

Fi=fi(m,n)

wi=ξi(m,n)

Φi=Ψi(m,n)

(9)

步骤4基于Matlab软件捕捉各组稳定图各自对应的稳定轴数据(即稳定图中的频率值)，并保存；假定第i幅稳定图中的稳定轴为wi=(w1,w2,…,wj)，其中j代表对应的稳定轴的总个数；

步骤5基于DBSCAN聚类算法对第i组参数结果和第i+1组参数结果进行聚类分析，具体实现步骤如下：

(1) 选取第i组参数结果和第i+1组参数结果中对应的wi和wi+1作为初始的对象pj，j代表两组参数结果中所有稳定轴的总和数；

(2) 根据置信区间[-0.95,0.95]定义各对象pj对应的领域阀值(NEps(pj))；

(3) 根据稳定图定阶原理确定系统真实阶次为NN，即各稳定轴上对应的最大阶次数为NN。由于在利用稳定图定阶时，人们常常根据稳定数大于0.8NN来辨识该稳定轴是否为系统的真实稳定轴[18]；于此，可定义DBSCAN聚类算法中的密度阀值MinPts=0.8NN，即当密度超过系统真实阶次的80%时，可定义为一簇；

步骤6基于如下步骤实现对N组模态参数结果的聚类分析：

(1) 基于步骤5能够实现对第1组参数结果(X1)和第2组参数结果(X2)进行聚类分析，得到聚类结果X1-2；

(2) 以X1-2和第3组参数结果(X3)为DBSCAN聚类的输入数据进行聚类，得到聚类结果X1-3；

(3) 基于相同的原理，实现对N组模态参数结果的聚类分析，得到最终的聚类结果X1-N；

(4) 统计X1-N中各聚类模态是由多少组相应模态聚类而成的；并将聚类组数大于0.8N的模态作为有效模态绘制于稳定图中。

综合上述，可将DBSCAN聚类算法与CDSI算法进行有效地结合，得到改进后的D-CDSI模态参数识别算法，该算法的一般流程图如图2所示。

图2 D-CDSI识别算法流程图Fig.2 Flowchart of D-CDSI

3 实际桥梁结构的模态参数识别

本章将以某大型试验桥梁结构为研究对象，首先分别运用CDSI算法和SSI算法识别各工况下的模态参数结果，其次综合DBSCAN聚类算法以实现模态参数的智能化筛选，最终将两种算法所得结果与MIDAS软件所得理论值作对比分析，以验证所提D-CDSI算法较SSI算法而言具有更好的识别精度。

3.1 试验桥的概况

试验室中的桥梁结构模型是按照实际桥梁结构尺寸的1/20来设计的，其桥跨尺寸如图3所示。其中主跨长度为1 900 cm，边跨长度650 cm，全长3 200 cm，主梁上共布置11处加速度传感器用于采集主梁的振动响应信号，采样频率为256 Hz，即每秒采集次数为256次。主墩采用门式双柱结构，每个主墩的顶点处共设计12根斜拉索，各边跨处各设计两处桥墩用于支撑桥面系。该桥梁结构各构件的材料，主体结构采用M15微粒混凝土，包括索塔和桥墩；混凝土中的受力钢筋采用直径6 mm的圆钢；采用5 mm厚的钢板模拟主梁；斜拉索为直径为10 mm的钢丝。

图3 桥梁模型(cm)Fig.3 Bridge model (cm)

3.2 地震波信号

通过承台下方的振动台来施加相应的地震波，地震波的方向为横桥向，峰值加速度在0.02g～0.1g之间，每组地震波持续时间为1 min。表1为各工况下地震波的相关参数表，图4为峰值加速度为0.06g地震波在前20 s对应的加速度响应时程曲线图。

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

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

3.3 参数结果理论值

利用MIDAS软件建立该斜拉桥的模型，如图5所示。其中主塔采用钢筋混凝土材质，弹模取2.06×105N/mm2，泊松比取0.17；主梁和承力索均采用钢材材质，弹模取3.43×104N/mm2，泊松比取0.3；桥墩墩顶均采用固结；全桥模型节点数共计771，单元共计735。对其进行相关的模型修正，得到其在特征值分析下的前5阶模态振型和固有频率值，图6为该斜拉桥前2阶的模态振型图，表2为前5阶的频率理论值。

图5 3D实体模型Fig.5 Bridge model in three dimensions

(a) 一阶竖向弯曲振型(f1=1.13 Hz)

(b) 二阶竖向弯曲振型(f2=2.46 Hz)图6 前2阶模态振型图(MIDAS)Fig.6 Modal pattern of the first two orders (MIDAS)

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

3.4 各工况下CDSI算法识别结果

试验过程共输入5种工况对应的地震波，且各工况对应的试验时间为1 min，经大量试验后发现，可将每种工况分为12组，即每次以5 s对应的输入信号和输出信号作为CDSI识别算法的输入。图7为W1工况下0～5 s和5～10 s对应的稳定图。

(a) 稳定图(时间：0～5 s)

(b) 稳定图(时间：5～10 s)图7 CDSI识别算法所得稳定图(W1工况)Fig.7 Stabilization diagram under CDSI (working condition: W1)

3.5 各工况下SSI算法识别结果

根据3.4节可知，CDSI识别算法需要以每5 s对应的输入信号和输出信号作为参数识别的输入数据。而采用SSI算法识别时，仅需要将每5 s对应的输出信号作为参数识别的输入即可，图8为SSI识别算法所得W1工况下0～5 s对应的稳定图。

图8 SSI识别算法所得稳定图(W1工况,时间：0～5 s)Fig.8 Stabilization diagram under SSI (W1, time:0～5 s)

3.6 稳定图聚类分析结果

基于图2所示D-CDSI聚类算法的步骤实现对W1工况下共12组模态参数结果的聚类分析，得到图9所示聚类稳定图，根据该图可清晰的看到相关的模态都聚集到了一起，即所提算法能够在一定程度上实现同类模态参数间的聚类。

图9 D-CDSI聚类稳定图(W1工况)Fig.9 D-CDSI Clustering Stabilization diagram (W1)

基于图2所示D-CDSI聚类算法的步骤，将CDSI算法更换为SSI算法可以得到D-SSI聚类算法以实现对W1工况下共12组模态参数结果的聚类分析，得到图10所示聚类稳定图。由图10可知，系统的第一阶模态存在严重的遗漏现象。

图10 D-SSI聚类稳定图(W1工况)Fig.10 D-SSI Clustering Stabilization diagram (W1)

基于图2所示算法对D-CDSI聚类所得5种工况下共60组模态参数结果进行DBSCAN聚类处理，并采用各稳定轴的平均值为最终的频率值绘制最终的聚类稳定图，结果如图11所示。

图11 D-CDSI聚类稳定图Fig.11 D-CDSI Clustering Stabilization diagram

同时基于图2所示算法对D-SSI聚类所得5种工况下共60组模态参数结果进行DBSCAN聚类处理，并采用各稳定轴的平均值为最终的频率值绘制最终的聚类稳定图，结果如图12所示。

图12 D-SSI聚类稳定图Fig.12 D-SSI Clustering Stabilization diagram

对D-SSI聚类算法结果、D-CDSI聚类算法结果以及MIDAS软件所得理论值作对比分析，结果见表3。对比表中数据可知：

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

(1) D-SSI聚类算法无法识别出该斜拉桥的第一阶真实模态频率值；

(2) 对比D-SSI聚类算法结果与MIDAS理论值可知，最大误差百分比为9.2%，最小误差百分比为3.4%；

(3) 对比D-CDSI聚类算法结果与MIDAS理论值可知，最大误差百分比为6.6%，最小误差百分比为1.1%；

(4) 对比D-CDSI聚类算法结果与D-CDSI聚类算法结果可知，最大误差百分比为8.2%，最小误差百分比为1.1%；

综合上述可知所提D-CDSI算法所得频率值的聚类结果与MIDAS理论值间具有更高的匹配度。

3.7 模态振型结果

桥梁结构的模态参数识别主要包括三方面，分别是识别是结构对应的固有频率值、阻尼比以及模态振型的识别。表3已经给出了固有频率值的识别结果，就频率值的数值大小而言，可知本文所提算法的频率值识别结果具有可靠性。

考虑到桥梁结构的阻尼比结果易受环境的影响，以致其结果值经常出现较大的偏差，所以本文未对阻尼比进行识别。

图13为该斜拉桥前三阶模态振型图，其中横坐标中的11个点代表主梁上传感器的布置位置，纵坐标的数值[-1,1]代表归一化处理后该桥梁结构在竖向上的结构变形量。点划线为本文算法识别结果；实线为MIDAS所得振型图。图中前三阶模态振型的模态置信准则MAC达到0.96，可知所提算法识别所得振型图具有较高的准确性。

(a) 1阶振型图

(b) 2阶振型图

4 结论

桥梁结构的模态参数识别作为桥梁健康检测系统中的主要环节之一，参数识别的精确程度直接影响着桥梁健康评估的准确程度。本文针对现阶段被广泛运用的确定-随机子空间算法存在的不足，提出了将基于密度的聚类算法嵌入到该识别算法中，以得到改进的D-CDSI算法，进而提高参数识别的效率。同时将所提算法运用于某大型斜拉桥的模态参数识别中，结果表明：

(1) 所提算法能够有效避免人工参与辨识稳定图而带来的弊端；

(2) 通过将识别所得固有频率值和MIDAS对应的理论值作对比，可知识别结果与理论值很接近，误差百分比在6.6%以内；

(3) 将识别所得模态振型图和MIDAS对应的理论振型图作对比，可知两振型图间具有很高的相似度。

综合上述可知，所提改进D-CDSI算法能够运用于该试验桥梁结构的模态参数识别中。但因缺少实桥的相关数据，暂时无法将所提算法运用于识别实际桥梁结构的模态参数结果，所以该算法是否能够很好地识别出实际桥梁结构的模态参数结果还需进一步实证。