基于模态综合法和模态叠加法的密集模态结构响应重构

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

(1.中南大学土木工程学院，湖南，长沙 410075；2.轨道交通工程结构防灾减灾湖南省重点实验室，湖南，长沙 410075；3.中南大学交通运输工程学院，湖南，长沙 410075；4.中国建筑第五工程局有限公司，湖南，长沙 410004)

土木工程结构的性能在其使用过程中会随着时间的推移而下降，为了监测结构的健康状况和进行寿命预估，有必要准确获得结构关键区域的系统响应[1-2]。结构动力响应重构在过去十几年中受到越来越多的关注。然而，鉴于工程结构的几何复杂性及部件多样性，某些位置(例如结构交界面、狭缝等)皆不易安装传感器，实际中并不总是能在期望的位置进行测量[3-4]。因此，预测关键点的动态响应成为进行结构健康监测的重要组成部分[5- 6]。

目前，响应重构方法主要可分两种：频域方法和时域方法。在频域方法中，RIBEIRO等[7]首先提出了传递率概念，定义为频域中已知响应函数和其他未知响应函数之间的传递率矩阵，通过传递率可以实现结构任意位置的响应重构；DEVRIENDT等[8]提出另一个传递率概念，定义为响应谱和单个参考响应谱之间的比值，在响应重构中仅应用于全局传递性；LAW 等[9- 11]基于频响函数将上述位移传递率推广至加速度传递率，完成了全模型或子结构上的响应重构。以上的频域方法使用频率响应域中的数据，而且通过逆傅里叶变换来将重构的频域响应转换回时域，会引入额外的计算成本和数值截断误差。在时域方法中，HE等[12-13]提出一种基于经验模态分解(EMD)的时域响应重构方法，将采集响应分解为由固有模态函数表示的模态响应，然后，利用这些模态响应计算待测位置的模态响应，最后，采用模态叠加法实现时域动态响应重构。为了确保通过EMD得到本征模态函数(IMF)的只包含一个频率成分，还加入了间歇性准则进行频率筛选。该方法在计算成本方面非常有效，并且适用于基于不同类型的传感器测量的动态响应重构，相比于频域方法有着更高的计算效率。

然而，基于EMD的响应重构方法应用于大型土木结构存在一定困难。一方面，该方法需要获取整个结构的全局刚度矩阵与质量矩阵，而大型土木工程结构通常维数较大，在提取刚度矩阵、质量矩阵时需要大量的计算资源和时间；另一方面，带间歇性准则的EMD方法不能有效地分离密集模态，从而难以实现存在密集模态结构的响应的重构。

针对以上两种限制，本文提出一种基于模态综合法和改进模态叠加法的响应重构方法。该方法结合了模态综合法[14-16]，将结构全模型经过两次坐标变换缩聚为自由度更少的超单元模型，并通过在超单元模型上计算的模态振型进行动态响应重构；当存在密集模态时，在超单元模型的基础上，将结构模态分为密集模态和剩余模态[17-18]，其中剩余模态的重构通过基于EMD的重构方法实现，密集模态的重构响应可由模态振型和剩余模态计算得到，最后通过模态响应叠加获得待测位置响应。

1 理论推导

1.1 子结构模态综合法

对于一般的土木工程结构而言，可以根据结构的连接特性、实际工程需要等因素划分成多个子结构。在有限元基础上，子结构的动力学运动方程可表示为：

式中：上标s表示第s个子结构；Ms、Cs及Ks分别为子结构的有限元模型的质量矩阵、阻尼矩阵及刚度矩阵；Xs(t) 、分别为自由度位移、速度及加速度；fs(t)为子结构所受外力，gs(t)为子结构的界面力。

根据Craig-Bampton 提出的固定界面模态综合法，式(1)中的自由度可分成内部组“i”和边界组“j”，其中，内部组中的自由度不与任何子结构共有，边界组必须包含实际的“边界”，即子结构与子结构之间共有的自由度，也可以将子结构内部任意自由度归到边界组中。将式(1)中的各矩阵的元素按如下形式重新排序：

提取第s个子结构的模态转换矩阵 Φs，该矩阵由两个子集组成，即固定界面选取的主模态集和全部界面坐标的约束模态集：

将式(4)代入式(1)，并令Xs=Φsqs，其中qs表示第s个子结构的模态坐标向量，则式(1)可表示为：

根据式(3)，模态坐标q和原坐标X之间的关系可具体表示为：

即经过第一次坐标变换后，边界组中自由度的模态坐标下的响应与对应的原结构响应相等。

