基于贝叶斯的钢筋混凝土梁模型修正方法

宋彦朋，陈 辉,，黄 斌

(1.武汉工程大学 邮电与信息工程学院，湖北 武汉 430073；2.武汉理工大学 土木工程与建筑学院，湖北 武汉 430070)

在土木工程结构中，计算模型通常采用有限元格式。土木工程结构尤其是城市道路桥梁结构，主要是由钢筋混凝土构成。混凝土或钢筋混凝土结构相对于钢结构而言，材料内部的非均匀性、各向异性、样本参数的个体间差异显著。因此有限元仿真模型与实际结构存在比钢结构更为显著的差异。为了减小差异并建立足够精确的有限元模型，有必要根据测量数据更新有限元模型的参数[1～3]。

使用测量数据来修改结构的计算模型是一个非常活跃的研究方向。传统的基于动力的模型修正方法主要是基于结构固有频率、振型、模态曲率、柔度矩阵等[4]。这些方法基本都是基于特征值方程的动力学反问题。在动力试验中往往只能得到前几阶模态数据，因此在基于有限元方法的模型修正中往往会求解不适定方程，造成病态解[5]。很多学者采用不同的方法来避免病态解。例如，Ren等[6]采用模态截断奇异值分解算法来修正自由模态的钢筋混凝土实验梁，并将修正后的模型作为基准模型用于该梁的损伤识别，取得了很好的效果；Hua等[7]将Tikhonov正则化用于平面数值桁架的有限元模型修正来避免修正方程求解的病态问题；李英超等[8]将奇异值截断正则化算法(Truncated Singular Value Decomposition，TSVD)用于三维导管架平台的缩尺简化模型的模型修正，很好地避免了多解问题。从以上方法中可以看出，正则化技术是处理方程病态问题的一种非常有效的手段，但正则化参数或者奇异值截断数以及测量误差对求解结果也有较大影响，另外，这些文章多半都是基于数值算例或者单一材料的实验算例，对钢筋混凝土结构的研究较少。另一方面，基于Beck等[9,10]提出的基于贝叶斯框架的模型修正或结构参数识别的方法也取得了巨大的发展，很多学者也将其应用到混凝土结构的有限元模型修正中。该方法最大的特点是通过随机抽取参数样本，正向计算结构响应并和测量结果比较来选择合理的待修正参数，从根本上避免逆向求解问题。易伟建等[11]运用贝叶斯方法，通过测量频率和振型MAC值构建目标函数，对一个四层混凝土框架结构进行了模型修正，并将修正后的模型作为基准模型用于该结构的损伤识别。房长宇等[12]用贝叶斯方法基于测量频率对预应力混凝土梁进行了模型修正。Wan等[13]基于测量频率，将贝叶斯方法用于实际混凝土人行天桥的模型修正中。从上述基于贝叶斯框架的模型修正方法中可以看出，该方法最重要的两个问题是目标函数的构建和抽样方法的选取。

考虑到混凝土材料的特殊性以及施工及环境的影响，混凝土结构各区域材料参数差异较其他材料显著，因此需要尽可能地对其作全局修正，但这会导致修正参数较多，传统的吉布斯(Gibbs)抽样方法或MH(Metropolis-Hasyings)抽样方法限制了该方法的使用；另外，要进行全局修正，通常需要对大小发生变化的修正参数进行定位，但和振型相关MAC值只能体现整体模态差异，无法体现局部差异，而振型是位置坐标的函数，因此在目标函数的构建上应该充分使用该测量信息。基于贝叶斯方法的最新进展，本文考虑采用贝叶斯模型修正方法对钢筋混凝土实验梁进行模型修正，以避免模型修正中的病态问题。同时基于测量频率和模态来构建目标函数，采用更为高效的DRAM(Delayed Rejection Adaptive Metropolis)抽样方法来提高抽样的遍历性。为了确认结构真实模型，文章通过对比实际物理特征对修正结果进行了确认。

1 贝叶斯模型修正的基本原理

基于贝叶斯方法的随机模型修正是结合先验知识(主观信息)和测试数据(客观信息)，采用MCMC(Markov Chain Monte Carlo)方法推断修正参数的后验概率分布。其表达式如下：

=cp(x/θ)π(θ)∝p(x/θ)π(θ)

(1)

