CN111522063A - 基于分频联合反演的叠前高分辨率流体因子反演方法 - Google Patents

基于分频联合反演的叠前高分辨率流体因子反演方法 Download PDF

Info

Publication number
CN111522063A
CN111522063A CN202010362773.XA CN202010362773A CN111522063A CN 111522063 A CN111522063 A CN 111522063A CN 202010362773 A CN202010362773 A CN 202010362773A CN 111522063 A CN111522063 A CN 111522063A
Authority
CN
China
Prior art keywords
inversion
seismic
model
prestack
angle
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
CN202010362773.XA
Other languages
English (en)
Other versions
CN111522063B (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.)
Hunan University of Science and Technology
Original Assignee
Hunan University of Science and Technology
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 Hunan University of Science and Technology filed Critical Hunan University of Science and Technology
Priority to CN202010362773.XA priority Critical patent/CN111522063B/zh
Publication of CN111522063A publication Critical patent/CN111522063A/zh
Application granted granted Critical
Publication of CN111522063B publication Critical patent/CN111522063B/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. analysis, for interpretation, for correction
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/362—Effecting static or dynamic corrections; Stacking
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/38—Seismology; Seismic or acoustic prospecting or detecting specially adapted for water-covered areas
    • G01V1/3808—Seismic data acquisition, e.g. survey design
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/10—Aspects of acoustic signal generation or detection
    • G01V2210/12—Signal generation
    • G01V2210/129—Source location
    • G01V2210/1297—Sea bed
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/50—Corrections or adjustments related to wave propagation
    • G01V2210/51—Migration
    • G01V2210/512—Pre-stack

Abstract

本发明提供一种基于分频联合反演的叠前高分辨率流体因子反演方法，该方法包括：基于孔弹性理论对传统形式的精确Zoeppritz方程进行改写得到包含流体因子项的新形式精确Zoeppritz方程；利用工区实际地震数据的分频处理结果、测井数据统计得到的多个储层参数的先验分布函数及推导的新形式Zoeppritz方程构建贝叶斯框架下的多波数据分频联合反演目标函数；利用泰勒级数展开对目标函数进行求解得到所述储层参数扰动量的求解表达式；对求解表达式进行迭代求解，利用求解结果对初始储层参数值进行迭代更新，直至达到迭代终止条件。本发明能够提高储层流体因子的预测精度。

Description

