CN111914475B - 一种加速刻画高斯水文地质参数场的贝叶斯逆模拟方法 - Google Patents

一种加速刻画高斯水文地质参数场的贝叶斯逆模拟方法 Download PDF

Info

Publication number
CN111914475B
CN111914475B CN202010602434.4A CN202010602434A CN111914475B CN 111914475 B CN111914475 B CN 111914475B CN 202010602434 A CN202010602434 A CN 202010602434A CN 111914475 B CN111914475 B CN 111914475B
Authority
CN
China
Prior art keywords
gaussian
sample
markov
distribution
chains
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202010602434.4A
Other languages
English (en)
Other versions
CN111914475A (zh
Inventor
徐腾
张文俊
陈沁宇
鲁春辉
谢一凡
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hohai University HHU
Original Assignee
Hohai University HHU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hohai University HHU filed Critical Hohai University HHU
Priority to CN202010602434.4A priority Critical patent/CN111914475B/zh
Publication of CN111914475A publication Critical patent/CN111914475A/zh
Application granted granted Critical
Publication of CN111914475B publication Critical patent/CN111914475B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/27—Design optimisation, verification or simulation using machine learning, e.g. artificial intelligence, neural networks, support vector machines [SVM] or training a model
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N7/00—Computing arrangements based on specific mathematical models
    • G06N7/01—Probabilistic graphical models, e.g. probabilistic networks

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • Software Systems (AREA)
  • Artificial Intelligence (AREA)
  • General Engineering & Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Algebra (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computing Systems (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Physics (AREA)
  • Mathematical Analysis (AREA)
  • Probability & Statistics with Applications (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Medical Informatics (AREA)
  • Computer Hardware Design (AREA)
  • Geometry (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种加速刻画高斯水文地质参数场的贝叶斯逆模拟方法，包括以下步骤：设立初始水文参数场，建立温度梯度和跳跃因子梯度；通过预条件克兰克尼科尔森建议分布得到当前建议采样；对于每条马尔科夫链，当满足接受概率时，下一个采样实现为当前建议采样，否者维持不变；对于任何一对不同温度的马尔科夫链，当样本个数和建议交换频率的商为整数时，且满足交换接受概率，冷链和热链之间的信息进行交换，否者维持不变；重复以上过程进行下一个高斯场样本采样，根据设定的马尔科夫链的长度，判断是否终止取样。本方法是一个渐进精确的采样方法，在处理地质统计学逆模拟问题中，比预条件克兰克尼科尔森马尔科夫蒙特卡洛方法效率更高。

Description

一种加速刻画高斯水文地质参数场的贝叶斯逆模拟方法
技术领域
本发明属于水文统计学技术领域，具体涉及一种加速刻画高斯水文地质参数场的贝叶斯逆模拟方法。
背景技术
地下水状态跟人类社会的可持续发展息息相关。一个值得信赖的地下水水流和溶质运移模拟在地下水水位预报，溶质运移预测，地下水资源管理以及地下水环境危险性评估方面是具有非常重要的意义。地下水水流和溶质运移模拟的质量好坏在很大程度上取决于地下水水文参数刻画的好坏(例如，水文渗透系数，孔隙度，等等)。然而，相对于含水层的尺度而言，在实际野外得到的这些水文参数的试验测量数据是非常稀少地。而稀少的水文地质参数测量数据是很难被用来直接刻画研究区域水文地质参数场。所以，如何利用这些稀少的野外水文地质参数来尽可能准确的刻画整个含水层的水文地质参数场是一个非常大的挑战。特别是对于大尺度的研究区而言，在最近几十年来，同化实时状态数据的随机逆模拟已经被证明是一个非常有用的能够刻画水文地质参数场的工具。
并且，综合以上研究还可发现，当今主流的反演非线性高斯参数场的逆模拟方法主要还是集中在集合卡尔曼滤波、集合平滑法、逆序贯模拟法或基于马尔科夫蒙特卡洛的贝叶斯方法当中。其中，尽管基于MCMC的方法计算相当耗时并且采样效率低下，但它是经典的可以处理复杂参数分布及非线性问题的逆模拟方法。反观其他方法，他们都依赖于某种线性化或一阶、二阶矩近似，虽然这些方法通常比马尔科夫蒙特卡洛计算快得多，但随着非线性问题的增加，反演的结果会变得越来越不精确并且会造成不稳定的近似。因此，对于需要处理高精度的问题，以及要严格检验近似估算方法的先进性的时候，渐近精确方法(如马尔科夫蒙特卡洛)就是唯一的解决方案。但是马尔科夫蒙特卡洛方法计算非常耗时，难以解决大尺度高维度问题。
发明内容
发明目的：为了克服背景技术的不足，本发明公开了一种加速刻画高斯水文地质参数场的贝叶斯逆模拟方法，该方法通过耦合并行退火技术，通过交换冷链和热链的采样信息从而使得目标链能够探索到更大的模型空间使得马尔科夫目标链能够探索到更大的模型空间，从而提高探索目标后验分布来提高具有高斯分布的水文地质参数场的反演效率。
技术方案：本发明的一种加速刻画高斯水文地质参数场的贝叶斯逆模拟方法，包括以下步骤：
S1、建立初始高斯参数场、设定平行退火技术中需要的温度梯度以及预条件克兰克尼科尔森建议分布计算中的跳跃梯度；
S2、通过预条件克兰克尼科尔森建议分布得到当前的建议样本，用于计算建议样本的似然度；
S3、根据建议分布和当前分布的似然度计算所有马尔科夫链的当前的接受概率，再根据Metroplis Hastings法则，判断新的样本是建议样本，还是当前样本；
S4、根据设定的交换频率，通过并行退火技术计算任意一对不同温度的马尔科夫链的交换接受概率，再根据Metroplis Hastings法则，判断冷链、热链是否交换，使得目标链能够尽快的探索整个模型空间；
S5、重复以上过程进行下一个高斯场样本采样，根据设定的马尔科夫链的长度，判断是否终止取样。
其中，S1具体为：
设立初始高斯场
Figure GDA0003045900430000021
服从均值为0，方差为C的高斯分布，即
Figure GDA0003045900430000022
设立一个温度梯度，T1<T2<…<Ti<…Tn，T1＝1，设立一个跳跃因子梯度，β1<β2<…<βi<…βn，βn<1，
式中n且为总的马尔科夫链的个数，i表示马尔科夫链的编号，i∈(1,n)，并且交换建议频率为d。
进一步的，S2具体为：
在第kth个样本迭代生成一个预条件克兰克尼科尔森建议分布
Figure GDA0003045900430000023
且对于所有马尔科夫链i＝1,...,n，
Figure GDA0003045900430000024
进一步的，S3具体为：
对于每个马尔科夫链i，当满足接受概率
Figure GDA0003045900430000025
时，则下一个样本
Figure GDA0003045900430000026
否者下一个样本维持不变，即
Figure GDA0003045900430000027
Figure GDA0003045900430000028
式中
Figure GDA0003045900430000029
为建议分布的对数似然度，
Figure GDA00030459004300000210
为当前分布的对数似然度。
进一步的，S4具体为：
对于任意一对不同温度的马尔科夫链，第ith和第jth条马尔科夫链，如果
Figure GDA00030459004300000211
且当满足交换接受概率
Figure GDA00030459004300000212
时，这两条链进行交换
Figure GDA00030459004300000213
否者维持不变，
Figure GDA0003045900430000031
有益效果：与现有技术相比，本发明的优点为：首先，本发明能够刻画水文地质参数高斯场，能够处理水文地质参数反演中的高维度线性及非线性高斯场问题；其次，本方法是一个渐进精确的采样方法，在处理地质统计学逆模拟问题中，比预条件克兰克尼科尔森马尔科夫蒙特卡洛方法效率更高。
附图说明
图1为实施例中观测井和抽水井分布示意图；
图2为本发明中贝叶斯逆模拟方法示意图；
图3为实施例中对数水力渗透系数参考场的频率直方图；
图4为实施例中pCN-MCMC和pCN-PT采样的对数水力渗透系数频率直方图；
图5为实施例中pCN-MCMC和pCN-PT采样的均方根误差；
图6为实施例中pCN-MCMC和pCN-PT潜在尺度因子R值的变化曲线。
具体实施方式
下面结合附图和实施例对本发明的技术方案作进一步的说明。
为了验证本发明是可以加速刻画高斯水文地质参数场的贝叶斯逆模拟方法，在一个饱和，稳定流的承压含水层中验证，该含水层空间二维为5000m*5000m，厚度为50m，且离散为100*100*1个单元格，每个单元格为50m*50m*50m，东、西边界为隔水边界且给定水位分别为20m和0m，南、北边界为隔水边界，并如图1所示，给定25口观测井和4口抽水井，圆圈代表观察井，星号代表抽水井，且抽水井#1、#2、#3、#4的抽水120m3/d，70m3/d，90m3/d，90m3/d。通过序贯高斯模拟[Deutsch and Journel,1998]，并依据均值为2.5ln[m/d]，标准差为2ln[m/d]，最大相关长度和最小相关长度分别为2000m和1500m，角度为135度的指数变异函数模型生成高维度的服从高斯分布的对数水文渗透系数参考场。
通过地下水流动模拟器MODFLOW[McDonald and Harbaugh，1988]计算如下地下水稳定流流动方程式：
Figure GDA0003045900430000032
式中：
Figure GDA0003045900430000033
·表示散度算子；K为水力渗透系数[m/d]；
Figure GDA0003045900430000034
为拉卜拉算子；H为水位值[m]；q为含水层单位体积的体积流量[1/d]。
pCN-PT：耦合并行退火技术的预条件克兰克尼科尔森马尔科夫蒙特卡洛方法
pCN-MCMC：传统的预条件克兰克尼科尔森马尔科夫蒙特卡洛方法
如图2所示的本发明的加速刻画高斯水文地质参数场的贝叶斯逆模拟方法(pCN-PT)，包括步骤如下：
S1、建立初始高斯参数场、设定平行退火技术中需要的温度梯度以及预条件克兰克尼科尔森建议分布计算中的跳跃梯度；
S2、通过预条件克兰克尼科尔森建议分布得到当前的建议样本，用于计算建议样本的似然度；
S3、根据建议分布和当前分布的似然度计算所有马尔科夫链的当前的接受概率，再根据Metroplis Hastings法则，判断新的样本是建议样本，还是当前样本；
S4、根据设定的交换频率，通过并行退火技术计算任意一对不同温度的马尔科夫链的交换接受概率，再根据Metroplis Hastings法则，判断冷链、热链是否交换，使得目标链能够尽快的探索整个模型空间；
S5、重复以上过程进行下一个高斯场样本采样，根据设定的马尔科夫链的长度，判断是否终止取样。
具体实施时：
(1)设立初始高斯场
Figure GDA0003045900430000041
服从均值为0，方差为C的高斯分布，即
Figure GDA0003045900430000042
设立一个温度梯度且T1<T2<…<Ti<…Tn，且T1＝1；设立一个跳跃因子梯度β1<β2<…<βi<…βn，且βn<1，式中n且为总的马尔科夫链的个数，i表示马尔科夫链的编号，i∈(1,n)，并且交换建议频率为d。
(2)在第kth个样本迭代生成一个预条件克兰克尼科尔森建议分布
Figure GDA0003045900430000043
且对于所有马尔科夫链i＝1,...,n，
Figure GDA0003045900430000044
(3)对于每个马尔科夫链i，当满足接受概率
Figure GDA0003045900430000045
时，则下一个样本
Figure GDA0003045900430000046
否者下一个样本维持不变，即
Figure GDA0003045900430000047
Figure GDA0003045900430000048
式中
Figure GDA0003045900430000049
为建议分布的对数似然度；
Figure GDA00030459004300000410
为当前分布的对数似然度。
(4)对于任意一对不同温度的马尔科夫链，例如，第ith和第jth条马尔科夫链，如果
Figure GDA00030459004300000411
且当满足交换接受概率
Figure GDA00030459004300000412
时，这两条链进行交换
Figure GDA00030459004300000413
否者维持不变
Figure GDA00030459004300000414
(5)重复以上过程进行下一个高斯场样本采样。
对于本发明的方法pCN-PT，设定20条并行的马尔科夫链，对于温度梯度为1.76i-1，跳跃因子成指数衰减且为0.7820-i+1。
pCN-MCMC的跳跃因子为0.005。采样的样本总数为800000个。
如图3所示，显示的是对数水力渗透系数参考场的频率直方图，如图4所示，显示的是通过pCN-MCMC和pCN-PT采样的对数水力渗透系数场的频率直方图(已扣除马尔科夫链的burn-in阶段)，可以看出这两种方法采样之后都能获得参考场的高斯分布特征且统计信息基本一致。如图5所示，显示的是pCN-MCMC和pCN-PT样本的均方根误差，可以看出通过pCN-PT方法计算的样本均分根误差能够更快的减小并且收敛到均方根误差均值。如图6所示，显示的是两种方法在随着马尔科夫链增长的潜在尺度因子(potential scale factor)R变化曲线(当R值小于1.2是表明采样样本达到稳定分布)，可以看出，虽然在这两种方法中，R值都随着马尔科夫链的增长而减小，但是pCN-PT方法采样计算的R值在采样结束前大部分都能小于1.2，而pCN-MCMC计算的R值在采样结束之前，大部分都还大于1.2，由此可知，发明的新方法pCN-PT能够更快的达到一个稳定分布。
通过以上分析发明的新方法pCN-PT能够更快速高效的处理高维度的非线性高斯场问题。

Claims (5)

1.一种加速刻画高斯水文地质参数场的贝叶斯逆模拟方法，其特征在于，包括以下步骤：
S1、建立初始地下含水层高斯渗透系数场、设定并 行退火技术中需要的温度梯度以及预条件克兰克尼科尔森建议分布计算中的跳跃梯度；
S2、通过预条件克兰克尼科尔森建议分布得到当前的高斯渗透系数场的建议样本，用于计算建议样本的似然度；
S3、根据渗透系数场的建议分布和当前分布的似然度计算所有马尔科夫链的当前的接受概率，再根据Metropolis Hastings法则，判断新的高斯渗透系数场样本是建议样本，还是当前样本；
S4、根据设定的交换频率，通过并行退火技术计算任意一对不同温度的马尔科夫链的交换接受概率，再根据Metropolis Hastings法则，判断冷链、热链是否交换，使得目标链能够尽快的探索整个地下水系统模型空间；
S5、重复以上过程进行下一个高斯场样本采样，根据设定的马尔科夫链的长度，判断是否终止取样，得到最终地下含水层高斯渗透系数场。
2.根据权利要求1所述的加速刻画高斯水文地质参数场的贝叶斯逆模拟方法，其特征在于，S1具体为：
设立初始地下含水层高斯渗透系数
Figure FDA0003111651280000011
服从均值为0，方差为C的高斯分布，即
Figure FDA0003111651280000012
设立一个温度梯度，T1<T2<…<Ti<…Tn，T1＝1，设立一个跳跃因子梯度，β1<β2<…<βi<…βn，βn<1，
式中n为总的马尔科夫链的个数，i表示马尔科夫链的编号，i∈(1,n)，并且交换建议频率为d。
3.根据权利要求2所述的加速刻画高斯水文地质参数场的贝叶斯逆模拟方法，其特征在于，S2具体为：
在第kth个高斯渗透系数样本迭代生成一个预条件克兰克尼科尔森建议分布
Figure FDA0003111651280000016
且对于所有马尔科夫链i＝1,...,n，
Figure FDA0003111651280000013
4.根据权利要求3所述的加速刻画高斯水文地质参数场的贝叶斯逆模拟方法，其特征在于，S3具体为：
对于每个马尔科夫链i，当满足接受概率
Figure FDA0003111651280000014
时，则下一个样本
Figure FDA0003111651280000015
否则下一个样本维持不变，即
Figure FDA0003111651280000021
Figure FDA0003111651280000022
式中
Figure FDA0003111651280000023
为建议分布的对数似然度，
Figure FDA0003111651280000024
为当前分布的对数似然度。
5.根据权利要求4所述的加速刻画高斯水文地质参数场的贝叶斯逆模拟方法，其特征在于，S4具体为：
对于任意一对不同温度的马尔科夫链，第ith和第jth条马尔科夫链，如果
Figure FDA0003111651280000025
且当满足交换接受概率
Figure FDA0003111651280000026
时，这两条链进行交换
Figure FDA0003111651280000027
否则维持不变，
Figure FDA0003111651280000028
CN202010602434.4A 2020-06-29 2020-06-29 一种加速刻画高斯水文地质参数场的贝叶斯逆模拟方法 Active CN111914475B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010602434.4A CN111914475B (zh) 2020-06-29 2020-06-29 一种加速刻画高斯水文地质参数场的贝叶斯逆模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010602434.4A CN111914475B (zh) 2020-06-29 2020-06-29 一种加速刻画高斯水文地质参数场的贝叶斯逆模拟方法

Publications (2)

Publication Number Publication Date
CN111914475A CN111914475A (zh) 2020-11-10
CN111914475B true CN111914475B (zh) 2021-09-07

Family

ID=73226697

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010602434.4A Active CN111914475B (zh) 2020-06-29 2020-06-29 一种加速刻画高斯水文地质参数场的贝叶斯逆模拟方法

Country Status (1)

Country Link
CN (1) CN111914475B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114486680B (zh) * 2022-01-24 2022-09-13 北京市生态环境保护科学研究院 深层非饱和带岩土渗透系数的确定方法、装置及电子设备

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101793977A (zh) * 2010-01-14 2010-08-04 南京大学 一种水文地质参数的估计方法
CN101819682A (zh) * 2010-04-09 2010-09-01 哈尔滨工程大学 基于马尔科夫链蒙特卡洛粒子滤波的目标跟踪方法
CN102999477B (zh) * 2012-12-21 2016-05-25 中国科学院计算机网络信息中心 一种基于mcmc的并行分类方法
CN107016205A (zh) * 2017-04-14 2017-08-04 宋凯 一种地下水数值模拟的多模型构建方法
CN111177646B (zh) * 2019-12-30 2023-07-28 武汉市陆刻科技有限公司 一种岩溶含水层渗透场反演优化方法

Also Published As

Publication number Publication date
CN111914475A (zh) 2020-11-10

Similar Documents

Publication Publication Date Title
CN111950140B (zh) 一种考虑多种不确定性的大坝渗流性态分析方法
Martino et al. A fast universal self-tuned sampler within Gibbs sampling
CN101356537A (zh) 使用集合卡尔曼滤波进行实时油藏模型更新的方法、系统和装置
WO2012112736A2 (en) System and method for uncertainty quantification in reservoir simulation
CN104933323B (zh) 融合产品成败型数据和故障时间数据的可靠性评估方法
CN111914475B (zh) 一种加速刻画高斯水文地质参数场的贝叶斯逆模拟方法
CN113360983A (zh) 一种边坡可靠度分析与风险评估方法
CN114611832B (zh) 一种基于贝叶斯多模型集对分析的海水入侵预测方法
Stramer et al. On simulated likelihood of discretely observed diffusion processes and comparison to closed-form approximation
CN118464716B (zh) 一种土壤环境污染在线检测方法及检测装置
Blair et al. A Bayesian approach to electrical conductivity relaxation and isotope exchange/secondary ion mass spectrometry
Zhang et al. Efficient surrogate modeling based on improved vision transformer neural network for history matching
CN117686442A (zh) 一种氯离子扩散浓度检测方法、系统、介质及设备
Kim et al. Jointly calibrating hydrologic model parameters and state adjustments
CN117829322A (zh) 基于周期性时间序列与多维度的关联型数据预测方法
Man et al. A higher‐order predictor–corrector scheme for two‐dimensional advection–diffusion equation
Zheng et al. Identification of the degradation coefficient for an anomalous diffusion process in hydrology
Su et al. Approximate analytical solution and parameter estimation for one‐dimensional horizontal absorption based on the van Genuchten–Mualem model
Tarkhamtham et al. High-order generalized maximum entropy estimator in kink regression model
CN113496070B (zh) 地层俘获截面曲线的处理方法、装置、设备及介质
CN114548578A (zh) 基于梯度增强多元回归模型的北斗监测滑坡位移量预测方法
CN109783875A (zh) 一种海洋重力数据的自我迭代更新优化算法
Dey et al. Modified estimation and confidence intervals of an asymmetric loss‐based process capability index C pm′ C^′_pm
Su et al. Sequential implicit sampling methods for Bayesian inverse problems
Naesseth et al. Capacity estimation of two-dimensional channels using sequential Monte Carlo

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant