基于反射率法的多波AVO储层弹性参数反演方法及系统

摘要

本发明提供一种基于反射率法的多波AVO储层弹性参数反演方法及系统，该方法包括：采集地震叠前道集、测井数据以及实际井旁角度域地震道集；基于地震叠前道集、测井数据以及实际井旁角度域地震道集确定地震子波以及振幅缩放因子；采集测井数据的统计模型参数；根据统计模型参数确定符合该工区的模型参数先验分布函数；采集地震构造解释资料；根据地震构造解释资料以及测井数据建立深度域的初始弹性参数模型；根据初始弹性参数模型确定PP波与PS波反演残差；构建最大后验概率意义下的反演目标函数；根据反演目标函数以及PP波与PS波反演残差确定最优弹性参数模型。该方案满足地震叠前反演识别油气储层，特别是页岩气储层表征的要求。

说明书

技术领域

本发明关于地球物理勘探技术领域，特别是关于油气田、页岩气地震勘探和储层参数 的处理技术，具体的讲是一种基于反射率法的多波AVO储层弹性参数反演方法及系统。

背景技术

地震反演是获得地下介质内部图像、对储层进行精细描述的有效方法，也是高分辨率 地震勘探的最终表现形式，地震数据反演很大程度上提高了储层表征的价值。

随着地震研究的重点由勘探逐渐向开发转移，由常规油气向非常规油气尤其是页岩气 的转移，通过地震反演等手段来揭示地下油气藏的精细分布特征对油气藏储层进行精雕细 刻的技术也受到越来越多的关注。从目前的研究来看，地震反演研究主要涉及基于 Zoeppritz方程的振幅随偏移距的变化(Amplitude Versus Offset，AVO)反演方法和全波形 反演方法。富含有机质高度非均质性页岩气薄层之间的层间多次波等干扰使得页岩气储层 层间反射往往十分微弱，变化规律复杂，加上页岩气储层的调谐作用也会改变地震振幅， 使这类储层的AVO研究具有很大的挑战性，尤其储层较薄时。

基于佐布里兹Zoeppritz方程或其近似公式的AVO反演方法利用叠前道集上振幅随偏 移距变化的信息来提取弹性参数，但这类方法在应用时不考虑多次波、转换波、球面扩散 和透射损失等波的传播效应，储层弹性参数反演结果误差大，不能页岩气储层表征精度要 求。而全波形反演虽然能够利用全波场的信息来预测模型参数，但其计算量巨大，在反演 尺度和计算效率上不能满足实际油藏储层精细表征要求，尤其对于实际三维大偏移距多波 多分量地震数据。

总之，目前基于地震数据的储层弹性参数反演方法研究存在的主要问题是：

1、传统的叠前AVO反演方法一般只考虑纵波PP波数据，PP波数据主要包含纵波速 度信息，但是对模型速度和密度不敏感，单独利用PP地震AVO数据反演横波速度信息和 密度具有较大的不确定性。

2、PP波与横波PS波旅行时间差异大，时间匹配困难大。

3、振幅的准确性对于叠前反演非常关键，但是传统的叠前AVO反演并不考虑透射损 失等波的传播效应，从而在深层位置模拟数据振幅与实际不匹配。

4、油气藏储层特别是页岩气储层最常见的结构是泥砂岩薄互层结构，这种地质条件 下层间多次波发育，各层反射波之间相互干涉叠加，极大地改变了叠前AVO响应特征， 基于Zoeppritz方程或其近似表达式的AVO反演方法不再适用于这类储层。

因此，如何开发出一种新的AVO储层反演处理方案，其能有效获得地下介质弹性参 数模型信息作为先验信息，以降低反演的不确定性，提高反演的精度是本领域亟待解决的 技术难题。

发明内容

为了克服现有技术存在的上述技术问题，本发明提供了一种基于反射率法的多波AVO 储层弹性参数反演方法及系统，通过确定振幅缩放因子、基于工区内所有测井数据统计模 型参数的先验信息、建立初始弹性参数模型、模拟角度域的多波地震叠前道集并计算反演 残差、更新参数模型直至反演残差达到最大迭代次数，进而输出最终模型参数，以满足地 震叠前反演识别油气储层，特别是页岩气储层表征的要求。

本发明的目的之一是，提供一种基于反射率法的多波AVO储层弹性参数反演方法， 包括：采集工区内的地震叠前道集、测井数据以及实际井旁角度域地震道集；基于所述的 地震叠前道集、测井数据以及实际井旁角度域地震道集确定地震子波以及振幅缩放因子； 采集工区内测井数据的统计模型参数；根据所述的统计模型参数确定符合该工区的模型参 数先验分布函数；采集地震构造解释资料；根据所述的地震构造解释资料以及测井数据建 立深度域的初始弹性参数模型；根据所述深度域的初始弹性参数模型确定PP波与PS波反 演残差；构建最大后验概率意义下的反演目标函数；根据所述的反演目标函数以及PP波 与PS波反演残差确定最优弹性参数模型。

在本发明的优选实施方式中，基于所述的地震叠前道集、测井数据以及实际井旁角度 域地震道集确定地震子波以及振幅缩放因子包括：根据所述的地震叠前道集以及测井数 据、采用统计方法提取依赖于入射角度的地震子波；将所述的测井数据作为输入模型；采 用反射法正演模拟角度域的PP波与PS波道集；将所述的PP波与PS波道集与所述的实际 井旁角度域地震道集进行对比，得到振幅缩放因子；将所述的振幅缩放因子应用于所述的 地震子波。

在本发明的优选实施方式中，根据所述的统计模型参数确定符合该工区的模型参数先 验分布函数包括：所述的统计模型参数满足三变量高斯分布；分析所述的测井数据，得到 所述统计模型参数的均值；确定所述统计模型参数的自相关系数以及互相关系数；根据所 述的自相关系数以及互相关系数构建三参数相关的协方差矩阵；根据所述的三参数相关的 协方差矩阵形成符合该工区的模型参数先验分布函数。

在本发明的优选实施方式中，根据所述的地震构造解释资料以及测井数据建立深度域 的初始弹性参数模型包括：根据所述的地震构造解释资料、基于沉积模式建立地质模型； 将所述的测井资料按照构造模式进行插值和外推，得到每条测线的时间域的初始弹性参数 模型，所述的初始弹性参数模型包括纵波速度模型、横波速度模型以及密度模型；通过时 深转换关系将所述时间域的初始弹性参数模型转换为深度域的初始弹性参数模型。

在本发明的优选实施方式中，根据所述深度域的初始弹性参数模型以及地震子波确定 PP波与PS波反演残差包括：将所述深度域的初始弹性参数模型作为输入；在频率域通过 传递矩阵算法确定每一个角度的PP波与PS波总反射率矩阵，所述PP波与PS波总反射率 矩阵包含反射系数与传播时间、层间多次波与透射损失；对所述PP波与PS波总反射率矩 阵进行傅里叶反变换；将傅里叶反变换后的PP波与PS波总反射率矩阵与地震子波进行褶 积，得到角度域的PP波与PS波角道集；采集实际角度域地震道集；将所述角度域的PP 波与PS波角道集与实际角度域地震道集作差，得到PP波与PS波反演残差。

在本发明的优选实施方式中，构建最大后验概率意义下的反演目标函数包括：利用所 述模型参数先验分布函数对反演结果进行约束；根据Bayesian原理、综合反演似然函数以 及先验分布函数得到后验概率分布函数；根据所述的后验概率分布函数构建反演目标函 数。

