一种堆石坝沉降序贯预测方法与流程

本发明专利涉及水利工程中堆石坝的沉降预测，特别涉及一种基于贝叶斯原理和监测数据的堆石坝沉降序贯预测方法。

背景技术：

在过去的几十年中，我国水电岩土工程快速发展，完成了一大批200m级乃至300m级的高堆石坝的施工建设，而堆石坝的动态施工和长期运行安全性面临重大挑战。堆石坝是众多坝型中最有竞争力的坝型(郦能惠等,2012)，但面板堆石坝坝体为散粒料填筑体，具有变形大且持续变形时间长的特点(钮新强等,2016)。堆石坝的行为响应(如沉降和孔压)决定着坝体施工和长期运行的安全性，是大坝安全监测系统中重要的组成部分。例如，堆石坝后期沉降变形可能直接影响到防渗混凝土面板与垫层料之间接触效果，甚至造成混凝土面板的脱空问题，影响到周边缝的防渗效果(李守巨等,2013)。堆石坝的沉降等是决定堆石坝安全性和经济性的关键，准确预测堆石坝的沉降对于指导施工、保证施工质量、防渗效果及坝体施工运行全周期安全性等具有重要意义。

目前在国内外学术研究和工程实践中，预测堆石坝沉降的方法主要有工程类比法、基于试验数据和计算模型的直接预测法、基于人工智能的反馈分析方法等。工程类比方法是确定堆石料参数并预测堆石坝沉降的简便方法，通过对比类似工程的堆石料参数以及相应的沉降值来进行，该方法直观简便计算量小；但该方法受经验性主导，估计的沉降量结果较为粗略，难以反映坝体分层分区填筑过程中的复杂机理，不能预测不同时段下的沉降，受到实际工程具体条件的影响导致预测准确性低。

基于试验数据和计算模型的直接预测法在实际工程中应用较为广泛，它通常的做法是基于原位和室内试验等数据，得到堆石料本构模型参数的取值，将其代入到堆石坝物理计算模型，直接预测堆石坝在不同时期的沉降。需注意的是，决定堆石坝预测值的关键在于堆石料自身的物理力学参数特征，以及能反映坝体沉降变形机理的物理模型。受堆石材料自身性质的离散性和变异性、测量误差、试验条件、缩尺效应等影响，如室内流变试验几个月的周期远远短于坝体的后期持续数年的变形，测定的物理力学参数与实际值存在一定差异，由此计算的堆石坝响应(如应力变形)与实测值差别较大(马洪琪,2011)。这意味着，在实际中仅依赖于勘测或试验数据的预测值常常偏离监测量。

基于人工智能算法(hashash等,2003；huang等,2015；li等,2016)，结合监测(peck,1969；2001)预警理念，反馈分析法可用于预测堆石坝的沉降等。人工智能算法在文献中多有报道，例如，马刚等(2012)采用基于粒子迁徙的粒子群算法和径向基函数神经网络的方法，结合实测位移资料，反演水布垭面板堆石坝静力本构模型参数和流变模型参数；倪沙沙和迟世春(2017)采用粒子群算法优化的径向基神经网络反演糯扎渡大坝的渗透系数；马洪琪(2011)总结了天生桥一级、洪家渡、三板溪、水布垭和巴贡5座面板堆石坝的反演结果，结果表明200m级面板堆石坝的计算结果与监测结果只在一定程度上吻合，尤其是水平位移的计算结果存在明显的偏差。基于人工智能的方法解单一，难以同时融合勘测和监测数据信息，无法考虑监测误差和堆石料本身固有的变异性。

实际中的堆石坝施工运行，及相关沉降预测问题至少有两大特点：首先，堆石坝是分层分仓填筑的，即堆石坝的施工是分步进行的，后期的监测也是随时间推进的(ladd,1991)，堆石坝预测方法最好可以随后期数据不断的加入逐渐更新堆石料参数，并动态预测堆石坝的沉降等。此外，在堆石坝的施工运行中，会监测到大量数据，如何充分融合这些监测数据信息，达到准确预测堆石坝的沉降也是一个实际难题。

技术实现要素：

本发明的目的就是针对现有技术的缺陷，提供一种堆石坝序贯预测方法，按照堆石坝施工和运行特征，分步分时段更新堆石料参数统计特征，实现了堆石坝沉降的序贯动态预测。

本发明提供了一种堆石坝沉降序贯预测方法，其特征在于包括以下步骤：

步骤1：确定堆石坝材料参数的先验分布；

步骤2：收集第i(i＝1,2，...，n)时段内的堆石坝监测数据，确定监测误差和模型误差的统计特征，构建第i时段的似然函数；

