基于时间序列模型的结构损伤识别方法

张玉建, 罗永峰, 郭小农, 刘 俊, 朱钊辰

(同济大学 土木工程学院, 上海 200092)

既有结构的定期检测和维护,是保障结构安全、正常服役的前提[1-3],而对结构损伤的准确识别则是进行检测和维护的基础.现有的损伤识别方法包括基于结构振动特性的识别方法和基于时频响应的识别方法[4-5].前者一般通过结构的模态参数进行损伤识别[6-13],然而,由于实际结构的高阶振型难以获取、高阶振型频率难以测量,此类方法的应用受到制约.后者是通过对结构的振动响应进行时域或频域处理,然后提取损伤特征量,进而评估结构的损伤状况.基于时频响应的识别方法因其具有数据样本易于获取、识别过程简单快速等特点,成为结构损伤识别领域的研究热点.

基于时频响应的识别方法需要对结构动态响应数据进行时域或频域处理.作为时域数据处理的重要方式,时间序列分析已被众多学者应用于结构损伤识别领域.刁延松等[14]利用计量经济学中的协整理论对时间序列的自回归模型(auto-regressive model, AR)系数进行处理,提取损伤指标并进行结构损伤识别.Noh等[15]采用NCREE标准模型建立AR模型,根据AR模型前三阶系数定义损伤因子并进行结构的损伤定位.上述方法仅利用模型的有限阶系数构造损伤特征参数,在数据处理过程中存在损伤信息遗漏问题,导致对部分损伤工况产生误判.Fugate等[16]基于AR模型残差的损伤识别方法,提出采用质量控制图形法监测结构状态.Sakellariou等[17]提出基于随机输出误差振动的地震激励下结构损伤检测与评估方法.这类方法仅考虑了单位置损伤识别,难以保证多位置损伤情况下的识别效果.Gul等[18-19]提出基于自回归滑动平均模型(auto-regressive moving average model, ARMA)的统计模式识别方法,以提高时间序列方法在损伤检测过程中的鲁棒性,并对不同结构的时间序列分析进行实测验证,但该方法仅验证了特定损伤程度下的损伤识别效果,难以保证不同损伤程度下的损伤识别效果.

分析现有研究文献可知,以时间序列分析作为时域数据处理方式的损伤识别方法,在应用过程中仍存在一些问题.产生这些问题的根本原因在于方法所建立的时间序列模型仅为统计意义上的数据序列规律总结,与结构本身的物理特性很难建立直接联系,进而导致数据处理过程中存在损伤信息遗漏,易产生误判现象,对多位置、多程度的复杂损伤工况难以进行有效处理,也难以实现结构损伤的定量识别.

针对上述问题,本文提出一种基于具有外部输入的自回归模型(auto-regressive model with eXogenous input, ARX)的损伤识别方法,通过建立考虑结构动力特性的ARX时间序列模型,对实测数据进行拟合,根据数据拟合度构造表征结构损伤的参数,即损伤因子,进而评估结构的损伤状态.以简支梁模型与桁架结构模型为例,通过一系列的损伤工况,分别对多位置和多程度的结构损伤进行识别,并分析激振位置和测量噪声对损伤识别结果的影响,证明了本文所提方法的有效性与适用性.

1 ARX模型基本原理

1.1 ARX模型一般表达式

时间序列是指按照时间顺序或其他物理量顺序排列的一组随机数据[20],而时间序列模型则是揭示数据内部排列规律的系统函数.时间序列分析需要承认随机数据之间的有序性和相关性,并通过数据内部的相互关系来辨识系统的变化规律.常用的时间序列模型包括具有外部输入的自回归滑动平均模型(auto-regressive moving average model with eXogenous input, ARMAX)、ARX模型、AR模型等.ARMAX模型的基本形式[21]如式(1)所示:

y(t)+a1y(t-Δt)+…+anay(t-naΔt)=

b1u(t)+b2u(t-Δt)+…+bnbu(t-nbΔt)+

e(t)+c1e(t-Δt)+…+cnce(t-ncΔt)

(1)

式中:y(t)、u(t)和e(t)分别为模型的输出、输入和白噪声干扰项,a1,…,ana、b1,…,bnb、c1,…,cnc均为待定系数;Δt表示时间序列数据点的时间间隔.

式(1)可采用更简洁的方式进行表达,即:

A(q)y(t)=B(q)u(t)+C(q)e(t)

(2)

式中:A(q)、B(q)和C(q)均为包含模型系数的多项式:

A(q)=1+a1q-1+…+anaq-na