在本发明的优选实施方式中，根据所述的反演目标函数以及PP波与PS波反演残差确 定最优弹性参数模型包括：确定所述反演目标函数的梯度以及Hessian矩阵；根据所述的 梯度、Hessian矩阵以及PP波与PS波反演残差，采用高斯-牛顿迭代法确定所述深度域的 初始弹性参数模型的更新量；将所述深度域的初始弹性参数模型的更新量作为输入；根据 所述深度域的初始弹性参数模型的更新量以及地震子波确定更新后的PP波与PS波反演残 差；根据更新后的PP波与PS波反演残差构建更新后的反演目标函数；根据更新后的反演 目标函数确定更新后的梯度以及Hessian矩阵；通过所述更新后的PP波与PS波反演残差 确定最大迭代次数；根据所述的最大迭代次数以及更新后的反演目标函数的梯度以及 Hessian矩阵确定最优弹性参数模型，所述的最优弹性参数模型包括纵波速度、横波速度和 密度：将所述的最优弹性参数模型作为地震三参数叠前AVO反演结果。

本发明的目的之一是，提供了一种基于反射率法的多波AVO储层弹性参数反演的系 统，包括：地震道集采集装置，用于采集工区内的地震叠前道集、测井数据以及实际井旁 角度域地震道集；地震子波确定装置，用于基于所述的地震叠前道集、测井数据以及实际 井旁角度域地震道集确定地震子波以及振幅缩放因子；统计模型参数采集装置，用于采集 工区内测井数据的统计模型参数；分布函数确定装置，用于根据所述的统计模型参数确定 符合该工区的模型参数先验分布函数；解释资料采集装置，用于采集地震构造解释资料； 弹性参数模型建立装置，用于根据所述的地震构造解释资料以及测井数据建立深度域的初 始弹性参数模型；反演残差确定装置，用于根据所述深度域的初始弹性参数模型确定PP 波与PS波反演残差；目标函数构建装置，用于构建最大后验概率意义下的反演目标函数； 最优弹性参数模型确定装置，用于根据所述的反演目标函数以及PP波与PS波反演残差确 定最优弹性参数模型。

在本发明的优选实施方式中，所述的地震子波确定装置包括：地震子波确定模块，用 于根据所述的地震叠前道集以及测井数据、采用统计方法提取依赖于入射角度的地震子 波；输入模型确定模块，用于将所述的测井数据作为输入模型；道集模拟模块，用于采用 反射法正演模拟角度域的PP波与PS波道集；振幅缩放因子确定模块，用于将所述的PP 波与PS波道集与所述的实际井旁角度域地震道集进行对比，得到振幅缩放因子；应用模 块，用于将所述的振幅缩放因子应用于所述的地震子波。

在本发明的优选实施方式中，所述的分布函数确定装置包括：高斯分布模块，所述的 统计模型参数满足三变量高斯分布；均值确定模块，用于分析所述的测井数据，得到所述 统计模型参数的均值；相关系数确定模块，用于确定所述统计模型参数的自相关系数以及 互相关系数；协方差矩阵构建模块，用于根据所述的自相关系数以及互相关系数构建三参 数相关的协方差矩阵；分布函数确定模块，用于根据所述的三参数相关的协方差矩阵形成 符合该工区的模型参数先验分布函数。

在本发明的优选实施方式中，所述的弹性参数模型建立装置包括：地质模型建立模块， 用于根据所述的地震构造解释资料、基于沉积模式建立地质模型；差值外推模块，用于将 所述的测井资料按照构造模式进行插值和外推，得到每条测线的时间域的初始弹性参数模 型，所述的初始弹性参数模型包括纵波速度模型、横波速度模型以及密度模型；转换模块， 用于通过时深转换关系将所述时间域的初始弹性参数模型转换为深度域的初始弹性参数 模型。

在本发明的优选实施方式中，所述的反演残差确定装置包括：输入确定模块，用于将 所述深度域的初始弹性参数模型作为输入；总反射率矩阵确定模块，用于在频率域通过传 递矩阵算法确定每一个角度的PP波与PS波总反射率矩阵，所述PP波与PS波总反射率矩 阵包含反射系数与传播时间、层间多次波与透射损失；傅里叶反变换模块，用于对所述PP 波与PS波总反射率矩阵进行傅里叶反变换；褶积模块，用于将傅里叶反变换后的PP波与 PS波总反射率矩阵与地震子波进行褶积，得到角度域的PP波与PS波角道集；采集模块， 用于采集实际角度域地震道集；反演残差确定模块，用于将所述角度域的PP波与PS波角 道集与实际角度域地震道集作差，得到PP波与PS波反演残差。

在本发明的优选实施方式中，所述的目标函数构建装置包括：约束模块，用于利用所 述模型参数先验分布函数对反演结果进行约束；分布函数确定模块，用于根据Bayesian原 理、综合反演似然函数以及先验分布函数得到后验概率分布函数；目标函数构建模块，用 于根据所述的后验概率分布函数构建反演目标函数。

在本发明的优选实施方式中，所述的最优弹性参数模型确定装置包括：梯度确定模块， 用于确定所述反演目标函数的梯度以及Hessian矩阵；更新量确定模块，用于根据所述的 梯度、Hessian矩阵以及PP波与PS波反演残差，采用高斯-牛顿迭代法确定所述深度域的 初始弹性参数模型的更新量；更新量输入模块，用于将所述深度域的初始弹性参数模型的 更新量作为输入；更新残差确定模块，用于根据所述深度域的初始弹性参数模型的更新量 以及地震子波确定更新后的PP波与PS波反演残差；更新目标函数构建模块，用于根据更 新后的PP波与PS波反演残差构建更新后的反演目标函数；更新梯度确定模块，用于根据 更新后的反演目标函数确定更新后的梯度以及Hessian矩阵；最大迭代次数确定模块，用 于通过所述更新后的PP波与PS波反演残差确定最大迭代次数；最优弹性参数模型确定模 块，用于根据所述的最大迭代次数以及更新后的反演目标函数的梯度以及Hessian矩阵确 定最优弹性参数模型，所述的最优弹性参数模型包括纵波速度、横波速度和密度：反演结 果确定模块，用于将所述的最优弹性参数模型作为地震三参数叠前AVO反演结果。

本发明的有益效果在于，提供了一种基于反射率法的多波AVO储层弹性参数反演方 法及系统，相对于Zoeppritz方程而言，该方案能够模拟水平层状介质情况下除直达波和折 射波以外的全波场响应，从而能够同时利用多种信息进行储层参数预测；该方案联合PP 波与PS波进行三参数反演，由于PS波地震数据包含更多的横波速度和密度信息，该方案 能有效提高反演的稳定性；该方案无需考虑纵横波时间匹配问题；层间多次波干扰地震数 据的分辨能力，而且难于在不损伤有效波的前提下去除，该方案能够有效利用层间多次波 进行模型参数预测，适用于层间多次波发育的储层；该方案考虑了透射损失，使深层反射 波的振幅与实际情况相符合，有利于合成记录与实际记录的振幅匹配；该方案通过三变量 高斯先验分布引入模型参数之间的相关性，在降低反演不确定性的同时，提高了反演的精 度。

为让本发明的上述和其他目的、特征和优点能更明显易懂，下文特举较佳实施例，并 配合所附图式，作详细说明如下。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技 术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明 的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根 据这些附图获得其他的附图。

图1为本发明实施例提供的一种基于反射率法的多波AVO储层弹性参数反演方法的 流程图；

图2为图1中的步骤S102的具体流程图；

图3为图1中的步骤S104的具体流程图；

图4为图1中的步骤S106的具体流程图；

图5为图1中的步骤S107的具体流程图；

图6为图1中的步骤S108的具体流程图；

