贝叶斯框架下大坝渗流参数反演组合代理模型研究

余红玲，王晓玲，王 成，曾拓程，余 佳，盖世聪

（1.天津大学 水利工程仿真与安全国家重点实验室，天津 300072；2.天津市政工程设计研究总院有限公司，天津 300051）

1 研究背景

大坝渗流性态分析对于大坝安全稳定运行具有重要意义。渗流参数（如坝体和坝基岩土体的渗透系数）是进行大坝渗流性态准确分析的关键参数。大坝运行期间获得的渗压、渗流量等监测资料是大坝渗流性态的最直接反映，故基于渗流监测数据进行参数反演分析是获取渗流参数的有效手段。在现有众多参数反演方法中，贝叶斯反演方法能够考虑参数反演过程中的不确定性，因此，其被广泛应用于渗流参数反演研究中［1-3］。然而，贝叶斯反演方法需要大量调用渗流正演模型，计算耗时，效率较低。因此，有必要采用有效方法提高贝叶斯反演方法的计算效率。

基于机器学习算法的代理模型作为一种提高贝叶斯反演计算效率的可行方法，近年来受到越来越多的关注。代理模型是通过数学模型逼近一组输入变量（独立变量）与输出变量（响应变量）的方法，通过寻求输入输出变量间响应关系，可代替真实系统快速给出所求解［4］。常用的代理模型包括响应面法、径向基函数（RBF）、人工神经网络（ANN）、Kriging 模型、支持向量回归（SVR）、多元自适应回归样条（MARS）等。然而，单一的代理模型并不总是能捕捉到输入输出变量间的复杂非线性关系［5］。采用单一模型，容易忽略模型结构和模型参数选择的不确定性，从而低估模型预测的不确定性［6-8］。这种预测不确定性可能会传播到贝叶斯反演方法中，并影响渗流参数反演结果的准确性。组合代理模型可通过综合单一模型的不同特性，最小化单一模型的不利影响来处理这种预测不确定性，提高模型整体预测精度［7-8］。

组合代理模型提供了更加可靠和准确的预测结果。然而，在组合代理模型中包含一个表现不佳的单一模型通常会导致组合模型的整体性能恶化。因此，需要谨慎地选择组成组合代理模型的各个单一模型，以减少性能不佳的单一模型的不良影响。SVR 通过引入核函数将低维空间中的数据映射到高维特征空间，且基于结构风险最小化原理，泛化能力高［9］；Kriging 具有逼近和映射复杂非线性函数的能力，局部拟合精度高［10］；多元自适应样条回归（MARS）是一种快速、灵活、自适应的非参数回归模型，能有效克服高维数据的不连续问题［11］。因此，考虑到单一代理模型的预测不确定性，本文选择SVR、Kriging 和MARS 这三种机器学习算法，通过对各单一模型进行加权求和构建组合代理模型。

确定各单一模型的权重系数是构建组合代理模型的关键一步，通过为各单一模型分配权重可进一步减少性能不佳的模型的不良影响。目前组合代理模型的权重确定方法较为丰富，其大致可以分为三类。第一类，将设计空间划分为多个子域，对每个子域使用一组权重因子，如Yin 等［12］将设计空间分成多个子域，通过计算样本的均方根误差给每个子域分配权重系数，从而构建一个全局的组合代理模型。第二类，将权重计算过程处理为一个优化问题，如Goel 等［13］和Cheng 等［14］采用一种基于广义交叉验证均方根误差的启发式策略来确定组合代理模型的权重；Viana 等［15］和Christelis 等［16］通过最小化均方根误差来获得权重系数；Ye 等［17］提出基于最小化局部均方根误差的权重系数确定方法；Wang 等［18］采用差分进化-禁忌搜索混合优化算法得到最优权重系数。第三类，利用迭代计算策略确定权重系数，如Liu 等［19］提出先采用广义交叉验证均方根误差确定初始权重，再通过迭代计算获得最佳权重的方法。

