陈 闯， 俞 鹏， 王银辉

(1.浙江大学 宁波理工学院，宁波 315100；2.浙江公路技师学院，杭州 310023)

结构损伤识别一直是结构安全状态评估的热点问题，随着结构健康监测技术的发展和应用，大量、不同类型的结构监测信息为准确、有效的识别结构损伤提供了数据支撑[1]。然而，对于整体结构而言，损伤时刻和损伤位置往往具有随机性，尤其对于早期、微小损伤以及远离监测位置的损伤，监测数据中包含的损伤信息少且信噪比低，这极大的增加了结构损伤识别的难度[2]。

目前，国内外学者基于结构监测数据已提出了一系列的损伤识别方法，主要是包括时域、频域以及时频域的损伤识别方法。基于马氏距离的时域损伤识别方法以其较好的适用性及敏感性在结构损伤识别领域已有较多应用。Mosavi等[3]以马氏距离为损伤位置评判指标提出一种在环境激励下钢梁桥损伤位置的识别方法，工字钢模型试验验证了该方法对损伤具有很好的敏感性。曹军宏等[4]提出一种基于独立分量分析结合马氏距离判断结构损伤的方法，钢框架模型的冲击试验结果验证了该方法的有效性。骆志高等[5]利用不同程度裂纹对应的马氏距离作为标准值，通过对比测试组马氏距离与标准值的差异判断损伤程度。Zhou等[6]通过马氏距离的最大值评判强迫振动下结构的损伤程度，并与频率响应函数为指标的损伤识别方法进行比较，证明了该方法的损伤识别能力。刘纲等[7]利用马氏距离作为损伤特征值用于评价自回归模型的损伤敏感向量，结果表明该方法具有较好的抗噪性。王清华等[8]采用电磁干扰技术结合马氏距离指标对木结构进行了损伤程度评估。相比于直接利用采集的数据构造马氏距离，Szymon等[9]利用结构的加速度响应构造的Hank矩阵计算马氏距离，并与子空间损伤识别方法进行了比较，证明了两种方法的高效性，并能够降低风、温度等周围环境对识别效果的影响。罗钧等[10]利用自回归滑动平均模型的系数向量构造马氏距离对剪切型框架结构进行了损伤定位研究。

综上所述，马氏距离作为结构损伤评判指标的适用性及有效性已被证实。然而，对于结构健康监测中信噪比低、包含的损伤信息较少情况下的损伤识别研究较少。针对以上情况，本文基于马氏距离累积量和经验模态分解提出了一种结构损伤识别定位的“两步法”，数值模拟和模型试验结果验证了该方法的有效性。

1 结构损伤识别“两步法”的基本原理

1.1 马氏距离基本原理

1936年印度统计学家Mahalanobis提出了马氏距离的基本原理[11]，其主要目的是比较两样本(集合)的差异性。马氏距离的特点是能够考虑样本中元素的大小和特性之间相关性，同时不区分数据类型，其计算公式为

(1)

1.2 经验模态分解(EMD)基本原理

(2)

式中，ci为第i阶IMF，rn(t)为残差项。

经验模态分解能够将非平稳信号按照频率高低排列成各阶本征模态函数，结构损伤导致的刚度降低将引起结构响应中某些频带能量的改变，从而造成频段能量的转移，这种频带能量的转移对结构的损伤具有较好敏感性[13-14]。

1.3 损伤识别向量的构造

结构发生损伤通常导致结构动力参数的改变，然而，较小的参数变化很难有效识别。根据结构损伤的“单一向”特性(即一旦损伤发生，产生的结构性裂缝只能一直存在且持续增加)，可通过损伤信息的“叠加累积”形成一种“放大效应”，使微小损伤有效识别，且信噪比低的信号有效提取损伤信息。本文通过将监测数据样本的马氏距离进行累积形成损伤信息的“放大效应”，从而实现损伤信息的累积。

假设监测数据样本长度为l，等分成k组，每组数据长度为j，即l=j×k。定义第i组马氏距离累积量(Mahalanobis distance cumulates, MDC)为

(3)

式中，xn为待测样本中的各类数据，μ为参考样本的均值估计，S为参考样本的协方差估计；MDCi为第i组马氏距离累积量。损伤识别向量Rk为各组MDC值组成的向量，表示为

(4)

1.4 结构损伤识别的“两步法”

本文提出的结构损伤识别“两步法”主要是针对由于损伤较小或传感器测点距离损伤较远而导致的监测信息中包含的损伤信息较少且信噪比低的情况，首先需要对原始的监测数据进一步处理，即对原始监测数据进行经验模态分解，进而在更精细化的频段内进行损伤识别。具体的步骤为

步骤1 利用损伤前后直接监测的物理参量(加速度、应变等)作为参考样本和待测样本计算MDC值构造损伤识别向量，并以参考样本损伤识别向量MDC值的均值估计作为阈值，即