图7为图1中的步骤S109的具体流程图；

图8为本发明实施例提供的一种基于反射率法的多波AVO储层弹性参数反演系统的 结构框图；

图9为本发明实施例提供的一种基于反射率法的多波AVO储层弹性参数反演系统中 地震子波确定装置的结构框图；

图10为本发明实施例提供的一种基于反射率法的多波AVO储层弹性参数反演系统中 的分布函数确定装置的结构框图；

图11为本发明实施例提供的一种基于反射率法的多波AVO储层弹性参数反演系统中 的弹性参数模型建立装置的结构框图；

图12为本发明实施例提供的一种基于反射率法的多波AVO储层弹性参数反演系统中 反演残差确定装置的结构框图；

图13为本发明实施例提供的一种基于反射率法的多波AVO储层弹性参数反演系统中 的目标函数构建装置的结构框图；

图14为本发明实施例提供的一种基于反射率法的多波AVO储层弹性参数反演系统中 的最优弹性参数模型确定装置的结构框图；

图15为本发明提供的具体实施例中基于反射率法的多波AVO反演方法的流程图；

图16为本发明提供的具体实施例中输入的叠前AVO PP波角道集地震数据示意图；

图17为本发明提供的具体实施例中输入的叠前AVO PS波角道集地震数据示意图；

图18为本发明提供的具体实施例中无噪声情况下采用基于反射率法的多波AVO反演 方法得到的纵波速度模型示意图；

图19为本发明提供的具体实施例中无噪声情况下采用基于反射率法的多波AVO反演 方法得到的横波速度模型示意图；

图20为本发明提供的具体实施例中无噪声情况下采用基于反射率法的多波AVO反演 方法得到的密度模型示意图；

图21为本发明提供的具体实施例中信噪比为4:1的情况下采用基于反射率法的多波 AVO反演方法得到的纵波速度模块示意图；

图22为本发明提供的具体实施例中信噪比为4:1的情况下采用基于反射率法的多波 AVO反演方法得到的横波速度模块示意图；

图23为本发明提供的具体实施例中信噪比为4:1的情况下采用基于反射率法的多波 AVO反演方法得到的密度模型示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地 描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本 发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实 施例，都属于本发明保护的范围。

本发明实施例提供一种基于反射率法的多波AVO储层弹性参数反演方法，以满足地震 叠前反演识别油气储层，特别是页岩气储层表征的要求。

图1为本发明实施例提供的一种基于反射率法的多波AVO储层弹性参数反演方法的流 程图，由图1可知，所述的方法包括：

S101：采集工区内的地震叠前道集、测井数据以及实际井旁角度域地震道集。

S102：基于所述的地震叠前道集、测井数据以及实际井旁角度域地震道集确定地震子 波以及振幅缩放因子。本发明假设反演之前地震子波已经得到，因而需要基于实际的地震 叠前道集和测井数据采取统计方法提取子波，子波在传播过程中受地层的影响会发生波形 或频率变化，提取依赖于入射角度的地震子波能有效提高振幅匹配程度。图2为步骤S102 的具体流程图，由图2可知，该步骤具体包括：

S201：根据所述的地震叠前道集以及测井数据、采用统计方法提取依赖于入射角度的 地震子波。在具体的实施方式中，基于实际的地震叠前道集和测井数据采取统计方法提取 若干个依赖于入射角度的地震子波。

S202：将所述的测井数据作为输入模型；

S203：采用反射法正演模拟角度域的PP波与PS波道集；

S204：将所述的PP波与PS波道集与所述的实际井旁角度域地震道集进行对比，得到 振幅缩放因子；

S205：将所述的振幅缩放因子应用于所述的地震子波。

在具体的实施方式中，基于测井数据和反射法正演模拟地震角道集并结合实际井旁地 震数据确定振幅缩放因子。实际的地震振幅往往是相对值，采用反射率正演模拟的地震数 据振幅与实际振幅存在一定数值差异。以测井数据为输入模型利用反射法正演模拟角度域 的PP波与PS波道集，与实际井旁角度域地震道集对比，计算振幅缩放因子，并应用于所 提取的地震子波，达到模拟记录与实际记录的振幅匹配。

在本发明的其他实施方式中，当地震数据信噪比较高时，为角道集的每一道使用统一 的振幅缩放因子，以保证振幅随偏移距的变化关系；当信噪比低时，可分近、中、远偏移 距分别计算振幅缩放因子，保证模拟记录与实际记录的最佳匹配，减少噪声对反演过程的 影响。

由图1可知，该方法还包括：

S103：采集工区内测井数据的统计模型参数。

S104：根据所述的统计模型参数确定符合该工区的模型参数先验分布函数。即该步骤 是基于工区内所有测井数据建立模型参数的先验信息，包括纵波速度、横波速度和密度的 均值，以及三参数相关的协方差矩阵。图3为步骤S104的具体流程图，由图3可知，该 步骤具体包括：

S301：所述的统计模型参数满足三变量高斯分布。为了降低地震反演的不确定性，提 高反演过程的稳定性，需要从其它途径获得地下介质地震弹性参数模型的信息作为先验信 息。本发明采用三变量高斯分布函数作为先验分布函数。假设统计模型参数服从三变量高 斯分布。

S302：分析所述的测井数据，得到所述统计模型参数的均值；

S303：确定所述统计模型参数的自相关系数以及互相关系数；

S304：根据所述的自相关系数以及互相关系数构建三参数相关的协方差矩阵；

S305：根据所述的三参数相关的协方差矩阵形成符合该工区的模型参数先验分布函 数。

也即，该步骤通过测井数据分析，获得各模型参数的均值，求取各参数的自相关系数 和互相关系数，构建三参数相关的协方差矩阵，形成符合该工区的模型参数先验分布函数。

在后续的反演目标函数中，三变量高斯分布函数相应的规则化表达式为：

$R\left(m\right)=\frac{1}{2}{(m-\mu )}^T{C}_m^{-1}(m-\mu )---\left(1\right)$

其中，m＝[vp,vs,den]^T为参数模型，μ和C_m分别为所统计的模型参数的均值和三参 数协方差矩阵。假设不同深度(时间)点的弹性参数互不相关，可以得到C_m的表达式为：

$C_m\left(\begin{array}{ccccccccc}{\sigma }_{\mathrm{pp}}^1& 0& ...& {\sigma }_{\mathrm{ps}}^1& 0& ...& {\sigma }_{\mathrm{pd}}^1& 0& ...\\ 0& {\sigma }_{\mathrm{pp}}^2& 0& ...& {\sigma }_{\mathrm{ps}}^2& 0& ...& {\sigma }_{\mathrm{pd}}^2& 0\\ & & & & .& & & & \\ & & & & .& & & & \\ & & & & .& & & & \\ {\sigma }_{\mathrm{ps}}^1& 0& ...& {\sigma }_{\mathrm{ss}}^1& 0& ...& {\sigma }_{\mathrm{sd}}^2& 0& ...\\ 0& {\sigma }_{\mathrm{ps}}^2& 0& ...& {\sigma }_{\mathrm{ss}}^2& 0& ...& {\sigma }_{\mathrm{sd}}^2& 0\\ & & & & .& & & & \\ & & & & .& & & & \\ & & & & .& & & & \\ {\sigma }_{\mathrm{pd}}^1& 0& ...& {\sigma }_{\mathrm{sd}}^1& & & {\sigma }_{\mathrm{dd}}^1& & \\ 0& {\sigma }_{\mathrm{pd}}^2& 0& ...& {\sigma }_{\mathrm{sd}}^2& 0& ...& {\sigma }_{\mathrm{dd}}^2& 0\\ & & & & .& & & & \\ & & & & .& & & & \\ & & & & .& & & & \end{array}\right)---\left(2\right)$