上述权重确定方法均能获得较为合理和准确的权重系数。本文基于差分进化自适应Metropolis（Differential Evolution Adaptive Metropolis，DREAMZS）算法提出一种新的权重确定方法。DREAMZS算法是近年来提出的一种马尔科夫链蒙特卡洛采样方法［20］，它基于统计抽样理论和随机模拟思想，通过统计分析模型输出参数的均值、方差等估计量，拟合输出参数的概率分布情况，从而定量描述参数的不确定性［21］。相比于其他采样方法，DREAMZS算法可同时运行多条平行的马尔科夫链，根据各条链当前位置和过去状态的样本采用差分进化算法自适应地调整推荐分布，从而实现每条链的随机采样。DREAMZS算法具有较高的求解效率，并可以保持平衡性和遍历性，现广泛应用于不确定性分析问题中［22-24］。因此，本文采用DREAMZS算法随机采样的优势计算模型权重系数的随机分布函数，在考虑不确定性的条件下获得模型的权重系数。

本文在贝叶斯框架下，建立大坝渗流参数反演的组合代理模型，避免单一代理模型预测精度较低的缺陷，且组合代理模型与计算耗时的渗流正演模型相比，能够显著提高参数反演的计算效率。本文所建模型为大坝渗流参数反演提供了一种新思路，并可以应用于其他相关领域的参数反演。

2 研究框架

本文提出的贝叶斯框架下大坝渗流参数反演组合代理模型主要包括4 个部分：①数据集生成；②组合代理模型建立；③贝叶斯反演；④案例分析。研究框架如图1 所示。

图1 研究框架

（1）数据集生成。根据参数先验分布和取值范围，采用拉丁超立方抽样（LHS）方法生成待反演渗流参数的样本点，将样本点输入到耦合VOF 法的三维非稳态水气两相流渗流模型［25］中进行模拟计算以获得监测点处的水头模拟值，待反演渗流参数样本点和水头模拟值构成数据集，进一步可划分为训练集和测试集。

（2）组合代理模型建立。根据训练集训练SVR、Kriging 和MARS 三种机器学习模型，避免单一模型低估预测不确定性的缺陷，提高预测精度；基于DREAMZS算法计算出各模型的权重系数，通过加权求和的方法建立组合代理模型。

（3）贝叶斯反演。设置水头监测值和待反演渗流参数的先验分布；由于求解待反演渗流参数后验分布的贝叶斯公式中的分母部分存在积分区间难以选择，求解困难的问题，故采用马尔科夫链蒙特卡罗（MCMC）采样方法求解待反演渗流参数的后验分布，其中，将（2）中建立的组合代理模型替代计算耗时的渗流正演模型，从而提高反演计算效率。

（4）案例分析。将本文模型应用于工程实例，并结合监测数据对计算结果进行对比分析，从而证实该模型的有效性和优越性。

3 贝叶斯框架下大坝渗流参数反演组合代理模型

3.1 贝叶斯理论框架 在贝叶斯方法中，待反演渗流参数被视为随机变量，用相应的分布函数来表征。根据贝叶斯原理，在观测数据的支持下，待反演渗流参数的后验分布p(θ |y )可以表示为

式中：θ为待反演的渗流参数；y 为渗流观测数据；p( θ )为渗流参数的先验分布，可通过专家经验或历史资料获得；p( y|θ )为似然函数，表示渗流模型预测值和观测值的似然度。渗流模型预测值和观测值的误差ε计算为

式中f（θ）为渗流模型的预测值。误差ε通常被称为观测误差，一般假设其服从均值为零，协方差矩阵为∑的多元正态分布ε～N（0，∑），其中可以根据监测仪器的精度确定，I表示单位矩阵。则高斯似然函数可表示为

式中： |∑|为∑的秩；n 为渗流观测点的维数。

式（1）中的分母为归一化常数，其积分区间比较难以选择，通常采用马尔可夫链蒙特卡洛（MCMC，Markov Chain Monte Carlo）方法进行抽样计算［26］，本文采用MCMC 方法中的DREAMZS算法进行抽样计算。

MCMC 方法需要抽取大量样本并把每一组样本代入渗流正演模型中进行计算，对于单次正演计算耗时较长的复杂渗流模型，计算效率较低。本文通过构建组合代理模型来替代贝叶斯反演求解过程中的渗流正演模型，从而提高贝叶斯反演的求解效率。

3.2 基于DREAMZS 算法的组合代理模型 组合代理模型可以综合各单一代理模型的优势，提高代理模型的整体预测精度，组合代理模型的预测值fen( θ )可以表示为