(5)

式中，mean(Rk)为Rk向量的均值估计。当采集到的物理参量信噪比较高，信号中包含了较多的损伤信息时，往往MDC值的变化已经很明显，则仅需第一步即可识别出裂缝引起的刚度降低。

步骤2 对原始结构监测数据进行EMD分解，获得各阶IMF，以损伤前后的各阶IMF作为参考样本和待测样本构造损伤识别向量，并建立各阶IMF样本的MDC值概率密度函数，以参考样本MDC值概率密度函数95%置信区间的上限值作为损伤识别的阈值，并通过待测样本与阈值之间的对比关系进一步识别结构损伤。

本文提出的结构损伤识别两步法，在数值模拟过程中，微小的刚度变化经过第一步方法即能够识别出损伤；然而，在模型试验和实际工程中由于噪声、温度等影响，为损伤的精确识别增加了难度，距离损伤位置较远处的传感器或者损伤较小的情况下的MDC值变化不明显，因此，在直接利用监测数据无法明确识别出损伤的情况下，则进入第二步基于原始监测数据EMD分解后，利用各阶IMF进行更精细化的损伤识别。

2 数值模拟算例

2.1 有限元模型

采用有限元分析软件ANSYS14.0建立跨径为10 m的简支梁模型，主梁采用BEAM4单元模拟，截面高和宽分别为0.6 m和0.3 m，惯性矩Izz和惯性矩Iyy分别为0.001 35 m4和0.005 4 m4，弹性模量为2.06 GPa，密度为7 850 kg/m3，泊松比为0.3。主梁共划分成25个单元和26个节点，每个单元长度为0.4 m。在节点1和节点26处设置支座，简支梁单元和节点编号如图1所示。

图1 简支梁有限元模型(cm)

在主梁13号节点施加F=100×Sin(2t)的正弦力激励，其振幅为100 kN，荷载激励时间步长为0.1 s，数据采集时长为500 s。目前，通常采用降低单元刚度的形式来模拟单元的损伤，本文通过降低9号单元刚度10%和20%两种工况模拟结构损伤。图2为3号节点和10号节点损伤前后的加速度时程曲线，可以看出，损伤前后加速度变化不大，单从加速度时程曲线上很难识别出结构损伤。

(a) 3号节点加速度时程曲线

(b) 10号节点加速度时程曲线

图2 损伤前后简支梁节点加速度时程曲线

Fig.2 Acceleration time-history curves of nodes of the simple support beam before and after damage

2.2 损伤位置识别

按照上文1.2节的方法对加速度数据构造马氏距离累积量MDC值，首先，分别以健康状态下的加速度数据作为参考样本，损伤后的加速度数据作为待测样本，将数据按照j=100进行分段，计算参考样本和待测样本的MDC值变化如图3所示。

(a) 损伤前后3号节点MDC值变化图

(b) 损伤前后10号节点MDC值变化图

Fig.3 MDC value variation node 3 and node 10 before and after damage

从图3中可以看出，当结构出现损伤后，3号和10号节点的MDC值明显的增大，通过比较健康状态与损伤状态的MDC值即可识别结构是否发生了损伤。随着单元损伤程度的加深(9号单元刚度降低10%到降低20%)，MDC值的整体增加愈加明显。所以，当采集数据包含了较多损伤信息或者其损伤较为明显时，其MDC值的变化较为明显，可以直接通过对比MDC值的变化即可识别出结构的损伤。图4为各个节点在损伤前后MDC值的均值，通过比较均值的变化可判断结构损伤的位置。

图4(a)表示损伤前后各个节点处MDC的均值，可以看出，结构损伤导致了MDC均值明显增大，且随着损伤的加深，MDC均值变化愈加明显。图4(b)为损伤前后各个节点相较于健康状态MDC值均值的变化情况，与图4(a)变化一致，发生损伤单元即9号节点处MDC均值的变化量较其他单元更为明显，形成了明显的峰值。因此，可通过各节点处MDC均值的变化情况来判断损伤的位置。

2.3 抗噪性分析

噪声是工程试验和现场监测过程中无法避免的，主要是由于环境激励、测试仪器与传感器本身等因素所导致的，结构在动力测试的过程中，一般会将环境等产生的噪声作为一种激励结构振动的输入，则其输出信号中包含了原测试信号与噪声。通常采用信号与噪声的相对比例，即信噪比描述监测数据中噪声对监测数据的影响。信噪比越小，说明输出信号中的噪音越大，反之亦然。现对10号节点健康状态、单元损伤10%与损伤20%情况下结构的加速度信号中加入5%高斯白噪声，并对其损伤识别情况进行分析，以验证该方法的抗噪性。图5分别为10号节点在加入了5%的高斯白噪声后的加速度时程损伤前后的对比图。可以看出，通过对比损伤指标向量MDC值均值变化依然可以较好的识别出结构是否发生损伤。

