基于解析二阶盲辨识的结构工作模态分析及真假模态区分

李涧鸣，包腾飞,3，周喜武，高瑾瑾，顾 昊

(1.河海大学水文水资源与水利工程科学国家重点实验室，南京 210098；2.河海大学水利水电学院，南京 210098；3.三峡大学水利与环境学院，宜昌 443002；4.江苏省水利工程科技咨询股份有限公司，南京 210029；5.水利部交通运输部国家能源局南京水利科学研究院，南京 210029)

模态参数识别是结构健康监测领域的关键问题之一。对于大型土木和水利工程结构，通常采用工作模态分析(Operational modal analysis, OMA)技术进行模态识别[1]。不同于试验模态分析，OMA 无需施加人工激励，仅需环境激励下结构的振动响应识别参数，能更真实反映运行状态下的结构特性。

为实现结构服役安全状态的长期、实时和在线监测，需开发具备自动识别能力的OMA 方法。现有研究大多基于随机子空间方法(Stochastic subspace identification, SSI)和稳定图进行自动分析[2-3]。由于SSI 系统阶次冗余，且结构在运行中易受到非平稳激励和噪声干扰，会不可避免地引入虚假模态，往往需要借助专家经验才能从稳定图中分辨出真实模态[4-6]。为了减少人为因素，实现自动识别，国内外学者致力于研究通过聚类、指标阈值和智能算法等方法自动区分虚假模态和真实模态。

MAGALHAES 等[7]利用谱系聚类对稳定极点分组，将超过一定数量极点的分组自动归为真实模态。REYNDERS 等[8]基于谱系聚类和K 均值聚类提出了一种三阶段聚类自动识别方法，并通过诸多模态指标来判断真实模态。郑沛娟等[9]引入图论聚类通过二次聚类法识别结构真实模态。张永祥等[10]构造了2 个不同维度的分块Hankel 矩阵并进行同阶极点匹配得到更清晰的稳定图，进而通过谱系聚类进行自动识别。祝青鑫等[11]综合应用了主成分分析、K 均值聚类和谱系聚类提出了一种模态参数自动识别方法。虽然通过聚类等方法可以一定程度实现SSI 自动识别，但并不能根本解决稳定图难以分辨的困难，尤其对于激励相对较弱的结构，稳定图中极点往往分布混乱，聚类结果很难保证准确可靠。

鉴于上述原因，本文另辟蹊径，提出了一种基于解析二阶盲辨识的模态参数识别方法，并基于模态指标和K 均值聚类提出了相应的真假模态区分方法，以期实现模态参数自动识别。该方法避免了引入稳定图造成的模态选择困难，能够根据结构振动响应有效分离出各阶模态响应，通过构造反映时、频域模态响应特征的模态特征指标区分真实模态和虚假模态，进而利用K 均值聚类自动筛选出真实模态，并采用频域参数拟合法估计模态参数。通过8 自由度数值算例和混凝土重力拱坝工程实例验证本文方法的准确性和优越性。

1 基于解析二阶盲辨识的模态参数识别

1.1 二阶盲辨识算法

二阶盲辨识(Second-order blind identification,SOBI)[12]是一种原本用于信号处理的盲源分离技术，PONCELET 等[13]成功将该技术引入结构模态分析中，并证明了其适用性。

盲源分离模型假定实测信号和未知源信号满足线性混合关系，即对于n通道混合实测信号x(t)=[x1(t),x2(t),···,xn(t)]T，m组独立的未知源信号s(t)=[s1(t),s2(t),···,sm(t)]T，满足如下关系：

式中，A称为混合矩阵。BSS 的目标是根据实测信号x(t)同时估计混合矩阵A和源信号s(t)，满足：

式中，B为A的(广义)逆矩阵，称为分离矩阵。

SOBI 的基本思想是寻求使得实测信号的时延协方差矩阵近似对角化的源信号。定义时延协方差矩阵R(τ)=E{x(t)x(t-τ)T}，其中，τ为时间延迟系数，1≤τ≤p。采用近似联合对角化[12]技术求得矩阵J，使得JR(τ)JT近似为对角矩阵。则问题为：

式中，off(·)为矩阵的非对角元。

实际计算中，需进行白化处理。引入线性变换：

式中，W为白化矩阵，可由主成分分析确定。对协方差矩阵R(0)=E{x(t)x(t)T}进行特征分解：

式中，E和D分别为R(0)的正交特征向量矩阵和特征值矩阵。则W计算如下：

容易证明，E{z(t)z(t)T}为单位矩阵。白化后的时延协方差矩阵为R(τ)=E{z(t)z(t-τ)T}，通过Jacobi算法可求出正交矩阵J，满足式(3)。则混合矩阵A和源信号s(t)为：

