CN114740528A - 一种超微分拉普拉斯块约束的叠前多波联合反演方法 - Google Patents

一种超微分拉普拉斯块约束的叠前多波联合反演方法 Download PDF

Info

Publication number
CN114740528A
CN114740528A CN202210322566.0A CN202210322566A CN114740528A CN 114740528 A CN114740528 A CN 114740528A CN 202210322566 A CN202210322566 A CN 202210322566A CN 114740528 A CN114740528 A CN 114740528A
Authority
CN
China
Prior art keywords
wave
inversion
seismic
ultramicro
distribution
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
CN202210322566.0A
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.)
CNOOC China Ltd
Original Assignee
CNOOC China 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 CNOOC China Ltd filed Critical CNOOC China Ltd
Priority to CN202210322566.0A priority Critical patent/CN114740528A/zh
Publication of CN114740528A publication Critical patent/CN114740528A/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/30—Analysis
    • 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
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/15—Correlation function computation including computation of convolution operations

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Remote Sensing (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geophysics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Computing Systems (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明提供一种超微分拉普拉斯块约束的叠前多波联合反演方法，先基于地震数据提取统计子波，并基于测井数据统计该工区纵波速度、横波速度以及密度的空间导数所满足的分布规律，确定微分拉普拉斯分布中参数的选取，之后建立符合沉积模式和地震相的初始纵横波速度和密度模型，然后通过褶积子波计算初始合成PP波与PS地震数据，并与真实地震数据进行比较计算残差，进而求取正演算子对模型参数的导数，之后基于贝叶斯原理，加入超微分拉普拉斯分布块约束，构建最大后验概率意义下的反演目标函数，获得反演结果，最后根据高斯牛顿法迭代，获取最终的弹性参数反演结果，指导地层解释和储层预测，有效提高地震弹性参数反演的精度和有效性。

Description

一种超微分拉普拉斯块约束的叠前多波联合反演方法
技术领域
本发明创造属于油气资源勘探中地震属性预测技术领域，尤其是涉及一种超微分拉普拉斯块约束的叠前多波联合反演方法。
背景技术
地层弹性参数是描述地层性质的重要属性，是地震振幅信息与储层物性参数之间的桥梁。精确的弹性参数（如纵波速度、横波速度和密度等）不仅可以用来进行地层解释，并且可以从中获取地层的物性参数（如孔隙度、渗透率和含水饱和度等），因此叠前弹性参数反演已经成为油气勘探中的一项常规的工作流程。正问题的物理完备性是反问题得以成功求解的必要条件。现有技术中地震叠前正演方法一般为基于Zoeppritz方程及其线性近似式的正演，包括Aki-Richard近似式，Hilterman近似式，Simth-Gidlow近似式，Fatti近似式以及Gray近似式等，这些线性近似式都是在不同的假设条件下提出的，且方程中所用的地层弹性参数都不尽相同，因此，利用不同的近似式可以直接反演不同的地层弹性参数。相比于精确Zoeppritz方程，基于其线性近似式的反演在数学上较为简单。由于线性近似，在所有的近似式中，反射系数与地层弹性参数之间均为线性关系，因此基于线性近似式的反演为线性反演方法；而精确Zoeppritz方程为非线性方程，因此基于精确式的反演明显较为复杂。起初，叠前反演技术主要是基于Zoeppritz方程近似式的反演，但由于其近似式的诸多假设，例如小角度入射（入射角不超过30°），界面两侧弹性参数变化不大等，基于近似式的反演在精确上明显不足，且容易产生假象，因此后来又有学者研究基于精确Zoeppritz方程的叠前反演方法。但精确Zoeppritz方程及其近似式，均是基于单一界面假设的，不考虑透射损失、转换波及层间多次波，这与地震波在实际地层中的传播情况明显不符，同时也导致反演结果中存在诸多的假象。
目前主要的叠前反演技术均是基于纯纵波（PP波）的反演。随着四分量采集技术的发展，转换波（PS波）信息越来越多的受到地球物理工作者的关注。相比于PP波，PS波穿过含流体储层,特别是气云区的能力更强，在地震剖面上能更好的成像，因此更加有利于储层的反演。此外，PS波较之PP波对横波速度和密度信息更加敏感，因此基于PP波和PS波的联合反演不仅可以获得准确的纵波速度信息，还可以更加准确的获取横波速度和密度信息。此外，由于正演算子的病态性，在反演中通常需要加入正则化项来使反演结果稳定。正则化项的选取对于反演结果的精度以及分辨率有着很大的影响，目前最普遍的正则化方式为高斯约束、柯西约束以及微分拉普拉斯约束。相比于高斯约束，后两者对于反演结果的分辨率有很大的提升。尽管如此，柯西约束和微分拉普拉斯约束均有其各自使用的弊端，这些弊端影响着最终反演结果的质量。因此现有技术中的地震叠前反演方法主要存在以下问题：1、正演算子不能正确表征地震波在地下的真实传播过程，这导致反演结果中存在着诸多假象；2、PP波在穿过含流体储层时衰减严重，特别是储层中存在气云区时，因而其地震剖面上对于储层区域的成像效果较差，影响反演结果，此外，PP波对于横波速度和密度的敏感性较低，导致其二者的反演精度较低；3、基于常规的正则化项的反演不能很好的刻画地层的块状性质，导致反演的分辨率下降，储层边界刻画困难。
发明内容
本发明创造要解决的问题是旨在克服上述现有技术中存在的缺陷，提出一种超微分拉普拉斯块约束的叠前多波联合反演方法，以有效提高地震弹性参数反演的精度和有效性。
为解决上述技术问题，本发明创造的技术方案是这样实现的：
一种超微分拉普拉斯块约束的叠前多波联合反演方法，包括如下步骤：
101、基于地震数据提取统计子波；基于测井数据统计该工区纵波速度、横波速度以及密度的空间导数所满足的分布规律，以确定微分拉普拉斯分布中参数的选取；根据贝叶斯理论，模型的先验信息表征了待反演变量的某些特征，因此，首先需要对待反演变量进行统计分析，求取合适的参数，以满足该参数下的分布函数可以包含所有的或者绝大多数的待反演参数。在对实际数据进行反演时，该过程往往需要对测井数据进行统计，通过井数据统计出满足该地区地质条件的分布规律；
102、结合测井资料、地震层位和地震构造解释资料，建立符合沉积模式和地震相的初始声波阻抗模型，具体的，用地震构造解释资料，基于沉积模式建立地质模型，并将测井资料，按构造模式进行插值和外推，得到初始的光滑的纵横波速度以及密度模型；
103、基于一维波动方程解析解计算地层反射系数，通过褶积子波计算初始合成PP波与PS地震数据，并与真实地震数据进行比较计算残差，以求取正演算子对模型参数的导数；对于模型资料来说，子波是人为给定的，因此直接利用该子波与上述计算的反射系数褶积即可。而对于实际地震资料，如步骤101所述，需要从原始数据上统计子波，通常需要在不同角度的数据上分别进行统计，以获取更加精确的地震子波，从而保证反演结果的准确性；
104、基于贝叶斯原理，加入超微分拉普拉斯分布块约束，构建最大后验概率意义下的反演目标函数，以获得反演结果；反问题的求解通常是不适定的，致使反演结果质量变差，因此，加入待反演参数的某种先验信息有利于促使解向先验信息的方向更新。对于叠前三参数反演来说，加入超微分拉普拉斯块约束可以提升反演结果的分辨率，使地层边界的刻画更为清晰；
105、根据高斯牛顿法迭代，获取最终的弹性参数反演结果，以用于指导地层解释和储层预测。
进一步，步骤102中，建立模型利用空间插值方法。
进一步，步骤102中建立模型的方法包括：
先利用散点插值的方法对各个层位的数据进行插值，完成地质层位建模；然后根据地质层位进行纵横波速度和密度横向插值，计算得到地下每个点的纵横波速度及密度值，完成初始弹性参数建模。
进一步，步骤103中，假设地下介质为N个各向同性层组成的水平层状介质，在第1层介质的顶界面布置一个点震源及若干个检波器，总反射率
Figure 100002_DEST_PATH_IMAGE001
可以通过一个6维向量
Figure 100002_DEST_PATH_IMAGE002
递归求得
Figure 100002_DEST_PATH_IMAGE003
， (1)
每一个向量代表的是第n层以下所有介质的总体响应
Figure 100002_DEST_PATH_IMAGE004
， (2)
其中
Figure 100002_DEST_PATH_IMAGE005
是第n层介质的传递矩阵，
Figure 100002_DEST_PATH_IMAGE006
描述了第n层介质对地震波相位的影响，
Figure 100002_DEST_PATH_IMAGE007
和
Figure 100002_DEST_PATH_IMAGE008
分别描述了第n个界面对下行波和上行波振幅的影响；
假设最后一层介质内只有透射波而没有反射波，因此所有波形的反射率为0，即
Figure 100002_DEST_PATH_IMAGE009
, (3)
从第N层至第1层递归使用式（2），可以得到第0层的向量
Figure 100002_DEST_PATH_IMAGE010
, (4)
进而可以得到所有地层的PP及PS波总反射率
Figure 100002_DEST_PATH_IMAGE011
， (5)
计算得到一定频率、慢度范围的
Figure 100002_DEST_PATH_IMAGE012
后，利用慢度积分和频率积分可以获得时间—空间域的地震道集
Figure 100002_DEST_PATH_IMAGE013
， (6)
其中，
Figure 100002_DEST_PATH_IMAGE014
地震子波的傅里叶变换，
Figure 100002_DEST_PATH_IMAGE015
是零阶贝赛尔函数，然后对p进行采样，之后对总反射率
Figure 100002_DEST_PATH_IMAGE016
进行频率积分
Figure 100002_DEST_PATH_IMAGE017
， (7)
即得到截距时间—射线参数域地震记录。
进一步，对于式（6）中，需对慢度p进行充分采样，以避免出现空间混叠效应。
进一步，为避免对p进行过度采样导致计算效率明显下降，允许直接对总反射率
Figure 605923DEST_PATH_IMAGE016
进行频率积分。
进一步，步骤104中获得反演目标函数的步骤包括：
先应用贝叶斯方法加入更加能刻画边界信息的超微分拉普拉斯分布块约束，应用的数学表达式为
Figure 100002_DEST_PATH_IMAGE018
， (8)
其中，
Figure 100002_DEST_PATH_IMAGE019
是后验概率分布；
Figure 100002_DEST_PATH_IMAGE020
是从模型参数空间映射到观测数据空间的似然函数；
Figure 100002_DEST_PATH_IMAGE021
是为先验分布；
Figure 100002_DEST_PATH_IMAGE022
为边缘概率密度，省略边缘概率密度项后，将贝叶斯方程表述为以下形式
Figure 100002_DEST_PATH_IMAGE023
， (9)
对于超微分拉普拉斯块约束，其先验分布为
Figure 100002_DEST_PATH_IMAGE024
， (10)
其中，
Figure 100002_DEST_PATH_IMAGE025
可表示为
Figure 100002_DEST_PATH_IMAGE026
， (11)
在上式中，μ为不同模型参数的均值，通常设定为初始模型，
Figure 100002_DEST_PATH_IMAGE027
为不同模型参数的缩放因子，D为一阶差分算子，可表示为
Figure 100002_DEST_PATH_IMAGE028
， (12)
Figure 100002_DEST_PATH_IMAGE029
可表示为
Figure 100002_DEST_PATH_IMAGE030
， (13)
而似然函数可表示为
Figure 100002_DEST_PATH_IMAGE031
， (14)
根据式（9），可以得到基于超微分拉普拉斯块约束的后验概率分布
Figure 100002_DEST_PATH_IMAGE032
， (15)
因而可以建立反演的目标函数为
Figure 100002_DEST_PATH_IMAGE033
， (16)
最终，基于超微分拉普拉斯块约束的叠前联合反演的目标函数可表示为
Figure DEST_PATH_IMAGE034
。 (17)
进一步，边缘概率密度是一个常数，它的定义规则是使后验概率分布的总和等于1。
进一步，在反演中，除了加入刻画边界信息的超微分拉普拉斯分布外，还需要加入高斯分布，以避免反演结果剧烈抖动。
进一步，在反演中，目标函数中除了有PP波的数据残差外，还应该包含PS波的数据残差项。
本发明创造具有的优点和积极效果是：
本发明方法基于一维波动方程解析解可以充分的考虑透射损失、转换波以及层间多次波，从而获得更加精确的反演结果。在传统的基于PP波反演的基础上加入了PS转换波信息，从而获得更加稳定且精确的横波速度和密度反演结果。并且基于贝叶斯理论，通过加入超微分拉普拉斯块约束，对模型参数的导数进行约束，从而对地层得边界信息得到一个更加准确的刻画，获得更高分辨率的反演结果。与传统的柯西约束及拉普拉斯分布相比，新的约束方法理论上对分辨率的提升效果更佳，从而为后续的储层解释以及物性参数反演提供准确的数据体。
附图说明
图1是本发明基于超微分拉普拉斯块约束的非线性叠前多波联合反演方法流程图；
图2是本发明实施例所使用的超微分拉普拉斯分布与传统的高斯、柯西及拉普拉斯分布的概率密度函数对比（2a）、以及超微分拉普拉斯分布中选取不同的参数p的概率密度函数对比（2b）；
图3是本发明方法实施例中输入的无噪叠前PP波（3a）和PS波（3b）地震数据；
图4是本发明实施例基于图3的无噪叠前数据所获取的反演结果，其中，虚线为输入的初始模型，实线为真实模型，点线为反演结果；
图5是本发明实施例输入的信噪比为10分贝的叠前地震记录；
图6是本发明实施例基于图5的无噪叠前数据所获取的反演结果，其中，虚线为输入的初始模型，实线为真实模型，点线为反演结果。
具体实施方式
需要说明的是，在不相冲突的情况下，本发明创造中的实施例及实施例中的特征可以相互组合。
在本发明创造的描述中，需要理解的是，术语“中心”、“纵向”、“横向”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明创造和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明创造的限制。此外，术语“第一”、“第二”等仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”等的特征可以明示或者隐含地包括一个或者更多个该特征。在本发明创造的描述中，除非另有说明，“多个”的含义是两个或两个以上。
在本发明创造的描述中，需要说明的是，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以通过具体情况理解上述术语在本发明创造中的具体含义。
针对现有地震叠前反演方法所存在的不稳定性及正演算子的不完备性等问题，下面对本发明创造的具体实施例做详细说明。如图1至6所示，本发明提供一种超微分拉普拉斯块约束的叠前多波联合反演方法：
本发明假设反演之前地震子波已知，因而，需要基于实际的地震叠前道集和测井数据采取统计方法提取子波。另外，加入的超微分拉普拉斯块分布需要包含待反演的模型参数。因此，在反演之前需要通过测井数据对纵横波速度以及密度信息进行统计分析，以确定合适的参数从而使获得的超微分拉普拉斯分布满足该工区地层弹性参数所服从的分布规律。
结合测井资料、地震层位和地震构造解释资料，建立符合沉积模式和地震相的初始模型：建立弹性参数模型主要利用空间插值方法，其技术流程是首先利用散点插值的方法对各个层位的数据进行插值，完成地质层位建模，然后根据地质层位进行纵横波速度和密度横向插值，即将测井信息进行横向插值，计算得到地下每个点的纵横波速度及密度值，完成初始弹性参数建模的任务。对于模型数据而言，低频模型（如图4和图6中的虚线所示）是真实弹性参数（如图4和图6中的实线所示）的低通滤波结果。
基于一维波动方程解析解计算地层反射系数，通过褶积子波计算初始合成PP波与PS地震数据。假设地下介质为N个各向同性层组成的水平层状介质，在第1层介质的顶界面布置一个点震源及若干个检波器。总反射率
Figure 352424DEST_PATH_IMAGE001
（即完全弹性波动方程在频率—慢度域的解）可以通过一个6维向量
Figure 559591DEST_PATH_IMAGE002
递归求得
Figure 892483DEST_PATH_IMAGE003
， (1)
每一个向量
Figure 100002_DEST_PATH_IMAGE035
代表的是第n层以下所有介质的总体响应
Figure 100002_DEST_PATH_IMAGE036
, (2)
其中
Figure 549598DEST_PATH_IMAGE005
是第n层介质的传递矩阵，
Figure DEST_PATH_IMAGE037
描述了第n层介质对地震波相位的影响，
Figure DEST_PATH_IMAGE038
和
Figure 981848DEST_PATH_IMAGE008
分别描述了第n个界面对下行波和上行波振幅的影响。
首先假设最后一层介质内只有透射波而没有反射波，因此所有波形的反射率为0，即
Figure DEST_PATH_IMAGE039
, (3)
从第N层至第1层递归使用式（2），可以得到第0层的向量
Figure DEST_PATH_IMAGE040
, (4)
根据上述对向量的定义，可以得到所有地层的PP及PS波总反射率，
Figure DEST_PATH_IMAGE041
， (5)
计算得到一定频率、慢度范围的
Figure 297816DEST_PATH_IMAGE012
后，利用慢度积分和频率积分可以获得时间—空间域的地震道集
Figure 699978DEST_PATH_IMAGE013
，(6)
其中，
Figure DEST_PATH_IMAGE042
地震子波的傅里叶变换，
Figure DEST_PATH_IMAGE043
是零阶贝赛尔函数。然而上式需对慢度p进行充分采样，否则将出现空间混叠效应。然后对p进行过度采样会导致计算效率明显下降，为了节约计算成本，直接对总反射率
Figure DEST_PATH_IMAGE044
进行频率积分
Figure DEST_PATH_IMAGE045
， (7)
即截距时间—射线参数域地震记录。
（4）基于贝叶斯原理，加入更加能刻画边界信息的超微分拉普拉斯分布块约束，构建最大后验概率意义下的反演目标函数。贝叶斯方法是用来计算条件概率的一种概率统计方法，数学上它的表达式为
Figure DEST_PATH_IMAGE046
， (8)
其中，
Figure DEST_PATH_IMAGE047
是后验概率分布；
Figure DEST_PATH_IMAGE048
是从模型参数空间映射到观测数据空间的似然函数；
Figure DEST_PATH_IMAGE049
是在观测之前对模型参数的先验了解，称为先验分布；
Figure DEST_PATH_IMAGE050
是一个常数，称为边缘概率密度，它的定义规则是使后验概率分布的总和等于1。一般省略边缘概率密度项，将贝叶斯方程表述为以下形式
Figure DEST_PATH_IMAGE051
， (9)
对于超微分拉普拉斯块约束，其先验分布为
Figure DEST_PATH_IMAGE052
, (10)
其中，
Figure DEST_PATH_IMAGE053
可表示为
Figure DEST_PATH_IMAGE054
, (11)
在上式中，μ为不同模型参数的均值，通常设定为初始模型，
Figure 916327DEST_PATH_IMAGE027
为不同模型参数的缩放因子。D为一阶差分算子，可表示为
Figure DEST_PATH_IMAGE055
， (12)
Figure DEST_PATH_IMAGE056
可表示为
Figure DEST_PATH_IMAGE057
， (13)
而似然函数可表示为
Figure DEST_PATH_IMAGE058
， (14)
根据式（9），可以得到基于超微分拉普拉斯块约束的后验概率分布
Figure DEST_PATH_IMAGE059
， (15)
因而可以建立反演的目标函数为
Figure DEST_PATH_IMAGE060
， (16)
通常在反演中，除了加入刻画边界信息的超微分拉普拉斯分布外，还需要加入高斯分布使得反演结果抖动不是很剧烈。此外，对于联合反演来说，目标函数中除了有PP波的数据残差外，还应该包含PS波的数据残差项，最终，基于超微分拉普拉斯块约束的叠前联合反演的目标函数可表示为
Figure DEST_PATH_IMAGE061
。(17)
根据高斯牛顿法，我们可以获得最终的纵横波速度及密度反演结果。高斯牛顿法的迭代公式为
Figure DEST_PATH_IMAGE062
， (18)
其中，
Figure DEST_PATH_IMAGE063
为目标函数关于模型参数的一阶导数（梯度），
Figure DEST_PATH_IMAGE064
为目标函数关于模型参数的二阶导数项，即海森矩阵，其二者可表示为
Figure DEST_PATH_IMAGE065
， (19)
Figure DEST_PATH_IMAGE066
， (20)
其中，海森矩阵中忽略了正演算子对模型参数的二阶导数项。
设定好迭代次数后，就可以获得最终的弹性参数反演结果（如图4和图6中的点线所示）。图4为本发明实施例基于图3的无噪地震数据进行反演获取的弹性参数反演结果，其中，虚线为输入的初始模型，实线为真实模型，点线为反演结果。图6为基于图5的含噪地震数据（信噪比为10）的反演结果，其中，虚线为输入的初始模型，实线为真实模型，点线为反演结果。
本发明是在研究了以下问题的基础之上提出的：（1）正演算子不能正确表征地震波在地下的真实传播过程，这导致反演结果中存在着诸多假象；（2）PP波在穿过含流体储层时衰减严重，特别是储层中存在气云区时，因而其地震剖面上对于储层区域的成像效果较差，影响反演结果，此外，PP波对于横波速度和密度的敏感性较低，导致其二者的反演精度较低；（3）基于常规的正则化项的反演不能很好的刻画地层的块状性质，导致反演的分辨率下降，储层边界刻画困难。本发明首先对一维波动方程解析解进行介绍，指出该正演方法相比于常规的基于Zoeppritz方程及其近似式的正演方法得优势。之后指出转换波数据对于储层（尤其是气云区）的成像效果要优于纵波数据，进而对基于一维波动方程解析解的PP和PS波联合反演方法进行介绍。最后，为了精确刻画储层边界，提高反演结果的分辨率以及稳定性，引入超微分拉普拉斯块约束，从分布图像上证明该约束较之传统的高斯、柯西以及拉普拉斯分布的优势，结合贝叶斯理论，通过高斯牛顿法获得最终的高分辨率反演结果，为地层解释及储层预测提供数据基础。
综上所述，本发明具有以下优点：
1、本技术方案能够提供一种高分辨率的叠前弹性参数反演方法。传统的基于Zoeppritz方程及其近似式的反演方法中无法考虑真实地层中地震波的复杂的传播规律，而本技术方案基于一维波动方程解析解可以充分的考虑透射损失、转换波以及层间多次波，从而获得更加精确的反演结果。
2、本技术方案在传统的基于PP波反演的基础上加入了PS转换波信息，从而获得更加稳定且精确的横波速度和密度反演结果。
3、基于贝叶斯理论，本技术方案通过加入超微分拉普拉斯块约束，对模型参数的导数进行约束，从而对地层得边界信息得到一个更加准确的刻画，获得更高分辨率的反演结果。与传统的柯西约束及拉普拉斯分布相比，新的约束方法理论上对分辨率的提升效果更佳，从而为后续的储层解释以及物性参数反演提供准确的数据体。
对于本领域技术人员而言，显然本发明创造不限于上述示范性实施例的细节，而且在不背离本发明创造的精神或基本特征的情况下，能够以其他的具体形式实现本发明创造。
因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明创造的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明创造内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。
此外，应当理解，虽然本说明书按照实施方式加以描述，但并非每个实施方式仅包含一个独立的技术方案，说明书的这种叙述方式仅仅是为清楚起见，本领域技术人员应当将说明书作为一个整体，各实施例中的技术方案也可以经适当组合，形成本领域技术人员可以理解的其他实施方式。

Claims (10)

1.一种超微分拉普拉斯块约束的叠前多波联合反演方法，其特征在于，包括如下步骤：
101、基于地震数据提取统计子波；基于测井数据统计该工区纵波速度、横波速度以及密度的空间导数所满足的分布规律，以确定微分拉普拉斯分布中参数的选取；
102、结合测井资料、地震层位和地震构造解释资料，建立符合沉积模式和地震相的初始纵横波速度和密度模型；
103、基于一维波动方程解析解计算地层反射系数，通过褶积子波计算初始合成PP波与PS地震数据，并与真实地震数据进行比较计算残差，以求取正演算子对模型参数的导数；
104、基于贝叶斯原理，加入超微分拉普拉斯分布块约束，构建最大后验概率意义下的反演目标函数，以获得反演结果；
105、根据高斯牛顿法迭代，获取最终的弹性参数反演结果，以用于指导地层解释和储层预测。
2.根据权利要求1所述的一种超微分拉普拉斯块约束的叠前多波联合反演方法，其特征在于：步骤102中，建立模型利用空间插值方法。
3.根据权利要求1或2所述的一种超微分拉普拉斯块约束的叠前多波联合反演方法，其特征在于，步骤102中建立模型的方法包括：
先利用散点插值的方法对各个层位的数据进行插值，完成地质层位建模；然后根据地质层位进行纵横波速度和密度横向插值，计算得到地下每个点的纵横波速度及密度值，完成初始弹性参数建模。
4.根据权利要求1所述的一种超微分拉普拉斯块约束的叠前多波联合反演方法，其特征在于：步骤103中，假设地下介质为N个各向同性层组成的水平层状介质，在第1层介质的顶界面布置一个点震源及若干个检波器，总反射率
Figure DEST_PATH_IMAGE001
可以通过一个6维向量
Figure DEST_PATH_IMAGE002
递归求得
Figure DEST_PATH_IMAGE003
， (1)
每一个向量
Figure DEST_PATH_IMAGE004
代表的是第n层以下所有介质的总体响应
Figure DEST_PATH_IMAGE005
， (2)
其中
Figure DEST_PATH_IMAGE006
是第n层介质的传递矩阵，
Figure DEST_PATH_IMAGE007
描述了第n层介质对地震波相位的影响，
Figure DEST_PATH_IMAGE008
和
Figure DEST_PATH_IMAGE009
分别描述了第n个界面对下行波和上行波振幅的影响；
假设最后一层介质内只有透射波而没有反射波，因此所有波形的反射率为0，即
Figure DEST_PATH_IMAGE010
, (3)
从第N层至第1层递归使用式（2），可以得到第0层的向量
Figure DEST_PATH_IMAGE011
, (4)
进而可以得到所有地层的PP及PS波总反射率
Figure DEST_PATH_IMAGE012
， (5)
计算得到一定频率、慢度范围的
Figure DEST_PATH_IMAGE013
后，利用慢度积分和频率积分可以获得时间—空间域的地震道集
Figure DEST_PATH_IMAGE014
， (6)
其中，
Figure DEST_PATH_IMAGE015
地震子波的傅里叶变换，
Figure DEST_PATH_IMAGE016
是零阶贝赛尔函数，然后对p进行采样，之后对总反射率
Figure DEST_PATH_IMAGE017
进行频率积分
Figure DEST_PATH_IMAGE018
， (7)
即得到截距时间—射线参数域地震记录。
5.根据权利要求4所述的一种超微分拉普拉斯块约束的叠前多波联合反演方法，其特征在于：对于式（6）中，需对慢度p进行充分采样，以避免出现空间混叠效应。
6.根据权利要求4所述的一种超微分拉普拉斯块约束的叠前多波联合反演方法，其特征在于：为避免对p进行过度采样导致计算效率明显下降，允许直接对总反射率
Figure DEST_PATH_IMAGE019
进行频率积分。
7.根据权利要求1所述的一种超微分拉普拉斯块约束的叠前多波联合反演方法，其特征在于，步骤104中获得反演目标函数的步骤包括：
先应用贝叶斯方法加入更加能刻画边界信息的超微分拉普拉斯分布块约束，应用的数学表达式为
Figure DEST_PATH_IMAGE020
， (8)
其中，
Figure DEST_PATH_IMAGE021
是后验概率分布；
Figure DEST_PATH_IMAGE022
是从模型参数空间映射到观测数据空间的似然函数；
Figure DEST_PATH_IMAGE023
是为先验分布；
Figure DEST_PATH_IMAGE024
为边缘概率密度，省略边缘概率密度项后，将贝叶斯方程表述为以下形式
Figure DEST_PATH_IMAGE025
， (9)
对于超微分拉普拉斯块约束，其先验分布为
Figure DEST_PATH_IMAGE026
， (10)
其中，
Figure DEST_PATH_IMAGE027
可表示为
Figure DEST_PATH_IMAGE028
， (11)
在上式中，μ为不同模型参数的均值，通常设定为初始模型，
Figure DEST_PATH_IMAGE029
为不同模型参数的缩放因子，D为一阶差分算子，可表示为
Figure DEST_PATH_IMAGE030
， (12)
Figure DEST_PATH_IMAGE031
可表示为
Figure DEST_PATH_IMAGE032
， (13)
而似然函数可表示为
Figure DEST_PATH_IMAGE033
， (14)
根据式（9），可以得到基于超微分拉普拉斯块约束的后验概率分布
Figure 961846DEST_PATH_IMAGE034
， (15)
因而可以建立反演的目标函数为
Figure DEST_PATH_IMAGE035
， (16)
最终，基于超微分拉普拉斯块约束的叠前联合反演的目标函数可表示为
Figure DEST_PATH_IMAGE036
(17)
8.根据权利要求7所述的一种超微分拉普拉斯块约束的叠前多波联合反演方法，其特征在于：边缘概率密度是一个常数，它的定义规则是使后验概率分布的总和等于1。
9.根据权利要求7所述的一种超微分拉普拉斯块约束的叠前多波联合反演方法，其特征在于：在反演中，除了加入刻画边界信息的超微分拉普拉斯分布外，还需要加入高斯分布，以避免反演结果剧烈抖动。
10.根据权利要求7所述的一种超微分拉普拉斯块约束的叠前多波联合反演方法，其特征在于：在反演中，目标函数中除了有PP波的数据残差外，还应该包含PS波的数据残差项。
CN202210322566.0A 2022-03-30 2022-03-30 一种超微分拉普拉斯块约束的叠前多波联合反演方法 Pending CN114740528A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210322566.0A CN114740528A (zh) 2022-03-30 2022-03-30 一种超微分拉普拉斯块约束的叠前多波联合反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210322566.0A CN114740528A (zh) 2022-03-30 2022-03-30 一种超微分拉普拉斯块约束的叠前多波联合反演方法

Publications (1)

Publication Number Publication Date
CN114740528A true CN114740528A (zh) 2022-07-12

Family

ID=82277406

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210322566.0A Pending CN114740528A (zh) 2022-03-30 2022-03-30 一种超微分拉普拉斯块约束的叠前多波联合反演方法

Country Status (1)

Country Link
CN (1) CN114740528A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117434606A (zh) * 2023-12-21 2024-01-23 中国海洋大学 基于改进的拉普拉斯滤波逆时偏移的地震剖面去噪方法

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117434606A (zh) * 2023-12-21 2024-01-23 中国海洋大学 基于改进的拉普拉斯滤波逆时偏移的地震剖面去噪方法
CN117434606B (zh) * 2023-12-21 2024-03-29 中国海洋大学 基于改进的拉普拉斯滤波逆时偏移的地震剖面去噪方法

Similar Documents

Publication Publication Date Title
CN111208561B (zh) 基于时变子波与曲波变换约束的地震声波阻抗反演方法
Velis Stochastic sparse-spike deconvolution
US9470811B2 (en) Creating a high resolution velocity model using seismic tomography and impedance inversion
US9784868B2 (en) Method and apparatus for deghosting seismic data
RU2460095C2 (ru) Преобразование радона волнового фронта
US7082368B2 (en) Seismic event correlation and Vp-Vs estimation
Sun et al. 2-D wavepath migration
US11733413B2 (en) Method and system for super resolution least-squares reverse time migration
FR2843202A1 (fr) Methode pour former un modele representatif de la distribution d'une grandeur physique dans une zone souterraine, affranchi de l'effet de bruits correles entachant des donnees d'exploration
CN103293552A (zh) 一种叠前地震资料的反演方法及系统
US10310117B2 (en) Efficient seismic attribute gather generation with data synthesis and expectation method
EA032186B1 (ru) Сейсмическая адаптивная фокусировка
US20180017692A1 (en) Device and method for estimating pre-stack wavelet model from seismic gathers
Liu et al. Effects of conjugate gradient methods and step-length formulas on the multiscale full waveform inversion in time domain: Numerical experiments
CN109143352B (zh) 一种各向异性介质地震反射特征方程建立方法
CN114740528A (zh) 一种超微分拉普拉斯块约束的叠前多波联合反演方法
CN114428343A (zh) 基于归一化互相关的Marchenko成像方法及系统
US8380440B2 (en) 3D residual binning and flatness error correction
WO2022256666A1 (en) Method and system for reflection-based travel time inversion using segment dynamic image warping
US12000971B2 (en) Method and system for seismic processing using virtual trace bins based on offset attributes and azimuthal attributes
CN112213774B (zh) 一种浅层q模型估计方法及装置
Dickens et al. Spatial resolution of diffraction tomography
CN109490964A (zh) 一种改进的高精度avo弹性参数快速反演方法
US20230184975A1 (en) Method and system for determining seismic velocities using global path tracing
US20240184007A1 (en) Method and system for determining attenuated seismic time

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