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

1.本发明创造属于油气资源勘探中地震属性预测技术领域，尤其是涉及一种超微分拉普拉斯块约束的叠前多波联合反演方法。


背景技术：

2.地层弹性参数是描述地层性质的重要属性，是地震振幅信息与储层物性参数之间的桥梁。精确的弹性参数（如纵波速度、横波速度和密度等）不仅可以用来进行地层解释，并且可以从中获取地层的物性参数（如孔隙度、渗透率和含水饱和度等），因此叠前弹性参数反演已经成为油气勘探中的一项常规的工作流程。正问题的物理完备性是反问题得以成功求解的必要条件。现有技术中地震叠前正演方法一般为基于zoeppritz方程及其线性近似式的正演，包括aki-richard近似式，hilterman近似式，simth-gidlow近似式，fatti近似式以及gray近似式等，这些线性近似式都是在不同的假设条件下提出的，且方程中所用的地层弹性参数都不尽相同，因此，利用不同的近似式可以直接反演不同的地层弹性参数。相比于精确zoeppritz方程，基于其线性近似式的反演在数学上较为简单。由于线性近似，在所有的近似式中，反射系数与地层弹性参数之间均为线性关系，因此基于线性近似式的反演为线性反演方法；而精确zoeppritz方程为非线性方程，因此基于精确式的反演明显较为复杂。起初，叠前反演技术主要是基于zoeppritz方程近似式的反演，但由于其近似式的诸多假设，例如小角度入射（入射角不超过30
°
），界面两侧弹性参数变化不大等，基于近似式的反演在精确上明显不足，且容易产生假象，因此后来又有学者研究基于精确zoeppritz方程的叠前反演方法。但精确zoeppritz方程及其近似式，均是基于单一界面假设的，不考虑透射损失、转换波及层间多次波，这与地震波在实际地层中的传播情况明显不符，同时也导致反演结果中存在诸多的假象。
3.目前主要的叠前反演技术均是基于纯纵波（pp波）的反演。随着四分量采集技术的发展，转换波（ps波）信息越来越多的受到地球物理工作者的关注。相比于pp波，ps波穿过含流体储层,特别是气云区的能力更强，在地震剖面上能更好的成像，因此更加有利于储层的反演。此外，ps波较之pp波对横波速度和密度信息更加敏感，因此基于pp波和ps波的联合反演不仅可以获得准确的纵波速度信息，还可以更加准确的获取横波速度和密度信息。此外，由于正演算子的病态性，在反演中通常需要加入正则化项来使反演结果稳定。正则化项的选取对于反演结果的精度以及分辨率有着很大的影响，目前最普遍的正则化方式为高斯约束、柯西约束以及微分拉普拉斯约束。相比于高斯约束，后两者对于反演结果的分辨率有很大的提升。尽管如此，柯西约束和微分拉普拉斯约束均有其各自使用的弊端，这些弊端影响着最终反演结果的质量。因此现有技术中的地震叠前反演方法主要存在以下问题：1、正演算子不能正确表征地震波在地下的真实传播过程，这导致反演结果中存在着诸多假象；2、pp波在穿过含流体储层时衰减严重，特别是储层中存在气云区时，因而其地震剖面上对于储层区域的成像效果较差，影响反演结果，此外，pp波对于横波速度和密度的敏感性较低，导致其二者的反演精度较低；3、基于常规的正则化项的反演不能很好的刻画地层的块状性
质，导致反演的分辨率下降，储层边界刻画困难。


技术实现要素：

