基于博弈K均值聚类的模态参数自动识别算法

郭 鹏, 李东升, 郭 鑫

(1. 汕头大学 智能制造技术教育部重点实验室 广东省结构安全与监测工程技术研究中心,广东 汕头 515063;2. 内蒙古大学 交通学院,呼和浩特 010030)

环境激励下的模态参数识别作为结构健康监测领域内重要的组成部分,受到工程界的广泛关注。由于模态参数识别理论大多要求为白噪声激励,而实测响应通常是在非白噪声激励下产生的,他们之间的差异往往导致虚假模态的产生。现有的运行模态分析技术在识别模态参数时,需要通过人为设定一系列阈值或参数来剔除虚假模态。然而,不同结构的参数设置情况往往是不同的,因此模态参数识别的结果会受到主观因素的影响。此外,在实际工程中,虚假模态的剔除需要经验丰富的专家消耗大量的时间和精力对识别结果进行人为定阶,模态参数不能实时自动识别,阻碍了结构健康监测的自动化。因此,开发和验证基于运行期间结构响应的模态参数自动识别算法显得至关重要。

随机子空间(stochastic subspace identification, SSI)可以非常有效地进行环境振动激励下的参数识别,具有操作方便,噪声稳定性好和识别精度高等特点,成为运行模态参数识别的主流方法。SSI方法主要有两种形式:一种是基于协方差的随机子空间(covariance-driven stochastic subspace identification, SSI-COV)算法;另一种是数据驱动的随机子空间算法。然而,不管哪种算法,都需要预先确定系统阶次,而实际待识别结构阶次未知,通常需要结合稳定图来确定模型阶次,定阶比较困难。为此,大量学者做了相关的研究。 Magalhaes等[1]以Infante D. Henrique桥为研究对象,利用SSI-COV结合层次聚类开发了一种模态参数自动识别方法,然而该方法需要人为设定阈值,具有一定的主观性。Zhang等[2]利用数据驱动的随机子空间算法,通过统计学的方法自动辨识了New Carquinez悬索桥的2年的模态参数,该方法对于虚假模态的讨论并不深入,且仍然需要人为选定各阶模态。为了解决模态参数识别中需要人为选定阶次的问题,Yang等[3]利用基于自然激励技术的特征系统实现算法结合模糊C均值聚类对某斜拉桥进行模态自动识别,该方法不需要用户设定参数,能够实时有效地对某混凝土斜拉桥的模态参数进行在线追踪辨识。尽管如此,Yang等提出的方法的前提是激励为白噪声,然而土木工程中并不能保证噪声是白色且平稳的,Qu等[4]指出了基于自然激励技术的特征系统实现算法的缺点,提出了一种基于虚拟频响函数结合特征系统实现的算法,实际结构表明其能够有效地识别各阶次频率,然而稳定图中仍然存在大量虚假模态。Reynders等[5]总结了量化模态参数不相似度的指标,并提出基于非层次聚类和层次聚类两阶段的稳定图自动分析方法,并对Z24桥进行了模态参数识别。考虑到Reynders等提出的方法至少需要一个用户定义的参数才能进行模态自动识别,且该参数需从实际数据中导出,Neu等[6]改进其算法,提出了三阶段聚类的稳定图分析方法,并将其应用于风洞试验数据中,得到了很好的结果。遗憾的是Neu等只在风洞试验进行测试,并未将其应用于实际工程中。与此同时,Sun等[7]引入了基于密度的聚类方法,进一步发展了Reynders等提出的两阶段聚类过程,并将其应用于某斜拉桥的模态参数自动辨识,但是该方法计算效率不高且结果易受噪声影响。Mao等[8]在克服密集模态的问题中,在Neu的基础上提出采用模态不相似度指标的第一主成分为聚类目标,以苏通大桥为研究对象,将20 d的数据用于分析,成功提取出模态参数,但仍存在聚类数及初值选取问题,其结果鲁棒性不高。Qu等[9]通过改进的频域分解法来区分密集模态,以含有密集模态的三自由度数值算例验证了该方法的有效性,由于其是在频域处理数据的,因此提出的方法精度依赖数据的频率分辨率。汤宝平等[10]将模态参数以及模态能量作为聚类因子计算其相似性,进而实现模态自动识别,然而其稳定图中仍然存在部分虚假模态。实际上,上述算法在剔除虚假模态时应用的都是传统聚类分析方法,并未考虑其物理意义。