(3)

B(q)=b1+b2q-1+…+bnbq-nb

(4)

C(q)=1+c1q-1+c2q-2+…+cncq-nc

(5)

式(3)~(5)中:q为后移算子,如时间t的变量xij(t)乘以k阶后移算子后可得xij(t-k),即:

xij(t)qk=xij(t-k)

(6)

根据各多项式阶数取值的不同,ARMAX模型可退化为常见的一些时间序列模型.例如,当nb和nc均为0时,式(1)将转化为AR模型;当na和nb均为0时,式(1)转化为滑动平均模型(moving average model, MA);仅当nc为0时,式(1)转化为ARX模型.ARX模型的基本形式为

y(t)+a1y(t-Δt)+…+anay(t-naΔt)=

b1u(t)+…+bnbu(t-nbΔt)+e(t)

(7)

当输入u(t)与其系数bt为向量时,该模型即为多输入单输出模型.以双输入单输出模型为例,此时bt为两元素的行向量,u(t)为两元素的列向量,模型的基本形式如式(8)所示:

y(t)+a1y(t-Δt)+…+anay(t-naΔt)=

b11u1(t)+…+b1nb1u1(t-nb1Δt)+

b21u2(t)+…+b2nb2u2(t-nb2Δt)+e(t)

(8)

式中:bij为第i个输入变量的第j阶待定参数.

1.2 考虑结构动力特性的ARX模型

多自由度体系粘滞阻尼系统的振动微分方程为

(9)

根据式(9),第i自由度的运动方程可表示为

(10)

式中:N为体系自由度数,i为自由度编号.

式(10)对x求二阶导数,可得:

(11)

式中:括号内数值表示对变量求导的阶数,后同.

当结构运动形式为自由振动时,式(11)中右侧荷载项为0,即:

(12)

(13)

(14)

式中:Δt为数据的时间间隔.

将式(13)和(14)代入式(12),整理可得:

(15)

考虑到质量矩阵与刚度矩阵的稀疏性,只有当两个自由度相互关联时,对应的矩阵元素才不为0.以N= 3为例,此时与自由度i的相邻自由度数为2,则式(15)可表示为

式中:

(17)

2 基于ARX模型的损伤识别方法

2.1 损伤识别机理

由1.2节分析可知,按式(16)建立的对应于自由度i的ARX模型包含结构的质量、刚度和阻尼等参数.在结构发生损伤后,伴随着结构动力特性的变化,由无损结构所建立的各节点ARX模型将无法反映损伤结构的节点加速度动态响应规律.因此,在建立考虑结构动力特性的ARX模型后,可利用此模型进行结构的损伤状态评估.

为了更好地理解ARX模型在损伤识别领域的应用机理,本节采用一个包含三个自由度的弹簧质点子系统进行说明,如图1所示.该系统各质点均沿水平方向自由振动,其中k、c分别表示质点间弹簧的刚度系数和阻尼系数,m1、m2、m3表示各质点的质量.

图1 弹簧质点系统

(18)

2.2 损伤特征因子

(19)

式中:p为计算阶次;Z表达为式(20):

(20)

式中:n为样本数量.

根据定义,可以得到拟合度ξ所在区间为(0,1].拟合度越大,代表待识别结构与未损伤结构越接近.无损伤情况下,yi将与Yi完全相等;而当结构出现损伤时,式(16)中模型各系数也随之变化,由无损伤情况下计算得到的ARX模型预测得到的加速度时程数据,将与实测的加速度时程数据产生偏差.因此,可根据拟合程度定义损伤因子(damage factor, DF)数值大小为1-ξ.

2.3 损伤识别流程

由2.1节分析可知,可将模型预测所得某一节点的加速度时程数据与实测加速度时程数据进行对比,根据两者拟合度的大小,评估该节点附近单元的损伤状态.由此,本文提出基于ARX时间序列模型的损伤识别方法,具体分析步骤如下:

(1) 针对损伤识别对象的结构形式,根据式(7)、(15),确定所采用ARX模型的各多项式阶数;

(4) 针对待测结构制定测量方案,布置加速度传感器,施加外界激励(如锤击),获取结构各节点实测加速度时程数据{Y1,Y2,…,Yi,…};

(5) 将Yi-1与Yi+1作为输入,代入对应于节点i的ARX模型,得到未损伤结构中节点i的加速度响应yi;

(6) 代入yi与Yi至式(19)、(20),计算节点i的损伤因子DF;

(7) 改变i,按步骤(5)~步骤(6)重复计算,最终得到各节点的损伤因子DF;