其中为第i个深度点纵波速度的自相关系数，为第i个深度点横波速度与密度的 互相关系数，其它类推。这些系数通过对工区内测井数据进行时间延迟统计得到。由于三 变量高斯分布通过协方差矩阵融合纵波速度和横波速度以及密度之间的相关性，降低了三 个属性参数之间的不确定性，因而三变量高斯分布可以有效地降低反演的不确定性。后续 反演过程中需要用到先验分布对模型参数的一阶导数、二阶导数：

$\left(\begin{array}{c}\frac{\partial R\left(m\right)}{\partial m}=C_m^{-1}(m-\mu )\\ \frac{{\partial }^{2}R\left(m\right)}{\partial m^2}=C_m^{-1}\end{array}\right)---\left(3\right)$

由图1可知，该方法还包括：

S105：采集地震构造解释资料。

S106：根据所述的地震构造解释资料以及测井数据建立深度域的初始弹性参数模型。 图4为步骤S106的具体流程图，由图4可知，该步骤具体包括：

S401：根据所述的地震构造解释资料、基于沉积模式建立地质模型；

S402：将所述的测井资料按照构造模式进行插值和外推，得到每条测线的时间域的初 始弹性参数模型，所述的初始弹性参数模型包括纵波速度模型、横波速度模型以及密度模 型；

S403：通过时深转换关系将所述时间域的初始弹性参数模型转换为深度域的初始弹性 参数模型。

建立弹性参数模型主要利用三维空间插值方法，其的技术流程是首先利用散点插值的 方法对各个层位的数据进行插值，完成地质层位建模，然后根据地质层位进行弹性参数横 向插值，即将测井信息进行横向插值，计算得到地下每个点上的弹性参数值，完成初始弹 性参数建模的任务。

由图1可知，该方法还包括：

S107：根据所述深度域的初始弹性参数模型确定PP波与PS波反演残差。图5为步骤 S107的具体流程图，由图5可知，该步骤具体包括：

S501：将所述深度域的初始弹性参数模型作为输入。

S502：在频率域通过传递矩阵算法确定每一个角度的PP波与PS波总反射率矩阵，所 述PP波与PS波总反射率矩阵包含所有层的反射系数与传播时间、层间多次波与透射损失；

S503：对所述PP波与PS波总反射率矩阵进行傅里叶反变换；

S504：将傅里叶反变换后的PP波与PS波总反射率矩阵与地震子波进行褶积，得到角 度域的PP波与PS波角道集；

S505：采集实际角度域地震道集；

S506：将所述角度域的PP波与PS波角道集与实际角度域地震道集作差，得到PP波 与PS波反演残差。如图16所示，为本发明提供的具体实施例中输入的叠前AVO PP波角 道集地震数据示意图，图17为本发明提供的具体实施例中输入的叠前AVO PS波角道集地 震数据示意图，图中纵轴表示时间，单位为毫秒，横轴表示入射角度。

根据Haskell-Dunkin矩阵法，假设一维介质共有N层，首先定义六元素的传递向量：

v_N＝[1 0 0 0 0 0]^T                 (4)

该向量表示半空间介质的初始向量，由第N层的向量可以通过矩阵连续相乘得到第0 层的向量，

v₀＝P₀P₁…P_N-1v_N           (5)

P_i为第i层的传播算子，由下式定义，

$P_i=\left(\begin{array}{ccc}{T}_i^{-1}& \mathrm{if}& i=N-1\\ E_i{T}_i^{-1}& \mathrm{if}& i=0\\ T_i{E}_i{T}_i^{-1}& & \mathrm{Orthers}\end{array}\right)---\left(6\right)$

其中，E_i代表了地震波在第i层内传播时的相位变化，为了模拟得到动校正后的AVO 道集，将各个角度在第i层内的传播时间设定为垂直旅行时：

$E_i=\left(\begin{array}{cccccc}{e}^{(-\mathrm{i\omega }{h}_i(1/v_p+1/v_s))}& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& e^{(-\mathrm{i\omega }{h}_i(1/v_p-1/v_s))}& 0& 0& 0\\ 0& 0& 0& e^{\left(\mathrm{i\omega }{h}_i(1/v_p-1/v_s)\right)}& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& e^{\left(\mathrm{i\omega }{h}_i(1/v_p+1/v_s)\right)}\end{array}\right)---\left(7\right)$

其中T_i代表了第i层界面对总反射率的贡献，具体表达式如下：

$T=\left(\begin{array}{cccccc}{t}_{11}& t_{12}& t_{13}& t_{14}& t_{15}& t_{16}\\ t_{21}& t_{22}& t_{23}& t_{24}& t_{25}& t_{26}\\ t_{31}& t_{32}& t_{33}& t_{34}& t_{35}& t_{36}\\ t_{41}& t_{42}& t_{43}& t_{44}& t_{45}& t_{46}\\ t_{51}& t_{52}& t_{53}& t_{54}& t_{55}& t_{56}\\ t_{61}& t_{62}& t_{63}& t_{64}& t_{65}& t_{66}\end{array}\right)---\left(8\right)$

t₁₁＝-(p²+q_pq_s)/μ＝t₁₆,t₁₂＝-2pq_p/μ,t₁₃＝-(p²-q_pq_s)/μ＝-t₁₄,

t₁₅＝-2pq_s/μ

t₂₁＝iq_s/β²＝-t₂₃＝-t₂₄＝-t₂₆,t₃₁＝-ip(Γ+2q_pq_s)＝t₃₆＝t₄₁＝t₄₆,t₃₂＝-4ip²q_p

t₃₃＝-ip(Γ-2q_pq_s)＝t₄₃＝-t₃₄＝-t₄₄

t₃₅＝-2iΓq_s

t₄₂＝-2iΓq_p

t₄₅＝-4ip²q_s

$t_{51}=-i{q}_p/v_s^2=t_{53}=t_{54}=-t_{56}$

t₆₁＝-μ(Γ²+4p²q_pq_s)＝t₆₆

t₆₂＝-4μΓpq_p

t₆₃＝-μ(Γ²-4p²q_pq_s)＝-t₆₄

t₆₅＝-4μΓpq_s

t₂₂＝t₂₅＝t₅₅＝t₅₂＝0

Γ＝2p²-1/v_s²

矩阵T^-1由矩阵T的各元素重排得到，

$T^{-1}=\left(\begin{array}{cccccc}{t}_{61}& t_{51}& t_{31}& t_{31}& t_{21}& t_{11}\\ -t_{65}& 0& -t_{45}& -t_{35}& 0& -t_{15}\\ -t_{63}& -t_{51}& -t_{33}& -t_{33}& t_{12}& {-t}_{13}\\ t_{63}& -t_{51}& t_{33}& t_{33}& t_{21}& t_{13}\\ -t_{62}& 0& -t_{42}& -t_{32}& 0& -t_{12}\\ t_{61}& -t_{51}& t_{31}& t_{31}& -t_{21}& t_{11}\end{array}\right)---\left(9\right)$

根据公式(5)，从第N层递归计算到第0层得到v₀后，即可求得频率-慢度域内的总 反射率函数，

$\left(\begin{array}{c}{r}_{\mathrm{PP}}=-\frac{v_0\left[4\right]}{v_0\left[1\right]}\\ r_{\mathrm{PS}}=-\frac{v_0\left[5\right]}{v_0\left[1\right]}\end{array}\right)---\left(10\right)$

对频率域的反射率进行傅里叶反变换，再与子波褶积即可获得角度域(或慢度域)的 PP波或PS波叠前AVO道集，

