CN112147683B - 基于岩石物理关系约束的叠前稀疏层反演方法及系统 - Google Patents

基于岩石物理关系约束的叠前稀疏层反演方法及系统 Download PDF

Info

Publication number
CN112147683B
CN112147683B CN201910571036.8A CN201910571036A CN112147683B CN 112147683 B CN112147683 B CN 112147683B CN 201910571036 A CN201910571036 A CN 201910571036A CN 112147683 B CN112147683 B CN 112147683B
Authority
CN
China
Prior art keywords
reflectivity
wave velocity
inversion
matrix
establishing
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
CN201910571036.8A
Other languages
English (en)
Other versions
CN112147683A (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.)
China Petroleum and Chemical Corp
Sinopec Exploration and Production Research Institute
Original Assignee
China Petroleum and Chemical Corp
Sinopec Exploration and Production Research Institute
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 China Petroleum and Chemical Corp, Sinopec Exploration and Production Research Institute filed Critical China Petroleum and Chemical Corp
Priority to CN201910571036.8A priority Critical patent/CN112147683B/zh
Publication of CN112147683A publication Critical patent/CN112147683A/zh
Application granted granted Critical
Publication of CN112147683B publication Critical patent/CN112147683B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/62—Physical property of subsurface
    • G01V2210/624—Reservoir parameters

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

公开了一种基于岩石物理关系约束的叠前稀疏层反演方法及系统。该方法可以包括：步骤1：读取角度叠加数据和角度子波，建立角度子波褶积矩阵；步骤2：根据测井数据，建立岩石物理关系；步骤3：建立基于岩石物理关系约束的AVO近似公式；步骤4：建立AVO正演方程组；步骤5：构建反射系数奇偶分解矩阵；步骤6：根据AVO正演方程组与反射系数奇偶分解矩阵，建立反演目标函数；步骤7：计算弹性参数反演结果并进行低频补偿，获得最终的弹性参数反演结果。本发明通过计算纵波速度反射率、纵横波速度比反射率以及密度反射率，结合岩石物理关系约束与反射系数奇偶分解理论，实现了关于纵横波速度比的稀疏层反演。

Description