(8) 根据DF数值对结构损伤状态进行评估.

3 数值算例

为验证所提方法的有效性与适用性,本文分别采用一简支梁结构和桁架结构作为数值算例进行模拟验算.

3.1 算例一

3.1.1数值模型

简支梁跨度L=10 m,采用箱型截面,截面尺寸为□300 mm×200 mm×10 mm×10 mm,弹性模量E=2.1×105MPa,泊松比ν=0.3,材料密度ρ=7 850 kg·m-3,结构模型见图2.采用ANSYS v17.0建立有限元模型,单元类型为Beam188单元,沿梁长度方向划分24个单元.阻尼类型采用Rayleigh阻尼,阻尼比取0.02.以瞬时冲击荷载作为激励,获取模型各节点竖向加速度响应作为仿真数据,加速度信号的采样频率为100 Hz.本文共选取2个位置(节点号13、19)作为激励点,并分别定义为①号激励点与②号激励点,冲击荷载大小为100 N,荷载作用时间0.000 1 s.图2所示模型上方数字为节点编号,下方数字为单元编号.

图2 简支梁模型

由于既有结构构件的损伤类型较多,本文仅讨论钢结构中常见的材料腐蚀、截面畸变、壁板变形等缺陷所引起的杆件刚度削减,不考虑服役过程中构件的质量损失与阻尼变化.文中构件的损伤程度即为构件刚度削减的相对值.考虑到结构的对称性,假定简支梁损伤出现在4处不同的位置,损伤位置为图2所示阴影区域,采用弹性模量折减的方式进行损伤模拟,弹性模量折减程度即为损伤程度.

需要说明的是,在利用相邻节点的加速度时程数据建立ARX模型预测某节点的加速度时,由于支座节点固定,不发生振动,因此,本文所提方法无法识别支座附近单元的损伤情况.

施加激励后,获取任意两个节点的加速度时程反应,并截取前2 s数据共200个数据点进行记录,如图3所示.其中,示例1对应节点远离激振位置,示例2对应节点靠近激振位置.由图3可知,本模型的自由振动在前2 s已基本衰减为0.此外,考虑到实际检测过程中可能出现人工操作引起的误差,需略去施加激励后一定时间段内的数据样本.因此,本文截取1~2 s时间区段的100个数据建立ARX模型,并进行损伤识别.在按式(19)计算损伤因子时,设定计算阶次p=2.

a 示例1

b 示例2

3.1.2损伤工况

定义损伤工况如表1所示,其中工况1~工况4为单位置损伤,工况5~工况8为多位置损伤.

表1 数值模拟损伤工况(算例一)

3.1.3测量噪声

实际工程中,实测输出响应信号通常会包含一定水平的测量误差.已有研究表明,测量误差具有一定同济规律,因此可以利用概率学理论和数理统计学方法对测量误差进行统计分析和数值模拟,从而弥补采用计算机数值模拟无法得到真实测量数据的不足.

本文采用文献[22]的方法,对测点的加速度响应施加正态分布的白噪声,作为测量噪声以引入测量误差,取噪声强度为n%,即测点加速度响应均方值R的n%,施加噪声后的响应为

at=at,0+rSN·R·n/100

(21)

式中:at,0和at分别为施加噪声前、后的加速度响应;rSN为符合标准正态分布的随机数.

本文数值算例中,噪声强度分别取0(无噪声)、1%、2%、4%.

3.1.4结果与分析

(1) 损伤位置识别

应用本文方法对工况1~工况8进行损伤识别.图4为选择①号激振点进行激励,在不考虑测量噪声的情况下4种单位置损伤工况的损伤识别结果,图中阴影部分表示损伤位置.图5为选择①号激振点进行激励,在不考虑测量噪声的情况下4种多位置损伤工况的损伤识别结果.以工况6为例,将识别结果与各位置单独发生损伤时的识别结果进行对比,如图6所示.

a 损伤工况1

b 损伤工况2

c 损伤工况3

d 损伤工况4

a 损伤工况5

b 损伤工况6

c 损伤工况7

d 损伤工况8

图6 识别结果对比(工况6,①号激振点,无噪声)

图4表明,对于单位置损伤工况,损伤杆件附近区域的DF值较大,并在与损伤杆件相连的节点达到峰值,由此可以对杆件损伤位置做出较为准确的判断.图5表明,多位置损伤工况下各损伤位置处的DF显著高于其他位置,因此可以较为准确地判断损伤位置.图6表明,多位置损伤识别结果近似为各单位置损伤识别结果的线性叠加,各位置损伤对识别结果的相互影响可忽略不计.

