CN112307536B - 一种大坝渗流参数反演方法 - Google Patents

Abstract

本发明公开一种大坝渗流参数反演方法，包括如下步骤：步骤1：采用拉丁超立方抽样方法从待反演渗流参数的取值范围内抽取样本点；步骤2：通过渗流数值模拟计算各样本点对应的响应值并构建训练样本和测试样本；步骤3：利用训练样本训练各代理模型；步骤4：基于DREAM算法确定各代理模型的权重系数，通过加权求和构建集合代理模型；步骤5：采用步骤4中所构建的集合代理模型替代计算耗时的大坝渗流数值模拟模型，基于贝叶斯反演方法计算获得待反演渗流参数，该方法准确高效地获得大坝渗流参数，从而为大坝渗流性态分析提供更可靠的支撑和依据。

Description

一种大坝渗流参数反演方法

技术领域

本发明属于水利水电工程中大坝渗流性态分析领域，具体涉及一种大坝渗流参数反演方法。

背景技术

渗流是影响大坝安全稳定的一个重要因素^[1]，大坝渗流性态分析研究为保障大坝渗流安全稳定提供了理论与技术支持。准确获取岩土体介质的渗流参数对大坝渗流性态进行科学分析的重要前提。在实际工程中，通常根据室内试验和原位试验获取渗流参数，然而由于大坝工程的复杂性和地质条件的不确定性，所得到的渗流参数和现实情况均有较大的偏差^[2]。将水头或渗流量等较容易观测的变量作为已知信息，对岩土体介质的渗流参数进行反演计算是一种有效手段。

目前常用的参数反演方法主要分为确定性反演方法和随机反演方法两大类。确定性反演方法发展历史较长，而且通常简单有效，在实际中得到了广泛应用^[3-5]。然而，确定性反演方法得出的反演参数仅仅是一组使模拟值与观测值间差值最小的最优参数组合，难以充分考虑参数反演过程中的不确定性，导致反演结果的可靠性得不到保障^[6]。为了量化参数反演过程中的不确定性，人们提出了随机反演方法。随机反演方法将待反演的参数作为随机变量，利用观测资料更新参数的概率分布，并采用最终的概率分布来描述这些随机变量，进而获得参数的最优解以及不确定性特征^[2]。贝叶斯反演方法是一种广泛使用的随机反演分析方法。当获得现场监测数据和待估参数的先验分布后，参数的后验分布可以通过贝叶斯规则反推出来^[7]。尽管贝叶斯反演方法具有严格的数学基础，但是它需要大量调用计算耗时的正演模型以使后验分布达到收敛，计算效率较低^[8]。

采用代理模型来替代计算耗时的数值模拟模型是提高贝叶斯反演效率的有效方法^[9]。常用的代理模型包括多项式响应面(PR)^[10]、人工神经网络(ANN)^[11]、径向基函数(RBF)^[12]、克里金(Kriging)^[13]、支持向量机(SVR)^[14]以及多元自适应回归样条(MARS)^[15]等。目前大多数研究都采用的是单个代理模型^[16-18]。然而，单个代理模型的预测性能不稳定，在某些样本输入情况下得到的预测结果可能会具有较大误差，将影响参数反演结果的准确性。

发明内容

针对目前渗流参数反演方法大多忽略了反演过程中的不确定性，而贝叶斯反演方法虽然能够考虑不确定性问题，但计算效率较低，且通常采用的单个代理模型的预测性能不稳定的不足，本发明提出一种大坝渗流参数反演方法，该方法准确高效地获得大坝渗流参数，从而为大坝渗流性态分析提供更可靠的支撑和依据。

为克服现有技术的不足，本发明采用的技术方案如下：

一种大坝渗流参数反演方法，包括如下步骤：

步骤1：采用拉丁超立方抽样方法从待反演渗流参数的取值范围内抽取样本点；

步骤2：通过渗流数值模拟计算各样本点对应的响应值并构建训练样本和测试样本；

步骤3：利用训练样本训练各代理模型；

步骤4：基于DREAM算法确定各代理模型的权重系数，通过加权求和构建集合代理模型；

步骤5：采用步骤4中所构建的集合代理模型替代计算耗时的大坝渗流数值模拟模型，基于贝叶斯反演方法计算获得待反演渗流参数。

