贝叶斯理论在既有砌体损伤定位中的应用研究

夏恒，崔鸿知

（长沙理工大学 土木工程学院，湖南长沙 410114）

0 引言

目前，部分既有砌体[1]属于历史文化保护遗产，其历史人文价值与日俱增，但材料老化、荷载作用、结构疲劳、超期服役等因素会导致砌体结构出现退化和损坏的现象，由于结构复杂、损伤位置隐蔽等情况，很难对结构进行准确的损伤定位，假如对其损伤放任不管，随着损伤的累积，将可能造成重大安全事故和经济损失，故对这些砌体结构进行损伤识别研究具有重要的意义。

在20 世纪50 年代，我国的传统检测方法（筒压法、扁顶法、原位单砖双剪法、原位轴压法、点荷法）属于定期的人工检测，在进行检测时还需预知结构损伤的大致位置，对于一些难以测量的部位，其人工检测的危险程度也会增加。到了20 世纪80 年代，随着科学技术的发展，基于动力学的损伤识别[2]方法逐渐成为主流，其原理是结构发生损伤时，相关参数会发生改变，通过在结构上布置相应的传感器可对结构的健康状况进行实时监测。

在基于动力学的损伤定位方法上， Pandey[3]提出绝对变化量指标，根据损伤前后应变模态差的绝对值进行损伤定位。董聪[4]根据结构应变的一阶变分关系，推导出了应变模态相较于位移振型和自振频率更能有效、准确地进行损伤定位分析，并且在进行试验后数据表明：损伤位置应变模态的变化十分明显，且变化值的大小与损伤程度成正比关系。Vanik[5]提出基于参数与贝叶斯理论相结合的在线健康监测方法。顾培因[6]提出一种在无法获取损伤前的模态数据下的一种直接指标法，该方法可以直接利用损伤后的应变模态参数进行损伤定位。殷志翔[7]利用应变模态差对索网结构实现了损伤识别。邹云峰[8]提出了基于应变模态响应重构的损伤识别方法。吴韬[8-10]等人将动力特性用于砌体结构的损伤识别。

研究发现，目前砌体结构的损伤定位研究大部分属于确定性分析，即损伤前的参数都有准确的数据。但在现实中，大部分的砌体结构因为各种原因，很难获取损伤前的参数。本文以湖南省长沙市湖南大学财经学院的省级保护性建筑砌体柱为例，考虑测量误差和模型误差的影响，以改良的应变模态差值公式嵌入到贝叶斯模型修正的框架中。该方法可以有效地处理损伤定位中的不确定性因素，能准确地识别结构损伤位置，并且还能将此次获得的后验信息作为下次更新的先验信息，从而实现长期的动态观测。

1 贝叶斯基本理论及MCMC 方法

1.1 贝叶斯基本理论

贝叶斯统计推断[11]方法是引入概率学的知识来做出判断，它是一种承认生产实践的继续性方法，将相关先验信息（如历史可靠数据）参与决策。推断基础来自于贝叶斯公式，其公式如下：

式（1）中P（θi）是事件θi 发生的先验分布，与事件D 无关，先验信息主要来源于可靠的历史数据和专家经验，假设先验分布为正态分布N（μ，λ2），可以利用先验信息确定μ 和λ2的值；P(D|θi)是指在事件θi 发生的基础上，事件D 的条件分布；是一个全概率数，起一个正则化因子的作用；P(θi|D)是在事件D 发生的基础上，事件θi 的分布，也被称为后验分布。一般情况下后验概率是难以表达的，对后验概率的计算可以采用MCMC-马尔科夫链蒙特卡罗法。

1.2 MCMC-马尔科夫链蒙特卡罗法

在马尔科夫链蒙特卡罗法中进行抽样是至关重要的，可以降低其样本的相关性。不同的抽样方法，将产生不同的马尔科夫链蒙特卡罗算法。1987 年，Durane[13]将马尔可夫链蒙特卡罗法与分子动力学相结合，提出了混合蒙特卡罗算法。混合蒙特卡罗算法结合了梅特罗波利斯-黑斯挺算法、吉布斯算法以及确定性分子动力学方法，可以用于高维复杂分布的抽样模拟，改进了建议概率分布的产生机制，避免陷入局部解。

混合蒙特卡罗法的算法流程[14]如下：

（1）链指标k=0，初始化位置向量x0和动力向量p0。

