基于经验模态分解和模型缩聚的动力响应重构方法研究

邹云峰，付正亿，何旭辉，卢玄东，阳劲松，周 帅

(1.中南大学土木工程学院，湖南，长沙 410075；2.中南大学交通运输工程学院，湖南，长沙 410075；3.中国建筑第五工程局有限公司，湖南，长沙 410004)

随着新技术、新材料在土木工程结构中的应用，人们对土木工程结构的安全性及可靠性提出更高要求。结合安装在工程结构上的传感器采集所得数据，可对工程结构进行健康监测及使用寿命评估[1-3]，但结构健康监测及使用寿命评估需要解决的关键问题是合理地监控关键位置的健康状态[4-6]，因此需要依赖于传感器系统。而大多数工程结构有限元模型自由度数目过多[7-10]，此外，鉴于工程结构的几何复杂性及部件多样性，某些位置(例如结构交界面、狭缝等)不易安装传感器，而这些位置往往是结构健康监测及使用寿命评估的关键部位[11-13]。

基于上述问题，近年来通过已安装传感器位置的响应对未安装传感器位置响应进行重构的方法得到广泛重视。Ribeiro 等[14]针对多自由度系统提出一种基于广义传递矩阵的响应重构法，在频域内建立两组响应的传递矩阵，进而实现响应待测点的重构工作。Law 等[15-17]基于频响函数将上述位移传递率推广至加速度传递率，完成了全模型或模型子结构上的响应重构。在进行子结构基础上的响应重构时，将子结构间的界面力视为外力并建立传递函数，从而完成响应重构。以上工作为频域上的响应重构，而完成频域基础上的响应重构工作时，需要多次运用傅里叶变换及逆变换，就计算效率而言，时域上的重构方法更具优越性。为此，Zhang 等[18-19]通过模态分析法，在时域基础上建立两组响应的模态转换矩阵，并结合归一化方法，完成了结构应变、位移响应重构。He 等[20-21]拓展了该重构方法，提出一种基于经验模态分解(Empirical Mode Decomposition,EMD)的响应重构方法，该方法通过间歇准则与经验模态分解法，得到测得响应的各单频时域响应，结合模态传递矩阵完成各单频响应的重构，再通过模态叠加法得到未安装传感器位置的动力响应。由于该方法不能完成应力及应变的重构，He 等[22]又基于欧拉伯努利梁理论，结合悬臂梁模拟及试验，进一步实现了应力及应变的重构。但现有的基于EMD 的时域响应重构方法，需要获取整个结构的全局刚度矩阵与质量矩阵。土木工程复杂结构的自由度数目较大，在提取刚度矩阵、质量矩阵时需要大量的计算资源与计算时间。同时，由于矩阵维数过大，也会导致计算传递矩阵时间过长，计算效率过低的问题，难以推广应用于工程实际中的实时监测[23-26]。

针对以上问题，本文提出基于经验模态分解和模型缩聚的动力响应重构法。不同于传统的基于带有间歇性准则的经验模态分解法，该响应重构在超单元的基础上完成。根据响应采集点及待测点所在位置划分子结构，建立各子结构的有限元模型，对于工程结构而言，一般选取低阶模态作为子结构的主模态，通过坐标变换，得到各子结构的超单元模型[27-28]。然后，进行二次坐标变换，将各子结构耦合为全模型，提取该模型的模态转换矩阵。结合带有间歇性准则的经验模态分解法，将采集点的响应转换成各个单频响应，通过模态转换矩阵即可得到待重构响应的各单频响应，并通过模态叠加法即可完成超单元基础上的重构，此时所得的重构响应即为待测点的响应。

1 基于经验模态分解和模型缩聚的动力响应重构理论推导

1.1 子结构划分

一般情况下，根据响应采集位置和待测点所处位置划分子结构，将这两个测点所处位置分别作为各自子结构的界面，由此形成子结构。通过有限元法提取子结构的刚度矩阵和质量矩阵，便于后期计算。

对于一般的子结构，其有限元基础上的动力学运动方程可表示为：

根据Craig-Bampton 模态综合法，提取子结构s的模态转换矩阵 Φs，结合布尔矩阵L，通过两次坐标变换，可得到该模型的超单元基础上的动力学运动方程式(2)。各子结构耦合为超单元模型后，各子结构间的界面力gs(t)耦合为0。

式中，各向量主要由两部分构成，各子结构界面自由度集及子结构内部单元各自由度主模态集。皆为超单元模型的质量矩阵、阻尼矩阵及刚度矩阵。模型的超单元位移响应可表示为：

式中：p为超单元模型的位移响应集；pin为第n个子结构的内部单元自由度的主模态位移响应集；pjn为第n个子结构的界面自由度的位移响应集。