为了解决上述方法中存在的聚类意义不明确的问题,本文在SSI-COV的框架上,结合虚假模态和结构模态的数量关系,分析了虚假模态产生的原因,提出了基于博弈K均值聚类的模态参数自动识别算法。该方法的核心思想是利用议价博弈理论对结构模态和虚假模态进行聚类,剔除虚假模态,再结合层次聚类将结构模态进一步细分为各阶模态。最后,通过数值算例和实际结构对本文提出的方法进行了较全面和严格的检验。

1 随机子空间算法和博弈K均值聚类算法基础

1.1 协方差驱动的随机子空间算法

一个动力学系统离散的状态空间方程[11]为

(1)

式中:xk为状态向量;Ad为离散系统矩阵;wk为过程噪声;yk为输出向量;Cd为离散输出矩阵;vk为测量噪声; 下标“k”为第k个时间步。

Hankel矩阵是随机子空间方法基础,构造Hankel矩阵有两种形式

(2)

式中:Yp为Hankel矩阵的前i行表示“过去”矩阵;Yf为Hankel矩阵的后i行表示“未来”矩阵;j为Hankel矩阵列数。

第二种形式的Hankel矩阵由Yf的第一行比Yp的最后一行延迟两个采样时间,形式如下

(3)

定义协方差矩阵,并假设过程噪声和观测噪声为零均值高斯白噪声,可得

(4)

利用式(4)计算得到的协方差,进一步构造Toeplitz矩阵

(5)

由第二种Hankel矩阵计算得到的Toeplitz矩阵为

(6)

将式(4)代入式(5),得

(7)

根据式(5)和式(6)之间的关系,结合式(7),可得

Θ2=OiAdΓi

(8)

对Toeplitz矩阵Θ进行奇异值分解求得可观矩阵Oi

(9)

根据式(7)和式(9)可得

(10)

(11)

最后,利用可观矩阵Oi的求得系统矩阵Ad

(12)

式中, (·)+为伪逆运算。对系统矩阵Ad进行特征值分解可以得到一个由极点组成的对角矩阵,从而能够进一步识别模态频率、模态阻尼比。在求系统的振型时,输出矩阵Cd为逆向可控矩阵Γi的前l行。

从式(2)可知,构造Hankel矩阵需要提前知道计算模型的阶次,目前仍然缺乏模型阶次确定的有效手段。为了回避模型阶次的问题,工程领域以稳定图法来替代,所以自动解释稳定图成为OMA(operational modal analysis)模态参数辨识的最主要问题。式(4)表明,随机子空间法要求激励满足零均值高斯白噪声,如激励为非白噪声,计算将会得到非系统阶次的模态,这是产生虚假模态的原因。

1.2 博弈K均值聚类算法

聚类分析是一种根据定义的逻辑规则对对象所拥有的特征进行分类或分组的技术。通常,处理后的数据具有高度的内部同质性和外部异质性。

给定n个样本X={x1,x2,…,xn},每个样本由一个特征向量Pi表示,特征向量的维数是m,K均值聚类的目标是将n个样本分到k个不同的簇中,并最小化样本和各聚类中心之间的平方差之和,根据这些假设,K均值聚类可以表示[12]为

(13)

式中:yik为指示函数, 当样本xi属于簇k时为1,否则为0;Ck为聚类中心;D(xi,Ck)为样本xi和Ck的距离函数。

K均值聚类是一种简单、快速的方法,然而该算法对初始聚类中心的选取较为敏感,且容易陷入局部最优解。为了使聚类中心之间的竞争能够吸引最大数量的样本加入对应的子集,Rezaee等[13]提出了博弈K均值聚类。该算法将聚类中心视为不同的参与者。这些参与者通过竞争将数量最多的样本吸引到其子集中。

根据Nash议价博弈理论[14],对于结构模态和虚假模态,可以看作两个参与者,其目标函数为

max(U-u0)(V-v0)

(14)

式中:U和V为结构和虚假模态的效用函数;u0和v0为效用崩溃点,即双方决定不议价的点。其表达式为

(15)

(16)

