基于HHT 的结构模态参数自动化识别方法和试验验证

2022-11-05何定桥

工程力学 2022年11期

何定桥，杨军

(清华大学土木工程安全与耐久教育部重点实验室，北京 100084)

结构模态参数是结构的固有属性，在地震荷载、风荷载、温度荷载的作用下结构出现损伤时模态参数会发生变化，可以在一定程度上反映结构的健康状态，因此，结构模态参数的识别是结构健康监测的重要内容[1-3]。

随着MEMS 芯片、无线通讯、人工智能和云计算等技术快速发展，基于无线智能传感器和云平台的结构健康监测系统开始出现。这类系统解决了传统建筑结构监测的传感器布线难、系统可维护性差、数据时效性差、适用场景有限等诸多问题[4-5]。通过无线智能传感器配合云平台可以实现建筑结构的全天候不间断监测。虽然结构模态参数识别方法已经非常成熟，种类繁多[6-7]，但与实验室中的模态参数识别方法不同，利用结构健康监测系统对结构进行长期监测的过程中识别算法不能出现主观的参数选择过程。基于HHT 的模态参数识别中多个步骤中均需要研究人员进行主观判断与筛选[8]， HHT 的第一步EMD 会产生虚假的IMF 分量，对虚假分量的识别与剔除往往依赖研究人员的主观判断[9]。真实的IMF 分量也会存在模态混叠现象，即单个IMF 包含结构多阶模态信息，如果直接对存在模态混叠现象的IMF 分量进行Hilbert 变换，会损失其中部分模态的有效信息。

本文基于HHT 提出了一种自动化结构模态参数识别方法，首先将K-L 散度与DNN 结合对 EMD产生的虚假IMF 分量进行了自动化识别与剔除，接着将SSA 与Butterworth 滤波器结合对产生模态混叠现象的IMF 进行了混叠模态分离，最后利用Hilbert 变换实现模态参数的识别。

1 理论基础

HHT 算法由HUANG 等[10]在1998 年提出，是一种分析非线性系统非平稳信号的方法，由经验模态分解EMD 和Hilbert 变换组成。

设x(t)为结构的监测数据，首先找到上包络线fmax(t)与下包络线fmin(t)，记m(t)为二者均值：

如果h1(t)满足以下两个条件，那么h1(t)为第1 个IMF：1)h1(t)极值点个数等于零点的个数或相差1；2)h1(t)上下包络线的均值为0。若不满足，则重复上述过程，假设经过k次重复得到的hk(t)符合以上两个条件，则hk(t)为第1 个IMF，记作c1(t)，得到第1 个IMF 后，将其从原信号中除去，得到剩余部分：

令x1(t)为新的x(t)，重复上述步骤，以此得到c1(t),c2(t), ···,cn(t)，直到余量函数为单调函数或者小于某一阈值时停止循环，此时便将原信号分解为一组IMF 与一个余量函数。理论情况下，分解得到的每个IMF 均包含结构的一阶模态信息，将IMF 进行Hilbert 变换构造复信号，再通过最小二乘法拟合复信号的幅频函数与相频函数即可得到结构的模态阻尼频率(小阻尼比下近似为模态频率)与模态阻尼比[11]。

传统HHT 结构模态参数识别无法直接应用于结构健康监测的系统，主要包括两个原因：1) HHT第1 步是将监测数据进行EMD 分解，得到多个IMF，部分IMF 包含结构的固有模态信息，为真实的EMD 分量，但部分IMF 不包含结构的固有模态信息，为虚假分量，对虚假分量的识别与剔除往往依赖研究人员的主观判断。2) 理想状态下单个IMF 应包含结构的单阶模态信息，但实际IMF可能会包含不同尺度分量，使后续的Hilbert 变换失去物理意义，因此，需要对存在模态混叠现象的IMF 进行分离。本文针对以上两个问题对传统HHT 结构模态参数识别方法进行了改进，提出了自动化的模态参数识别方法，使得基于HHT 的结构模态参数识别不包含人为干预的主观过程，可以有效应用于建筑结构的全天候不间断长期监测，整体流程可以用图1 表示。