1.2 子结构的动力响应重构理论

假设子结构界面自由度m的响应Xm(t)已知，需要重构另一子结构边界自由度r的响应Xr(t)。首先，结合坐标转换将已知响应Xm(t)转换为pm(t)，由于边界自由度为约束模态，则Xm(t)等于pm(t)。基于带有间歇准则的EMD 方法，将pm(t)分解为各单频的时域响应。其中，带有间歇准则的经验模态分解法(EMD)通过以下步骤实现。

首先，确定输入信号pm(t)的极大值及极小值点；然后，由各极值点拟合上包络线及下包络线，计算上下包络线的平均值m1(t)；再由式h1(t)计算第一个分量，h1(t)可表示为：

h1(t)=pm(t)-m1(t)(4)

判断h1(t)是否满足IMF 条件，若否，则将h1(t)视为新的信号重新筛；建立信号h1(t)的包络线，计算上下包络线的平均值m11，得到分量h11(t)，h11(t)可表示为：

重复k次以上筛选过程，直到分量h1k为本征模态函数(IMF)，h1k可表示为：

设数据中筛选得到的第一个IMF 分量f1(t)=h1k，把剩余信号pm(t)-f1(t)重复以上筛选步骤获得第二个IMFf2(t)。再将剩余信号不断循环，依次得到其他剩余的IMF，直到最后剩余信号r(t)为单调函数，则停止筛选，剩余信号r(t)可表示为：

原始信号pm(t)可表示为：

通过EMD 分解得到的IMF 可能含有多个频率成分，为了得到单频的IMF，需要设置带通滤波器。对于大型复杂结构，结构动力响应重构的阶数一般较多，通过傅里叶变换，可得到每阶模态的频率ωi，由此确定带通滤波器的滤波区间[ωiLωiH]，其中ωiL＜ωi＜ωiH。将pm(t)通过滤波器后，再由EMD 分解得到各个单频的IMF，则时域信号可表示为：

式中：di(t) 为第i阶模态响应；si(t)为其余非模态响应的IMF。

界面自由度m的响应pm(t)及自由度r的响应pr(t)可表示为：

超单元的r自由度响应与其对应的有限元响应关系可表示为：

式中，Xr表示超单元的r自由度响应对应的有限元响应。

具体重构过程如图1 所示。

图1 重构流程图Fig.1 Flowchart of the response reconstruction method

2 数值模拟分析

2.1 位移响应重构

通过悬臂梁仿真算例验证本文提出的基于经验模态分解和模型缩聚的动力响应重构法的有效性。如图2 所示，该悬臂梁由铝材制成，铝材弹性模量为69 600 MPa，密度2730 kg/m3，长宽高分别为5 m、0.5 m 和0.05 m。采用ANSYS 建立有限元模型，单元类型为beam3 单元，共有40 个单元，40 个节点，120 个自由度。以自由端为位移响应采集点，跨中节点为位移响应待测点，如图3 所示，以1 单元～20 单元为子结构1，21 单元～40 单元为子结构2。

图3 子结构划分示意图Fig.3 The schematic diagram of the substructures

该重构方法与外界激励无关，将竖向随机荷载均布于悬臂梁上，该随机荷载由6 阶低通巴特沃斯滤波器过滤的白噪声模拟。每阶模态阻尼设为1%，采样频率为1 kHz。为评估该重构方法的有效性，将均方误差(MSE)及相关系数(correlation coefficient)设为评估重构效果的指标，其中MSE可表示为：

式中，z(t)、k和分别表示待测点的重构响应、矩阵z(t)的列数及待测点响应的理论值。

如图4(a)所示为模拟生成的60 s 响应采集点的位移响应，图4(b)所示为该位移响应的频谱图，由该频谱图可知已识别得到的悬臂梁固有频率为1.63 Hz 和10.38 Hz。根据该结构的频率可得表1所示的带通滤波器。如图4(c)和图4(d)所示，结合带有间歇准则的经验模态分解法(EMD)，可将响应采集点的位移响应分解为多个单频时域信号。

图4 响应采集点的位移Fig.4 The displacement measurement data of the response gauge location

表1 带通滤波器的频率范围Table 1 the frequency range of the band-pass filter

由式(13)可得重构所得的响应待测点的位移响应，如图5 所示为响应理论值与重构值的对比图，其中图5(a)表示60 s 时域的对比图，为更清晰地反映重构效果，如图5(b)所示，任取其中10 s的对比结果，由图可知，重构所得的响应与理论值吻合较好。且响应采集点的理论值与重构值的相关系数为0.986(接近1)，均方误差(MSE)为1.28×10-5远小于3%)，以上表明了该重构方法的有效性及精确性。就计算效率而言，结合模型缩聚，悬臂梁有限元模型转换为超单元模型，刚度矩阵和质量矩阵大小由原来的120×120 降低至10×10，由此，计算量减小，计算效率得到提高。基于有限元基础上的经验模态分解法完成该悬臂梁响应重构所需时间为80.42 s，而基于经验模态分解和模型缩聚的动力响应重构方法完成该悬臂梁响应重构所需时间为6.89 s，由此可知，计算效率得以提高。