式中：x为观测信息；π(θ)为待修参数θ的先验分布；p(x/θ)为在θ给定下的条件分布，通常称为似然函数；c是一个不依赖于θ的常数因子。

在土木工程或者机械结构中，一般选用结构的弹性模量、刚度、密度、质量或某些几何尺寸作为修正参数，标记为θ，观测值记为DN，N为观测次数。将上述贝叶斯原理引入到有限元模型修正中，可将修正参数的后验概率分布写为：

p(θ/DN)=cp(DN/θ)p(θ)

(2)

式中：p(θ/DN)为结构后验概率；p(DN/θ)为似然函数；p(θ)为待修正参数θ的先验分布；μ0为先验参数均值，covσ0为先验参数的协方差矩阵；y为测量响应；y(θ)为经模型仿真计算(如有限元计算)得到的响应；covy为测量信息协方差矩阵；c为与θ无关的常数。根据贝叶斯假设采用无偏广义先验，p(θ)值应为1，待修正参数的后验概率可以写为：

(3)

式中：J(θ)为目标函数，可表示为：

(4)

2 目标函数的建立

结构的响应通常为频率和振型，在基于贝叶斯方法的模型修正方法中，常用的是频率、振型和MAC值或者它们的组合来构造目标函数。实际模态测试中，频率的测试较为容易和精确，而且不会受测点缺失的影响；但模态振型的测试中仅前几阶低阶模态能较为准确测量到。同时在测试过程中也经常会出现部分传感器信号缺失或者故障的情况发生，这会造成模态信息不完整。因此，本文从实际模态测试角度构造目标函数以修正结构参数。

考虑每次测量的频率与振型和抽样计算结果的误差，可将目标函数(4)进一步定义为：

(5)

将式(5)表示的目标函数J(θ)代入式(3)中可以得到待修正参数的后验概率密度的表达式。为了找到修正参数的后验概率密度函数的最大值，可用MCMC方法近似计算后验概率密度函数的分布，即采用MCMC抽样使目标函数J(θ)最小，最终可得到最优修正参数。

3 计算后验概率密度的DRAM抽样方法

在采用MCMC抽样计算待修正参数后验概率密度的过程中，产生Markov链的常用方法是Gibbs 抽样[14]和MH抽样[15]。这两种方法的不足是当参数较多时，采样过程平滑段出现逐渐增多，且参数之间易产生相互影响最终难以收敛。近年来，发展了一种新的抽样方法即DRAM抽样，它是一种更高效的抽样方法[16]，本文采用了这种抽样方法。该方法的抽样步骤如下：

(1)假设在t时刻，Markov链当前值为θt，然后根据建议分布q(θt,C0)随机抽取候选样本θt+1，其中C0为初始抽样方差。

(2)根据式(3)，计算接受概率p。

(3)从[0，1]均匀分布中随机产生一个变量u，当p>u时，接受候选样本θt+1，即θt+1=θt；当p≤u时，拒绝候选样本θt+1。

(4)根据二次建议分布q2(θt,C0)随机抽取候选样本θt+1，根据式(3)，计算接受概率p2，然后执行步骤(3)。

(5)当执行到第t+1步，且t+1≥N0(非适应抽样次数)时，更改建议分布的协方差矩阵为Ct+1，从(1)开始执行。其中任意抽样时刻t建议分布的协方差矩阵Ct满足：

(6)重复步骤(1)～(5)，直至Markov链趋于平稳后迭代终止。然后计算从收敛开始后修正参数样本的均值和方差等统计特征。

4 基于贝叶斯的钢筋混凝土梁有限元模型修正

4.1 修正参数的选择

(6)

式中：Δ[K]i为结构第i个单元的刚度变化量；αi为结构第i个单元的修正系数。为确保本文提出方法的有效性，在数值算例中验证了几种实际测量中常见工况下的模型修正效果。

4.2 基于贝叶斯的钢筋混凝土梁有限元模型修正

对钢筋混凝土梁的初始有限元模型进行计算，得到其自振频率和位移振型。同时对实际模型进行模态测试，分析并得到有效的测量信息。将初始模型自由度和测量自由度进行匹配，得到匹配后的初始有限元模型对应的频率和位移振型。将匹配后的所有振型进行质量归一化。最后将测量频率和计算频率以及归一化以后的位移模态代入式(3)中，用DRAM方法进行抽样并计算得到实际结构修正参数的后验概率，最终可以得到修正参数的优值。