2 自动化模态参数识别

2.1 数值模型介绍

采用一6 层平面混凝土框架作为数值算例，模型使用SAP2000 商业有限元软件建立，层高3 m，柱距6 m。梁均为高0.6 m，宽0.25 m 的矩形梁，柱为0.4 m×0.4 m 矩形截面柱，该模型为线弹性模型，楼面荷载3 kN/m2，混凝土采用C30，模型外观如图2 所示。模型前4 阶频率分别为6.27 Hz、18.38 Hz、29.27 Hz 和38.25 Hz。

2.2 EMD 虚假分量识别与剔除

经EMD 分解后得到的IMF 并不都反映结构的真实模态，由于采样率不足以及样条插值选取不适当会产生EMD 虚假分量，尤其是低频分量，导致分解结果不满足能量守恒。目前常见的EMD虚假分量识别方法包括相关系数法，然而IMF 分量间的相关系数较为接近，且判别过程较为主观，容易造成误判。PENG 等[12]提出通过WPT 将原信号分解为窄带信号再进行EMD；祁泉泉等[13]提出可以采用滤波器对原信号进行滤波再进行EMD分解；黄迪山[14]提出考虑虚假模态导致的能量不守恒。以上方法需预估结构自振频率的区间，若结构发生损伤自振频率变化时，结构自振频率会超出事先设定的截止频率，不适于结构的长期监测。

2.2.1 K-L 散度的应用

K-L 散度是一种量化两种概率分布之间差异的方式，又称为相对熵，可以用来衡量IMF 分量与原信号的“关系远近”。若EMD 的一个分量为真实分量，则其与原信号较为相似，两者之间的K-L 散度较小，反之K-L 散度较大[15]。

对数值模型从底部输入白噪声激励，时长6 s，PGA 为40 mm/s2，采样频率为500 Hz。数值模型在荷载激励下第6 层加速度响应经过EMD 分解得到的IMF 与原信号的K-L 散度如图3 所示。IMF1-3与原信号的K-L 散度分别为0.0274、0.0772、0.0823，均小于0.1，其余的IMF 为EMD 虚假分量，其与原信号的K-L 散度均大于0.5。

为了更好分离真实和虚假EMD 分量，真实分量的K-L 散度应该足够小，而虚假EMD 分量的K-L 散度应该足够大。在K-L 散度计算的核密度估计中，窗宽h的选择对计算K-L 散度精度的影响很大，窗宽过大或过小都会导致真实和虚假IMF的K-L 散度接近，从而影响分离效果。宋娜等[16]利用遗传算法寻找最优的h最大化分离度，但是文中预设了原信号的真实概率密度函数已知，这一点在自动化模态参数识别中难以符合要求，且K-L 散度阈值的仍然是主观确定的。对于不同的信号，因为采样频率、信号质量、噪声强度等影响，区分真实或虚假EMD 分量的阈值不同，在处理不同数据时需要对阈值进行调整，这一点在模态自动化识别中常常难以实现。

2.2.2 DNN 的应用

深度神经网络可以有效地应用于EMD 虚假分量的识别。将IMF 在不同窗宽下的K-L 散度及相关数据作为输入特征，将IMF 是否为虚假分量作为输出结果对DNN 进行训练。当获得新的信号，对信号进行EMD 分解后，DNN 可以自动识别该信号EMD 分解得到的IMF 中虚假分量，而无需人工选择K-L 散度的窗宽h和阈值[17]。

DNN 训练及验证样本为人工模拟结构的加速度信号经过EMD 后得到的IMF。仿真信号的确定需要参考待识别结构的动力特性与传感器特性，仿真信号的生成公式为：