图5 响应采集点的理论值与重构值对比图Fig.5 The comparison of the theoretical results and the reconstructed results

2.2 噪声对重构精度的影响

结合图2 所示的悬臂梁数值模拟算例分析四个不同噪声等级(噪声等级分别为0%、5%、10%、20%)对悬臂梁响应重构精度的影响。为模拟传感器的测量噪声，生成位移响应后，将噪声单元加入各位移响应。其中，噪声单元指高斯脉冲过程中的加速度的均方根占最大均方根百分比。通过噪声等级来表征噪声单元，如：10%的噪声等级表示该噪声为高斯脉冲过程中的加速度的均方根占最大均方根百分比为10%。基于该重构方法，通过悬臂梁自由端的响应重构悬臂梁跨中节点的响应。如图6 所示，为四个不同噪声等级下待测点的响应重构值与理论值对比图。定义两个指标以衡量噪声等级对重构精度的影响，即均方误差和相关系数。表2、表3 分别表示不同噪声等级下的响应重构均方误差和重构响应的理论值与重构值的相关系数。图7、图8 分别为重构响应的理论值与重构值的均方误差、相关系数随噪声等级的变化曲线图。

表2 不同噪声等级下的响应重构均方误差Table 2 Mean square errors of the four cases

表3 不同噪声等级下的响应重构相关系数Table 3 Correlation coefficients the four cases

图6 四个不同噪声等级下待测点的响应重构结果Fig.6 The reconstructed results at four different noise levels

图7 不同噪声等级下的均方误差Fig.7 Reconstruction performance measured in MSE under different levels

图8 不同噪声等级下的相关系数Fig.8 Reconstruction performance measured in correlation coefficient under different levels

由不同噪声等级下的重构对比图、MSE(皆小于3%)及相关系数(皆接近于1)可知，不同噪声等级下，响应重构的精度皆能得到保证。但随着噪声等级增加，MSE 随之增大，相关系数随之减小，EMD 的边界效应越来越显著，由此可知，噪声等级对重构精确度存在一定的影响。

2.3 主模态数量对重构精度的影响

结合图2 所示的悬臂梁数值模拟算例，分析不同主模态数量对悬臂梁响应重构精度的影响，态)下的重构效果，以研究主模态数量对重构精度的影响。结合该重构方法，通过悬臂梁自由端的响应重构悬臂梁跨中节点的响应。如图9 所示，为4 个工况下待测点的响应重构值与理论值对比图，仍以均方误差(MSE)及相关系数来衡量主模态数量对重构精度的影响。如表4、表5 所示，分别为不同主模态数量下响应重构的均方误差(MSE)和相关系数。由图9、表4 及表5 可知，4 种工况的均方误差(MSE)均较小(小于3%)及相关系数讨论4 种工况(工况1：1 阶主模态；工况2：2 阶主模态；工况3：3 阶主模态；工况4：4 阶主模均接近1。由此可得，主模态的数量对该重构方法的精度影响较小，主要是由于子结构界面的模态转换矩阵由约束模态计算所得，与主模态的相关参数并无太大关系。

表4 各工况下的响应重构均方误差Table 4 Mean square errors of the four cases

表5 各工况下的响应重构相关系数Table 5 Correlation coefficient of the four cases

图9 各工况下待测点的响应重构Fig.9 The comparison of the theoretical and reconstructed results in the different conditions

3 结论

本文的研究工作将有限元基础上的响应重构拓展到超单元，将模态综合法和带有间歇准则的经验模态分解法结合以解决大型复杂工程结构的重构问题。得到的主要结论如下：

(1)该重构方法无需考虑各子结构边界条件，有效地降低了有限元模型相关参数数学矩阵的维度，可以更好地适应大型复杂工程结构的动力响应重构，且较大地提高了响应重构的计算效率。

(2)基于悬臂梁仿真算例，通过该响应重构方法得到的待测点响应重构值与理论值较接近，其相关系数均接近于1 且均方误差(MSE)小于3%，由此验证了该方法的有效性。

(3)噪声及主模态数量影响研究表明，随着噪声等级的提高，该重构方法的精度能得到保证，但仍对重构精度有一定影响；就主模态数量对重构精度影响而言，改变主模态数量对重构精度影响不大，由于子结构界面的模态转换矩阵由约束模态计算得到，与主模态相关参数并无太大关系。