基于自相关峰值变化率的古木结构损伤识别

2020-08-27侯艳芳胡卫兵

结构工程师 2020年3期

侯艳芳胡卫兵杨佳

（1.西安建筑科技大学土木工程学院，西安710055；2.陕西工业职业技术学院土木工程学院，咸阳712000）

0 引言

中国古建筑是人类文化遗产中独具特点又弥足珍贵的一支，蕴含了浩瀚无涯的精神文化，它既能客观地反映历史的发展脉络及其面貌，又能将许多文字无法记录的史料保存至今。中国古建筑内涵深厚而形神皆备，它具有许多特质，在技术方面善于运用木结构，将木材技术发挥到极致，为世界其他文明所罕见［1］。在千百年的使用过程中，由风雨侵蚀、自然灾害、材料老化和长期荷载作用等因素共同作用，中国古木建筑已存在不同程度的损伤和破坏。近现代交通（地铁等）所产生的高频微幅振动，在原有损伤的基础上不断累积损伤并使抗力衰减，从而在自然灾害发生甚至是在正常环境作用下结构发生损坏的可能性大幅提高。一座古建筑就是一本立体的且无法复制的史书，一旦消失其损失无法弥补。若能利用有效的方法和手段对古木建筑进行损伤识别和评定，进而及时控制和修复损伤，避免结构的损坏，对古建筑的修缮和保护具有重要意义。

传统的检测方法通常采用人工激振，但是该方法需要较大的激振能量，容易引起古建筑木结构的损伤，反而对研究起到反作用。因此，利用环境激励（车辆、风、机械振动等）下的结构振动响应数据进行损伤检测成为结构健康监测中的热点［2］。

于哲峰［3-5］等提出了互相关函数幅值向量的概念，以某测点k为参考点，将其他各测点与参考点之间的互相关函数的最大值rkl组成一个幅值向量，进行了随机激励下基于互相关函数幅值向量的结构损伤检测的实验研究及损伤定位方法研究。杨智春等［6］提出了基于互相关函数幅值向量和连续小波变换联合检测损伤的方法。党晓娟等［7-8］提出了基于互相关函数幅值向量与统计学原理相结合的损伤统计检测方法，通过对一个复合材料层合悬臂板进行分层损伤检测的数值模拟，验证了该方法的有效性；并将互相关函数幅值向量与离散小波分解结合起来，提出了一种损伤检测新方法。

雷家艳［9］以相邻两测点的动力响应时程为基础，定义了新的互相关函数幅值向量和损伤因子，并通过试验数据分析验证了该方法的可行性。由于没有了参考点的牵制，对损伤的局部反映更直观。

本文根据相关函数的原理，计算各振动测点的自相关函数，并将所有的峰值组成自相关函数幅值向量；计算两个完好结构的自相关函数幅值向量之间的互相关系数，以此作为结构损伤识别的基准值。计算结构损伤前后各测点自相关函数峰值的变化率，确定结构的损伤位置。通过算例分析，验证了该方法的可行性。

1 相关函数的基本理论[10-12]

1.1 相关与相关函数

相关是信号分析中一个非常重要的概念。相关函数有互相关函数和自相关函数，分别表示两个时间序列之间或同一个时间序列在任意两个不同时刻取值之间的相关程度，即：互相关函数描述的是两个信号x（t）和y（t）在任意两个不同时刻t1、t2的取值之间的相关程度；自相关函数表示信号x（t）在某一时刻的数据与另一时刻的数据之间的依赖关系。

对于连续型随机过程，信号x（t）与y（t）的互相关函数记为

式中，T为所取信号的时间过程，τ表示两个信号取值的时间差，称为滞后。

若随机变量是离散型随机序列，则两个信号之间的互相关函数表示为

式中，N为随机序列的离散值数；τ为时间间隔。

如果y（t）=x（t），表示信号x（t）与自己的时延信号x（t+τ）的相似程度，也即是自相关函数。

连续型随机变量的自相关函数为

离散型随机序列的自相关函数为

φ（τ）的大小可以定量的衡量两组信号之间的相似程度。φ（τ）的数值不仅与两个信号自身的特点有关，还与两组信号之间的时间间隔τ有关。φ（τ）的最大值所对应的τ=τi值表示，当两个信号的时延为τi时，它们之间最相似。

1.2 自相关函数幅值向量及互相关系数

根据相关函数的概念，随机过程两个状态之间的相关程度，总不会大于在同一个点上的状态与其自身的相关程度，所以自身与自身强相关。因此，求随机振动结构各测点的自相关函数能更明确反映出结构损伤前后各测点峰值的变化程度。

自相关函数在τ=0 时取得最大值。这是因为此时x（t+τ）就是与x（t）自身求相关，相关程度最大。所以，φxx（τ）的最大值为φxx（0）。

对于随机振动的结构，若同时测得n个测点的动力响应（加速度、速度、位移等），则各测点动力响应的自相关函数为

当τ=0 时，记各测点xi的自相关函数φxi（τ）的最大值为：φxi（0）=δii，由所有的最大值组成的自相关函数幅值向量为