基于分频联合反演的叠前高分辨率流体因子反演方法
技术领域
本发明涉及油气田地震勘探开发和储层参数预测方法，特别是一种基于分频联合反演的叠前高分辨率流体因子反演方法。
背景技术
储层流体识别是油气藏勘探开发的重要环节之一，其直接关系着致密砂岩储层和页岩储层等非常规储层的“甜点”预测。可靠的流体识别结果可以提高勘探成功率、降低投入成本，对油气藏勘探开发与储层评价具有重要意义。通常，利用岩芯和测井数据分析储层的含流体情况最为直接而且准确，但这必须以大量的岩芯样品和测井数据分析结果为基础，并且局部的岩芯和单口井的分析结果不能客观反映周围岩石的流体分布情况。当储层岩芯样品和测井数据有限时，基于岩芯和测井数据分析的储层流体分布结果横向分辨率低，难于实现整个工区储层流体分布评价。流体因子是储层流体识别的重要表征参数，在油气藏勘探开发中发挥着重要作用，而获取流体因子的有效途径之一就是叠前地震反演。地震数据具有很好的空间连续性，在储层描述与表征中起着非常重要的作用。高精度反演方法使反演结果具有重要价值，可以提高地震数据解释的可信度。利用反演结果可以更好地实现储层空间描述，很大程度地提高储层表征的价值。目前，流体因子主要由基于精确Zoeppritz方程线性近似公式的反演方法估计得到，但近似公式假设条件多，且在强阻抗差、大入射角等情况下，存在较大的计算误差，使得流体因子反演结果的精度不能满足实际储层精细表征要求。此外，大量学者研究指出，通过谱分解技术分离得到的地震数据中的高频信息可以对较薄的互层结构进行有效识别，但常规反演方法对这部分高频地震信息的利用并不充分。
总之，目前基于叠前地震数据的流体因子地震估计方法研究存在以下问题：1、现有包含流体因子项的精确Zoeppritz方程近似公式假设条件多，计算精度有限，极大地限制了现有叠前流体因子反演方法的适用性和精度。2、常规叠前AVA流体因子反演方法并不能直接反演得到流体因子项，而是通过反演得到其对应的反射率然后间接求解得到，这一过程降低了最终的估计精度。3、地震数据中包含的高频信息对薄储层具有较好的识别效果，但现有反演方法对这部分信息的利用不充分。4、储层流体通常具有明显的分层特征，如何更加清晰的刻画流体边界有待进一步研究。
发明内容
本发明的目的在于针对现有技术中存在的上述缺陷，提供一种基于分频联合反演的叠前高分辨率流体因子反演方法。本发明能提供一种高可靠性的储层流体地震识别方法，有效地提高现有叠前AVA流体因子反演方法的精度与稳定性，满足地震叠前反演识别油气储层、特别是复杂油气藏流体识别的要求。
本发明基于分频联合反演的叠前高分辨率流体因子反演方法，包括如下步骤：
(1)采用广义S正反变换对叠前角道集数据进行分频处理，得到不同频段对应的叠前角道集数据；
(2)基于不同频段对应的叠前角道集数据提取角度依赖的子波，基于精确Zoeppritz方程和测井数据正演合成地震角道集并结合实际井旁角道集数据确定振幅缩放因子；
(3)基于工区内所有测井数据统计模型参数的先验信息及其均值，包括流体因子、剪切模量和密度，并计算包含流体因子、剪切模量和密度统计相关性的协方差矩阵；
(4)利用地震构造解释资料和工区所有测井数据，基于沉积模式建立时间域的流体因子、剪切模量和密度参数的初始模型；
(5)基于孔弹性理论，对传统形式的精确Zoeppritz方程进行改写，得到包含流体因子、剪切模量、密度以及干岩石纵横波速度比的平方
Figure BDA0002475663090000031
的新形式FMR-Zoeppritz方程；其中，VP、VS是干岩石的纵、横波速度；
(6)首先由测井数据统计分析得到固定常数
Figure BDA0002475663090000032
的值，然后基于时间域的初始模型和新形式的FMR-Zoeppritz方程正演模拟角度域的多波地震叠前道集，最后由正演模拟记录和实际记录直接计算叠前反演残差；
(7)基于Bayesian原理、先验信息和新形式的FMR-Zoeppritz方程正演算子构建最大后验概率意义下的反演目标函数，借助泰勒展开对目标函数进行求解，得到流体因子、剪切模量和密度参数扰动量的求解表达式；
(8)基于上述参数扰动量的求解表达式和反演残差，采用迭代重加权最小二乘算法求解模型参数的扰动量，并利用求解得到的扰动量对初始模型进行更新；同时利用井旁道反演结果选择合理的权重因子，重复以上步骤(6)、(7)和(8)，并通过反演残差控制最大迭代次数，获得最优储层参数，包括流体因子、剪切模量和密度。
本发明方法具有如下有益效果：
(1)本发明能够提供一种高可靠性的储层流体地震识别方法，相比于基于岩芯和测井数据分析估算工区储层含流体情况的方法而言，能提供横向连续性更好且真实反映地下流体分布情况的流体因子数据体。
(2)本发明能够直接通过叠前地震数据反演得到流体因子，很好地消除了间接求解对流体因子估计精度的影响。
(3)本发明所述的叠前AVA反演基于新推导的FMR-Zoeppritz方程，相对于已有的线性近似方程而言，新的FMR-Zoeppritz方程在强阻抗差以及大入射角等条件下均具有非常高的计算精度，很好地克服了现有近似公式由于假设条件太多而导致的适用性差、反演精度低等问题。
(4)本发明引入了服从微分拉普拉斯分布的垂向块约束项，能够更好地刻画流体边界，提高流体信息预测结果的精度和分辨率。
(5)本发明能够提供一种高精度的分频联合反演方法，能够充分利用地震数据所包含的高频信息，进一步提升流体因子反演结果的纵向分辨率。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单的介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其它的附图。
图1是本发明基于分频联合反演的叠前高分辨率流体因子反演方法的流程框图。
图2是本发明实施例输入原始叠前AVA(Amplitude versus angle，振幅随入射角度的变化)PP波角道集地震数据。
图3至图5是本发明实施例以图1所示角道集数据为输入采用广义S正反变换进行分频处理得到的分频段叠前AVA(Amplitude versus angle，振幅随入射角度的变化)PP波角道集地震数据；其中，图3是低频，图4是中频，图5是高频。
图6是本发明实施例以图1所示原始叠前角道集作为输入地震数据并采用基于FMR-Zoeppritz方程的反演方法反演得到的流体因子(左)、剪切模量(中)、密度(右)。
图7是本发明实施例以图3所示不同频段对应的叠前角道集作为联合输入地震数据并采用基于FMR-Zoeppritz方程的分频联合反演方法反演得到的流体因子(左)、剪切模量(中)、密度(右)。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整的描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。
参见图1，是本发明基于分频联合反演的叠前高分辨率流体因子反演方法的流程框图，具体包括如下步骤：
101、采用广义S正反变换对叠前角道集数据进行分频处理，得到不同频段对应的叠前角道集数据：采用广义S变换将原始叠前角道集数据变换至频率域，结合测井数据分析地震数据的最大分辨能力，确定频段划分原则，根据这一原则对频率域的地震数据进行频段划分，最后将不同频段对应的地震数据采用广义S反变换转换至时间域，从而得到不同频段对应的时间域叠前角道集数据。
102、基于不同频段对应的叠前角道集数据提取角度依赖的子波，基于精确Zoeppritz方程和测井数据正演合成地震角道集并结合实际井旁角道集数据确定振幅缩放因子：首先，基于分解得到的不同频段对应的地震叠前角道集和测井数据采用统计方法提取若干个依赖于入射角度的地震子波；以测井数据为输入模型利用精确Zoeppritz方程正演模拟角度域的地震道集，与实际井旁角度域地震道集对比，计算振幅缩放因子，并应用于所提取的地震子波。
103、基于工区内所有测井数据统计模型参数的先验信息及其均值，包括流体因子、剪切模量和密度，并计算包含流体因子、剪切模量和密度统计相关性的协方差矩阵：假设模型参数服从高斯分布的同时引入服从微分拉普拉斯分布的垂向块约束项，通过测井数据分析获得需要的模型参数，求取各参数的自相关系数和互相关系数，构建协方差矩阵，形成符合该工区的模型参数先验分布函数。
104、利用地震构造解释资料(层位、断层等信息)和工区所有测井数据，基于沉积模式建立时间域的流体因子、剪切模量和密度参数的初始模型：利用地震构造解释资料(层位、断层等信息)，基于沉积模式建立地质模型，并将测井资料，按构造模式进行插值和外推，得到每条测线的初始参数模型，包括流体因子、剪切模量和密度的初始模型。
105、基于孔弹性理论，对传统形式的精确Zoeppritz方程进行改写，得到包含流体因子、剪切模量、密度以及干岩石纵横波速度比的平方
Figure BDA0002475663090000061
的新形式FMR-Zoeppritz方程：首先基于孔弹性理论，利用流体因子、剪切模量、纵横波速度以及密度之间的关系，对传统形式的Zoeppritz方程进行改写，最终得到关于流体因子、剪切模量、密度以及干岩石纵横波速度比的平方
Figure BDA0002475663090000062
的新形式FMR-Zoeppritz方程；其中，VP、VS是干岩石的纵、横波速度。
106、首先由测井数据统计分析得到固定常数
Figure BDA0002475663090000063
的值，然后基于时间域的初始模型和新形式的FMR-Zoeppritz方程正演模拟角度域的多波地震叠前道集，最后由正演模拟记录和实际记录直接计算叠前反演残差：首先利用测井数据进行统计分析，估计得到固定常数
Figure BDA0002475663090000064
的值，然后以时间域的初始模型为输入，直接利用FMR-Zoeppritz方程计算每一个角度的地震反射系数向量，将前述获取的子波与反射系数进行褶积，获得角度域的叠前角道集，并与实际角度域地震道集作差，获得叠前地震反演残差。
107、基于Bayesian原理、先验信息和新形式的FMR-Zoeppritz方程正演算子构建最大后验概率意义下的反演目标函数，借助泰勒展开对目标函数进行求解，得到流体因子、剪切模量和密度参数扰动量的求解表达式：利用以上所述先验模型对反演结果进行约束，并假设似然函数服从Gaussian(高斯)分布，先验模型服从Gaussian(高斯)与Differentiable Laplace(微分拉普拉斯)的联合分布；将似然函数和先验模型代入贝叶斯公式中，根据后验概率确定反演的目标函数，由于模型参数不好直接求解，因此考虑借助于泰勒级数展开对目标函数进行简化，变成关于模型参数扰动量的函数，然后将目标函数对扰动量进行求导并令导数等于零，得到扰动量的迭代求解公式；其中泰勒展开需要求取简化方程正演算子关于模型参数的一阶偏导数，该导数可以通过解析或数值方法求得。
108、基于上述参数扰动量的求解表达式和反演残差，采用迭代重加权最小二乘算法求解模型参数的扰动量，并利用求解得到的扰动量对初始模型进行更新，同时利用井旁道反演结果选择合理的权重因子，重复以上步骤106、107和108，并通过反演残差控制最大迭代次数，获得最优储层参数，包括流体因子、剪切模量和密度：以反演残差为输入，采用迭代重加权最小二乘算法对上述非线性的参数扰动量的求解表达式进行进一步求解，得到模型参数更新量。利用得到的更新量对初始模型进行更新，同时利用井旁道反演结果选择合理的权重因子，重复以上步骤106、107和108，并通过反演残差控制最大迭代次数确定最优的参数模型作为地震叠前AVA反演的最终结果，包括流体因子、剪切模量和密度。
简而言之，本发明方法的工作步骤如下：利用广义S正反变化对原始地震数据进行分频处理→基于分频地震数据提取角度依赖的子波，通过测井数据正演模拟和井旁道地震数据确定振幅缩放因子→基于工区内所有测井数据统计模型参数的先验信息→利用地震构造解释资料和测井数据，基于沉积模式建立储层参数的初始模型→基于孔弹性理论推导新形式的FMR-Zoeppritz方程，利用FMR-Zoeppritz方程正演模拟角度域的地震叠前道集并计算反演残差→借助泰勒级数展开将目标函数改写成关于模型参数扰动量的函数，然后对扰动量求导并令导数等于零得到扰动量的迭代求解公式→利用迭代重加权最小二乘算法对扰动量进行求解，并利用求解得到的扰动量对模型参数进行更新→重复以上工作步骤直至反演残差达到要求或达到最大迭代次数→输出最终模型参数，包括流体因子、剪切模量、密度模型。
下面结合具体实施例对本发明方法的工作步骤详细叙述如下：
(1)如图2所示，是本发明实施例输入的原始叠前AVA角道集地震数据，图中纵轴表示时间，单位为毫秒，横轴表示入射角度。利用广义S变换将原始叠前角道集数据变换至频率域，结合测井数据分析地震数据的最大分辨能力，确定频段划分原则，根据这一原则对频率域的地震数据进行频段划分，最后将不同频段对应的地震数据利用广义S反变换转换至时间域，从而得到不同频段对应的时间域叠前角道集数据。如图3、图4、图5所示，是本发明实施例获得的分频叠前AVA角道集地震数据，图3为低频段对应的角道集数据，图4为中频段对应的角道集数据，图5为高频段对应的角道集数据，图中纵轴表示时间，单位为秒，横轴表示入射角度。
(2)本发明假设反演之前地震子波已知，因而，需要基于实际的地震叠前道集和测井数据采取统计方法提取子波，子波在传播过程中受地层的影响会发生波形或频率变化，提取依赖于入射角度的地震子波能有效提高振幅匹配程度。
实际的地震振幅往往是相对值，采用精确Zoeppritz方程正演模拟的地震数据振幅与实际振幅存在一定数值差异。以测井数据为输入模型利用精确Zoeppritz方程正演模拟角度域的叠前地震道集，与实际井旁角度域地震道集对比，计算振幅缩放因子，并应用于所提取的地震子波，达到模拟记录与实际记录的振幅匹配。当地震数据信噪比较高时，为角道集的每一道使用统一的振幅缩放因子，以保证振幅随偏移距的变化关系；当信噪比低时，可分近、中、远偏移距分别计算振幅缩放因子，保证模拟记录与实际记录的最佳匹配，减少噪声对反演过程的影响。
(3)建立模型参数的先验模型。为了降低地震反演的不确定性，提高反演过程的稳定性，需要从其它途径获得地下介质地震弹性参数模型的信息作为先验信息。本发明采用高斯分布函数和微分拉普拉斯分布函数作为先验分布函数，基于工区内所有测井数据得到的各模型参数及其均值，求取各参数的自相关系数和互相关系数，构建协方差矩阵，形成符合该工区的模型参数先验分布函数，其表达式为：
Figure BDA0002475663090000091
其中，Cm为包含模型参数统计相关性的块对角矩阵；m为模型参数向量；μ为模型参数的均值向量(不同的模型参数需要分别求取)；D是一阶微分算子；分母中的ks,s＝1,2,3是尺度参数，对于不同的模型参数其值也会不一样。
(4)利用地震构造解释资料(层位、断层等信息)，基于沉积模式建立地质模型，并将测井资料，按构造模式进行插值和外推，得到每条测线的参数初始模型，包括流体因子初始模型、剪切模量初始模型和密度初始模型。
建立弹性参数模型主要利用三维空间插值方法，其技术流程是首先利用散点插值的方法对各个层位的数据进行插值，完成地质层位建模，然后根据地质层位进行弹性参数横向插值，即将测井信息进行横向插值，计算得到地下每个点上的弹性参数值，完成初始弹性参数建模的任务。
(5)根据孔弹性理论，对传统形式的Zoeppritz方程进行改写，得到关于流体因子、剪切模量、密度以及干岩石纵横波速度比的平方
Figure BDA0002475663090000101
的新形式FMR-Zoeppritz方程，其推导过程如下：
Russell流体因子f表达式如下：
Figure BDA0002475663090000102
其中，VP、VS和ρ是饱和岩石的纵、横波速度和密度；
Figure BDA0002475663090000103
表示干岩石纵横波速度比的平方。
Zoeppritz方程描述的平面波入射到两个半无限均匀各向同性弹性介质的分界面时的反射和透射情况。其表达式如下：
Figure BDA0002475663090000104
根据Snell定律有：
Figure BDA0002475663090000105
将等式(4)代入等式(3)中，将原式中的反射角
Figure BDA0002475663090000106
和透射角θ2、
Figure BDA0002475663090000107
进行替换，得到如下表达式：
Figure BDA0002475663090000111
在各向同性介质中，纵波速度和横波速度与弹性模量存在如下关系：
Figure BDA0002475663090000112
其中，M和μ分别表示饱和岩石的纵波模量和横波模量。
将公式(6)代入公式(2)可得：
Figure BDA0002475663090000113
然后将公式(6)和公式(7)代入公式(5)可得新形式的FMR-Zoeppritz方程：
Figure BDA0002475663090000114
基于新形式的FMR-Zoeppritz方程的正演过程可表示如下：
Figure BDA0002475663090000115
式(9)中，g为基于FMR-Zoeppritz方程的正演算子；W表示子波矩阵；r为反射系数序列；下标pp和ps代表不同的反射波模式。
(6)首先由测井数据统计分析得到固定常数
Figure BDA0002475663090000116
的值，然后基于时间域的初始模型和新形式的FMR-Zoeppritz方程正演模拟角度域的多波地震叠前道集，最后由正演模拟记录和实际记录直接计算叠前反演残差：首先利用测井数据进行统计分析，估计得到固定常数
Figure BDA0002475663090000121
的值，然后以时间域的初始模型为输入，直接利用FMR-Zoeppritz方程计算每一个角度的地震反射系数向量，将前述获取的子波与反射系数进行褶积，获得角度域的叠前角道集，并与实际角度域地震道集作差(如图3、图4、图5所示，是本发明实施例输入的实际分频叠前AVA角道集地震数据，图3为低频段对应的角道集数据，图4为中频段对应的角道集数据，图5为高频段对应的角道集数据，图中纵轴表示时间，单位为秒，横轴表示入射角度)，获得叠前地震反演残差。
(7)基于Bayesian原理、先验信息和新形式的FMR-Zoeppritz方程正演算子构建最大后验概率意义下的反演目标函数，借助泰勒展开对目标函数进行求解，得到流体因子、剪切模量和密度参数扰动量的求解表达式：利用以上所述先验模型对反演结果进行约束，并假设似然函数服从Gaussian(高斯)分布，先验模型服从Gaussian(高斯)与DifferentiableLaplace(微分拉普拉斯)的联合分布。将似然函数和先验模型代入贝叶斯公式中，根据后验概率确定反演的目标函数，由于模型参数不好直接求解，因此考虑借助于泰勒级数展开对目标函数进行简化，变成关于模型参数扰动量的函数，然后将目标函数对扰动量进行求导并令导数等于零，得到扰动量的迭代求解公式。其中泰勒展开需要求取FMR-Zoeppritz方程正演算子关于流体因子、剪切模量和密度参数的一阶偏导数，该导数可以通过解析或数值方法求得。设模型参数为mT＝(m1,m2,",mn)T，观测地震数据为dT＝(d1,d2,",dn)T，由贝叶斯理论可以得到，在已知叠前地震数据的情况下，反演地下介质弹性参数的问题可以归结为求解一个后验概率函数：
Figure BDA0002475663090000122
其中，P(d)＝∫p(d|m)p(m)dm为归一化因子，可以看作是常数；P(d|m)为似然函数；P(m)为先验概率分布。假设似然函数满足高斯分布，则有：
Figure BDA0002475663090000131
式(11)中，g为正演算子，CD是噪声的协方差矩阵。
对于分频联合反演而言，利用广义S正反变换对某一个角度对应的地震数据进行分解得到不同频段的地震数据在尺度上相互独立，假设不同频段地震数据包含的噪声数据均服从高斯分布，则分频数据多波联合反演的似然函数可表示为：
Figure BDA0002475663090000132
出于实际问题考虑，一般假设噪声相对于观测数据而言是不相关的，那么噪声的协方差矩阵就会简化为一个对角矩阵，在这种情况下，噪声的协方差矩阵可表示为，CD＝σn 2I，其中，σn 2是噪声的方差，I是Nd×Nd的单位阵，Nd是观测数据的长度。则多波地震数据后验概率分布可以表示为：
Figure BDA0002475663090000141
通常，上述后验概率分布函数最大值对应的解即为我们所求的最优解。从式(13)可以看出，求解最大后验概率解的问题转化为求解如下目标函数的最小值的问题：
Figure BDA0002475663090000142
其中，ηi,i＝2,...l控制PP波数据不同频率成分的权重，αi,i＝1,2,...l控制PS转换波不同频率成分的权重，β＝σpp控制先验信息的权重。
上述目标函数是复杂的非线性函数，可借助泰勒级数展开对其进行求解。首先，将基于新形式FMR-Zoeppritz方程的正演算子g(m)在初始模型参数m0处进行泰勒展开，如下：
Figure BDA0002475663090000151
对上式取一阶近似并带入公式(14)，同时对约束项进行改写可得：
Figure BDA0002475663090000152
将上面得到的目标函数对弹性参数扰动量Δm进行求导并令导数等于零得：
Figure BDA0002475663090000153
其中：
Figure BDA0002475663090000154
Figure BDA0002475663090000155
Figure BDA0002475663090000156
由于等式两边都含有Δm，因此采用迭代重加权最小二乘算法进行求解。这样我们就可以得到模型参数的更新量Δm0，利用这一更新量可以对初始模型m0进行更新：
m1＝m0+Δm0 (18)；
利用m1取代公式(17)中的m0，我们可以计算得到新的模型更新量Δm1。这样模型参数就形成了一个持续迭代更新的过程，迭代更新公式如下：
mj+1＝mj+λjΔmj,j＝0,1,2... (19)；
其中，λj是j次迭代的步长，
Figure BDA0002475663090000161
其中正演算子关于模型参数的一阶导数
Figure BDA0002475663090000162
和
Figure BDA0002475663090000163
可通过解析的方
法求得：
Figure BDA0002475663090000164
其中，W是子波矩阵，rPP(m)和rPS(m)分别为纵波反射系数序列和转换横波反射系数序列。对于其中某一个角度θi有：
Figure BDA0002475663090000171
将等式(21)推广到n个输入角度的情况，G(m)关于模型参数的一阶偏导数就会变成一个(n×N)×3N矩阵，然后再推广到纵横波联合反演的情况，
Figure BDA0002475663090000172
最终就变成了一个(2n×N)×3N的矩阵。
将公式(8)所示的FMR-Zoeppritz方程的写成矩阵形式如下：
AR＝b (22)；
将公式(22)两边同时对界面两侧的流体因子、剪切模量和密度求偏导：
Figure BDA0002475663090000173
整理便可得到某一界面反射系数和透射系数关于界面两侧的流体因子、剪切模量和密度的偏导数：
Figure BDA0002475663090000174
将上式计算得到结果代入公式(21)中，便可得到矩阵
Figure BDA0002475663090000175
的值，进而求得模型参数的更新量，实现模型参数的迭代更新。
(8)以反演残差为输入，采用迭代重加权最小二乘算法对上述非线性的参数扰动量的求解表达式进行进一步求解，得到模型参数更新量，利用得到的更新量对初始模型进行更新。
(9)同时利用井旁道反演结果选择合理的权重因子，重复以上步骤(6)、(7)和(8)，并通过反演残差控制最大迭代次数确定最优的参数模型作为地震叠前AVA反演的最终结果，包括流体因子、剪切模量和密度。如图7所示，是本发明实施例利用基于FMR-Zoeppritz方程的分频联合反演方法反演得到的流体因子(左)、剪切模量(中)和密度(右)模型图，图中纵轴表示时间，单位秒，横轴自左至右表示流体因子(单位：Gpa)、剪切模量(单位：Gpa)和密度(单位：克/立方厘米))。图6是本发明实施例以图1所示原始叠前角道集为输入利用基于FMR-Zoeppritz方程的反演方法反演得到流体因子(左)、剪切模量(中)和密度(右)模型图。通过对比分析可以看出，基于新FMR-Zoeppritz方程的反演算法能够高精度的估计流体因子，并且分频联合反演能够进一步提升流体因子的估计精度和分辨率。
本发明由于采取以上技术方案，具有以下优点：1、本发明技术方案能够提供一种高可靠性的储层流体地震识别方法，相比于基于岩芯和测井数据分析估算工区储层含流体情况的方法而言，能提供横向连续性更好且真实反映地下流体分布情况的流体因子数据体。2、本发明技术方案可以直接通过叠前地震数据反演得到能够流体因子，很好的消除了间接求解对流体因子估计精度的影响。3、本发明技术方案所述的叠前AVA反演基于新推导的FMR-Zoeppritz方程，相对于已有的线性近似方程而言，新的FMR-Zoeppritz方程在强阻抗差以及大入射角等条件下均具有非常高的计算精度，很好的克服了现有近似公式由于假设条件太多而导致的适用性差、反演精度低等问题。4、本发明技术方案引入了服从微分拉普拉斯分布的垂向块约束项，能够更好的刻画流体边界，提高流体信息预测结果的精度和分辨率。5、本发明技术方案能够提供一种高精度的分频联合反演方法，能够充分利用地震数据所包含的高频信息，进一步提升流体因子反演结果的纵向分辨率。