4.3 修正结果验证

模型修正的目的是使得修正后结构的计算响应与测量结果一致，同时修正后的参数必须能和真实结构的物理特征对应。因此，在用本文提出的方法得到钢筋混凝土梁的修正参数后，首先应将修正参数代入修正后的有限元模型中计算响应，并和测量结果进行对比以评估修正结果的有效性；同时，为了确保修正结果具有真实意义，还应该对钢筋混凝土梁构造明显的物理特征(如可见的细小裂缝)，来进一步验证修正结果是否具有真实意义。

5 数值算例

取一悬臂梁，如图1所示，其密度ρ=7.8×103 kg/m3，初始弹性模量E=2.10×l011Pa。将梁划分为七个单元，假定悬臂梁弹性模量发生了变化，变化后各单元弹性模量真实值如表1所示。假定待修正参数为各单元弹性模量，相对于初始弹性模量E的变化量用参数αi(i=1,2,…,7)表示。实际模态测试中，测量信息往往有限，同时高阶模态测量误差较大。能有效使用的往往是少数低阶模态。因此，需要就测量模态阶数对参数识别结果的影响进行分析。为此，分四种工况讨论，每种工况的测量信息为：(1)结构前二阶频率和前二阶位移振型；(2)结构前三阶频率和前三阶位移振型；(3)结构前四阶频率和前四阶位移振型；(4)结构前四阶频率。

图1 悬臂梁有限元模型/mm

表1 各单元弹性模量真值

分别将以上四种工况下的测量频率和位移模态信息代入目标函数进行参数识别。仿真测量误差假定为变异系数为0.01的正态分布随机变量。四种工况下各识别参数样本的Markov链如图2所示。在工况1～3下，悬臂梁中待识别参数的均值如图3所示。

图2 四种工况下采用DRAM抽样的各单元弹模变化参数样本的Markov链

图3 悬臂梁中待识别参数的均值

从图2a～2c中可以看出，不同测量信息下，三种基于频率与振型识别结果中，所有单元的修正参数均在真值附近波动，并呈现稳定的收敛现象，说明本方法识别结果具有很强的稳定性；另外，还发现可靠的测量信息越丰富，Markov链的稳定性越高。另外，图2d的结果说明，仅基于前四阶测量频率的Markov链并不收敛，说明仅依赖于前几阶低阶频率的贝叶斯方法无法有效识别较多结构参数的变化，可能需要更多高阶频率参与修正，或采用多层MCMC方法。综合整个仿真算例来看，本文提出的方法在测量的少数低阶频率与振型信息下，能有效修正结构模型。

6 混凝土梁的模型修正试验

通过对一实际钢筋混凝土梁进行模型修正，以验证本文提出方法的有效性。试验梁几何尺寸和配筋等情况如图4a所示。实际模型如图4b所示。混凝土的力学参数为：密度ρ= 2400 kg/m3，弹性模量E=26 GPa，剪切弹性模量G= 0.4E，泊松比μ=0.33。充分考虑梁的剪切效应，选用Timoshenko梁单元类型梁单元，梁被分成110个单元用于有限元计算以保证计算结果收敛。修正过程中，考虑到修正参数过多容易使抽样结果发散，因此将修正模型作如下处理：将110个单元中的1～10单元作为一个修正参数(单元)，以此类推。因此总共修正参数(单元)为11个。

图4 钢筋混凝土试验梁尺寸和动力测试

为了准确地验证修正结果的物理意义，将试验梁缓慢加载至轻微裂缝出现后停止加载，使之产生明显的裂缝这一物理特征，但同时要保证试验梁处于正常使用极限范围之内。在裂缝处做好红色标记以便于和修正结果比较。最终在单元4，6，7，8处有不同程度的裂缝产生，如图4b所示，为非对称裂缝。其中4，6单元的裂缝深度在梁深度的一半附近，7单元裂缝达到裂缝深度40%附近，基本在跨中，8单元裂缝深度接近梁深度的20%，各裂缝宽度均不超过0.5 mm；将试验梁放置在软橡胶垫块上进行自由模态测试。模态试验中，采用六个加速度传感器分三批对12个测点进行测量，传感器布置如图4b所示，加速度响应图和频谱图如图4c，4d所示。采用工作模态分析(Operational Modal Analysis，OMA)识别结构频率和模态，所得结果如表2和图5所示。各阶频率的测量误差变异系数为0.02，在1 h内总共进行了12次模态测试。将12次测量模态各测点变异系数中的最大值作为振型的变异系数，为0.02。测量误差均假定为零均值的正态分布。其中振型已进行质量归一化。为了提高修正精度和计算效率，认为梁两端的单元刚度基本一致，跨中出现裂缝的单元刚度变化较大。因此，考虑到试验梁质量保持不变，将1～3单元和9～11单元的刚度相对于初始有限元刚度的变化量分别作为两个待修正参数，中间的4～8单元的刚度相对于初始有限元刚度的变化量作为另外五个待修正参数，总共7个待修正参数。取测量结果中的前四阶频率和位移振型代入目标函数对待修正参数进行抽样。值得注意的是在式(4)中调整目标函数中频率与位移振型的权重比例为10000以使两者数量级一致。

