CN111177646A - 一种岩溶含水层渗透场反演优化方法 - Google Patents
一种岩溶含水层渗透场反演优化方法 Download PDFInfo
- Publication number
- CN111177646A CN111177646A CN201911391154.7A CN201911391154A CN111177646A CN 111177646 A CN111177646 A CN 111177646A CN 201911391154 A CN201911391154 A CN 201911391154A CN 111177646 A CN111177646 A CN 111177646A
- Authority
- CN
- China
- Prior art keywords
- probability distribution
- distribution function
- permeability
- function
- reference probability
- 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.)
- Granted
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/15—Correlation function computation including computation of convolution operations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A10/00—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE at coastal zones; at river basins
- Y02A10/40—Controlling or monitoring, e.g. of flood or hurricane; Forecasting, e.g. risk assessment or mapping
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Computational Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Algebra (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Computing Systems (AREA)
- Complex Calculations (AREA)
- Amplifiers (AREA)
- Feedback Control In General (AREA)
Abstract
本发明提出了一种岩溶含水层渗透场反演优化方法,在传统的Metropolis‑Hasting算法的基础上,通过改进参考概率分布对岩溶含水层渗透场进行反演,提高了反演模型结果的精度;本发明可以在保证马尔科夫链收敛性的同时提高抽样的全局能力和多链抽样的效率,并获取高精度的反演模型结果。
Description
技术领域
本发明涉及水文地质学参数优化(即反演)领域,提出了一种新型的岩溶含水层渗透场反演优化方法。
背景技术
随着环境水文地质学正问题研究日益成熟,其相关的反问题也日益受到重视。水文地质学反问题是根据基于观测数据对水力传导率分布进行估计。对于岩溶含水层,由于内部岩溶管道、裂缝存在不同程度的发育,含水层的水力学特性会相差几个数量级,水力参数评估十分困难。
贝叶斯推理是处理非线性、不确定性系统反问题一种行之有效的概率方法。该方法通过对未知参数的后验概率函数进行抽样来获得参数的估计分布。MCMC(马尔科夫链蒙特卡罗)算法是一种基于贝叶斯推理的标准抽样方法,但是由于高维反演问题的特强非线性以及MCMC算法中计算参数设置的人为性,抽样结果是否具有对后验概率分布的代表性与马尔科夫(马尔可夫)链设计者的经验有很大关系。基于Metropolis-Hastings算法的MCMC抽样方法采用多链搜索的方法对其进行了改进。但是,对于高维度非常复杂的参数空间,基于Metropolis-Hasting算法的实际随机游走往往不能对其充分搜索,从而陷入一个局部高密度区域。因此需要在传统的Metropolis-Hasting算法的基础上进行改进,提高反演模型结果的精度。
发明内容
有鉴于此,本发明提出了一种岩溶含水层渗透场反演优化方法,能在保证马尔科夫链收敛性的同时提高抽样的全局能力和多链抽样的效率,解决了现有反演模型结果精度不够准确的问题。
本发明的技术方案是这样实现的:本发明提供了一种岩溶含水层渗透场反演优化方法,包括以下步骤:
S1,从岩溶含水层渗透场模型获取渗透率数据(如量级和大概空间分布等),通过贝叶斯推理,根据该渗透数据获取后验概率密度函数,通过将后验概率密度函数局部近似为高斯分布,获取初始参考概率分布函数;
S2,通过低秩矩阵近似法,根据渗透率数据对初始参考概率分布函数中的Hessian矩阵局部进行改写,获取改写后的参考概率分布函数;
S3,获取参考概率分布函数的方差数据,通过截断牛顿法将方差数据引入改写后的参考概率分布函数中,获取最终参考概率分布函数;
S4,通过最终参考概率分布函数对岩溶含水层渗透场模型进行反演。
在以上技术方案的基础上,优选的,步骤S1中,从岩溶含水层渗透场模型获取渗透率数据(如量级和大概空间分布等),通过贝叶斯推理,包括以下步骤:
获取渗透数据,所述渗透数据包括:水头观测结果hobs=[h1,h2,...,hM]以及对渗透率s=[log(P1),log(P2),...,log(PN)],其中,M表示测量点的数目,N表示中待拟合点的数目,P代表渗透率值;
通过贝叶斯推理,后验概率密度函数∏post(s)正比于先验概率密度函数Πprior(s)和似然函数Π(hobs|s)的积,即:Πpost(s)∝Πprior(s)∏(hobs|s);
通过L2准则可以将后验概率密度函数进一步表示为:
其中,Γnoise以及Γprior都为协方差矩阵。
进一步优选的,步骤S1中,通过将后验概率密度函数局部近似为高斯分布,获取初始参考概率分布函数,包括以下步骤:
将后验概率密度函数局部近似为高斯分布,获取基于局部梯度信息的初始参考概率分布函数:
进一步优选的,局部梯度g可以表示为:
进一步优选的,步骤S2中,通过低秩矩阵近似法,根据渗透数据对初始参考概率分布函数中的Hessian矩阵局部进行改写,包括以下步骤:
在统计设置中,Hessian矩阵由两部分组成:
其中,Hmisfit表示参数字段的子空间,该参数字段受渗透数据的约束,通过低秩矩阵近似法,根据参数字段对初始参考概率分布函数中的Hessian矩阵局部进行改写。
进一步优选的,改写后的参考概率分布函数为:
H=L-T(LTHmisfitL+I)L-1;
其中,L是先验协方差矩阵Γprior的Cholesky因式分解,T代表转置,I=LTL。
进一步优选的,获取参考概率分布函数的方差数据,通过截断牛顿法将方差数据引入改写后的参考概率分布函数中,获取最终参考概率分布函数包括以下步骤:最终参考概率分布函数为:
其中,σ2为参考概率分布函数的方差。
本发明的一种岩溶含水层渗透场反演优化方法相对于现有技术具有以下有益效果:
(1)通过低秩矩阵近似法以及截断牛顿法对参考概率分布函数进行改进,通过改进后的参考概率分布函数对岩溶含水层渗透场进行反演,提高反演模型的精度;
(2)通过低秩矩阵近似法以及截断牛顿法对参考概率分布函数进行改进,使改进后的参考概率分布函数保证马尔科夫链收敛性的同时提高抽样的全局能力和多链抽样的效率。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其他的附图。
图1为本发明一种岩溶含水层渗透场反演优化方法的流程图。
具体实施方式
下面将结合本发明实施方式,对本发明实施方式中的技术方案进行清楚、完整地描述,显然,所描述的实施方式仅仅是本发明一部分实施方式,而不是全部的实施方式。基于本发明中的实施方式,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施方式,都属于本发明保护的范围。
如图1所示,本发明的一种岩溶含水层渗透场反演优化方法,包括以下步骤:
S1,从岩溶含水层渗透场模型获取渗透率数据(如量级和大概空间分布等),通过贝叶斯推理,根据该渗透数据获取后验概率密度函数,通过将后验概率密度函数局部近似为高斯分布,获取初始参考概率分布函数;具体包括:
获取含水层渗透属性数据,所述含水层渗透属性数据包括:水头观测结果hobs=[h1,h2,...,hM]以及对渗透率s=[log(P1),log(P2),...,log(PN)],其中,M表示测量点的数目,N表示中待拟合点的数目,P代表渗透率值;
通过贝叶斯推理,后验概率密度函数Πpost(s)正比于先验概率密度函数Πprior(s)和似然函数Π(hobs|s)的积,即:∏post(s)∝∏prior(s)∏(hobs|s);
通过L2准则可以将后验概率密度函数进一步表示为:
其中,Γnoise以及Γprior都为协方差矩阵;
将后验概率密度函数局部近似为高斯分布,获取基于局部梯度信息的初始参考概率分布函数:
局部梯度g可以表示为:
应当理解的是,贝叶斯推断(Bayesian inference)是推论统计的一种方法。这种方法使用贝叶斯定理,在有更多证据及信息时,更新特定假设的概率。贝叶斯推断是统计学(特别是数理统计学)中很重要的技巧之一。贝叶斯更新(Bayesian updating)在序列分析中格外的重要。贝叶斯推断应用在许多的领域中,包括科学、工程学、哲学、医学、体育运动、法律等。在决策论的哲学中,贝叶斯推断和主观概率有密切关系,常常称为贝叶斯概率。
应当理解的是,假设测量误差是可以累计的,表现为均值为0、协方差矩阵为Γnoise的高斯分布。同时假设先验概率密度函数为均值为协方差矩阵为Γprior的高斯分布,则后验概率密度函数可以用L2准则进一步表示为:
该后验概率密度函数即为反问题的统计解。为了获得参数的估计值,需要对后验概率分布进行抽样。目前,马尔可夫链蒙特卡罗方法(MCMC)是后验概率分布最有效的方法。Metropolis-Hasting算法是应用最为广泛的MCMC抽样方法,该方法将先验模型看做一个黑箱,并且采用高斯参考概率分布函数。
应当理解的是,反问题的统计解为当后验概率密度函数最大化时,即指数项最小化。其中,如果第一项为误差协方差倒数加权的数据失配项,第二项为先验概率函数协方差倒数加权的确定性优化的正则化项,这相当于,这相当于将确定性逆问题的最小二乘成本函数最小化。此外,确定性代价函数的Hessian矩阵的逆是后验密度协方差的近似值。
S2,通过低秩矩阵近似法,根据渗透率数据对初始参考概率分布函数中的Hessian矩阵局部进行改写,获取改写后的参考概率分布函数:
在统计设置中,Hessian矩阵由两部分组成:
其中,Hmisfit表示参数字段的子空间,该参数字段受渗透数据的约束,通过低秩矩阵近似法,根据参数字段对初始参考概率分布函数中的Hessian矩阵局部进行改写;
改写后的参考概率分布函数为:
H=L-T(LTHmisfitL+I)L-1;
其中,L是先验协方差矩阵Γprior的Cholesky因式分解,T代表转置,I=LTL;
应当理解的是,从参考概率分布函数q(sk,y)中抽取样本,需要正定的Hessian矩阵然而,该条件仅在最小二乘代价函数的局部最小值附近得到满足,在参数空间的任意位置处,不能保证Hessian矩阵是正定的。因此,为了确保马尔可夫链的收敛,需要对局部Hessian进行适当的修改。另一方面,从确定性反问题中获得的证据表明,等式的数据拟合项的Hessian,是一个紧凑的算子。实际上,对于许多欠定的逆问题,在给定的反演步骤中仅通知参数字段的有限数量的元素。因此,用其低级变量替换局部Hessian矩阵是合理的。
应当理解的是,对于Vr仅包含主要的特征向量,而Dr的对角元素代表相对应的主要特征值,下标r代表的秩的阶数,其通过截断分解来确定,小于系统设定阈值α的特征值被认为是主要受先验约束,所以这些特征值会被消灭,因此这种截断仍能够确保得到的低秩Hessian矩阵是正定的。
设置初始值T0,对于每个节点k=0,...,N-1循环进行:
1.计算▽log(π(Tk))和H(Tk);
2.从初始参考概率分布函数q(sk,y)中抽取一个样本y;
3.计算∏(y),q(Tk,y)和q(y,Tk);
5.用α(Tk,y)接收和设置Tk+1=y,否则和设置Tk+1=Tk;
结束循环。
S3,获取参考概率分布函数的方差数据,通过截断牛顿法将方差数据引入改写后的参考概率分布函数中,获取最终参考概率分布函数:
其中,σ2为参考概率分布函数的方差。
应当理解的是,随机牛顿方法是基于目标概率密度是高斯分布的基本假设得出的。当目标密度接近高斯时,初始参考概率分布函数的确定项和随机项都很好地工作。然而,对于先验模型是非线性且参数空间维度大的水力学问题,目标密度很可能与高斯不同。在这种情况下,如果仍然应用完全牛顿步骤马尔可夫链可能被困在低概率区域。因此截断牛顿步骤是合理的。另一方面,如果目标密度不是高斯分布,则需要调整随机项,否则,该初始参考概率分布函数可能会被拒绝。
应当理解的是,对于大规模问题,在过渡阶段,即在提议开始探索静止密度分布之前,提议方差最优地缩放为σ2=l2N-1/2作为参数的数量接近无限,并且在静止阶段中按比例σ2=l2N-1/3。这里,l是缩放常数,N表示未知参数的数量。实现时,动态调整的值以确保MCMC链的良好收敛:如果样本被拒绝,则l在下一次迭代中减半,直到链演变为新点;一旦生成了接受的样本,则将l设置回初始值。
S4,通过最终参考概率分布函数对岩溶含水层渗透场模型进行反演。
应当理解的是,通过最终参考概率分布函数对岩溶含水层渗透场模型进行反演,可以提高反演模型结果的精确度。
以上所述仅为本发明的较佳实施方式而已,并不用以限制本发明,凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。
Claims (8)
1.一种岩溶含水层渗透场反演优化方法,其特征在于:
S1,从岩溶含水层渗透场模型获取渗透率数据(如量级和大概空间分布等),通过贝叶斯推理,根据该渗透数据获取后验概率密度函数,通过将后验概率密度函数局部近似为高斯分布,获取初始参考概率分布函数;
S2,通过低秩矩阵近似法,根据渗透率数据对初始参考概率分布函数中的Hessian矩阵局部进行改写,获取改写后的参考概率分布函数;
S3,获取参考概率分布函数的方差数据,通过截断牛顿法将方差数据引入改写后的参考概率分布函数中,获取最终参考概率分布函数;
S4,通过最终参考概率分布函数对岩溶含水层渗透场模型进行反演。
2.如权利要求1所述的一种岩溶含水层渗透场反演优化方法,其特征在于:步骤S1中,从岩溶含水层渗透场模型获取渗透率数据(如量级和大概空间分布等),通过贝叶斯推理,包括以下步骤:
获取渗透数据,所述渗透数据包括:水头观测结果hobs=[h1,h2,...,hM]以及对渗透率s=[log(P1),log(P2),...,log(PN)],其中,M表示测量点的数目,N表示中待拟合点的数目,P代表渗透率值;
通过贝叶斯推理,后验概率密度函数Πpost(s)正比于先验概率密度函数Πprior(s)和似然函数Π(hobs|s)的积,即:Πpost(s)∝Πprior(s)Π(hobs|s);
通过L2准则可以将后验概率密度函数进一步表示为:
其中,Γnoise以及Γprior为协方差矩阵。
6.如权利要求5所述的一种岩溶含水层渗透场反演优化方法,其特征在于:改写后的参考概率分布函数为:
H=L-T(LTHmisfitL+I)L-1;
其中,L是先验协方差矩阵Γprior的Cholesky因式分解,T代表转置,I=LTL。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201911391154.7A CN111177646B (zh) | 2019-12-30 | 2019-12-30 | 一种岩溶含水层渗透场反演优化方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201911391154.7A CN111177646B (zh) | 2019-12-30 | 2019-12-30 | 一种岩溶含水层渗透场反演优化方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN111177646A true CN111177646A (zh) | 2020-05-19 |
CN111177646B CN111177646B (zh) | 2023-07-28 |
Family
ID=70647420
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201911391154.7A Active CN111177646B (zh) | 2019-12-30 | 2019-12-30 | 一种岩溶含水层渗透场反演优化方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN111177646B (zh) |
Cited By (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111914475A (zh) * | 2020-06-29 | 2020-11-10 | 河海大学 | 一种加速刻画高斯水文地质参数场的贝叶斯逆模拟方法 |
CN112307536A (zh) * | 2020-09-18 | 2021-02-02 | 天津大学 | 一种大坝渗流参数反演方法 |
CN118050300A (zh) * | 2024-04-16 | 2024-05-17 | 河北天辰仪器设备有限公司 | 一种土工布智能垂直渗透系数测定方法及测定装置 |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US7254091B1 (en) * | 2006-06-08 | 2007-08-07 | Bhp Billiton Innovation Pty Ltd. | Method for estimating and/or reducing uncertainty in reservoir models of potential petroleum reservoirs |
US20130191095A1 (en) * | 2010-09-27 | 2013-07-25 | Total Sa | Simulation of a geological phenomenon |
CN106153859A (zh) * | 2016-06-22 | 2016-11-23 | 福建工程学院 | 一种盐岩动水溶蚀的试验装置和计算方法 |
US20170140079A1 (en) * | 2014-06-30 | 2017-05-18 | Cgg Services Sas | Ensemble-based multi-scale history-matching device and method for reservoir characterization |
-
2019
- 2019-12-30 CN CN201911391154.7A patent/CN111177646B/zh active Active
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US7254091B1 (en) * | 2006-06-08 | 2007-08-07 | Bhp Billiton Innovation Pty Ltd. | Method for estimating and/or reducing uncertainty in reservoir models of potential petroleum reservoirs |
US20130191095A1 (en) * | 2010-09-27 | 2013-07-25 | Total Sa | Simulation of a geological phenomenon |
US20170140079A1 (en) * | 2014-06-30 | 2017-05-18 | Cgg Services Sas | Ensemble-based multi-scale history-matching device and method for reservoir characterization |
CN106153859A (zh) * | 2016-06-22 | 2016-11-23 | 福建工程学院 | 一种盐岩动水溶蚀的试验装置和计算方法 |
Non-Patent Citations (1)
Title |
---|
马琦琦: "《改进的贝叶斯迭代反演方法及其在白云岩致密储层识别的应用》" * |
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111914475A (zh) * | 2020-06-29 | 2020-11-10 | 河海大学 | 一种加速刻画高斯水文地质参数场的贝叶斯逆模拟方法 |
CN112307536A (zh) * | 2020-09-18 | 2021-02-02 | 天津大学 | 一种大坝渗流参数反演方法 |
CN112307536B (zh) * | 2020-09-18 | 2022-11-22 | 天津大学 | 一种大坝渗流参数反演方法 |
CN118050300A (zh) * | 2024-04-16 | 2024-05-17 | 河北天辰仪器设备有限公司 | 一种土工布智能垂直渗透系数测定方法及测定装置 |
Also Published As
Publication number | Publication date |
---|---|
CN111177646B (zh) | 2023-07-28 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
US10852439B1 (en) | Global ionospheric total electron content prediction system | |
CN108761574B (zh) | 基于多源信息融合的降雨量估算方法 | |
Moradkhani et al. | Evolution of ensemble data assimilation for uncertainty quantification using the particle filter‐Markov chain Monte Carlo method | |
Pathiraja et al. | Data‐driven model uncertainty estimation in hydrologic data assimilation | |
Jeremiah et al. | Bayesian calibration and uncertainty analysis of hydrological models: A comparison of adaptive Metropolis and sequential Monte Carlo samplers | |
US9317929B2 (en) | Decomposition apparatus and method for refining composition of mixed pixels in remote sensing images | |
Welling et al. | Bayesian learning via stochastic gradient Langevin dynamics | |
CN111177646A (zh) | 一种岩溶含水层渗透场反演优化方法 | |
Kang et al. | Bayesian inference for the spatial random effects model | |
Rings et al. | Bayesian model averaging using particle filtering and Gaussian mixture modeling: Theory, concepts, and simulation experiments | |
CN110232471B (zh) | 一种降水传感网节点布局优化方法及装置 | |
CN106203511A (zh) | 一种图像相似块评估方法 | |
Elmi et al. | Spaceborne river discharge from a nonparametric stochastic quantile mapping function | |
Wei et al. | A comparative assessment of multisensor data merging and fusion algorithms for high-resolution surface reflectance data | |
Levy et al. | Using deep generative neural networks to account for model errors in Markov chain Monte Carlo inversion | |
Dehkordi et al. | Machine learning-based estimation of suspended sediment concentration along Missouri River using remote sensing imageries in Google Earth Engine | |
Denes et al. | Persistence and material coherence of a mesoscale ocean eddy | |
Reggiani et al. | A Bayesian processor of uncertainty for precipitation forecasting using multiple predictors and censoring | |
CN115825388A (zh) | 一种重金属估算模型的训练方法、估算方法、装置及设备 | |
CN116089832A (zh) | 重力卫星地下水储量的降尺度方法、装置和计算机设备 | |
Buehner et al. | Sea ice data assimilation | |
CN112287607B (zh) | 一种多尺度电离层tec预测方法 | |
Hase et al. | Atmospheric inverse modeling via sparse reconstruction | |
Barillec et al. | Data assimilation for precipitation nowcasting using Bayesian inference | |
Laine et al. | Aerosol model selection and uncertainty modelling by adaptive MCMC technique |
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 |