各子结构进行第一次坐标变换之后，整合所有子结构，则整个模型的动力学运动方程可表示为：

其中：

由于各个子结构之间存在共有自由度，因此模态坐标q中的元素不是独立的，需要经过第二次坐标变换消除非独立坐标。结合布尔矩阵L可将各子结构耦合成全模型的超单元模型，当各子结构耦合为一个超单元模型之后，界面力为零，即LTg(t)=0，令q=Lp，其中p表示独立坐标，并代入式(8)，则整个超单元模型的动力学运动方程可表示为：

式(10)中各矩阵可表示为：

第二次坐标变换只消除耦合模型中的非独立坐标或改变各个坐标的排列顺序，并不会改变元素大小，因此根据式(7)，可以得到以下关系：

基于超单元模型的响应重构需要求解模型的超单元模态振型矩阵，如下所示：

1.2 基于改进模态叠加法的子结构重构理论

假设结构第w阶和第v阶模态为密集模态，下标a和r分别表示密集模态和剩余模态，下标m和u分别表示响应采集处自由度和待测响应自由度。将响应采集自由度和待测自由度分别置于边界组中，结合式(11)，超单元基础上的已知响应集与未知响应集如下所示：

式中：Da1(t),Da2(t)···Daa(t)为模态坐标中的密集模态响应；Dr1(t),Dr2(t)···Drr(t)为模态坐标中的剩余模态响应。由式(17)得模态坐标下的密集模态响应为：

式中：上标+表示矩阵的广义逆。在这里，已知响应的数目应不小于密集模态的总数。结合式(16)和式(18)，将各个模态组分进行叠加可得未知响应向量：

定义如下四个矩阵：

式中：下标b表示第b自由度(响应采集量中的某个响应自由度)；下标r1、r2…rr表示剩余的r个模态；为第b个自由度的剩余模态响应。式(20)可表示为如下形式：

若结构的模态中存在密集模态，在响应的傅里叶谱上可能难以辨别其中某一组密集模态具体包含了哪几阶模态，这时候可以选择将该组密集模态可能包含的模态都视为密集模态来计算，并根据式(25)重构。该重构方法适用于应力、应变、位移、加速度等动力响应的重构。

2 数值模拟算例

通过输电塔案例验证所提出的响应重构方法对大型结构的适用性。输电塔的有限元模型由ANSYSAPDL建立，如图1所示，其中单元类型为beam188，共有3123个单元，2764个节点，16 584个自由度，结构底部通过四个固定约束与地面连结。模型建成后，在ANSYS内提取刚度矩阵和质量矩阵并导入MATLAB中，结构响应的模拟和重构都在MATLAB中实现。

图1 输电塔模型Fig.1 The model of power transm ission tower

为了获取结构响应，假设在结构中设置了6个响应采集点，如图1所示，这些采集点位于离地高度较小的位置，也是在实际中便于安装传感器的位置。此外，响应重构点Loc.R 的位置也在图1中标出。将由15阶低通巴特沃斯滤波器过滤的白噪声模拟的随机荷载施加于输电塔的多个节点上，每阶模态阻尼设为1%，采样频率为1 kHz。在超单元模型上根据式(13)计算模态频率矩阵，得出前11阶结构固有频率分别为8.14 Hz、17.89 Hz、36.31 Hz、38.55 Hz、39.38 Hz、45.29 Hz、73.52 Hz、86.41 Hz、87.32 Hz、87.78 Hz 和90.40 Hz，图2是响应采集位置Loc.K 的响应傅里叶谱。通过对六个响应采集点的响应傅里叶谱进行分析，第11阶模态(频率为90.4Hz)及其后的模态并没有明显的幅值，为了简化计算，在本案例中仅考虑对响应有较大贡献的前10阶模态。根据超单元模型计算的固有频率值，将第1阶、2阶、6阶、7阶模态视为单阶剩余模态；第3阶～5阶模态和第8阶～10阶模态内的频率间隔很小，很难使用传统的EMD方法对其进行准确的分离，因此将第3阶～5阶模态和第8阶～10阶模态归入密集模态组。

图2 Loc.K 的响应傅里叶谱Fig.2 Spectrum of the response of Loc.K