$\left(\begin{array}{c}{G}_{\mathrm{pp}}\left(m\right)=W*r_{\mathrm{pp}}\left(m\right)\\ G_{\mathrm{ps}}\left(m\right)=W*r_{\mathrm{ps}}\left(m\right)\end{array}\right)---\left(11\right)$

该步骤基于深度域的初始弹性参数模型和反射率法正演模拟角度域的多波地震叠前 道集，由正演模拟记录和实际记录直接计算PP波与PS波反演残差，与Zoeppritz方程相 比，该步骤可以很方便地考虑层间多次波与透射损失，并且不需要匹配所模拟的PP波与 PS波旅行时间。

由图1可知，该方法还包括：

S108：构建最大后验概率意义下的反演目标函数。即基于Bayesian原理、先验信息和 反射率法正演算子构建最大后验概率意义下的反演目标函数。图6为步骤S108的具体流 程图，由图6可知，该步骤具体包括：

S601：利用所述模型参数先验分布函数对反演结果PP波与PS波反演残差进行约束；

S602：根据贝叶斯Bayesian原理、综合反演似然函数以及先验分布函数得到后验概率 分布函数。在本发明中，假设地震数据噪声和模型空间满足高斯Gauss分布，则反演似然 函数和先验概率分布满足Gauss分布。且后验概率分布函数满足Gauss分布。

S603：根据所述的后验概率分布函数构建反演目标函数。

由图1可知，该方法还包括：

S109：根据所述的反演目标函数以及PP波与PS波反演残差确定最优弹性参数模型。 图7为步骤S109的具体流程图，由图7可知，该步骤具体包括：

S701：确定所述反演目标函数的梯度以及海森Hessian矩阵。即对目标函数求梯度及 Hessian矩阵，其中梯度的求取需要求解基于反射率法的正演算子对模型参数的一阶导数， 该导数可以通过解析或数值方法求得。为了达到计算精度与效率的平衡，采用主项近似构 建Hessian矩阵。

S702：根据所述的梯度、Hessian矩阵以及PP波与PS波反演残差，采用高斯-牛顿迭 代法确定所述深度域的初始弹性参数模型的更新量。

设模型参数为m^T＝(m₁,m₂,…,m_n)^T，观测地震数据为d^T＝(d₁,d₂,…,d_n)^T，由贝叶斯 理论可以得到，在已知叠前地震数据的情况下，反演地下介质弹性参数的问题可以归结为 求解一个后验概率函数，

$P\left(m\right|d)=\frac{P\left(d\right|m)P\left(m\right)}{P\left(d\right)}---\left(12\right)$

其中P(d)＝∫p(d|m)p(m)dm为归一化因子_，可以看做是常数，P(d|m)为似然函数， P(m)为先验概率分布。假设先验分布和似然函数都满足高斯分布，

$\left(\begin{array}{c}P\left(m\right)=P_{\mathrm{om}}\exp (-\frac{1}{2}R\left(m\right))\\ P\left(d\right|m)={\left({\left(2\pi \right)}^n\det C_D\right)}^{1/2}\exp (-\frac{1}{2}{(d-G\left(m\right))}^T{C}_D^{-1}(d-G\left(m\right)))\end{array}\right)---\left(13\right)$

则多波地震数据后验概率分布可以表示为：

$\left(\begin{array}{c}P\left(m\right|d)\propto P\left(d\right|m)P\left(m\right)\\ \propto {\left({\left(2\pi \right)}^{\mathrm{npp}}\det C_{\mathrm{DPP}}\right)}^{1/2}\exp (-\frac{1}{2}{(d_{\mathrm{PP}}-G_{\mathrm{PP}}\left(m\right))}^T{C}_{\mathrm{DPP}}^{-1}(d_{\mathrm{PP}}-G_{\mathrm{PP}}\left(m\right)))\\ *{\left({\left(2\pi \right)}^{\mathrm{nps}}\det C_{\mathrm{DPS}}\right)}^{1/2}\exp (-\frac{1}{2}{(d_{\mathrm{PS}}-G_{\mathrm{PS}}\left(m\right))}^T{C}_{\mathrm{DPS}}^{-1}(d_{\mathrm{PS}}-G_{\mathrm{PS}}\left(m\right)))\\ *\exp (-\frac{1}{2}R\left(m\right))\end{array}\right)---\left(14\right)$

一般可通过求解上式的最大值作为最优解，等价于求解下式目标函数的最小值所对应 的解，

$\left(\begin{array}{c}S\left(m\right)=\frac{1}{2}{(d_{\mathrm{PP}}-G_{\mathrm{PP}}\left(m\right))}^T{C}_{\mathrm{DPP}}^{-1}(d_{\mathrm{DPP}}-G_{\mathrm{PP}}\left(m\right))\\ +\frac{1}{2}{(d_{\mathrm{PS}}-G_{\mathrm{PS}}\left(m\right))}^T{C}_{\mathrm{DPS}}^T(d_{\mathrm{PS}}-G_{\mathrm{PS}}\left(m\right))+\frac{1}{2}R\left(m\right)\end{array}\right)---\left(15\right)$

可通过高斯-牛顿方法迭代求解该目标函数，首先计算目标函数关于模型参数的一阶偏 导数。若假设地震数据中随机噪声之间相互独立，并服从高斯分布，PP波地震数据的噪音 方差为σ_pp，PS波地震数据的噪音方差为σ_ps，则目标函数中的协方差矩阵和可 以简化对角线为常数的对角阵。目标函数的一阶梯度可表示为，

$\left(\begin{array}{c}\frac{\partial S}{\partial m}={\left(G_{\mathrm{pp}}^r\right)}^T(G_{\mathrm{pp}}\left(m\right)-d_{\mathrm{pp}})+\alpha {\left(G_{\mathrm{ps}}^r\right)}^T(G_{\mathrm{ps}}\left(m\right)-d_{\mathrm{ps}})+\beta \frac{\partial R\left(m\right)}{\partial m}={\gamma }_n\\ G_{\mathrm{pp}}^r=\frac{\delta G_{\mathrm{pp}}^r}{\mathrm{\delta m}},G_{\mathrm{ps}}^r=\frac{\delta G_{\mathrm{ps}}^r}{\mathrm{\delta m}}\end{array}\right)---\left(16\right)$

其中α＝σ_pp/σ_ps控制纵横波数据的使用比重，β＝σ_pp控制先验信息的权重。正演 算子关于模型参数的一阶导数与可通过解析的方法求得，

$\left(\begin{array}{c}{G}_{\mathrm{pp}}^r=W*\frac{\delta r_{\mathrm{pp}}}{\mathrm{\delta m}}\\ G_{\mathrm{ps}}^r=W*\frac{\delta r_{\mathrm{ps}}}{\mathrm{\delta m}}\end{array}\right)---\left(17\right)$

令代表第i个地层的第i₀个弹性参数。i₀＝1,2,3分别代表纵波速度、横波速度和密 度，则反射率对模型参数可表示为，

$\left(\begin{array}{c}\frac{\delta r_{\mathrm{pp}}}{\delta m_i^{i0}}=\frac{\frac{\partial v_0\left[4\right]}{\partial m_i^{i0}}{v}_0\left[1\right]-v_0\left[4\right]\frac{\partial v_0\left[1\right]}{\partial m_i^{i0}}}{{\left(v_0\right[1\left]\right)}^2}\\ \frac{\delta r_{\mathrm{ps}}}{\delta m_i^{i0}}=\frac{\frac{\partial v_0\left[5\right]}{\partial m_i^{i0}}{v}_0\left[1\right]-v_0\left[5\right]\frac{\partial v_0\left[1\right]}{\partial m_i^{i0}}}{{\left(v_0\right[1\left]\right)}^2}\end{array}\right)---\left(18\right)$