(2) 损伤程度识别

为考察不同杆件损伤程度的识别效果,本节以工况2和工况6为例,在原有20%损伤程度的基础上新增5组不同程度的损伤情况,损伤识别结果如图7所示,损伤位置处相邻节点的DF值随损伤程度变化曲线如图8所示.

a 损伤工况2

b 损伤工况6

图7表明,不同损伤程度下各节点DF分布形式基本相同,即使在较小损伤程度的情况下,仍然能够较为准确地确定损伤位置.图8表明,损伤位置相邻节点的DF值与损伤程度高度线性相关.因此,在进行结构损伤程度识别时,可将某一特定的损伤程度设定为基准损伤工况,并通过有限元数值模拟手段计算该损伤程度下节点的DF数值.根据DF与损伤程度的线性关系,通过对比实测结构与基准损伤工况下的DF数值,即可对损伤程度进行较为准确的判断.

a 损伤工况2

b 损伤工况6

(3) 激振位置对损伤识别结果的影响

为考察激励位置对杆件损伤识别效果的影响,本节以工况2和工况6为例,分别以激励点①和激励点②作为激励位置,应用本文方法对该简支梁结构进行损伤识别,得到不同激励位置的结构损伤识别结果如图9所示.

a 损伤工况2

b 损伤工况6

图9表明,不同激励位置的各节点DF值具有微小差异,但对损伤识别结果影响较小,可忽略不计.实际工程应用中,可根据不同激励位置所求得的DF数值进行综合判断,从而获得较为准确的识别结果.

(4) 测量噪声对损伤识别结果的影响

为考察测量噪声水平对杆件损伤识别效果的影响,本节以工况2和工况6为例,引入三组不同程度的测量噪声并将其识别结果与无噪声情况下的识别结果进行对比,对比结果如图10所示.

a 损伤工况2

b 损伤工况6

图10 不同噪声水平下多位置损伤识别结果(工况2、工况6,①号激振点)

Fig.10 Multi-position damage identification results under different noise levels(condition 2 & 6, excitation point 1)

图10表明,测量噪声对单位置及多位置损伤识别结果均具有一定程度的影响.随着测量噪声强度等级的不断提高,未损伤区域的节点DF变化显著,因此,可能存在高噪声环境下损伤位置的误判情况.不同噪声等级下损伤位置节点DF基本保持不变,因此,在准确识别损伤位置的情况下,测量噪声对损伤程度识别结果的影响较小,可忽略不计.

3.2 算例二

3.2.1数值模型

由1.1及1.2节理论推导部分可以看出,本文所提方法基于双输入单输出ARX模型,适用于梁式结构受力体系的损伤识别.桁架结构依据其实际受力形式可等效为普通梁结构,如图11所示,因此可用本文所提方法对其进行损伤识别.

图11 桁架结构模型

桁架结构为上弦支撑的平面静定钢桁架结构[12-13],跨度L=18 m,高度h=2 m.上下弦杆长度均为2 m,斜腹杆长度为2.24 m,采用圆钢管截面,截面面积为5×10-4m2,弹性模量E=2.1×1011Pa,泊松比ν=0.3,材料密度ρ=7 850 kg·m-3,结构模型如图11所示.采用ANSYS v17.0建立有限元模型,单元类型为Link1单元.阻尼类型采用Rayleigh阻尼,阻尼比取0.02.以瞬时冲击荷载作为激励,获取模型各节点竖向加速度响应作为仿真数据,加速度信号的采样频率为100 Hz.本文选取节点13作为激励点,冲击荷载大小为100 N,荷载作用时间0.000 1 s.图11所示模型节点位置处数字为节点编号,单元旁数字为单元编号.

考虑到结构的对称性,假定桁架结构损伤出现在4处不同的位置,损伤位置为图11所示阴影区域(或称加粗区域),采用弹性模量折减的方式进行损伤模拟,弹性模量折减程度即为损伤程度.在识别过程中,上下弦杆的损伤可通过等效梁结构相应单元中部节点的DF值进行判断,如单元11可通过5号节点进行判断;斜腹杆的损伤可通过等效梁结构相应单元两侧节点的DF值进行判断,如单元24可通过7、8号节点进行判断.

3.2.2损伤工况

定义损伤工况如表2所示,其中工况1、2为单位置损伤,工况3、4为多位置损伤.

表2 数值模拟损伤工况(算例二)

3.2.3结果与分析