按照基于模态综合法和模态叠加法的密集模态结构响应重构方法计算Loc.R 的动力响应，重构的基本流程如图3所示。图4所示为子结构划分示意图，上半部分结构为子结构1，下半部分结构为子结构2，每个子结构缩聚后均保留20个主模态。本算例中超单元模型的边界组自由度包括子结构1和子结构2的界面自由度，以及输电塔底部多个自由度(包含所有响应采集点)和重构点对应的自由度，压缩后的超单元模型共包含148个自由度，在超单元模型上的重构过程如下：

图3 重构流程图Fig.3 Flow chart of the response reconstruction method

图4 子结构划分示意Fig.4 Division of two substructure

1)由式(13)计算超单元模型的振型 Φ～；

2)根据式(21)～式(23)计算Tuam、Subr和Smbr。在该案例中，对应于下标m的响应采集位置包括Loc.K、Loc.L、Loc.M、Loc.N、Loc.O、Loc.P；对应于下标m的响应重构位置仅有Loc.R；对应于下标a的密集模态包括第3阶、4阶、5阶、8阶、9阶、10 阶模态；对应于下标r的剩余模态包括第1阶、2阶、6阶、7阶模态；

3)使用EMD方法提取Loc.K 的剩余模态响应，由此构建矩阵式(24)中的矩阵EMD方法的性能取决于带通滤波器分离不同模式的能力。表1 给出了提取剩余模态响应的带通滤波器的频率范围；

表1 使用EMD方法的每个带通滤波器频率范围Table1 Frequency ranges for each band-pass filter using the EMD method

4)根据式(25)重构Loc.R 的响应。

重构的理论值与重构值对比如图5所示。为了定量评价重构的精度，定义一种相对误差(RE)用以反映重构值偏离理论值的程度，以及平均绝对误差(MAE)用以评价重构值的平均偏移量，如下所示：

式中：n为时间序列长度；p(t)为 重构值；为理论值。图5中重构值和理论值的RE和MAE分别为1.129%和0.103%，重构精度较高。

图5 理论响应与重构响应Fig.5 Theoretical response and reconstructed response

接下来使用传统的EMD 方法进行重构，即按照式(26)的重构方法获取待测位置的响应。选取Loc.K 的响应重构Loc.R的响应，表2给出了每个带通滤波器的频率范围。Loc.R 的重构响应和理论响应对比如图6所示，RE和MAE分别为67.254%和0.815%，重构误差很大，说明在该案例中传统的EMD方法不能准确实现待测位置响应的重构。

表2 使用传统EMD 方法的每个带通滤波器频率范围Table2 Frequency ranges for each band-pass filter using traditional EMDmethod

图6 理论响应与重构响应(使用传统的EMD方法)Fig.6 Theoretical response and reconstructed response(using traditional EMDmethod)

此外，大型结构的参数矩阵(刚度矩阵和质量矩阵)存在着大量的数据冗余，即包含很多零元素，在重构全过程中会耗费大量时间和内存。表3给出了原模型和子结构的数据量对比，可以看出，划分为两个子结构之后，参与运算的数据大约减少一半，有效缓解了计算机内存不足和计算时间冗长的问题。

表3 原模型和超单元模型的运算数据量Table3 Calculation data of original model and super-element model

2.1 主模态数量的影响

根据式(3)可知，在坐标变换过程中人为选择了保留的主模态数量，该数目通常远小于原结构的自由度数量，以此达到模型缩聚的效果。显然，不同的主模态数量会改变最终超单元模型的振型计算值，从而影响重构结果。为了验证主模态数量对所提出的响应重构方法的影响，对比每个子结构选取的主模态数量分别为2、4、6、10和20时的重构结果。待测位置响应的重构值和理论值对比如图7所示，此外还给出了各主模态数量下的重构精度，如表4所示。

表4 不同主模态数量下的重构精度Table4 Reconstruction accuracy under different number of main modes

图7 不同主模态数量下的理论值与重构值Fig.7 Theoretical value and reconstructed value under different number of main modes

2.2 子结构划分方式的影响

表5 不同子结构划分方式下的重构精度Table5 Reconstruction accuracy under different different substructure division

图8 子结构划分方式1Fig.8 Substructure division 1

图10 子结构划分方式3Fig.10 Substructure division 3

图11 不同子结构划分方式下的理论值与重构值Fig.11 Theoreticalvalue and reconstructed value under different substructure division

图9 子结构划分方式2 Fig.9 Substructure division 2

相对于直接提取原模型的参数矩阵而言，通过分别提取各子结构的参数矩阵再进行耦合的方法，在保证重构精度的前提下，可以减少运算的数据量。表6给出了不同子结构划分下的重构运算数据量对比，可以看出，一定范围内子结构数量越多，运算数据总量越少，对计算效率的提升也越明显。但是随着子结构数量的增大，将需要更多人力设置相应的界面耦合。因此在实际中应综合考虑计算效率和人力投入来划分合适的子结构数量。