基于岩石物理关系约束的叠前稀疏层反演方法及系统
技术领域
本发明涉及油气勘探领域，更具体地，涉及一种基于岩石物理关系约束的叠前稀疏层反演方法及系统。
背景技术
叠前地震数据中包含了丰富的偏移距信息，振幅随偏移距/入射角的变化揭示了地下介质的岩性变化和孔隙内流体成分变化。因此利用叠前AVA同步反演可以从角度部分叠加地震数据中提取多种岩石弹性参数，其中纵横波速度比对储层岩性及孔隙内流体的变化更为敏感，是地震解释人员较常应用的“烃类指示因子”。然而常规叠前AVA同步反演无法直接获取纵横波速度比，弹性参数的间接转换通常会引入累计误差，另外纵横波速度比相对纵波阻抗\纵波速度有更大的不确定性，其反演精度对道集质量以及入射角范围也更为敏感。常规叠前AVA同步反演，其本质还是稀疏脉冲反演，因此垂向分辨率较低。
前人对叠前地震反演问题进行了深入的研究。叠前AVA反演最早可追溯至Smith提出的加权叠加方法，该方法属于带限反演，未考虑到子波的带限效应，并且反演结果依然是弹性参数反射率。将贝叶斯理论引入到叠前地震反演中，通过假设地震数据的似然函数和模型参数的先验分布均服从多变量高斯分布，给出了模型参数后验分布的均值和方差的解析解，并指出了确定性反演的解为模型参数后验分布的期望；而随机反演的解，可以通过MCMC等技术从后验概率中抽样实现。提出“长尾巴”分布相对高斯分布可以进一步提高反演的垂向分辨率，并借助于参数协方差矩阵的去相关技术去除三参数的相关性，提高了叠前反演的不适定性，针对长尾巴分布中的Lp范数分布、单变量柯西分布和Huber分布分别实现了叠前二项和三项反演。通过引入岩石物理经验公式作为约束，进一步提高了叠前反演的稳定性，实现了基于角道集的叠前三参数同步反演算法，该算法为HRS软件的叠前反演模块的核心算法。进一步发展贝叶斯叠前反演，通过引入岩石物理关系约束以及点约束，提高了反演结果的鲁棒性。为了进一步提高反演结果的垂向分辨率，提出基于基追踪优化算法的稀疏层反演，并将稀疏层反演从叠后反演扩展至叠前反演，然而该方法并未考虑到由弹性参数相关而引入的不适定性。上述叠前反演算法并未考虑到由弹性参数间接转换而造成纵横波速度比存在累计误差，以及纵横波速度比反演的垂向分辨率较低等问题。因此，有必要开发一种基于岩石物理关系约束的叠前稀疏层反演方法及系统。
公开于本发明背景技术部分的信息仅仅旨在加深对本发明的一般背景技术的理解，而不应当被视为承认或以任何形式暗示该信息构成已为本领域技术人员所公知的现有技术。
发明内容
本发明提出了一种基于岩石物理关系约束的叠前稀疏层反演方法及系统，其能够通过计算纵波速度反射率、纵横波速度比反射率以及密度反射率，结合岩石物理关系约束与反射系数奇偶分解理论，实现了关于纵横波速度比的稀疏层反演。
根据本发明的一方面，提出了一种基于岩石物理关系约束的叠前稀疏层反演方法。所述方法可以包括：步骤1：读取角度叠加数据和角度子波，并根据所述角度子波建立角度子波褶积矩阵；步骤2：根据测井数据，建立岩石物理关系；步骤3：根据所述岩石物理关系与Aki-Richards近似式，建立基于岩石物理关系约束的AVO近似公式；步骤4：根据所述角度叠加数据、所述角度子波褶积矩阵以及所述AVO近似公式，建立AVO正演方程组；步骤5：确定目的层段的薄层最大时间厚度，根据反射系数奇偶分解理论，构建反射系数奇偶分解矩阵；步骤6：根据所述AVO正演方程组与所述反射系数奇偶分解矩阵，并在最小化弹性参数奇偶分量系数的1范数的约束下，建立反演目标函数；步骤7：根据所述反演目标函数，计算弹性参数反演结果并进行低频补偿，获得最终的弹性参数反演结果。
优选地，所述步骤2包括：步骤21：根据测井数据，通过拟合纵波速度曲线与横波速度曲线的线性关系，获得Castagna经验公式的斜率Cαβ；步骤22：根据测井数据，通过拟合纵波速度曲线与密度曲线的指数关系，获得Gardnar经验公式的指数Cαρ；步骤23：根据所述Castagna经验公式的斜率，并结合弹性参数反射率之间的转换关系，利用公式(1)建立纵波速度反射率与纵横波速度比反射率之间的岩石物理关系：
Figure BDA0002110886280000031
其中，Rα表示纵波速度反射率，Rγ表示纵横波速度比反射率，ΔRα表示纵波速度反射率的扰动量，
Figure BDA0002110886280000032
表示背景纵横波速度比；步骤24：根据所Gardnar经验公式的指数与所述Castagna经验公式的斜率，并结合弹性参数反射率之间的转换关系，利用公式(2)建立密度反射率与纵横波速度比反射率之间的岩石物理关系：
Figure BDA0002110886280000033
其中，Rρ表示纵波速度反射率，ΔRρ表示纵波速度反射率的扰动量。
优选地，所述AVO近似公式为：
Rpp(θ)＝A(θ)Rγ+B(θ)ΔRα+C(θ)ΔRρ (3)
其中，
Figure BDA0002110886280000034
Figure BDA0002110886280000041
Rpp(θ)表示入射角为θ的角度反射系数。
优选地，所述AVO正演方程组为：
d＝Gr (4)
其中，
Figure BDA0002110886280000042
d(θi)表示入射角为θi的角度叠加数据，其中i＝1,2,…K，K表示入射角的个数；
Figure BDA0002110886280000043
Rγ、ΔRα和ΔRρ分别表示纵横波速度比反射率向量、纵波速度反射率扰动量向量和密度反射率扰动量向量；
Figure BDA0002110886280000044
W(θi)表示入射角为θi的角度子波褶积矩阵，其中i＝1,2,…K；A(θi)、B(θi)、C(θi)均为对角阵，其中i＝1,2,…K；其对角阵元素为：
Figure BDA0002110886280000045
其中，j＝1,2,…N，N表示反射系数的采样点个数。
优选地，所述反射系数奇偶分解矩阵为：
D＝[De Do] (5)
其中，De表示反射系数偶分量矩阵，
Figure BDA0002110886280000051
Do表示反射系数奇分量矩阵，
Figure BDA0002110886280000052
其中，M表示薄层最大时间厚度对应的采样点个数。
优选地，所述反演目标函数为：
Figure BDA0002110886280000053
其中，x表示待反演参数，
Figure BDA0002110886280000054
xγ、xα和xρ分别表示纵横波速度比反射率的奇偶分量系数、纵波速度反射率扰动量的奇偶分量系数以及密度反射率扰动量的奇偶分量系数；
Figure BDA0002110886280000055
λ表示稀疏约束权重；上标T表示矩阵转置。
优选地，所述步骤7包括：步骤71：根据所述反演目标函数，利用迭代重加权最小二乘算法进行求解，计算待反演参数；步骤72：根据所述待反演参数，计算弹性参数反射率；步骤73：根据弹性参数反射率，计算弹性参数，并进行低频补偿。
优选地，所述步骤71中的所述迭代重加权最小二乘算法的计算步骤为：步骤711：设置待反演参数x的初始值为零向量；步骤712：计算非一致性加权矩阵Q，Q为对角矩阵，其对角线元素为：
Figure BDA0002110886280000061
其中，ε表示稳定系数，xl表示x的第l个元素，l＝1,2,...,3N；步骤713：通过结合Cholesky分解和上下三角矩阵分解，求解线性方程组
Figure BDA0002110886280000062
得到待反演参数x；步骤714：重复步骤712-步骤713，直至达到规定的迭代次数，输出待反演参数x。
根据本发明的另一方面，提出了一种基于岩石物理关系约束的叠前稀疏层反演系统，其特征在于，该系统包括：存储器，存储有计算机可执行指令；处理器，所述处理器运行所述存储器中的计算机可执行指令，执行以下步骤：步骤1：读取角度叠加数据和角度子波，并根据所述角度子波建立角度子波褶积矩阵；步骤2：根据测井数据，建立岩石物理关系；步骤3：根据所述岩石物理关系与Aki-Richards近似式，建立基于岩石物理关系约束的AVO近似公式；步骤4：根据所述角度叠加数据、所述角度子波褶积矩阵以及所述AVO近似公式，建立AVO正演方程组；步骤5：确定目的层段的薄层最大时间厚度，根据反射系数奇偶分解理论，构建反射系数奇偶分解矩阵；步骤6：根据所述AVO正演方程组与所述反射系数奇偶分解矩阵，并在最小化弹性参数奇偶分量系数的1范数的约束下，建立反演目标函数；步骤7：根据所述反演目标函数，计算弹性参数反演结果并进行低频补偿，获得最终的弹性参数反演结果。
优选地，所述步骤2包括：步骤21：根据测井数据，通过拟合纵波速度曲线与横波速度曲线的线性关系，获得Castagna经验公式的斜率Cαβ；步骤22：根据测井数据，通过拟合纵波速度曲线与密度曲线的指数关系，获得Gardnar经验公式的指数Cαρ；步骤23：根据所述Castagna经验公式的斜率，并结合弹性参数反射率之间的转换关系，利用公式(1)建立纵波速度反射率与纵横波速度比反射率之间的岩石物理关系：
Figure BDA0002110886280000071
其中，Rα表示纵波速度反射率，Rγ表示纵横波速度比反射率，ΔRα表示纵波速度反射率的扰动量，
Figure BDA0002110886280000072
表示背景纵横波速度比；步骤24：根据所Gardnar经验公式的指数与所述Castagna经验公式的斜率，并结合弹性参数反射率之间的转换关系，利用公式(2)建立密度反射率与纵横波速度比反射率之间的岩石物理关系：
Figure BDA0002110886280000073
其中，Rρ表示纵波速度反射率，ΔRρ表示纵波速度反射率的扰动量。
优选地，所述AVO近似公式为：
Rpp(θ)＝A(θ)Rγ+B(θ)ΔRα+C(θ)ΔRρ (3)
其中，
Figure BDA0002110886280000074
Figure BDA0002110886280000075
Rpp(θ)表示入射角为θ的角度反射系数。
优选地，所述AVO正演方程组为：
d＝Gr (4)
其中，
Figure BDA0002110886280000076
d(θi)表示入射角为θi的角度叠加数据，其中i＝1,2,…K，K表示入射角的个数；
Figure BDA0002110886280000081
Rγ、ΔRα和ΔRρ分别表示纵横波速度比反射率向量、纵波速度反射率扰动量向量和密度反射率扰动量向量；
Figure BDA0002110886280000082
W(θi)表示入射角为θi的角度子波褶积矩阵，其中i＝1,2,…K；A(θi)、B(θi)、C(θi)均为对角阵，其中i＝1,2,…K；其对角阵元素为：
Figure BDA0002110886280000083
其中，j＝1,2,…N，N表示反射系数的采样点个数。
优选地，所述反射系数奇偶分解矩阵为：
D＝[De Do] (5)
其中，De表示反射系数偶分量矩阵，
Figure BDA0002110886280000084
Do表示反射系数奇分量矩阵，
Figure BDA0002110886280000091
其中，M表示薄层最大时间厚度对应的采样点个数。
优选地，所述反演目标函数为：
Figure BDA0002110886280000092
其中，x表示待反演参数，
Figure BDA0002110886280000093
xγ、xα和xρ分别表示纵横波速度比反射率的奇偶分量系数、纵波速度反射率扰动量的奇偶分量系数以及密度反射率扰动量的奇偶分量系数；
Figure BDA0002110886280000094
λ表示稀疏约束权重；上标T表示矩阵转置。
优选地，所述步骤7包括：步骤71：根据所述反演目标函数，利用迭代重加权最小二乘算法进行求解，计算待反演参数；步骤72：根据所述待反演参数，计算弹性参数反射率；步骤73：根据弹性参数反射率，计算弹性参数，并进行低频补偿。
优选地，所述步骤71中的所述迭代重加权最小二乘算法的计算步骤为：步骤711：设置待反演参数x的初始值为零向量；步骤712：计算非一致性加权矩阵Q，Q为对角矩阵，其对角线元素为：
Figure BDA0002110886280000095
其中，ε表示稳定系数，xl表示x的第l个元素，l＝1,2,...,3N；步骤713：通过结合Cholesky分解和上下三角矩阵分解，求解线性方程组
Figure BDA0002110886280000101
得到待反演参数x；步骤714：重复步骤712-步骤713，直至达到规定的迭代次数，输出待反演参数x。
本发明的方法和装置具有其它的特性和优点，这些特性和优点从并入本文中的附图和随后的具体实施方式中将是显而易见的，或者将在并入本文中的附图和随后的具体实施方式中进行详细陈述，这些附图和具体实施方式共同用于解释本发明的特定原理。
附图说明
通过结合附图对本发明示例性实施例进行更详细的描述，本发明的上述以及其它目的、特征和优势将变得更加明显，其中，在本发明示例性实施例中，相同的参考标号通常代表相同部件。
图1示出了根据本发明的基于岩石物理关系约束的叠前稀疏层反演方法的步骤的流程图。
图2a、图2b、图2c分别示出了根据本发明的一个实施例的中心角分别8度、18度、28度的角度部分叠加数据的示意图。
图3示出了根据图2a、图2b、图2c的中心角分别8度、18度、28度的角度子波的示意图。
图4a示出了根据本发明的一个实施例的纵波速度与横波速度的交会图及其拟合关系的示意图。
图4b示出了根据本发明的一个实施例的纵波速度与密度的交会图及其拟合关系的示意图。
图5a、图5b、图5c分别示出了根据本发明的一个实施例的纵横波速度比反演剖面、纵波速度反演剖面、密度反演剖面的示意图。
图6a、图6b、图6c分别示出了根据本发明的一个实施例的井旁道的纵横波速度比反演结果、纵波速度反演结果、密度反演结果的示意图。
具体实施方式
下面将参照附图更详细地描述本发明。虽然附图中显示了本发明的优选实施例，然而应该理解，可以以各种形式实现本发明而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了使本发明更加透彻和完整，并且能够将本发明的范围完整地传达给本领域的技术人员。
图1示出了根据本发明的基于岩石物理关系约束的叠前稀疏层反演方法的步骤的流程图。
在该实施例中，根据本发明的基于岩石物理关系约束的叠前稀疏层反演方法可以包括：步骤1：读取角度叠加数据和角度子波，并根据角度子波建立角度子波褶积矩阵；步骤2：根据测井数据，建立岩石物理关系；步骤3：根据岩石物理关系与Aki-Richards近似式，建立基于岩石物理关系约束的AVO近似公式；步骤4：根据角度叠加数据、角度子波褶积矩阵以及AVO近似公式，建立AVO正演方程组；步骤5：确定目的层段的薄层最大时间厚度，根据反射系数奇偶分解理论，构建反射系数奇偶分解矩阵；步骤6：根据AVO正演方程组与反射系数奇偶分解矩阵，并在最小化弹性参数奇偶分量系数的1范数的约束下，建立反演目标函数；步骤7：根据反演目标函数，计算弹性参数反演结果并进行低频补偿，获得最终的弹性参数反演结果。
在一个示例中，步骤2包括：步骤21：根据测井数据，通过拟合纵波速度曲线与横波速度曲线的线性关系，获得Castagna经验公式的斜率Cαβ；步骤22：根据测井数据，通过拟合纵波速度曲线与密度曲线的指数关系，获得Gardnar经验公式的指数Cαρ；步骤23：根据Castagna经验公式的斜率，并结合弹性参数反射率之间的转换关系，利用公式(1)建立纵波速度反射率与纵横波速度比反射率之间的岩石物理关系：
Figure BDA0002110886280000121
其中，Rα表示纵波速度反射率，Rγ表示纵横波速度比反射率，ΔRα表示纵波速度反射率的扰动量，
Figure BDA0002110886280000122
表示背景纵横波速度比；步骤24：根据所Gardnar经验公式的指数与Castagna经验公式的斜率，并结合弹性参数反射率之间的转换关系，利用公式(2)建立密度反射率与纵横波速度比反射率之间的岩石物理关系：
Figure BDA0002110886280000123
其中，Rρ表示纵波速度反射率，ΔRρ表示纵波速度反射率的扰动量。
在一个示例中，AVO近似公式为：
Rpp(θ)＝A(θ)Rγ+B(θ)ΔRα+C(θ)ΔRρ (3)
其中，
Figure BDA0002110886280000124
Figure BDA0002110886280000125
Rpp(θ)表示入射角为θ的角度反射系数。
在一个示例中，AVO正演方程组为：
d＝Gr (4)
其中，
Figure BDA0002110886280000126
d(θi)表示入射角为θi的角度叠加数据，其中i＝1,2,…K，K表示入射角的个数；
Figure BDA0002110886280000127
Rγ、ΔRα和ΔRρ分别表示纵横波速度比反射率向量、纵波速度反射率扰动量向量和密度反射率扰动量向量；
Figure BDA0002110886280000131
W(θi)表示入射角为θi的角度子波褶积矩阵，其中i＝1,2,…K；A(θi)、B(θi)、C(θi)均为对角阵，其中i＝1,2,…K；其对角阵元素为：
Figure BDA0002110886280000132
其中，j＝1,2,…N，N表示反射系数的采样点个数。
在一个示例中，反射系数奇偶分解矩阵为：
D＝[De Do] (5)
其中，De表示反射系数偶分量矩阵，
Figure BDA0002110886280000133
Do表示反射系数奇分量矩阵，
Figure BDA0002110886280000141
其中，M表示薄层最大时间厚度对应的采样点个数。
在一个示例中，反演目标函数为：
Figure BDA0002110886280000142
其中，x表示待反演参数，
Figure BDA0002110886280000143
xγ、xα和xρ分别表示纵横波速度比反射率的奇偶分量系数、纵波速度反射率扰动量的奇偶分量系数以及密度反射率扰动量的奇偶分量系数；
Figure BDA0002110886280000144
λ表示稀疏约束权重；上标T表示矩阵转置。
在一个示例中，步骤7包括：步骤71：根据反演目标函数，利用迭代重加权最小二乘算法进行求解，计算待反演参数；步骤72：根据待反演参数，计算弹性参数反射率；步骤73：根据弹性参数反射率，计算弹性参数，并进行低频补偿。
在一个示例中，步骤71中的迭代重加权最小二乘算法的计算步骤为：步骤711：设置待反演参数x的初始值为零向量；步骤712：计算非一致性加权矩阵Q，Q为对角矩阵，其对角线元素为：
Figure BDA0002110886280000145
其中，ε表示稳定系数，是一个很小的数，一般可取10-9，xl表示x的第l个元素，l＝1,2,...,3N；步骤713：通过结合Cholesky分解和上下三角矩阵分解，求解线性方程组
Figure BDA0002110886280000151
得到待反演参数x；步骤714：重复步骤712-步骤713，直至达到规定的迭代次数，输出待反演参数x。
具体地，根据本发明的基于岩石物理关系约束的叠前稀疏层反演方法可以包括：
步骤1：读取角度叠加数据和角度子波，并根据角度子波建立角度子波褶积矩阵。
步骤2包括：步骤21：根据测井数据，通过拟合纵波速度曲线与横波速度曲线的线性关系，获得Castagna经验公式的斜率Cαβ，其中，Castagna经验公式为：
β≈Cαβα+Dαβ (7)
其中，β表示横波速度，α表示纵波速度，Dαβ表示Castagna公式的截距；步骤22：根据测井数据，通过拟合纵波速度曲线与密度曲线的指数关系，获得Gardnar经验公式的指数Cαρ，Gardnar经验公式为：
Figure BDA0002110886280000152
其中，ρ表示密度，Dαρ表示Gardnar公式的系数；步骤23：对Castagna经验公式左右两边进行微分，并同时除以β得到：
Figure BDA0002110886280000153
考虑到纵横波速度比的计算公式为
Figure BDA0002110886280000154
对该公式左右两边进行微分，并同时除以γ得到：
Figure BDA0002110886280000155
综合公式(9)和公式(10)，得到：
Figure BDA0002110886280000161
其中，
Figure BDA0002110886280000162
表示纵波速度反射率，可以用Rα表示，
Figure BDA0002110886280000163
表示纵横波速度比反射率，可以用Rγ表示，同时引入纵波速度反射率的扰动量ΔRα，将约等式修正为等式，即得到公式(1)为纵波速度反射率与纵横波速度比反射率之间的岩石物理关系；步骤24：对Gardnar经验公式左右两边进行微分，可以得到：
Figure BDA0002110886280000166
将公式(12)代入公式(11)可以得到：
Figure BDA0002110886280000164
通过引入密度反射率的扰动量ΔRρ，将约等式修正为等式，即得到公式(2)为密度反射率与纵横波速度比反射率之间的岩石物理关系。
步骤3：利用Aki-Richards近似式可以有效描述角度反射系数与弹性参数反射率的关系，Aki-Richards近似式为：
Figure BDA0002110886280000165
将公式(1)与公式(2)的岩石物理关系式，代入公式(14)即可得到基于岩石物理关系约束的AVO近似公式为公式(3)。
步骤4：根据角度叠加数据、角度子波褶积矩阵以及AVO近似公式，建立AVO正演方程组为公式(4)。
步骤5：确定目的层段的薄层最大时间厚度，根据反射系数奇偶分解理论，构建反射系数奇偶分解矩阵为公式(5)。
步骤6：根据AVO正演方程组与反射系数奇偶分解矩阵，并在最小化弹性参数奇偶分量系数的1范数的约束下，建立反演目标函数为公式(6)。
步骤7包括：步骤71：根据反演目标函数，利用迭代重加权最小二乘算法进行求解，计算待反演参数，其中，步骤71中的迭代重加权最小二乘算法的计算步骤为：步骤711：设置待反演参数x的初始值为零向量；步骤712：计算非一致性加权矩阵Q；步骤713：通过结合Cholesky分解和上下三角矩阵分解，求解线性方程组
Figure BDA0002110886280000171
得到待反演参数x；步骤714：重复步骤712-步骤713，直至达到规定的迭代次数，输出待反演参数x；步骤72：根据待反演参数，通过公式(15)计算弹性参数反射率：
Figure BDA0002110886280000172
其中，I为N行N列的单位矩阵，Γ和Λ均为对角矩阵，其对角线元素分别为：
Figure BDA0002110886280000173
步骤73：根据弹性参数反射率，计算弹性参数，并进行低频补偿。
本方法通过计算纵波速度反射率、纵横波速度比反射率以及密度反射率，结合岩石物理关系约束与反射系数奇偶分解理论，实现了关于纵横波速度比的稀疏层反演。
应用示例
为便于理解本发明实施例的方案及其效果，以下给出一个具体应用示例。本领域技术人员应理解，该示例仅为了便于理解本发明，其任何具体细节并非意在以任何方式限制本发明。
图2a、图2b、图2c分别示出了根据本发明的一个实施例的中心角分别8度、18度、28度的角度部分叠加数据的示意图。
图3示出了根据图2a、图2b、图2c的中心角分别8度、18度、28度的角度子波的示意图。
图4a示出了根据本发明的一个实施例的纵波速度与横波速度的交会图及其拟合关系的示意图。
图4b示出了根据本发明的一个实施例的纵波速度与密度的交会图及其拟合关系的示意图。
根据测井数据，通过拟合纵波速度与横波速度的线性关系，获得Castagna经验公式的斜率Cαβ＝0.4557，如图4a所示；通过拟合纵波速度与密度的指数关系，获得Gardnar经验公式的指数Cαρ＝0.2320，如图4b所示。进而根据公式(1)建立纵横波速度比反射率与纵波速度反射率的岩石物理关系，根据公式(2)建立纵横波速度比反射率与密度反射率的岩石物理关系。
根据岩石物理关系，并结合Aki-Richards近似公式，通过公式(3)构建基于岩石物理关系约束的AVO近似公式；根据角度叠加数据、角度子波褶积矩阵以及基于岩石物理关系约束的AVO近似公式，通过公式(4)建立AVO正演方程组；确定目的层段的薄层最大时间厚度，构建反射系数奇偶分解矩阵为公式(5)；根据AVO正演方程组与反射系数奇偶分解矩阵，构建反演目标函数为公式(6)；根据反演目标函数，利用迭代重加权最小二乘算法获取弹性参数反射率的奇偶分量系数，然后通过公式(10)计算弹性参数反射率，最后通过道积分获得弹性参数，并对弹性参数进行低频补偿。
图5a、图5b、图5c分别示出了根据本发明的一个实施例的纵横波速度比反演剖面、纵波速度反演剖面、密度反演剖面的示意图。由图可以看出反演结果具有较高的垂向分辨率，特别是在1600ms附近，展示了较为丰富的储层细节信息。
图6a、图6b、图6c分别示出了根据本发明的一个实施例的井旁道的纵横波速度比反演结果、纵波速度反演结果、密度反演结果的示意图，验证井位于CDP301处，黑线为实测井曲线，灰线为井旁道反演结果，从图中可以看出井旁道反演结果与实测曲线反演结果吻合度较高，并且具有较高的垂向分辨率。
综上所述，本发明通过计算纵波速度反射率、纵横波速度比反射率以及密度反射率，结合岩石物理关系约束与反射系数奇偶分解理论，实现了关于纵横波速度比的稀疏层反演。
本领域技术人员应理解，上面对本发明的实施例的描述的目的仅为了示例性地说明本发明的实施例的有益效果，并不意在将本发明的实施例限制于所给出的任何示例。
根据本发明的实施例，提供了一种基于岩石物理关系约束的叠前稀疏层反演系统，其特征在于，该系统包括：存储器，存储有计算机可执行指令；处理器，所述处理器运行所述存储器中的计算机可执行指令，执行以下步骤：步骤1：读取角度叠加数据和角度子波，并根据角度子波建立角度子波褶积矩阵；步骤2：根据测井数据，建立岩石物理关系；步骤3：根据岩石物理关系与Aki-Richards近似式，建立基于岩石物理关系约束的AVO近似公式；步骤4：根据角度叠加数据、角度子波褶积矩阵以及AVO近似公式，建立AVO正演方程组；步骤5：确定目的层段的薄层最大时间厚度，根据反射系数奇偶分解理论，构建反射系数奇偶分解矩阵；步骤6：根据AVO正演方程组与反射系数奇偶分解矩阵，并在最小化弹性参数奇偶分量系数的1范数的约束下，建立反演目标函数；步骤7：根据反演目标函数，计算弹性参数反演结果并进行低频补偿，获得最终的弹性参数反演结果。
在一个示例中，步骤2包括：步骤21：根据测井数据，通过拟合纵波速度曲线与横波速度曲线的线性关系，获得Castagna经验公式的斜率Cαβ；步骤22：根据测井数据，通过拟合纵波速度曲线与密度曲线的指数关系，获得Gardnar经验公式的指数Cαρ；步骤23：根据Castagna经验公式的斜率，并结合弹性参数反射率之间的转换关系，利用公式(1)建立纵波速度反射率与纵横波速度比反射率之间的岩石物理关系：
Figure BDA0002110886280000201
其中，Rα表示纵波速度反射率，Rγ表示纵横波速度比反射率，ΔRα表示纵波速度反射率的扰动量，
Figure BDA0002110886280000202
表示背景纵横波速度比；步骤24：根据所Gardnar经验公式的指数与Castagna经验公式的斜率，并结合弹性参数反射率之间的转换关系，利用公式(2)建立密度反射率与纵横波速度比反射率之间的岩石物理关系：
Figure BDA0002110886280000203
其中，Rρ表示纵波速度反射率，ΔRρ表示纵波速度反射率的扰动量。
在一个示例中，AVO近似公式为：
Rpp(θ)＝A(θ)Rγ+B(θ)ΔRα+C(θ)ΔRρ(3)
其中，
Figure BDA0002110886280000204
Figure BDA0002110886280000205
Rpp(θ)表示入射角为θ的角度反射系数。
在一个示例中，AVO正演方程组为：
d＝Gr (4)
其中，
Figure BDA0002110886280000206
d(θi)表示入射角为θi的角度叠加数据，其中i＝1,2,…K，K表示入射角的个数；
Figure BDA0002110886280000207
Rγ、ΔRα和ΔRρ分别表示纵横波速度比反射率向量、纵波速度反射率扰动量向量和密度反射率扰动量向量；
Figure BDA0002110886280000211
W(θi)表示入射角为θi的角度子波褶积矩阵，其中i＝1,2,…K；A(θi)、B(θi)、C(θi)均为对角阵，其中i＝1,2,…K；其对角阵元素为：
Figure BDA0002110886280000212
其中，j＝1,2,…N，N表示反射系数的采样点个数。
在一个示例中，反射系数奇偶分解矩阵为：
D＝[De Do] (5)
其中，De表示反射系数偶分量矩阵，
Figure BDA0002110886280000213
Do表示反射系数奇分量矩阵，
Figure BDA0002110886280000221
其中，M表示薄层最大时间厚度对应的采样点个数。
在一个示例中，反演目标函数为：
Figure BDA0002110886280000222
其中，x表示待反演参数，
Figure BDA0002110886280000223
xγ、xα和xρ分别表示纵横波速度比反射率的奇偶分量系数、纵波速度反射率扰动量的奇偶分量系数以及密度反射率扰动量的奇偶分量系数；
Figure BDA0002110886280000224
λ表示稀疏约束权重；上标T表示矩阵转置。
在一个示例中，步骤7包括：步骤71：根据反演目标函数，利用迭代重加权最小二乘算法进行求解，计算待反演参数；步骤72：根据待反演参数，计算弹性参数反射率；步骤73：根据弹性参数反射率，计算弹性参数，并进行低频补偿。
在一个示例中，步骤71中的迭代重加权最小二乘算法的计算步骤为：步骤711：设置待反演参数x的初始值为零向量；步骤712：计算非一致性加权矩阵Q，Q为对角矩阵，其对角线元素为：
Figure BDA0002110886280000225
其中，ε表示稳定系数，xl表示x的第l个元素，l＝1,2,...,3N；步骤713：通过结合Cholesky分解和上下三角矩阵分解，求解线性方程组
Figure BDA0002110886280000231
得到待反演参数x；步骤714：重复步骤712-步骤713，直至达到规定的迭代次数，输出待反演参数x。
本系统通过计算纵波速度反射率、纵横波速度比反射率以及密度反射率，结合岩石物理关系约束与反射系数奇偶分解理论，实现了关于纵横波速度比的稀疏层反演。
本领域技术人员应理解，上面对本发明的实施例的描述的目的仅为了示例性地说明本发明的实施例的有益效果，并不意在将本发明的实施例限制于所给出的任何示例。
以上已经描述了本发明的各实施例，上述说明是示例性的，并非穷尽性的，并且也不限于所披露的各实施例。在不偏离所说明的各实施例的范围和精神的情况下，对于本技术领域的普通技术人员来说许多修改和变更都是显而易见的。