根据公式(5)，通过偏微分公式推导，可以得到v₀各个分量关于模型参数的一阶导 数，从而得到正演算子关于模型参数的一阶导数。

目标函数的二阶导数(Hessian矩阵)可近似为，

$\frac{{\partial }^{2}S\left(m\right)}{\partial m^2}={\left(G_{\mathrm{pp}}^r\right)}^T{G}_{\mathrm{pp}}^r+\alpha {\left(G_{\mathrm{ps}}^r\right)}^T{G}_{\mathrm{ps}}^r+\beta \frac{{\partial }^{2}R\left(m\right)}{\partial m^2}=H---\left(19\right)$

最终可以得到模型参数的更新公式：

$m_{n+1}=m_n-{\eta }_n{H}_n^{-1}{\gamma }_n---\left(20\right)$

S703：将所述深度域的初始弹性参数模型的更新量作为输入；

S704：根据所述深度域的初始弹性参数模型的更新量以及地震子波确定更新后的PP 波与PS波反演残差；

S705：根据更新后的PP波与PS波反演残差构建更新后的反演目标函数；

S706：根据更新后的反演目标函数确定更新后的梯度以及Hessian矩阵；

S707：通过所述更新后的PP波与PS波反演残差确定最大迭代次数；

S708：根据所述的最大迭代次数以及更新后的反演目标函数确定最优弹性参数模型， 所述的最优弹性参数模型包括纵波速度、横波速度和密度：

S709：将所述的最优弹性参数模型作为地震三参数叠前AVO反演结果。

也即，本发明根据目标函数的梯度、Hessian矩阵和反演残差，采用高斯-牛顿迭代法 计算模型参数的更新量，将深度域的初始弹性参数模型的更新量作为输入重复迭代以上步 骤S704、S705，并通过反演残差控制最大迭代次数，获得最优模型参数，包括纵波速度、 横波速度和密度。图18为本发明提供的具体实施例中无噪声情况下采用基于反射率法的 多波AVO反演方法得到的纵波速度模型示意图，图19为本发明提供的具体实施例中无噪 声情况下采用基于反射率法的多波AVO反演方法得到的横波速度模型示意图，图20为本 发明提供的具体实施例中无噪声情况下采用基于反射率法的多波AVO反演方法得到的密 度模型示意图。图中纵轴表示时间，单位毫秒，横轴自左至右表示纵波速度(单位：千米 /秒)、横波速度(单位：千米/秒)和密度(单位：克/立方厘米)。

基于反射率法的多波AVO反演方法不仅对纵波速度、横波速度有精确的预测结果， 由于引用了PS波横波进行联合反演，其对密度模型也预测准确。图21为本发明提供的具 体实施例中信噪比为4:1的情况下采用基于反射率法的多波AVO反演方法得到的纵波速度 模块示意图，图22为本发明提供的具体实施例中信噪比为4:1的情况下采用基于反射率法 的多波AVO反演方法得到的横波速度模块示意图，图23为本发明提供的具体实施例中信 噪比为4:1的情况下采用基于反射率法的多波AVO反演方法得到的密度模型示意图，三变 量高斯先验约束对保持反演过程稳定起到了关键作用。

如上所示，即为本发明提供的一种基于反射率法的多波AVO储层弹性参数反演方法， 本发明具体采取以下工作步骤来实现上述技术方案：基于地震数据提取角度依赖的子波， 通过测井数据正演模拟和井旁道地震数据确定振幅缩放因子→基于工区内所有测井数 据统计模型参数的先验信息→利用地震构造解释资料和测井数据，基于沉积模式建立初 始弹性参数模型→利用反射率法正演模拟角度域的多波地震叠前道集并计算反演残差 →计算反演目标函数关于模型参数的梯度和Hessian矩阵→根据高斯-牛顿迭代公式更 新参数模型→重复以上工作步骤直至反演残差达到要求或达到最大迭代次数→输出 最终模型参数，包括纵波速度、横波速度、密度模型，用于油气藏储层预测。

针对常规的AVO反演方法所存在的不稳定性等上述问题，本发明还提供一种基于反 射率法的多波AVO储层弹性参数反演系统。本发明是在研究了以下问题的基础之上提出 的：(1)PP波地震数据对横波速度及密度不敏感，反演存在较大不确定性，PS波数据包 含更多的横波速度和密度信息，能有效提高三参数反演的稳定性；(2)多次波、透射损失 等波的传播效应对地震振幅变化影响较大，在反演过程中考虑这类因素有助于提高反演的 准确度；(3)反射率法能够正演模拟得到水平层状介质情况下地震波场的完全解，从而方 便我们在反演过程中考虑转换波、多次波、透射损失等的影响，特别是对于转换波，无需 考虑PP波与PS波的时间匹配问题；(4)引入地下介质的先验模型参数信息，并提供三参 数之间的相关性，对于降低反演的不确定性，提高反演的精度至关重要。本发明首先利用 贝叶斯理论建立多波AVO非线性反演的理论框架，并通过测井数据分析获得模型参数的 先验分布及三参数之间相关性，以对反演过程进行约束；采用反射率法合成地震记录并与 实际地震数据对比，通过高斯-牛顿迭代方法求解模型参数，实现多波AVO地震叠前反演， 预测得到储层纵横波速度及密度信息。

图8为本发明实施例提供的一种基于反射率法的多波AVO储层弹性参数反演系统的结 构框图，由图8可知，所述的系统包括：

地震道集采集装置101，用于采集工区内的地震叠前道集、测井数据以及实际井旁角 度域地震道集。

地震子波确定装置102，用于基于所述的地震叠前道集、测井数据以及实际井旁角度 域地震道集确定地震子波以及振幅缩放因子。本发明假设反演之前地震子波已经得到，因 而需要基于实际的地震叠前道集和测井数据采取统计方法提取子波，子波在传播过程中受 地层的影响会发生波形或频率变化，提取依赖于入射角度的地震子波能有效提高振幅匹配 程度。图9为地震子波确定装置的结构框图，由图9可知，地震子波确定装置102具体包 括：

地震子波确定模块201，用于根据所述的地震叠前道集以及测井数据、采用统计方法 提取依赖于入射角度的地震子波。在具体的实施方式中，基于实际的地震叠前道集和测井 数据采取统计方法提取若干个依赖于入射角度的地震子波。

输入模型确定模块202，用于将所述的测井数据作为输入模型；

道集模拟模块203，用于采用反射法正演模拟角度域的PP波与PS波道集；

振幅缩放因子确定模块204，用于将所述的PP波与PS波道集与所述的实际井旁角度 域地震道集进行对比，得到振幅缩放因子；

应用模块205，用于将所述的振幅缩放因子应用于所述的地震子波。

在具体的实施方式中，基于测井数据和反射法正演模拟地震角道集并结合实际井旁地 震数据确定振幅缩放因子。实际的地震振幅往往是相对值，采用反射率正演模拟的地震数 据振幅与实际振幅存在一定数值差异。以测井数据为输入模型利用反射法正演模拟角度域 的PP波与PS波道集，与实际井旁角度域地震道集对比，计算振幅缩放因子，并应用于所 提取的地震子波，达到模拟记录与实际记录的振幅匹配。

在本发明的其他实施方式中，当地震数据信噪比较高时，为角道集的每一道使用统一 的振幅缩放因子，以保证振幅随偏移距的变化关系；当信噪比低时，可分近、中、远偏移 距分别计算振幅缩放因子，保证模拟记录与实际记录的最佳匹配，减少噪声对反演过程的 影响。