（2）链指标k=k+1，从正态分布N（0，M）中随机产生一个新的动力向量pt。

（3）赋予蛙跳算法初值[x(0),p(0)=(xk-1,pt)]并允许L 时间步以获得一个建议状态（x*，p*）=[x(t+Lε)，pt+Lε)]。

（4）计算接受概率α=min[1,exp(-△H)],其中△H=H（ x*，p*）-H（xk-1,pt）。如果接受建议状态（xk，pk）=（x*，p*），否则（xk，pk）= （xk-1，pt）。

（5）重复步骤（2）—（4），直至到达预设的链长为止。

再根据后验样本可计算后验分布的各阶矩，进行相应的统计推断。

2 有限元模型修正方法

部分既有砌体由于年代久远，无法准确获取结构损伤前的参数，只能根据现场测量数据和相关历史资料初步建立无损的有限元模型，但所建立的模型与实际情况存在差异。为了修正有限元模型，本文把结构自身物理参数，如材料密度、弹性模量、以及泊松比等作为修正参数，以实际结构的测量值与有限元模型的动力响应特征值的相关信息来构建目标函数[15]，以相关信息差值作为修正依据。

式（2）中：y（r）表示为现实结构与有限元模型动力响应特征值的差值，可以是应力、应变、频率、振型、阻尼、应变模态等；ME表示为现实结构测量的数据 ；MA（r）表示为有限元模型的特征值，其中r 表示为修正的参数向量。有限元模型修正流程如图1 所示。

图1 有限元模型修正流程图

3 基于贝叶斯理论的损伤定位方法

3.1 基于应变模态的损伤识别原理

在基于动力学的结构损伤定位中，用于损伤定位的参数需要满足两个条件：一是对局部损伤敏感；二是位置坐标的单调函数[4]。

在文献[4]中证实了应变模态相较于其他参数更能有效、准确解决损伤定位问题。

由动力学方程：

文献[4]推导出结构应变的一阶变分函数为：

式（4）中{Δε}为应变的变化；[Δψ]为应变模态的变化；[Уr]结构自振频率矩阵；[φ]为位移模态矩阵；[ψ]为应变模态矩阵；[ΔУr]结构自振频率变化；[Δφ]位移模态的变化；其中{Δε}和[Δψ]的变化在位置上存在明确的对应关系，是具有进行损伤识别的能力。

损伤定位的公式以应变模态差的绝对值为依据，其公式为：

式子（5）中[ψI]和[ψD]表示为无损结构和有损结构的应变模态矩阵；i 表示坐标位置；j 表示模态阶数。在损伤位置处，应变模态差Δψ会发生突变，Δψ的大小表示损伤的严重程度，Δψ 恒大于0。

本文综合考虑不确定性因素带来的影响，实际的应变模态差的绝对值调整为：

式（6）中,f1、f2分别是测量误差和模型误差，假设f1、f2都服从正态分布，则其均值分别为μ1、μ2，方差分别为σ21、σ22。

3.2 似然函数与后验分布

由改良的定位公式（6），可推出似然函数表示为：

确定似然函数和先验分布后，后验分布[16]则为：

式（9）中，P(D)为边缘似然函数，其值为固定常数。其中参数θ 中的每一个θj，其边缘概率密可以通过积分获得为：

3.3 损伤识别流程

本文将砌体结构识别流程分为三个主要部分：（1）先基于历史资料初步建立有限元模型，再激励砌体结构获取相关参数数据，进行动力特性分析，最后以相关信息差值对有限元模型进行修正；（2）将现场收集的数据作为先验信息，再通过马尔科夫链产生候选样本数据，依据混合蒙特卡罗法计算接受概率，确定样本集；（3）通过贝叶斯方法，计算应变模态的后验分布，判断是否收敛。具体流程图如图2 所示。

图2 损伤识别流程图

4 砌体结构损伤定位研究

4.1 有限元模型材料参数和边界条件

烧结砖砌体柱来源于湖南省长沙市湖南大学财经学院的省级保护性建筑的第一层柱体，结构尺寸为720 mm×470 mm×25 440 mm，如图3 所示。实体结构多处带有明显裂缝，且砌块间的砂浆风化程度严重，几乎散失了强度。损伤前的有限元模型依据现场测量数据、相关资料和GB 50003—2011《砌体结构设计规范》确定烧结砖为MU15、砂浆为M7.5，砌块密度设为1 800 kg/m3、弹性模量设为1 600（ff 为砌体抗压强度，单位为MPa），计算为3 312 MPa、体积剪切模量为1 577 MPa、剪切模量1 440 MPa、线性膨胀系数设为5×10-6/℃、收缩率设为-0.1 mm/m、摩擦系数设为0.70、热导率设为1.4 W/(m·K)、泊松比设为0.15。