对于N自由度线性系统，由模态叠加原理可知，结构振动响应y(t)与模态响应q(t)可以通过线性变换建立联系，即：

式中， Ψ为模态向量组成的模态矩阵。

可以看出，式(9)与式(1)具有相似的形式，模态矩阵 Ψ对应于混合矩阵A，模态响应q(t)对应于源信号s(t)。因此，可以通过SOBI 算法得到Ψ和q(t)，进而识别模态参数。

1.2 解析二阶盲辨识算法

由结构动力学可知，式(9)中只有当阻尼满足经典阻尼假定(如比例阻尼)的情况下 Ψ才为N阶实矩阵，故SOBI 理论上仅适用于阻尼简化的实模态情形。在振动测试中往往需考虑更一般的复模态情形，此时Ψ ∈CN×2N，q(t)∈C2N×1。因此，有必要将SOBI 推广以处理复模态问题。

由于模态向量与模态响应以共轭复数对的形式出现，则式(9)可表示为：

式中， (·)∗为复共轭。

与p(t)实部pR(t)对应的解析信号表示为：

与式(11)类似，x(t)对应的解析信号为：

则由式(12)可得：

相应地，可将式(1)推广为如下形式：

则1.1 节中的SOBI 可以自然推广为复数形式的解析SOBI。即根据式(15)，由白化处理后的yˆ(t)得到时延协方差矩阵，再由近似联合对角化求出酉矩阵Jc使得时延协方差矩阵近似对角化。类似式(7)和式(8)得到混合矩阵Ac与解析源信号sˆ(t)。相关过程与1.1 节完全类似，故不再赘述。

1.3 频域参数拟合法估计模态参数

利用式(14)和式(15)的相似关系，可由解析SOBI 估计模态矩阵 Φ和与模态响应q(t)相关联的解析信号pˆR(t)。由矩阵 Φ得到结构振型，由可通过频域参数拟合法估计固有频率和阻尼比。

根据振动力学，N自由度系统的运动方程为：

式中：u(t)为位移向量；f(t)为外部激励向量；M、C、K分别为质量矩阵、阻尼矩阵和刚度矩阵。系统的频响函数H(iω)具有如下形式：

式中：ωnr第r阶自振圆频率； ζr为阻尼比。

在频响函数的第r阶共振峰附近，忽略共轭项，则式(17)可作以下近似：

假设激励为白噪声，则pR(t)的自功率谱密度函数P(ω)可作为H(ω)的估计。由式(19)可得：

P(ω)谱峰附近取l个频率点ω1,ω2,···,ωl，则：

由式(21)可通过最小二乘法估计λr，则由式(18)

可得固有频率fr和阻尼比 ζr为：

2 基于模态特征指标和K 均值聚类的真假模态自动区分

2.1 模态特征指标

由于结构被充分激励的模态数目难以事先确定，在解析SOBI 分析中假定源信号的数目与实测信号维数相同，计算中会引入多余的虚假模态。另外，实测信号中还存在系统和测量噪声、非平稳激励以及滤波等引入的虚假模态。为实现模态参数自动识别，需要对结构真实模态和虚假模态进行自动区分，最后通过频域参数拟合法识别真实模态参数。通过比较模态表现出的不同模态响应特性，分别从时域和频域角度出发，提出了总包络长度指标和谱熵指标两个模态特征指标以区分真假模态。

2.1.1 总包络长度指标

在时域中，真实模态的模态响应通常具有更短的包络线[15]。根据解析信号理论，对于模态j，模态响应曲线qj(t)的包络线可表示为。由勾股定理， dt时段内的包络线长度近似为dl=则总包络线长度为：

2.1.2 谱熵指标

在频域中，真实模态的自功率谱密度函数通常更稳定，更有规律性。可以将其类比成概率密度函数。对于模态j，模态响应曲线qj(t)的自功率谱密度函数表示为Pj(ω)，为满足概率积分为1 的条件，可对其进行归一化：

在信息论中可通过信息熵度量信源的不确定性。则类似定义谱熵描述功率谱密度的不确定性：

2.2 K 均值聚类区分真假模态

对于解析SOBI 计算的各模态响应，分别计算总包络长度和谱熵两个模态特征指标，并归一化到[0,1]区间以均衡权重。对于模态j，其总包络长度和谱熵指标值构成二维数据点。本文采用聚类数为2 的K 均值聚类算法(K-means Clustering)自动将模态特征指标构成的二维数据点分成两个簇，分别代表真实模态和虚假模态构成的集合。K 均值聚类算法的目标是寻求使得各簇中所有数据点与最近的聚类中心的距离之和最小的最优聚类方案，可表示如下：

式中：S1和S2分别为真实和虚假模态构成的簇；c1和c2分别为S1和S2对应的聚类中心点；//·//2为L2 范数，表征欧式距离。算法步骤如下：