式中,β1和β2取决于结构模态和虚假模态的数量,在稳定图中,结构模态的极点数Np=OmaxL,虚假模态的极点数=Omax(Omax+1)/2-Np,Omax为人为给定的系统的最大阶次,L为研究频率范围内的模态数。d1和d2分别为待分类的某阶模态与结构模态聚类中心、虚假模态聚类中心最近的模态特征指标。

(17)

(18)

式中:α1和α2为常数,从统计学角度出发确定其值,本文取0.7;D1和D2分别为结构模态特征指标和虚假模态特征指标与对应聚类中心的集合,令

(19)

(20)

则,对于区分结构模态和虚假模态问题,博弈K均值聚类算法的目标函数可写成

(21)

为不失一般性,GBK聚类的目标函数为

(22)

其中

(23)

(24)

2 两阶段模态自动识别算法

为使剔除虚假模态算法具有通用性,本章以稳定图为基础,提出了两阶段模态自动识别算法。在第一阶段,应用SSI得到结构在各阶次下的模态参数特征,如频率、阻尼比、振型等特征,并对这些特征进行博弈K均值聚类,将其划分为结构模态和虚假模态。在第二阶段中,把剔除掉虚假模态的极点,应用层次聚类进一步划分为各阶结构模态,从而实现各阶模态的分离。

2.1 阶段一:基于GBK聚类区分结构模态和虚假模态

GBK聚类是基于样本集合划分的聚类算法,利用博弈论的思想将结构模态和虚假模态看作两个竞争者,通过目标函数最大化选取最优的划分。

通常情况下,在不考虑气动阻尼时,工程结构的阻尼比小于0.2。为了减少计算量,首先,剔除阻尼比在[0, 20%]以外的极点。考虑到结构模态在相邻阶次下表现为高相似性,构造相邻阶次模态偏差指标作为聚类的特征。模态的相似度可以用频率偏差、阻尼比偏差和模态置信准则(modal assurance criterion, MAC)偏差来描述,然而实际工程中布设的测点是有限的,采用MAC对振型相似度进行计算时,可能出现与不同频率相关联的振型系数的MAC值很大的情况,这就是空间混叠现象。针对这个问题,采用观测相似度(modal observability correlation, MOC)[15]来替代MACMOC(CMO)定义如下

(25)

在剔除虚假模态阶段需要考虑的第二个问题是,在实际工程中,结构的阻尼比并不是比例阻尼的形式,因此计算的振型在复平面中表现为不共线的复向量,如图1所示,并且结构模态的振型不共线程度是低于虚假模态的。在图1中,实线表示振型平均相位角。平均相位偏差(mean phase deviation, MPD)可以用于量化振型的不共线程度,MPD(DMP)定义如下

图1 复模态及平均相位Fig.1 Complex mode shapes and mean phase

(26)

(27)

式中:φji为振型;θ为振型平均相位角。

因此,聚类特征向量在本文可以定义为

Pi=[dfi,i+1dξi,i+1dCMO,i,i+1dDMP,i,i+1]

(28)

式中:dfi,i+1,dξi,i+1分别为相邻极点的频率偏差和阻尼比偏差;dCMO,dDMP分别为相邻极点的观测相似度和平均相位偏差,dCMO=1-CMO,dDMP=1-DMP。

考虑到各类偏差指标的概率分布可能不同,其中频率偏差远小于振型和阻尼比的偏差,直接将特征向量进行聚类会导致聚类效果不理想,因此,为了消除各指标偏差程度不一致对聚类结果的影响,需要在聚类前进行Box-Cox变换[16]将特征向量进行正态化。Box-Cox变换的一般形式为

(29)

式中,γ为变换参数,γ的估计一般有两种方法,分别是最大似然估计法和Bayes方法。

已知极点为结构模态和虚假模态的集合,因此将聚类的簇设置为2,设置聚类中心,利用式(17)和式(18)计算崩溃点d1,d2以及相关的参数。对聚类特征向量按照式(19)和式(20)进行迭代,获得新的聚类中心、距离函数,直到式(21)的目标函数不再增大;终止迭代找到最佳解,得到两个簇,选择样本特征小的作为结构模态,另一簇为虚假模态。

2.2 阶段二:利用层次聚类提取各阶模态

在2.1节中,应用GBK聚类将SSI-COV得到的极点聚类成结构模态和虚假模态。本节利用层次聚类将前一阶段提取出的结构模态进一步划分成各阶次模态,层次聚类提取各阶模态的流程分为以下4步。

