基于MCMC法的混凝土坝坝体坝基变形模量随机反演

程 井，李培聪，李同春，袁 平

(1. 河海大学 水利水电学院，江苏 南京 210098；2. 中冶长天国际工程有限责任公司，湖南 长沙 410007)

坝体、坝基物理力学参数的选定是坝工设计和分析坝体应力、变形以及裂缝形成机理的基础。传统反演方法主要建立在确定性分析基础上。然而，由于混凝土坝工作环境的复杂多变性，坝体变形是一个随机变量，变形过程是一个随机过程，观测值只是一个样本序列或样本值的实现，即观测的变形具有不确定性，只能从概率意义来研究变形的均值与方差[1]。因此，考虑不确定性的随机反演方法更符合实际情况。

目前，混凝土重力坝主要通过对原型观测资料分析，建立确定性数学模型，反求大坝的材料参数以及某些结构特性[2]。近年来，部分学者开始考虑水工结构反演中的不确定性问题。杨杰等[3]提出基于最大熵原理的贝叶斯不确定性反分析法；苏怀智等[4]在传统的确定性模型中引入区间数因子，考虑了不确定性因素对反演结果的影响；杜永峰等[5]基于模糊理论提出考虑测量数据不确定性的结构物理参数识别方法；雷鹏等[6-7]进一步研究了参数不确定性区间反演。贝叶斯方法是一种考虑不确定性因素的有效反演方法，马尔可夫链蒙特卡罗方法(MCMC)是目前贝叶斯方法的标准抽样方法，在岩土工程、水资源优化以及土壤学等学科已有学者进行相关研究[8-10]。然而，对于许多复杂问题，似然函数难以获得，因此，Marjoram等[11]首次提出了无似然函数的马尔可夫链蒙特卡罗方法，在此基础上，Daniel等[12]通过调节马尔可夫链的容差，提高该方法的计算效率，简化了许多复杂的工程问题。

本文基于Bayesian理论，考虑参数先验信息分布及观测值的不确定性，采用无似然函数的马尔可夫链蒙特卡罗（MCMC without likelihoods）方法进行弹性模量的随机反演，研究参数后验分布的统计特性，并分别对观测值出现波动情况（人为观测误差及拟合中出现误差）与观测值进行比较，评价模型误差和预测精度。

1 基于现代贝叶斯理论的参数随机反演

贝叶斯原理中，令X为一连续型分布的随机变量，其先验概率密度函数为π(X)，则后验概率密度函数f(X∣D)为：

式中：X为随机变量（待反演参数）；D为观测样本；f (X∣D) 为随机变量的后验分布；f (D∣X) 为似然函数，数值上等于D在X上的条件分布；π(X) 为随机变量（参数X）的先验分布。

1.1 马尔可夫链蒙特卡罗方法（MCMC）

对于复杂问题，若先验分布和后验分布不满足共轭分布条件，常规贝叶斯方法难以得到后验分布估计值（均值、标准差）的解析式，因此，需要借助于数值方法或近似方法进行模拟[8]。为解决该问题，采用马尔可夫链蒙特卡罗方法计算后验分布。而MCMC方法的基本思想则是通过随机抽样建立一个最终平稳分布为所求后验分布的马尔可夫链，通过马尔可夫链得到后验分布的样本，进而得到对应的期望值和标准差。常用的MCMC方法有Gibbs抽样、Metropolis抽样等基本算法。采用的无似然函数MCMC算法具体步骤[11]如下：

(1) 初始化马尔可夫链初始状态X0=x0；

(2) 在i时刻（第i次循环），马尔可夫链状态Xi=xi，由转移核q(Xi→Xi+1)采样Xi+1；

(3)根据Xi+1及已知先验计算模型M，产生随机样本D′；

(4) 如果D′=D（已知抽样样本），进入下一步，否则从第（2）步重新开始；

(5) 计算接受率h：

(6) 若满足接收条件，则接收Xi+1，否则从第(2)步重新开始。

抽样过程中，抽样终止以马尔可夫链收敛为前提，因此，马尔可夫链收敛与否对模型参数后验估计具有重要影响。BGR诊断法[13]是马尔可夫链收敛性诊断常用方法之一，该方法基于区间统计特性及长度进行诊断，定义诊断指标为：

式中：L 为 马尔可夫链总的序列区间长度或统计特性；l为马尔可夫链单链区间长度均值或统计特性均值。若马尔可夫链收敛，值应接近于1。

