CN114428303A - 一种基于高精度非线性反演算法的高分辨分频联合反演方法 - Google Patents

一种基于高精度非线性反演算法的高分辨分频联合反演方法 Download PDF

Info

Publication number
CN114428303A
CN114428303A CN202011060018.2A CN202011060018A CN114428303A CN 114428303 A CN114428303 A CN 114428303A CN 202011060018 A CN202011060018 A CN 202011060018A CN 114428303 A CN114428303 A CN 114428303A
Authority
CN
China
Prior art keywords
inversion
work area
seismic
data
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.)
Pending
Application number
CN202011060018.2A
Other languages
English (en)
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 CN202011060018.2A priority Critical patent/CN114428303A/zh
Publication of CN114428303A publication Critical patent/CN114428303A/zh
Pending legal-status Critical Current

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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • 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/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
    • 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
    • G01V2210/6242—Elastic parameters, e.g. Young, Lamé or Poisson
    • 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/63—Seismic attributes, e.g. amplitude, polarity, instant phase
    • 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/66—Subsurface modeling

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

一种基于高精度非线性反演算法的高分辨分频联合反演方法
技术领域
本发明涉及油气田、页岩气地震勘探和储层参数预测方法，特别是一种基于高精度非线性反演算法的高分辨分频联合反演方法。
背景技术
页岩气作为全球重要接替资源，其勘探开发技术成为目前人们的研究热点和难点。随着页岩储层勘探的不断深入，其对叠前AVA反演结果的分辨率要求也是越来越高，从储层结构上看，以页岩为主的非常规储层大多呈现互层结构，比如盐间页岩油储层，这类互层沉积特征使得有限带宽的地震数据不能很好的对其进行分辨。针对，这一问题，相关学者开展大量研究，杜炳毅等基于Rüger及崔杰的近似反射系数公式，选取多角度纵波和转换横波叠前道集，以岩石物理模型和测井资料作为约束实现了VTI介质纵横波联合反演。侯栋甲等基于Rüger近似公式实现了基于贝叶斯理论的VTI介质多波叠前联合反演。Zhang等基于改进的Rüger近似公式，实现了AVO梯度和截距以及各向异性纵波速度的同时反演，并提出了VTI介质的两步各向异性AVO反演方法，先反演出各向异性纵波速度，然后再利用远偏移距数据估计出Thomsen参数。雍杨等推导了VTI介质的多波AVA方程，并实现了基于该方程的多波联合弹性参数和各向异性参数反演。Zhang和Li推导得到了包含各向异性参数二阶项的VTI介质反射系数近似表达式并实现了同时反演，并测试了VTI介质弹性参数和各向异性参数同时反演方法在富含粘土的页岩地层的应用效果。但这些近似公式推导过程中的假设条件使得它们在强阻抗、强各向异性以及大角度等条件下，存在较大的误差，极大地限制了AVA反演的精度。钟峙等利用差分法实现了基于精确VTI反射系数方程的叠前AVA岩性参数反演。Amit和Subhashis利用遗传算法实现了基于VTI精确反射系数方程多波叠前联合反演。但这些方法效率低，目前难以应用于实际生产。大量学者研究指出，通过谱分解技术分离得到的地震数据中的高频信息可以对这类较薄的互层结构进行有效识别，但常规反演方法对这部分高频地震信息的利用并不充分。现有生产中的分频联合反演方法只是简单的将不同频段的地震数据与常规线性反演相结合，由于受到线性反演方法的精度和分辨率的影响，使得最终结合的效果远不能满足现阶段页岩储层勘探开发的高分辨率需求。此外，现阶段虽然有学者实现了基于近似公式的分频联合反演，如张繁昌等基于Aki-Richards近似方程实现了叠前AVA多频信息同时反演，但没有对具体的频段划分原则，以及高频地震信息的权重因子调节方法进行系统的研究和阐述。不合理的频段划分和权重因子选择会是反演结果的失真，从而导致最终的反演结果失去参考价值。
发明内容
针对上述问题，本发明实施例提供一种基于高精度非线性反演算法的高分辨分频联合反演方法，以提高现有叠前AVA反演方法的分辨率，满足叠前地震反演识别油气储层，特别是具有互层结构特征的盐间页岩储层表征的要求。所述方法包括以下步骤：
S101，利用工区测井数据正演分析确定工区地震数据的频段划分原则，并基于该原则对工区地震数据进行分频处理，以获得不同频段的实际地震角道集数据；
S102，利用不同频段的实际地震角道集数据提取角度依赖与频率依赖的地震子波；
S103，基于工区测井数据确定工区模型参数及其均值，并工区模型参数及其均值基于计算包含地层参数统计相关性的协方差矩阵，以形成该工区的模型参数先验分布函数；
S104，利用工区地震构造解释资料和测井数据，建立时间域的初始地层弹性参数模型；
S105，基于时间域的初始地层弹性参数模型和角度依赖与频率依赖的地震子波正演模拟不同频段的多波地震叠前角道集，以获得不同频段的模拟地震角道集数据，通过比较不同频段的模拟地震角道集数据与不同频段的实际地震角道集数据来确定叠前反演残差；
S106，根据贝叶斯理论，利用高精度正演方程和不同频段的实际地震角道集数据，通过引入该工区的模型参数先验分布函数构建最大后验概率意义下的反演目标函数，并对反演目标函数进行求解，以获得模型参数扰动量的求解表达式；
S107，基于模型参数扰动量的求解表达式和叠前反演残差，采用迭代重加权最小二乘算法计算模型参数的扰动量，同时利用工区井旁道反演结果选择相应的权重因子，重复迭代以上步骤S104至S106，并通过叠前反演残差的大小来控制最大迭代次数，从而获得高分辨率地层参数反演结果。
根据本发明的一个实施例，上述步骤S101包括以下内容：
S101.1，利用工区测井数据正演分析确定工区地震数据的频段划分原则；
S101.2，利用广义S变换将工区地震数据从时间域变换到频率域；
S101.3，按照确定好的频段划分原则对变换到频率域的工区地震数据进行频段划分，并将分离开来的不同频段的地震数据从频率域反变换回时间域。
根据本发明的一个实施例，上述步骤S102包括以下内容：
利用统计方法，针对不同频段的地震数据，从每一个频段的地震数据提取若干依赖于角度的地震子波，将各个角度和频率的子波与振幅缩放因子相乘，以得到一组角度依赖与频率依赖的地震子波集；
其中，以工区测井数据为输入，利用精确反射系数方程正演模拟角度道集，将模拟的角度道集与工区实际井旁角度域地震道集进行对比，确定所述振幅缩放因子。
根据本发明的一个实施例，上述步骤S103包括以下内容：
S103.1，假设工区模型参数服从柯西分布，基于工区测井数据确定工区模型参数及其均值，并计算各模型参数的自相关系数和互相关系数；
S103.2，基于各模型参数的自相关系数和互相关系数构建包含地层参数统计相关性的协方差矩阵，以形成该工区的模型参数先验分布函数。
根据本发明的一个实施例，上述步骤S104包括以下内容：
S104.1，利用工区地震构造解释资料，根据沉积模式建立工区地质模型；
S104.2，将工区测井资料，按照工区地质模型的构造模式进行插值和外推，得到地下每个点上的弹性参数值，从而获得该工区的时间域的初始地层弹性参数模型。
根据本发明的一个实施例，上述步骤S105包括以下内容：
S105.1，以时间域的初始地层弹性参数模型为输入，利用高计算精度的正演方程计算每一个对应角度的反射系数向量，并将频率依赖的地震子波与反射系数向量进行褶积，从而获得不同频段的模拟地震角道集数据；
S105.2，将不同频段的模拟地震角道集数据与不同频段的实际地震角道集数据作差，获得叠前反演残差。
根据本发明的一个实施例，上述步骤S106中，所述反演目标函数为AVA非线性反演目标函数。
根据本发明的一个实施例，上述步骤S106中，采用广义线性反演方法对AVA非线性反演目标函数进行求解。
此外，本发明还提供一种存储介质，其中存储有计算机程序，其特征在于，所述计算机程序被处理器执行时，实现如上所述的基于高精度非线性反演算法的高分辨分频联合反演方法的步骤。
此外，本发明还提供一种计算机设备，包括存储器和处理器，其中所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时，实现如上所述的基于高精度非线性反演算法的高分辨分频联合反演方法的步骤。
本发明利用以上技术方案，其具有以下优点：
1、本发明的技术方案能够提供一种高分辨率的储层参数估计方法，相比于常规反演方法而言能够提供可靠且分辨率更高的地层参数反演结果；
2、本发明的技术方案能够提供一种用于分频联合反演的频段划分原则，使最终的分频联合反演结果更加合理；
3、本发明的技术方案所述的基于高精度正演算子的分频联合反演方法很好的融合了高精度分线性反演方法与分频联合反演方法的优势，能够进一步有效的提升反演结果的分辨率。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其它的附图。
图1是本发明实施例基于高精度非线性反演算法的高分辨分频联合反演方法的流程图；
图2是本发明实施例输入的叠前AVA(Amplitude versus angle，振幅随入射角度的变化)PP波角道集地震数据；
图3是本发明实施例输入的分频叠前AVA(Amplitude versus angle，振幅随入射角度的变化)PP波角道集地震数据，其中，0-10Hz(上)、10-30Hz(中)、30-40Hz(下)；
图4是本发明实施例以图2所示角道集作为输入地震数据并采用基于精确Zoeppritz方程反演方法反演得到的纵波速度(a)、横波速度(b)、密度(c)；
图5是本发明实施例以图3所示不同频段的角道集作为联合输入地震数据并结合精确Zoeppritz方程反演得到的纵波速度(a)、横波速度(b)、密度(c)；
具体实施方式
针对传统叠前AVA反演对高频地震信息利用不充分，反演结果对具有互层结构等特征的盐见页岩油储层分辨能力弱等问题，本发明的目的是提供一种基于高精度非线性反演算法的高分辨分频联合反演方法，本发明是在研究了以下问题的基础之上提出的：1、有限带宽的地震观测数据分辨能力有限，使得基于该数据的反演结果不能对具有互层结构特征的盐间页岩油层进行有效识别；2、地震数据中包含的高频信息对薄互层具有较好的识别效果，但现有反演方法对这部分信息的利用不充分；3、缺乏合理的频段划分原则，可能造成反演结果失真；4、缺乏有效的权重因子调节机制，可能会造成反演结果突变。本发明首先对实际数据进行分频处理，得到不同频段对应的角道集数据，然后将这些角道集数据作为输入与基于高精度正演算子的非线性反演方法进行结合，最终实现高分辨分频联合反演，合理有效地提升储层参数反演结果的分辨率。
如图1所示，为本发明实施例一种基于高精度非线性反演算法的高分辨分频联合反演方法的工作流程图，具体包括：
101、利用测井数据正演分析确定地震数据频段划分原则，并基于该原则对地震数据进行分频处理：首先利用测井数据正演分析确定地震数据频段划分原则，然后利用广义S变换将时间域的角道集数据变换到频率域，最后按照确定好的分频原则对地震数据进行频段划分并将分离开来的不同频段对应的地震数据反变换回时间域；
S变换可以定义为：
Figure BDA0002712053790000061
其中是地震信号，t和τ是时间，τ表示窗函数的中心点，通过改变它可以控制窗函数的位置，f是频率。在S变换中，基本小波以及高斯窗函数的表达式如下：
Figure BDA0002712053790000062
Figure BDA0002712053790000063
从上述公式可以看出，S变换的时窗宽度与频率的变化趋势相反，也就是说低频段对应的时窗较宽，获得的频率分辨率较高；而高频段对应的时窗相对较窄，获得的时间分辨率较高。S变换不仅有效地克服了短时傅里叶变换窗口固定不变的缺点，同时又引进了小波变换的多分辨率分析，S变换可以说是既结合了它们的优点又很好地克服了它们存在的不足，因此得到了广泛应用。
S变换的窗函数是固定的，其大小不能根据实际需要进行调整，缺乏灵活性，限制了它的适用性。对传统的S变换进行了推广，提出现有的广义S变换。广义S变换的主要思想是通过替换S变换中的基本小波和窗函数，演变出许多不同形式的广义S变换。通过用变量f/γGS代替(2)式中的f对传统形式的S变换进行推广，提出了一种广义S变换。其对应的基函数如下所示：
Figure BDA0002712053790000064
其中，参数γGS为常数，且γGS>0。当γGS取值为1时即是标准的S变换。通过控制这个尺度参数的大小就能达到控制窗口宽度的目的，从而改变频谱的时间分辨率或频率分辨率。由于时间分辨率与频率分辨率具有对立性，因此在利用广义S变换进行实际数据处理应用时需要选择合理的γGS值。广义S变换与傅立叶谱有直接联系，其反变换可以通过求傅立叶反变换实现：
Figure BDA0002712053790000065
102、利用得到的分频地震数据提取角度依赖与频率依赖的地震子波，基于测井数据和精确反射系数方程正演模拟地震角道集并结合实际井旁地震数据确定振幅缩放因子：首先，基于角道集数据和测井数据采取统计方法可以提取若干个依赖于入射角度的地震子波，每一个频段的地震数据均可提取若干角度依赖的地震子波，最终可以得到一组角度依赖与频率依赖的地震子波集；以测井数据为输入模型利用精确反射系数方程正演模拟角度道集，与实际井旁角度域地震道集对比，计算振幅缩放因子，并应用于所提取的地震子波。
本发明假设反演之前地震子波已知，因而，需要基于实际的地震叠前道集和测井数据采取统计方法提取子波。因为基于实际的地震叠前道集和测井数据采取统计方法提取子波的振幅会被归一化或标准化，丢失原本的振幅信息。用这样的子波进行正演得到的地震记录和实际的地震记录在振幅上会较大差距，之后反演步骤中计算的残差会非常大，无法得到正确的迭代更新量。所以为了获取与真实地震资料同一量级的合成地震记录，应当将各个角度和频率的子波与振幅缩放因子相乘，将子波调整到真实振幅量级，这样正演出的地震记录与真实地震记录的差值(残差)，才可以用于计算每一步迭代反演的真实更新量。
103、基于工区内所有测井数据提取得到模型参数及其均值，并计算包含地层参数统计相关性的协方差矩阵：假设模型参数服从柯西分布，利用测井数据求取各参数的自相关系数和互相关系数，构建协方差矩阵，形成符合该工区的模型参数先验分布函数。
为了降低地震反演的不确定性，提高反演过程的稳定性，需要从其它途径获得地下介质地震弹性参数模型的信息作为先验信息。本发明采用柯西分布函数作为先验分布函数，基于工区内所有测井数据得到的各模型参数及其均值，求取各参数的自相关系数和互相关系数，构建协方差矩阵，形成符合该工区的模型参数先验分布函数。在后续的反演目标函数中先验模型相应的规则化表达式为：
Figure BDA0002712053790000071
其中，Φ为包含弹性参数和各向异性参数统计相关性的协方差矩阵，N为模型参数的长度，μ为模型参数的均值向量(不同的模型参数需要分别求取)。
104、利用地震构造解释资料和测井数据，基于沉积模式建立时间域的初始地层参数模型：利用地震构造解释资料，基于沉积模式建立地质模型，并将测井资料，按构造模式进行插值和外推，得到每条测线的初始弹性参数模型，包括常规弹性参数和各向异性参数模型。
建立弹性参数模型主要利用三维空间插值方法，其的技术流程是首先利用散点插值的方法对各个层位的数据进行插值，完成地质层位建模，然后根据地质层位进行弹性参数横向插值，即将测井信息进行横向插值，计算得到地下每个点上的弹性参数值，完成初始弹性参数建模的任务。
105、基于时间域的初始弹性参数模型和角度依赖与频率依赖的地震子波正演模拟不同频段的多波地震叠前角道集，由正演模拟记录和不同频率段的实际记录直接计算对应的叠前反演残差：以时间域的初始模型为输入，利用高计算精度的正演方程计算每一个对应角度的反射系数向量，将频率依赖的地震子波与之进行褶积，获得不同频段对应的合成角度道集，并与不同频段对应的实际角道集数据作差，获得叠前地震反演残差。
VTI介质高精度正演方程推导基于Thomsen参数，Thomsen参数是一组新的各向异性参数，利用两个垂向速度和三个无量纲的各向异性参数对各向异性介质五个独立的刚度张量进行替换，它们之间的关系如下：
Figure BDA0002712053790000081
Figure BDA0002712053790000082
其中，vP0和vS0分别表示垂直方向(对称轴方向)的纵波速度和横波速度，ρ表示密度，γ,ε和δ是Thomsen各向异性参数。Thomsen参数具有明确的物理含义：γ表示的是横向和垂向SH波速度的差异，ε表示的是横向和垂向P波速度的差异，δ描述的是P波速度在近垂直入射时随相位角的变化，由于在VTI介质中SH波从P波和SV波中完全解耦，所以只考虑P波和SV波的反射和透射系数。
基于弱各向异性假设，VTI介质中的P波相速度和刚度c13可以利用各向异性参数ε和δ线性表示为：
VP(θ)≈vP0(1+δsin2θcos2θ+εsin4θ) (7)
c13≈c33(1+δ)-2c55 (8)根据公式(7)，可以得到形式相对简单的水平慢度：
Figure BDA0002712053790000091
将公式(6)和(8)代入VTI精确反射系数方程可得：
Figure BDA0002712053790000092
其中
Figure BDA0002712053790000093
Figure BDA0002712053790000094
Figure BDA0002712053790000095
Figure BDA0002712053790000096
Figure BDA0002712053790000097
Figure BDA0002712053790000098
Figure BDA0002712053790000099
Figure BDA00027120537900000910
其中
Figure BDA0002712053790000101
Figure BDA0002712053790000102
Figure BDA0002712053790000103
Figure BDA0002712053790000104
在上述VTI反射和透射系数方程中，不同波模式的垂直慢度可表示为：
Figure BDA0002712053790000105
Figure BDA0002712053790000106
其中
Figure BDA0002712053790000107
Figure BDA0002712053790000108
Figure BDA0002712053790000109
其中，qα是平面P波的垂直慢度，qβ是对应的SV波的垂直慢度，cij是刚度张量，ρ是密度，p＝sin(θ)/VP(θ)是水平慢度，VP(θ)是P波相速度，θ是相位角。在上述方程中，α表示P波模式，β表示对应的SV波模式。
将公式(6)和(8)代入公式(14)中，将独立的刚度张量利用Thomsen各向异性参数进行替换，然后通过略去高阶项运算可得：
Figure BDA0002712053790000111
Figure BDA0002712053790000112
Figure BDA0002712053790000113
Figure BDA0002712053790000114
将公式(15)中的K1和公式(16)代入公式(13)中可得：
Figure BDA0002712053790000115
Figure BDA0002712053790000116
经过上述推导就得到了形式简单的水平慢度和垂直慢度，使得VTI介质的反射透射系数方程形式得到极大的简化。公式(10)是VTI介质反射和透射系数表达式的隐式表示形式。公式(9)，(10)，(11)和(17)是公式(10)所示VTI介质反射和透射系数隐式表示形式中间变量的显示表达式。将这五个公式联立，便构成本专利所使用的高精度正演算子。
106、基于Bayesian理论，将尺度独立的分频地震数据联合构建得到最大后验概率意义下的反演目标函数，利用广义线性反演思想对反演目标函数进行求解，得到模型参数扰动量的求解表达式：假设地震数据噪声服从Gaussian分布，先验模型服Cauchy分布，然后基于Bayesian原理，将尺度相互独立的分频地震数据联合构建得到最大后验概率意义下的AVA非线性联合反演目标函数。然后采用广义线性反演方法对AVA非线性联合反演目标函数进行求解，求解过程中需要计算正演算子关于模型参数的一阶偏导数，该导数可以通过解析或数值方法求得。
根据实际分辨能力的需要，利用傅里叶变换或广义S变换可将地震数据分成多个频段(dpp1,dpp2,...dppl,dps1,dps2,...dpsl)。基于贝叶斯理论，利用高精度的正演方程和不同频段的地震数据构建AVA非线性联合反演目标函数，并引入合理的先验模型。利用稳定的高精度非线性目标函数求解算法对目标函数进行求解，提高反演结果的分辨率。
在分频联合反演中，利用傅里叶变换或广义S变换对某一个角度对应的地震数据进行分解得到不同频段的地震数据在尺度上相互独立，假设不同频段地震数据包含的噪声数据均服从高斯分布，则分频数据纵横波联合反演的似然函数可表示为：
Figure BDA0002712053790000121
其中，G是高精度正演算子(各向同性介质：精确Zoeppritz方程；VTI介质：精确方程或其它高精度近似公式)，l是所分的频段数。
基于贝叶斯理论，引入合理的先验模型便可得到分频联合反演方法的目标函数：
Figure BDA0002712053790000122
其中，μi,i＝2,...l控制PP波数据不同频率成分的权重，αi,i＝1,2,...l控制PS转换波不同频率成分的权重。同理，上述目标函数也可以利用广义线性反演与迭代重加权最小二乘组合算法进行求解。
采用广义线性反演的思想对上述目标函数进行求解，可得
Δmk＝(Hk)-1γk (20)
其中
Figure BDA0002712053790000131
Figure BDA0002712053790000132
最终，根据上述推导，可以得到模型参数的更新迭代公式：
mk+1＝mk+ηkΔmk,k＝0,1,2... (21)
其中，ηk是第k次迭代的步长因子。
107、根据模型参数扰动量的求解表达式公式和反演残差，采用迭代重加权最小二乘算法计算模型参数的扰动量，同时利用井旁道反演结果选择合理的权重因子，重复迭代以上步骤104、105和106，并通过反演残差控制最大迭代次数，获得高分辨率地层参数反演结果：上述步骤得到的模型参数扰动量的求解表达式是非线性，需要进一步借助迭代重加权最小二乘算法进行求解，同时基于井旁道反演结果选择合理的权重因子。利用以上所述加上扰动量后的参数模型为输入，重复迭代以上所述步骤104、105和106，并通过反演残差和最大迭代次数确定最优的参数模型作为最终的反演结果。从图4和图5所示反演结果对比可以看出，本发明能够进一步有效的提升反演结果的精度和有效性，充分验证了其可行性和有效性。
实施例三
此外，本实施例还提供一种存储介质，其中存储有计算机程序，其特征在于，所述计算机程序被处理器执行时，实现如上所述的基于高精度非线性反演算法的高分辨分频联合反演方法的步骤。
实施例四
另外，本实施例还提供一种计算机设备，包括存储器和处理器，其中所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时，实现如上所述的基于高精度非线性反演算法的高分辨分频联合反演方法的步骤。
需要说明的是，本发明实施例的方法可以由单个设备执行，例如一台计算机或服务器等。本实施例的方法也可以应用于分布式场景下，由多台设备相互配合来完成。在这种分布式场景的情况下，这多台设备中的一台设备可以只执行本发明实施例的方法中的某一个或多个步骤，这多台设备相互之间会进行交互以完成所述的方法。
流程图中或在此以其他方式描述的任何过程或方法描述可以被理解为，表示包括一个或更多个用于实现特定逻辑功能或过程的步骤的可执行指令的代码的模块、片段或部分，并且本发明的优选实施方式的范围包括另外的实现，其中可以不按所示出或讨论的顺序，包括根据所涉及的功能按基本同时的方式或按相反的顺序，来执行功能，这应被本发明的实施例所属技术领域的技术人员所理解。
应当理解，本发明的各部分可以用硬件、软件、固件或它们的组合来实现。在上述实施方式中，多个步骤或方法可以用存储在存储器中且由合适的指令执行系统执行的软件或固件来实现。例如，如果用硬件来实现，和在另一实施方式中一样，可用本领域公知的下列技术中的任一项或他们的组合来实现：具有用于对数据信号实现逻辑功能的逻辑门电路的离散逻辑电路，具有合适的组合逻辑门电路的专用集成电路，可编程门阵列(PGA)，现场可编程门阵列(FPGA)等。
本技术领域的普通技术人员可以理解实现上述实施例方法携带的全部或部分步骤是可以通过程序来指令相关的硬件完成，所述的程序可以存储于一种计算机可读存储介质中，该程序在执行时，包括方法实施例的步骤之一或其组合。
此外，在本发明各个实施例中的各功能单元可以集成在一个处理模块中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个模块中。上述集成的模块既可以采用硬件的形式实现，也可以采用软件功能模块的形式实现。所述集成的模块如果以软件功能模块的形式实现并作为独立的产品销售或使用时，也可以存储在一个计算机可读取存储介质中。
上述提到的存储介质可以是只读存储器，磁盘或光盘等。
总而言之，随着油气资源勘探的不断推进，其对地震反演结果分辨能力的要求也越来越高。基于近似公式的分频联合反演方法以及基于高精度正演方程的非线性反演方法都可以较好地改善反演结果的精度和分辨率，但它们的反演结果还是很难对具有互层结构特征的盐见页岩油层进行有效识别。本技术将上述方法的优势进行融合，提出了高分辨率分频联合反演。首先利用测井数据正演分析确定合理的分频原则并基于该原则对地震数据进行分频处理，然后将分频地震数据与基于高精度正演算子的非线性反演方法相结合，并基于测井数据确定高频信息的权重因子，最终实现基于高精度非线性反演算法的高分辨率分频联合反演。单井数据测试结果充分验证了本发明所述技术方案的可行性和有效性。
虽然本发明所公开的实施方式如上，但所述的内容只是为了便于理解本发明而采用的实施方式，并非用以限定本发明。任何本发明所属技术领域内的技术人员，在不脱离本发明所公开的精神和范围的前提下，可以在实施的形式上及细节上作任何的修改与变化，但本发明的保护范围，仍须以所附的权利要求书所界定的范围为准。