步骤3：联合第i时段的先验分布和步骤2的似然函数，构建第i时段的贝叶斯方程；

步骤4：利用马尔可夫链蒙特卡洛模拟法求解第i时段的贝叶斯方程获得第i时段的堆石料参数后验分布；

步骤5：基于第i步更新后的堆石料参数后验分布预测堆石坝的沉降；

步骤6：将第i时段得到的堆石料参数后验分布作为下一时段的先验分布；i＝i+1，并重复步骤2-5；如果i＝n，停止对堆石料参数的更新及对堆石坝沉降的预测。

上述技术方案中，步骤1包括以下步骤：

(a)收集勘测数据，包括原位和室内试验数据、堆石料参数相关案例文献资料，结合类似堆石坝工程参数取值，确定堆石坝参数的可能范围；

(b)采用统计和检验方法，分析堆石料参数的均值、标准差、分位值，并识别堆石料参数的概率分布类型；

(c)基于堆石坝参数的可能范围、均值、标准差、分位值、分布类型，确定堆石料参数的先验概率密度。

上述技术方案中，步骤2根据下式构建第i时段的似然函数：

式中，堆石料参数为ф，其代入到堆石坝计算模型中，即得到沉降计算量f(φ)。di＝[di,1,di,2,...,di,m]^t为在第i层中收集到的q个监测数据，di,j为第i层更新中第j(j＝1,...,q)个监测数据。是主对角线元素为的对角矩阵，其中的元素σi,j＝ci,err×di,j，ci,err为监测数据di,j的监测误差的变异系数；是矩阵ri的行列式。第1步更新中用到的监测数据记为dm1，第2步中的监测数据为dm2，第i-1步中的监测数据为dmi-1。

上述技术方案中，步骤3中基于下式构建第i时段的贝叶斯方程：

p(φ|dm1,dm2,...,dmi)＝ki…k2k1

p(dmi|φ，dm1,dm2,...,dmi-1)…p(dm2|φ，dm1)p(dm1|φ)p0(φ)

式中，p0(φ)为堆石料参数的先验分布，p(dm1|φ)为第1阶段的似然函数，p(dm2|φ)为第2阶段的似然函数，p(dmi|φ，dm1,dm2,...,dmi-1)为第i阶段的似然函数，ki…k2k1为常数。其中，从第1至第i(i≥2)层更新，融入的监测数据量逐渐增大；在经历了第i(i＝1,2,...,n)层更新后，得到的堆石料参数融合了先验和第1至i阶段的监测数据信息，将更新后的参数代入计算堆石坝沉降的确定性物理或统计模型中，即得到第i层时的预测量。

上述技术方案中，步骤4中利用mcmc方法进行模拟抽样，待堆石料参数的马尔可夫链收敛且后验分布符合目标分布后，抽取马尔科夫链后25％的样本，即为融合了先验信息和监测数据信息后的堆石料参数样本。

上述技术方案中，步骤5包括将堆石料参数样本代入当前步的堆石坝计算模型中，得到堆石坝沉降的所有可能值及统计信息，包括预测的沉降均值、标准差、分位值。

本发明按照堆石坝施工特征和运行全过程，将整个周期分成若干时间段；每一时段中，基于贝叶斯原理融合勘测和监测数据信息，得到更新后的堆石料参数统计特征；将前一时段更新后的后验信息作为当前时段施工运行中堆石料的先验信息，进一步融合当前时段新收集到的监测数据，实现对堆石料参数的再次更新；按照不同时段，逐步更新参数并进行沉降预测；在每一时段中，将更新的参数代入到堆石坝物理计算模型，从而准确预测每一时段的后续沉降量。该方法由于符合堆石坝工程施工和运行上的时间推进特征，逐步更新堆石料参数，逐步预测堆石坝的沉降等，这里称之为堆石坝的序贯更新方法。本发明在堆石坝的施工运行的全周期中，融合大量监测数据及勘测数据信息；按照堆石坝施工和运行特征，分步分时段更新堆石料参数统计特征，实现了堆石坝沉降的序贯动态预测。对于动态指导堆石坝的分层施工和长期运行，保证堆石坝的长期安全性等具有指导意义。

附图说明

图1是基于贝叶斯原理和监测数据的堆石坝沉降序贯预测方法的实现步骤；

图2是基于贝叶斯原理和监测数据的堆石坝沉降序贯预测方法框架。

具体实施方式

下面结合附图和具体实施例对本发明作进一步的详细说明，便于清楚地了解本发明，但它们不对本发明构成限定。

如图1所示，本发明提供了一种堆石料参数的序贯更新及堆石坝沉降动态序贯预测的理论框架：