4.本发明创造要解决的问题是旨在克服上述现有技术中存在的缺陷，提出一种超微分拉普拉斯块约束的叠前多波联合反演方法，以有效提高地震弹性参数反演的精度和有效性。
5.为解决上述技术问题，本发明创造的技术方案是这样实现的：一种超微分拉普拉斯块约束的叠前多波联合反演方法，包括如下步骤：101、基于地震数据提取统计子波；基于测井数据统计该工区纵波速度、横波速度以及密度的空间导数所满足的分布规律，以确定微分拉普拉斯分布中参数的选取；根据贝叶斯理论，模型的先验信息表征了待反演变量的某些特征，因此，首先需要对待反演变量进行统计分析，求取合适的参数，以满足该参数下的分布函数可以包含所有的或者绝大多数的待反演参数。在对实际数据进行反演时，该过程往往需要对测井数据进行统计，通过井数据统计出满足该地区地质条件的分布规律；102、结合测井资料、地震层位和地震构造解释资料，建立符合沉积模式和地震相的初始声波阻抗模型，具体的，用地震构造解释资料，基于沉积模式建立地质模型，并将测井资料，按构造模式进行插值和外推，得到初始的光滑的纵横波速度以及密度模型；103、基于一维波动方程解析解计算地层反射系数，通过褶积子波计算初始合成pp波与ps地震数据，并与真实地震数据进行比较计算残差，以求取正演算子对模型参数的导数；对于模型资料来说，子波是人为给定的，因此直接利用该子波与上述计算的反射系数褶积即可。而对于实际地震资料，如步骤101所述，需要从原始数据上统计子波，通常需要在不同角度的数据上分别进行统计，以获取更加精确的地震子波，从而保证反演结果的准确性；104、基于贝叶斯原理，加入超微分拉普拉斯分布块约束，构建最大后验概率意义下的反演目标函数，以获得反演结果；反问题的求解通常是不适定的，致使反演结果质量变差，因此，加入待反演参数的某种先验信息有利于促使解向先验信息的方向更新。对于叠前三参数反演来说，加入超微分拉普拉斯块约束可以提升反演结果的分辨率，使地层边界的刻画更为清晰；105、根据高斯牛顿法迭代，获取最终的弹性参数反演结果，以用于指导地层解释和储层预测。
6.进一步，步骤102中，建立模型利用空间插值方法。
7.进一步，步骤102中建立模型的方法包括：先利用散点插值的方法对各个层位的数据进行插值，完成地质层位建模；然后根据地质层位进行纵横波速度和密度横向插值，计算得到地下每个点的纵横波速度及密度值，完成初始弹性参数建模。
8.进一步，步骤103中，假设地下介质为n个各向同性层组成的水平层状介质，在第1层介质的顶界面布置一个点震源及若干个检波器，总反射率可以通过一个6维向量递归求得
，
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)每一个向量代表的是第n层以下所有介质的总体响应，
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)其中是第n层介质的传递矩阵，描述了第n层介质对地震波相位的影响，和分别描述了第n个界面对下行波和上行波振幅的影响；假设最后一层介质内只有透射波而没有反射波，因此所有波形的反射率为0，即,
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)从第n层至第1层递归使用式（2），可以得到第0层的向量,
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)进而可以得到所有地层的pp及ps波总反射率，
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)计算得到一定频率、慢度范围的后，利用慢度积分和频率积分可以获得时间—空间域的地震道集，
ꢀꢀ
(6)其中，地震子波的傅里叶变换，是零阶贝赛尔函数，然后对p进行采样，之后对总反射率进行频率积分，
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(7)即得到截距时间—射线参数域地震记录。
9.进一步，对于式（6）中，需对慢度p进行充分采样，以避免出现空间混叠效应。
10.进一步，为避免对p进行过度采样导致计算效率明显下降，允许直接对总反射率进行频率积分。
11.进一步，步骤104中获得反演目标函数的步骤包括：先应用贝叶斯方法加入更加能刻画边界信息的超微分拉普拉斯分布块约束，应用
的数学表达式为，
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(8)其中，是后验概率分布；是从模型参数空间映射到观测数据空间的似然函数；是为先验分布；为边缘概率密度，省略边缘概率密度项后，将贝叶斯方程表述为以下形式，
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)对于超微分拉普拉斯块约束，其先验分布为，
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(10)其中，可表示为，
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(11)在上式中，μ为不同模型参数的均值，通常设定为初始模型，为不同模型参数的缩放因子，d为一阶差分算子，可表示为，
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(12)可表示为，
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(13)而似然函数可表示为，
ꢀꢀꢀꢀꢀꢀ
(14)根据式（9），可以得到基于超微分拉普拉斯块约束的后验概率分布
，(15)因而可以建立反演的目标函数为，(16)最终，基于超微分拉普拉斯块约束的叠前联合反演的目标函数可表示为。(17)进一步，边缘概率密度是一个常数，它的定义规则是使后验概率分布的总和等于1。
12.进一步，在反演中，除了加入刻画边界信息的超微分拉普拉斯分布外，还需要加入高斯分布，以避免反演结果剧烈抖动。
13.进一步，在反演中，目标函数中除了有pp波的数据残差外，还应该包含ps波的数据残差项。
14.本发明创造具有的优点和积极效果是：本发明方法基于一维波动方程解析解可以充分的考虑透射损失、转换波以及层间多次波，从而获得更加精确的反演结果。在传统的基于pp波反演的基础上加入了ps转换波信息，从而获得更加稳定且精确的横波速度和密度反演结果。并且基于贝叶斯理论，通过加入超微分拉普拉斯块约束，对模型参数的导数进行约束，从而对地层得边界信息得到一个更加准确的刻画，获得更高分辨率的反演结果。与传统的柯西约束及拉普拉斯分布相比，新的约束方法理论上对分辨率的提升效果更佳，从而为后续的储层解释以及物性参数反演提供准确的数据体。
附图说明
15.图1是本发明基于超微分拉普拉斯块约束的非线性叠前多波联合反演方法流程图；图2是本发明实施例所使用的超微分拉普拉斯分布与传统的高斯、柯西及拉普拉斯分布的概率密度函数对比（2a）、以及超微分拉普拉斯分布中选取不同的参数p的概率密度函数对比（2b）；图3是本发明方法实施例中输入的无噪叠前pp波（3a）和ps波（3b）地震数据；图4是本发明实施例基于图3的无噪叠前数据所获取的反演结果，其中，虚线为输入的初始模型，实线为真实模型，点线为反演结果；图5是本发明实施例输入的信噪比为10分贝的叠前地震记录；图6是本发明实施例基于图5的无噪叠前数据所获取的反演结果，其中，虚线为输入的初始模型，实线为真实模型，点线为反演结果。
具体实施方式
16.需要说明的是，在不相冲突的情况下，本发明创造中的实施例及实施例中的特征可以相互组合。
17.在本发明创造的描述中，需要理解的是，术语“中心”、“纵向”、“横向”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明创造和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明创造的限制。此外，术语“第一”、“第二”等仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”等的特征可以明示或者隐含地包括一个或者更多个该特征。在本发明创造的描述中，除非另有说明，“多个”的含义是两个或两个以上。
18.在本发明创造的描述中，需要说明的是，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是固定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以通过具体情况理解上述术语在本发明创造中的具体含义。
19.针对现有地震叠前反演方法所存在的不稳定性及正演算子的不完备性等问题，下面对本发明创造的具体实施例做详细说明。如图1至6所示，本发明提供一种超微分拉普拉斯块约束的叠前多波联合反演方法：本发明假设反演之前地震子波已知，因而，需要基于实际的地震叠前道集和测井数据采取统计方法提取子波。另外，加入的超微分拉普拉斯块分布需要包含待反演的模型参数。因此，在反演之前需要通过测井数据对纵横波速度以及密度信息进行统计分析，以确定合适的参数从而使获得的超微分拉普拉斯分布满足该工区地层弹性参数所服从的分布规律。
20.结合测井资料、地震层位和地震构造解释资料，建立符合沉积模式和地震相的初始模型：建立弹性参数模型主要利用空间插值方法，其技术流程是首先利用散点插值的方法对各个层位的数据进行插值，完成地质层位建模，然后根据地质层位进行纵横波速度和密度横向插值，即将测井信息进行横向插值，计算得到地下每个点的纵横波速度及密度值，完成初始弹性参数建模的任务。对于模型数据而言，低频模型（如图4和图6中的虚线所示）是真实弹性参数（如图4和图6中的实线所示）的低通滤波结果。
21.基于一维波动方程解析解计算地层反射系数，通过褶积子波计算初始合成pp波与ps地震数据。假设地下介质为n个各向同性层组成的水平层状介质，在第1层介质的顶界面布置一个点震源及若干个检波器。总反射率（即完全弹性波动方程在频率—慢度域的解）可以通过一个6维向量递归求得，
ꢀꢀ
(1)每一个向量代表的是第n层以下所有介质的总体响应
,(2)其中是第n层介质的传递矩阵，描述了第n层介质对地震波相位的影响，和分别描述了第n个界面对下行波和上行波振幅的影响。
22.首先假设最后一层介质内只有透射波而没有反射波，因此所有波形的反射率为0，即,(3)从第n层至第1层递归使用式（2），可以得到第0层的向量,(4)根据上述对向量的定义，可以得到所有地层的pp及ps波总反射率，，(5)计算得到一定频率、慢度范围的后，利用慢度积分和频率积分可以获得时间—空间域的地震道集，(6)其中，地震子波的傅里叶变换，是零阶贝赛尔函数。然而上式需对慢度p进行充分采样，否则将出现空间混叠效应。然后对p进行过度采样会导致计算效率明显下降，为了节约计算成本，直接对总反射率进行频率积分，(7)即截距时间—射线参数域地震记录。
23.（4）基于贝叶斯原理，加入更加能刻画边界信息的超微分拉普拉斯分布块约束，构建最大后验概率意义下的反演目标函数。贝叶斯方法是用来计算条件概率的一种概率统计方法，数学上它的表达式为，(8)
其中，是后验概率分布；是从模型参数空间映射到观测数据空间的似然函数；是在观测之前对模型参数的先验了解，称为先验分布；是一个常数，称为边缘概率密度，它的定义规则是使后验概率分布的总和等于1。一般省略边缘概率密度项，将贝叶斯方程表述为以下形式，(9)对于超微分拉普拉斯块约束，其先验分布为,(10)其中，可表示为,(11)在上式中，μ为不同模型参数的均值，通常设定为初始模型，为不同模型参数的缩放因子。d为一阶差分算子，可表示为，(12)可表示为，(13)而似然函数可表示为，(14)根据式（9），可以得到基于超微分拉普拉斯块约束的后验概率分布，
ꢀꢀꢀꢀ
(15)因而可以建立反演的目标函数为
，(16)通常在反演中，除了加入刻画边界信息的超微分拉普拉斯分布外，还需要加入高斯分布使得反演结果抖动不是很剧烈。此外，对于联合反演来说，目标函数中除了有pp波的数据残差外，还应该包含ps波的数据残差项，最终，基于超微分拉普拉斯块约束的叠前联合反演的目标函数可表示为。(17)根据高斯牛顿法，我们可以获得最终的纵横波速度及密度反演结果。高斯牛顿法的迭代公式为，(18)其中，为目标函数关于模型参数的一阶导数（梯度），为目标函数关于模型参数的二阶导数项，即海森矩阵，其二者可表示为，(19)，(20)其中，海森矩阵中忽略了正演算子对模型参数的二阶导数项。
24.设定好迭代次数后，就可以获得最终的弹性参数反演结果（如图4和图6中的点线所示）。图4为本发明实施例基于图3的无噪地震数据进行反演获取的弹性参数反演结果，其中，虚线为输入的初始模型，实线为真实模型，点线为反演结果。图6为基于图5的含噪地震数据（信噪比为10）的反演结果，其中，虚线为输入的初始模型，实线为真实模型，点线为反
演结果。
25.本发明是在研究了以下问题的基础之上提出的：（1）正演算子不能正确表征地震波在地下的真实传播过程，这导致反演结果中存在着诸多假象；（2）pp波在穿过含流体储层时衰减严重，特别是储层中存在气云区时，因而其地震剖面上对于储层区域的成像效果较差，影响反演结果，此外，pp波对于横波速度和密度的敏感性较低，导致其二者的反演精度较低；（3）基于常规的正则化项的反演不能很好的刻画地层的块状性质，导致反演的分辨率下降，储层边界刻画困难。本发明首先对一维波动方程解析解进行介绍，指出该正演方法相比于常规的基于zoeppritz方程及其近似式的正演方法得优势。之后指出转换波数据对于储层（尤其是气云区）的成像效果要优于纵波数据，进而对基于一维波动方程解析解的pp和ps波联合反演方法进行介绍。最后，为了精确刻画储层边界，提高反演结果的分辨率以及稳定性，引入超微分拉普拉斯块约束，从分布图像上证明该约束较之传统的高斯、柯西以及拉普拉斯分布的优势，结合贝叶斯理论，通过高斯牛顿法获得最终的高分辨率反演结果，为地层解释及储层预测提供数据基础。
26.综上所述，本发明具有以下优点：1、本技术方案能够提供一种高分辨率的叠前弹性参数反演方法。传统的基于zoeppritz方程及其近似式的反演方法中无法考虑真实地层中地震波的复杂的传播规律，而本技术方案基于一维波动方程解析解可以充分的考虑透射损失、转换波以及层间多次波，从而获得更加精确的反演结果。
27.2、本技术方案在传统的基于pp波反演的基础上加入了ps转换波信息，从而获得更加稳定且精确的横波速度和密度反演结果。
28.3、基于贝叶斯理论，本技术方案通过加入超微分拉普拉斯块约束，对模型参数的导数进行约束，从而对地层得边界信息得到一个更加准确的刻画，获得更高分辨率的反演结果。与传统的柯西约束及拉普拉斯分布相比，新的约束方法理论上对分辨率的提升效果更佳，从而为后续的储层解释以及物性参数反演提供准确的数据体。
29.对于本领域技术人员而言，显然本发明创造不限于上述示范性实施例的细节，而且在不背离本发明创造的精神或基本特征的情况下，能够以其他的具体形式实现本发明创造。
30.因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明创造的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明创造内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。
31.此外，应当理解，虽然本说明书按照实施方式加以描述，但并非每个实施方式仅包含一个独立的技术方案，说明书的这种叙述方式仅仅是为清楚起见，本领域技术人员应当将说明书作为一个整体，各实施例中的技术方案也可以经适当组合，形成本领域技术人员可以理解的其他实施方式。