该信号由n条频率为f1,f2, ···,fn幅值为A1,A2, ···,An的正弦信号与白噪声叠加而成。假设一条仿真信号经过EMD 后得到S个IMF。每个IMF 样本共提取21 个特征，第k个IMF 的特征1-10 为该IMF 与原信号在不同窗宽下的K-L 散度值：

式中：N为信号的长度；特征1 为窗宽最小(N/200)时该IMF 与原信号的K-L 散度；特征10 为窗宽最大(N/20)时该IMF 与原信号的K-L 散度。第k个IMF 的特征11-20 为该IMF 与原仿真信号在K-L 散度值的归一化排名：

第k个IMF 的特征21 为该IMF 与原信号在不同窗宽下的K-L 散度值的标准差(也即特征1-10的标准差)：

随机生成500 条仿真信号，经过EMD 后得到4309 条IMF 信号。DNN 采用有监督式机器学习，对每一条IMF 标记是否为虚假EMD 分量。将IMF 进行傅里叶变换并提取主频fIMF，与原信号包含的频率成分f1,f2, ···,fn进行对比，若fIMF与原信号包含的任一频率成分相差不超过5%，则将该IMF 标记为EMD 真实分量，否则标记为虚假分量。

经过试验，确定DNN 网络层数为6 层，包含1 个输入层、4 个隐藏层、1 个输出层，输入层神经元数量与样本的特征数量一致为21 个，隐藏层神经元数量根据经验公式估计为10 个，判别是否是虚假EMD 分量属于分类问题，因此，输出层神经元数量为1 个。选择Sigmoid 函数作为输出层的激活函数，隐藏层的激活函数选择ReLU。学习率和丢失率的设定通常需要经过多次试验确定，模型学习率LR 设置为0.1。丢失率的设置是为了简化网络结构，防止模型过拟合，提高容错率选择为5%。损失函数选择二元交叉熵函数，梯度下降算法选择小批量梯度下降法。

将70%的数据划分为训练集，30%的数据划分为验证集对模型进行训练，模型训练用时23.89 s。如图4 所示，模型在训练集上损失函数由0.628 下降到0.304，准确率达87.20%；在验证集上损失函数由0.576 下降到0.304，准确率达87.63%，在训练集和验证集上均表现出很高的准确度。

2.2.3 模型验证

提取数值模型在白噪声荷载激励下的加速度数据作为响应信号。为模拟实际情况，在信号中部分时间片段中加入干扰噪声。提取数值模型2、4、6、8、10 节点数据，添加5%、10%、20%、30%的噪声干扰，共生成20 条信号。将原信号进行EMD分解得到161 条IMF，使用DNN 对IMF 识别结果如表1。

表1 EMD 虚假分量识别准确率Table 1 EMD false component identification accuracy

在5%、10%、20%、30%的噪声干扰水平下，识别准确率达85.7%、81%、68.3%、60.5%，整体识别准确率达73.3%。在间段噪声干扰下，DNN 仍可准确识别EMD 的虚假分量。

2.3 混叠模态分离

结构监测信号EMD 后得到的IMF 通常不仅仅包含单阶模态信息，而有可能包含多阶模态信息。模态混叠产生的原因主要包括[18]：

1) 输入信号中混有间断的高频小幅值信号；

2) 信号中存在间段的噪声；

3) 信号中存在频率相近的分量。

模态混叠问题的解决方法包括集总经验模态分解[19](EEMD)、掩膜信号法[20]、独立分量分析方法[21](ICA)、解相关算法等。但是以上方法都具有一定不足之处，EEMD 需要预设信号的信噪比，这一点在模态自动化识别中难以实现；掩膜信号法中对不同的信号掩膜信号幅度和频率难以确定；ICA 的幅度不确定性会限制对信号能量信息的判断；解相关算法中的相关系数阈值需要主观确定，不适用于模态自动化识别。

2.3.1 分离步骤