Claims (8)

1.一种基于岩石物理关系约束的叠前稀疏层反演方法，其特征在于，包括：
步骤1：读取角度叠加数据和角度子波，并根据所述角度子波建立角度子波褶积矩阵；
步骤2：根据测井数据，建立岩石物理关系；
步骤3：根据所述岩石物理关系与Aki-Richards近似式，建立基于岩石物理关系约束的AVO近似公式；
步骤4：根据所述角度叠加数据、所述角度子波褶积矩阵以及所述AVO近似公式，建立AVO正演方程组；
步骤5：确定目的层段的薄层最大时间厚度，根据反射系数奇偶分解理论，构建反射系数奇偶分解矩阵；
步骤6：根据所述AVO正演方程组与所述反射系数奇偶分解矩阵，并在最小化弹性参数奇偶分量系数的1范数的约束下，建立反演目标函数；
步骤7：根据所述反演目标函数，计算弹性参数反演结果并进行低频补偿，获得最终的弹性参数反演结果；
其中，所述步骤2包括：
步骤21：根据测井数据，通过拟合纵波速度曲线与横波速度曲线的线性关系，获得Castagna经验公式的斜率Cαβ；
步骤22：根据测井数据，通过拟合纵波速度曲线与密度曲线的指数关系，获得Gardnar经验公式的指数Cαρ；
步骤23：根据所述Castagna经验公式的斜率，并结合弹性参数反射率之间的转换关系，利用公式(1)建立纵波速度反射率与纵横波速度比反射率之间的岩石物理关系：
Figure FDA0003507507790000011
其中，Rα表示纵波速度反射率，Rγ表示纵横波速度比反射率，ΔRα表示纵波速度反射率的扰动量，
Figure FDA0003507507790000021
表示背景纵横波速度比；
步骤24：根据所Gardnar经验公式的指数与所述Castagna经验公式的斜率，并结合弹性参数反射率之间的转换关系，利用公式(2)建立密度反射率与纵横波速度比反射率之间的岩石物理关系：
Figure FDA0003507507790000022
其中，Rρ表示密度反射率，ΔRρ表示密度反射率的扰动量。
2.根据权利要求1所述的基于岩石物理关系约束的叠前稀疏层反演方法，其中，所述AVO近似公式为：
Rpp(θ)＝A(θ)Rγ+B(θ)ΔRα+C(θ)ΔRρ (3)
其中，
Figure FDA0003507507790000023
Figure FDA0003507507790000024
Rpp(θ)表示入射角为θ的角度反射系数。
3.根据权利要求1所述的基于岩石物理关系约束的叠前稀疏层反演方法，其中，所述AVO正演方程组为：
d＝Gr (4)
其中，
Figure FDA0003507507790000025
d(θi)表示入射角为θi的角度叠加数据，其中i＝1,2,…K，K表示入射角的个数；
Figure FDA0003507507790000026
Rγ、ΔRα和ΔRρ分别表示纵横波速度比反射率向量、纵波速度反射率扰动量向量和密度反射率扰动量向量；
Figure FDA0003507507790000031
W(θi)表示入射角为θi的角度子波褶积矩阵，其中i＝1,2,…K；A(θi)、B(θi)、C(θi)均为对角阵，其中i＝1,2,…K；其对角阵元素为：
Figure FDA0003507507790000032
其中，j＝1,2,…N，N表示反射系数的采样点个数。
4.根据权利要求3所述的基于岩石物理关系约束的叠前稀疏层反演方法，其中，所述反射系数奇偶分解矩阵为：
D＝[De Do] (5)
其中，De表示反射系数偶分量矩阵，
Figure FDA0003507507790000033
Do表示反射系数奇分量矩阵，
Figure FDA0003507507790000041
其中，M表示薄层最大时间厚度对应的采样点个数。
5.根据权利要求4所述的基于岩石物理关系约束的叠前稀疏层反演方法，其中，所述反演目标函数为：
Figure FDA0003507507790000042
其中，x表示待反演参数，
Figure FDA0003507507790000043
xγ、xα和xρ分别表示纵横波速度比反射率的奇偶分量系数、纵波速度反射率扰动量的奇偶分量系数以及密度反射率扰动量的奇偶分量系数；
Figure FDA0003507507790000044
λ表示稀疏约束权重；上标T表示矩阵转置。
6.根据权利要求5所述的基于岩石物理关系约束的叠前稀疏层反演方法，其中，所述步骤7包括：
步骤71：根据所述反演目标函数，利用迭代重加权最小二乘算法进行求解，计算待反演参数；
步骤72：根据所述待反演参数，计算弹性参数反射率；
步骤73：根据弹性参数反射率，计算弹性参数，并进行低频补偿。
7.根据权利要求6所述的基于岩石物理关系约束的叠前稀疏层反演方法，其中，所述步骤71中的所述迭代重加权最小二乘算法的计算步骤为：
步骤711：设置待反演参数x的初始值为零向量；
步骤712：计算非一致性加权矩阵Q，Q为对角矩阵，其对角线元素为：
Figure FDA0003507507790000051
其中，ε表示稳定系数，xl表示x的第l个元素，l＝1,2,...,3N；
步骤713：通过结合Cholesky分解和上下三角矩阵分解，求解线性方程组
Figure FDA0003507507790000052
得到待反演参数x；
步骤714：重复步骤712-步骤713，直至达到规定的迭代次数，输出待反演参数x。
8.一种基于岩石物理关系约束的叠前稀疏层反演系统，其特征在于，该系统包括：
存储器，存储有计算机可执行指令；
处理器，所述处理器运行所述存储器中的计算机可执行指令，执行以下步骤：
步骤1：读取角度叠加数据和角度子波，并根据所述角度子波建立角度子波褶积矩阵；
步骤2：根据测井数据，建立岩石物理关系；
步骤3：根据所述岩石物理关系与Aki-Richards近似式，建立基于岩石物理关系约束的AVO近似公式；
步骤4：根据所述角度叠加数据、所述角度子波褶积矩阵以及所述AVO近似公式，建立AVO正演方程组；
步骤5：确定目的层段的薄层最大时间厚度，根据反射系数奇偶分解理论，构建反射系数奇偶分解矩阵；
步骤6：根据所述AVO正演方程组与所述反射系数奇偶分解矩阵，并在最小化弹性参数奇偶分量系数的1范数的约束下，建立反演目标函数；
步骤7：根据所述反演目标函数，计算弹性参数反演结果并进行低频补偿，获得最终的弹性参数反演结果；
其中，所述步骤2包括：
步骤21：根据测井数据，通过拟合纵波速度曲线与横波速度曲线的线性关系，获得Castagna经验公式的斜率Cαβ；
步骤22：根据测井数据，通过拟合纵波速度曲线与密度曲线的指数关系，获得Gardnar经验公式的指数Cαρ；
步骤23：根据所述Castagna经验公式的斜率，并结合弹性参数反射率之间的转换关系，利用公式(1)建立纵波速度反射率与纵横波速度比反射率之间的岩石物理关系：
Figure FDA0003507507790000061
其中，Rα表示纵波速度反射率，Rγ表示纵横波速度比反射率，ΔRα表示纵波速度反射率的扰动量，
Figure FDA0003507507790000062
表示背景纵横波速度比；
步骤24：根据所Gardnar经验公式的指数与所述Castagna经验公式的斜率，并结合弹性参数反射率之间的转换关系，利用公式(2)建立密度反射率与纵横波速度比反射率之间的岩石物理关系：
Figure FDA0003507507790000063
其中，Rρ表示密度反射率，ΔRρ表示密度反射率的扰动量。
CN201910571036.8A 2019-06-28 2019-06-28 基于岩石物理关系约束的叠前稀疏层反演方法及系统 Active CN112147683B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910571036.8A CN112147683B (zh) 2019-06-28 2019-06-28 基于岩石物理关系约束的叠前稀疏层反演方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910571036.8A CN112147683B (zh) 2019-06-28 2019-06-28 基于岩石物理关系约束的叠前稀疏层反演方法及系统

