CN114609669A - 基于方位弹性阻抗的hti型裂缝储层参数预测方法及系统 - Google Patents

基于方位弹性阻抗的hti型裂缝储层参数预测方法及系统 Download PDF

Info

Publication number
CN114609669A
CN114609669A CN202210500271.8A CN202210500271A CN114609669A CN 114609669 A CN114609669 A CN 114609669A CN 202210500271 A CN202210500271 A CN 202210500271A CN 114609669 A CN114609669 A CN 114609669A
Authority
CN
China
Prior art keywords
hti
elastic impedance
fracture
reservoir
parameters
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
Application number
CN202210500271.8A
Other languages
English (en)
Other versions
CN114609669B (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.)
Qingdao Youjin Energy Technology Co ltd
Original Assignee
Qingdao Youjin Energy Technology 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 Qingdao Youjin Energy Technology Co ltd filed Critical Qingdao Youjin Energy Technology Co ltd
Priority to CN202210500271.8A priority Critical patent/CN114609669B/zh
Publication of CN114609669A publication Critical patent/CN114609669A/zh
Application granted granted Critical
Publication of CN114609669B publication Critical patent/CN114609669B/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

本发明属于地球物理勘测技术领域，提供了一种基于方位弹性阻抗的HTI型裂缝储层参数预测方法及系统，包括以下步骤：获取HTI型裂缝储层的第一刚性矩阵和第二刚性矩阵，计算准裂缝密度参数和PP波反射系数方程；对所得到的PP波反射系数方程进行反演计算，构建方位弹性阻抗方程；选择出以反演目标函数最小的方位弹性阻抗，使用被选择出的方位弹性阻抗预测HTI型裂缝储层参数。本发明通过参数化HTI型裂缝储层，计算准裂缝密度参数、垂直裂缝储层PP波反射系数方程和方位弹性阻抗方程，构建反演目标函数，预测脆性指数，提高了HTI型裂缝储层参数的预测精度。

Description