SSA 可以根据时间序列构造出轨迹矩阵，并对轨迹矩阵进行分解、重构，从而提取出代表原时间序列不同成分的信号，分离出包含多阶模态信息的IMF 中的单阶模态信息。再利用Butterworth阻带滤波器将该阶模态信息从原IMF 信号中过滤，从余量函数中继续提取剩余的模态信息。直到所有的单阶模态均被有效分离。

具体的分离步骤如图5 所示：

1) 获取结构在白噪声荷载激励下加速度响应，进行EMD 后去除虚假分量，得到反应结构真实模态信息的IMF，假设该IMF 为第i个EMD 分量，记为IMFi，计算IMF 自功率谱；

2) 根据自功率谱的峰值判断IMFi包含的模态数量；

3) 若IMFi包含的模态信息多于1 阶(假设包含n阶)，则利用SSA 从IMFi中分离出自功率谱幅值最大的模态，记为IMFi.1；

4) 对IMFi使用Butterworth 阻带滤波器，过滤IMFi.1，得到ri.1，利用SSA 从ri.1中分离出自功率谱幅值最大的模态，记为IMFi.2；

5) 对IMFi使用Butterworth 阻带滤波器，过滤IMFi.2，得到ri.2；

6) 重复步骤4)～步骤5)，直至得到IMFi.1,IMFi.2, ···, IMFi.n。

2.3.2 模型验证

对数值模型节点11 在白噪声荷载激励下的加速度时程进行EMD 分解，去除EMD 虚假分量，得到前3 阶IMF 为EMD 的真实分量，原信号及前3 阶IMF 的自相关函数如图6 所示。

对上述自相关函数作FFT，得到它们的自功率谱如图7 所示，原信号包含了结构的前4 阶模态信息，理想状态下IMF1～IMF4应该分别包含结构的第4 阶～第1 阶模态信息，实际IMF1包含了结构的第2 阶、第3 阶和第4 阶模态信息，IMF2包含了结构的第1 阶和第2 阶模态信息，IMF3包含了结构的1 阶模态信息。

IMF1中第4 阶模态的自功率谱幅值最大。将IMF1进行SSA，IMF1经过SVD 分解后得到1200个特征值，将特征值排序后选择前80%的特征值进行重构得到IMF1.1，IMF1.1反应了结构的第4 阶模态信息，IMF1.1及其自功率谱如图8 所示。

使用阻带滤波器在IMF1中过滤IMF1.1的模态信息，可以得到r1.1。r1.1包含结构的第2 阶和第3 阶模态信息，按照相同的方法利用SSA 分离可以得到IMF1.2和IMF1.3。IMF1.2和IMF1.3的自相关函数和自功率谱如图9和图10 所示。

至此包含结构第2 阶～第4 阶模态信息的固有模态函数IMF1已经被分离为了IMF1.1、IMF1.2和IMF1.3，均只包含结构的单阶模态信息。包含2 阶模态信息的IMF2混叠模态分离方式与IMF1完全相同，在此不再赘述。

2.4 Hilbert 变换

最终得到的固有模态函数包括IMF3(第1 阶模态)、IMF2.2(第1 阶模态)、IMF2.1(第2 阶模态)、IMF1.2(第2 阶模态)、IMF1.3(第3 阶模态)、IMF1.1(第4 阶模态)。对IMF3、IMF1.2、IMF1.3、IMF1.1进行Hilbert 变换，并通过最小二乘法拟合幅频函数与相频函数，便可得到结构的模态频率与模态阻尼比。模型对于频率的识别精度较高，结构前4 阶频率识别最大误差为1.8%，最小为0.3%。模型对高阶模态阻尼比识别存在一定误差。模态识别结果及与SAP2000 有限元模型计算结果对比如表2 所示。

表2 模态参数识别结果Table 2 Modal parameter identification results

3 试验验证

3.1 振动台试验介绍