步骤1：对于各数据点组成的集合{p1,p2,···,pm}，随机选择2 个初始聚类中心点c1、c2。

步骤2：分别计算各数据点pj到各聚类中心ck的欧式距离//pj-ck//2，寻找与各数据点距离最小的聚类中心，并将数据点分配到对应的簇中。

步骤3：采用平均的方法重新计算分配完成后的新簇的聚类中心，计算公式为：

式中，nk为Sk中数据点的个数。

通过K 均值聚类可自动提取真实模态，依次采用1.3 节所述的频域参数拟合法即可估计出各真实模态对应的模态参数。

3 8 自由度质量-弹簧-阻尼模型数值算例

为验证本文提出方法的有效性，并检验其自动区分真假模态的能力，采用MATLAB 平台建立了一个8 自由度质量-弹簧-阻尼模型，其质量、刚度和阻尼矩阵分别为：

模型考虑非经典阻尼，阻尼矩阵可由满足瑞利阻尼的部分添加非经典项形成。其中瑞利阻尼部分的第1 阶和第8 阶模态的阻尼比设定为2%。对模型的各自由度施加高斯白噪声以模拟环境激励。采用50 Hz 采样频率进行信号采集，各自由度采集60 min 时长的加速度响应。对各加速度信号均添加信噪比为10 dB 的高斯白噪声以模拟真实测试环境中的测量噪声。结合模型的固有频率特点，对信号进行重采样，奈奎斯特频率为6.25 Hz。经过重采样，上述加速度响应信号中仅包含模型前4 阶模态信息。采用解析SOBI 进行模态识别时，假设源信号数等于信号通道数，因此会引入4 阶虚假模态。如2.1 节所述，对解析SOBI 识别的8 个模态响应分别计算总包络长度和谱熵指标，并采用K 均值聚类方法区分真假模态，如图1 所示。可以看出，K 均值聚类算法成功将各模态分成两个簇，其中4 阶真实模态(“+”表示)构成的簇比4 阶虚假模态(“×”表示)构成的簇的聚类中心更靠近坐标原点，可实现真假模态自动区分。

图1 模态特征指标的K 均值聚类结果Fig.1 The result of K-means clustering for modal metrices

根据K 均值聚类结果可得4 阶真实模态对应的模态矩阵及模态响应相关联的解析信号。由频域参数拟合法识别各阶固有频率和阻尼比，见表1。表1 中同时列出了相应理论值，并计算相对误差。固有频率和阻尼比的相对误差绝对值分别小于0.04%和6%。作为对比，采用增强频域分解法(Enhance frequency domain decomposition, EFDD)[16]对固有频率和阻尼比进行识别，见表2。可以看出，本文方法对固有频率和阻尼比的识别结果均具有较高精度，其中固有频率的精度与EFDD 相近，而阻尼比的精度比EFDD 更高。

表1 本文方法识别的模态参数与理论值的比较Table 1 Comparisons between modal parameters identified by the proposed method and theoretical values

表2 EFDD 识别的模态参数与理论值的比较Table 2 Comparisons between modal parameters identified by EFDD and theoretical values

由模态矩阵得到结构振型，可通过模态置信准则(Modal assurance criterion, MAC)评价振型识别结果。第r阶和第s阶模态的模态向量vr和vs之间的MAC 值定义为：

式中：(·)H为共轭转置；MAC 值越接近1，则振型越相关；反之越接近0，则越不相关。分别计算本文方法识别的各阶振型与相应理论值之间的MAC值，见表1。相应地，EFDD 识别的各阶振型与理论值之间MAC 见表2。可以看出，两种方法识别的各阶MAC 值均大于99.9%，具有较高精度。

4 龙羊峡大坝工程实例

4.1 工程概况

龙羊峡水利枢纽位于黄河上游，由混凝土重力拱坝主坝、两岸重力墩、两岸重力副坝、泄水建筑物、引水建筑物和水电站厂房等组成。挡水前缘总长度1226 m，其中主坝长396 m。最大坝高178 m，最大底宽80 m。建基标高2432 m，水库正常蓄水位2600 m。工程处于高地震烈度区，基本烈度为VII 度，大坝按VIII 度地震设防。大坝于1986 年布置了强震监测系统，并于2009 年进行了技术改造，实现了震动实时监测能力。监测系统采用TDA-23 型三向加速度传感器和TDA324 FA/CA 数字采集器，共布置13 个测点，采样频率为200 Hz。

4.2 大坝模态参数识别

采用主坝2600 m 廊道中布置的5 个加速度传感器共15 个监测通道的监测数据进行大坝模态识别，如图2 所示。采用2016 年6 月6 日11 点至12 点的60 min 加速度响应数据进行分析，如图3 所示。