表6 不同子结构划分的重构运算数据量。Table 6 Calculation data of reconstruction under differentsubstructure division

2.3 噪声等级的影响

现实中响应采集不可避免地存在噪声，为了所提出的响应重构方法在测量噪声环境下的适用性，研究在同一组外部激励下，不同等级的噪声(5%、10%、15%和20%)对重构的影响。测量噪声用高斯脉冲过程来表示，其均方根设置为应变最大均方根的百分比，这里使用噪声等级来表征噪声程度，即添加的噪声是响应最大均方根值等同于噪声等级的高斯过程。重构结果如图12和表7所示，可以看出，重构精度会受到噪声的影响，随着噪声的增大而降低，而且重构响应上的锯齿也趋于明显。

表7 不同噪声等级下的重构精度Table7 Reconstruction accuracy under different noise level

图12 不同噪声等级下的理论值与重构值Fig.12 Theoreticalvalue and reconstructed value under different noise level

EMD方法本身有着一定的降噪性能，而在本文所提出的重构方法中，仅用EMD提取剩余的单频模态响应，对于密集模态响应则通过已知响应、提取的剩余模态响应和振型矩阵外推至响应重构点，这样一来，已知响应中的噪声将被保留至响应重构点，对重构结果产生影响。因此，对于密集模态成分主导的结构响应，该方法难以在高噪声的情况下重构出准确的结果。对于该案例而言，在不大于15%的噪声下进行重构能得到较为准确的结果。

2.4 阻尼的影响

在前文的论述中，该案例采用1%模态阻尼以模拟阻尼效应，而实际工程结构可能存在更高的阻尼，因此接下来讨论不同的阻尼对重构结果的影响。对于同一组外部激励，考虑三种模态阻尼(2%、4%和6%)，不同阻尼下的重构结果如图13和表8所示所示。

图13 不同阻尼下的理论值与重构值Fig.13 Theoreticalvalue and reconstructed value under different damp

表8 不同阻尼下的重构精度Table8 Reconstruction accuracy under different damp

从重构结果可以看出，RE随着阻尼的增大而明显增大，说明在高阻尼下重构精度相对较低。本文所提出的方法将结构响应频率的模态振型和模态响应一一对应实现重构，其基本思想是：结构响应可分解为多个模态响应的叠加，每个模态响应对应着一个固有频率。图14 给出了不同阻尼下的Loc.K 的响应傅里叶谱，可以看出，随着阻尼的增高，各个模态的能量趋于分散，难以与单一频率相对应，因此高阻尼下重构精度差。对于该输电塔结构而言，应保证重构在模态阻尼不大于4%下的情况下进行。

图14 不同阻尼下Loc.K 的响应傅里叶谱Fig.14 Spectrum of the response of Loc.K under different damp

3 结论

本文提出了一种基于模态综合法和模态叠加法的密集模态结构响应重构方法，以输电塔结构为例开展了响应重构数值计算，探讨了主模态数量、子结构划分方式、噪声等级以及阻尼对重构结果的影响，得到主要结论如下：

(1)该方法结合了模态综合法，通过模型缩聚，参与计算的刚度和质量矩阵数据量大大减少，由此改善了计算机内存不足的问题，从而更加适用于大型结构。

(2)鉴于EMD法难以分离密集模态，本文提出将密集模态视作一个整体，并由剩余模态替代的方法，改善了传统EMD方法不能分离频率间隔较小的模态而无法准确重构未知响应这一不足，且适用于于应力、应变、位移、加速度等各种动力响应的重构。

(3)以输电塔为例的数值计算表明，当主模态数量很小时，重构精度较差，随着主模态数量增加，重构精度趋于稳定，可以考虑保留两倍于结构响应贡献较大的模态阶数的主模态数量；子结构划分方式对该重构精度的影响极小，应综合考虑计算效率和人力投入划分适当的子结构数量。

(4)测量噪声等级越高，重构精度越低，对于密集模态主导且噪声等级较高的情况，难以得到理想的重构结果，作者在今后的研究中也将寻求适合的滤波降噪方法应用于该重构方法；结构阻尼对重构有较大影响，这是因为高阻尼下各个模态的能量趋于分散，难以和单一频率相对应，应保证案例的输电塔结构的重构在模态阻尼不大于4%下的情况下进行。