振动台试验由清华大学土木系纪晓东课题组设计并在中国建筑科学研究院地震模拟振动台上完成[22]。试验结构为3 层钢筋混凝土框架剪力墙模型，长边框剪方向为2 跨，长4.7 m，短边框架方向为1 跨，宽3 m，层高为2 m，总高6 m。模型总重47.4 t，材料为C30 混凝土与HRB400 钢筋。模型每层均有附加配重，第1 层、第2 层、第3层顶板配重分别为8.94 t、8.94 t、9.36 t。模型通过基础梁与振动台连接，模型外观如图11 所示。

如表3 所示，加载分为7 个大组共25 个工况。第1 组为双向白噪声加载，第2 组～第7 组加载顺序均为框剪方向的单向地震波加载、双向白噪声加载、框架方向的单向地震波加载、双向白噪声加载。第2 组～第5 组输入地震波为JMA Kobe NS，加速度峰值分别为70 cm/s2、200 cm/s2、400 cm/s2、620 cm/s2，第6 组、第7 组输入地震波为JR Takatori 00，加速度峰值分别为400 cm/s2、7620 cm/s2，白噪声加载加速度峰值均为50 cm/s2。

表3 试验工况加载次序Table 3 Test loading sequence

3.2 EMD 虚假分量识别

试验所用加速度传感器采样频率为200 Hz，根据Nyquist 定理，可以准确识别到的最高频率为100 Hz，通过Opensees 有限元计算显示该结构的第6 阶频率为12.82 Hz。目标是分别识别结构x和y方向的前3 阶频率，因为需要识别结构在地震荷载输入发生损伤后的模态参数，结构损伤刚度下降后周期会延长，自振频率会降低，模态将会在0 Hz～10 Hz 内更加密集，所以设定DNN 训练所用的数值仿真信号包含6 阶信号，最高频率为30 Hz。每条仿真信号采样时间为10 s，随机生成1000 条仿真信号，经过EMD 后得到6652 条IMF 数据。

(3)相关手续、权责不清楚。在存货模式中，物流等监管其中角色重要。但是法律上的职责界定存在模糊。具体职责确定，风险、收益等的情况。现实中存在监管企业承担了过大的风险，但是收益过低的情况。造成一种不匹配情景。

将每条IMF 进行傅里叶变换，若IMF 的主频与仿真信号包含的频率成分相差5%以内，则将该IMF 标记为EMD 真实分量，否则标记为EMD 虚假分量。神经网络特征的选取与本文第2 节相同：

1) 特征1-10 为该IMF 与原信号在不同窗宽下的K-L 散度值；

2) 特征11-20 为该IMF 与原仿真信号在K-L散度值的归一化排名；

3) 特征21 为IMF 与原信号在不同窗宽下的K-L 散度值的标准差。

根据IMF 的特征与标签建立训练数据集。神经网络超参数需要不断调整优化，最终参数选择如表4 所示。

表4 DNN 模型超参数选取Table 4 DNN model parameter selection

使用70%的数据作为训练集，30%的数据作为验证集对模型进行训练，模型训练批次为20 次(迭代次数20 次)，总用时27.48 s。如图12，在训练集上损失函数由0.631 下降到0.416，准确率达到82.42%；在验证集上损失函数由0.592 下降到0.413，准确率达到81.96%，模型在训练集和验证集上均表现出很高的准确度。

提取表3 中工况1～工况7 中每个工况最后一次白噪声激励下结构第3 层加速度信号，信号时长取10 s，对信号进行EMD 分解后提取IMF 的主频，使用DNN 进行分类，标记为“虚假”或“真实”模态，IMF 的主频及DNN 识别结果如表5 所示。

3.3 混叠模态分离

对工况1-1 下结构第3 层x方向加速度响应进行EMD，可以得到6 条IMF，其中IMF4～IMF7为EMD 虚假分量， IMF1～IMF3为EMD 真实分量。由图13(a)中原信号的自功率谱可以看出原信号包含了结构x方向的前4 阶模态信息。