基于方位弹性阻抗的HTI型裂缝储层参数预测方法及系统
技术领域
本发明属于地球物理勘测技术领域，具体涉及一种基于方位弹性阻抗的HTI型裂缝储层参数预测方法及系统。
背景技术
本部分的陈述仅仅是提供了与本发明相关的背景技术信息，不必然构成在先技术。
我国页岩气储层地质的地球物理条件差异较大，勘探开发难度较大。工程甜点评价是页岩等非常规油气储层勘探开发的重要环节，其中脆性和裂缝密度预测是工程甜点评价中非常关键的参数。AVAZ（Amplitude Variation with Azimuth）反演技术是研究裂缝储层的重要手段。
裂缝是非常规（致密油气、页岩油气等）储层重要的油气存储空间与运移通道，脆性可有效指导储层开发方案的设计。稳定而有效的刻画地下储层脆性和裂缝特征，可提高裂缝油气藏的勘探与开发效率。各向同性均匀岩石中发育一组垂直定向硬币型裂缝时，可将其等效为HTI介质，地震AVAZ反演方法是垂直裂缝储层的有效预测方法。但针对垂直裂缝储层脆性和裂缝参数预测的模型参数化和反问题的稳定求解依旧存在较大的挑战。
脆性是岩石的一种固有性质，表现为岩石受到外力时产生较小的形变甚至不发生形变就可以发生破碎。目前，国内外学者研究储层脆性的方法主要包括：脆性矿物含量评价法，应力-应变曲线评价法和弹性参数评价法。针对脆性地震评价一般采用弹性参数评价法，该方法中，弹性力学参数主要包括岩石的杨氏模量和泊松比；杨氏模量与岩石的刚度有关，泊松比能够表征岩石的塑性。脆性地震预测方法通过综合利用地震数据与测井资料，通过地震反演方法提取出杨氏模量和泊松比这两个弹性参数，从而实现脆性甜点评价。国内外许多专家分析了杨氏模量和泊松比与岩石脆性之间的关系，利用这两个参数提出了多种脆性指数表达方法；如何获得能够直接反映杨氏模量等力学参数的各向异性反射系数方程是裂缝储层脆性评价的重点所在。
各向异性反射方程的推导是构建五维地震正演算子和裂缝参数等工程甜点地震评价的基础；但是现有技术中反射系数方程所包含的待反演参数较多，严重影响地震预测的精确度；因此，降低表征地震反射系数方程的参数个数，是提高裂缝储层反演结果合理性的重要手段。
反问题求解算法的选取是储层脆性及裂缝特征精确刻画的关键。最小二乘方法是最基础的反问题求解手段，它对正演算子合理性具有较高的要求，同时，对地震资料中的噪音较为敏感，通过加入约束项可在一定程度上提高待反演参数的预测准确度；但是，难以量化预测结果的不确定性。
发明内容
为了解决上述问题，本发明提出了一种基于方位弹性阻抗的HTI型裂缝储层参数预测方法及系统，通过参数化HTI型裂缝储层，计算准裂缝密度参数、垂直裂缝储层PP波反射系数方程和方位弹性阻抗方程，构建反演目标函数，预测脆性指数，提高了HTI型裂缝储层参数的预测精度。
根据一些实施例，本发明的第一方案提供了一种基于方位弹性阻抗的HTI型裂缝储层参数预测方法，采用如下技术方案：
一种基于方位弹性阻抗的HTI型裂缝储层参数预测方法，包括以下步骤：
获取HTI型裂缝储层的第一刚性矩阵
Figure 194876DEST_PATH_IMAGE001
和第二刚性矩阵
Figure 917367DEST_PATH_IMAGE002
，计算准裂缝密度参数和PP波反射系数方程；
对所得到的PP波反射系数方程进行反演计算，构建方位弹性阻抗方程；
选择出以反演目标函数最小的方位弹性阻抗，使用被选择出的方位弹性阻抗预测HTI型裂缝储层参数。
作为进一步的技术限定，所预测的HTI型裂缝储层参数包括准裂缝密度参数和脆性指数。
进一步的，所述脆性指数至少包括全工区参数的最小值和最大值。
作为进一步的技术限定，所述第一刚性矩阵为背景柔度矩阵与裂缝柔度矩阵之和；所述第一刚性矩阵与所述第二刚性矩阵互为逆矩阵。
进一步的，所述背景柔度矩阵与所述HTI型裂缝储层的杨氏模量、泊松比和剪切模量相关；所述裂缝柔度矩阵与所述HTI型裂缝储层的扰动裂缝的法向弱度和切向弱度相关。
进一步的，对所述第二刚性矩阵进行线性处理，结合散射理论，得到PP波反射系数。
进一步的，所述PP波反射系数与所述HTI型裂缝储层的质量密度、PP波在所述HTI型裂缝储层上的入射角，以及观测方位角相关。
进一步的，基于所得到的PP波反射系数构建所述PP波反射系数方程，通过约束稀疏脉冲反演法对所述PP波反射系数方程进行反演计算，得到每个方位的弹性阻抗，结合道积分思想，构建方位弹性阻抗方程。
进一步的，所述反演目标函数采用约束稀疏脉冲反演法的目标函数，即为：
Figure 539978DEST_PATH_IMAGE003
，
其中，下标m代表第m个采样点，
Figure 672145DEST_PATH_IMAGE004
是入射角为θ，方位角为ϕ的弹性阻抗反射系数；
Figure 12996DEST_PATH_IMAGE005
和
Figure 521600DEST_PATH_IMAGE006
分别代表实际和合成地震记录；
Figure 998718DEST_PATH_IMAGE007
是通过测井数据估计的弹性阻抗梯度模型约束；
Figure 567365DEST_PATH_IMAGE008
和
Figure 129933DEST_PATH_IMAGE009
分别代表p阶和q阶范数；常数
Figure 170789DEST_PATH_IMAGE010
和
Figure 767993DEST_PATH_IMAGE011
是通过工区先验认识给定的权重系数。
根据一些实施例，本发明的第二方案提供了一种基于方位弹性阻抗的HTI型裂缝储层参数预测系统，采用如下技术方案：
一种基于方位弹性阻抗的HTI型裂缝储层参数预测系统，包括：
获取模块，被配置为获取HTI型裂缝储层的第一刚性矩阵和第二刚性矩阵，计算准裂缝密度参数和PP波反射系数方程；
构建方程模块，被配置为对所得到的PP波反射系数方程进行反演计算，构建方位弹性阻抗方程；
参数预测模块，被配置为选择出以反演目标函数最小的方位弹性阻抗，使用被选择出的方位弹性阻抗预测HTI型裂缝储层参数。
与现有技术相比，本发明的有益效果为：
本发明通过参数化HTI型裂缝储层，计算准裂缝密度参数、垂直裂缝储层PP波反射系数方程和方位弹性阻抗方程，构建反演目标函数，预测脆性指数，提高了HTI型裂缝储层参数的预测精度。
附图说明
构成本发明的一部分的说明书附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。
图1是本发明实施例一种的基于方位弹性阻抗的HTI型裂缝储层参数预测方法的流程图；
图2是本发明实施例一中的各向同性背景中发育一组垂直定向分布的裂缝所形成的等效HTI介质模型示意图；
图3是本发明实施例二中的基于方位弹性阻抗的HTI型裂缝储层参数预测系统的结构框图。
具体实施方式
下面结合附图与实施例对本发明作进一步说明。
应该指出，以下详细说明都是例示性的，旨在对本发明提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本发明所属技术领域的普通技术人员通常理解的相同含义。
需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本发明的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。
在不冲突的情况下，本发明中的实施例及实施例中的特征可以相互组合。
实施例一
本发明实施例一介绍了一种基于方位弹性阻抗的HTI型裂缝储层参数预测方法。
本实施例采用Schoenberg线性滑动理论和Hudson理论参数化HTI型裂缝储层，基于Born逆散射理论和稳相法推导杨氏模量、泊松比、质量密度和裂缝密度四参数直接表征的垂直裂缝储层PP波反射系数方程和方位弹性阻抗方程；再结合贝叶斯反演方法和低频模型正则化约束，构建反演目标函数，形成脆性指数和准裂缝密度参数预测技术；最后结合实际案例应用证明了本实施例所提出方法的有效性。
如图1所示的一种基于方位弹性阻抗的HTI型裂缝储层参数预测方法，包括以下步骤：
获取HTI型裂缝储层的第一刚性矩阵和第二刚性矩阵，计算准裂缝密度参数和PP波反射系数方程；
对所得到的PP波反射系数方程进行反演计算，构建方位弹性阻抗方程；
选择出以反演目标函数最小的方位弹性阻抗，使用被选择出的方位弹性阻抗预测HTI型裂缝储层参数。
当各向同性背景岩石中发育一组垂直定向排列的硬币型裂缝时，可将其看作具有水平对称轴的横向各向同性介质（horizontal transversely isotropic medium，即HTI介质），如图2所示。Schoenberg线性滑动理论指出，定向裂缝介质的整体等效柔度矩阵
Figure 241962DEST_PATH_IMAGE012
可表示为背景柔度矩阵
Figure 291826DEST_PATH_IMAGE013
与裂缝柔度矩阵
Figure 142233DEST_PATH_IMAGE014
之和，同时，介质柔度矩阵是其刚度矩阵的逆矩阵，则垂直裂缝诱导的HTI介质有效刚度矩阵
Figure 593943DEST_PATH_IMAGE015
可表示为：
Figure 504392DEST_PATH_IMAGE016
（1）
对于大多数储层，裂缝密度均较小，对公式（1）进行泰勒级数的展开，仅仅保留一阶项，得
Figure 41553DEST_PATH_IMAGE017
（2）
其中，各向同性背景的柔度矩阵可表示为：
Figure 424212DEST_PATH_IMAGE018
（3）
当垂直定向裂缝组法向与x轴一致时，对于小裂缝密度情况，裂缝柔度矩阵可近似表达为
Figure 730428DEST_PATH_IMAGE019
（4）
其中，
Figure 811779DEST_PATH_IMAGE020
，
Figure 836236DEST_PATH_IMAGE021
和
Figure 28445DEST_PATH_IMAGE022
分别表示背景介质的杨氏模量，泊松比和剪切模量；
Figure 189168DEST_PATH_IMAGE023
是背景介质刚度矩阵第ii个元素，
Figure 441420DEST_PATH_IMAGE024
和
Figure 635729DEST_PATH_IMAGE025
分别表示扰动裂缝的法向弱度和切向弱度。
裂缝法向和切向弱度参数是裂缝密度的线性函数，对于含气储层，即有：
Figure 864585DEST_PATH_IMAGE026
，
Figure 646858DEST_PATH_IMAGE027
（5）
将公式（3）、公式（4）和公式（5）带入公式（2），可得裂缝诱导的HTI介质等效刚度矩阵为：
Figure 381596DEST_PATH_IMAGE028
（6）
其中，
Figure 115066DEST_PATH_IMAGE029
Figure 914657DEST_PATH_IMAGE030
，
Figure 49972DEST_PATH_IMAGE031
，
Figure 644026DEST_PATH_IMAGE032
。
假设某一反射界面将具有相同裂缝倾向的两个HTI介质划分开，且裂缝密度较小，界面上下弹性差异较小，则包含
Figure 864792DEST_PATH_IMAGE033
，
Figure 766994DEST_PATH_IMAGE034
，
Figure 22395DEST_PATH_IMAGE035
和
Figure 787351DEST_PATH_IMAGE036
的项可以忽略，则该界面上下介质刚度参数扰动量可近似表示为：
Figure 495412DEST_PATH_IMAGE037
（7）
其中，
Figure 371227DEST_PATH_IMAGE038
,
Figure 481134DEST_PATH_IMAGE039
,
Figure 416991DEST_PATH_IMAGE040
Figure 612349DEST_PATH_IMAGE041
,
Figure 285995DEST_PATH_IMAGE042
，其中，
Figure 250409DEST_PATH_IMAGE043
，
Figure 357168DEST_PATH_IMAGE044
，
Figure 774242DEST_PATH_IMAGE045
。
根据散射理论，可以得到P-P波反射系数
Figure 257439DEST_PATH_IMAGE046
和散射函数的线性近似关系为：
Figure 76359DEST_PATH_IMAGE047
（8）
其中，
Figure 354019DEST_PATH_IMAGE048
表示PP波散射函数在端点
Figure 523969DEST_PATH_IMAGE049
处的函数值，即有
Figure 539418DEST_PATH_IMAGE050
，（i, j, k, l=1,2,3）（9）
其中，重复的角标指的是Einstein求和；
Figure 275161DEST_PATH_IMAGE051
代表储层质量密度，
Figure 989302DEST_PATH_IMAGE052
是克罗内克函数。
Figure 380969DEST_PATH_IMAGE053
，
Figure 940388DEST_PATH_IMAGE054
，
Figure 733901DEST_PATH_IMAGE055
和
Figure 353363DEST_PATH_IMAGE056
分别代表入射波和反射波的慢度和极化向量，即
Figure 232326DEST_PATH_IMAGE057
（10）
其中，
Figure 855157DEST_PATH_IMAGE058
和
Figure 503176DEST_PATH_IMAGE059
分别代表入射角和观测方位角。
假设参考介质是均匀各向同性的，则公式（9）可表示为
Figure 293539DEST_PATH_IMAGE060
（11）
其中，
Figure 659799DEST_PATH_IMAGE061
，
Figure 826600DEST_PATH_IMAGE062
（12）
角标(m, n)和角标(i, j, k, l)的关系为
Figure 329126DEST_PATH_IMAGE063
,
Figure 555970DEST_PATH_IMAGE064
。
将公式（10）带入公式（12）可得
Figure 143946DEST_PATH_IMAGE065
。将差分矩阵公式（7）、公式（12）、公式（A-1）及公式（11）带入公式（8）中，则得到HTI型裂缝储层P波入射，P波反射的反射系数方程为：
Figure 132016DEST_PATH_IMAGE066
（13）
其中：
Figure 754628DEST_PATH_IMAGE067
，
Figure 886794DEST_PATH_IMAGE068
，
Figure 227645DEST_PATH_IMAGE069
，
Figure 736249DEST_PATH_IMAGE070
为了消除背景各向同性反射系数与裂隙诱导方位各向异性反射系数间的表达差异，根据Mesdag和Quevedo给出的准弱度参数概念，提出一种新的裂隙密度参数化方法，即
Figure 213367DEST_PATH_IMAGE071
，其中，
Figure 516435DEST_PATH_IMAGE072
代表裂缝密度在全工区的平均值；可将qe称为准裂缝密度参数，则公式（13）可变换为：
Figure 79003DEST_PATH_IMAGE073
（14）
弹性阻抗是一种层属性，而且弹性阻抗反演方法比AVO同步反演方法的稳定性更强；因此，本实施例研究了基于HTI介质弹性阻抗的脆性参数和裂缝参数预测方法。根据Connolly的弹性阻抗概念，可扩展得到方位弹性阻抗的定义，即
Figure 385439DEST_PATH_IMAGE074
（15）
其中，EI和
Figure 717063DEST_PATH_IMAGE075
分别代表方位弹性阻抗及其上下层的平均值。
假设上下层介质性质差异较小，则
Figure 456611DEST_PATH_IMAGE076
，然后将公式（14）带入公式（15）可得：
Figure 506476DEST_PATH_IMAGE077
（16）
对公式（16）等号两边同时求积分，并进行标准化处理，可得：
Figure 356882DEST_PATH_IMAGE078
（17）
其中，
Figure 808592DEST_PATH_IMAGE079
代表相应属性在全工区的平均值。
在本实施例中，采用约束稀疏脉冲反演方法反演每个方位的弹性阻抗。
约束稀疏脉冲反演是综合地震资料、测井数据和解释层位等基础数据，以地震解释层位作为地层控制格架，以测井数据为约束条件，从井点出发，由井旁道向外推内插完成地震反演。此方法保持了地震资料原有的横向分辨率，同时纵向上与测井数据吻合度较高。其主要过程可以分为以下三步：
（1）通过最大似然反褶积求得一个具有稀疏特性的反射系数序列；
（2）通过最大似然反演得到相对弹性阻抗；
（3）将测井数据中的低频趋势补充到相对弹性阻抗中，得到地震资料的绝对弹性阻抗，完成整个反演过程；
其中，约束稀疏脉冲反演的目标函数可写为：
Figure 719041DEST_PATH_IMAGE080
（18）
其中，下标m代表第m个采样点，
Figure 990623DEST_PATH_IMAGE081
是入射角为θ，方位角为ϕ的弹性阻抗反射系数；
Figure 638861DEST_PATH_IMAGE082
和
Figure 945077DEST_PATH_IMAGE083
分别代表实际和合成地震记录；
Figure 26428DEST_PATH_IMAGE084
是通过测井数据估计的弹性阻抗梯度模型约束；
Figure 50885DEST_PATH_IMAGE085
和
Figure 243094DEST_PATH_IMAGE086
分别代表p阶和q阶范数；常数
Figure 403817DEST_PATH_IMAGE087
和
Figure 656069DEST_PATH_IMAGE088
是通过工区先验认识给定的权重系数，用于调节公式（18）中所包含的三项对反演结果的贡献度。
非线性反演算法，通过寻找公式（18）的全局最小值即可获得最优弹性阻抗反射系数估计，再利用道积分思想便可求得最终的方位弹性阻抗
Figure 902242DEST_PATH_IMAGE089
，即
Figure 829966DEST_PATH_IMAGE090
，其中，t0代表初始时刻，ti代表第i时刻。
当目标储层拥有N个反射界面，采用P个方位角，每个方位角M个入射角的地震记录，公式（17）可表达为矩阵形式：
Figure 110775DEST_PATH_IMAGE091
（19）
其中,
Figure 533928DEST_PATH_IMAGE092
Figure 267398DEST_PATH_IMAGE093
,其中，
Figure 801410DEST_PATH_IMAGE094
，
Figure 936725DEST_PATH_IMAGE095
，
Figure 530779DEST_PATH_IMAGE096
，
Figure 751545DEST_PATH_IMAGE097
，
Figure 293777DEST_PATH_IMAGE098
，
Figure 50642DEST_PATH_IMAGE099
，
Figure 314134DEST_PATH_IMAGE100
，
Figure 258081DEST_PATH_IMAGE101
，
Figure 898010DEST_PATH_IMAGE102
，
Figure 509382DEST_PATH_IMAGE103
，
Figure 943774DEST_PATH_IMAGE104
，
Figure 369159DEST_PATH_IMAGE105
，
Figure 547199DEST_PATH_IMAGE106
，
Figure 13078DEST_PATH_IMAGE107
，
Figure 618371DEST_PATH_IMAGE108
，
Figure 802490DEST_PATH_IMAGE109
Figure 518642DEST_PATH_IMAGE110
。
由于模型参数间的相关性会导致反演结果的不稳定，利用协方差矩阵奇异值分解思想消除模型参数的相关性。
假设模型参数的协方差矩阵可表示为：
Figure 42290DEST_PATH_IMAGE111
；利用奇异值分解方法，协方差矩阵
Figure 818485DEST_PATH_IMAGE112
可写为
Figure 218461DEST_PATH_IMAGE113
（20）
其中，
Figure 738304DEST_PATH_IMAGE114
是由特征值组成的对角矩阵，
Figure 178775DEST_PATH_IMAGE115
是由特征向量组成的单位正交阵。
因此，N个界面的去相关矩阵
Figure 125872DEST_PATH_IMAGE116
可表示为：
Figure 19003DEST_PATH_IMAGE117
（21）
其中，
Figure 342537DEST_PATH_IMAGE118
代表Kronecker积，
Figure 637515DEST_PATH_IMAGE119
是L阶单位矩阵。所以，去相关后的正演算子可写为：
Figure 755512DEST_PATH_IMAGE120
（22）
柯西分布表征变量具有稀疏分布特征，基于Bayesian反演理论，将柯西分布作为先验概率密度函数
Figure 130081DEST_PATH_IMAGE121
来提高反演结果的分辨率，高斯分布作为似然函数
Figure 257306DEST_PATH_IMAGE122
，则后验概率密度函数
Figure 406790DEST_PATH_IMAGE123
可写为：
Figure 695689DEST_PATH_IMAGE124
（ 23）
其中，
Figure 563413DEST_PATH_IMAGE125
代表噪音方差，
Figure 228749DEST_PATH_IMAGE126
代表模型参数的方差；
Figure 232740DEST_PATH_IMAGE127
是指数函数。
针对公式（23）的最大后验概率估计，可将其等效写为：
Figure 958119DEST_PATH_IMAGE128
（ 24）
通过融合从测井、地质等方法获得的初始模型低频信息，最终的目标函数可写为：
Figure 65138DEST_PATH_IMAGE129
（25）
式中，
Figure 534166DEST_PATH_IMAGE130
，
Figure 392662DEST_PATH_IMAGE131
，
Figure 288943DEST_PATH_IMAGE132
表示模型参数的初值；
Figure 865680DEST_PATH_IMAGE133
代表模型参数的约束系数，控制低频约束模型对反演结果的贡献度；对公式（25）求解，可得
Figure 872819DEST_PATH_IMAGE134
（26）
式中，
Figure 851402DEST_PATH_IMAGE135
，
Figure 918584DEST_PATH_IMAGE136
。
其中，
Figure 976758DEST_PATH_IMAGE137
称为柯西稀疏矩阵。
采用迭代重加权最小二乘算法即可求解公式（26），获得杨氏模量、泊松比、质量密度和准裂缝密度的估计；准裂缝密度可直接指示裂缝发育情况，将估计的杨氏模量和泊松比带入Rickman脆性指数表达式，即可实现岩石脆性的描述，Rickman脆性指数BI可表示为：
Figure 787588DEST_PATH_IMAGE138
（27）
其中，下标min和max分别代表相应参数在全工区的最小值和最大值，可通过测井数据等先验信息获得。除此之外，本实施例中所提出的方法应用到实际工区时，需要根据工区先验信息，给出整个目标区的裂缝方位分布，并将方位地震数据中蕴含的裂缝方位全部校正到0度方向。
以中国某裂缝型页岩气工区为例进行本实施例中所介绍方法的合理性验证。
该工区的宽方位地震数据包括六个方位扇，即0度-30度、30度-60度、60度-90度、90度-120度、120度-150度及150度-180度，叠加后将每个角度范围的中间值作为该叠加数据的方位角。每个方位角对应的地震数据包括3个入射角，即10度，19度和29度。根据工区前期的勘探认识，采用方位角为15度和105度的地震数据进行方位弹性阻抗反演，最后采用本实施例提出的反演方法实现准裂缝密度及脆性参数的预测。
反演的方位弹性阻抗在入射角和方位角域均存在差异，这为基于方位弹性阻抗的脆性参数和准裂缝密度参数估计奠定基础。在有效融合杨氏模量、泊松比、密度、准裂缝密度参数的低频先验信息后，进而可估计储层脆性参数。
目的层的准裂缝密度参数较大，表明裂缝较为发育，但脆性预测值较低。横向对比可发现，目的层整体表现出较好的工程甜点属性，即高脆性、高裂缝密度特征，可以沿层进行水平钻井作业。后期水平井钻进显示，目的层裂缝较为发育，脆性指数较高，是良好的油气储层；预测结果与实钻结果较为吻合，证明了本实施例所介绍方法的合理性。
本实施例推导了岩石力学参数和准裂缝密度参数表征的HTI介质四参数反射系数方程，提出了一个基于贝叶斯反演框架的脆性参数和准裂缝密度参数反演方法；基于Schoenberg线性滑动理论和Born逆散射理论，在界面弱对比和裂缝密度较小的假设下，推导了杨氏模量、泊松比、质量密度和准裂缝密度表征的垂直裂缝储层PP波反射系数方程，并进一步给出了方位弹性阻抗方程，以提高脆性和准裂缝密度参数估计精度；然后结合Bayesian反演框架，提出了裂缝型含气储层脆性和裂缝密度直接反演方法，工区案例证明所提方法是可行的。
本实施例所介绍的方法主要针对含气裂缝储层，但如果采用饱含水时的裂缝弱度参数表达式（即公式5）来重新参数化刚度矩阵，即公式（6），便可将本实施例中所介绍的方法扩展到含油水裂缝储层的工程甜点评价中。
实施例二
本发明实施例二介绍了一种基于方位弹性阻抗的HTI型裂缝储层参数预测系统。
如图3所示的一种基于方位弹性阻抗的HTI型裂缝储层参数预测系统。
获取模块，被配置为获取HTI型裂缝储层的第一刚性矩阵和第二刚性矩阵，计算准裂缝密度参数和PP波反射系数方程；
构建方程模块，被配置为对所得到的PP波反射系数方程进行反演计算，构建方位弹性阻抗方程；
参数预测模块，被配置为选择出以反演目标函数最小的方位弹性阻抗，使用被选择出的方位弹性阻抗预测HTI型裂缝储层参数。
详细步骤与实施例一提供的基于方位弹性阻抗的HTI型裂缝储层参数预测方法相同，在此不再赘述。
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种基于方位弹性阻抗的HTI型裂缝储层参数预测方法，其特征在于，包括以下步骤：
获取HTI型裂缝储层的第一刚性矩阵和第二刚性矩阵，计算准裂缝密度参数和PP波反射系数方程；
对所得到的PP波反射系数方程进行反演计算，构建方位弹性阻抗方程；
选择出以反演目标函数最小的方位弹性阻抗，使用被选择出的方位弹性阻抗预测HTI型裂缝储层参数。
2.如权利要求1中所述的一种基于方位弹性阻抗的HTI型裂缝储层参数预测方法，其特征在于，所述第一刚性矩阵为背景柔度矩阵与裂缝柔度矩阵之和；所述第一刚性矩阵与所述第二刚性矩阵互为逆矩阵。
3.如权利要求2中所述的一种基于方位弹性阻抗的HTI型裂缝储层参数预测方法，其特征在于，所述背景柔度矩阵与所述HTI型裂缝储层的杨氏模量、泊松比和剪切模量相关；所述裂缝柔度矩阵与所述HTI型裂缝储层的扰动裂缝的法向弱度和切向弱度相关。
4.如权利要求3中所述的一种基于方位弹性阻抗的HTI型裂缝储层参数预测方法，其特征在于，对所述第二刚性矩阵进行线性处理，结合Born近似和稳定相位法，得到PP波反射系数。
5.如权利要求4中所述的一种基于方位弹性阻抗的HTI型裂缝储层参数预测方法，其特征在于，所述PP波反射系数与所述HTI型裂缝储层的质量密度、PP波在所述HTI型裂缝储层上的入射角，以及观测方位角相关。
6.如权利要求5中所述的一种基于方位弹性阻抗的HTI型裂缝储层参数预测方法，其特征在于，基于所得到的PP波反射系数构建所述PP波反射系数方程，通过约束稀疏脉冲反演法对所述PP波反射系数方程进行反演计算，得到每个方位的弹性阻抗，结合道积分思想，构建方位弹性阻抗方程。
7.如权利要求6中所述的一种基于方位弹性阻抗的HTI型裂缝储层参数预测方法，其特征在于，所述反演目标函数采用约束稀疏脉冲反演法的目标函数，即为：
Figure 685299DEST_PATH_IMAGE001
，
其中，下标m代表第m个采样点，
Figure 564262DEST_PATH_IMAGE002
是入射角为θ，方位角为ϕ的弹性阻抗反射系数；
Figure 927373DEST_PATH_IMAGE003
和
Figure 575392DEST_PATH_IMAGE004
分别代表实际和合成地震记录；
Figure 631335DEST_PATH_IMAGE005
是通过测井数据估计的弹性阻抗梯度模型约束；
Figure 732015DEST_PATH_IMAGE006
和
Figure 898816DEST_PATH_IMAGE007
分别代表p阶和q阶范数；常数
Figure 666921DEST_PATH_IMAGE008
和
Figure 622326DEST_PATH_IMAGE009
是通过工区先验认识给定的权重系数。
8.如权利要求1中所述的一种基于方位弹性阻抗的HTI型裂缝储层参数预测方法，其特征在于，所预测的HTI型裂缝储层参数包括准裂缝密度参数和脆性指数。
9.如权利要求1中所述的一种基于方位弹性阻抗的HTI型裂缝储层参数预测方法，其特征在于，所述脆性指数至少包括全工区参数的最小值和最大值。
10.一种基于方位弹性阻抗的HTI型裂缝储层参数预测系统，其特征在于，包括：
获取模块，被配置为获取HTI型裂缝储层的第一刚性矩阵和第二刚性矩阵，计算准裂缝密度参数和PP波反射系数方程；
构建方程模块，被配置为对所得到的PP波反射系数方程进行反演计算，构建方位弹性阻抗方程；
参数预测模块，被配置为选择出以反演目标函数最小的方位弹性阻抗，使用被选择出的方位弹性阻抗预测HTI型裂缝储层参数。
CN202210500271.8A 2022-05-10 2022-05-10 基于方位弹性阻抗的hti型裂缝储层参数预测方法及系统 Active CN114609669B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210500271.8A CN114609669B (zh) 2022-05-10 2022-05-10 基于方位弹性阻抗的hti型裂缝储层参数预测方法及系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210500271.8A CN114609669B (zh) 2022-05-10 2022-05-10 基于方位弹性阻抗的hti型裂缝储层参数预测方法及系统