进一步地，步骤1中，所述的采用拉丁超立方抽样方法从待反演渗流参数的取值范围内抽取样本点包括：

根据工程经验以及室内试验和原位试验数据确定待反演渗流参数的取值范围和服从的分布函数类型，采用拉丁超立方抽样方法从待反演渗流参数的取值范围内抽取样本点。

进一步地，步骤2中，所述的通过渗流数值模拟计算各样本点对应的响应值并构建训练样本和测试样本包括：

将抽取的样本点逐个输入到大坝渗流数值模拟模型中进行模拟计算，得出样本点对应的响应值，将样本点与响应值组成样本对，取设定比例的样本对为训练样本，剩下的样本对为测试样本。

进一步地，步骤3中，所述的利用训练样本训练各代理模型包括：

利用步骤2中获得的训练样本训练各代理模型，并利用测试样本测试各代理模型的预测性能。

进一步地，步骤4中，所述的基于DREAM算法确定各代理模型的权重系数，通过加权求和构建集合代理模型包括：

DREAM算法是近年来提出的一种具有较高求解效率的多条马尔科夫链并行采样的方法，各条平行链在搜索过程中通过引入差分进化遗传算法自适应地调整搜索步长和方向，有效搜索多个全局最优区域，从而提高样本的多样性,目前已被成功用于不确定性分析问题中。采用DREAM算法确定各代理模型的权重系数的计算步骤如下:

在[0,1]范围内随机生成各代理模型的初始权重作为先验信息，随机产生权重系数的N个初始样本(即为N条马尔科夫链的起点)，计算各样本的后验概率值,开始进行N条平行链的进化；

针对第j(j＝1,2,…,N)条马尔科夫链，产生候选样本w_j,t+1,计算如公式(1)所示,

