CN113433587B - 岩石物理反演方法及系统 - Google Patents

岩石物理反演方法及系统 Download PDF

Info

Publication number
CN113433587B
CN113433587B CN202010207180.6A CN202010207180A CN113433587B CN 113433587 B CN113433587 B CN 113433587B CN 202010207180 A CN202010207180 A CN 202010207180A CN 113433587 B CN113433587 B CN 113433587B
Authority
CN
China
Prior art keywords
parameters
rock
probability density
density
bulk modulus
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
CN202010207180.6A
Other languages
English (en)
Other versions
CN113433587A (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.)
Petrochina Co Ltd
Original Assignee
Petrochina Co Ltd
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 Petrochina Co Ltd filed Critical Petrochina Co Ltd
Priority to CN202010207180.6A priority Critical patent/CN113433587B/zh
Publication of CN113433587A publication Critical patent/CN113433587A/zh
Application granted granted Critical
Publication of CN113433587B publication Critical patent/CN113433587B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30Analysis
    • G01V1/306Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/282Application of seismic models, synthetic seismograms
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30Analysis
    • G01V1/303Analysis for determining velocity profiles or travel times
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/60Analysis
    • G01V2210/62Physical property of subsurface
    • G01V2210/622Velocity, density or impedance
    • G01V2210/6222Velocity; travel time
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/60Analysis
    • G01V2210/62Physical property of subsurface
    • G01V2210/622Velocity, density or impedance
    • G01V2210/6224Density
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/60Analysis
    • G01V2210/62Physical property of subsurface
    • G01V2210/624Reservoir parameters
    • G01V2210/6244Porosity

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明提供了一种岩石物理反演方法,所述方法包含:获取待反演地区的地震弹性参数和储层物理参数;通过Gassmann方程和Xu‑White模型建立所述地震弹性参数和所述储层物理参数之间的关系方程;于贝叶斯理论框架下利用ADMM‑MCMC算法迭代所述关系方程中的物性参数计算其对应的后验概率密度数据;根据所述后验概率密度数据中最大后验密度值对所述待反演地区进行反演获得反演结果。

Description