由图8可知，该系统还包括：

统计模型参数采集装置103，用于采集工区内测井数据的统计模型参数。

分布函数确定装置104，用于根据所述的统计模型参数确定符合该工区的模型参数先 验分布函数。即该步骤是基于工区内所有测井数据建立模型参数的先验信息，包括纵波速 度、横波速度和密度的均值，以及三参数相关的协方差矩阵。图10为分布函数确定装置 的结构框图，由图10可知，分布函数确定装置104具体包括：

高斯分布模块301，所述的统计模型参数满足三变量高斯分布。为了降低地震反演的 不确定性，提高反演过程的稳定性，需要从其它途径获得地下介质地震弹性参数模型的信 息作为先验信息。本发明采用三变量高斯分布函数作为先验分布函数。假设统计模型参数 服从三变量高斯分布。

均值确定模块302，用于分析所述的测井数据，得到所述统计模型参数的均值；

相关系数确定模块303，用于确定所述统计模型参数的自相关系数以及互相关系数；

协方差矩阵构建模块304，用于根据所述的自相关系数以及互相关系数构建三参数相 关的协方差矩阵；

分布函数确定模块305，用于根据所述的三参数相关的协方差矩阵形成符合该工区的 模型参数先验分布函数。

也即，该步骤通过测井数据分析，获得各模型参数的均值，求取各参数的自相关系数 和互相关系数，构建三参数相关的协方差矩阵，形成符合该工区的模型参数先验分布函数。

在后续的反演目标函数中，三变量高斯分布函数相应的规则化表达式为公式(1)。

其中，m＝[vp,vs,den]^T为参数模型，μ和C_m分别为所统计的模型参数的均值和三参 数协方差矩阵。假设不同深度(时间)点的弹性参数互不相关，可以得到C_m的表达式为公式 (2)。

其中为第i个深度点纵波速度的自相关系数，为第i个深度点横波速度与密度的 互相关系数，其它类推。这些系数通过对工区内测井数据进行时间延迟统计得到。由于三 变量高斯分布通过协方差矩阵融合纵波速度和横波速度以及密度之间的相关性，降低了三 个属性参数之间的不确定性，因而三变量高斯分布可以有效地降低反演的不确定性。后续 反演过程中需要用到先验分布对模型参数的一阶导数、二阶导数如公式(3)。

由图8可知，该系统还包括：

解释资料采集装置105，用于采集地震构造解释资料。

弹性参数模型建立装置106，用于根据所述的地震构造解释资料以及测井数据建立深 度域的初始弹性参数模型。图11为弹性参数模型建立装置的结构框图，由图11可知，弹 性参数模型建立装置106包括：

地质模型建立模块401，用于根据所述的地震构造解释资料、基于沉积模式建立地质 模型；

差值外推模块402，用于将所述的测井资料按照构造模式进行插值和外推，得到每条 测线的时间域的初始弹性参数模型，所述的初始弹性参数模型包括纵波速度模型、横波速 度模型以及密度模型；

转换模块403，用于通过时深转换关系将所述时间域的初始弹性参数模型转换为深度 域的初始弹性参数模型。

建立弹性参数模型主要利用三维空间插值方法，其的技术流程是首先利用散点插值的 方法对各个层位的数据进行插值，完成地质层位建模，然后根据地质层位进行弹性参数横 向插值，即将测井信息进行横向插值，计算得到地下每个点上的弹性参数值，完成初始弹 性参数建模的任务。

由图8可知，该系统还包括：

反演残差确定装置107，用于根据所述深度域的初始弹性参数模型确定PP波与PS波 反演残差。图12为反演残差确定装置的结构框图，由图12可知，反演残差确定装置107 具体包括：

输入确定模块501，用于将所述深度域的初始弹性参数模型作为输入。

总反射率矩阵确定模块502，用于在频率域通过传递矩阵算法确定每一个角度的PP波 与PS波总反射率矩阵，所述PP波与PS波总反射率矩阵包含所有层的反射系数与传播时 间、层间多次波与透射损失；

傅里叶反变换模块503，用于对所述PP波与PS波总反射率矩阵进行傅里叶反变换；

褶积模块504，用于将傅里叶反变换后的PP波与PS波总反射率矩阵与地震子波进行 褶积，得到角度域的PP波与PS波角道集；

采集模块505，用于采集实际角度域地震道集；

反演残差确定模块506，用于将所述角度域的PP波与PS波角道集与实际角度域地震 道集作差，得到PP波与PS波反演残差。如图16所示，为本发明提供的具体实施例中输 入的叠前AVO PP波角道集地震数据示意图，图17为本发明提供的具体实施例中输入的叠 前AVO PS波角道集地震数据示意图，图中纵轴表示时间，单位为毫秒，横轴表示入射角 度。

根据Haskell-Dunkin矩阵法，假设一维介质共有N层，首先定义六元素的传递向量如 公式(4)。该向量表示半空间介质的初始向量，由第N层的向量可以通过矩阵连续相乘 得到第0层的向量，如公式(5)。P_i为第i层的传播算子，由公式(6)定义。

其中，E_i代表了地震波在第i层内传播时的相位变化，为了模拟得到动校正后的AVO 道集，将各个角度在第i层内的传播时间设定为垂直旅行时，如公式(7)。其中T_i代表了 第i层界面对总反射率的贡献，具体表达式如公式(8)。

矩阵T^-1由矩阵T的各元素重排得到，如公式(9)。

根据公式(5)，从第N层递归计算到第0层得到v₀后，即可求得频率-慢度域内的总 反射率函数，如公式(10)。

对频率域的反射率进行傅里叶反变换，再与子波褶积即可获得角度域(或慢度域)的 PP波或PS波叠前AVO道集，如公式(11)。

该步骤基于深度域的初始弹性参数模型和反射率法正演模拟角度域的多波地震叠前 道集，由正演模拟记录和实际记录直接计算PP波与PS波反演残差，与Zoeppritz方程相 比，该步骤可以很方便地考虑层间多次波与透射损失，并且不需要匹配所模拟的PP波与 PS波旅行时间。

由图8可知，该系统还包括：

目标函数构建装置108，用于构建最大后验概率意义下的反演目标函数。即基于 Bayesian原理、先验信息和反射率法正演算子构建最大后验概率意义下的反演目标函数。 图13为目标函数构建装置的结构框图，由图13可知，目标函数构建装置108具体包括：

约束模块601，用于利用所述模型参数先验分布函数对反演结果进行约束；

分布函数确定模块602，用于根据Bayesian原理、综合反演似然函数以及先验分布函 数得到后验概率分布函数。在本发明中，假设地震数据噪声和模型空间满足Gauss(高斯) 分布，则反演似然函数和先验概率分布满足Gauss分布。且后验概率分布函数满足Gauss 分布。

目标函数构建模块603，用于根据所述的后验概率分布函数构建反演目标函数。

由图8可知，该系统还包括：

最优弹性参数模型确定装置109，用于根据所述的反演目标函数以及PP波与PS波反 演残差确定最优弹性参数模型。图14为最优弹性参数模型确定装置的结构框图，由图14 可知，该最优弹性参数模型确定装置109具体包括：

梯度确定模块701，用于确定所述反演目标函数的梯度以及Hessian矩阵。即对目标函 数求梯度及Hessian矩阵，其中梯度的求取需要求解基于反射率法的正演算子对模型参数 的一阶导数，该导数可以通过解析或数值方法求得。为了达到计算精度与效率的平衡，采 用主项近似构建Hessian矩阵。