Claims (10)

1.一种基于高精度非线性反演算法的高分辨分频联合反演方法，其特征在于，包括以下步骤：
S101，利用工区测井数据正演分析确定工区地震数据的频段划分原则，并基于该原则对工区地震数据进行分频处理，以获得不同频段的实际地震角道集数据；
S102，利用不同频段的实际地震角道集数据提取角度依赖与频率依赖的地震子波；
S103，基于工区测井数据确定工区模型参数及其均值，并工区模型参数及其均值基于计算包含地层参数统计相关性的协方差矩阵，以形成该工区的模型参数先验分布函数；
S104，利用工区地震构造解释资料和测井数据，建立时间域的初始地层弹性参数模型；
S105，基于时间域的初始地层弹性参数模型和角度依赖与频率依赖的地震子波正演模拟不同频段的多波地震叠前角道集，以获得不同频段的模拟地震角道集数据，通过比较不同频段的模拟地震角道集数据与不同频段的实际地震角道集数据来确定叠前反演残差；
S106，根据贝叶斯理论，利用高精度正演方程和不同频段的实际地震角道集数据，通过引入该工区的模型参数先验分布函数构建最大后验概率意义下的反演目标函数，并对反演目标函数进行求解，以获得模型参数扰动量的求解表达式；
S107，基于模型参数扰动量的求解表达式和叠前反演残差，采用迭代重加权最小二乘算法计算模型参数的扰动量，同时利用工区井旁道反演结果选择相应的权重因子，重复迭代以上步骤S104至S106，并通过叠前反演残差的大小来控制最大迭代次数，从而获得高分辨率地层参数反演结果。
2.如权利要求1所述的基于高精度非线性反演算法的高分辨分频联合反演方法，其特征在于，步骤S101包括：
S101.1，利用工区测井数据正演分析确定工区地震数据的频段划分原则；
S101.2，利用广义S变换将工区地震数据从时间域变换到频率域；
S101.3，按照确定好的频段划分原则对变换到频率域的工区地震数据进行频段划分，并将分离开来的不同频段的地震数据从频率域反变换回时间域。
3.如权利要求1所述的基于高精度非线性反演算法的高分辨分频联合反演方法，其特征在于，步骤S102包括：
利用统计方法，针对不同频段的地震数据，从每一个频段的地震数据提取若干依赖于角度的地震子波，将各个角度和频率的子波与振幅缩放因子相乘，以得到一组角度依赖与频率依赖的地震子波集；
其中，以工区测井数据为输入，利用精确反射系数方程正演模拟角度道集，将模拟的角度道集与工区实际井旁角度域地震道集进行对比，确定所述振幅缩放因子。
4.如权利要求1所述的基于高精度非线性反演算法的高分辨分频联合反演方法，其特征在于，步骤S103包括：
S103.1，假设工区模型参数服从柯西分布，基于工区测井数据确定工区模型参数及其均值，并计算各模型参数的自相关系数和互相关系数；
S103.2，基于各模型参数的自相关系数和互相关系数构建包含地层参数统计相关性的协方差矩阵，以形成该工区的模型参数先验分布函数。
5.如权利要求1所述的基于高精度非线性反演算法的高分辨分频联合反演方法，其特征在于，步骤S104包括：
S104.1，利用工区地震构造解释资料，根据沉积模式建立工区地质模型；
S104.2，将工区测井资料，按照工区地质模型的构造模式进行插值和外推，得到地下每个点上的弹性参数值，从而获得该工区的时间域的初始地层弹性参数模型。
6.如权利要求1所述的基于高精度非线性反演算法的高分辨分频联合反演方法，其特征在于，步骤S105包括：
S105.1，以时间域的初始地层弹性参数模型为输入，利用高计算精度的正演方程计算每一个对应角度的反射系数向量，并将频率依赖的地震子波与反射系数向量进行褶积，从而获得不同频段的模拟地震角道集数据；
S105.2，将不同频段的模拟地震角道集数据与不同频段的实际地震角道集数据作差，获得叠前反演残差。
7.如权利要求1所述的基于高精度非线性反演算法的高分辨分频联合反演方法，其特征在于，步骤S106中，所述反演目标函数为AVA非线性反演目标函数。
8.如权利要求7所述的基于高精度非线性反演算法的高分辨分频联合反演方法，其特征在于，步骤S106中，采用广义线性反演方法对AVA非线性反演目标函数进行求解。
9.一种存储介质，其中存储有计算机程序，其特征在于，所述计算机程序被处理器执行时，实现如权利要求1至8中任一项所述的基于高精度非线性反演算法的高分辨分频联合反演方法的步骤。
10.一种计算机设备，包括存储器和处理器，其中所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时，实现如权利要求1至8中任一项所述的基于高精度非线性反演算法的高分辨分频联合反演方法的步骤。
CN202011060018.2A 2020-09-30 2020-09-30 一种基于高精度非线性反演算法的高分辨分频联合反演方法 Pending CN114428303A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011060018.2A CN114428303A (zh) 2020-09-30 2020-09-30 一种基于高精度非线性反演算法的高分辨分频联合反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011060018.2A CN114428303A (zh) 2020-09-30 2020-09-30 一种基于高精度非线性反演算法的高分辨分频联合反演方法