表2 测量频率 Hz

图5 混凝土梁初始有限元模型和实际测量的振型

图6为本文方法得到的Markov链，从图中可以看出，抽样超过3000次后各单元待修正参数开始收敛，因此取燃烧期为3000。根据收敛后的Markov链进行最大后验估计，得到修正参数均值和标准差，如图7所示。对比图7中的刚度变化和图4b中静力加载后记录的裂缝位置和深度，可以发现：刚度下降的位置和裂缝位置对应，裂缝深度占梁高度的比例和刚度下降量比较接近，单元4中裂缝深度为梁高的50%，识别结果显示刚度下降约为40%；单元5中裂缝深度为梁高的50%，但裂缝宽度很小，识别结果显示刚度下降接近10%，单元6中裂缝深度为梁高的60%，识别结果显示刚度下降量接近57%，单元7中裂缝深度为梁高的40%，识别刚度下降量接近38%；单元8中裂缝深度为梁高的15%，识别刚度下降量接近14%；同时注意到梁首端三个单元刚度略有下降，末端三个单元的刚度均有所增大。整个修正结果和试验梁的实际刚度变化位置和大小情况基本吻合。

图6 试验混凝土梁的DRAM抽样结果

图7 混凝土梁的贝叶斯模型修正结果

为了进一步对修正结果进行验证，将修正后的参数代入新的有限元模型中进行计算，得到修正后的结构响应。图8比较了修正前后试验梁的前六阶计算频率相对于测量频率的误差。从图8中可以看出，修正前结构频率和测量频率之间的最大误差接近26%，最小频率误差也接近10%；修正后试验梁前四阶计算频率和测量频率之间的相对误差很小，未参与修正过程的第五、六两阶频率相对于测量结果也有大幅下降。这进一步说明修正结果是有效的。

图8 修正前后频率误差

图9比较了修正前后的前四阶计算位移振型和初始有限元振型。从图9中可以看出，测量的振型和初始有限元振型有一定的误差，并且随着振型阶次的提高，测量振型和初始值之间的误差越来越大。这再次显示了仿真和试验的区别；图10是修正前后前四阶振型MAC值，可以发现，经过修正后的计算振型和测量结果更加接近，前四阶振型MAC值均更接近1。另外值得注意的是，从图9e可以看出，第五阶测量振型明显偏离实际情况的允许范围，这时因为高阶的振型由于仪器精度等原因的限制，会产生脱离实际的误差，因此若作为测量信息使用则会产生错误的修正结果。

图9 初始有限元和修正前后的前五阶位移振型

图10 修正前后模态MAC值

7 结 论

本文针对钢筋混凝土梁模型修正中的方程病态问题及模态测试中测点信息不完整和模态信息不完备的实际情况，采用贝叶斯方法进行模型修正，目标函数的构建充分考虑了模态测试的特点，抽样方法采用高效的DRAM抽样方法。数值算例结果表明，基于DRAM抽样和模态测量数据的贝叶斯模型修正方法能有效修正结构有限元模型参数，即使在仅用前两阶模态的情况下也能进行准确修正，说明该方法具有很好的鲁棒性和灵活性。模态试验结果显示：该方法得到的修正结果能使试验梁结构修正后的前六阶计算频率与测量值吻合；同时使前四阶振型的MAC值也进一步提高并非常接近1；且静力加载导致的裂缝位置和程度与修正的位置和大小有很好的对应关系，说明了该方法能有效应用于实际钢筋混凝土梁结构的模型修正。