(a) 各节点MDC均值图

(b) 各节点MDC均值相对变化图

Fig.4 MDC average values and relative variation of MDC average values of each node of the simply supported beam

3 模型试验方法验证

结构在实际监测环境中，往往会受到如风、温度等外界环境的影响；加之，结构损伤具有随机性，传感器与损伤的相对位置往往决定了是否能够收集到足够的损伤信息；此外，初始损伤信噪比低，也极大的增加了损伤识别的难度。本节对工字钢进行了损伤识别的模型试验，从梁底向上切割工字钢梁人为制造损伤，距离损伤远近不同的位置分别布置了应变传感器，以验证监测位置距离损伤较远、信噪比较低情况下的损伤识别效果。

3.1 工字钢模型及损伤试验

选取一根14#工字钢进行模型试验研究，其横截面尺寸为140 mm×80 mm×5.5 mm(高度×宽度×腹板厚度)，质量密度为16.89 kg/m，标准跨径为2.3 m，计算跨径为2.0 m。

(a) 加入5%白噪声后损伤前后加速度对比图

(b) 10号节点MDC值变化图

图5 简支梁10号节点抗噪性分析

Fig.5 Anti-noise analysis of the node 10 of the beam

工字钢底部等间距布置5个断面的应变测点，每个断面2个应变测点，测点编号为y1～y10；布置环境温度传感器T1及工字钢表面温度传感器T2和T3，传感器布置及现场试验采集设备如图6和图7所示。利用激振器在工字钢接近跨中位置进行强迫振动激励，采样频率50 Hz，采集健康状态下24 h的动态应变数据。在工字钢底部向上切割5 mm和10 mm损伤，分别定义为损伤工况a和损伤工况b，并继续进行强迫振动，保持50 Hz的采样频率，分别各采集24 h应变数据。

图6 应变片与温度传感器布置图

图7 工字钢试验加载及数据采集

3.2 趋势项分析

通过24 h的动态应变数据采集，获得在健康状态下工字钢底部的应变时程曲线如图8所示，可以看出，应变时程具有明显的趋势项，采用滑动平均值法[15]获得整体趋势项及剔除趋势项后应变。图9为提取的应变趋势项，与图10测试得到温度变化对比可以看出，其应变趋势项与温度变化具有较高的相关性，表1列出了应变趋势项与温度的相关性分析结果，可以看出，各测点的应变趋势项与工字钢表面温度的相关系数均较高。因此，可以认为应变趋势项为温度趋势项，在数据分析时将温度趋势剔除，只分析在外荷载激励下的受迫振动数据。

图8 应变时程曲线

3.3 损伤识别分析

以健康状态下剔除趋势项后的应变数据为参考样本，损伤工况a和损伤工况b下剔除趋势项后的应变数据计为待测样本。参考样本和待测样本均按照j=10 000进行分段，分别计算损伤识别向量的MDC值，各个位置损伤前后MDC值的变化情况如图11所示。

图9 应变趋势项时程曲线

图10 大气温度及工字钢表面温度时程曲线

Fig.10 Time-history curves of the atmospheric temperature and the surface temperature of thei-beam

表1 应变趋势项与温度相关系数

Tab.1 Correlation coefficients between the strain trend items and the temperature

相关系数应变测点y1y3y5y7y9温度测点T1T2T30.710.930.850.700.940.850.690.940.850.710.920.830.710.930.85

通过图11可以看出，测点y3和y5的MDC值相较于其他传感器变化更为明显，MDC值均值也随着损伤的发生而逐渐增大，并且随着损伤的加重，其MDC均值的增幅更大，可判断损伤位置在测点y3和y5之间，这与实际损伤的损伤位置相吻合。因此，可直接通过损伤识别过程中的第一步，即可识别出损伤及损伤位置。

3.4 经验模态分解的应用

为了进一步验证本文提出方法在距离损伤位置较远处和直接利用MDC值很难判断是否损伤发生情况下的识别效果，选取测点y1和y3两位置处的应变结果进行分析。通过G.Rilling的EMD工具箱对各测点的应变数据进行分解得到各阶的本征模态函数，测点y3健康状态下的EMD分解后第3阶本征模态函数如图12所示。

(a) 测点y1损伤前后MDC值变化

(b) 测点y3损伤前后MDC值变化

(c) 测点y5损伤前后MDC值变化

(d) 测点y7损伤前后MDC值变化

(e) 测点y9损伤前后MDC值变化

图11 工字钢部分应变测点损伤前后MDC值变化

Fig.11 MDC values of partial testing strain points of thei-beam before and after damage