在后面的计算中，为表达方便，自相关函数幅值向量简记为A。

表示两组信号相关程度的一个重要特征数值是互相关系数：

互相关系数ρxy表达的是两组信号或向量之间的相关程度，其取值范围为［0，1］。当x（t）=y（t），即信号与自身相关时，ρxy=1；当x（t）和y（t）不相关时，ρxy=0。

互相关系数ρ越大，相关程度就越高；反之，ρ值越小，也就说明两组向量之间的相关性越小。因此，可以根据损伤结构与完好结构的ACFA之间互相关系数的大小判别结构是否发生了损伤，即可以将互相关系数作为结构损伤识别的损伤因子。

1.3 自相关函数幅值向量的固定形态

设x（t）是遍历过程的一个样本函数，它是定义在（-∞＜t＜∞）区间内的一个非周期函数，一般不满足绝对可积条件，进行傅里叶变换时需要引入截尾函数，即x（t）是指在有限时间区间上的截尾函数。x（t）和x（t+τ）的傅里叶变换为

它们的傅里叶逆变换为

其中，X*(ω)是X(ω)的共轭。

设结构所受激励荷载f(t)的傅里叶变换为

结构频率响应函数的第j行为：

则有

将式（13）代入式（10），得

可见，在F(ω)一定的情况下，自相关函数φxx(τ)只与各激励点之间的频响函数H(ω)有关。因此，自相关函数幅值向量中的各分量元素具有固定的比例关系，也即是具有固定的形态。

2 木结构古建筑的结构特性

中国古建筑木结构有别于现代的钢结构、钢筋混凝土结构。尤其在结构抗震上，中国木结构古建筑体现出了“以柔克刚、摩擦耗能、滑移隔震”的抗震原则，是高度有机的结构体系。这与现代结构大多采用“加强结构、硬抗地震”的抗震原则有很大的不同［13］。

木结构古建筑在结构上最显著的一个特点是梁柱等构件之间采用榫卯节点连接。榫头和卯口间因摩擦滑移而具有耗能作用，这是木结构古建筑具有较好抗震性能的一个重要原因［14］。

方东平等［15］采用杆件单元和半刚性节点单元建立木结构古建筑的有限元模型，通过定义和引入反映古建筑木结构榫卯节点特性的半刚性节点单元，对榫卯节点的力学性能作了定量研究，证实了木结构古建筑的榫卯节点是半刚性的。

薛建阳［16］将榫卯的力学模型比拟为变刚度弹簧单元，进行了动力时程分析。结果表明：由于大震下柱底与础石的滑移和榫卯与斗拱的摩擦阻尼作用，大幅度减少了上部结构的地震响应。因此，中国古建筑木结构具有良好的耗能与减震性能。

3 算例分析

3.1 有限元模型的建立

以西安钟楼为工程背景，建立木结构古建筑有限元模型时，选取一榀木框架进行损伤识别分析。木梁长4 m，截面尺寸为300 mm×700 mm；木柱高6 m，截面直径为500 mm，梁柱采用Beam188单元。基于榫卯节点的半刚性特点，采用Combin14 弹簧单元模拟梁柱的榫卯连接。每个节点共有六个自由度，其中三个平动自由度UX、UY、UZ 和三个转动自由度ROTX、ROTY、ROTZ。榫卯节点刚度采用文献［17］中根据钟楼监测数据得到的近似值：1×1010kN·m/rad。木材的弹性模量为1×1010N/m2，泊松比为0.25，密度为410 kg/m3，采用Rayleigh定义的黏性比例阻尼。

木结构框架有限元模型及节点编号如图1所示。

图1 有限元模型（a）及节点编号（b）Fig.1 Finite element model and nodes number

韩广森［18］利用列车荷载经验分析模型，计算出了西安地铁2 号、6 号线单线运行时的作用荷载，得到了不同速度下的荷载值，并通过MATLAB得到了激励的时程曲线和频谱曲线。将经过土层传播折减后的激励作用到钟楼的有限元模型上，并对2 号、6 号单线运行下的钟楼响应进行了分析。

本文在建立木框架有限元模型进行分析时，随机激励选择西安地铁单线运行时，速度分别为V1=20 km/h 和V2=40 km/h 时的荷载值作为随机激励样本1 和随机激励样本2 施加在柱脚节点上。激励时长取T=12 s，时间步长为Δt=0.01 s。

3.2 结构的损伤识别

选取木框架上节点编号为1-9 的一跨梁进行分析。梁长4 m，网格划分时单元大小取0.5 m，共8 个单元。由于建立有限元模型时梁两端边单元的节点与梁柱的榫卯节点重合，本文选取梁上中间6 个单元进行损伤识别与定位分析，单元编号为①-⑥，对应的节点编号为1-7。

首先对完好结构分别施加随机激励1 和随机激励2，提取完好结构在两个随机激励下的两组加速度动力响应。计算完好结构在两个随机激励下的两组自相关函数幅值向量分别为A1和A2，经归一化后的两条曲线如图2所示。