图2 主坝2600 m 廊道加速度传感器布置图Fig.2 Positions of accelerometers in the 2600 m dam gallery

图3 大坝60 min 加速度时程曲线Fig.3 Dam acceleration time series in 60 min

对信号进行重采样，奈奎斯特频率为10 Hz。采用本文提出的方法进行模态识别，假设源信号数等于信号通道数15。对解析SOBI 识别的15 个模态响应分别计算总包络长度和谱熵指标，并采用K 均值聚类方法区分真假模态，如图4 所示。可以看出，K 均值聚类可自动区分真假模态。根据K 均值聚类结果可得5 阶真实模态对应的模态矩阵及模态响应相关联的解析信号pˆR(t)。由解析信号实部pR(t)计算得到各阶模态响应的自相关函数和自功率谱密度函数，分别如图5 和图6 所示。为简化起见，各模态响应自相关函数幅值归一化到[-1, 1]区间。可以看出，各阶模态时域和频域均符合结构响应特征，解析SOBI 具有较好的模态分离效果。由模态矩阵得到大坝振型如图7 所示，相应描述见表2。可以看出，第1 阶、第3 阶和第4 阶振型为对称振型，第2 阶和第5 阶振型为反对称振型。

图4 模态特征指标的K 均值聚类结果Fig.4 The result of K-means clustering for modal metrices

图5 各阶模态响应的自相关函数Fig.5 Autocorrelation functions of modal responses for each mode

图6 各阶模态响应的自功率谱密度函数Fig.6 Auto power spectral density functions of modal responses for each mode

图7 本文方法的振型识别结果Fig.7 Mode shapes results identified by the proposed method

根据频域参数拟合法，由pR(t)计算得到自功率谱密度函数P(ω)，各谱峰附近取7 个频率点ω1,ω2,···,ω7，并相应得到P(ω1),P(ω2),···,P(ω7)，通过式(21)求得各阶模态复极点 λr，见表3。最后由式(22)和式(23)得到各阶固有频率和阻尼比，见表4。作为对比，采用EFDD 对固有频率和阻尼比进行识别，见表4。由于第3 阶和第4 阶模态固有频率较接近，EFDD 难以识别第4 阶模态，其余模态参数与本文方法识别的结果相近。计算本文方法和EFDD 识别的各阶振型之间的MAC 值，见表4。可以看出，除了第3 阶和第4 阶密集模态外，其余模态振型的MAC 值均大于97%。可见，采用本文方法能有效识别结构的模态参数，尤其对传统方法较难识别的密集模态具有更好的识别效果。

表3 频域参数拟合法计算表Table 3 Calculation table of frequency domain parameter fitting method

表4 本文方法和EFDD 方法识别的模态参数Table 4 Modal parameters identified by the proposed method and EFDD

由SSI 方法得到稳定图如图8 所示。稳定图中的极点尤其是第3 阶和第4 阶密集模态附近分布不清晰，较难分辨真实和虚假模态，更难以进行自动识别。而采用本文方法避免了引入稳定图，各阶模态尤其是密集模态分离效果较好，较易区分真假模态，识别过程清晰，更适合进行模态参数自动识别。

图8 随机子空间的稳定图Fig.8 Stabilization diagram of SSI

5 结论

本文基于解析二阶盲辨识提出了一种结构模态参数识别方法，并相应提出了真假模态自动区分方法。通过解析二阶盲辨识分离模态响应，然后构造时、频域模态特征指标表征模态特性，并利用K 均值聚类区分真假模态，采用频域参数拟合法估计模态参数。通过数值算例和龙羊峡混凝土大坝工程实例验证了方法的有效性，得到以下结论：

(1) 解析二阶盲辨识具有较好的模态分离效果。分离出的各阶模态响应符合结构响应特征。

(2) 本文提出的总包络长度和谱熵两个模态特征指标可以有效区分真实和虚假模态，通过K 均值聚类可以自动筛选出真实模态。

(3) 本文方法对模态参数的识别结果具有较高精度。数值算例中识别的固有频率和阻尼比与理论值的相对误差绝对值分别小于0.04%和6%。识别的各阶振型与理论值之间的MAC 值均大于99.9%。本文方法识别的频率和振型结果精度与EFDD 接近，而阻尼比结果精度比EFDD 更高。工程实例中识别的大坝固有频率和阻尼比除第4 阶模态外均与EFDD 方法识别的结果相近，除第3 阶、第4阶模态外，本文方法与EFDD 识别振型的MAC 值均大于97%。

(4) 本文方法对传统方法较难识别的密集模态具有更好的识别效果。且可以克服由于引入稳定图造成的极点分布不清晰，模态难以分辨的问题。更适合进行模态参数自动识别。