IMF1主要包含了结构的第2 阶、第3 阶和第4 阶信息(图13(b))，IMF2和IMF3包含了结构的第2 阶和第1 阶信息(图13(c)、图13(d))。

使用本文提出的方法对IMF1进行混叠模态分离。如图14 所示，IMF1.1反应了结构的第3 阶模态信息，IMF1.2反应了结构的第2 阶模态信息，IMF1.3反应了结构的第4 阶模态信息。

因篇幅限制，本节只给出工况1 的x方向数据混叠模态分离步骤，工况1 的y方向、工况2-7的x与y方向计算过程与工况1 的x方向数据相同，过程将不再赘述。

3.4 Hilbert 变换

对IMF3、IMF1.2、IMF1.1和IMF1.3进行Hilbert变换构造复信号，得到幅值函数与相位函数，通过最小二乘法拟合幅值函数与相位函数，如图15～图18 所示，便可得到结构x方向的模态频率与模态阻尼比。

通过最小二乘法拟合幅值函数与相位函数，便可得到试验结果x方向的模态频率与模态阻尼比。第1 阶～第4 阶频率识别结果为1.64 Hz、7.10 Hz、13.49 Hz 与16.21 Hz，第1 阶～第4 阶阻尼比识别结果为10.6%、5.2%、3.3%、5.8%。

为了评价模态参数识别的准确性，将自动化模态参数识别方法的结果分别与有限元模态分析结果和频响函数法识别的模态参数进行对比。选用两种方法作为对比的原因是：1)结构的有限元建模比较精确的前提下模态分析结果可以看作模态参数的理论值；2)频响函数法利用了输入激励的数据，识别结果可以看作试验结构真实的模态参数。

振动台试验可以提取结构基础的加速度数据作为结构的输入激励并计算频响函数，通过实频曲线与虚频曲线可以识别结构的频率和阻尼比。图19 显示了频响函数法识别的结构频率和阻尼比。

表6 显示了自动化模态参数识别方法与其他两种方法得到模态参数的结果。自动化模态参数可以准确识别结构的前4 阶模态，与有限元模型的结果相比，前4 阶频率识别的误差分别为21.9%、0.1%、3.8%、12.8%；与频响函数法对比，前4 阶频率识别的误差分别为8.9%、0.0%、2.2%、7.4%，前4 阶阻尼比识别的误差分别为6.1%、0.9%、1.1%、4.1%，自动化模态参数识别方法在频率和阻尼比两方面识别均具有较高的精度。

表6 模态参数识别结果对比Table 6 Comparison of modal parameters identification results

3.5 自动化模态参数识别的应用

将经过验证的自动化模态参数识别方法应用于试验结构在工况1～工况7(共25 次加载)中地震激励输入作用下的结构，对发生损伤后的结构进行模态参数识别。图20～图22 中，加载次序为奇数时，输入激励为白噪声激励，白噪声为双向白噪声，双向PGA 均为50 cm/s2；加载次序为偶数时，输入激励为地震荷载，输入激励选用的地震波、PGA、加载方向均在图20～图22 中标出，也可见表3 中的试验工况。

在结构输入地震荷载后，前3 阶频率均出现下降。第1 次加载时(工况1-1)，结构完整，x方向前3 阶频率为1.64 Hz、7.10 Hz、13.49 Hz；y方向前3 阶频率为0.92 Hz、3.67 Hz、6.17 Hz。第25 次加载时(工况7-25)，结构x方向前3 阶频率为0.63 Hz、3.41 Hz、9.44 Hz，相较结构完整时分别下降62%、52%、30%；y方向前3 阶频率为0.55 Hz、1.49 Hz、3.34 Hz，相较结构完整时分别下降40%、59%、46%。结构构件在地震荷载下的损伤导致构件刚度以及结构的总体刚度发生下降，进而导致自振频率下降。