更新量确定模块702，用于根据所述的梯度、Hessian矩阵以及PP波与PS波反演残差， 采用高斯-牛顿迭代法确定所述深度域的初始弹性参数模型的更新量。

设模型参数为m^T＝(m₁,m₂,…,m_n)^T，观测地震数据为d^T＝(d₁,d₂,…,d_n)^T，由贝叶斯 理论可以得到，在已知叠前地震数据的情况下，反演地下介质弹性参数的问题可以归结为 求解一个后验概率函数，如公式(12)。

其中P(d)＝∫p(d|m)p(m)dm为归一化因子，可以看做是常数，P(d|m)为似然函数， P(m)为先验概率分布。假设先验分布和似然函数都满足高斯分布，如公式(13)。则多 波地震数据后验概率分布可以表示为如公式(14)。

一般可通过求解上式的最大值作为最优解，等价于求解下式目标函数的最小值所对应 的解，如公式(15)。

可通过高斯-牛顿方法迭代求解该目标函数，首先计算目标函数关于模型参数的一阶偏 导数。若假设地震数据中随机噪声之间相互独立，并服从高斯分布，PP波地震数据的噪音 方差为σ_pp，PS波地震数据的噪音方差为σ_ps，则目标函数中的协方差矩阵和可 以简化对角线为常数的对角阵。目标函数的一阶梯度可表示为如公式(16)。

其中α＝σ_pp/σ_ps控制纵横波数据的使用比重，β＝σ_pp控制先验信息的权重。正演 算子关于模型参数的一阶导数与可通过解析的方法求得，如公式(17)。

令代表第i个地层的第i₀个弹性参数。i₀＝1,2,3分别代表纵波速度、横波速度和密 度，则反射率对模型参数可表示为如公式(18)。

根据公式(5)，通过偏微分公式推导，可以得到v₀各个分量关于模型参数的一阶导 数，从而得到正演算子关于模型参数的一阶导数。

目标函数的二阶导数(Hessian矩阵)可近似为如公式(19)。最终可以得到模型参数 的更新公式如公式(20)。

更新量输入模块703，用于将所述深度域的初始弹性参数模型的更新量作为输入；

更新残差确定模块704，用于根据所述深度域的初始弹性参数模型的更新量以及地震 子波确定更新后的PP波与PS波反演残差；

更新目标函数构建模块705，用于根据更新后的PP波与PS波反演残差构建更新后的 反演目标函数；

更新梯度确定模块706，用于根据更新后的反演目标函数确定更新后的梯度以及 Hessian矩阵；

最大迭代次数确定模块707，用于通过所述更新后的PP波与PS波反演残差确定最大 迭代次数；

最优弹性参数模型确定模块708，用于根据所述的最大迭代次数以及更新后的反演目 标函数确定最优弹性参数模型，所述的最优弹性参数模型包括纵波速度、横波速度和密度：

反演结果确定模块709，用于将所述的最优弹性参数模型作为地震三参数叠前AVO反 演结果。

也即，本发明根据目标函数的梯度、Hessian矩阵和反演残差，采用高斯-牛顿迭代法 计算模型参数的更新量，将深度域的初始弹性参数模型的更新量作为输入重复迭代以上模 块更新残差确定模块704、更新目标函数构建模块705，并通过反演残差控制最大迭代次 数，获得最优模型参数，包括纵波速度、横波速度和密度。图18为本发明提供的具体实 施例中无噪声情况下采用基于反射率法的多波AVO反演方法得到的纵波速度模型示意图， 图19为本发明提供的具体实施例中无噪声情况下采用基于反射率法的多波AVO反演方法 得到的横波速度模型示意图，图20为本发明提供的具体实施例中无噪声情况下采用基于 反射率法的多波AVO反演方法得到的密度模型示意图。图中纵轴表示时间，单位毫秒， 横轴自左至右表示纵波速度(单位：千米/秒)、横波速度(单位：千米/秒)和密度(单 位：克/立方厘米)。

基于反射率法的多波AVO反演方法不仅对纵波速度、横波速度有精确的预测结果， 由于引用了PS波横波进行联合反演，其对密度模型也预测准确。图21为本发明提供的具 体实施例中信噪比为4:1的情况下采用基于反射率法的多波AVO反演方法得到的纵波速度 模块示意图，图22为本发明提供的具体实施例中信噪比为4:1的情况下采用基于反射率法 的多波AVO反演方法得到的横波速度模块示意图，图23为本发明提供的具体实施例中信 噪比为4:1的情况下采用基于反射率法的多波AVO反演方法得到的密度模型示意图，三变 量高斯先验约束对保持反演过程稳定起到了关键作用。

综上所述，本发明提供了一种基于反射率法的多波AVO储层弹性参数反演方法及系 统，图15为本发明提供的具体实施例中基于反射率法的多波AVO反演方法的流程图，由 图15可知，其主要包括：

1、基于地震数据提取角度依赖的子波；基于测井数据和反射率法正演模拟地震角道 集并结合实际井旁地震数据确定振幅缩放因子；

2、基于工区内所有测井数据统计模型参数的先验信息，包括纵波速度、横波速度和 密度的均值，以及三参数相关的协方差矩阵；

3、利用地震构造解释资料和测井数据，基于沉积模式建立深度域的初始弹性参数模 型；

4、基于深度域的初始弹性参数模型和反射率法正演模拟角度域的多波地震叠前道集， 由正演模拟记录和实际记录直接计算PP波与PS波反演残差；

5、基于Bayesian原理、先验信息和反射率法正演算子构建最大后验概率意义下的反 演目标函数，并计算目标函数的梯度和Hessian矩阵；

6、根据目标函数的梯度、Hessian矩阵和反演残差，采用高斯-牛顿迭代法计算模型参 数的更新量，重复迭代以上步骤(4)、(5)和(6)，并通过反演残差控制最大迭代次 数，获得最优弹性参数模型，包括纵波速度、横波速度和密度。

本发明由于采取以上技术方案，其具有以下优点：

1、本技术方案所述的地震叠前AVO反演基于反射率法，相对于Zoeppritz方程而言， 反射率法能够模拟水平层状介质情况下除直达波和折射波以外的全波场响应，从而使得我 们能够同时利用多种信息进行储层参数预测；

2、本技术方案联合PP波与PS波进行三参数反演，由于PS波地震数据包含更多的横 波速度和密度信息，该方案能有效提高反演的稳定性；

3、本技术方案无需考虑纵横波时间匹配问题；

4、层间多次波干扰地震数据的分辨能力，而且难于在不损伤有效波的前提下去除， 本技术方案能够有效利用层间多次波进行模型参数预测，适用于层间多次波发育的储层；

5、本技术方案考虑了透射损失，使深层反射波的振幅与实际情况相符合，有利于合 成记录与实际记录的匹配；

6、本技术方案通过三变量高斯先验分布引入模型参数之间的相关性，在降低反演不 确定性的同时，提高了反演的精度。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，可以通过计 算机程序来指令相关的硬件来完成，所述的程序可存储于一般计算机可读取存储介质中， 该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、 光盘、只读存储记忆体(Read-Only Memory，ROM)或随机存储记忆体(Random Access  Memory，RAM)等。

本领域技术人员还可以了解到本发明实施例列出的各种功能是通过硬件还是软件来 实现取决于特定的应用和整个系统的设计要求。本领域技术人员可以对于每种特定的应 用，可以使用各种方法实现所述的功能，但这种实现不应被理解为超出本发明实施例保护 的范围。

本发明中应用了具体实施例对本发明的原理及实施方式进行了阐述，以上实施例的说 明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依 据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内 容不应理解为对本发明的限制。