图2 A1和A2的ACFA曲线Fig.2 ACFA curves of A1and A2

由图可以看出，完好结构在两个不同随机激励下的两条自相关函数幅值向量曲线几乎吻合。这说明对于随机振动的结构，在结构完好状态下，各测点动力响应的自相关函数幅值向量具有固定的形态。因此，可以采用随机振动结构各测点动力响应的自相关函数幅值向量判别结构的损伤。计算图中两条自相关函数幅值向量曲线之间的互相关系数，如表1所示。

表1 A1和A2之间的互相关系数Table 1 Cross correlation coefficient of A1 and A2

将A1和A2之间的互相关系数ρ12=0.999 997 7作为损伤识别的基准值。若损伤结构的ACFA 与完好结构的ACFA 之间的互相关系数明显小于该基准值，则认为结构发生了损伤。

通过降低梁单元的刚度模拟结构的损伤。以④单元为例，将单元④的刚度降低30%，提取结构损伤状态下梁上各测点的加速度响应，并计算各测点加速度响应自相关函数的最大值δii，组成损伤结构的自相关函数幅值向量，记为Ad。损伤结构Ad与完好结构A1、A2的自相关函数幅值向量曲线如图3所示。

图3 损伤结构与完好结构的ACFA曲线Fig.3 ACFA curves of A1、A2 and Ad

从图中可以看出，单元④损伤时的自相关函数幅值向量曲线与完好结构A1、A2的曲线相比有明显的波动。完好结构的自相关函数幅值向量曲线取为A1，即将A1记为Au。单元④损伤时，经计算，损伤结构Ad与完好结构Au之间的互相关系数如表2所示。

表2 Ad与Au之间的互相关系数Table 2 Cross correlation coefficient of Ad and Au

取ρdu=0.999 02。同样的方法，可以计算出其他各单元分别发生损伤时，与完好结构之间的互相关系数，如表3所示。

表3 损伤结构与完好结构之间的互相关系数Table 3 Cross correlation coefficient of Ad and Au

将表中各单元损伤时的互相关系数与基准值做对比，可以发现，损伤结构与完好结构之间的互相关系数有明显的降低，由此可以判断结构发生了损伤。因此，使用结构损伤前后各测点自相关函数幅值向量之间的互相关系数作为损伤因子可以有效判别结构的损伤。

3.3 损伤定位分析

由前述损伤因子的变化可知，当结构发生损伤时，损伤结构与完好结构之间的相关性会降低，也即是损伤结构的自相关函数幅值向量曲线会发生波动。分别将各单元的刚度降低30%，绘制各单元发生损伤时与完好结构的自相关函数幅值向量曲线，如图4所示。

从图中可以看出，当结构发生损伤时，自相关函数幅值向量的曲线会发生波动。梁两端边单元损伤时，各测点自相关函数峰值的波动较小；中间单元损伤时，梁上各测点自相关函数的峰值波动幅度较大，从而导致结构损伤前后自相关函数幅值向量之间的相关性降低。

为了衡量损伤前后曲线的波动程度，计算各测点损伤前后峰值的变化率，即：

当各单元分别发生损伤时，梁上各测点自相关函数峰值损伤前后的变化率如图5所示。

从图中可以看出，当梁上某一单元发生损伤时，该单元的两个节点的峰值变化率最大。如单元②损伤时，损伤前后相邻的两个节点2和3的峰值变化率最大，此时对应着单元②两端的两个节点。其他单元分别损伤时，也满足这种对应关系。所以，可以利用结构损伤前后各测点峰值变化率的最大绝对值来确定结构的损伤位置，即结构损伤后峰值变化率最大的相邻两个节点的区域对应着结构的损伤单元。

由于信号自身与自身的强相关性，使用各测点的自相关函数幅值向量比两点之间的互相关函数幅值向量能更明确表达结构损伤前后峰值的变化程度，不需要再对各测点的相对变化进行二次差分，进而避免了由于二次差分造成的节点与单元的位置偏移，使损伤前后各测点峰值的相对变化更加突出，识别结果更加突出明显。

3.4 测量噪声对损伤识别的影响

为了模拟实际工程中噪声对测量信号的影响，向节点加速度中分别加入信噪比SNR=20 db、10 db、0 db 的高斯白噪声。分别计算不同噪声水平影响下，结构损伤前后各测点自相关函数峰值的变化率。为了方便对比，将无噪声和三组不同噪声水平的自相关函数峰值的变化率归一化，如图6所示。

从图中可以看出，单元③损伤时，无噪声情况下，单元的两个节点3和4的峰值变化率明显高于其他各测点，能够准确进行损伤定位。随着信噪比的降低，其他各测点的峰值变化率逐渐接近损伤单元的两个节点峰值变化率的数值，但当SNR=0 dB时，仍能对损伤单元③进行准确定位。因此，对于受随机振动的结构，采用各测点损伤前后自相关函数峰值的变化率作为损伤定位指标具有较好的抗噪声干扰能力。