将1,2,...,n时间阶段的监测数据分别记为dm1,dm2...,dmn，通过融合众多监测数据，从而更新堆石料参数的概率模型。由于堆石坝的施工及运行是动态进行的，所以随时间逐步融合不同阶段的数据。

堆石料参数的序贯更新始于先验信息(即prior1)和监测数据dm1的融合。在堆石坝序贯预测的框架下，在第1层更新中，prior1和监测数据dm1的信息分别通过先验分布p0(ф)和似然函数p(dm1|ф)来定量表示。当同时拥有p0(ф)和p(dm1|ф)时，即可基于贝叶斯原理获取堆石料参数在第1层更新中的后验分布：

p(φ|dm1)＝k1p(dm1|φ)p0(φ)(1)

式中，k1＝p(dm1)^-1是第1层更新中被均一化后的常数，其大小与ф无关。p(φ|dm1)反映的是融合prior1和监测数据dm1的信息后堆石料参数的认知水平，即为堆石料参数的后验信息1。

在第2层更新中，堆石料参数的后验信息1作为该层的先验信息(即为prior2)；通过进一步获取的监测数据dm2，实现对prior2的更新。在第2层中，堆石料参数的后验分布p(φ|dm1,dm2)定量地表示为：

p(φ|dm1,dm2)＝k2p(dm2|φ，dm1)p(φ|dm1)(2)

式中，k2＝p(dm2|dm1)^-1是第2层更新中被均一化后的常数。p(dm2|φ，dm1)是在给定数据dm1及参数ф信息时得到dm2的概率；p(φ|dm1)是第2层更新中堆石料参数的先验分布，该分布也是式(1)表示的第1层更新后得到的后验分布。

类似地，在第3层更新中，式(2)中的后验分布p(φ|dm1,dm2)将作为本层中的先验分布。通过进一步融入监测数据dm3的信息，更新先验分布p(φ|dm1,dm2)，得到第3层的后验分布：

p(φ|dm1,dm2,dm3)＝k3p(dm3|φ，dm1,dm2)p(φ|dm1,dm2)(3)

式中，k3＝p(dm3|dm1,dm2)^-1是第3层更新中被均一化后的常数。p(dm3|φ，dm1,dm2)是在给定数据dm3时反映参数ф信息的似然函数；p(φ|dm1,dm2)是第3层更新中堆石料参数的先验分布，该分布也是式(2)表示的第2层更新后得到的后验分布。最终，联合式(1)-(3)，在贝叶斯序贯更新的框架下，得到经过3层更新后得到的堆石料参数的后验分布：

式(4)表示的经过第3层更新得到的堆石料参数的后验分布；在实际的堆石坝施工和长期运行过程中，全周期中有多个时间阶段，阶段总数记为n；在经历了n层贝叶斯更新后得到的堆石料参数的后验分布即为：

值得注意的是，从第1至第n(n≥2)层更新，融入的监测数据量逐渐增大；在经历了第n层更新后，得到的堆石料参数融合了先验和第1至n阶段的监测数据信息，将更新后的参数代入计算堆石坝沉降的确定性物理或统计模型f(φ)中，即得到第n层时的预测量。

在堆石坝分层建筑和长期运行中，随时间逐步收集堆石坝第1至当时第i时间区段(i＝1,2,…,n)内的监测数据(dm1，dm2，…,dmi,…)，融合先验信息，得到堆石料材料参数的后验分布，获得堆石坝在当前时间段之后的沉降预测量。

本发明提供了一种堆石坝沉降序贯预测方法，其特征在于包括以下步骤：

步骤1，确定堆石坝材料参数的先验分布：

(a)收集勘测数据，包括原位和室内试验数据、堆石料参数相关案例文献资料，结合类似堆石坝工程参数取值，确定堆石坝参数的可能范围；

(b)采用统计和检验方法，分析堆石料参数的均值、标准差、分位值，并识别堆石料参数的概率分布类型；

(c)基于堆石坝参数的可能范围、均值、标准差、分位值、分布类型，确定堆石料参数的先验概率密度。

步骤2：收集第i(i＝1,2，...，n)时段内的堆石坝监测数据，确定监测误差和模型误差的统计特征，构建第i时段的似然函数；

步骤3：联合第i时段的先验分布和步骤2的似然函数，构建第i时段的贝叶斯方程；基于下式构建第i时段的贝叶斯方程：

p(φ|dm1,dm2,...,dmi)＝ki…k2k1

p(dmi|φ，dm1,dm2,...,dmi-1)…p(dm2|φ，dm1)p(dm1|φ)p0(φ)

