CN111239805B - 基于反射率法的块约束时移地震差异反演方法及系统 - Google Patents
基于反射率法的块约束时移地震差异反演方法及系统 Download PDFInfo
- Publication number
- CN111239805B CN111239805B CN202010091778.3A CN202010091778A CN111239805B CN 111239805 B CN111239805 B CN 111239805B CN 202010091778 A CN202010091778 A CN 202010091778A CN 111239805 B CN111239805 B CN 111239805B
- Authority
- CN
- China
- Prior art keywords
- inversion
- difference
- seismic
- time
- lapse seismic
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 139
- 238000002310 reflectometry Methods 0.000 title claims abstract description 57
- 238000004364 calculation method Methods 0.000 claims abstract description 27
- 238000004088 simulation Methods 0.000 claims abstract description 12
- 238000012545 processing Methods 0.000 claims abstract description 11
- 238000006243 chemical reaction Methods 0.000 claims abstract description 8
- 238000009499 grossing Methods 0.000 claims abstract description 8
- 230000006870 function Effects 0.000 claims description 32
- 238000009826 distribution Methods 0.000 claims description 18
- 230000008569 process Effects 0.000 claims description 14
- 239000011159 matrix material Substances 0.000 claims description 11
- 238000004590 computer program Methods 0.000 claims description 10
- 238000003860 storage Methods 0.000 claims description 9
- 238000012360 testing method Methods 0.000 claims description 6
- 239000000126 substance Substances 0.000 claims description 5
- 238000005315 distribution function Methods 0.000 claims description 4
- 238000000605 extraction Methods 0.000 claims description 3
- 238000005457 optimization Methods 0.000 claims description 3
- 238000001228 spectrum Methods 0.000 claims description 2
- 230000008859 change Effects 0.000 abstract description 4
- 239000010410 layer Substances 0.000 description 14
- 239000003921 oil Substances 0.000 description 11
- 230000000694 effects Effects 0.000 description 10
- 238000012544 monitoring process Methods 0.000 description 9
- 230000004044 response Effects 0.000 description 7
- 238000010586 diagram Methods 0.000 description 6
- 239000004743 Polypropylene Substances 0.000 description 5
- 230000005540 biological transmission Effects 0.000 description 5
- 238000004519 manufacturing process Methods 0.000 description 5
- 238000005516 engineering process Methods 0.000 description 4
- 238000002347 injection Methods 0.000 description 4
- 239000007924 injection Substances 0.000 description 4
- 239000011229 interlayer Substances 0.000 description 4
- 239000000243 solution Substances 0.000 description 4
- 230000008901 benefit Effects 0.000 description 3
- 238000009795 derivation Methods 0.000 description 3
- 238000012546 transfer Methods 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 230000007704 transition Effects 0.000 description 2
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 2
- 238000010521 absorption reaction Methods 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 230000000903 blocking effect Effects 0.000 description 1
- 238000004422 calculation algorithm Methods 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 239000010779 crude oil Substances 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 239000006185 dispersion Substances 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 239000000463 material Substances 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000005065 mining Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000000704 physical effect Effects 0.000 description 1
- -1 polypropylene Polymers 0.000 description 1
- 229920001155 polypropylene Polymers 0.000 description 1
- 238000011084 recovery Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 238000000638 solvent extraction Methods 0.000 description 1
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/282—Application of seismic models, synthetic seismograms
-
- 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/303—Analysis for determining velocity profiles or travel times
-
- 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/622—Velocity, density or impedance
- G01V2210/6222—Velocity; travel time
-
- 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/622—Velocity, density or impedance
- G01V2210/6224—Density
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
本发明公开了一种基于反射率法的块约束时移地震差异反演方法及系统，该方法包括：利用反射率法得到时移地震差异反演中的正演算子；利用测井数据来正演模拟角度道集进行井震标定，与井旁地震道进行相关对比，确定时深转换关系，拾取层位信息；基于线性褶积模型对实际差异地震记录分角度提取子波；结合测井数据建立全频段插值结果模型，并通过平滑处理获取差异反演的低频初始模型；基于贝叶斯反演推理框架构建差异反演的目标函数，得到正演模拟记录和实际记录的残差，并在贝叶斯框架中加入优化的超拉普拉斯块约束，建立差异反演公式；进行时移地震差异反演，经过迭代计算，当得到的计算结果变化幅度小于一阈值时终止迭代，得到时移地震反演结果。
Description
技术领域
本发明涉及时移地震油藏监测技术领域，尤指一种基于反射率法的块约束时移地震差异反演方法、系统、计算机设备及计算机可读存储介质，用于时移地震中的高精度、高分辨率储层差异弹性参数估计。
背景技术
时移地震油藏监测技术已经被成功地应用于监测油气生产过程中的油藏变化，从而用以寻找死油区、确定新井位和优化注采方案，提高了油藏采收率。但是，油藏在生产过程中会伴随注水或注气开采，一方面由于原油的采出和注入水或气，改变了油藏的含油饱和度而产生了地震响应的变化，另一方面，油藏在开采过程中其压力系统也发生了变化，这可能引起可监测的地震振幅差异。
在时移地震油藏监测技术中，反演是获取储层差异弹性参数的重要手段，其主要任务是直接利用差异地震数据估计差异弹性参数，为储层建模提供基础约束信息，进而获取储层岩相、物性分布信息。大多数的时移地震都是运用分别反演的方法，简单易求。但往往由于不同年份时移地震数据间的反演没有进行耦合，而导致不正确的解释结果。与分别反演方法相比，差异反演在反演过程中直接利用差异地震数据进行计算，能够提供必要的耦合条件，减少计算量。
在时移地震差异反演过程中，由于叠后数据往往假设地震数据采集是自激自收的，振幅不随偏移距变化。实际上这种假设与实际多次覆盖观测系统不符，而且损失了隐藏在叠前资料中的丰富信息。与叠后差异地震数据相比，叠前差异地震数据携带了丰富的地下介质的岩性信息和流体信息。因此为了能够更好地表征油藏在开采过程中的变化情况，需要运用叠前差异地震数据进行计算。
在现有技术中，常规的叠前时移地震差异反演的技术存在以下缺陷：
1、该方法往往运用Zoeppritz方程的近似式进行计算，简化了求解的过程，但该方法是小角度线性近似的，因此在大角度(例如大于30°)的情况下，计算的反射系数误差较大。
2、该方法的计算过程基于单界面假设仅计算一次反射波信息，忽略透射损失、层间多次波和层间转换波的影响，这些会导致最终差异反演的结果存在一定误差。
3、该方法基于弹性假设，忽略衰减与频散影响，实际地下介质并非完全弹性，地震波的衰减效应将会导致子波振幅减小以及波形畸变，降低地震记录的分辨率。
4、反演获得的差异弹性参数往往是块状的，该方法对块状边界的刻画不够清晰，在高斯以及拉普拉斯约束下很难得到稀疏性较高的差异反演结果。
5、调谐干涉影响反演结果的分辨率，由于传统方法基于单界面假设推导，故反演分辨率受到界面调谐影响；即使在时移地震差异反演过程中运用精确Zoeppritz方程进行计算，也会由于为考虑透射损失较多次波信息，而导致差异反演中的正演算子无法准确求取，从而影响差异反演的结果。
综上来看，在现有的时移地震差异反演过程中，由于差异弹性参数的特征、小角度线性近似、不完全波形模拟以及忽略衰减效应等情况，会导致常规差异反演结果准确度不够、分辨率降低；因此，亟需一种可以解决上述问题的时移地震差异反演方法，实现更为准确以及分辨率更高的反演结果，满足实际油藏地震储层表征需求。
发明内容
为解决上述问题，本发明提出了一种基于反射率法的块约束时移地震差异反演方法，可以基于反射率法(一维波动方程解析解)，进行时移地震差异反演中正演算子的求取，并利用测井数据来正演模拟角度道集进行井震标定，进一步拾取层位信息、提取子波、建立初始模型，基于贝叶斯差异反演推理框架构建目标函数，在贝叶斯框架中加入优化的超拉普拉斯块约束，来更加精细地刻画差异弹性参数的块状边界，从而获取有效准确的反演结果，为实际油藏地震储层表征提供有力的技术支持。
在本发明一实施例中，提出了一种基于反射率法的块约束时移地震差异反演方法，该方法包括：
获取测井数据；
根据所述测井数据，利用反射率法得到时移地震差异反演中的正演算子；
根据所述时移地震差异反演中的正演算子，利用所述测井数据来正演模拟角度道集进行井震标定，与井旁地震道进行相关对比，确定时深转换关系，拾取层位信息；
基于线性褶积模型对实际差异地震记录分角度提取子波；
根据拾取的所述层位信息及提取的所述子波，结合测井数据建立全频段插值结果模型，并通过平滑处理获取差异反演的低频初始模型；
根据所述低频初始模型，基于贝叶斯反演推理框架构建差异反演的目标函数，得到正演模拟记录和实际记录的残差，并在贝叶斯框架中加入优化的超拉普拉斯块约束，建立差异反演公式；
利用所述差异反演公式进行时移地震差异反演，经过迭代计算，当得到的计算结果变化幅度小于一阈值时终止迭代，得到时移地震反演结果。
在本发明另一实施例中，还提出了一种基于反射率法的块约束时移地震差异反演系统，该系统包括：
数据获取模块，用于获取测井数据；
正演算子计算模块，用于根据所述测井数据，利用反射率法得到时移地震差异反演中的正演算子；
层位信息拾取模块，用于根据所述时移地震差异反演中的正演算子，利用所述测井数据来正演模拟角度道集进行井震标定，与井旁地震道进行相关对比，确定时深转换关系，拾取层位信息；
子波提取模块，用于基于线性褶积模型对实际差异地震记录分角度提取子波；
模型建立模块，用于根据拾取的所述层位信息及提取的所述子波，结合测井数据建立全频段插值结果模型，并通过平滑处理获取差异反演的低频初始模型；
差异反演公式建立模块，用于根据所述低频初始模型，基于贝叶斯反演推理框架构建差异反演的目标函数，得到正演模拟记录和实际记录的残差，并在贝叶斯框架中加入优化的超拉普拉斯块约束，建立差异反演公式；
差异反演迭代计算模块，用于利用所述差异反演公式进行时移地震差异反演，经过迭代计算，当得到的计算结果变化幅度小于一阈值时终止迭代，得到时移地震反演结果。
在本发明另一实施例中，还提出了一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现基于反射率法的块约束时移地震差异反演方法。
在本发明另一实施例中，还提出了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现基于反射率法的块约束时移地震差异反演方法。
本发明提出的基于反射率法的块约束时移地震差异反演方法、系统、计算机设备及计算机可读存储介质，基于反射率法，进行时移地震差异反演中正演算子的求取，并利用测井数据来正演模拟角度道集进行井震标定，进一步拾取层位信息、提取子波、建立初始模型，基于贝叶斯差异反演推理框架构建目标函数，在贝叶斯框架中加入优化的超拉普拉斯块约束，来更加精细地刻画差异弹性参数的块状边界，从而得到更为准确以及分辨率更高的差异反演结果，为实际油藏地震储层表征提供有力的数据支持。
附图说明
图1是本发明一实施例的基于优化的超拉普拉斯块约束的反射率法时移地震差异反演方法流程示意图。
图2是本发明一实施例的基于优化的超拉普拉斯块约束的反射率法时移地震差异反演系统架构示意图。
图3是本发明一实施例的计算机设备结构示意图。
图4A至图4C分别是基础数据的正演结果、监测数据的正演结果及差异数据的正演结果示意图。
图5A及图5B分别为高斯分布及优化的超拉普拉斯分布的度量示意图。
图6A及图6B分别为基于高斯约束的反射率法时移地震差异反演方法在信噪比为1000和8时得到的纵波速度、横波速度及密度示意图。
图7A及图7B分别为基于优化的超拉普拉斯块约束的反射率法时移地震差异反演方法在信噪比为1000和8时得到的纵波速度、横波速度及密度示意图。
具体实施方式
下面将参考若干示例性实施方式来描述本发明的原理和精神。应当理解，给出这些实施方式仅仅是为了使本领域技术人员能够更好地理解进而实现本发明，而并非以任何方式限制本发明的范围。相反，提供这些实施方式是为了使本公开更加透彻和完整，并且能够将本公开的范围完整地传达给本领域的技术人员。
本领域技术人员知道，本发明的实施方式可以实现为一种系统、装置、设备、方法或计算机程序产品。因此，本公开可以具体实现为以下形式，即：完全的硬件、完全的软件(包括固件、驻留软件、微代码等)，或者硬件和软件结合的形式。
根据本发明的实施方式，提出了一种基于优化的超拉普拉斯块约束的反射率法时移地震差异反演方法、系统、计算机设备及计算机存储介质。其中，该方法主要是基于反射率法，进行时移地震差异反演中正演算子的求取，并利用测井数据来正演模拟角度道集进行井震标定，进一步拾取层位信息、提取子波、建立初始模型，基于贝叶斯差异反演推理框架构建目标函数，在贝叶斯框架中加入优化的超拉普拉斯块约束，来更加精细地刻画差异弹性参数的块状边界，从而获取有效准确的反演结果。
下面参考本发明的若干代表性实施方式，详细阐释本发明的原理和精神。
图1是本发明一实施例的基于优化的超拉普拉斯块约束的反射率法时移地震差异反演方法流程示意图。如图1所示，该方法包括：
步骤S1，获取测井数据；
步骤S2，根据所述测井数据，利用反射率法(一维波动方程解析解)得到时移地震差异反演中的正演算子；
步骤S3，根据所述时移地震差异反演中的正演算子，利用所述测井数据来正演模拟角度道集进行井震标定，与井旁地震道进行相关对比，确定时深转换关系，拾取层位信息；由于考虑透射效应、层间多次、衰减效应等影响，合成地震记录与井旁道地震记录更加匹配。
步骤S4，基于线性褶积模型对实际差异地震记录分角度提取子波；
步骤S5，根据拾取的所述层位信息及提取的所述子波，结合测井数据建立全频段插值结果模型，并通过平滑处理获取差异反演的低频初始模型；
步骤S6，根据所述低频初始模型，基于贝叶斯反演推理框架构建差异反演的目标函数，得到正演模拟记录和实际记录的残差，并在贝叶斯框架中加入优化的超拉普拉斯块约束，建立差异反演公式；
步骤S7，利用所述差异反演公式进行时移地震差异反演，经过迭代计算，当得到的计算结果变化幅度小于一阈值时终止迭代，得到时移地震反演结果。
在一实施例中，该方法还包括：
步骤S8，根据所述时移地震反演结果进行井旁地震道反演测试，对比所述时移地震反演结果与测井结果，根据对比结果优化反演参数；其中，所述反演参数至少包括：最大迭代次数、角度范围及超参数；
步骤S9，根据优化后的反演参数对工区地震数据进行并行化时移地震反演，得到工区的时移地震反演结果。
在一实施例中，步骤S4具体的处理过程包括：基于线性褶积模型，假设地下反射系数是具有白噪谱的随机序列，利用统计学原理提取子波，由于子波随角度是变化的，因此需要对实际差异地震记录分角度提取子波，井旁地震模拟匹配确定振幅缩放因子，以差异测井数据为输入模型利用反射率法模拟角度域的PP波道集，与实际井旁角度域地震道集对比，计算振幅缩放因子，并应用于所提取的地震子波。
在一实施例中，步骤S5具体的处理过程包括：
步骤S51，根据所述低频初始模型提取模型参数的均值，假设模型参数服从高斯分布，得到弹性参数的自相关系数以及互相关系数；
步骤S52，根据所述弹性参数的自相关系数以及互相关系数得到所述弹性参数的协方差矩阵，建立服从该工区统计特征的模型参数的先验分布函数。
在一实施例中，步骤S6具体的处理过程包括：
步骤S61，基于贝叶斯反演推理框架，根据所述先验分布函数构建差异反演的目标函数，计算正演模拟记录和实际记录的残差；
步骤S62，根据所述残差，在正演过程中运用泰勒展开略去高阶项，得到反射率法的差异正演线性近似公式；
步骤S63，根据所述反射率法的差异正演线性近似公式，在贝叶斯框架中加入优化的超拉普拉斯块约束，建立差异反演公式。
其中，建立的差异反演公式为：
其中，
为差异反演结果；
W为差异反演中的正演算子；
Ce为噪音的协方差矩阵；
CΔm为差异弹性参数的协方差矩阵；
D为一阶微分算子；
Δd为基础数据与监测数据之差；
β为约束项的权重；
μ为差异数据的均值向量；
进一步的，
算子为：
其中，κj(j＝1,2,3)分别表示纵波速度、横波速度及密度的尺度因子；
e为常数因子；
Δm为基础数据与监测数据对应的弹性参数之差；
p值根据不同井数据特征进行确定。
为了对上述基于优化的超拉普拉斯块约束的时移地震差异反演方法进行更为清楚的解释，下面结合一个具体的实施例来进行说明。
针对导致常规时移地震差异反演结果准确度不够、分辨率降低等问题，本发明提出了基于优化的超拉普拉斯块约束的反射率法时移地震差异反演方法，可得到更为准确以及分辨率更高的差异反演结果，时移地震中基础数据与监测数据的正演过程如下：
首先，时移地震中基础数据与监测数据的正演过程如下：
d1＝G(m1)+e1＝R(m1)*S+e1； (1)
d2＝G(m2)+e2＝R(m2)*S+e2； (2)
其中，d1和d2分别是基础数据和监测数据；
m1和m2分别是两期数据对应的弹性参数；
e1和e2分别是两期数据所含噪音；
G是正演算子；
R是反射系数；
S是地震记录的子波。
为了便于求解计算，需将反射率法的正演过程进行线性近似，来获取差异正演公式。在假设两期数据变化较小的前提下，油藏开采前后才会具有相同的地质背景；可以用同一个模型来获取背景响应，所以利用泰勒展开能够将式1和式2改写成：
其中，式3及式4中的(“≈”号后)第一项表示的是两期数据中相同的背景响应，中间项相当于两期数据在这一相同背景响应下的不同扰动，最后一项表示数据中的噪声。将两期正演结果相减，得到差异地震记录(式5)，用W表示差异反演中的正演算子：
运用反射率法求解正演算子时，需要假设地下介质是水平层状的，且有N层。反射率法主要是从底界面开始，运用传递算子来逐层向上递推，获得顶界面以下的总体反射相应。用向量v来表示界面之下频率慢度域的总反射响应：
v＝[Δ -RSPΔ -RSSΔ RPPΔ RPSΔ |R|Δ]T； (6)
其中，Δ表示缩放因子；
R表示频率慢度域(ω-p)反射系数，其中第一个下标表示入射波类型，第二个下标是反射波类型；如，SP表示入射波为S波(横波)，反射波为P波(纵波)；
|R|表示反射系数的行列式。
接着运用传递算子Q来逐层向上递推，得到上层的总反射响应：
其中，Qn表示第n层的传递算子；
En表示第n层介质的相位变化；
介质的底界面以下是弹性半空间，因此在该处只有透射传播没有反射，所以运用式7，从底界面逐层向上递推得到顶界面的总响应v0：
v0＝Q0Q1…QN-1vN； (8)
根据式8计算得到的向量v0，可以获得PP波在频率慢度域(ω-p)的反射系数：
为了求解差异正演中的算子W，需要求反射系数的一阶偏导：
其中，
分别表示第n层的纵波速度、横波速度和密度。进一步的，式10中的部分偏导如下：
在进行叠前差异反演时，需要运用角道集数据，因此需要将单期数据的正演算子G及差异数据的正演算子W转至时间角度域(τ-θ)。
首先，运用
将慢度域反射系数转至角度域，接着再进行快速傅里叶逆变换(IFFT)：
贝叶斯框架方法中，差异弹性参数Δm的后验分布如下：
其中，P(Δd|Δm)是从弹性参数到观测数据的似然函数；
P(Δm)是先验分布，可以引入与弹性参数相关的信息来进行约束反演，提高稳定性及分辨率；
P(Δd)是边缘概率密度，是一个常数。
假设正演式5中的噪音项e服从零均值高斯分布，那么差异反演中的似然函数可以表示为：
其中，Ce是噪音的协方差矩阵，假设数据中的噪音是不相关的，可以将Ce转化为对角矩阵(Ce＝σ2I)，其中σ2是噪音的方差，I是单位阵，大小为Nd×Nd，Nd是观测数据的长度。
由于时移地震差异数据中存在分块现象，为了能够更好地刻画差异数据的地层边界，需要在先验约束中再加入有着长尾巴特征的块约束项。因此，先验约束由两部分组成：1、包含先验低频趋势和不同模型参数之间协方差的高斯分布项PG(Δm)；2、具有长尾巴分布特征的块约束项PB(Δm)。其中，高斯分布项可以表示为：
其中，CΔm是差异弹性参数的协方差矩阵；
N是模型参数的长度；
μ是差异数据的均值向量。
根据式13，求出差异弹性参数的后验分布后，对其求负对数，并舍去常数项，可以得到目标函数J(Δm)。因为要使后验概率最大，所以目标函数需要达到最小，才能求出恰当的差异弹性参数：
min(J(Δm))＝min(L(Δd|Δm)+β(LG(Δm)+LB(Δm)))； (16)
其中，β为约束项的权重；似然函数的负对数L(Δd|Δm)，高斯分布项的负对数LG(Δm)和块约束项的负对数LB(Δm)可以表示为：
其中，D为一阶微分算子，可以表示为
κj(j＝1,2,3)分别表示纵波速度、横波速度及密度的尺度因子，通过测试求得尺度因子的值；
块约束函数B(x)为正则化项。
将目标函数关于差异弹性参数Δm求导，并令导数为零(式18)，可以求出最小目标函数所对应的差异弹性参数：
其中，
算子是块约束函数B(x)的部分导数项。
式19是求解贝叶斯框架下最优反演结果的方程，为了获得最终的反演结果，需要应用迭代重加权方法。当求解的结果不再发生大幅度变化时，迭代将会终止，即可得到最终的反演结果。
在上述反演过程中引入了正则化项：块约束函数B(x)，能够增强梯度的稀疏性，可以更好地刻画地层边界。因此，块约束函数的选取尤为重要，需要满足三个特性：较好的稀疏性，求解便利以及稳定性好。
超拉普拉斯约束
是在拉普拉斯分布的基础上，通过自定义α值，来选取具有较好稀疏性的结果，曾被应用在弱故障特征增强领域，并取得了很好的效果。
将超拉普拉斯约束经过负对数求取后，由于α选取部分数值时，函数具有非凸性的，导致后续的反演求解较为困难。在此基础上进行了改进，提出了优化的超拉普拉斯约束来解决一阶求导的问题：
其中，e是常数因子，一般取很小的数值，例如可以取0.0001。
式21中的p值选取，对差异弹性参数的反演有一定影响，反演前首先需要对不同的p值进行测试，选取最适合工区数据的参数值。为了突出差异数据的边界，需要在反演结果的垂向梯度上加强稀疏性，因此需要对差异弹性参数进行一阶偏导的概率密度分布统计，从而确定合适的p值。
综合上述的分析，可以将优化的超拉普拉斯分布(式21)引入到块约束函数B(x)。在目标函数中对块约束函数进行求导，来得到算子
(式22)：
其中，p值得选取需要根据不同井数据特征来确定，针对本示例中给出的模型，可以统计出最合适的值是p＝0.2。
进一步，可以结合步骤S7，利用所述差异反演公式(式19)进行时移地震差异反演，经过迭代计算，当得到的计算结果变化幅度小于一阈值时终止迭代，得到时移地震反演结果。
另外，结合步骤S8，可以根据所述时移地震反演结果进行井旁地震道反演测试，对比所述时移地震反演结果与测井结果，根据对比结果优化反演参数；其中，所述反演参数至少包括：最大迭代次数、角度范围及超参数；
最后，可以根据优化后的反演参数向整个工区进行推广，对工区地震数据进行并行化时移地震反演，得到整个工区的时移地震反演结果。
本发明由于采取以上技术方案，至少具有以下优点：
1、相比于现有的叠前AVO技术，本技术方案能够提供一种高可靠性、高分辨率的弹性参数估计方法。
2、本技术方案提供一个解析求解相应正演方程导数矩阵方法，不再使用数值求解方法，使得基于该类方法反演方法效率大幅度提高。
3、本技术方案对AVO反演方法的补充和扩展，利用相对成熟的AVO角道集优化技术作为算法的输入支持的同时，通过波动方程解析解计算全波场信息，对应的反演结果具有较高的分辨率，可满足储层预测的精度要求。
4、本技术方案不再基于线性近似，适用于任何角度范围道集，为大角度地震数据反演提供坚实的理论基础。
5、本技术方案是模拟全波场信息，包括一次反射波之外的透射损失、层间多次波、层间转换波等信息，充分考虑层厚影响，消除薄层的干涉作用，进而可获取高精度的反演结果。
6、本技术方案将衰减补偿与叠前波形反演方法结合，在反演过程中考虑吸收衰减效应，进而获取更为准确以及分辨率更高的反演结果。
需要说明的是，尽管在上述实施例及附图中以特定顺序描述了本发明方法的操作，但是，这并非要求或者暗示必须按照该特定顺序来执行这些操作，或是必须执行全部所示的操作才能实现期望的结果。附加地或备选地，可以省略某些步骤，将多个步骤合并为一个步骤执行，和/或将一个步骤分解为多个步骤执行。
在介绍了本发明示例性实施方式的方法之后，接下来，参考图2对本发明示例性实施方式的基于反射率法的块约束时移地震差异反演系统进行介绍。
基于反射率法的块约束时移地震差异反演系统的实施可以参见上述方法的实施，重复之处不再赘述。以下所使用的术语“模块”或者“单元”，可以是实现预定功能的软件和/或硬件的组合。尽管以下实施例所描述的装置较佳地以软件来实现，但是硬件，或者软件和硬件的组合的实现也是可能并被构想的。
基于同一发明构思，本发明还提出了一种基于反射率法的块约束时移地震差异反演系统，如图2所示，该系统包括：
数据获取模块201，用于获取测井数据；
正演算子计算模块202，用于根据所述测井数据，利用反射率法得到时移地震差异反演中的正演算子；
层位信息拾取模块203，用于根据所述时移地震差异反演中的正演算子，利用所述测井数据来正演模拟角度道集进行井震标定，与井旁地震道进行相关对比，确定时深转换关系，拾取层位信息；
子波提取模块204，用于基于线性褶积模型对实际差异地震记录分角度提取子波；
模型建立模块205，用于根据拾取的所述层位信息及提取的所述子波，结合测井数据建立全频段插值结果模型，并通过平滑处理获取差异反演的低频初始模型；
差异反演公式建立模块206，用于根据所述低频初始模型，基于贝叶斯反演推理框架构建差异反演的目标函数，得到正演模拟记录和实际记录的残差，并在贝叶斯框架中加入优化的超拉普拉斯块约束，建立差异反演公式；
差异反演迭代计算模块207，用于利用所述差异反演公式进行时移地震差异反演，经过迭代计算，当得到的计算结果变化幅度小于一阈值时终止迭代，得到时移地震反演结果；
进一步，该系统还可以包括：参数优化模块208，用于根据所述时移地震反演结果进行井旁地震道反演测试，对比所述时移地震反演结果与测井结果，根据对比结果优化反演参数；其中，所述反演参数至少包括：最大迭代次数、角度范围及超参数；
最后，通过差异反演迭代计算模块207，根据优化后的反演参数对工区地震数据进行并行化时移地震反演，得到整个工区的时移地震反演结果。
应当注意，尽管在上文详细描述中提及了基于反射率法的块约束时移地震差异反演系统的若干模块，但是这种划分仅仅是示例性的并非强制性的。实际上，根据本发明的实施方式，上文描述的两个或更多模块的特征和功能可以在一个模块中具体化。反之，上文描述的一个模块的特征和功能可以进一步划分为由多个模块来具体化。
基于前述发明构思，如图3所示，本发明还提出了一种计算机设备300，包括存储器310、处理器320及存储在存储器310上并可在处理器320上运行的计算机程序330，所述处理器320执行所述计算机程序330时实现前述基于优化的超拉普拉斯块约束的时移地震差异反演方法。
基于前述发明构思，本发明还提出了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现基于反射率法的块约束时移地震差异反演方法。
为了对上述基于反射率法的块约束时移地震差异反演方法、系统、计算机设备及计算机可读存储介质进行更为清楚的解释，下面结合一个具体的实施例来进行说明，然而值得注意的是该实施例仅是为了更好地说明本发明，并不构成对本发明不当的限定。
将基于精确Zoeppritz方程的时移地震正演模型与(本发明利用的)反射率法的时移地震正演模型进行对比，其对比结果如图4A至图4C所示。其中，图4A是基础数据的正演结果，图4B是监测数据的正演结果，图4C是差异数据的正演结果。在每幅图中包含三部分，左侧部分是反射率法的正演结果，中间部分是精确Zoeppritz方法的正演结果，右侧部分是两种方法正演结果的差异。通过对比可知，相较于精确Zoeppritz方程的时移地震正演模型，可能会导致差异反演中的正演算子无法准确求取，从而影响差异反演的结果，而利用反射率法的时移地震正演模型则可避免该些问题，能够准确求取差异反演中的正演算子。
进一步的结合图5A及图5B所示，分别为高斯分布以及优化的超拉普拉斯分布的度量图。结合图5A及图5B，以信噪比1000和8为例，分别进行基于高斯约束的反射率法时移地震差异反演及基于优化的超拉普拉斯块约束的反射率法时移地震差异反演。
如图6A及图6B所示，分别为基于高斯约束的反射率法时移地震差异反演方法在信噪比为1000和8时得到的纵波速度(Vp)、横波速度(Vs)、密度。其中，图6A是信噪比为1000时的反演结果，图6B是信噪比为8时的反演结果；图中，线601是真实差异井数据，线602是差异反演结果，虚线603是初始模型；
如图7A及图7B所示，分别为基于优化的超拉普拉斯块约束的反射率法时移地震差异反演方法在信噪比为1000和8时得到的纵波速度(Vp)、横波速度(Vs)、密度。其中，图7A是信噪比为1000时的反演结果，图7B是信噪比为8时的反演结果；图中，线701是真实差异井数据，线702是差异反演结果，虚线703是初始模型。
通过对比上述两种方法得到的反演结果可以看出，图6A及图6B得到的反演结果准确度不够、分辨率较低；再结合图7A及图7B来看，利用本发明提出的基于优化的超拉普拉斯块约束的反射率法时移地震差异反演方法，能够获得准确以及分辨率更高的反演结果，满足实际油藏地震储层表征需求。
综上来看，本发明提出的基于反射率法的块约束时移地震差异反演方法、系统、计算机设备及计算机可读存储介质基于反射率法，进行时移地震差异反演中正演算子的求取，并利用测井数据来正演模拟角度道集进行井震标定，进一步拾取层位信息、提取子波、建立初始模型，基于贝叶斯差异反演推理框架构建目标函数，在贝叶斯框架中加入优化的超拉普拉斯块约束，来更加精细地刻画差异弹性参数的块状边界，从而得到更为准确以及分辨率更高的差异反演结果，为实际油藏地震储层表征提供有力的数据支持。
虽然已经参考若干具体实施方式描述了本发明的精神和原理，但是应该理解，本发明并不限于所公开的具体实施方式，对各方面的划分也不意味着这些方面中的特征不能组合以进行受益，这种划分仅是为了表述的方便。本发明旨在涵盖所附权利要求的精神和范围内所包括的各种修改和等同布置。
Claims (10)
1.一种基于反射率法的块约束时移地震差异反演方法，其特征在于，该方法包括：
获取测井数据；
根据所述测井数据，利用反射率法得到时移地震差异反演中的正演算子；
根据所述时移地震差异反演中的正演算子，利用所述测井数据来正演模拟角度道集进行井震标定，与井旁地震道进行相关对比，确定时深转换关系，拾取层位信息；
基于线性褶积模型对实际差异地震记录分角度提取子波；
根据拾取的所述层位信息及提取的所述子波，结合测井数据建立全频段插值结果模型，并通过平滑处理获取差异反演的低频初始模型；
根据所述低频初始模型，基于贝叶斯反演推理框架构建差异反演的目标函数，得到正演模拟记录和实际记录的残差，并在贝叶斯框架中加入优化的超拉普拉斯块约束，建立差异反演公式；
利用所述差异反演公式进行时移地震差异反演，经过迭代计算，当得到的计算结果变化幅度小于一阈值时终止迭代，得到时移地震反演结果。
2.根据权利要求1所述的基于反射率法的块约束时移地震差异反演方法，其特征在于，该方法还包括：
根据所述时移地震反演结果进行井旁地震道反演测试，对比所述时移地震反演结果与测井结果，根据对比结果优化反演参数；其中，所述反演参数至少包括：最大迭代次数、角度范围及超参数；
根据优化后的反演参数对工区地震数据进行并行化时移地震反演，得到工区的时移地震反演结果。
3.根据权利要求1所述的基于反射率法的块约束时移地震差异反演方法，其特征在于，基于线性褶积模型对实际差异地震记录分角度提取子波，包括：
基于线性褶积模型，假设地下反射系数是具有白噪谱的随机序列，利用统计学原理对实际差异地震记录分角度提取子波。
4.根据权利要求1所述的基于反射率法的块约束时移地震差异反演方法，其特征在于，根据拾取的所述层位信息及提取的所述子波，结合测井数据建立全频段插值结果模型，并通过平滑处理获取差异反演的低频初始模型，包括：
根据所述低频初始模型提取模型参数的均值，假设模型参数服从高斯分布，得到弹性参数的自相关系数以及互相关系数；
根据所述弹性参数的自相关系数以及互相关系数得到所述弹性参数的协方差矩阵，建立服从工区统计特征的模型参数的先验分布函数。
5.根据权利要求4所述的基于反射率法的块约束时移地震差异反演方法，其特征在于，根据所述低频初始模型，基于贝叶斯反演推理框架构建差异反演的目标函数，得到正演模拟记录和实际记录的残差，并在贝叶斯框架中加入优化的超拉普拉斯块约束，建立差异反演公式，包括：
基于贝叶斯反演推理框架，根据所述先验分布函数构建差异反演的目标函数，计算正演模拟记录和实际记录的残差；
根据所述残差，在正演过程中运用泰勒展开略去高阶项，得到反射率法的差异正演线性近似公式；
根据所述反射率法的差异正演线性近似公式，在贝叶斯框架中加入优化的超拉普拉斯块约束，建立差异反演公式。
6.根据权利要求5所述的基于反射率法的块约束时移地震差异反演方法，其特征在于，根据所述反射率法的差异正演线性近似公式，在贝叶斯框架中加入优化的超拉普拉斯块约束，建立差异反演公式，其中，
差异反演公式为：
其中，
为差异反演结果；
W为差异反演中的正演算子；
Ce为噪音的协方差矩阵；
CΔm为差异弹性参数的协方差矩阵；
D为一阶微分算子；
Δd为基础数据与监测数据之差；
β为约束项的权重；
μ为差异数据的均值向量；
其中，
算子为：
其中，κj(j＝1,2,3)分别表示纵波速度、横波速度及密度的尺度因子；
e为常数因子；
Δm为基础数据与监测数据对应的弹性参数之差；
p值根据不同井数据特征进行确定。
7.一种基于反射率法的块约束时移地震差异反演系统，其特征在于，该系统包括：
数据获取模块，用于获取测井数据；
正演算子计算模块，用于根据所述测井数据，利用反射率法得到时移地震差异反演中的正演算子；
层位信息拾取模块，用于根据所述时移地震差异反演中的正演算子，利用所述测井数据来正演模拟角度道集进行井震标定，与井旁地震道进行相关对比，确定时深转换关系，拾取层位信息；
子波提取模块，用于基于线性褶积模型对实际差异地震记录分角度提取子波；
模型建立模块，用于根据拾取的所述层位信息及提取的所述子波，结合测井数据建立全频段插值结果模型，并通过平滑处理获取差异反演的低频初始模型；
差异反演公式建立模块，用于根据所述低频初始模型，基于贝叶斯反演推理框架构建差异反演的目标函数，得到正演模拟记录和实际记录的残差，并在贝叶斯框架中加入优化的超拉普拉斯块约束，建立差异反演公式；
差异反演迭代计算模块，用于利用所述差异反演公式进行时移地震差异反演，经过迭代计算，当得到的计算结果变化幅度小于一阈值时终止迭代，得到时移地震反演结果。
8.根据权利要求7所述的基于反射率法的块约束时移地震差异反演系统，其特征在于，该系统还包括：
参数优化模块，用于根据所述时移地震反演结果进行井旁地震道反演测试，对比所述时移地震反演结果与测井结果，根据对比结果优化反演参数；其中，所述反演参数至少包括：最大迭代次数、角度范围及超参数；
所述差异反演迭代计算模块，还用于根据优化后的反演参数对工区地震数据进行并行化时移地震反演，得到工区的时移地震反演结果。
9.一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1至6任一所述方法。
10.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现权利要求1至6任一所述方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010091778.3A CN111239805B (zh) | 2020-02-13 | 2020-02-13 | 基于反射率法的块约束时移地震差异反演方法及系统 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010091778.3A CN111239805B (zh) | 2020-02-13 | 2020-02-13 | 基于反射率法的块约束时移地震差异反演方法及系统 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111239805A CN111239805A (zh) | 2020-06-05 |
| CN111239805B true CN111239805B (zh) | 2021-02-05 |
Family
ID=70873107
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010091778.3A Active CN111239805B (zh) | 2020-02-13 | 2020-02-13 | 基于反射率法的块约束时移地震差异反演方法及系统 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111239805B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114063148B (zh) * | 2020-07-31 | 2023-09-26 | 中国石油天然气股份有限公司 | 一种基于贝叶斯判别的折射波初至优选方法及系统 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102508294A (zh) * | 2011-10-20 | 2012-06-20 | 西北大学 | 一种利用时移地震勘探资料进行差异avo分析的方法 |
| CN106646603A (zh) * | 2017-01-04 | 2017-05-10 | 中海石油(中国)有限公司 | 一种实际时移地震资料处理差异的可靠性判断方法 |
| CN107092029A (zh) * | 2017-04-26 | 2017-08-25 | 中国石油大学(北京) | 一种地震反演方法和装置 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7636275B2 (en) * | 2007-02-06 | 2009-12-22 | Conocophillips Company | Direct time lapse inversion of seismic data |
-
2020
- 2020-02-13 CN CN202010091778.3A patent/CN111239805B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102508294A (zh) * | 2011-10-20 | 2012-06-20 | 西北大学 | 一种利用时移地震勘探资料进行差异avo分析的方法 |
| CN106646603A (zh) * | 2017-01-04 | 2017-05-10 | 中海石油(中国)有限公司 | 一种实际时移地震资料处理差异的可靠性判断方法 |
| CN107092029A (zh) * | 2017-04-26 | 2017-08-25 | 中国石油大学(北京) | 一种地震反演方法和装置 |
Non-Patent Citations (1)
| Title |
|---|
| 基于精确Zoeppritz方程的时移地震差异反演;周林;《2019年中国地球科学联合学术年会论文集》;20191027;第1142-1144页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111239805A (zh) | 2020-06-05 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US11231516B2 (en) | Direct migration of simultaneous-source survey data | |
| US7082368B2 (en) | Seismic event correlation and Vp-Vs estimation | |
| CN103293552B (zh) | 一种叠前地震资料的反演方法及系统 | |
| US9702993B2 (en) | Multi-parameter inversion through offset dependent elastic FWI | |
| US9702997B2 (en) | Wave-equation migration velocity analysis using image warping | |
| US20170176617A1 (en) | Automated near surface analysis by surface-consistent refraction methods | |
| US11294087B2 (en) | Directional Q compensation with sparsity constraints and preconditioning | |
| EP2728382B1 (en) | Spatial expansion seismic data processing method and apparatus | |
| US10310117B2 (en) | Efficient seismic attribute gather generation with data synthesis and expectation method | |
| US10838092B2 (en) | Estimating multiple subsurface parameters by cascaded inversion of wavefield components | |
| CN111522063B (zh) | 基于分频联合反演的叠前高分辨率流体因子反演方法 | |
| CN111025387B (zh) | 一种页岩储层的叠前地震多参数反演方法 | |
| Ma et al. | Multichannel block sparse Bayesian learning reflectivity inversion with lp-norm criterion-based Q estimation | |
| Djebbi et al. | Frequency domain multiparameter acoustic inversion for transversely isotropic media with a vertical axis of symmetry | |
| CN111025388B (zh) | 一种多波联合的叠前波形反演方法 | |
| CN110687597B (zh) | 一种基于联合字典的波阻抗反演方法 | |
| CN111239805B (zh) | 基于反射率法的块约束时移地震差异反演方法及系统 | |
| Liu et al. | Absolute acoustic impedance inversion using convolutional neural networks with transfer learning | |
| CN116088048A (zh) | 包含不确定性分析的各向异性介质多参数全波形反演方法 | |
| CN110988991B (zh) | 一种弹性参数反演方法、装置及系统 | |
| CN112363222A (zh) | 叠后自适应宽带约束波阻抗反演方法及装置 | |
| CN111077573A (zh) | 一种确定地层弹性参数的方法、装置及系统 | |
| Przebindowska | Acoustic full waveform inversion of marine reflection seismic data | |
| CN115453620B (zh) | 一种基于非稳态反演的avo校正方法 | |
| CN115480311A (zh) | 一种多道叠前波形反演方法及装置 |
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 |