图3 烧结砖砌体柱

有限元模型采用ANSYS Workbench 建模，单块砖的模型为240 mm×115 mm×53 mm，边界条件设置为：上部承受荷载为均布荷载，方向垂直向下无偏心，根据GB 50009—2012《建筑结构荷载规范》和设计资料，其均布荷载值设为20 kN/m2，底部采用 Fix supporied（固定约束），加速度设为9.8 m/s2，方向竖直向下。由于模型整齐，单元类型采用六面体SOLID186 单元，单元尺寸为50 mm，网格划分采用程序自动控制。有限元模型如图4 所示。

图4 有限元模型及检测点布置图

4.2 试验方案和数据结果

在实体结构上取五点分别为a、b、c、d、e，这五点处于结构正面，如图4 所示。其中a、b、c 三处是带有明显裂缝的位置，d、e 两处无裂缝，通过比较现实结构和有限元模型的应变模态差值判断结构损伤位置，差值越大说明损伤越严重。对实体结构a、b、c、d、e 测量点分别测量3 次，以其平均值作为先验样本（表1），再经过混合蒙特卡罗法计算，设置样本接受的概率，经过多次循环后，使样本值向目标值收敛。由于前面的样本值高度依赖于初始值，故舍弃前200 个样本。本文分析前四阶模态下应变模态差值的后验分布。

表1 实际结构检测点应变模态

有限元模型在修正参数后得到无损伤前应变模态云图（图5），检测点a、b、c、d、e 的应变模态如表2 所示，实际结构和有限元模型检测点的应变模态如图6 所示。

表2 有限元模型检测点应变模态

图5 前四阶应变模态云图

图6 实际结构与有限元模型检测点应变模态

从图6 中可以明显发现：测量点处应变模态值的大小与结构损伤程度成正相关。测量点（a、b、c）在一阶模态和四阶模态下的应变模态值明显高于二阶模态和三阶模态下的应变模态值，其应变模态差值在一阶模态和四阶模态下有着更为明显的峰值，故在损伤定位时，针对此结构采用一阶模态和四阶模态能更准确地识别损伤位置。在检测点（d、e）处，各阶模态上的应变模态差值相较于测量点（a、b、c），应变模态差值保持在稳定低值范围内。将五点的应变模态差值作为先验信息带入计算模型中，一阶模态下测量点处的应变模态差值统计结果如表3 所示。

表3 测量点处一阶应变模态差值

考虑测量误差和模型误差的影响（式6），其均值和标准差在不同的采样点取值不同，将修正的应变模态差值的均值带入马尔科夫链中，生成样本集，通过贝叶斯算法确定后验分布图（图7），具体数值如表4 所示。

表4 测量点处一阶应变模态差值的修正值

图7 各测量点应变模态差值的后验分布图

从图7 和表4 中可以发现经过贝叶斯推断后，其应变模态差值的均值和标准差有相对明显的修正，在带有裂缝的a、b、c点应变模态差值与没带裂缝d、e 点应变模态差值存在明显不同，在带有裂缝处的检测点其应变模态差值远远高于没有裂缝处的检测点，如图8 所示。通过应变模态差值的大小可以判断出a、b、c 点相较于d、e 点是存在损伤，并且其损伤的程度与应变模态差值成正相关。

图8 应变模态差值折线图

5 结论

贝叶斯理论属于一种概率统计理论，该方法可借助先验信息（包括可靠的历史数据以及专家经验），再基于总体信息和样本信息得出合理推断。本文将该方法与有限元方法共同运用于砌体柱的损伤定位研究，结论如下：

（1）砌体柱的损伤前有限元模型可以通过现场测量数据和相关资料以及经验进行初步的建立，利用实际结构与有限元模型的动力特征值的差值作为修正依据，以贝叶斯算法作为优化方法，做到信息差最小，完成有限元模型的建立。

（2）贝叶斯推断和有限元相结合的方法在结构损伤定位研究中具有可行性，贝叶斯理论可以有效地处理损伤定位中不确定性因素的影响，从而得出更为合理的推断。