式中：θ为输入变量，本研究中为待反演的渗流参数；M 为单一模型的个数；wi为各单一模型的权重系数； fi（θ）为各单一模型的预测值。

SVR、Kriging 和MARS 的详细介绍可分别参考相关文献［27-29］。确定权重系数是构建组合代理模型的关键一步，为了获得更加可靠的权重系数，本文采用DREAMZS算法来确定各模型的权重系数。

DREAMZS算法是由Vrugt 等［20］提出的一种多马尔科夫链的并行采样方法，它能自适应地调整建议分布的步长和方向。使用DREAMZS算法确定权重系数的计算步骤如下：

（1）初始化t= 0。

（2）在［0，1］内随机生成初始权重wj以及N 个初始样本（即N 条马尔科夫链的初始值）。

（3）生成候选样本wj，t+1。在DREAMZS算法中，候选样本是基于差分进化和snooker 更新生成的，分别如下式所示，

式中：wj，t为第j 条马尔科夫链第t 代样本；In为n 阶单位矩阵；e 服从均匀分布Un（－b，b），ε服从正态分布Nn（0，b′），b 和b′为自定义的极小值；δ为用于产生候选样本的平行链对数；γ（δ，n）为比例因子，一般定义为γ=2.38/（2δndef）0.5， wr1（m）和wr2（k）为之前的样本，r1（m）≠r2（k），（m=1，2，…，δ；k=1，2，…，δ）；γs为随机变量，w、wR1和wR2为三个之前的样本， D 为wR1和wR2的函数，D 和γs通常用于确定跳跃的距离。在本研究中，马尔科夫链的更新混合了90%的差分进化和10%的snooker 更新。

（4）基于交叉概率pc确定是否接受候选样本i=1，2，…，n），

式中u 为根据均匀分布U［0，1］产生的随机数。

（5）计算新候选样本的后验概率密度值和接受概率α（wj，t， wj，t+1），

式中：p( wj，t+1| fe)为新候选样本的后验概率密度值；fe为集合代理模型的预测响应值。

根据接受概率α( wj，t，wj，t+1)判断是否接受新候选样本wj，t+1。如果α( wj，t，wj，t+1)≥u，则接受新候选样本；否则不接受，并令wj，t+1=wj，t。

（6）采用Inter Quartile Range（IQR）方法［30］移除无用链。

（7）根据定量收敛判断指标—比例得分因子SR判断采样序列的收敛性［31］。

当SR＜1.2 时，表示算法收敛至稳定的后验分布，可停止迭代计算；否则重复步骤（3）—（7），继续进化平行链。当采样序列收敛后，对稳定的后验样本进行统计分析可获得模型权重系数的随机分布函数。

4 案例分析

以我国西南的L 水电站为研究对象，L 水电站为大渡河干流水电梯级开发的第12 级电站。水电站正常蓄水位1378.00 m，总库容2.4 亿m3，装机容量920 MW。拦河大坝为黏土心墙堆石坝，坝顶高程1385.50 m，坝顶宽12.0 m，长526.70 m，最大坝高79.5 m。

坝址区河床覆盖层深厚，一般为120.0～130.0 m，最大厚度148.6 m，按其物质组成、结构特征、成因和分布情况等自下而上分为4 大层7 个亚层，典型剖面地质图如图2 所示。第①层漂（块）卵（碎）砾石层，系冰水堆积（，分布于坝址区河谷底部。第②层系冰缘泥石流、冲积混合堆积（），主要分布于河床中下部及上坝址右岸谷坡，根据其物质组成及结构特征，可分为②-1 漂（块）卵（碎）砾石层、②-2 碎（卵）砾石土层和②-3 粉细砂及粉土层。第③层系冲、洪积堆积（Q4al+pl），分布于坝址区Ⅰ级阶地和河谷右岸，按其物质组成可分为③-1 含漂（块）卵（碎）砾石土层和③-2 砾石砂层。第④层为冲积（Q4al）堆积之漂卵砾石层，分布于坝址区现代河床及漫滩。

图2 L 水电站典型剖面地质图

坝基土体局部有架空现象，且局部有砂层或粉土层透镜体分布，因粗粒含量较高，土体具中等～强透水性，存在沿坝基覆盖层向下游渗漏的可能性。两岸坝肩强卸荷岩体结构松弛，具强透水性，弱卸荷岩体中等透水，故存在沿坝肩岩体向下游绕渗的问题。