Claims (9)

1.一种基于分频联合反演的叠前高分辨率流体因子反演方法，其特征在于包括如下步骤：
(1)采用广义S正反变换对叠前角道集数据进行分频处理，得到不同频段对应的叠前角道集数据；
(2)基于不同频段对应的叠前角道集数据提取角度依赖的子波，基于精确Zoeppritz方程和测井数据正演合成地震角道集并结合实际井旁角道集数据确定振幅缩放因子；
(3)基于工区内所有测井数据统计模型参数的先验信息及其均值，包括流体因子、剪切模量和密度，并计算包含流体因子、剪切模量和密度统计相关性的协方差矩阵；
(4)利用地震构造解释资料和工区所有测井数据，基于沉积模式建立时间域的流体因子、剪切模量和密度参数的初始模型；
(5)基于孔弹性理论，对传统形式的精确Zoeppritz方程进行改写，得到包含流体因子、剪切模量、密度以及干岩石纵横波速度比的平方
Figure FDA0002475663080000011
的新形式FMR-Zoeppritz方程；其中，VP、VS是干岩石的纵、横波速度；
(6)首先由测井数据统计分析得到固定常数
Figure FDA0002475663080000012
的值，然后基于时间域的初始模型和新形式的FMR-Zoeppritz方程正演模拟角度域的多波地震叠前道集，最后由正演模拟记录和实际记录直接计算叠前反演残差；
(7)基于Bayesian原理、先验信息和新形式的FMR-Zoeppritz方程正演算子构建最大后验概率意义下的反演目标函数，借助泰勒展开对目标函数进行求解，得到流体因子、剪切模量和密度参数扰动量的求解表达式；
(8)基于上述参数扰动量的求解表达式和反演残差，采用迭代重加权最小二乘算法求解模型参数的扰动量，并利用求解得到的扰动量对初始模型进行更新；同时利用井旁道反演结果选择合理的权重因子，重复以上步骤(6)、(7)和(8)，并通过反演残差控制最大迭代次数，获得最优储层参数，包括流体因子、剪切模量和密度。
2.根据权利要求1所述基于分频联合反演的叠前高分辨率流体因子反演方法，其特征在于：步骤(1)具体是，采用广义S变换将原始叠前角道集数据变换至频率域，结合测井数据分析地震数据的最大分辨能力，确定频段划分原则，根据这一原则对频率域的地震数据进行频段划分，最后将不同频段对应的地震数据采用广义S反变换转换至时间域，从而得到不同频段对应的时间域叠前角道集数据。
3.根据权利要求1所述基于分频联合反演的叠前高分辨率流体因子反演方法，其特征在于：步骤(2)具体是，首先，基于分解得到的不同频段对应的地震叠前角道集和测井数据采用统计方法提取若干个依赖于入射角度的地震子波；以测井数据为输入模型利用精确Zoeppritz方程正演模拟角度域的地震道集，与实际井旁角度域地震道集对比，计算振幅缩放因子，并应用于所提取的地震子波。
4.根据权利要求1所述基于分频联合反演的叠前高分辨率流体因子反演方法，其特征在于：步骤(3)具体是，假设模型参数服从高斯分布的同时引入服从微分拉普拉斯分布的垂向块约束项，通过测井数据分析获得需要的模型参数，求取各参数的自相关系数和互相关系数，构建协方差矩阵，形成符合该工区的模型参数先验分布函数。
5.根据权利要求1所述基于分频联合反演的叠前高分辨率流体因子反演方法，其特征在于：步骤(4)具体是，利用地震构造解释资料，基于沉积模式建立地质模型，并将测井资料，按构造模式进行插值和外推，得到每条测线的初始参数模型，包括流体因子、剪切模量和密度的初始模型。
6.根据权利要求1所述基于分频联合反演的叠前高分辨率流体因子反演方法，其特征在于：步骤(5)具体是，首先基于孔弹性理论，利用流体因子、剪切模量、纵横波速度以及密度之间的关系，对传统形式的Zoeppritz方程进行改写，最终得到关于流体因子、剪切模量、密度以及干岩石纵横波速度比的平方
Figure FDA0002475663080000031
的新形式FMR-Zoeppritz方程。
7.根据权利要求1所述基于分频联合反演的叠前高分辨率流体因子反演方法，其特征在于：步骤(6)具体是，首先利用测井数据进行统计分析，估计得到固定常数
Figure FDA0002475663080000032
的值，然后以时间域的初始模型为输入，直接利用FMR-Zoeppritz方程计算每一个角度的地震反射系数向量，将前述获取的子波与反射系数进行褶积，获得角度域的叠前角道集，并与实际角度域地震道集作差，获得叠前地震反演残差。
8.根据权利要求1所述基于分频联合反演的叠前高分辨率流体因子反演方法，其特征在于：步骤(7)具体是，利用以上所述先验模型对反演结果进行约束，并假设似然函数服从Gaussian即高斯分布，先验模型服从Gaussian即高斯与Differentiable Laplace即微分拉普拉斯的联合分布；将似然函数和先验模型代入贝叶斯公式中，根据后验概率确定反演的目标函数，由于模型参数不好直接求解，因此考虑借助于泰勒级数展开对目标函数进行简化，变成关于模型参数扰动量的函数，然后将目标函数对扰动量进行求导并令导数等于零，得到扰动量的迭代求解公式；其中泰勒展开需要求取简化方程正演算子关于模型参数的一阶偏导数，该导数可以通过解析或数值方法求得。
9.根据权利要求1所述基于分频联合反演的叠前高分辨率流体因子反演方法，其特征在于：步骤(8)具体是，以反演残差为输入，采用迭代重加权最小二乘算法对非线性的参数扰动量的求解表达式进行进一步求解，得到模型参数更新量，利用得到的更新量对初始模型进行更新；同时利用井旁道反演结果选择合理的权重因子，重复步骤(6)、步骤(7)和步骤(8)，并通过反演残差控制最大迭代次数确定最优的参数模型作为地震叠前AVA反演的最终结果，包括流体因子、剪切模量和密度。
CN202010362773.XA 2020-04-30 2020-04-30 基于分频联合反演的叠前高分辨率流体因子反演方法 Active CN111522063B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010362773.XA CN111522063B (zh) 2020-04-30 2020-04-30 基于分频联合反演的叠前高分辨率流体因子反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010362773.XA CN111522063B (zh) 2020-04-30 2020-04-30 基于分频联合反演的叠前高分辨率流体因子反演方法