其中,w_j,t为第j条马尔科夫链第t代的样本；I_n为n阶单位矩阵；e、ε为随机产生的很小的数，且ε服从U_n(-b,b)，ε服从N_n(0,b')，b、b'为自定义的极小值；δ代表用于产生候选样本的平行链对数；γ(δ,n)为比例因子，一般定义为

r₁(m)，r₂(k)为随机选取的平行链编号，并满足r₁(m),r₂(k)∈{1,2,…,N}，r₁(m)≠r₂(k),(m＝1,2,…,δ；k＝1,2,…,δ)；

根据交叉概率p_c确定是否接受候选样本分量

定义交叉概率p_c∈[0,1],n_def＝n,如果u≤1-p_c，则不接受，

反之，则接受.u为根据0-1均匀分布产生的随机数；

计算新候选样本w_j,t+1的后验概率值，并计算接受概率α(w_j,t,w_j,t+1),计算如公式(2)所示,

根据α(w_j,t,w_j,t+1)判断是否接受新候选样本w_j,t+1,如果α(w_j,t,w_j,t+1)≥u,则接受新候选样本,否则不接受,并令w_j,t+1＝w_j,t；

根据Inter－Quartile－Range(IQR)方法统计去除无用链；

采用比例得分因子S_R判断采样过程的收敛性,比例得分因子S_R的计算如公式(3)所示，

其中，g为每条马尔科夫链的进化代数，q为用于评价的马尔科夫链条数，B/g为q条马尔科夫链平均值的方差,W为q条马尔科夫链方差的平均值；

当S_R＜1.2时,表示权重系数收敛于稳定的后验分布，则计算结束，取权重系数后验分布的均值为各代理模型的最终权重w_j，否则重复公式(1)到公式(2)之间的步骤，继续进化平行链；

根据获得的各代理模型的权重系数w_j，通过加权求和构建集合代理模型，计算如公式(4)所示，

其中，f_en(x)为集合代理模型的预测值；f_j(x)为第j(j＝1,2,…,k)个代理模型的预测值；w_j为第j(j＝1,2,…,k)个代理模型的权重系数，且满足

进一步地，步骤5中，所述的采用步骤4中所构建的集合代理模型替代计算耗时的大坝渗流数值模拟模型，基于贝叶斯反演方法计算获得待反演渗流参数包括：

采用步骤4中所构建的集合代理模型f_en(x)替代计算耗时的渗流数值模拟模型，根据渗流实测数据，通过贝叶斯反演方法计算待反演渗流参数的后验分布，计算如公式(5)所示，

p(x|y_ob)∝p(y_ob|x)p(x) (5)

其中,p(x)为待反演渗流参数的先验分布函数；y_ob为监测值，监测误差服从均值为零，方差为

的正态分布；p(y_ob|x)为似然函数，用来衡量集合代理模型预测值f_en(x)和监测值y_ob之间的拟合程度，设定集合代理模型预测值f_en(x)和监测值y_ob之间的偏差为δ＝y_ob-f_en(x)，且服从一个协方差矩阵为Σ的多元正态分布，则p(y_ob|x)的计算如公式(6)所示，

其中,

为

的行列式；I为n×n的单位矩阵；

取待反演渗流参数后验分布的最大后验概率估计值为待反演渗流参数的反演值。

有益效果

1、本发明提出的基于DREAM算法集合代理模型的大坝渗流参数反演方法,能够考虑参数反演过程中的不确定性,准确高效地获得大坝渗流参数，从而为大坝渗流性态分析提供更可靠的支撑和依据。

2、采用本发明提出的基DREAM算法的集合代理模型替代计算耗时的渗流数值模型,能够大幅度提高原有贝叶斯反演方法的计算效率，降低计算成本；采用本发明提出的基于DREAM算法的集合代理模型能够解决单一代理模型预测性能不稳定的问题,从而提高参数反演结果的准确性。

3、本发明提出的基于DREAM算法的集合代理模型不仅可以用于参数反演领域，还可以广泛应用于参数优化、参数敏感性分析等领域，可推广性较强。

附图说明

图1为本发明实施例的涉及基于DREAM算法集合代理模型的大坝渗流参数反演方法的流程图；

图2为本发明实施例的大坝渗流结构示意图；

图3为本发明实施例的坝基渗透系数反演结果图；

图4为本发明实施例的单一代理模型(RBF、Kriging、MARS)以及基于DREAM算法集合代理模型的预测性能对比图(以水头测点3为例)。

具体实施方式：

以下结合实施案例和附图，对本发明的技术和方法进行详细描述，以下案例和附图用于说明本发明的构成，但不是用来限定本发明的范围。

为了便于理解本发明，下文将结合说明书附图和实施例对本发明作更详细地描述，但本发明的保护范围并不限于以下具体的实施例。

参照附图1，对本发明的具体步骤作进一步的详细描述:

步骤1：采用拉丁超立方抽样方法从待反演渗流参数的取值范围内抽取样本点；

步骤2：通过渗流数值模拟计算各样本点对应的响应值并构建训练样本和测试样本；

步骤3：利用训练样本训练各代理模型；

步骤4：基于DREAM算法确定各代理模型的权重系数，通过加权求和构建集合代理模型；

步骤5：采用步骤4中所构建的集合代理模型替代计算耗时的大坝渗流数值模拟模型，基于贝叶斯反演方法计算获得待反演渗流参数。

进一步地，步骤1中，所述的采用拉丁超立方抽样方法从待反演渗流参数的取值范围内抽取样本点包括：

根据工程经验以及室内试验和原位试验数据确定待反演渗流参数的取值范围和服从的分布函数类型，采用拉丁超立方抽样方法从待反演渗流参数的取值范围内抽取样本点。

进一步地，步骤2中，所述的通过渗流数值模拟计算各样本点对应的响应值并构建训练样本和测试样本包括：

将抽取的样本点逐个输入到大坝渗流数值模拟模型中进行模拟计算，得出样本点对应的响应值，将样本点与响应值组成样本对，取设定比例的样本对为训练样本，剩下的样本对为测试样本。

进一步地，步骤3中，所述的利用训练样本训练各代理模型包括：

利用步骤2中获得的训练样本训练各代理模型，并利用测试样本测试各代理模型的预测性能。

进一步地，步骤4中，所述的基于DREAM算法确定各代理模型的权重系数，通过加权求和构建集合代理模型包括：

DREAM算法是近年来提出的一种具有较高求解效率的多条马尔科夫链并行采样的方法，各条平行链在搜索过程中通过引入差分进化遗传算法自适应地调整搜索步长和方向，有效搜索多个全局最优区域，从而提高样本的多样性,目前已被成功用于不确定性分析问题中。采用DREAM算法确定各代理模型的权重系数的计算步骤如下:

在[0,1]范围内随机生成各代理模型的初始权重作为先验信息，随机产生权重系数的N个初始样本(即为N条马尔科夫链的起点)，计算各样本的后验概率值,开始进行N条平行链的进化；

针对第j(j＝1,2,…,N)条马尔科夫链，产生候选样本w_j,t+1,计算如公式(1)所示,

其中,w_j,t为第j条马尔科夫链第t代的样本；I_n为n阶单位矩阵；e、ε为随机产生的很小的数，且ε服从U_n(-b,b)，ε服从N_n(0,b')，b、b'为自定义的极小值；δ代表用于产生候选样本的平行链对数；γ(δ,n)为比例因子，一般定义为

r₁(m)，r₂(k)为随机选取的平行链编号，并满足r₁(m),r₂(k)∈{1,2,…,N}，r₁(m)≠r₂(k),(m＝1,2,…,δ；k＝1,2,…,δ)；

根据交叉概率p_c确定是否接受候选样本分量

定义交叉概率p_c∈[0,1],n_def＝n,如果u≤1-p_c，则不接受，

反之，则接受.u为根据0-1均匀分布产生的随机数；

计算新候选样本w_j,t+1的后验概率值，并计算接受概率α(w_j,t,w_j,t+1),计算如公式(2)所示,

根据α(w_j,t,w_j,t+1)判断是否接受新候选样本w_j,t+1,如果α(w_j,t,w_j,t+1)≥u,则接受新候选样本,否则不接受,并令w_j,t+1＝w_j,t；

根据Inter－Quartile－Range(IQR)方法统计去除无用链；

采用比例得分因子S_R判断采样过程的收敛性,比例得分因子S_R的计算如公式(3)所示，

其中，g为每条马尔科夫链的进化代数，q为用于评价的马尔科夫链条数，B/g为q条马尔科夫链平均值的方差,W为q条马尔科夫链方差的平均值；

当S_R＜1.2时,表示权重系数收敛于稳定的后验分布，则计算结束，取权重系数后验分布的均值为各代理模型的最终权重w_j，否则重复公式(1)到公式(2)之间的步骤，继续进化平行链；

根据获得的各代理模型的权重系数w_j，通过加权求和构建集合代理模型，计算如公式(4)所示，

其中，f_en(x)为集合代理模型的预测值；f_j(x)为第j(j＝1,2,…,k)个代理模型的预测值；w_j为第j(j＝1,2,…,k)个代理模型的权重系数，且满足

进一步地，步骤5中，所述的采用步骤4中所构建的集合代理模型替代计算耗时的大坝渗流数值模拟模型，基于贝叶斯反演方法计算获得待反演渗流参数包括：

采用步骤4中所构建的集合代理模型f_en(x)替代计算耗时的渗流数值模拟模型，根据渗流实测数据，通过贝叶斯反演方法计算待反演渗流参数的后验分布，计算如公式(5)所示，

p(x|y_ob)∝p(y_ob|x)p(x) (5)

其中,p(x)为待反演渗流参数的先验分布函数；y_ob为监测值，监测误差服从均值为零，方差为

的正态分布；p(y_ob|x)为似然函数，用来衡量集合代理模型预测值f_en(x)和监测值y_ob之间的拟合程度，设定集合代理模型预测值f_en(x)和监测值y_ob之间的偏差为δ＝y_ob-f_en(x)，且服从一个协方差矩阵为Σ的多元正态分布，则p(y_ob|x)的计算如公式(6)所示，

其中,

为

的行列式；I为n×n的单位矩阵；

取待反演渗流参数后验分布的最大后验概率估计值为待反演渗流参数的反演值。

下面是一个采用上述方法的实施例。

参照附图2,由于本发明实施例中以坝基渗透系数反演来证明本发明提出方法的有效性和优越性，为简化计算，将土石坝坝体和防渗墙均设为不透水体。坝基包括地层A、地层B、地层C三个地层，各地层均设为各向同性的均匀多孔介质。各地层的渗透系数设为待反演渗流参数，其取值范围分别设为K_A∈[0,2]×10^-3m/s、K_B∈[0.2,1.8]×10^-4m/s、K_C∈[0.5,1.5]×10^-5m/s，且其先验分布均设为均匀分布。设地层A、地层B、地层C的渗透系数的“真值”分别为1×10^-3m/s、1×10^-4m/s、1×10^-5m/s，上游水位H₁＝20m，下游水位H₂＝0m，通过渗流数值模拟，得到3个水头监测点的模拟值以及坝基下游渗流量的模拟值，并将其作为相应的“实际监测值”。

采用拉丁超立方抽样方法从这三个待反演渗透系数的取值范围内抽取50组渗透系数样本点。将这50组渗透系数样本点输入到大坝渗流数值模拟模型中进行模拟计算，得出各组渗透系数样本点对应的各监测点响应值，将渗透系数样本点与监测点响应值组成50组样本对。取50组样本对中的40组样本对作为训练样本，分别训练各个代理模型，将剩余的10组样本对作为检验代理模型预测性能的测试样本。

以40组训练样本训练径向基函数(RBF)、克里金(Kriging)和多元自适应回归样条(MARS)这三种代理模型，基于DREAM算法求得单个代理模型的权重系数分别为：水头测点1(w_RBF＝0.4650,w_Kriging＝0.4265,w_MARS＝0.1081)、水头测点2(w_RBF＝0.4910,w_Kriging＝0.4392,w_MARS＝0.0689)、水头测点3(w_RBF＝0.5761,w_Kriging＝0.3316,w_MARS＝0.0913)、坝基下游渗流量(w_RBF＝0.4918,w_Kriging＝0.5062,w_MARS＝0.0001)。各单个代理模型通过加权求和构建最终的集合代理模型f_en(x)。

将基于DREAM算法的集合代理模型f_en(x)替代计算耗时的渗流数值模型，通过贝叶斯反演方法获得待反演渗透系数的后验分布，如附图3所示，取该后验分布的最大后验概率估计值为待反演渗透系数的取值，即地层A、地层B、地层C的渗透系数反演值分别为1.04x10^-3m/s、0.99x10^-4m/s、1.17x10^-5m/s。反演值很接近“真值”，故本发明所提方法是有效的。

为了进一步验证本发明提出的基于DREAM算法的集合代理模型的优越性，基于10组测试样本数据，采用相关系数R²、均方根误差RMSE、平均绝对偏差MAE和中位数绝对偏差MAD这四种统计指标，对比分析三种单一代理模型(RBF、Kriging和MARS)和基于DREAM算法的集合代理模型的预测性能，如附图4所示。在附图4中可以看出，相比于其他几种代理模型，本发明提出的基于DREAM算法的集合代理模型具有较高的相关系数R²，而均方根误差RMSE、平均绝对偏差MAE和中位数绝对偏差MAD较低，故可说明本发明所提出的基于DREAM算法的集合代理模型具有较优越的预测性能。

本发明并不限于上文描述的实施方式。以上对具体实施方式的描述旨在描述和说明本发明的技术方案，上述的具体实施方式仅仅是示意性的，并不是限制性的。在不脱离本发明宗旨和权利要求所保护的范围情况下，本领域的普通技术人员在本发明的启示下还可做出很多形式的具体变换，这些均属于本发明的保护范围之内。

Claims (4)

1.一种大坝渗流参数反演方法，其特征在于：包括如下步骤：

步骤1：采用拉丁超立方抽样方法从待反演渗流参数的取值范围内抽取样本点；

步骤2：通过渗流数值模拟计算各样本点对应的响应值并构建训练样本和测试样本；

步骤3：利用训练样本训练各代理模型；

步骤4：基于DREAM算法确定各代理模型的权重系数，通过加权求和构建集合代理模型；

步骤5：采用步骤4中所构建的集合代理模型替代计算耗时的大坝渗流数值模拟模型，基于贝叶斯反演方法计算获得待反演渗流参数； 其中：

所述步骤4包括如下步骤：

在[0,1]范围内随机生成各代理模型的初始权重作为先验信息，随机产生权重系数的N个初始样本，计算各样本的后验概率值,开始进行N条平行链的进化；

针对第j(j＝1,2,…,N)条马尔科夫链，产生候选样本w_j,t+1,计算如公式(1)所示,

其中,w_j,t为第j条马尔科夫链第t代的样本；I_n为n阶单位矩阵；e、ε为随机产生的很小的数，且ε服从U_n(-b,b)，ε服从N_n(0,b')，b、b'为自定义的极小值；δ代表用于产生候选样本的平行链对数；γ(δ,n)为比例因子，一般定义为

r₁(m)，r₂(k)为随机选取的平行链编号，并满足r₁(m),r₂(k)∈{1,2,…,N}，r₁(m)≠r₂(k),(m＝1,2,…,δ；k＝1,2,…,δ)；

根据交叉概率p_c确定是否接受候选样本分量

定义交叉概率p_c∈[0,1],n_def＝n,如果u≤1-p_c，则不接受，

反之，则接受； u为根据0-1均匀分布产生的随机数；

计算新候选样本w_j,t+1的后验概率值，并计算接受概率α(w_j,t,w_j,t+1),计算如公式(2)所示,

根据α(w_j,t,w_j,t+1)判断是否接受新候选样本w_j,t+1,如果α(w_j,t,w_j,t+1)≥u,则接受新候选样本,否则不接受,并令w_j,t+1＝w_j,t；

根据Inter－Quartile－Range(IQR)方法统计去除无用链；

采用比例得分因子S_R判断采样过程的收敛性,比例得分因子S_R的计算如公式(3)所示，

其中，g为每条马尔科夫链的进化代数，q为用于评价的马尔科夫链条数，B/g为q条马尔科夫链平均值的方差,W为q条马尔科夫链方差的平均值；

当S_R<1.2时,表示权重系数收敛于稳定的后验分布，则计算结束，取权重系数后验分布的均值为各代理模型的最终权重w_j，否则重复公式(1)到公式(2)之间的步骤，继续进化平行链；

根据获得的各代理模型的权重系数w_j，通过加权求和构建集合代理模型，计算如公式(4)所示，

其中，f_en(x)为集合代理模型的预测值；f_j(x)为第j(j＝1,2,…,k)个代理模型的预测值；w_j为第j(j＝1,2,…,k)个代理模型的权重系数，且满足

其中：

所述步骤5包括如下步骤：

采用步骤4中所构建的集合代理模型替代计算耗时的渗流数值模拟模型，根据渗流实测数据，通过贝叶斯反演方法计算待反演渗流参数的后验分布，计算如公式(5)所示，

p(x|y_ob)∝p(y_ob|x)p(x) (5)

其中,p(x)为待反演渗流参数的先验分布函数；y_ob为监测值，监测误差服从均值为零，方差为

的正态分布；p(y_ob|x)为似然函数，用来衡量集合代理模型预测值f_en(x)和监测值y_ob之间的拟合程度，设定集合代理模型预测值f_en(x)和监测值y_ob之间的偏差为δ＝y_ob-f_en(x)，且服从一个协方差矩阵为Σ的多元正态分布，则p(y_ob|x)的计算如公式(6)所示，

其中,

为

的行列式；I为n×n的单位矩阵；

取待反演渗流参数后验分布的最大后验概率估计值为待反演渗流参数的反演值。

2.根据权利要求1所述的一种大坝渗流参数反演方法，其特征在于：所述步骤1包括如下步骤：

根据工程经验以及室内试验和原位试验数据确定待反演渗流参数的取值范围和服从的分布函数类型，采用拉丁超立方抽样方法从待反演渗流参数的取值范围内抽取样本点。

3.根据权利要求1所述的一种大坝渗流参数反演方法，其特征在于：所述步骤2包括如下步骤：

将抽取的样本点逐个输入到大坝渗流数值模拟模型中进行模拟计算，得出样本点对应的响应值，将样本点与响应值组成样本对，取设定比例的样本对为训练样本，剩下的样本对为测试样本。

4.根据权利要求1所述的一种大坝渗流参数反演方法，其特征在于：所述步骤3包括如下步骤：

利用步骤2中获得的训练样本训练各代理模型，并利用测试样本测试各代理模型的预测性能。

Images