由于③-2 层分布范围较小，故将其与③-1 层合并为一层进行研究。为了便于表示，本文分别以K1、K2、K3、K4、K5、K6表示①层、②-1 层、②-2 层、②-3 层、③层和④层的渗透系数。本文在2018年8月9日的上下游水位条件下（上游水位1368.92 m，下游水位1306.01 m），根据坝后8 个长观孔水头观测点（本文以M1、M2、M3、M4、M5、M6、M7、M8表示）在2018年8月9日的实测数据（分别为1327.17、1340.17、1316.10、1341.10、1329.20、1340.60、1335.08、1334.86 m），基于所提出的贝叶斯框架下大坝渗流参数反演组合代理模型对L 水电站坝基河床覆盖层的渗透系数开展反演研究。

4.1 组合代理模型构建 根据现场勘探数据和室内试验数据，获得覆盖层岩体的渗透系数的取值范围，如表1 所示。假设6 个覆盖层的渗透系数服从均匀分布，采用LHS 方法抽取100 个样本点，将100 个样本点依次输入到耦合VOF 法的三维非稳态水气两相流渗流模型［25］中进行模拟计算，得到100组长观孔水头模拟值，从而生成构建组合代理模型的数据集。随机选取80 组样本作为训练集训练SVR、Kriging 和MARS 这三种代理模型，剩下的20 组样本作为测试集测试代理模型的预测性能。采用DREAMZS算法计算组合代理模型的权重系数，设置3 条平行链，迭代次数设为15 000 次，对迭代稳定后的6000 次样本进行统计分析，计算结果如图3 和表2 所示。

表1 各覆盖层的渗透系数取值范围 （单位：cm/s）

图3 中对角线上的直方图为权重系数的后验分布图，每个直方图的峰值都比较明显，表示权重系数能被较好的识别；下三角上的图是3 个权重系数的二维联合分布图，红色数字标识表示相关系数，上三角上的图是3 个权重系数的二维散点图，分析可知W1和W2，W1和W3间的负相关性较高，相关系数分别为-0.46 和-0.82，而W2和W3间的相关性较小，相关系数仅为-0.13。表2 展示了各测点组合代理模型的权重系数分布函数的最大后验概率估计值、均值和标准差，可以明显看出，对于每个测点，SVR 的权重均值都比Kriging 和MARS 的大，说明SVR 的预测性能最佳。同时，每个测点上SVR的权重标准差均比Kriging 和MARS 的稍大，说明SVR 权重系数的不确定性也比Kriging 和MARS 的稍大。取权重系数分布函数的最大后验概率估计值为权重系数的最终取值，通过加权求和，从而实现组合代理模型的构建。

表2 各测点上组合代理模型的权重系数计算结果

图3 组合代理模型的权重系数计算结果（以测点M1为例，W1、W2和W3分别表示SVR、Kriging 和MARS 的权重系数）

4.2 贝叶斯反演分析 将4.1 节中建立的组合代理模型耦合到贝叶斯参数反演方法中，然后对6 个覆盖层的渗透系数进行反演分析。待反演渗透系数的先验分布设置为均匀分布，其取值范围如表1 所示。然后，采用DREAMZS算法评估待反演渗透系数的后验概率密度函数。其中，DREAMZS算法的平行链设置为4 条，迭代次数设为10 000 次。前6000 次迭代作为“燃烧期”，取后4000 次迭代计算的样本评估待反演渗透系数的后验概率密度函数。

6 个覆盖层的渗透系数反演结果如图4 所示。可以看出K1、K3、K4的后验分布直方图具有明显的峰值，表示长观孔水头监测数据的支持显著减少了这3 个渗透系数的不确定性。而K2、K5、K6的后验分布直方图比较扁平且峰值不明显，说明这几个参数的识别性较差，不确定性较大。本文取渗透系数的最大后验估计值为反演结果，如表3 所示。

图4 渗透系数反演计算结果

表3 渗透系数反演结果 （单位：cm/s）

5 讨论