其中，从第1至第i(i≥2)层更新，融入的监测数据量逐渐增大；在经历了第i(i＝1,2,...,n)层更新后，得到的堆石料参数融合了先验和第1至i阶段的监测数据信息，将更新后的参数代入计算堆石坝沉降的确定性物理或统计模型中，即得到第i层的预测量。

步骤4：利用马尔可夫链蒙特卡洛模拟法求解第i时段的贝叶斯方程获得第i时段的后验分布。利用mcmc方法进行模拟抽样，待堆石料参数的马尔可夫链收敛且后验分布符合目标分布后，抽取马尔科夫链后25％的样本，即为融合了先验信息和监测数据信息后的堆石料参数样本。在堆石坝的沉降预测中，由于本构模型中的参数较多，给堆石料参数的更新带来高维问题。在高维问题中ф分量间的相关性较为复杂，采用传统的直接积分法可能会导致马尔可夫链难以收敛或需要耗费大量的计算成本。马尔科夫链蒙特卡洛(mcmc)模拟是解决该问题的最有代表性的方法之一，文献中有诸多方法可以实现该目的，差分演化算法是其中一种，它改进产生建议样本的算法：让产生建议样本的过程遵从参数间的相关性。该方法是使用多链进行mcmc模拟的算法之一，它不是像常规metropolis和metropolis-hastings算法一样通过单链从单个初始值开始模拟，而是从一组初始值开始，并从每个初始值生成一个样本链形成多链。多链允许通过来自其它链样本间的相关性来决定当前链中的建议样本，从而有效解决建议分布不满足目标分布中参数间相关性而导致难以收敛或耗时较长的问题

步骤5：基于第i步更新后的堆石料参数预测堆石坝的沉降。将该样本代入当前步的堆石坝计算模型中，得到堆石坝沉降的所有可能值及统计信息，包括预测的沉降均值、标准差、分位值等。

步骤6：将第i时段得到的堆石料参数后验分布作为下一时段的先验分布；i＝i+1，并重复步骤2-5；如果i＝n，停止对堆石料参数的更新及对堆石坝沉降的预测。

上述技术方案中，步骤2中将每一层更新得到的参数ф代入到堆石坝计算模型中，即可得到沉降计算量f(φ)。但计算模型存在模型误差，任何的监测数据d都存在监测误差，模型计算值f(φ)将与监测值之间存在差异，该差值记为e，e同时包含监测误差和模型误差：

d＝f(φ)+e(6)

假设e＝[e1,e2,…,em]服从均值为零的正态分布，且各元素之间相关独立，m是当前更新层监测数据的数量。

数学上，第1层中的似然函数p(dm1|φ)表示的是e的概率密度分布；实际中表示的是在给定参数ф时观到dm1的概率。第1层更新中的似然函数表示为：

式中，d1＝[d1,1,d1,2,...,d1,m]^t为在第1层中收集到的m个监测数据，d1,j为第1层更新中第j(j＝1,...,m)个监测数据。是主对角线元素为的对角矩阵，其中的元素σ1,j＝c1,err×d1,j，c1,err为监测数据d1,j的监测误差的变异系数；是矩阵r1的行列式。

堆石料参数的第2层更新中，同样基于式(6)来推导似然函数p(dm2|φ，dm1)。类似地，p(dm2|φ，dm1)表达的是e的概率密度分布，表示为：

式中，p是第2层中监测数据d2＝[d2,1,d2,2,...,d2,p]^t的数据量。

类似地，第i(i＝1,2,…,n)层更新中，似然函数p(dmi|φ，dm1,dm2,...,dmi-1)表示为：

式中，q是第i层中监测数据di＝[di,1,di,2,...,di,q]^t的数据量。

将(7)-(9)代入到式(5)中，即得到第i层更新后堆石料参数的后验分布，将该更新后的参数代入了堆石坝沉降计算模型f(φ)中，即可得到第i时间段后的预测值。

本发明所提序贯预测方法符合堆石坝施工和运行特征：堆石坝的施工和运行全周期可以分成不同的时间段，在不同的时间段内可以收集到监测数据；在序贯预测的框架下，可以逐步融入监测数据，逐步增强对堆石料参数特征的认识，并在不同时间段内预测堆石坝在后续时间段内的沉降值，从而为下一步的施工或运行提供定性的参考依据。本发明能够通过逐步融入监测数据，提高预测沉降的准确性；相比于传统的人工智能方法，所提方法能有效考虑监测误差等因素，能够提供在非完备信息下预测的沉降的所有可能性，并通过概率密度分布的形式进行呈现。

本说明书未作详细描述的内容属于本领域专业技术人员公知的现有技术。