1.2 观测数据（样本）不确定性分析

坝体位移是一个随机变量，变形过程是一个随机过程，观测值只是一个样本序列或样本值的实现，因此具有不确定性。根据文献[2]大坝位移统计模型表达式为：

式中：δH，δT，δθ分别为水压分量、温度分量、时效分量。考虑到混凝土坝工作环境的复杂多变性以及各种不确定因素（如测量误差、仪器精度、拟合误差等），假设对于某个区间内的每个值有：

式（6）表明，大坝位移由两部分组成，一部分是由水压分量δH、温度分量δT和时效分量δθ等组成，另一部分ε ～N(0, σ2)为不确定性因素造成的随机误差。则对于水压分量有：

式中：H为上游水位测值与坝底高程之差；ai为水压因子回归系数；εH为水压分量中不确定性因素造成的随机误差。

1.3 坝体坝基变形模量随机反演

采用MCMC without likelihoods算法进行参数随机反演，具体步骤：

（1）根据观测坝体位移变形数据统计回归分析[2]，得到不同水位下实测位移统计模型水压分量D以及不确定性因素造成的随机误差 εH～N(0,σH2)；

（2）确定坝体、坝基变形模量X的先验分布π(X)；

（3）由于观测的变形具有不确定性，由统计模型分析得到的坝体位变形与实际观测值必然存在差距，因此，抽样过程中必须考虑不确定性因素造成的随机误差εH～N(0,σH2)；

（4）以后验概率密度函数f (X∣D)为目标函数，采用MCMC without likelihoods算法产生坝体、坝基变形模量随机样本，其中通过随机误差εH～N(0,σH2)筛选相应的有效样本，2.4节通过人为调整随机误差，进一步研究了随机误差对坝体、坝基变形模量分布的影响；

（5）将有效的马尔可夫链样本作为后验概率密度函数f (X∣D)的样本，计算后验分布f (X∣D)的估计值（均值，标准差）。

2 工程算例

2.1 工程概况

龙滩水电站前期工程于2006年10月开始下闸蓄水，2008年12月初步达到正常蓄水位。大坝坝顶高程382 m，建基面高程216.43 m，坝高165.57 m，坝顶宽度14 m。以11号坝段为典型坝段，根据实测坝体水平位移反演坝体、坝基变形模量。如图1所示，二维有限元模型范围为：向上、下游各延伸250 m、坝基建基面以下延伸235 m。有限元节点布置时，尽可能考虑将位移测点安排在单元节点上，采用四节点单元，共2 383个节点，2 244个单元。

图 1 11 号坝段有限元模型及垂线测点Fig. 1 FEM mesh and vertical monitoring points of 11# dam section

参数反演时，选取正垂线测点 PL11-1（高程 342.00～379.00 m），PL11-2-2（高程 270.00～342.00 m）与PL11-3（高程222.75～270.00 m）2010年5月28日至2011年8月15日的实测顺河向位移进行分析，依据式（4）分离坝体位移水压、温度及时变等主要分量，各垂线的复相关系数R分别为0.964，0.951和0.977，反演分析中输入不同水位下实测位移统计模型水压分量D，即对模型依次施加水位为335，340，345，355和360 m的上游静水压力，结果如表1所示。

表 1 不同水位下实测位移统计模型水压分量DTab. 1 Water pressure component D given by statistical model at different water levels

观测值只是一个样本序列或样本值的实现，也就是说观测的变形具有不确定性。因此，观测值必须考虑各种不确定因素。由图2(a)所示，拟合值与实测值存在误差，误差来源于上述各种不确定因素。根据1.2节所述引入水压分量的随机误差 εH～N(0,σH2)，图 2(b)，(c)和 (d)分别为 222.75，270.00和 342.00 m 高程拟合总体位移及其分量。

图 2 龙滩观测位移拟合值与实测值Fig. 2 Fitted values and measured values of observation displacements of Longtan Dam

2.2 单参数反演分析

假定坝基变形模量为确定量，坝体弹性模量为待反演随机变量。采用无似然函数的MCMC算法来反演单个参数，初步认为坝体弹性模量服从区间为[10，60]的均匀分布，模拟时马氏链样本数取20 000个。