5.1 组合代理模型和单一代理模型的预测性能评估 为了评估各单一代理模型（SVR、Kriging、MARS）和组合代理模型（Ensemble）的预测性能，根据测试集计算了各测点的四种性能指标，即决定系数（R）、均方根误差（RMSE）、平均绝对误差（MAE）和平均绝对百分比误差（MAPE），如图5 所示。正常情况下，当一个模型的R 值大于0.8 时，可认为这个模型是准确的。在图5 中，除了测点M1 中MARS 的R 值小于0.8 外，其余各测点单一模型和组合代理模型的R 指标均高于0.8。故从R 值上看，各模型的预测结果是准确的。通过进一步的比较可以发现，在各测点处，SVR 的预测性能均优于Kriging 和MARS，这与4.1 节中的分析一致。

图5 各模型的四种预测性能指标计算结果

除M1和M8，其余6 个测点的组合代理模型的R 值均比最佳的SVR 高，而组合代理模型的MAE 值和MAPE 值均比最佳的SVR 低，且组合代理模型的RMSE 值均比各单一代理模型的RMSE 值低。在M8处，组合代理模型的RMSE 值几乎与SVR 的RMSE 值相同。在M1处，组合代理模型的各项性能指标仅次于SVR，而优于Kriging 和MARS。综上所述，相比于单个代理模型，组合代理模型具有更优越的预测性能。

5.2 大坝渗流参数反演结果讨论 在普通计算机上执行一次渗流模拟模型至少需要耗费4 小时，而执行一次组合代理模型计算只需要几秒钟，这极大地提高了贝叶斯反演的计算效率。如果忽略掉运行组合代理模型和贝叶斯反演的时间，时间主要耗费在生成组合代理模型的100 个样本上。本文设置DREAMZS算法的迭代次数为10 000 次，需要调用10 000 次组合代理模型，则基于组合代理模型的贝叶斯反演方法相比于传统的贝叶斯方法将提高10 000/100=100 倍的计算效率。对于更加复杂耗时的渗流模拟，本文所提的组合代理模型提高计算效率的效果更明显。

为了进一步验证所提方法的有效性，将基于单一代理模型（SVR、Kriging 和MARS）和组合代理模型（Ensemble）的渗透系数反演值分别输入到渗流数值模型中进行模拟计算，并对长观孔测点的水头模拟值和实测值进行比较分析。

图6 展示了各监测点上基于组合代理模型和单一代理模型的反演水头与实测水头的对比情况。可以看出，各代理模型反演水头的变化趋势与实测值基本一致。对于测点M1、M2、M3、M6和M7，基于组合代理模型反演结果的绝对偏差均小于基于单一代理模型反演结果的绝对偏差。而对于M4、M5和M8，基于组合代理模型反演结果的绝对偏差均接近于基于最佳单一代理模型反演结果的绝对偏差。进一步地，可计算出8 个测点基于组合代理模型以及基于单一代理模型（SVR、Kriging 和MARS）的反演水头与实测水头的平均绝对误差，分别为5.08 m、5.89 m、6.30 m 和5.79 m。相比于单一代理模型SVR、Kriging 和MARS，基于本文所提组合代理模型的反演结果更准确，精度分别提高了13.78%、19.34%和12.27%。

图6 基于组合代理模型和单一代理模型的反演水头与实测水头的对比

6 结论

渗流参数反演是大坝渗流性态分析的重要环节。针对现有的贝叶斯反演方法大多采用单一机器学习算法来替代计算耗时的渗流正演模型，精度较低的问题，本文提出一种贝叶斯框架下大坝渗流参数反演组合代理模型，为大坝渗流参数反演研究提供一种新思路。主要取得以下成果：

（1）以SVR、Kriging 和MARS 三种机器学习算法构建组合代理模型，综合各模型的预测优势，提高了模型的整体预测精度，解决了单一代理模型精度较低的问题。

（2）提出基于DREAMZS算法的模型权重系数计算方法，利用DREAMZS算法并行采样的优势计算模型权重系数的随机分布函数，在考虑不确定性的条件下得到了更加可靠的权重系数。

（3）案例分析表明，本文所提的组合代理模型相比于SVR、Kriging 和MARS 具有更优越的预测性能，且在显著提高贝叶斯反演计算效率的同时，基于组合代理模型的贝叶斯反演方法相比于基于SVR、Kriging 和MARS 的贝叶斯反演方法能够获得更准确的反演结果，其平均精度分别提高了13.78%、19.34%和12.27%。