应用本文方法对工况1~工况4进行损伤识别,识别结果如图12所示.为考察不同杆件损伤程度的识别效果,本节以工况2为例,在原有20%损伤程度的基础上新增4组不同程度的损伤情况,损伤识别结果如图13所示,损伤位置处相邻节点的DF值随损伤程度变化曲线如图14所示.

a 损伤工况1

b 损伤工况2

c 损伤工况3

d 损伤工况4

图13 不同损伤程度下的单位置损伤识别结果

图14 损伤因子-损伤程度曲线

图12表明,各损伤工况下均可较为准确地识别损伤位置.值得注意的是,当损伤杆件为上下弦杆时,损伤区域中部节点的DF数值较大.同时可以看出,对于上下弦杆,若以原结构的单元相邻节点的DF值进行判断,识别效果不尽理想,原因在于在等效梁结构中,原结构弦杆的损伤相当于两个相邻单元同时发生损伤,此时两单元公共节点的DF值会高于端部的两节点DF值.因此,以中部节点的DF数值作为弦杆的损伤特征指标会获得更好的识别效果.

图13表明,不同损伤程度下各节点DF分布形式基本相同,即使在较小损伤程度的情况下,仍然能够大致确定损伤位置,但其准确性会有所降低.图14表明,与简支梁结构相同,桁架结构损伤区域内节点的DF值与损伤程度高度线性相关.因此,在进行结构损伤程度识别时,可将某一特定的损伤程度设定为基准损伤工况,并通过有限元数值模拟手段计算该损伤程度下节点的DF数值.根据DF与损伤程度的线性关系,通过对比实测结构与基准损伤工况下的DF数值,即可对损伤程度进行较为准确的判断.

在利用本文所提方法对桁架结构进行损伤识别时,激振位置与测量噪声对识别结果的影响程度仍需通过数值模拟进行评判.根据数值模拟结果,与算例一相同的是,激振位置和测量噪声两种因素对于桁架结构损伤识别结果准确性与精度的影响仍然较小,因此可以忽略.限于文献篇幅,本文不再赘述.

3.2.4方法对比

现有的桁架结构损伤识别方法主要为基于结构振动特性的识别方法,如应变模态法[9]、柔度阵法[10]、单元模态应变差法[12]、柔度曲率幅值突变系数法[13]等.柔度曲率幅值突变系数法通过计算结构各阶模态柔度曲率幅值突变系数,对结构的损伤状况进行综合判断.本节以柔度曲率幅值突变系数法为例,分别应用本文所提方法与此方法对算例二的桁架结构进行损伤识别,并对识别结果进行整理对比.以工况1为例,两种方法的损伤识别结果对比如图15所示,其中图15a为仅利用一阶柔度曲率幅值突变系数进行对比;图15b为利用前三阶柔度曲率幅值突变系数的平均值进行对比.

a 一阶模态

b 前三阶模态

图15表明,应用现有损伤识别方法进行损伤位置识别时,未损伤区域处的节点柔度曲率幅值突变系数较大(如节点18、19),因此在实际工程应用中可能产生误判现象.此外,对于多位置损伤情况的识别以及损伤程度的判定,现有方法均不能较好地解决.本文所提方法在一定程度上提高了识别的精度,并在多位置损伤工况下以及损伤程度判定方面具有一定优势.

4 结论

本文基于ARX时间序列模型,提出一种新的损伤识别方法.根据时间序列模型基本方程推导出双输入单输出ARX模型一般表达式,通过联立多自由度体系动力方程确定模型阶数,根据未损伤结构的加速度时程数据,建立考虑结构动力特性的各节点ARX时间序列模型,预测未损伤结构的节点加速度时程数据,通过对比理论预测数据和实测数据,根据拟合程度计算损伤因子DF,进而进行结构损伤评估.数值算例分析结果表明,根据本文方法计算得到的各节点DF数值分布可用于损伤位置的识别定位,并可根据损伤位置相邻节点的DF数值大小识别损伤程度.该方法具有以下优点:

(1) 应用本文所提方法建立的ARX模型符合动态数据的实际变化规律,模型拟合程度较高,拟合模型所需的样本数量较少;

(2) 可仅针对疑似损伤区域或重点检测区域进行损伤识别,无需对整体结构进行模态参数识别,所需传感器数量较少;

(3) 实际工程应用中,可通过有限元数值模拟计算基准损伤工况的节点DF数值,并将实测结构所得DF与其进行对比,即可较为准确地判断损伤程度;

(4) 多位置损伤工况下仍具有较好的识别效果.