Publications (2)

Publication Number Publication Date
CN111522063A true CN111522063A (zh) 2020-08-11
CN111522063B CN111522063B (zh) 2022-04-12

Family

ID=71908449

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010362773.XA Active CN111522063B (zh) 2020-04-30 2020-04-30 基于分频联合反演的叠前高分辨率流体因子反演方法

Country Status (1)

Country Link
CN (1) CN111522063B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113156510A (zh) * 2021-04-27 2021-07-23 中南大学 一种页岩储层脆性和各向异性参数预测方法及系统
CN113156509A (zh) * 2021-04-25 2021-07-23 中南大学 基于饱和介质精确Zeoppritz方程的地震振幅反演方法及系统
CN113253347A (zh) * 2021-05-14 2021-08-13 中南大学 基于vti介质的页岩储层avo反演表征方法及系统
CN114428303A (zh) * 2020-09-30 2022-05-03 中国石油化工股份有限公司 一种基于高精度非线性反演算法的高分辨分频联合反演方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20040199330A1 (en) * 2003-04-01 2004-10-07 Conocophillips Company Simultaneous inversion for source wavelet and AVO parameters from prestack seismic data
CN103257361A (zh) * 2013-05-24 2013-08-21 中国石油天然气集团公司 基于Zoeppritz方程近似式的油气预测方法及系统
CN104597490A (zh) * 2015-01-28 2015-05-06 中国石油大学(北京) 基于精确Zoeppritz方程的多波AVO储层弹性参数反演方法
WO2018201114A1 (en) * 2017-04-28 2018-11-01 Pioneer Natural Resources Usa, Inc. High resolution seismic data derived from pre-stack inversion and machine learning
CN109765613A (zh) * 2019-03-04 2019-05-17 太原理工大学 基于最速梯降叠前精确方程流体反演的页岩气识别方法
US20200041692A1 (en) * 2018-07-31 2020-02-06 Jan Schmedes Detecting Fluid Types Using Petrophysical Inversion

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20040199330A1 (en) * 2003-04-01 2004-10-07 Conocophillips Company Simultaneous inversion for source wavelet and AVO parameters from prestack seismic data
CN103257361A (zh) * 2013-05-24 2013-08-21 中国石油天然气集团公司 基于Zoeppritz方程近似式的油气预测方法及系统
CN104597490A (zh) * 2015-01-28 2015-05-06 中国石油大学(北京) 基于精确Zoeppritz方程的多波AVO储层弹性参数反演方法
WO2018201114A1 (en) * 2017-04-28 2018-11-01 Pioneer Natural Resources Usa, Inc. High resolution seismic data derived from pre-stack inversion and machine learning
US20200041692A1 (en) * 2018-07-31 2020-02-06 Jan Schmedes Detecting Fluid Types Using Petrophysical Inversion
CN109765613A (zh) * 2019-03-04 2019-05-17 太原理工大学 基于最速梯降叠前精确方程流体反演的页岩气识别方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
LIN ZHOU 等: "Prestack amplitude versus angle inversion for Young"s modulus and Poisson"s ratio based on the exact Zoeppritz equations", 《GEOPHYSICAL PROSPECTING》 *
ZHAOYUN ZONG 等: "Direct inversion for a fluid factor and its application in heterogeneous", 《GEOPHYSICAL PROSPECTING》 *
冉然 等: "基于Zoeppritz方程的纵横波模量反演", 《物探与化探》 *
周林 等: "基于精确Ｚ ｏｅ ｐｐ ｒｉｔｚ方程的非线性ＡＶＯ三参数反演", 《地球物理学报》 *

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114428303A (zh) * 2020-09-30 2022-05-03 中国石油化工股份有限公司 一种基于高精度非线性反演算法的高分辨分频联合反演方法
CN113156509A (zh) * 2021-04-25 2021-07-23 中南大学 基于饱和介质精确Zeoppritz方程的地震振幅反演方法及系统
CN113156509B (zh) * 2021-04-25 2022-07-08 中南大学 基于饱和介质精确Zoeppritz方程的地震振幅反演方法及系统
CN113156510A (zh) * 2021-04-27 2021-07-23 中南大学 一种页岩储层脆性和各向异性参数预测方法及系统
CN113156510B (zh) * 2021-04-27 2022-07-01 中南大学 一种页岩储层脆性和各向异性参数预测方法及系统
CN113253347A (zh) * 2021-05-14 2021-08-13 中南大学 基于vti介质的页岩储层avo反演表征方法及系统
CN113253347B (zh) * 2021-05-14 2022-07-01 中南大学 基于vti介质的页岩储层avo反演表征方法及系统