步骤1初始时每个极点为一个类,计算各极点的距离,聚类距离公式为

(30)

式中,f为频率。

步骤2依据Single-linkage准则计算这两个类之间的距离,找到最小距离,并由此得到拥有最小距离的两个类的编号,将这两个类合并为一个新的类,根据式(30)重新计算每个类之间的距离。Single linkage准则定义的两个类之间的距离为两个类之间距离最近的两个点之间的距离。

步骤3重复步骤2,直到最终只剩下一个类,得到层次聚类树状图。

步骤4将步骤3生成的树状图用高度为ed的类间距离进行切割,把独立出来的每个分支设置为各阶模态。类间距离ed的计算公式为

ed=μp+2σp

(31)

式中,μp和σp分别为结构模态聚类特征距离的样本均值和标准差。换句话说,μp和σp为Pi(1)+1-Pi(3)求得的距离样本的均值和标准差。所提出方法的流程图如图2所示。

图2 两阶段模态自动识别流程图Fig.2 Flow chart of the two-stage automatic modal identification

3 数值算例

通过一个简单的六自由度质量-弹簧线性结构数值模型,验证该方法的有效性,如图3所示。该模型中各质量块均为2 kg,各弹簧刚度均为100 N/m,采用Rayleigh阻尼,系统的质量矩阵M,刚度矩阵K和阻尼矩阵C如下所示

图3 六自由度质量-弹簧线性结构Fig.3 Six degree-of-freedom linear mass-spring structure

在6个自由度结构上同时作用零均值Gaussian白噪声激励,利用Newmark法计算各测点(x1～x6)结构加速度响应。积分时间间隔为0.01 s,总时长为3 600 s。利用SSI-COV对数据进行分析,将模型最大阶次设置为100。

六自由度质量-弹簧结构的加速度功率谱图,如图4所示。由图4可知,各加速度响应的功率谱峰值在3 Hz以内,考虑到采样频率为100 Hz,为了减小高频噪声对识别结果的影响,预先对采集的加速度信号进行了低通滤波处理。

图4 六自由度结构加速度功率谱Fig.4 Power spectra of six degree-of-freedom structure

设定最大模型阶次为100,以最小模型阶次为2开始,频率-阶次关系如图5所示。其中,实线为加速度功率谱曲线,曲线峰值点对应的频率值为该结构的各阶次频率。离散点为利用SSI-COV方法不经过任何处理识别出的频率点。由图5可知:随着模型阶次的增加,图5中出现了与功率谱峰值点并不对应的频率点,这些频率点与系统无关,为虚假模态。当模型阶次大于20阶时,出现了大量与结构模态非常接近的虚假模态,如1.1 Hz,2.3 Hz,2.4 Hz附近的频率。如果没有功率谱峰值作为对比,很容易被误判为结构模态。

图5 频率-模型阶次图Fig.5 Frequency-model order diagram

将这些极点进行第一阶段聚类,剔除虚假模态,聚类特征为该模态与其最邻近模态之间频率、阻尼比、观测相似度及平均相位偏差。对剔除完虚假模态的剩余模态按照模型的阶次进行第二阶段聚类,通过计算类间距离将这些极点进一步划分成各阶次结构模态,得到模态参数。

对模态偏差指标特征向量进行GBK聚类并剔除虚假模态后的稳定图,如图6所示。由图6可知,剔除虚假模态后的稳定图比未处理的稳定图少了许多虚假模态的干扰,并且各阶频率清晰可见,均位于功率谱曲线的峰值点处。将层次聚类应用于剔除虚假模态后的极点数据后,自动提取出结构的各阶模态。

图6 两阶段聚类后稳定图Fig.6 Stabilization diagram after two-stage clustering

为了检验该算法的识别精度,将结果提出的方法和理论解进行比较。理论值和本文方法识别结果之间的对比,如表1所示。振型如图7所示,图7中星号为各测点计算的振型,实线两端为理论值,可以看出本文提出的方法所识别的结果与理论值非常接近,频率的识别误差最大仅为0.205%,说明本文能够准确识别出结构的模态参数,并且全程不需要人为干预,实现了模态参数自动识别。

表1 六自由度模型模态参数Tab.1 Modal parameters of the six degree of freedom model