图13为y3应变剔除趋势项后信号构造了损伤识别向量MDC值的概率密度函数，分别采用Burr分布、Normal分布和Logistic分布拟合概率密度函数得到的结果如表2所示。采用极大似然法对3种概率密度函数拟合结果对比发现，Burr分布拟合的极大似然值最大，因此，可采用Burr分布对MDC值进行分布拟合，其95%置信区间为[963.3,1 045.5]，并将其置信区间上限作为损伤识别向量的阈值。

表2 概率密度函数拟合结果对比

Tab.2 Comparison results of the probability density function fitting

拟合方式Burr分布Normal分布Logistic分布极大似然值-4 375.67-4 478.33-4 429.21

图12 健康状态下测点y3经验模态分解的IMF3

图13 健康状态下测点y3 MDC值概率密度函数拟合图

Fig.13 Probability density function fitting of MDC values of testing pointy3in health condition

对测点y1与y3的应变数据利用EMD分解得到各阶IMF，再分别对各阶IMF构造损伤识别向量以判断损伤情况，因篇幅有限，只列出前四阶IMF构造的MDC值概率密度函数拟合结果如图14和图15所示。

(a) IMF1损伤前后MDC值概率密度函数拟合

(b) IMF2损伤前后MDC值概率密度函数拟合

(c) IMF3损伤前后MDC值概率密度函数拟合

(d) IMF4损伤前后MDC值概率密度函数拟合

图14 测点y1前四阶IMF损伤前后MDC值概率密度函数拟合

Fig.14 Probability density function fitting of MDC values of first four order IMFs of the testing pointy1before and after damage

图14仅列出了测点y1经EMD分解后前四阶IMF损伤前后的MDC值概率密度函数拟合结果。未发生损伤状态下，各阶IMF损伤识别向量MDC值的概率密度函数应较为一致；一旦结构发生损伤，某些频率IMF构造的损伤识别向量MDC值的概率密度函数面积将显著超过阈值，在统计学角度上表明超过阈值的MDC值的概率显著增加，超过阈值面积越大，损伤程度也越大。本例中测点y1利用分解的IMF3、IMF4、IMF6和IMF7信号构造损伤识别向量后，都较好的识别出结构是否发生损伤，其他阶IMF的MDC值变化不明显，且IMF3和IMF4相对于其他阶本征模态函数而言，具有更好的识别效果，其MDC值的概率密度函数随着损伤程度的增加变化更为明显。这是由于结构损伤导致了结构整体刚度的改变，在输入能量基本不变的情况下，各阶IMF频带能量发生了重分布，各阶频率对于损伤的敏感程度也各不相同。图15是y3应变经EMD分解后损伤前后前四阶IMF的MDC值概率密度函数拟合结果，测点y3的损伤识别情况与测点y1的较为相似，其IMF3、IMF4和IMF8信号构造损伤识别向量后，可以较好的识别出结构的损伤情况，而IMF3和IMF4相较于其他几阶IMF，更好的识别出结构是否发生损伤，且其对损伤程度敏感度明显大于其他阶本征模态函数。

(a) IMF1损伤前后MDC值概率密度函数拟合

(b) IMF2损伤前后MDC值概率密度函数拟合

(c) IMF3损伤前后MDC值概率密度函数拟合

(d) IMF4损伤前后MDC值概率密度函数拟合

图15 测点y3前四阶IMF损伤前后MDC值概率密度函数拟合

Fig.15 Probability density function fitting of MDC values of first four order IMFs of the testing pointy3before and after damage

4 结 论

本文提出了一种基于马氏距离累积量和经验模态分解的结构损伤识别“两步法”，首先，直接利用原始监测数据的马氏距离累积量构造损伤识别向量，并通过马氏距离累积量的均值作为阈值进行初步判断；在损伤较小、信噪比低的情况下，通过EMD方法将原始监测数据按频率高低分解成各阶IMF再构造损伤识别向量，并对其进行概率密度函数拟合，以概率密度函数95%置信区间上限值作为阈值进一步识别损伤。通过有限元模型和工字钢模型试验对该方法进行了验证，得出主要结论如下：

(1) 利用原始监测数据的马氏距离损伤累积量均值作为阈值，在损伤较为明显，监测数据中包含了较多损伤信息时可直接识别出损伤及损伤位置；随着损伤程度的增加，损伤识别向量MDC值的变化愈加明显。

(2) 温度趋势项是监测数据中最主要的影响因素，在利用该方法进行实际结构损伤识别时，必须先剔除温度趋势项的影响。

(3) 结构损伤造成了结构动力响应中各阶频率能量的改变，损伤识别向量MDC值的概率密度函数能够有效反映出结构刚度降低引起的频率能量改变，但损伤前后频率能量的重分布模式有待于进一步研究。