CN104614763A - 基于反射率法的多波avo储层弹性参数反演方法及系统 - Google Patents

Abstract

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

Description

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

技术领域

本发明关于地球物理勘探技术领域，特别是关于油气田、页岩气地震勘探和储层参数的处理技术，具体的讲是一种基于反射率法的多波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-\mathrm{\μ})}^T{C}_m^{-1}(m-\mathrm{\μ})---\left(1\right)$

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

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

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

$\left\{\begin{array}{c}\frac{\∂R\left(m\right)}{\∂m}=C_m^{-1}(m-\mathrm{\μ})\\ \frac{{\∂}^{2}R\left(m\right)}{\∂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\ω}{h}_i(1/v_p+1/v_s))}& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0\\ 0& 0& e^{(-\mathrm{i\ω}{h}_i(1/v_p-1/v_s))}& 0& 0& 0\\ 0& 0& 0& e^{\left(\mathrm{i\ω}{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\ω}{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₀后，即可求得频率-慢度域内的总反射率函数，

$\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}---\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)为先验概率分布。假设先验分布和似然函数都满足高斯分布，

$\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\mathrm{\π}\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}---\left(13\right)$

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

$\begin{array}{c}P\left(m\right|d)\∝P\left(d\right|m)P\left(m\right)\\ \∝{\left({\left(2\mathrm{\π}\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\mathrm{\π}\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}---\left(14\right)$

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

$\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}---\left(15\right)$

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

$\left\{\begin{array}{c}\frac{\∂S}{\∂m}={\left(G_{\mathrm{pp}}^r\right)}^T(G_{\mathrm{pp}}\left(m\right)-d_{\mathrm{pp}})+\mathrm{\α}{\left(G_{\mathrm{ps}}^r\right)}^T(G_{\mathrm{ps}}\left(m\right)-d_{\mathrm{ps}})+\mathrm{\β}\frac{\∂R\left(m\right)}{\∂m}={\mathrm{\γ}}_n\\ G_{\mathrm{pp}}^r=\frac{\mathrm{\δ}{G}_{\mathrm{pp}}^r}{\mathrm{\δm}},G_{\mathrm{ps}}^r=\frac{\mathrm{\δ}{G}_{\mathrm{ps}}^r}{\mathrm{\δm}}\end{array}\right.---\left(16\right)$

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

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

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

$\begin{array}{c}\frac{\mathrm{\δ}{r}_{\mathrm{pp}}}{\mathrm{\δ}{m}_i^{i0}}=\frac{\frac{\∂v_0\left[4\right]}{\∂m_i^{i0}}{v}_0\left[1\right]-v_0\left[4\right]\frac{\∂v_0\left[1\right]}{\∂m_i^{i0}}}{{\left(v_0\right[1\left]\right)}^2}\\ \frac{\mathrm{\δ}{r}_{\mathrm{ps}}}{\mathrm{\δ}{m}_i^{i0}}=\frac{\frac{\∂v_0\left[5\right]}{\∂m_i^{i0}}{v}_0\left[1\right]-v_0\left[5\right]\frac{\∂v_0\left[1\right]}{\∂m_i^{i0}}}{{\left(v_0\right[1\left]\right)}^2}\end{array}---\left(18\right)$

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

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

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

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

$m_{n+1}=m_n-{\mathrm{\η}}_n{H}_n^{-1}{\mathrm{\γ}}_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 AccessMemory，RAM)等。

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

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

Claims (14)

1.一种基于反射率法的多波AVO储层弹性参数反演方法，其特征是，所述的方法包括：

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

基于所述的地震叠前道集、测井数据以及实际井旁角度域地震道集确定地震子波以及振幅缩放因子；

采集工区内测井数据的统计模型参数；

根据所述的统计模型参数确定符合该工区的模型参数先验分布函数；

采集地震构造解释资料；

根据所述的地震构造解释资料以及测井数据建立深度域的初始弹性参数模型；

根据所述深度域的初始弹性参数模型确定PP波与PS波反演残差；

构建最大后验概率意义下的反演目标函数；

根据所述的反演目标函数以及PP波与PS波反演残差确定最优弹性参数模型。

2.根据权利要求1所述的方法，其特征是，基于所述的地震叠前道集、测井数据以及实际井旁角度域地震道集确定地震子波以及振幅缩放因子包括：

根据所述的地震叠前道集以及测井数据、采用统计方法提取依赖于入射角度的地震子波；

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

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

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

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

3.根据权利要求1或2所述的方法，其特征是，根据所述的统计模型参数确定符合该工区的模型参数先验分布函数包括：

所述的统计模型参数满足三变量高斯分布；

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

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

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

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

4.根据权利要求3所述的方法，其特征是，根据所述的地震构造解释资料以及测井数据建立深度域的初始弹性参数模型包括：

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

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

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

5.根据权利要求1或4所述的方法，其特征是，根据所述深度域的初始弹性参数模型以及地震子波确定PP波与PS波反演残差包括：

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

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

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

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

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

将所述角度域的PP波与PS波角道集与实际角度域地震道集作差，得到PP波与PS波反演残差。

6.根据权利要求5所述的方法，其特征是，构建最大后验概率意义下的反演目标函数包括：

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

根据Bayesian原理、综合反演似然函数以及先验分布函数得到后验概率分布函数；

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

7.根据权利要求6所述的方法，其特征是，根据所述的反演目标函数以及PP波与PS波反演残差确定最优弹性参数模型包括：

确定所述反演目标函数的梯度以及Hessian矩阵；

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

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

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

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

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

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

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

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

8.一种基于反射率法的多波AVO储层弹性参数反演系统，其特征是，所述的系统包括：

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

地震子波确定装置，用于基于所述的地震叠前道集、测井数据以及实际井旁角度域地震道集确定地震子波以及振幅缩放因子；

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

分布函数确定装置，用于根据所述的统计模型参数确定符合该工区的模型参数先验分布函数；

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

弹性参数模型建立装置，用于根据所述的地震构造解释资料以及测井数据建立深度域的初始弹性参数模型；

反演残差确定装置，用于根据所述深度域的初始弹性参数模型确定PP波与PS波反演残差；

目标函数构建装置，用于构建最大后验概率意义下的反演目标函数；

最优弹性参数模型确定装置，用于根据所述的反演目标函数以及PP波与PS波反演残差确定最优弹性参数模型。

9.根据权利要求8所述的系统，其特征是，所述的地震子波确定装置包括：

地震子波确定模块，用于根据所述的地震叠前道集以及测井数据、采用统计方法提取依赖于入射角度的地震子波；

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

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

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

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

10.根据权利要求8或9所述的系统，其特征是，所述的分布函数确定装置包括：

高斯分布模块，所述的统计模型参数满足三变量高斯分布；

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

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

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

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

11.根据权利要求10所述的系统，其特征是，所述的弹性参数模型建立装置包括：

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

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

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

12.根据权利要求8或11所述的系统，其特征是，所述的反演残差确定装置包括：

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

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

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

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

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

反演残差确定模块，用于将所述角度域的PP波与PS波角道集与实际角度域地震道集作差，得到PP波与PS波反演残差。

13.根据权利要求12所述的系统，其特征是，所述的目标函数构建装置包括：

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

分布函数确定模块，用于根据Bayesian原理、综合反演似然函数以及先验分布函数得到后验概率分布函数；

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

14.根据权利要求13所述的系统，其特征是，所述的最优弹性参数模型确定装置包括：

梯度确定模块，用于确定所述反演目标函数的梯度以及Hessian矩阵；

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

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

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

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

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

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

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

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