图7 六自由度结构的振型系数Fig.7 Mode shapes of six degree-of-freedom structure

4 工程实例

为进一步验证所提出的基于GBK的模态自动识别算法的有效性,将利用某大跨径悬索桥的实际监测数据进行模态参数识别。该桥梁在2020年发生了一次涡激振动,专家组仔细分析原因,初步认为是更换吊杆的工人在桥面两侧设置水马改变了该桥梁的的气动外形,进而在当时的风况下的导致了涡激振动。但即使水马被拆除,振动仍未立即消失。专家进一步分析,认为是长时间大幅振动导致桥梁的结构阻尼显著降低。加速度传感器布置图,如图8所示。桥面共布设了14个竖向加速度传感器(桥面两侧各7个)和7个横向加速度传感器。各加速度传感器以200 Hz的采样频率采集桥梁振动信号。

本文以竖向加速度传感器作为分析样本,首先利用SSI-COV识别模态,然后利用两阶段聚类法剔除虚假模态,以获得各阶次模态的模态参数。

采集时间为1 h样本的7通道竖向加速度响应时程曲线及功率谱曲线,如图9、图10所示。

图9 竖向加速度响应Fig.9 Vertical acceleration responses

图10 竖向加速度功率谱Fig.10 Power spectrum of vertical acceleration

由图10可知,各加速度响应功率谱的主要峰值在0.6 Hz以内,采样频率远大于桥梁结构的主要频带范围。因此应对采集到的信号进行降采样,并且为了减小高频噪声对识别结果的影响,预先对采集的加速度信号进行了低通滤波处理。

传统只考虑阻尼偏差和频率偏差的方法剔除虚假模态得到的结果,如图11所示。由图11可知,随着频率的增加,依然存在与频谱图峰值点并不对应的频率点,如阶次为36左右,在0.1 Hz附近存在离散的未剔除干净的虚假模态,0.36 Hz附近也存在类似稳定的极点,实际上,根据功率谱所示,这些频率点与结构无关,为虚假模态,因此仅考虑模态偏差并不能完全剔除虚假模态。

图11 传统两阶段聚类后稳定图Fig.11 Stabilization diagram after two-stage clustering

结合MOC偏差和MPD偏差,重构聚类特征向量,并应用本文第2章所提出的两阶段聚类方法,结果如图12所示,环境噪声及模型误差导致的虚假模态被有效剔除,两阶段聚类后从竖向加速度数据中识别的前6阶频率和阻尼比如表2所示。

表2 某大跨桥梁模态参数Tab.2 Modal parameters of a long-span suspension bridge

图12 改进两阶段聚类后稳定图Fig.12 Stabilization diagram after two-stage clustering

由于缺少桥梁结构模态参数的理论解,无法验证辨识参数的精度,参考章关永等[17]的研究,与本文识别出的频率基本一致。由于测点数量有限,振型并不能完整描述,因此没有对振型进行计算。相比于1999年测试的结果,频率变化不大,然而阻尼比显著降低,这表明随着桥梁运营时间的增加,结构频率的变化远小于阻尼比的变化,因此在桥梁长期运营中,更应该关注结构阻尼比特性的变化。

5 结 论

本文提出了基于GBK聚类的两阶段模态参数自动识别算法,首先通过基于GBK聚类区分结构模态和虚假模态,然后利用层次聚类识别模态的阶次。通过数值模型进行了验证,并将其应用于某大跨桥梁实际监测数据的模态识别。

得出以下结论:

(1) 针对稳定图中提取结构模态时存在经验阈值这一问题,将博弈聚类技术应用于模态特征指标,利用结构模态和虚假模态之间的数量关系,可以将表征结构模态的特征和表征虚假模态的特征自动分为两类。本文所提出的两阶段聚类方法,可以让稳定图中的极点只出现在结构模态中。实桥算例对比分析表明,该方法能够很好的剔除虚假模态,并识别出结构模态参数。

(2) 在长期的结构健康监测中,结构频率的变化远小于阻尼比的变化,在大跨桥梁结构健康监测系统中,需要重点关注对结构阻尼比的监测和分析,及时采取措施,如增加阻尼器,以提高结构抗风的能力。

综上所述,本文提出的算法能够实现实际桥梁结构模态参数自动化识别,且精度高。并且通过工程实例说明了阻尼比监测的重要性。