Publications (2)

Publication Number Publication Date
CN112147683A CN112147683A (zh) 2020-12-29
CN112147683B true CN112147683B (zh) 2022-06-21

Family

ID=73869254

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910571036.8A Active CN112147683B (zh) 2019-06-28 2019-06-28 基于岩石物理关系约束的叠前稀疏层反演方法及系统

Country Status (1)

Country Link
CN (1) CN112147683B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113985480B (zh) * 2021-11-08 2024-01-26 西南石油大学 一种基于角度校正的avo反演方法及装置

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102169190A (zh) * 2011-01-06 2011-08-31 中国科学院地质与地球物理研究所 一种调制补充子空间的井约束叠前弹性参数反演方法
CN103245970A (zh) * 2012-02-08 2013-08-14 中国石油化工股份有限公司 叠前地震宽角度反演方法
CN104422956A (zh) * 2013-08-22 2015-03-18 中国石油化工股份有限公司 一种基于稀疏脉冲反演的高精度地震谱分解方法
CN104898167A (zh) * 2015-06-18 2015-09-09 中国海洋石油总公司 一种稀疏正则化的岩石物理弹性参数提取方法
WO2018027003A1 (en) * 2016-08-05 2018-02-08 Chevron U.S.A. Inc. System and method for petro-elastic modeling

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106842310B (zh) * 2015-12-04 2020-10-13 中国石油化工股份有限公司 叠前地震四参数同步反演方法
US10928536B2 (en) * 2017-12-07 2021-02-23 Saudi Arabian Oil Company Mapping chemostratigraphic signatures of a reservoir with rock physics and seismic inversion

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102169190A (zh) * 2011-01-06 2011-08-31 中国科学院地质与地球物理研究所 一种调制补充子空间的井约束叠前弹性参数反演方法
CN103245970A (zh) * 2012-02-08 2013-08-14 中国石油化工股份有限公司 叠前地震宽角度反演方法
CN104422956A (zh) * 2013-08-22 2015-03-18 中国石油化工股份有限公司 一种基于稀疏脉冲反演的高精度地震谱分解方法
CN104898167A (zh) * 2015-06-18 2015-09-09 中国海洋石油总公司 一种稀疏正则化的岩石物理弹性参数提取方法
WO2018027003A1 (en) * 2016-08-05 2018-02-08 Chevron U.S.A. Inc. System and method for petro-elastic modeling

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
基于低频软约束-基追踪薄层反演方法研究;李治昊等;《地球物理学进展》;20181231;第33卷(第6期);第2403-2408页 *
基于低频软约束的叠前AVA稀疏层反演;张丰麒等;《石油地球物理勘探》;20170831;第52卷(第4期);第770-782页 *