Also Published As

Publication number Publication date
CN111522063B (zh) 2022-04-12

Similar Documents

Publication Publication Date Title
CN111522063B (zh) 基于分频联合反演的叠前高分辨率流体因子反演方法
CN111208561B (zh) 基于时变子波与曲波变换约束的地震声波阻抗反演方法
EP1746443B1 (en) Method of estimating elastic parameters and rock composition of underground formations using seismic data
US7082368B2 (en) Seismic event correlation and Vp-Vs estimation
CA3111394C (en) Reservoir characterization utilizing inversion of resampled seismic data
CA2293627A1 (en) Method for characterizing subsurface petrophysical properties using linear shape attributes
CN113031068B (zh) 一种基于反射系数精确式的基追踪叠前地震反演方法
CN114994758B (zh) 碳酸盐岩断控储层的波阻抗提取与结构表征方法和系统
CN113156510B (zh) 一种页岩储层脆性和各向异性参数预测方法及系统
CN113156500A (zh) 数据驱动的快速构造约束叠前地震多道反演方法
CN106842291B (zh) 一种基于叠前地震射线阻抗反演的不整合圈闭储层岩性预测方法
CN111077578B (zh) 岩层分布预测方法和装置
Eikrem et al. Bayesian estimation of reservoir properties—effects of uncertainty quantification of 4D seismic data
US20230140656A1 (en) Method and system for determining seismic processing parameters using machine learning
CN113589365B (zh) 基于时频域信息的储层尖灭线描述方法
CN112147677B (zh) 油气储层参数标签数据生成方法及装置
CN111239805B (zh) 基于反射率法的块约束时移地震差异反演方法及系统
CN111308550B (zh) 一种页岩vti储层的各向异性参数多波联合反演方法
CN111897004B (zh) 一种基于大数据分析技术的测井预测方法
RU2764378C1 (ru) Способ повышения разрешающей способности данных сейсморазведки и прогнозирования геологического строения в межскважинном пространстве на основе метода спектральной инверсии
Ng et al. Estimation of facies probabilities on the Snorre field using geostatistical AVA inversion
CN114428303A (zh) 一种基于高精度非线性反演算法的高分辨分频联合反演方法
CN115480311A (zh) 一种多道叠前波形反演方法及装置
Soares Conditioning reservoir models on vast data sets
CN118033733A (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