20世纪50年代开始，钢筋混凝土结构及构件的抗剪问题即成为工程界研究的热点和难点，通过大量试验数据对影响因素、破坏机理、分析模型等进行了系统、深入的研究，成果被纳入相关设计规程[1-5]，但因试验数据、模型结构、求解过程、影响参数等的不确定性，混凝土构件的受剪计算模型争议众多。因此，深入分析、解决不确定性的来源，对于从根本上提高结果的可靠性，降低主观判断误差，改进、优化分析模型具有重要的意义。
深受弯构件斜截面的剪切问题具有典型代表性，已有研究成果形成了大量的先验信息和计算模型，在此基础上保证相当的可靠度形成了相应设计方法，但承载力的预测与其准确值差距显著，即对剪切机理的认识仍不完善，各国规范及典型机理模型大多难以准确考虑不确定因素的影响。近年来，区别于传统频率统计方法，考虑研究对象的先验信息和样本空间的贝叶斯后验概率统计方法被引入进行承载力预测与分析[6-8]。
结合先验分布和样本信息，采用贝叶斯推断方法，通过计算所需后验分布函数的积分可完成对模型参数的估计，得到模型参数的均值和方差[9-10]。但在实际应用中，参数的后验分布大多为高维、复杂的分布函数，难以直接进行计算且结果偏差较大。Markov Chain Monte Carlo[11-14](MCMC)方法能有效地解决该积分计算困难的问题，其将高维积分的计算问题用随机模拟的方式来解决，MCMC方法的基本原理即从所需的后验分布中抽样，产生一系列后验样本，最后计算后验样本的统计值，为贝叶斯方法的后验参数分析提供了新思路。因此，本文在贝叶斯理论的基础上，根据大量的钢筋混凝土深受弯构件受剪试验研究，利用MCMC方法，考虑各种不确定性影响因素，编辑R语言，建立基于MCMC方法的深受弯构件受剪承载力概率模型，根据预定的置信水平确定深受弯构件的受剪承载力特征值，从而为深受弯构件受剪承载力的概率极限状态设计和概率安全性评估提供依据，为受剪问题的分析提供新方法。
1 概率模型的建立
1.1 MCMC方法简介
MCMC[15]方法的核心思想是通过建立一条平稳分布为π(x)的马尔科夫链，再对π(x)进行抽样，最后基于这些样本作各种统计推断。根据随机过程理论[16]，马尔科夫链是一个随机变量序列{X(0),X(1),…}，X(t+1)由条件分布π(x|X(t))产生，且X(t+1)的产生只与X(t)有关，与以前的状态{X(0),X(1),…,X(t-1)}无关，而本文中的后验样本刚好满足这一性质。若该链满足遍历性、不可约性，则无论初始值X(0)取什么，X(t)都会收敛到平稳分布，即马尔科夫链有收敛解。
则MCMC方法的基本步骤如下所示。
1) 在B上构造一条能收敛到平稳分布π(x)的马尔科夫链；
2) 从B中的某一点X(0)开始，用上述产生的马尔科夫链进行模拟，通过对后验分布π(x)抽样可得一系列后验样本：X(0),…,X(n)；
3) 进行蒙特卡洛积分，对抽样得到的样本进行积分：
根据大数定理，此积分可写成以下形式：
最后可得任意函数f(x)的后验期望和方差估计值为：
式(2)便是蒙特卡洛积分，n为通过模拟产生的总样本数，m为马尔科夫链达到平稳时的样本数。
MCMC随机模拟方法包括Gibbs[17]抽样方法与Metropolis-Hastings[18]算法，本文综合应用了上述两种抽样算法。
1.2 后验分布密度函数的产生
根据贝叶斯理论[19]，综合先验信息、总体信息与样本信息，取无信息先验分布
可得后验分布函数为：
根据Gibbs抽样方法可得式(3)的满条件分布如下所示：
1.3 概率模型
根据贝叶斯理论[20]，可基于已有受剪承载力计算模型建立深受弯构件受剪承载力概率模型：
式中：X表示钢筋混凝土深受弯构件受剪承载力影响因素的向量形式；Θ=(θ,σ)为拟合试验数据时所需要的模型参数；Vd为已有深受弯构件受剪承载力计算公式；γ(X,θ)为偏差修正函数；ε为正态随机变量；σ表示模型进行修正后仍存在的误差。对概率模型式(6)取对数可得：
式中，偏差修正项γ(X,θ)函数可以表示为：
通过对产生的后验分布进行随机模拟，可得到概率模型参数（θ,σ）的估计值，代入式(7)中可得概率模型。本文选取的深受弯构件已有受剪承载力计算模型分别为我国《混凝土结构设计规范》
美国ACI 318规范、加拿大
规范及欧洲规范EC2。
2 基于MCMC方法的抽样计算
2.1 试验数据整理
本文收集到的国内外深受弯构件受剪试验数据为645组，具体见参考文献[21]，根据理论和需要，式(8)中的hi(x)可确定为：h1(x)=ln2为修正常数项，h2(x)=ln(fcu/f y)，考虑混凝土抗压强度fcu和纵筋屈服强度fy的影响；h3(x)=ln(a/h0)考虑了剪跨比a/h0的影响；h4(x)=ln(l0/h)考虑了跨高比(l0/h)的影响；h5(x)=ln(b/h)考虑了构件截面尺寸(b×h)的影响，h6(x)=ln(eρv )和h7(x)=ln(eρh)分别考虑了竖向和水平向腹筋的(ρv,ρh)的影响，h8(x)=ln(ρ)考虑了纵向配筋率(ρ)的影响。
2.2 抽样过程
根据式(4)对参数θ 进行Gibbs抽样，X和Y为数据，θ 为未知参数，且θ=(θ(1),θ(2),…,θ(8))，选取参数θ 的一个初始值θ(0)，若第i次迭代开始时θ的值是θ(i-1)，则迭代算法如下所示。
1) 从满足条件分布
中抽取一个样本
2) 从满足条件分布
中抽取一个样本
3) 从满足条件分布
中抽取一个样本
4) 对i=1, 2, …, n重复以上各步，进而可得后验样本
，再通过式(1)和式(2)计算后验样本的统计值。
根据式(5)对参数σ进行Metropolis-Hastings算法抽样，先找到一个和后验分布相近的分布函数，称为建议分布密度函数：
从马尔科夫链中选取一个初始值σ（0），设第i次迭代开始时σ 值是σ（i-1)，则迭代过程如下所示。
1) 首先从建议分布密度函数q(σ;σ(i-1))中抽取一个待选样本σ*。
2) 计算接受概率：
3) 以概率α(σ;σ(i -1),σ*)接受σ(i)=σ*，以概率1-α(σ;σ(i -1),σ*)接受σ(i)=σ(i -1)。
4) 重复步骤1) ～步骤3) n次，进而可得后验样本σ(1),σ(2),…,σ(n)，再根据式(1)和式(3)计算后验样本的统计值。
上述两种抽样过程都通过R语言在计算机中实现。根据整理好的试验数据，计算出式(7)中的参数
然后利用R语言[22]软件对参数(θ,σ2)进行马尔科夫链-蒙特卡洛随机模拟，具体抽样流程图如图1所示。
图1 抽样流程图
Fig.1 The flow chart of sampling
2.3 计算结果分析
在随机模拟过程中，输入多个初始值，产生一条马尔科夫迭代链，在运行n=50000次迭代分析后，去除前m=2000个值，根据式(1)、式(2)可以得出待估参数(θ,σ2)，的估计值，图2为迭代形成马尔科夫链的基于我国规范的模拟轨迹图和后验密度函数图，在图2(a)～图2(h)中左边均为待估参数的模拟轨迹图，右边均为参数θ的后验密度函数图。
以图2(a)中左图为例，横坐标代表迭代次数，本文取50000次，纵坐标为参数θ1的估计值，迭代开始时先任取一初值，在经过50000次迭代算法后，参数θ的值会趋于一个固定值，表明该马尔科夫链最终收敛；以图2(a)中右图为例，横坐标表示估计值θ1的取值范围，模拟环境是后验样本N=5000个，估计值θ1的频率宽度为0.01962，纵坐标表示其累积概率，由图可知该后验分布密度函数曲线接近于正态分布的函数曲线，参数估计值会在最大概率处取到它的固定值，其他图为参数θ2~θ8的模拟结果，图中内容与θ1一致。由图可知，利用R语言模拟出来的马尔科夫链在运行50000次迭代分析后，后验样本会趋于一个固定值，表明该链最终达到收敛，该值就是待估参数的统计值。
图2 模型参数的模拟轨迹图和后验密度图
Fig.2 Simulated trajectory map and posterior density map of model parameters
表1~表4为基于四国规范的概率模型参数(θ,σ2)的运行结果，其中包括了(θ,σ2)的均值、方差、MCMC误差和参数估计值的2.5%、97.5%分位数。由表可知，由MCMC方法产生的模型参数估计值的方差都在0.05左右，MCMC误差不到0.1%，表明了此随机模拟方法能精确地求解待估参数的估计值，并且具有较高的可信度。
表1 基于GB5100―2100的受剪承载力概率模型参数运行结果
Table 1 Computational results of probabilistic shear strength model parameters based on GB 50010―2010
模型参数 θ1 θ2 θ3 θ4 θ5 θ6 θ7 θ8 σ2均值 0.7037 0.09419 -0.34357 0.12155 0.13143 -0.01479 -0.19380 0.32987 0.07543方差 0.05379 0.02788 0.0314 0.03791 0.02687 0.0274 0.01917 0.02169 0.0043 MCMC误差 0.00084 0.00039 0.00044 0.00054 0.00038 0.00039 0.00027 0.00031 0.00006 2.5%分位数 0.5123 0.04023 -0.4051 0.04753 0.07796 -0.06879 -0.23102 0.28807 0.06763 97.5%分位数 0.8974 0.1486 -0.28324 0.19541 0.18313 0.03901 -0.15564 0.37263 0.08438
表2 基于ACI 318的受剪承载力概率模型参数运行结果
Table 2 Computational results of probabilistic shear strength model parameters based on ACI 318
模型参数 θ1 θ2 θ3 θ4 θ5 θ6 θ7 θ8 σ2均值 -0.6492 -0.25013 -0.07604 -0.35147 -0.19265 -0.11734 0.02937 0.18213 0.07478方差 0.05181 0.02776 0.03127 0.03775 0.02676 0.02728 0.01909 0.0216 0.00426 MCMC误差 0.00087 0.00039 0.00044 0.00053 0.00038 0.00039 0.00027 0.00031 0.00006 2.5%分位数 -0.5034 -0.30387 -0.13731 -0.42518 -0.24589 0.06357 -0.00769 0.1405 0.06705 97.5%分位数 -0.7529 -0.19596 -0.01597 -0.27793 -0.14116 0.17091 0.06736 0.2247 0.08366
表3 基于CSA A23.3―04的受剪承载力概率模型参数运行结果
Table 3 Computational results of probabilistic shear strength model parameters based on CSA A23.3―04
模型参数 θ1 θ2 θ3 θ4 θ5 θ6 θ7 θ8 σ2均值 0.5087 0.02378 0.83417 -0.00789 -0.13415 0.11199 0.0098 -0.06441 0.08964方差 0.07582 0.03098 0.03459 0.04077 0.0292 0.03023 0.02096 0.02366 0.00507 MCMC误差 0.00087 0.00039 0.00044 0.00053 0.00038 0.00039 0.00027 0.00031 0.00006 2.5%分位数 0.3145 -0.03627 0.76639 -0.0881 -0.19146 0.05273 -0.03144 -0.11066 0.08024 97.5%分位数 0.6772 0.08477 0.9025 0.07238 -0.07679 0.1716 0.05105 -0.01799 0.10011
表4 基于EC2的受剪承载力概率模型参数运行结果
Table 4 Computational results of probabilistic shear strength model parameters based on EC2
模型参数 θ1 θ2 θ3 θ4 θ5 θ6 θ7 θ8 σ2均值 -0.6619 -0.25046 -0.07786 -0.35023 -0.19097 0.11819 0.02969 0.18146 0.07458方差 0.05242 0.02772 0.03123 0.0377 0.02672 0.02724 0.0191 0.02157 0.00425 MCMC误差 0.00087 0.00039 0.00044 0.00053 0.00038 0.00039 0.00027 0.00031 0.00006 2.5%分位数 -0.4723 -0.30412 -0.13904 -0.42384 -0.24412 0.06449 -0.00731 0.13989 0.06687 97.5%分位数 -0.8576 -0.19636 -0.01787 -0.27679 -0.13956 0.17168 0.06763 0.22397 0.08343
图3为对模型参数抽样产生的后验频数直方分布图，横坐标表示待估参数具有一定组距的取值范围，纵坐标表示其样本数。由图可见，模型参数在进行50000次迭代分析后，其估计值的后验频数大多集中在模拟结果值附近区域，与模拟轨迹图和后验密度函数图中模型参数θ的值基本一致，说明马尔科夫链蒙特卡洛方法在对参数进行抽样模拟和迭代分析时有较高的可信度，并且有效解决了运用一般数理统计方法估计模型参数时所存在的系统误差问题。根据后验样本的2.5%和97.5%分位数，以文献[20]中编号为5～30的一组深受弯构件为例，可得受剪承载力概率模型计算值具有置信水平为50%、80%和95%的特征值分别为162.0 kN、189.3 kN和222.6 kN。图4为由本文方法得到的深受弯构件受剪承载力的累计分布密度函数曲线，由图可见，基于MCMC方法建立的概率模型能确定不同置信水平下深受弯构件的受剪承载力特征值。将参数的模拟结果代入式(7)，并进行指数运算，整理可得基于中国规范的深受弯构件受剪承载力概率模型见式(14)，根据各参数影响显著性采用文献[9]中建议方法，剔除对计算结果影响较小的参数，即得简化概率模型，见式(15)，同理可得其他规范计算结果，见式(16)～式(18)。
图3 模型参数的后验频数直方图
Fig.3 Posterior frequency histogram of model parameters
图4 不同置信水平下的受剪承载力特征值
Fig.4 Characteristic values of shear strength at different confidence levels
美国ACI 318规范：
加拿大CSA规范：
欧洲EC2规范：
3 试验验证
四国规范受剪承载力计算值和概率模型计算值与试验值的比值的统计结果见表5，由表可见：概率模型计算值的均值更接近1.0，标准差较小，即基于四国规范所得得概率模型计算值较规范计算模型更接近试验破坏值，且离散性较小，精度较高。图5为四国规范剪力计算结果与概率模型剪力计算结果相对于混凝土强度的偏差对比，Vtest/VMCMC表示深受弯构件受剪承载力试验值和基于MCMC方法所得受剪承载力概率模型计算值的比值。由图可见：概率模型计算值与四国规范剪力计算值分布基本一致，但概率模型计算值更接近试验破坏值且分布更集中，离散性较小，如图5(a)中，Vtest/VGB主要分布在1.060～1.743之间，Vtest/VMCMC主要分布在0.739～1.337，说明经过MCMC方法模拟的计算值更接近于1.0，减小了由影响因素带来的偏差。进一步证明了MCMC方法的有效性，概率模型的离散性和随机性均显著减小。
图6为基于中国规范的深受弯构件受剪承载力概率模型计算值与中国规范以及试验结果的样本分布图，由图可见：中国规范的计算结果略为保守，且低于试验值，而概率模型计算值较接近试验值，说明经MCMC方法修正后的受剪概率模型计算精度较高，也进一步说明了本文建议方法的准确性和合理性。
表5 深受弯构件受剪承载力的模型计算值与试验值比值
Table 5 Comparison of biases between shear models and tests
计算模型 GB 50010-2010 ACI 318 CSAA 23.3-04 EC2规范计算值 MCMC.GB 规范计算值 MCMC.ACI规范计算值 MCMC.CSA 规范计算值 MCMC.EC2均值μ 1.401 1.038 1.254 1.113 1.754 1.046 1.254 1.047标准差σ 0.342 0.289 0.451 0.328 0.918 0.326 0.451 0.292
图5 概率模型计算值与规范计算值对比
Fig.5 Comparison of biases between shear models based on MCMC and the four codes
图6 基于GB50010―2010的剪力值对比图
Fig6.Comparison of biases between models and tests based on GB 50010―2010
4 结论
(1) 基于Bayesian-MCMC方法，结合645组试验数据和贝叶斯理论，利用R语言对深受弯构件受剪承载力概率模型参数进行随机模拟，建立概率模型，能够根据预定的置信水平确定深受弯构件的受剪承载力特征值。结果表明：概率模型计算结果和试验结果吻合良好，与四国规范抗剪承载力计算结果相比，更接近试验值，且精度比较高，离散性小。
(2) 基于MCMC方法得到的受剪承载力概率模型是在进行了50000次迭代分析而产生的结果，能合理地解释影响参数的不确定性，可信度较高。建议方法在土木工程领域有一定的推广性，能够为该类构件的承载力分析预测提供新思路、新方法。
[1]GB 50010―2010, 混凝土结构设计规范[S].北京: 中国建筑工业出版社, 2010.GB 50010―2010, Code for design of concrete structures[S].Beijing: China Architectural ＆ Building Press,2010.(in Chinese)
[2]GB 50011―2010, 建筑抗震设计规范[S].北京: 中国建筑工业出版社, 2010.GB 50011―2010, Code for seismic design of buildings[S].Beijing: China Architectural ＆ Building Press,2010.(in Chinese)
[3]ACI Committee 318, Building code requirement for structural concrete (ACI318―14) and commentary(ACI318R-14) [S].Farmington Hius, MI, USA:American Concrete Institute, 2005.
[4]CSA A23.3―04, Design of Concrete Structures [S].Mississauga, Canadian: Canadian Standards Association,2004.
[5]EN 1992-1-1: 2004, Europe Code2.Design of concrete structures [S].London, UK: British Standards Institution,2004.
[6]吴涛, 刘喜, 邢国华.基于贝叶斯理论的钢筋混凝土柱受剪承载力计算[J].工程力学, 2013, 30(5): 195―201.Wu Tao, Liu Xi, Xing Guohua.Study on the shear capacity of reinforced concrete column based on Bayesian theory [J].Engineering Mechanics, 2013,30(5): 195―201.(in Chinese)
[7]吴涛, 黄凯, 刘喜, 等.钢筋混凝土深受弯构件受剪承载力概率模型研究[J].工程力学, 2015, 32(11): 210―217.Wu Tao, Huang Kai, Liu Xi, et al.Study on probabilistic shear strength model for deep flexural member [J].Engineering Mechanics, 2015, 32(11): 210―217.(in Chinese)
[8]Kim J, LaFave J M.Probabilistic joint shear strength models for design of RC beam-column connections [J].ACI Structural Journal, 2008, 105(6): 770―780.
[9]Song J, Kang W H, Kim K S, et al.Probabilistic shear strength models for reinforced concrete beams without shear reinforcement [J].Structural Engineering &Mechanics, 2010, 11(1): 15―38.
[10]Gardoni P, Kiureghian A D, Mosalam K M.Probabilistic capacity models and fragility estimates for reinforced concrete columns based on experimental observations[J].Journal of Engineering Mechanics, 2015, 128(10):1024―1038.
[11]余波, 陈冰, 吴然立.剪切型钢筋混凝土柱抗剪承载力计算的概率模型[J].工程力学, 2017, 34(7): 136―145 Yu Bo, Cheng Bing, Wu Ranli.Probabilistic model for shear strength of shear-critical reinforced concrete column [J].Engineering Mechanics, 2017, 34(7): 136―145.(in Chinese)
[12]吴然立.锈蚀钢筋混凝土柱的抗剪机理及概率抗剪承载力研究[D].南宁: 广西大学, 2016.Wu Ranli.Study shear mechanism and probabilistic shear strength of corroded reinforced concrete columns[D].Nanning: Guangxi University, 2016.(in Chinese)
[13]Beck J L, Au S K.Bayesian updating of structural models and reliability using markov chain monte carlo simulation [J].Journal of Engineering Mechanics, 2002,128(4): 380―391.
[14]Muto M, Beck J L.Bayesian updating and model class selection for hysteretic structural models using stochastic simulation [J].Journal of Vibration & Control, 2008,14(1/2): 7―34.
[15]张忠诚.MCMC方法及应用[D].武汉: 湖北大学, 2008.Zhang Zhongcheng.MCMC method and its applications[D].Wuhan: Hubei University, 2008.(in Chinese)
[16]龚光鲁, 钱敏平.应用随机过程教程[M].北京: 清华大学出版社, 2004: 15―30.Gong Guanglu, Qian Minping.Applying the random process tutorial [M].Beijing: Tsinghua University Press,2004: 15―30 (in Chinese)
[17]赵琪.基于MCMC方法的贝叶斯统计推断[J].中国科技信息, 2012(10): 64―65.Zhao Qi.Bayesian inference based on the MCMC method [J].China Science and Technology Information,2012(10): 64―65.(in Chinese)
[18]Chib S, Greenberg E.Understanding the Metropolis-Hastings algorithm [J].American Statistician, 1995,49(4): 327―335.
[19]朱慧明, 韩玉启.贝叶斯多元统计推断理论[M].北京:科学出版社, 2006: 30―42.Zhu Huiming, Han Yuqi.Bayesian multivariate statistical inference theory [M].Beijing: Science Press, 2006: 30―42.(in Chinese)
[20]Chetchotisak P, Teerawong J, Yindeesuk S, et al.New strut-and-tie-models for shear strength prediction and design of RC deep beams [J].Computers and Concrete,2014, 14(1): 19―40.
[21]刘喜, 王洁, 吴涛, 等.深受弯构件抗剪影响参数显著性分析[J].工程力学, 2017, 34(7): 126―135.Liu Xi, Wang Jie, Wu Tao, et al.Significance analysis of influence parameters on shear capacity of deep flexural member [J].Engineering Mechanics, 2017, 34(7): 126―135.(in Chinese)
[22]韩明.贝叶斯统计学及其应用[M].上海: 同济大学出版社, 2015: 5―15.Han Ming.Bayesian statistics and its applications [M].Shanghai: Tongji University Press, 2015: 5―15.(in Chinese)