Also Published As

Publication number Publication date
CN112147683A (zh) 2020-12-29

Similar Documents

Publication Publication Date Title
Das et al. Convolutional neural network for seismic impedance inversion
US10310113B2 (en) Q-compensated full wavefield inversion
WO2018129844A1 (zh) 一种地震绕射波分离方法和装置
Korenaga et al. Seismic tomography of Shatsky Rise by adaptive importance sampling
US9291736B2 (en) Surface-consistent amplitude and deconvolution simultaneous joined inversion
US9348049B2 (en) Simultaneous joint estimation of the P-P and P-S residual statics
Yablokov et al. An artificial neural network approach for the inversion of surface wave dispersion curves
US11733413B2 (en) Method and system for super resolution least-squares reverse time migration
CN106556861B (zh) 一种基于全方位地震资料的方位avo反演方法
US20070073488A1 (en) Handling of static corrections in multiple prediction
CN103293552A (zh) 一种叠前地震资料的反演方法及系统
CN110895348B (zh) 一种地震弹性阻抗低频信息提取方法、系统及存储介质
CN109212599A (zh) 一种地震数据的叠前同步反演方法
Ball et al. A joint M onte C arlo analysis of seafloor compliance, R ayleigh wave dispersion and receiver functions at ocean bottom seismic stations offshore N ew Z ealand
Thiel et al. Comparison of acoustic and elastic full‐waveform inversion of 2D towed‐streamer data in the presence of salt
CN116047583A (zh) 基于深度卷积神经网络的自适应波阻抗反演方法及系统
CN112147683B (zh) 基于岩石物理关系约束的叠前稀疏层反演方法及系统
CN113970787B (zh) 物性参数反演方法、装置、计算机设备和存储介质
CN108957545A (zh) 气枪阵列子波方向性反褶积方法及系统
CN115327625B (zh) 一种储层岩性识别方法
CN112147681A (zh) 基于γ_Zoeppritz方程的叠前反演方法及系统
CN114325832A (zh) 一种裂缝参数和弹性参数同步反演方法及系统
Zand et al. Integrated algorithm for high‐resolution crustal‐scale imaging using complementary OBS and streamer data
Pendrel et al. Direct modeling of reservoir properties from seismic partial-angle stacks
Pérez et al. Estimating sparse-spike attributes from ava data using a fast iterative shrinkage-thresholding algorithm and least squares

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