Publications (2)

Publication Number Publication Date
CN114609669A true CN114609669A (zh) 2022-06-10
CN114609669B CN114609669B (zh) 2022-08-23

Family

ID=81869337

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210500271.8A Active CN114609669B (zh) 2022-05-10 2022-05-10 基于方位弹性阻抗的hti型裂缝储层参数预测方法及系统

Country Status (1)

Country Link
CN (1) CN114609669B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117331123A (zh) * 2023-11-06 2024-01-02 成都理工大学 缝洞密度反演方法及缝洞储层预测系统

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2103964A1 (fr) * 2008-03-21 2009-09-23 Ifp Procédé d'estimation de la densité de fractures dans un milieu rocheux
EP2506039A2 (en) * 2011-03-28 2012-10-03 Conocophillips Company Methods for Seismic Fracture Parameter Estimation and Gas Filled Fracture Identification From Vertical Well Log Data
CN104007462A (zh) * 2014-04-16 2014-08-27 孙赞东 基于衰减各向异性的裂缝预测方法
US20170192118A1 (en) * 2016-01-05 2017-07-06 Schlumerger Technology Corporation Amplitude Inversion on Partitioned Depth Image Gathers Using Point Spread Functions
CN107462937A (zh) * 2017-10-09 2017-12-12 中国石油集团川庆钻探工程有限公司地球物理勘探公司 页岩气储层压后裂缝密度三维地质模型构建方法
CN111897006A (zh) * 2020-07-15 2020-11-06 中国石油大学(华东) 一种基于方位弹性阻抗差异奇异值分解的裂缝密度及方向预测方法及系统与应用
US20210255359A1 (en) * 2020-02-19 2021-08-19 Manzar Fawad Method for estimating rock brittleness from well-log data
CN114048627A (zh) * 2021-11-24 2022-02-15 中国石油大学(华东) 基于贝叶斯反演的页岩储层裂缝及脆性预测方法及系统

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP2103964A1 (fr) * 2008-03-21 2009-09-23 Ifp Procédé d'estimation de la densité de fractures dans un milieu rocheux
EP2506039A2 (en) * 2011-03-28 2012-10-03 Conocophillips Company Methods for Seismic Fracture Parameter Estimation and Gas Filled Fracture Identification From Vertical Well Log Data
CN104007462A (zh) * 2014-04-16 2014-08-27 孙赞东 基于衰减各向异性的裂缝预测方法
US20170192118A1 (en) * 2016-01-05 2017-07-06 Schlumerger Technology Corporation Amplitude Inversion on Partitioned Depth Image Gathers Using Point Spread Functions
CN107462937A (zh) * 2017-10-09 2017-12-12 中国石油集团川庆钻探工程有限公司地球物理勘探公司 页岩气储层压后裂缝密度三维地质模型构建方法
US20210255359A1 (en) * 2020-02-19 2021-08-19 Manzar Fawad Method for estimating rock brittleness from well-log data
CN111897006A (zh) * 2020-07-15 2020-11-06 中国石油大学(华东) 一种基于方位弹性阻抗差异奇异值分解的裂缝密度及方向预测方法及系统与应用
CN114048627A (zh) * 2021-11-24 2022-02-15 中国石油大学(华东) 基于贝叶斯反演的页岩储层裂缝及脆性预测方法及系统

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
JIANLU WU ET AL.: "Effective and efficient approaches for calculating seismic ray velocity and attenuation in viscoelastic anisotropic media", 《GEOPHYSICS》 *
PETER R. MESDAG ET AL.: "Quantitative inversion of azimuthal anisotropy parameters from isotropic techniques", 《THE LEADING EDGE》 *
范白涛等: "基于Born逆散射理论的反射系数参数化及弱度参数反演方法研究与应用", 《地球物理学进展》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117331123A (zh) * 2023-11-06 2024-01-02 成都理工大学 缝洞密度反演方法及缝洞储层预测系统
CN117331123B (zh) * 2023-11-06 2024-04-02 成都理工大学 缝洞密度反演方法及缝洞储层预测系统