Publications (1)

Publication Number Publication Date
CN114428303A true CN114428303A (zh) 2022-05-03

Family

ID=81310038

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011060018.2A Pending CN114428303A (zh) 2020-09-30 2020-09-30 一种基于高精度非线性反演算法的高分辨分频联合反演方法

Country Status (1)

Country Link
CN (1) CN114428303A (zh)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2016041189A1 (zh) * 2014-09-19 2016-03-24 杨顺伟 一种评价页岩气储层及寻找甜点区的方法
WO2018107904A1 (zh) * 2016-12-12 2018-06-21 中国石油大学(华东) 一种精确反演杨氏模量和泊松比的方法
CN111522063A (zh) * 2020-04-30 2020-08-11 湖南科技大学 基于分频联合反演的叠前高分辨率流体因子反演方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2016041189A1 (zh) * 2014-09-19 2016-03-24 杨顺伟 一种评价页岩气储层及寻找甜点区的方法
WO2018107904A1 (zh) * 2016-12-12 2018-06-21 中国石油大学(华东) 一种精确反演杨氏模量和泊松比的方法
CN111522063A (zh) * 2020-04-30 2020-08-11 湖南科技大学 基于分频联合反演的叠前高分辨率流体因子反演方法

Similar Documents

Publication Publication Date Title
US11815642B2 (en) Elastic full wavefield inversion with refined anisotropy and VP/VS models
CN104597490B (zh) 基于精确Zoeppritz方程的多波AVO储层弹性参数反演方法
US8296069B2 (en) Pseudo-analytical method for the solution of wave equations
US8335651B2 (en) Estimation of propagation angles of seismic waves in geology with application to determination of propagation velocity and angle-domain imaging
US8082107B2 (en) Methods and computer-readable medium to implement computing the propagation velocity of seismic waves
US7082368B2 (en) Seismic event correlation and Vp-Vs estimation
US11493658B2 (en) Computer-implemented method and system employing nonlinear direct prestack seismic inversion for poisson impedance
US11041971B2 (en) Full wavefield inversion with an image-gather-flatness constraint
Mazzotti et al. Two-grid genetic algorithm full-waveform inversion
CN111522063B (zh) 基于分频联合反演的叠前高分辨率流体因子反演方法
EP2165221A1 (en) Method for velocity analysis using waveform inversion in laplace domain for geophysical imaging
CN103293552A (zh) 一种叠前地震资料的反演方法及系统
CN111025387B (zh) 一种页岩储层的叠前地震多参数反演方法
US20120095690A1 (en) Methods and computer-readable medium to implement inversion of angle gathers for rock physics reflectivity attributes
US11733413B2 (en) Method and system for super resolution least-squares reverse time migration
CN110187384B (zh) 贝叶斯时移地震差异反演方法及装置
US12032111B2 (en) Method and system for faster seismic imaging using machine learning
CN113031068B (zh) 一种基于反射系数精确式的基追踪叠前地震反演方法
WO2017136133A1 (en) Efficient seismic attribute gather generation with data synthesis and expectation method
RU2570827C2 (ru) Гибридный способ для полноволновой инверсии с использованием способа одновременных и последовательных источников
CN116088048A (zh) 包含不确定性分析的各向异性介质多参数全波形反演方法
Sun et al. Intelligent AVA inversion using a convolution neural network trained with pseudo-well datasets
CN108508481B (zh) 一种纵波转换波地震数据时间匹配的方法、装置及系统
US20210255345A1 (en) Velocity Tomography Using Time Lags of Wave Equation Migration
CN111175821B (zh) 一种vti介质的各向异性参数分步反演方法

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