岩石物理反演方法及系统
技术领域
本发明涉及石油地球物理勘探领域,尤指一种岩石物理反演方法及系统。
背景技术
定量地震解释是地震解释的最终目标。地震反演是利用地震数据获得地震弹性参数的重要手段,利用地震弹性参数可以进行储层(岩性)预测、物性预测以及含油气性检测等。特别是在物性参数预测时,通常只能利用地震弹性参数与储层物性参数之间的统计经验关系拟合计算。不同于常规的经验关系拟合,岩石物理反演是利用岩石物理理论模型来构建地震弹性参数与储层物性参数的理论关系。因此岩石物理反演是进行定量地震解释的重要途径。
岩石物理反演主要包括两方面内容:一个是利用岩石物理模型来构建岩石物理反演的反问题方程;另一个就是选取何种方法来求解岩石物理反演的反问题方程。通常岩石物理模型建立的地震弹性参数与储层物性参数之间的关系非常复杂,往往是非线性的。所以很多专家学者都利用了统计学反演方法来求解非线性岩石物理反问题。Mukerji等(2001)首次提出了将岩石物理模型与地质统计学方法结合来反演储层物性参数。Bachrach等(2006)建立了弹性参数与物性参数之间的统计关系,利用贝叶斯估计方法反演孔隙度和含水饱和度。Spikes等(2007)利用接触胶结模型进行岩石物理建模,利用测井数据与地震数据约束的概率化反演方法预测物性参数。Bosch(2009)利用Wyllie时间平均方程构建物性参数与弹性参数之间的关系,利用马尔可夫链蒙特卡罗MCMC方法反演物性参数。邓继新等(2009)提出了在贝叶斯框架下利用统计岩石物理模型反演孔隙度与含气饱和度的方法。印兴耀等(2014)利用弹性阻抗与储层物性参数之间的统计关系,在贝叶斯框架下实现储层物性参数反演。林恬等(2017)利用包含物模型构建岩石物理模型,提出了约束最小二乘与信赖域的储层参数反演方法。Li等(2018)利用等效介质理论构建三维岩石物理模板,利用最小二乘原理寻找三维模版网格点方法预测物性参数。杨培杰(2018)利用Simon模型建立了砂泥岩储层物性参数和地震弹性参数之间的定量关系,基于贝叶斯理论利用蒙特卡罗与遗传算法反演储层孔隙度和含水饱和度。张冰等(2018)提出了各向异性页岩储层统计岩石物理反演方法,在贝叶斯框架下利用模拟退火优化粒子群算法(SA-PSO)来估计储层物性参数以及各向异性参数的空间分布。由于岩石物理反问题很多时候都是非线性的,因此岩石物理反演方法主要存在的问题就是如何选取合适的反演算法求解非线性的岩石物理反演问题,提高反演问题的稳定性和精度。
发明内容
本发明目的在于利用地震弹性参数获取储层物性参数;通过岩石物理模型建立地震弹性参数与储层物性参数之间的理论关系,将交替方向乘子ADMM算法引入到马尔科夫链蒙特卡罗法MCMC算法中,提出了一种基于ADMM-MCMC算法的求解非线性岩石物理反演问题的岩石物理反演方法及系统。
为达上述目的,本发明所提供的岩石物理反演方法具体包含:获取待反演地区的地震弹性参数和储层物理参数;通过Gassmann方程和Xu-White模型建立所述地震弹性参数和所述储层物理参数之间的关系方程;于贝叶斯理论框架下利用ADMM-MCMC算法迭代所述关系方程中的物性参数计算其对应的后验概率密度数据;根据所述后验概率密度数据中最大后验密度值对所述待反演地区进行反演获得反演结果。
在上述岩石物理反演方法中,优选的,通过Gassmann方程和Xu-White模型建立所述地震弹性参数和所述储层物理参数之间的关系方程还包含:根据所述地震弹性参数和所述储层物理参数获得饱和流体的体积模量和密度;通过Voigt-Reuss-Hill平均方法计算获得待反演地区的岩石基质的体积模量、剪切模量和密度;根据所述岩石基质的体积模量和剪切模量,通过Xu-White模型获得岩石骨架的体积模量和剪切模量;根据岩石骨架的体积模量和剪切模量与饱和流体的体积模量,通过Gassmann方程获得饱和岩石的体积模量和剪切模量;根据饱和流体的密度、岩石基质的密度、饱和岩石的体积模量和剪切模量计算获得饱和岩石的纵波速度、横波速度和密度。
在上述岩石物理反演方法中,优选的,于贝叶斯理论框架下利用ADMM-MCMC算法迭代所述关系方程中的物性参数计算其对应的后验概率密度数据包含:通过饱和岩石的纵波速度、横波速度和密度构建观测弹性参数,通过以下关系方程计算后验概率密度数据:
在上式中,m为物性参数,p(m)为物性参数m的先验概率密度,p(dobs|m)为模型弹性参数fRPM(m)与观测弹性参数dobs之间的拟合程度,p(m|dobs)为后验密度值,fRPM(m)为利用岩石物理模型由物性参数计算得到模型弹性参数,dobs为观测弹性参数。
在上述岩石物理反演方法中,优选的,所述先验概率密度通过以下公式计算获得:
在上式中,p(m)为物性参数m的先验概率密度,φ为孔隙度、Sw为饱和度,Vsh为泥质含量,p(φ)为孔隙度的先验概率密度,p(Vsh)为泥质含量的先验概率密度,p(Sw)为饱和度的先验概率密度,μφ为孔隙度均值,为泥质含量均值,/>为饱和度均值,σφ为孔隙度方差,/>为泥质含量方差,/>为饱和度方差,exp为以自然指数e为底的指数函数,n为待估计样点个数,∑为对所有样点求和运算,π为圆周率。
在上述岩石物理反演方法中,优选的,所述拟合程度通过以下公式计算获得:
在上式中,p(dobs|m)为模型弹性参数fRPM(m)与观测弹性参数dobs之间的拟合程度,dobs为观测弹性参数,fRPM(m)为利用岩石物理模型由物性参数计算得到模型弹性参数,σε为误差。
在上述岩石物理反演方法中,优选的,于贝叶斯理论框架下利用ADMM-MCMC算法迭代所述关系方程中的物性参数计算其对应的后验概率密度数据包含:利用ADMM算法构造马尔科夫链,通过马尔科夫链替代MCMC算法中的马尔科夫链后,对所述关系方程中的物性参数进行逐个迭代计算,获得对应的后验概率密度数据。
本发明还提供一种岩石物理反演系统,所述系统包含采集模块、构建模块、计算模块和反演模块;所述采集模块用于获取待反演地区的地震弹性参数和储层物理参数;所述构建模块用于通过Gassmann方程和Xu-White模型建立所述地震弹性参数和所述储层物理参数之间的关系方程;所述计算模块用于于贝叶斯理论框架下利用ADMM-MCMC算法迭代所述关系方程中的物性参数计算其对应的后验概率密度数据;所述反演模块用于根据所述后验概率密度数据中最大后验密度值对所述待反演地区进行反演获得反演结果。
在上述岩石物理反演系统中,优选的,所述构建模块还包含分析单元,所述分析单元用于:根据所述地震弹性参数和所述储层物理参数获得饱和流体的体积模量和密度;通过Voigt-Reuss-Hill平均方法计算获得待反演地区的岩石基质的体积模量、剪切模量和密度;根据所述岩石基质的体积模量和剪切模量,通过Xu-White模型获得岩石骨架的体积模量和剪切模量;根据岩石骨架的体积模量和剪切模量与饱和流体的体积模量,通过Gassmann方程获得饱和岩石的体积模量和剪切模量;根据饱和流体的密度、岩石基质的密度、饱和岩石的体积模量和剪切模量计算获得饱和岩石的纵波速度、横波速度和密度。
在上述岩石物理反演系统中,优选的,所述计算模块还包含:利用ADMM算法构造马尔科夫链,通过马尔科夫链替代MCMC算法中的马尔科夫链后,对所述关系方程中的物性参数进行逐个迭代计算,获得对应的后验概率密度数据。
本发明还提供一种计算机设备,包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序,所述处理器执行所述计算机程序时实现上述方法。
本发明还提供一种计算机可读存储介质,所述计算机可读存储介质存储有执行上述方法的计算机程序。
本发明的有益技术效果在于:基于Gassmann方程和Xu-White模型岩石物理模型建立地震弹性参数与储层物性参数之间的非线性关系。然后在贝叶斯理论框架下,提出了一种基于ADMM-MCMC算法的求解非线性岩石物理反演问题的方法。利用ADMM算法逐个迭代待反演的物性参数,更加符合岩石物理反演,因为孔隙度、饱和度和泥质含量都是彼此独立的参数,这样可以尽量保证寻找每个参数的最优解,而不是所有参数组合的最优解。再基于对迭代生成的平稳分布马尔科夫链进行统计分析,就可以获取待反演参数最大后验密度,进而获得物性参数的最优解及其不确定性分析。实际测井数据和地震数据应用表明,基于ADMM-MCMC的非线性岩石物理反演方法可以获得较为准确的物性参数。
附图说明
此处所说明的附图用来提供对本发明的进一步理解,构成本申请的一部分,并不构成对本发明的限定。在附图中:
图1为本发明一实施例所提供的岩石物理反演方法的流程示意图;
图2为本发明一实施例所提供的参数获取流程示意图;
图3为本发明一实施例所提供的岩石物理反演系统的结构示意图;
图4为本发明一实施例所提供的岩石物理反演方法实际测井数据;
图5为本发明一实施例所提供的对测井曲线进行非线性岩石物理反演示意图;
图6A至图6C为本发明一实施例所提供的地震弹性参数示意图;
图7A至图7C为本发明一实施例所提供的反演结果示意图;
图8为本发明一实施例所提供的电子设备的结构示意图。
具体实施方式
以下将结合附图及实施例来详细说明本发明的实施方式,借此对本发明如何应用技术手段来解决技术问题,并达成技术效果的实现过程能充分理解并据以实施。需要说明的是,只要不构成冲突,本发明中的各个实施例及各实施例中的各个特征可以相互结合,所形成的技术方案均在本发明的保护范围之内。
另外,在附图的流程图示出的步骤可以在诸如一组计算机可执行指令的计算机系统中执行,并且,虽然在流程图中示出了逻辑顺序,但是在某些情况下,可以以不同于此处的顺序执行所示出或描述的步骤。
请参考图1所示,本发明所提供的岩石物理反演方法具体包含:S101获取待反演地区的地震弹性参数和储层物理参数;S102通过Gassmann方程和Xu-White模型建立所述地震弹性参数和所述储层物理参数之间的关系方程;S103于贝叶斯理论框架下利用ADMM-MCMC算法迭代所述关系方程中的物性参数计算其对应的后验概率密度数据;S104根据所述后验概率密度数据中最大后验密度值对所述待反演地区进行反演获得反演结果。具体的,在实际工作中,上述方法主要方法如下:首先利用Gassmann方程和Xu-White模型构建泥质砂岩的岩石物理模型,建立地震弹性参数与储层物性参数之间的非线性关系。再将交替方向乘子ADMM算法引入到马尔科夫链蒙特卡罗法MCMC算法中,提出了一种基于ADMM-MCMC算法的求解非线性岩石物理反演问题的方法。在贝叶斯理论框架下,利用ADMM算法逐个迭代待反演参数,再对迭代生成的平稳分布马尔科夫链进行统计分析,就可以获取待反演参数的最大后验密度,进而获得待反演参数的最优解及其不确定性分析。
请参考图2所示,在本发明一实施例中,通过Gassmann方程和Xu-White模型建立所述地震弹性参数和所述储层物理参数之间的关系方程还包含:S201根据所述地震弹性参数和所述储层物理参数获得饱和流体的体积模量和密度;S202通过Voigt-Reuss-Hill平均方法计算获得待反演地区的岩石基质的体积模量、剪切模量和密度;S203根据所述岩石基质的体积模量和剪切模量,通过Xu-White模型获得岩石骨架的体积模量和剪切模量;根据岩石骨架的体积模量和剪切模量与饱和流体的体积模量,通过Gassmann方程获得饱和岩石的体积模量和剪切模量;S204根据饱和流体的密度、岩石基质的密度、饱和岩石的体积模量和剪切模量计算获得饱和岩石的纵波速度、横波速度和密度。具体的,请参考如下所示:
(1)利用Gassmann方程和Xu-White模型进行岩石物理建模,具体建模过程如下:
①岩石基质的体积模量Km、剪切模量μm和密度ρm可以利用Voigt-Reuss-Hill平均方法计算:
ρm=Vshρc+(1-Vshq (3)
式中,Vsh为泥质含量,Kc和Kq分别为泥质和砂质矿物的体积模量,μc和μq分别为泥质和砂质矿物的剪切模量,ρc和ρq分别为泥质和砂质矿物的密度。
②饱和流体的体积模量Kfl和密度ρfl计算公式为:
ρfl=(1-Swhc+Swρw (5)
式中,Sw为含水饱和度,Khc和Kw分别为油或气等碳氢化合物和水的体积模量,ρhc和ρw分别为油或气等碳氢化合物和水的密度。
③岩石骨架的体积模量Kdry和剪切模量μdry可以利用Xu-White模型计算:
其中,
式中,αq为砂质矿物的孔隙纵横比,αc为泥质矿物的孔隙纵横比,φq为砂质矿物的孔隙度,且φq=Vshφ,φc为泥质矿物的孔隙度,且φc=(1-Vsh)φ,满足φqc=φ,Tijij(α)和Tiijj(α)为孔隙纵横比的函数。
④饱和岩石的体积模量Ksat和剪切模量μsat可以用Gassmann方程计算:
μsat=μdry (9)
⑤饱和岩石的纵波速度VP、横波速度VS和密度ρ计算公式为
ρ=ρm(1-φ)+ρflφ (10)
(2)以此,通过上述获得的参数即可构建岩石物理反演的反问题方程,如下:
dobs=fRPM(m)+ε, (13)
其中,dobs为观测弹性参数,dobs=[VP,VS,ρ],即纵波速度VP、横波速度VS和密度ρ等,m为待估计的物性参数,m=[φ,Vsh,Sw],即孔隙度φ、饱和度Sw和泥质含量Vsh等,fRPM(m)为利用岩石物理模型由物性参数计算得到模型弹性参数;ε为观测弹性参数dobs与模型弹性参数fRPM(m)之间的误差。
根据上述反问题方程,在贝叶斯理论框架下可进行转换,具体的,在本发明一实施例中,上述于贝叶斯理论框架下利用ADMM-MCMC算法迭代所述关系方程中的物性参数计算其对应的后验概率密度数据可进一步包含:通过饱和岩石的纵波速度、横波速度和密度构建观测弹性参数,通过以下关系方程计算后验概率密度数据:
在上式中,p(m)为物性参数m的先验概率密度,即未获得观测弹性参数dobs时物性参数m的分布规律,一般根据已有的信息或者经验知识确定的先验信息,p(dobs|m)为似然函数,表示模型弹性参数fRPM(m)与观测弹性参数dobs之间的拟合程度,p(m|dobs)为后验概率密度,表示在获得观测弹性参数dobs时物性参数m的分布规律。
假设物性参数m=[φ,Vsh,Sw]之间的相互独立,且服从高斯分布,其均值分别为μφ和/>方差分别为σφ,/>和/>则先验概率密度p(m)可通过以下公式计算获得:
在上式中,p(m)为物性参数m的先验概率密度,n为待估计参数的采样点数,,φ为孔隙度、Sw为饱和度,Vsh为泥质含量,p(φ)为孔隙度的先验概率密度,p(Vsh)为泥质含量的先验概率密度,p(Sw)为饱和度的先验概率密度,μφ为孔隙度均值,为泥质含量均值,/>为饱和度均值,σφ为孔隙度方差,/>为泥质含量方差,/>为饱和度方差,exp为以自然指数e为底的指数函数,n为待估计样点个数,∑为对所有样点求和运算,π为圆周率。
似然函数p(dobs|m)表示模型弹性参数fRPM(m)与观测弹性参数dobs之间的拟合程度,似然函数的值越大表示二者拟合的程度越好,反之则差。假设误差ε满足均值为0,方差为σε的高斯分布,可以得到似然函数p(dobs|m)为:
在上式中,p(dobs|m)为模型弹性参数fRPM(m)与观测弹性参数dobs之间的拟合程度,dobs为观测弹性参数,fRPM(m)为利用岩石物理模型由物性参数计算得到模型弹性参数,σε为误差。确定了先验概率密度p(m)和似然函数p(dobs|m)之后,后验概率密度p(m|dobs),也就是岩石物理反演的反问题的解,利用式(14)就可以求解出来。如果后验概率密度的形式非常简单,直接利用式(14)求解其解析解;但大多数情况下岩石物理模型往往是非线性的,因此后验概率密度形式比较复杂,没有显性的表达式。这个时候可以通过对后验概率密度进行抽样统计,对其进行间接表达。
在本发明一实施例中,上述步骤S103的于贝叶斯理论框架下利用ADMM-MCMC算法迭代所述关系方程中的物性参数计算其对应的后验概率密度数据可包含:利用ADMM算法构造马尔科夫链,通过马尔科夫链替代MCMC算法中的马尔科夫链后,对所述关系方程中的物性参数进行逐个迭代计算,获得对应的后验概率密度数据。
具体的,马尔科夫链蒙特卡罗法MCMC(Markov Chains Monte Carlo)就是一种启发式蒙特卡洛法,被广泛应用于非线性反演中。在MCMC算法中通常利用Metropolis-Hasting算法构建马尔可夫链,在反问题的求解中,目标分布(即平稳分布)一般为似然函数L(x),x为待反演参数,马尔科夫链在某时刻的状态为xt,下一时刻的潜在转移点为x*,根据Metropolis-Hasting算法:
①如果L(x*)≥L(xt),则接受建议,转移到x*,即xt+1=x*
②如果L(x*)<L(xt),根据接受概率进行判断,来决定是否接受x*
在Metropolis-Hasting算法中,马尔科夫链中的所有待反演参数同时迭代。如果待反演参数为多个参数(x,y,z),例如本文中岩石物理反问题中m=[φ,Vsh,Sw],那么需要接收概率同时判断的值,这里假设马尔科夫链在某时刻的状态为(xt,yt,zt),下一时刻的潜在转移点为(x*,y*,z*)。
为了进一步优化MCMC算法,本发明提出了ADMM-MCMC算法,即将交替方向乘子ADMM(Alternating Direction Method of Multipliers)算法引入到MCMC算法中。交替方向乘子ADMM算法是一种解决可分解凸优化问题的方法,主要应用于求解大规模系统的优化问题[31]。ADMM算法可以将原问题的目标函数等价分解成若干个可求解的子问题,然后并行求解每一个子问题,最后协调子问题的解得到原问题的全局解。这里主要利用ADMM方法来用于构造马尔科夫链,不同于Metropolis-Hasting算法中的所有元素同时迭代,这里马尔科夫链元素逐个迭代。
对于需要同时反演多个参数,假设马尔科夫链在某时刻的状态为(xt,yt,zt),下一时刻的潜在转移点为(x*,y*,z*),根据ADMM算法所有元素可以逐个迭代,如下一时刻迭代点的选择按照如下步骤:①首先迭代xt+1,xt+1的接收概率为②然后再迭代yt+1,yt+1的接收概率为/>③再迭代zt+1,zt+1的接收概率为为并且xt+1,yt+1,zt+1的迭代顺序可以任意交换。
使用ADMM算法逐个迭代待反演参数更加符合实际意义,特别是在岩石物理反演中,孔隙度、饱和度和泥质含量都是彼此独立的参数,这样可以尽量保证寻找每个参数的最优解,而不是所有参数组合的最优解。
综上所述,利用式(1)-(12)就可以建立地震弹性参数(纵波速度VP、横波速度VS和密度ρ等)与储层物性参数(孔隙度φ、饱和度Sw和泥质含量Vsh等)的理论关系,即式(13)中的fRPM(m)。再在贝叶斯理论框架下利用ADMM-MCMC算法,从待估计的未知物性参数空间中抽取大量样本并计算其对应的后验概率密度,寻找最大概率密度值并确定其不确定性,就可以得到每个物性参数的估计值。
请参考图3所示,本发明还提供一种岩石物理反演系统,所述系统包含采集模块、构建模块、计算模块和反演模块;所述采集模块用于获取待反演地区的地震弹性参数和储层物理参数;所述构建模块用于通过Gassmann方程和Xu-White模型建立所述地震弹性参数和所述储层物理参数之间的关系方程;所述计算模块用于于贝叶斯理论框架下利用ADMM-MCMC算法迭代所述关系方程中的物性参数计算其对应的后验概率密度数据;所述反演模块用于根据所述后验概率密度数据中最大后验密度值对所述待反演地区进行反演获得反演结果。
在上述岩石物理反演系统中,所述构建模块还可包含分析单元,所述分析单元用于:根据所述地震弹性参数和所述储层物理参数获得饱和流体的体积模量和密度;通过Voigt-Reuss-Hill平均方法计算获得待反演地区的岩石基质的体积模量、剪切模量和密度;根据所述岩石基质的体积模量和剪切模量,通过Xu-White模型获得岩石骨架的体积模量和剪切模量;根据岩石骨架的体积模量和剪切模量与饱和流体的体积模量,通过Gassmann方程获得饱和岩石的体积模量和剪切模量;根据饱和流体的密度、岩石基质的密度、饱和岩石的体积模量和剪切模量计算获得饱和岩石的纵波速度、横波速度和密度。
在上述岩石物理反演系统中,所述计算模块还包含:利用ADMM算法构造马尔科夫链,通过马尔科夫链替代MCMC算法中的马尔科夫链后,对所述关系方程中的物性参数进行逐个迭代计算,获得对应的后验概率密度数据。
本发明还提供一种计算机设备,包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序,所述处理器执行所述计算机程序时实现上述方法。
本发明还提供一种计算机可读存储介质,所述计算机可读存储介质存储有执行上述方法的计算机程序。
本发明的有益技术效果在于:基于Gassmann方程和Xu-White模型岩石物理模型建立地震弹性参数与储层物性参数之间的非线性关系。然后在贝叶斯理论框架下,提出了一种基于ADMM-MCMC算法的求解非线性岩石物理反演问题的方法。利用ADMM算法逐个迭代待反演的物性参数,更加符合岩石物理反演,因为孔隙度、饱和度和泥质含量都是彼此独立的参数,这样可以尽量保证寻找每个参数的最优解,而不是所有参数组合的最优解。再基于对迭代生成的平稳分布马尔科夫链进行统计分析,就可以获取待反演参数最大后验密度,进而获得物性参数的最优解及其不确定性分析。实际测井数据和地震数据应用表明,基于ADMM-MCMC的非线性岩石物理反演方法可以获得较为准确的物性参数。
将基于ADMM-MCMC的非线性岩石物理反演方法应用于实际测井数据;请参考图4所示,A井位于中国东部海上油田含气砂岩储层,测井曲线包括孔隙度、泥质含量和含水饱和度等物性参数,以及纵波速度、横波速度和密度等弹性参数。测井数据经过井震标定后由深度域转换成时间域,并按照2ms进行重采样。主力产气层为2760~2800ms之间的一套厚砂岩储层,其上覆2718ms左右有一层较薄的泥岩层,下覆2822~2852ms之间有一套厚泥岩层,中间夹杂两套较薄的砂岩储层;图4中从左到右为孔隙度数据、粘土含量数据、含水饱和度数据、纵波速度数据、横波速度数据和密度数据。
然后再对测井曲线进行非线性岩石物理反演,如图5所示。即利用图4中纵波速度、横波速度和密度来反演图4中的孔隙度、泥质含量和含水饱和度。选择孔隙度、泥质含量和含水饱和度测井曲线的50次3点平滑结果作为反演的初始值。图5中点虚线代表反演初始值,虚线代表实际测井数据,点线代表反演得到5条马尔科夫链结果,实线代表对5条马尔科夫链结果取平均作为非线性岩石物理反演的结果。由图中可以看到,非线性岩石物理反演结果与实际测井曲线基本一致,虽然幅度没有完全吻合,但整体的变化趋势是相同的。这里给出了非线性岩石物理反演结果与实际测井曲线之间的平均误差,其中孔隙度的平均误差约为16%,含水饱和度的平均误差约为17%,泥质含量的平均误差约为14%。图5中从左到右为孔隙度数据、含水饱和度数据和泥质含量数据。
在对测井曲线进行岩石物理反演之后,再对地震数据进行岩石物理反演。选择一条过A井的二维地震剖面,由AVO叠前地震反演分别获得了纵波速度、横波速度和密度等弹性参数的二维剖面,如图6A至图6C所示;图6A至图6C中箭头所指位置为A井。由二维弹性参数剖面上可以看到,无论是纵波速度、横波速度还是密度在2800ms左右都有一套较厚的低值区域,这是主要的目的层段;图6A至图6C中从上到下分别为纵波速度数据、横波速度数据和密度数据。
利用图6A至图6C中的弹性参数进行非线性岩石物理反演,获取物性参数剖面,如图7A至图7C所示,从上到下分别展示了反演得到的孔隙度、泥质含量和含气饱和度;由反演结果可以看到,主要目的层段反演得到的物性参数与井上的物性参数吻合较好,并且能够获取主要目的层段的物性参数空间展布情况。
如图8所示,该电子设备600还可以包括:通信模块110、输入单元120、音频处理单元130、显示器160、电源170。值得注意的是,电子设备600也并不是必须要包括图8中所示的所有部件;此外,电子设备600还可以包括图8中没有示出的部件,可以参考现有技术。
如图8所示,中央处理器100有时也称为控制器或操作控件,可以包括微处理器或其他处理器装置和/或逻辑装置,该中央处理器100接收输入并控制电子设备600的各个部件的操作。
其中,存储器140,例如可以是缓存器、闪存、硬驱、可移动介质、易失性存储器、非易失性存储器或其它合适装置中的一种或更多种。可储存上述与失败有关的信息,此外还可存储执行有关信息的程序。并且中央处理器100可执行该存储器140存储的该程序,以实现信息存储或处理等。
输入单元120向中央处理器100提供输入。该输入单元120例如为按键或触摸输入装置。电源170用于向电子设备600提供电力。显示器160用于进行图像和文字等显示对象的显示。该显示器例如可为LCD显示器,但并不限于此。
该存储器140可以是固态存储器,例如,只读存储器(ROM)、随机存取存储器(RAM)、SIM卡等。还可以是这样的存储器,其即使在断电时也保存信息,可被选择性地擦除且设有更多数据,该存储器的示例有时被称为EPROM等。存储器140还可以是某种其它类型的装置。存储器140包括缓冲存储器141(有时被称为缓冲器)。存储器140可以包括应用/功能存储部142,该应用/功能存储部142用于存储应用程序和功能程序或用于通过中央处理器100执行电子设备600的操作的流程。
存储器140还可以包括数据存储部143,该数据存储部143用于存储数据,例如联系人、数字数据、图片、声音和/或任何其他由电子设备使用的数据。存储器140的驱动程序存储部144可以包括电子设备的用于通信功能和/或用于执行电子设备的其他功能(如消息传送应用、通讯录应用等)的各种驱动程序。
通信模块110即为经由天线111发送和接收信号的发送机/接收机110。通信模块(发送机/接收机)110耦合到中央处理器100,以提供输入信号和接收输出信号,这可以和常规移动通信终端的情况相同。
基于不同的通信技术,在同一电子设备中,可以设置有多个通信模块110,如蜂窝网络模块、蓝牙模块和/或无线局域网模块等。通信模块(发送机/接收机)110还经由音频处理器130耦合到扬声器131和麦克风132,以经由扬声器131提供音频输出,并接收来自麦克风132的音频输入,从而实现通常的电信功能。音频处理器130可以包括任何合适的缓冲器、解码器、放大器等。另外,音频处理器130还耦合到中央处理器100,从而使得可以通过麦克风132能够在本机上录音,且使得可以通过扬声器131来播放本机上存储的声音。
本领域内的技术人员应明白,本发明的实施例可提供为方法、系统、或计算机程序产品。因此,本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且,本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器,使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中,使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品,该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上,使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理,从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。
以上所述的具体实施例,对本发明的目的、技术方案和有益效果进行了进一步详细说明,所应理解的是,以上所述仅为本发明的具体实施例而已,并不用于限定本发明的保护范围,凡在本发明的精神和原则之内,所做的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。

Claims (6)

1.一种岩石物理反演方法,其特征在于,所述方法包含:
获取待反演地区的地震弹性参数和储层物理参数;
通过Gassmann方程和Xu-White模型建立所述地震弹性参数和所述储层物理参数之间的关系方程;
于贝叶斯理论框架下利用ADMM-MCMC算法迭代所述关系方程中的物性参数,计算物性参数对应的后验概率密度数据;
根据所述后验概率密度数据中最大后验密度值对所述待反演地区进行反演获得反演结果;
通过Gassmann方程和Xu-White模型建立所述地震弹性参数和所述储层物理参数之间的关系方程包含:
根据所述地震弹性参数和所述储层物理参数获得饱和流体的体积模量和密度;
通过Voigt-Reuss-Hill平均方法计算获得待反演地区的岩石基质的体积模量、剪切模量和密度;
根据所述岩石基质的体积模量和剪切模量,通过Xu-White模型获得岩石骨架的体积模量和剪切模量;
根据岩石骨架的体积模量和剪切模量与饱和流体的体积模量,通过Gassmann方程获得饱和岩石的体积模量和剪切模量;
根据饱和流体的密度、岩石基质的密度、饱和岩石的体积模量和剪切模量计算获得饱和岩石的纵波速度、横波速度和密度;
于贝叶斯理论框架下利用ADMM-MCMC算法迭代所述关系方程中的物性参数,计算其对应的后验概率密度数据包含:
通过饱和岩石的纵波速度、横波速度和密度构建观测弹性参数,通过以下关系方程计算后验概率密度数据:
在上式中,m为物性参数,p(m)为物性参数m的先验概率密度,p(dobs|m)为模型弹性参数fRPM(m)与观测弹性参数dobs之间的拟合程度,p(m|dobs)为后验概率密度值,fRPM(m)为利用岩石物理模型由物性参数计算得到的模型弹性参数,dobs为观测弹性参数;
于贝叶斯理论框架下利用ADMM-MCMC算法迭代所述关系方程中的物性参数,计算其对应的后验概率密度数据还包含:通过用ADMM算法构造的马尔科夫链替代MCMC算法中的马尔科夫链后,对所述关系方程中的物性参数进行逐个迭代计算,获得对应的后验概率密度数据。
2.根据权利要求1所述的岩石物理反演方法,其特征在于,所述先验概率密度通过以下公式计算获得:
在上式中,p(m)为物性参数m的先验概率密度,m=[f,Vsh,Sw],f为孔隙度、Sw为饱和度,Vsh为泥质含量,p(f)为孔隙度的先验概率密度,p(Vsh)为泥质含量的先验概率密度,p(Sw)为饱和度的先验概率密度,mf为孔隙度均值,mVsh为泥质含量均值,mSw为饱和度均值,sf为孔隙度方差,为泥质含量方差,sSw为饱和度方差,exp为以自然指数e为底的指数函数,n为待估计样点个数,∑为对所有样点求和运算,p为圆周率。
3.根据权利要求1所述的岩石物理反演方法,其特征在于,所述拟合程度通过以下公式计算获得:
在上式中,p(dobs|m)为模型弹性参数fRPM(m)与观测弹性参数dobs之间的拟合程度,dobs为观测弹性参数,fRPM(m)为利用岩石物理模型由物性参数计算得到的模型弹性参数,sε为误差。
4.一种适用于权利要求1至3中任一项所述的岩石物理反演方法的岩石物理反演系统,其特征在于,所述系统包含采集模块、构建模块、计算模块和反演模块;
所述采集模块用于获取待反演地区的地震弹性参数和储层物理参数;
所述构建模块用于通过Gassmann方程和Xu-White模型建立所述地震弹性参数和所述储层物理参数之间的关系方程;
所述计算模块用于于贝叶斯理论框架下利用ADMM-MCMC算法迭代所述关系方程中的物性参数,计算其对应的后验概率密度数据;
所述反演模块用于根据所述后验概率密度数据中最大后验密度值对所述待反演地区进行反演获得反演结果;
所述构建模块还包含分析单元,所述分析单元用于:
根据所述地震弹性参数和所述储层物理参数获得饱和流体的体积模量和密度;
通过Voigt-Reuss-Hill平均方法计算获得待反演地区的岩石基质的体积模量、剪切模量和密度;
根据所述岩石基质的体积模量和剪切模量,通过Xu-White模型获得岩石骨架的体积模量和剪切模量;
根据岩石骨架的体积模量和剪切模量与饱和流体的体积模量,通过Gassmann方程获得饱和岩石的体积模量和剪切模量;
根据饱和流体的密度、岩石基质的密度、饱和岩石的体积模量和剪切模量计算获得饱和岩石的纵波速度、横波速度和密度;
所述计算模块还包含:通过用ADMM算法构造的马尔科夫链替代MCMC算法中的马尔科夫链后,对所述关系方程中的物性参数进行逐个迭代计算,获得对应的后验概率密度数据。
5.一种计算机设备,包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序,其特征在于,所述处理器执行所述计算机程序时实现权利要求1至3任一所述方法。
6.一种计算机可读存储介质,其特征在于,所述计算机可读存储介质存储有执行权利要求1至3任一所述方法的计算机程序。
CN202010207180.6A 2020-03-23 2020-03-23 岩石物理反演方法及系统 Active CN113433587B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010207180.6A CN113433587B (zh) 2020-03-23 2020-03-23 岩石物理反演方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010207180.6A CN113433587B (zh) 2020-03-23 2020-03-23 岩石物理反演方法及系统

Publications (2)

Publication Number Publication Date
CN113433587A CN113433587A (zh) 2021-09-24
CN113433587B true CN113433587B (zh) 2023-09-26

Family

ID=77752382

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010207180.6A Active CN113433587B (zh) 2020-03-23 2020-03-23 岩石物理反演方法及系统

Country Status (1)

Country Link
CN (1) CN113433587B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115586573B (zh) * 2022-09-15 2023-06-09 河海大学 一种致密砂岩储层的动态约束物性参数地震反演方法
CN116165701B (zh) * 2022-12-08 2024-02-20 中国矿业大学 一种碳酸盐岩储层的孔隙结构参数叠前地震直接反演方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104516017A (zh) * 2013-09-29 2015-04-15 中国石油化工股份有限公司 一种碳酸盐岩岩石物理参数地震反演方法
CN106168676A (zh) * 2015-05-22 2016-11-30 中国石油化工股份有限公司 基于地震资料的地层岩性和流体识别方法和装置
CN107290782A (zh) * 2016-03-30 2017-10-24 中国石油化工股份有限公司 储层孔隙度、含水饱和度和泥质含量参数同时反演新方法
CN110133720A (zh) * 2019-06-04 2019-08-16 南京信息工程大学 一种基于统计岩石物理模型的横波速度预测方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7286939B2 (en) * 2003-10-28 2007-10-23 Westerngeco, L.L.C. Method for estimating porosity and saturation in a subsurface reservoir

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104516017A (zh) * 2013-09-29 2015-04-15 中国石油化工股份有限公司 一种碳酸盐岩岩石物理参数地震反演方法
CN106168676A (zh) * 2015-05-22 2016-11-30 中国石油化工股份有限公司 基于地震资料的地层岩性和流体识别方法和装置
CN107290782A (zh) * 2016-03-30 2017-10-24 中国石油化工股份有限公司 储层孔隙度、含水饱和度和泥质含量参数同时反演新方法
CN110133720A (zh) * 2019-06-04 2019-08-16 南京信息工程大学 一种基于统计岩石物理模型的横波速度预测方法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
Multimodal Markov chain Monte Carlo method for nonlinear petrophysical seismic inversion;Leandro Passos de Figueiredo et al.;《Geophysics》;第84卷(第5期);第M1-M13页 *
Rock physics inversion based on an optimized MCMC method;Zhang Jia-Jia et al.;《APPLIED GEOPHYSICS》;第18卷;第288-298页 *
基于贝叶斯分类的储层物性参数联合反演方法;胡华锋 等;《石油物探》;第51卷(第3期);第225-232页 *
岩石物理驱动的储层裂缝参数与物性参数概率地震反演方法;潘新朋 等;《地球物理学报》;第61卷(第2期);第683-696页 *
贝叶斯网络结构稀疏学习研究进展;郭珉 等;《模式识别与人工智能》;第29卷(第10期);第907-723页 *

Also Published As

Publication number Publication date
CN113433587A (zh) 2021-09-24

Similar Documents

Publication Publication Date Title
CN110031896B (zh) 基于多点地质统计学先验信息的地震随机反演方法及装置
CA2580570C (en) Integrated anisotropic rock physics model
Jullum et al. A Gaussian-based framework for local Bayesian inversion of geophysical data to rock properties
CN109490963B (zh) 裂缝储层岩石物理建模方法及系统
US11010969B1 (en) Generation of subsurface representations using layer-space
US10984590B1 (en) Generation of subsurface representations using layer-space
EP2856373A1 (en) System and method for predicting rock strength
Eikrem et al. Iterated extended Kalman filter method for time‐lapse seismic full‐waveform inversion
Karimpouli et al. Estimating 3D elastic moduli of rock from 2D thin-section images using differential effective medium theory
CN113433587B (zh) 岩石物理反演方法及系统
CN110887772B (zh) 一种碳酸盐岩储层渗透率识别方法、系统及装置
CN109407150A (zh) 基于统计岩石物理的页岩储层可压裂性解释方法及系统
US20220091290A1 (en) Systems and methods for generating elastic property data as a function of position and time in a subsurface volume of interest
CN107797139A (zh) 页岩储层游离气含气量地震预测方法及系统
US20220091300A1 (en) Systems and methods for generating subsurface property data as a function of position and time in a subsurface volume of interest
US20220091291A1 (en) Systems and methods for generating subsurface data as a function of position and time in a subsurface volume of interest
US20210341642A1 (en) Nested model simulations to generate subsurface representations
Zhang et al. Bayesian variational time-lapse full waveform inversion
CN114460639A (zh) 页岩油储层渗透率的预测方法及装置
Pyrcz et al. Representative input parameters for geostatistical simulation
CN113671572B (zh) 基于室内砂箱的地震数据成像方法及装置
CN112649867B (zh) 虚拟井构建方法及系统
CN112462421A (zh) 储层信息预测方法、装置、电子设备及存储介质
Carrasquero et al. A Seismically-Driven Integrated Reservoir Characterisation for Modelling Purposes: A Green Field Case Study
Bakshevskaya et al. Methods of modeling hydraulic heterogeneity of sedimentary formations

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