图3为无似然函数MCMC抽样样本。图4为正态分布Q-Q检验图，其中横坐标为分位数，纵坐标为样本值；图中点近似为一条直线，说明样本基本服从正态分布。图5为BGR诊断法马氏链收敛性诊断图，通过比较可以看出样本数量达到7 500时，诊断指标稳定在1附近，所以在以下分析中有效链长度取为10 000～20 000。图6比较了参数的先验分布、后验分布概率密度曲线及其正态拟合曲线，结果表明：（1）坝体弹性模量均值由先验分布的45 GPa变为后验分布均值32.09 GPa，基本符合龙滩重力坝实际坝体弹性模量；（2）弹性模量后验分布区间比先验分布明显收窄，坝体弹性模量反演结果的标准差为3.81 GPa，变异系数为0.119。

图 3 MCMC抽样样本Fig. 3 Samples in MCMC

图 4 正态分布Q-Q检验图Fig. 4 Q-Q diagram of normal distribution

图 5 马氏链收敛性诊断Fig. 5 Convergence diagnosis of Markov chains

图 6 弹性模量后验分布概率密度分布曲线Fig. 6 Probability density curves of posterior distribution

2.3 多参数反演分析

假定坝体弹性模量及坝基变形模量均为待反演随机变量，采用无似然函数的MCMC算法反演多个参数。反演时，初步认为两参数均服从区间为[10，60]的均匀分布，模拟时马氏链样本数取20 000个。

图7为坝体弹性模量、坝基变形模量MCMC抽样样本，其Q-Q检验图表明两个变量的马氏链样本均服从正态分布。图8比较了两个参数的先验分布、后验分布概率密度曲线及其正态拟合曲线，结果表明：（1）坝体弹性模量均值由先验分布的45 GPa变为后验分布均值33.49 GPa，坝基变形模量的均值由先验分布35 GPa变为后验分布均值29.53 GPa，基本符合龙滩重力坝实际情况；（2）坝体弹性模量、坝基变形模量标准差分别为3.78 和1.90 GPa，变异系数分别为0.11和0.06，标准差与变异系数均明显缩小，说明参数的不确定性显著减小。

图 7 MCMC抽样样本Fig. 7 Samples in MCMC

图 8 参数后验分布概率密度分布曲线Fig. 8 Probability density curves of posterior distribution

采用经典的最小二乘法反演坝体、坝基变形模量，与本文提出的反演方法进行对比，在同一坝段，同样观测数据情况下，两种反演方法的结果比较如表2所示。结果表明，两种算法都能较好地给出参数的反演结果，最小二乘法尽管误差平方和更小，但缺少对不确定因素的考虑，只给出了一个确定值的结果；而MCMC算法充分考虑了不确定因素的影响，反演结果能更好地反映坝体、坝基变形模量的分布规律。

表 2 两种反演方法结果Tab. 2 Comparison between two inversion methods

2.4 坝体坝基变形模量敏感性分析

位移的不确定性决定了坝体弹性模量及坝基变形模量反演结果的不确定性。通过模拟不同的σH进行坝体弹性模量及坝基变形模量的敏感性分析。正垂线PL11-1（高程342.00 m）与PL11-3（高程222.75 m）测点实测样本标准差分别为σH= 0.12 mm和σH= 0.24 mm。图9为不同的观测样本标准差（随机误差）与坝体弹性模量、坝基变形模量标准差的关系，易发现随着观测样本随机性（波动）的增大，反演的坝体、坝基变形模量变异系数相应增大。

图 9 观测样本标准差与坝体弹性模量、坝基变形模量变异系数关系Fig. 9 Relationships between standard deviations of observation samples and variation coefficients of elastic modulus of dam body and deformation modulus of dam foundation

3 结 语

针对混凝土坝材料力学参数反演中存在大量不确定性的问题，提出了高混凝土重力坝坝体弹性模量与坝基变形模量的MCMC随机反演方法，给出了无似然函数的马尔可夫链蒙特卡罗法在工程中的实际应用，并研究了坝体变形不确定性对反演结果的影响。研究结果表明：（1）该随机反演方法可以得出所需反演参数（坝体弹性模量、坝基变形模量）的分布；（2）反演参数的后验分布变异性较先验分布有明显降低；（3）在考虑变形观测值离散性变化情况下，后验分布变异性与观测值离散性呈正相关关系。

MCMC随机反演方法计算简便、高效，具有很强的工程实用性，为混凝土坝参数随机反演提供了重要计算方法。但有时为了得到较为精确的参数尾部分布，需要的抽样数较大，这就导致有限元计算的总耗时巨大，因此需要结合有限元原理对该算法进行优化。