Also Published As

Publication number Publication date
CN114609669B (zh) 2022-08-23

Similar Documents

Publication Publication Date Title
EP1820137B1 (en) Integrated anisotropic rock physics model
MX2012010473A (es) Metodos y sistemas para realizar inversion elastica simultanea azimutal.
CN105842732A (zh) 多道稀疏反射系数的反演方法及系统
MXPA06012781A (es) Inversion de forma de onda llena en pila previa en 3-d.
CN105452903A (zh) 用于估计地层各向异性常数的模态色散的多频率反转
Wang et al. Data-driven S-wave velocity prediction method via a deep-learning-based deep convolutional gated recurrent unit fusion network
CA2991904A1 (en) Method and apparatus for unambiguously estimating seismic anisotropy parameters
CN112965103B (zh) 一种多重孔隙储层叠前地震概率化多道反演方法
CN105138731A (zh) 一种水合物分解引起海底斜坡不稳定性评价系统及方法
US11692972B2 (en) Systems and methods for ultrasonic characterization of permafrost, frozen soil and saturated soil samples
Jin et al. Reflection and transmission approximations for weak contrast orthorhombic media
CN114609669B (zh) 基于方位弹性阻抗的hti型裂缝储层参数预测方法及系统
Rwechungura et al. Application of particle swarm optimization for parameter estimation integrating production and time lapse seismic data
CN115993649A (zh) 基于等效方位杨氏模量的裂缝参数预测方法及系统
Chaikine et al. A new machine learning procedure to generate highly accurate synthetic shear sonic logs in unconventional reservoirs
Yuan et al. Quantitative uncertainty evaluation of seismic facies classification: A case study from northeast China
Bazulin et al. Determination of the elastic parameters of a VTI medium from sonic logging data using deep learning
US20230125277A1 (en) Integration of upholes with inversion-based velocity modeling
US11947062B2 (en) Velocity tomography using time lags of wave equation migration
CN115586572A (zh) 一种孔隙参数与储层参数的地震岩石物理解析反演方法
US10261205B2 (en) System for processing seismic data based upon volatility measurement model and related methods
Hajian et al. Genetic algorithm with applications in geophysics
Bayuk Why anisotropy is important for location of microearthquake events in shale?
CN111077573A (zh) 一种确定地层弹性参数的方法、装置及系统
US20230408718A1 (en) Method and system for generalizable deep learning framework for seismic velocity estimation robust to survey configuration

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