CN102508293B - 一种叠前反演的薄层含油气性识别方法 - Google Patents

Abstract

本发明提供一种叠前反演的薄层含油气性识别方法，包括：1)利用地震构造解释资料和测井数据，基于沉积模式建立初始弹性参数模型；2)基于Bayesian原理的确定性反演得到随机反演的初始模型和模型约束范围；3)采用岩石物理诊断与回归分析建立的三参数相关的Monte Carlo模拟进行随机抽样获得测井尺度模型；4)利用Hashin-Shtrikman边界平均方法进行计算获得地震尺度下弹性参数模型；5)正演模拟角度域叠前地震记录，和实际地震记录计算价值函数，并计算弹性参数模型的接受概率，并按概率确定新的测井尺度弹性参数模型；6)重复迭代以上步骤3)、4)和5)，确定最优的测井尺度弹性参数模型作为地震三参数叠前反演最终结果；7)采用三参数三维空间交汇方法进行油藏薄层含油气性识别。

Description

一种叠前反演的薄层含油气性识别方法

技术领域

本发明涉及油气田地震勘探和油气藏预测与描述方法，特别是涉及一种叠前反演的薄层含油气性识别方法。

背景技术

地震反演是获得地下介质内部图像、对储层进行精细描述的有效方法，也是高分辨率地震勘探的最终表现形式。地震数据反演很大程度上提高了储层表征的价值。随着油气藏勘探程度的增加，在勘探开发程度都比较高的油气区，发现大型的构造油气藏已经是非常困难的事情。而油气藏勘探开发工作者更感兴趣的工作是深化对已勘探、已开发油气藏的认识，对油气藏进行更为精细的描述与表征，进而寻找规模较小的隐蔽型油气藏和发现剩余油分布规律。随着地震研究的重点由勘探逐渐向开发转移，通过地震反演等手段来揭示地下油气藏的精细分布特征，对油气藏进行精雕细刻的技术也受到越来越多的关注。但从目前研究来看，尽管与原始地震道相比，地震反演获得了更高分辨率的成像，但由于地震数据是带限的，反演结果缺乏低频和高频成分，其分辨率无法达到油藏工程师的要求。通过引入正则化方法可以部分地解决这些问题，但即使应用像稀疏约束和块状约束的正则化方法，反演结果也不能达到测井数据的频带宽度。而随机反演综合了先验模型，而先验模型控制了地震数据缺失的低频和高频成分。然而弹性性质的高频变化不会影响地震振幅，因此模型参数高频部分的估计具有高度非唯一性。

总之，目前基于地震叠前反演识别薄层含油气性方法研究存在的主要问题是：1、非均质泥砂岩储层结构变化复杂，泥岩薄层、砂岩薄层厚度超出了基于地震资料进行识别的能力范围，目前的地震反演方法不能实现油藏中单一薄层精细描述，更难以基于地震数据实现单一薄层的定量分析与评价。2、测井数据与地震数据间的测量尺度差异阻碍两者的联合。测井数据、实验室岩芯测量数据具有较高的纵向分辨率，是基于点测量的弹性数据，而地震数据是低频信号，反映了波长内的综合弹性信息。因此，基于测井数据、实验室岩芯测量数据建立的油藏参数与岩石弹性参数变化关系不能直接应用于常规地震波阻抗反演数据解释与反演过程参数约束。3、不同尺度的地球物理数据无法直接进行有效的相互转化，导致不同尺度的地球物理数据不能直接联合获得高分辨率的弹性参数反演结果，从而无法识别和定量分析厚度小于地震分辨率的薄层储层。4、随机反演具有很强的非唯一性且计算量巨大，而降低反演结果多解性的约束方法不够理想，需要进一步提高。

发明内容

本发明实施例提供一种叠前反演的薄层含油气性识别方法，以满足基于地震叠前反演识别薄层含油气性要求。

为了达到上述技术目的，本发明实施例提供了一种叠前反演的薄层含油气性识别方法，所述方法包括：

(1)利用地震构造解释资料和测井数据，基于沉积模式建立初始弹性参数模型；

(2)基于Bayesian原理的确定性反演得到随机反演的初始模型和模型约束范围；

(3)采用岩石物理诊断与回归分析建立的纵波速度、横波速度和密度三参数相关的Monte Carlo模拟进行随机抽样获得测井尺度模型；

(4)利用Hashin-Shtrikman边界平均方法对随机抽样得到测井尺度模型进行计算获得地震尺度下弹性参数模型；

(5)基于地震尺度下弹性参数模型和完整的Zoeppritz非线性方程正演模拟角度域叠前地震记录，由正演模拟角度域叠前地震记录和实际地震记录计算价值函数，并计算弹性参数模型的接受概率，并按概率确定新的测井尺度弹性参数模型；

(6)重复迭代以上步骤(3)、(4)和(5)，并通过所述价值函数和最大迭代次数确定最优的测井尺度弹性参数模型作为地震三参数叠前反演最终结果；

(7)采用纵波速度、横波速度和密度三参数三维空间交汇方法进行油藏薄层含油气性识别。

可选的，在本发明的一实施例中，所述初始弹性参数模型包括：纵波速度模型、横波速度模型和密度模型。

可选的，在本发明的一实施例中，所述基于Bayesian原理的确定性反演得到随机反演的初始模型和模型约束范围，包括：

假设地震数据噪声和模型空间满足高斯分布，则反演似然函数和先验概率分布满足高斯分布，根据Bayesian原理，综合反演似然函数和先验分布函数得到后验概率分布函数，且后验概率分布函数满足高斯分布；

对后验概率分布函数求梯度，并令梯度为零，得到线性方程组；

利用共轭梯度法求解所述线性方程组得到最大后验概率解，即确定性叠前地震反演得到的弹性参数模型。

可选的，在本发明的一实施例中，所述测井尺度模型包括；纵波速度模型、横波速度模型和密度模型。

可选的，在本发明的一实施例中，所述地震尺度下弹性参数模型包括：纵波速度模型、横波速度模型和密度模型。

可选的，在本发明的一实施例中，所述采用纵波速度、横波速度和密度三参数三维空间交汇方法进行油藏薄层含油气性识别，包括：将所述随机反演得到的地震三参数叠前反演最终结果，包括纵波速度、横波速度和密度，进行三维空间交汇分析，并将交汇结果与基于实际测井纵波速度、横波速度、密度数据和含水饱和度建立的含油气性识别图版叠合，根据交汇图叠合结果，确定储层岩性与含油气特征，以实现油藏薄层含油气性识别。

本发明由于采取以上技术方案，其具有以下优点：1、本技术方案考虑了测井尺度模型与地震尺度模型差异，并将测井尺度弹性参数模型到地震尺度弹性参数模型转换应用为地震随机反演，提高模拟数据与实际数据一致性并获得了测井尺度弹性参数模型，有效提高反演分辨率同时，使利用基于测井数据建立交汇图版进行薄层岩性与含油气性识别变得可行。2、本技术方案将岩石物理诊断与回归分析结果应用于约束随机抽样建立弹性参数模型，有效降低了模型抽样随机性，提高计算效果，并具有明确的物理意义。3、本技术方案采用三参数反演与三参数交汇分析方法进行薄层含油气性识别，相对于双参数反演与交汇分析具有更高的识别精度。4、本技术方案采用了确定性反演与随机反演相结合的方法进行叠前地震反演，有效提高了反演精度，降低了计算量。

附图说明

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

图1是本发明实施例叠前反演的薄层含油气性识别方法流程图；

图2是本发明实施例输入的叠前AVO(Amplitude Versus Offset，振幅随偏移距的变化)角道集地震数据；

图3是本发明实施例基于确定性反演得到的纵波速度(左)、横波速度(中)和密度(右)模型；

图4是本发明实施例基于测井数据进行岩石物理诊断结果：纵波速度与模型数据对比(上)，横波速度与模型数据对比(下)；

图5是本发明实施例基于确定性反演结果通过随机抽样得到测井尺度纵波速度(左)、横波速度(中)和密度(右)模型；

图6是本发明实施例基于测井尺度弹性参数模型计算得到的地震尺度纵波速度(左)、横波速度(右)模型；

图7是本发明实施例基于快速模拟退火随机反演得到的测井尺度纵波速度(左)、横波速度(中)和密度(右)模型；

图8是本发明实施例确定性反演结果和随机反演结果纵波速度模型对比，图中左图为0-500ms反演结果，右图为左图中250ms-400m放大的结果；

图9是本发明实施例确定性反演结果和随机反演结果横波速度模型对比，图中左图为0-500ms反演结果，右图为左图中250ms-400m放大的结果；

图10是本发明实施例确定性反演结果和随机反演结果密度模型对比，图中左图为0-500ms反演结果，右图为左图中250ms-400m放大的结果；

图11是本发明实施例基于实际测井数据纵波速度、横波速度、密度和含水饱和度建立的三维空间交汇分析图版，用于薄层岩性识别与含油气性分析，色棒颜色表示储层含水饱和度。

具体实施方式

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

针对常规确定性反演结果无法满足基于地震叠前反演识别薄层含油气性要求等上述问题，本发明的目的是提供一种基于地震三参数叠前反演的薄层含油气性识别方法。本发明是在研究了以下问题基础上提出的：(1)低频地震信号对地下介质测量结果与高频的声波测井对地下介质测量结果存在差异；(2)基于测井声波数据直接合成地震记录与实际地震记录没有良好的对应关系；(3)基于测井尺度测量数据计算地震尺度测量数据后再合成地震记录与实际地震记录有更好的一致性；(4)基于三弹性参数交汇分析识别储层含油气性较双弹性参数交汇分析更有效、精度更高。本发明利用确定性地震反演作为初始模型输入和约束的条件，并在初始模型基础上，采用岩石物理诊断与回归分析建立的三参数相关的Monte Carlo(蒙特卡罗，或称计算机随机模拟方法)模拟进行随机抽样获得测井尺度高分辨率模型，并利用Hashin-Shtrikman边界理论将测井尺度模型转化为地震尺度模型；采用完全Zoeppritz(佐普里兹)方程合成地震记录并与实际地震数据对比，实现叠前地震随机反演，并将反演结果与利用测井数据建立图版叠合，实现薄层含油气性识别。

如图1所示，为本发明实施例一种基于地震三参数叠前反演的薄层含油气性识别方法，具体包括：

101、利用地震构造解释资料和测井数据，基于沉积模式建立初始弹性参数模型：利用地震构造解释资料，基于沉积模式建立地质模型，并将测井资料，按构造模式进行插值和外推，得到每条测线的初始弹性参数模型，包括纵波速度模型、横波速度模型和密度模型，初始弹性参数模型作为确定性反演时的初始输入模型和约束条件。

102、基于Bayesian原理的确定性反演得到随机反演的初始模型和模型约束范围：利用以上所述初始模型和角度域叠前地震道集，基于共轭梯度方法和Bayesian(贝叶斯)原理进行确定性叠前地震反演，得到纵波速度、横波速度与密度确定性反演模型。确定性反演中利用以上所述先验模型对反演结果进行约束，并假设地震数据噪声和模型空间满足Gauss(高斯)分布，则反演似然函数和先验概率分布满足Gauss分布。根据Bayesian原理，综合反演似然函数和先验分布函数得到后验概率分布函数，且后验概率分布函数满足Gauss分布。对后验概率分布函数求梯度，并令梯度为零，得到线性方程组。利用共轭梯度法求解该线性方程组得到最大后验概率解，即确定性叠前地震反演得到的弹性参数模型。

103、采用岩石物理诊断与回归分析建立的纵波速度、横波速度和密度三参数相关的Monte Carlo模拟进行随机抽样获得测井尺度模型：利用以上所述确定性叠前地震反演得到的弹性参数模型作为先验信息和初始解，通过建立蒙特卡洛马尔科夫链(Monte CarloMarkov Chain)对模型空间进行分步随机抽样，即先采用随机抽样方法抽样得到纵波速度的实现，并通过基于测井岩石物理诊断与回归分析建立的纵波速度与横波速度和密度关系进行预测再加一个随机项来生成横波速度和密度模型的实现，从而建立测井尺度弹性参数模型，包括纵波速度模型、横波速度模型和密度模型。与三个参数分别独立进行随机抽样相比本发明考虑了纵波速度与横波速度和密度之间的相关性，也降低了模型空间，提高了计算效率。

104、利用Hashin-Shtrikman边界平均方法对随机抽样得到测井尺度模型进行计算获得地震尺度下弹性参数模型：利用以上所述建立测井尺度弹性参数模型，包括纵波速度模型、横波速度模型和密度模型，通过广义的Hashin-Shtrikman边界理论方法得到地震尺度下弹性参数模型，包括纵波速度模型、横波速度模型和密度模型，实现测井尺度高频模型到地震尺度低频模型转换。

105、基于地震尺度下弹性参数模型和完整的Zoeppritz非线性方程正演模拟角度域叠前地震记录，由正演模拟角度域叠前地震记录和实际地震记录计算价值函数，并计算弹性参数模型的接受概率，并按概率确定新的测井尺度弹性参数模型：利用以上所述建立的地震尺度弹性参数模型，包括纵波速度模型、横波速度模型和密度模型，基于完全的Zoeppritz非线性方程正演模拟角度域叠前地震记录，并利用模拟地震记录和实际地震记录计算价值函数，即通过L2模表征合成角度域叠前地震记录和观测叠前角度域地震记录之间的差异以及反演模型偏离初始模型的程度。从而计算弹性参数模型的接受概率，并按概率确定新的测井尺度弹性参数模型。

106、重复迭代以上步骤103、104和105，并通过所述价值函数和最大迭代次数确定最优的测井尺度弹性参数模型作为地震三参数叠前反演最终结果：利用以上所述确定的新的测井尺度弹性参数模型，重复迭代以上所述步骤103、104和105，并通过以上所述价值函数和最大迭代次数确定最优的测井尺度弹性参数模型作为地震三参数叠前反演最终结果，包括纵波速度、横波速度和密度模型。

107、采用纵波速度、横波速度和密度三参数三维空间交汇方法进行油藏薄层含油气性识别：将以上所述随机反演得到的测井尺度模型最终结果，包括纵波速度、横波速度和密度，进行三维空间交汇分析，并将交汇结果与基于实际测井纵波速度、横波速度、密度数据和含水饱和度建立的含油气性识别图版叠合，根据交汇图叠合结果，确定储层岩性与含油气特征，实现油藏薄层含油气性识别。

本发明具体采取以下工作步骤来实现上述技术方案：利用地震构造解释资料和测井数据，基于沉积模式建立初始弹性参数模型→基于Bayesian原理的确定性反演得到随机反演的初始模型和模型约束范围→采用岩石物理诊断与回归分析建立的三参数相关的MonteCarlo模拟进行随机抽样获得测井尺度高分辨率模型→利用Hashin-Shtrikman边界平均方法对随机抽样得到测井尺度模型进行计算获得地震尺度下弹性参数模型→基于地震尺度下弹性参数模型和完整的Zoeppritz非线性方程正演模拟叠前角度域地震记录，由正演模拟地震记录和实际地震记录计算价值函数，并计算新解的接受概率，并按概率接受新解→重复以上步骤直到满足最大迭代次数或者价值函数满足收敛条件为止，并输出测井尺度下弹性参数反演结果→采用三参数三维空间交汇方法进行油藏薄层含油气性识别。技术方案与工作步骤详细叙述如下：

(1)利用地震构造解释资料，基于沉积模式建立地质模型，并将测井资料，按构造模式进行插值和外推，得到每条测线的初始弹性参数模型，包括纵波速度模型、横波速度模型和密度模型，初始弹性参数模型作为确定性反演时的初始输入模型和约束条件。

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

(2)利用以上所述初始模型和角度域叠前地震道集(如图2所示，是本发明实施例输入的叠前AVO角道集地震数据，图中纵轴表示时间，单位毫秒，横轴表示地震道数)，基于共轭梯度方法和Bayesian原理进行确定性叠前地震反演，得到纵波速度、横波速度与密度确定性反演模型(如图3所示，是本发明实施例基于确定性反演得到的纵波速度(左)、横波速度(中)和密度(右)模型，图中纵轴表示时间，单位毫秒，横轴自左至右表示纵波速度(单位：千米/秒)、横波速度(单位：千米/秒)和密度(单位：克/立方厘米))。确定性反演中以上所述先验模型对反演结果进行约束，并假设地震数据噪声和模型空间满足Gauss分布，则反演似然函数和先验概率分布满足Gauss分布。根据Bayesian原理，综合反演似然函数和先验分布函数得到后验概率分布函数，且后验概率分布函数满足Gauss分布。对后验概率分布函数求梯度，并令梯度为零，得到线性方程组。利用共轭梯度法求解该线性方程组得到最大后验概率解，即确定性叠前地震反演得到的弹性参数模型。具体步骤包括：

考虑两个事件，A和B可以定义它们的事件发生概率分别为P(A)和P(B)，而A和B的联合概率为：

P(A，B)＝P(A|B)P(B)＝P(B|A)P(A)(1)

由公式(1)可以得到Bayesian理论公式：

$P\left(B\right|A)=\frac{P\left(A\right|B)P\left(B\right)}{P\left(A\right)}---\left(2\right)$

对于叠前地震反演来讲，事件B为模型参数m^T＝(m₁，m₂，…，m_m)^T，事件A为叠前观测地震数据d^T＝(D₁，d₂，…，d_M)^T。由Bayesian理论可以得到，在已知叠前地震数据的情况下，反演地下介质弹性参数的反演问题可以归结为求解下式

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

在这里P(d)＝∫p(d|m)p(m)dm为归一化因子，可以作为常数，它保证后验概率分布函数P(m |d)的积分和为1，P(d |m)为似然函数，P(m)为先验概率分布。地质模型与叠前地震数据的关系可以表示为，

d＝G(m)+n    (4)

G(m)为地震正演算子。由于Zoeppritz方程为强非线性方程，不利于直接反演，本发明采用了以相对波阻抗表示的近似方程：

$R_{\mathrm{PP}}\left(\mathrm{\θ}\right)=\frac{1}{2}(1+{\tan }^2\mathrm{\θ})\frac{{\mathrm{\ΔI}}_p}{I_p}-4{\left(\frac{V_s}{V_p}\right)}^2{\sin }^2\mathrm{\θ}\frac{\mathrm{\Δ}{I}_s}{I_s}-[\frac{1}{2}{\tan }^2\mathrm{\θ}-2{\left(\frac{V_s}{V_p}\right)}^2{\sin }^2\mathrm{\θ}]\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}---\left(5\right)$

(5)式中的R_PP(θ)即为本发明中应用于确定性反演的地震正演算子G(m)。假定叠前地震记录中的噪音n满足Gauss分布，那么似然函数满足Gauss分布可以得到：

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

假设先验概率分布函数同样满足Gauss分布即：

$P\left(m\right)={\left({\left(2\mathrm{\π}\right)}^n\det C_M\right)}^{1/2}\exp (-\frac{1}{2}{(m-m_{\mathrm{prior}})}^T{C}_M^{-1}(m-m_{\mathrm{prior}}))---\left(7\right)$

则后验概率分布函数可以表示为：

$P\left(m\right|d)\∝P\left(d\right|m)P\left(m\right)$

$\∝\exp (-\frac{1}{2}{(d-G\left(m\right))}^T{C}_D^{-1}(d-G\left(m\right)))\exp (-\frac{1}{2}{(m-m_{\mathrm{prior}})}^T{C}_M^{-1}(m-m_{\mathrm{prior}}))---\left(8\right)$

m_prior为以上所述建立的初始弹性参数模型，它对最终反演结果有较大的影响。

通过Monte Carlo的方法求解式(8)可以得到模型参数的后验概率分布函数，用于评价反演模型参数的不确定性。但Monte Carlo的反演方法计算量比较大，可以取式(8)的最大值作为最优解，也就是模型参数的最大后验概率解。求解(8)的最大后验概率解也就等价于求解下式的极小值：

$O\left(m\right)=\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))}_{+}\frac{1}{2}{(m-m_{\mathrm{prior}})}^T{C}_M^{-1}\left({m-m}_{\mathrm{prior}}\right)---\left(9\right)$

进一步假设随机噪声之间相互独立，并服从方差为σ_n的高斯分布，则目标函数中的协方差矩阵可以简化为常数。而对于模型参数向量来说，当反射系数之间不相关且服从零均值、方差为σ_m的高斯分布时还可以对目标函数O(m)进行简化，即：

$O\left(m\right)=\frac{1}{2}{(d_{\mathrm{PP}}-G_{\mathrm{PP}}\left(m\right))}^T{(d_{\mathrm{PP}}-G_{\mathrm{PP}}\left(m\right))}_{+}\frac{\mathrm{\μ}}{2}{(m-m_{\mathrm{prior}})}^T\left({m-m}_{\mathrm{prior}}\right)---\left(10\right)$

对目标函数求梯度并令梯度为零，可以将上式的非线性问题转化为线性问题即：

$\▿O\left(m\right)=(G_{\mathrm{PP}}^T{G}_{\mathrm{PP}}+\mathrm{\μI})m-G_{\mathrm{PP}}^T{d}_{\mathrm{PP}}=0---\left(11\right)$

令 $G=G_{\mathrm{PP}}^T{G}_{\mathrm{PP}}+\mathrm{\μI},$ $d=G_{\mathrm{PP}}^T{d}_{\mathrm{PP}},$ (10)式可以简化为，

Gm＝d    (12)

用共轭梯度法求解式(12)即可由叠前地震数据反演得到最大后验概率解，实现确定性反演，反演结果可以作为下一步随机反演的先验信息和初始解。

(3)利用以上所述确定性叠前地震反演得到的弹性参数模型作为先验信息和初始解，通过建立蒙特卡洛马尔科夫链(Monte Carlo Markov Chain)对模型空间进行分步随机抽样，即先采用随机抽样方法抽样得到纵波速度的实现，并通过基于测井岩石物理诊断(如图4所示，是本发明实施例基于测井数据进行岩石物理诊断结果：纵波速度与模型数据(孔隙度)对比(上)，横波速度与模型数据(孔隙度)对比(下)，色棒颜色表示测井曲线伽马值)与回归分析建立的纵波速度与横波速度和密度关系进行预测并加一个随机项来生成横波速度和密度的实现，从而建立测井尺度弹性参数模型，包括纵波速度模型、横波速度模型和密度模型(如图5所示，是本发明实施例基于确定性反演结果通过随机抽样得到测井尺度纵波速度(左)、横波速度(中)和密度(右)模型，图中纵轴表示时间，单位毫秒，横轴自左至右表示纵波速度(单位：千米/秒)、横波速度(单位：千米/秒)和密度(单位：克/立方厘米))。与三个参数分别独立进行随机抽样相比本发明考虑了纵波速度与横波速度和密度之间的相关性，也降低了模型空间，提高了计算效率。

首先通过均匀分布随机抽样uⁱ，且uⁱ满足uⁱ∈U[0，1]。将随机抽样得到的uⁱ代入下式

由公式(13)计算得到yⁱ后，代入下式得，

${\mathrm{\α}}_{k+1}^i={\mathrm{\α}}_k^i+y^i(B_i-A_i)---\left(14\right)$

为第i个变量第k迭代的取值，A_i和B_i分别为纵波速度第i个变量的最小取值和最大取值。A_i和B_i的取值根据以上所述步骤(2)中反演得到的结果分别加减一个弹性参数的最大变化量Δ得到。当时，表明随机抽样成功，否则，重新抽样，从而完成了纵波速度的随机抽样。

由岩石物理诊断和回归分析以及其它先验信息建立纵波速度v_p和横波速度v_s的关系式v_s＝f_s(v_p)，并估算剩余方差σ_s；建立纵波速度α和密度ρ的关系式ρ＝f_ρ(α)，并估算剩余方差σ_ρ。通过关系式和

对横波速度

和密度

进行随机抽样，直到满足

和这样就完成了对横波速度

和密度

的随机抽样，从而形成了测井尺度下的弹性参数模型，包括纵波速度模型、横波速度模型和密度模型，作为下一步地震尺度弹性参数模型计算的输入模型。

(4)利用以上所述建立测井尺度弹性参数模型，包括纵波速度模型、横波速度模型和密度模型，通过广义的Hashin-Shtrikman边界理论方法得到地震尺度下弹性参数模型，包括纵波速度模型、横波速度模型和密度模型(如图6所示，是本发明实施例基于测井尺度弹性参数模型计算得到的地震尺度纵波速度(左)、横波速度(右)模型，图中纵轴表示时间，单位毫秒，横轴自左至右表示纵波速度(单位：千米/秒)、横波速度(单位：千米/秒))，实现测井尺度精确模型到地震尺度低频模型转换。

测井尺度的数据是单点数据分辨率较高，而地震波在地下介质中传播时是一种低频三维体波，因此其包含的弹性参数信息是一定三维空间内岩石弹性参数的综合，而非仅仅一个点，因此地震数据其分比率相对较低。随机抽样得到测井尺度下的纵波速度α_i、横波速度β_i和密度ρ_i后，为了匹配测井尺度模型和地震尺度下的模型，需要对测井尺度下的模型进行计算获得地震尺度下的低频模型。本发明中利用广义Hashin-Shtrikman上下边界平均的方法来对两种尺度的数据进行匹配，即利用测井尺度模型计算地震尺度模型。下面给出广义的Hashin-Shtrikman边界理论方程：

K_low≤K_eff≤K_up    (15)

G_low≤G_eff≤G_up    (16)

上式中K_eff和G_eff分别为等效介质的体积模量和剪切模量；K_low和K_up分别为等效介质体积模量的上下限，且满足

$K_{\mathrm{low}}={\left[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{f_i}{K_i+\frac{4}{3}{G}_{\min }}\right]}^{-1}-\frac{4}{3}{G}_{\min };$ $K_{\mathrm{up}}={\left[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{f_i}{K_i+\frac{4}{3}{G}_{\max }}\right]}^{-1}-\frac{4}{3}{G}_{\max };---\left(17\right)$

ρ_i和α_i分别为半个波长内第i个点的密度和纵波速度。

G_low和G_up分别为等效介质体积模量的上下限，且满足

$G_{\mathrm{low}}={\left[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{f_i}{G_i+\mathrm{\ζ}(K_{\min },G_{\min })}\right]}^{-1}-\mathrm{\ζ}(K_{\min },G_{\min });---\left(18\right)$

$G_{\mathrm{up}}={\left[\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{f_i}{G_i+\mathrm{\ζ}(K_{\max },G_{\max })}\right]}^{-1}-\mathrm{\ζ}(K_{\max },G_{\max });---\left(19\right)$

其中

v_si为一定三维空间非均匀储层内第i个点的横波速度。实际数据计算与分析表明，利用广义Hashin-Shtrikman边界理论计算的地震尺度下弹性参数上、下限非常接近，因此可以通过上、下边界平均得到地震尺度下非均匀岩石储层的等效弹性参数，计算精度较高。

$K_{\mathrm{mix}}=\frac{1}{2}(K_{\mathrm{up}}+K_{\mathrm{low}})---\left(20\right)$

$G_{\mathrm{mix}}=\frac{1}{2}(G_{\mathrm{up}}+G_{\mathrm{low}})---\left(21\right)$

而对于地震尺度下非均匀岩石的密度计算采用体积密度平均方法：

${\mathrm{\ρ}}_i^{\mathrm{new}}=\frac{1}{n}\underset{j=i-n/2}{\overset{i+n/2}{\mathrm{\Σ}}}{\mathrm{\ρ}}_j---\left(22\right)$

上式中ρ_j为一定三维空间非均匀储层内第j个点的密度。进而可以由下式计算得到一定三维空间非均匀储层内纵、横波速度的平均值：

${\mathrm{\α}}_i^{\mathrm{new}}=\sqrt{(K_{\mathrm{mix}}+\frac{4}{3}{G}_{\mathrm{mix}})/{\mathrm{\ρ}}_i^{\mathrm{new}}}---\left(23\right)$

${\mathrm{\β}}_i^{\mathrm{new}}=\sqrt{G_{\mathrm{mix}}/{\mathrm{\ρ}}_i^{\mathrm{new}}}---\left(24\right)$

根据式(15)-(24)即可利用测井尺度下的随机抽样模型，计算得到地震尺度下的弹性参数模型。为基于地震数据反演测井尺度弹性参数奠定基础。

(5)利用以上所述建立的地震尺度弹性参数模型，包括纵波速度模型、横波速度模型和密度模型，基于完全的Zoeppritz非线性方程正演模拟角度域叠前地震记录，并利用模拟地震记录和实际地震记录计算价值函数，即通过L2模表征合成角度域叠前地震记录和观测叠前角度域地震记录之间的差异以及反演模型偏离初始模型的程度。从而计算弹性参数模型的接受概率，并按概率确定新的测井尺度弹性参数模型。

在确定性反演中，正演采用了Zoeppritz方程的线性近似式，但为了更加精确地反演弹性参数，在随机反演中，本发明采用完全的Zoeppritz方程进行叠前地震记录正演模拟。通过完全的Zoeppritz方程计算不同角度的纵波反射系数，再利用褶积模型，将不同角度的反射系数与子波进行褶积运算，得到角度域的叠前地震记录。解析求解Zoeppritz方程可以得到纵波反射系数的解析表达式：

$R_{\mathrm{PP}}\left(\mathrm{\θ}\right)=-1+\frac{2C\cos \mathrm{\θ}}{D}---\left(25\right)$

上式(25)中各参数详细表达式由式(26)-(37)给出

$C=\frac{{\mathrm{\γ}}^2{X}^2{F}_2}{q_2}+\frac{\mathrm{\γ\η}{F}_4}{q_1}+4Q^2{\sin }^2\mathrm{\θ}{F}_2{F}_3{F}_4---\left(26\right)$

$F_2=\sqrt{1-\frac{{\sin }^2\mathrm{\θ}}{q_1^2}}---\left(27\right)$

$F_3=\sqrt{1-{\mathrm{\γ}}^2{\sin }^2\mathrm{\θ}}---\left(28\right)$

$F_4=\sqrt{1-\frac{{\mathrm{\γ}}^2{\sin }^2\mathrm{\θ}}{q_2^2}}---\left(29\right)$

$Q=\frac{{\mathrm{\γ}}^2\mathrm{\η}}{q_2^2}-\frac{1}{q_1^2}---\left(30\right)$

X＝η-2Q sin²θ(31)

Y＝1+2Q sin²θ (32)

Z＝η-1-2Q sin²θ(33)

$T=\frac{{\mathrm{\γ}}^2\mathrm{XZ}}{q_2}+2{\mathrm{QYF}}_3{F}_4---\left(34\right)$

D＝D₁+D₂(35)

$D_1=\frac{{\left(\mathrm{\γ}X\sin \mathrm{\θ}\right)}^2}{q_1{q}_2}+\frac{{\left(\mathrm{\γX}\right)}^2{F}_1{F}_2}{q_2}+\frac{Y^2{F}_3{F}_4}{q_1}---\left(36\right)$

$D_2=\mathrm{\γ\η}(\frac{F_1{F}_4}{q_1}+\frac{F_2{F}_3}{q_2})+{4Q}^2{\sin }^2\mathrm{\θ}{F}_1{F}_2{F}_3{F}_4---\left(37\right)$

基于以上所述的地震尺度下弹性参数模型利用公式(25)计算得到不同角度的反射系数，再与利用测井数据和实际待反演的角度域叠前地震记录提取的子波进行褶积运算生成角度域叠前合成地震记录。计算公式如下，

$d_{\mathrm{syn}}(t,\mathrm{\θ})=w\left(t\right)\⊗R_{\mathrm{pp}}(t,\mathrm{\θ})---\left(38\right)$

本发明采用快速模拟退火算法(Very Fast Simulated Annealing)进行全局寻优，在优化过程中，每迭代一次需要计算一次价值函数，价值函数相当于退火过程中金属的能量。通过L2模表征合成角度域叠前地震记录和观测叠前角度域地震记录之间的差异，并且同时适当考虑反演模型偏离初始模型的程度对结果的影响，反演模型偏离初始模型的程度同样采用L2模来计算。价值函数计算公式如下：

E＝||d_syn-d_obs||₂+α||m-m_lowf||₂(39)

其中d_syn为当前模型的合成角度域叠前地震记录，d_obs为观测角度域叠前地震数据，m为当前模型，m_lowf为利用确定性反演得到的反演结果，α为权值，用于调整随机反演结果对初始模型的依赖程度。当α选取较大的值时，反演结果偏离初始模型较小，反之，当α选取较小的值时，反演结果允许较大地偏离初始模型。

在利用快速模拟退火算法进行全局寻优过程中，新解按概率接受，如果新解使得价值函数降低，即新解的Gibbs概率大于当前解的Gibbs概率，则无条件接受新解；如新解的Gibbs概率小于当前解的Gibbs概率，则按转移概率来接受新解。当前解的Gibbs概率计算方程如下式所示：

$\mathrm{Gibbs}\left(m_i\right)=\frac{\exp \left(\frac{-E\left(m_i\right)}{T_k}\right)}{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}\exp \left(\frac{-E\left(m_j\right)}{T_k}\right)}---\left(40\right)$

如果ΔE＝E(m^new)-E(m^old)≥0，接受新解m^new为当前解。如果ΔE＝E(m^new)-E(m^old)≥0，按如下式所示的转移概率接受新解m^new为当前解：

$h\left(\mathrm{\ΔE}\right)=\frac{\mathrm{Gibbs}\left(m^{\mathrm{new}}\right)}{\mathrm{Gibbs}\left(m^{\mathrm{new}}\right)+\mathrm{Gibbs}\left(m^{\mathrm{old}}\right)}=\frac{1}{1+\exp (\mathrm{\ΔE}/T_k)}\≈\exp \left(\frac{-\mathrm{\ΔE}}{T_k}\right)---\left(41\right)$

首先通过均匀分布随机抽样uⁱ，且uⁱ满足uⁱ∈U[0，1]如果h(ΔE)＞uⁱ则接受新解，否则，放弃新解，保持当前解不变。确定新解即地震尺度下弹性模型后，利用新解对应的测井尺度下弹性参数模型作为初始输入，进行后面的迭代。

(6)利用以上所述确定的新的测井尺度弹性参数模型，重复迭代以上所述步骤(3)、步骤(4)和步骤(5)，并通过以上所述价值函数和最大迭代次数确定最优的测井尺度弹性参数模型作为地震三参数叠前随机反演最终结果，包括纵波速度、横波速度和密度(如图7所示，是本发明实施例基于快速模拟退火随机反演得到的测井尺度纵波速度(左)、横波速度(中)和密度(右)模型，图中纵轴表示时间，单位毫秒，横轴自左至右表示纵波速度(单位：千米/秒)、横波速度(单位：千米/秒)和密度(单位：克/立方厘米))。

利用快速模拟退火算法进行全局寻优过程中，每次迭代需要同时进行降温处理。降温策略是利用快速模拟退火算法进行随机反演的关键步骤之一，降温过快，会导致解很快收敛到局部极值，而降温过快又会导致计算量过大，浪费大量机时。本发明采用以下三种有效的降温策略。

第一种策略，T_i为第i个变量的温度，并且由下式给出，

T_i(k)＝T_0iexp(-c_ik^1/D)(42)

k为迭代次数，T_0i为第i个变量的初始温度，D为变量的维数，c_i为第i个变量的给定常数。本发明对于所有变量给定相同的初始温度和相同的给定常数，即T_0i≡T₀，c_i≡c。

第二种策略，T_i为第i个变量的温度，并且由下式给出，

$T\left(k\right)=\frac{T_0}{\ln (k+1)}---\left(43\right)$

第三种策略，T_i为第i个变量的温度，并且由下式给出，

$T\left(k\right)=\frac{T_0}{k+1}---\left(44\right)$

三种降温策略均是在迭代过程中自适应调整温度，因而算法具有自适应性，避免了人为选取参数的随意性。

重复迭代以上所述步骤(3)、步骤(4)和步骤(5)，直到迭代次数达到预先设定最大值或当满足价值函数E(m^new)＜ε终止迭代，并将所有迭代中价值函数最小时对应的测井尺度弹性参数模型作为最终反演结果，用于下一步弹性参数交汇分析，进行薄层含油气性识别。基于随机反演得到的反演结果较确定性反演结果具有更高的分辨率(纵波速度、横波速度和密度对比结果分布如图8所示，是本发明实施例确定性反演结果和随机反演结果纵波速度模型对比，图中左图为0-500ms反演结果，右图为左图中250ms-400m放大的结果，图中纵轴表示时间，单位毫秒，横轴表示纵波速度(单位：千米/秒)；如图9所示，是本发明实施例确定性反演结果和随机反演结果横波速度模型对比，图中左图为0-500ms反演结果，右图为左图中250ms-400m放大的结果，图中纵轴表示时间，单位毫秒，横轴表示横波速度(单位：千米/秒)；如图10所示，是本发明实施例确定性反演结果和随机反演结果密度模型对比，图中左图为0-500ms反演结果，右图为左图中250ms-400m放大的结果，图中纵轴表示时间，单位毫秒，横轴表示密度(单位：克/立方厘米))，能与测井数据分析结果较好匹配，克服测量尺度对薄层含油气性识别的影响。

(7)将以上所述测井尺度模型反演最终结果，包括纵波速度、横波速度和密度，进行三维空间交汇分析，并将交汇结果与基于实际测井纵波速度、横波速度、密度数据和含水饱和度建立的含油气性识别图版叠合(如图11所示，是本发明实施例基于实际测井数据纵波速度、横波速度、密度和含水饱和度建立的三维空间交汇分析图版，用于薄层岩性识别与含油气性分析，色棒颜色表示储层含水饱和度)，根据交汇图叠合结果，确定储层岩性与含油气特征，实现油藏薄层含油气性识别。

基于实际测井数据进行纵波速度、横波速度和密度进行三维空间三参数交汇分析，并用颜色代表储层含水饱和度，建立基于三参数交汇分析储层含油气性识别图版，在实际测井数据含水饱和度不变化时，可以利用Gassmann方程进行流体替换得到不同含水饱和度时弹性参数。将以上所述测井尺度模型最终反演结果，包括纵波速度、横波速度和密度，同样进行三维空间交汇分析，并将交汇分析结果与以上所述建立图版叠合，根据图版分析储层岩性与含油气性，实现油藏薄层含油气性识别。

本发明由于采取以上技术方案，其具有以下优点：1、本技术方案考虑了测井尺度模型与地震尺度模型差异，并将测井尺度弹性参数模型到地震尺度弹性参数模型转换应用为地震随机反演，提高模拟数据与实际数据一致性并获得了测井尺度弹性参数模型，有效提高反演分辨率同时，使利用基于测井数据建立交汇图版进行薄层岩性与含油气性识别变得可行。2、本技术方案将岩石物理诊断与回归分析结果应用于约束随机抽样建立弹性参数模型，有效降低了模型抽样随机性，提高计算效果，并具有明确的物理意义。3、本技术方案采用三参数反演与三参数交汇分析方法进行薄层含油气性识别，相对于双参数反演与交汇分析具有更高的识别精度。4、本技术方案采用了确定性反演与随机反演相结合的方法进行叠前地震反演，有效提高了反演精度，降低了计算量。

在一个或多个示例性的设计中，本发明实施例所描述的上述功能可以在硬件、软件、固件或这三者的任意组合来实现。如果在软件中实现，这些功能可以存储与电脑可读的媒介上，或以一个或多个指令或代码形式传输于电脑可读的媒介上。电脑可读媒介包括电脑存储媒介和便于使得让电脑程序从一个地方转移到其它地方的通信媒介。存储媒介可以是任何通用或特殊电脑可以接入访问的可用媒体。例如，这样的电脑可读媒体可以包括但不限于RAM、ROM、EEPROM、CD-ROM或其它光盘存储、磁盘存储或其它磁性存储装置，或其它任何可以用于承载或存储以指令或数据结构和其它可被通用或特殊电脑、或通用或特殊处理器读取形式的程序代码的媒介。此外，任何连接都可以被适当地定义为电脑可读媒介，例如，如果软件是从一个网站站点、服务器或其它远程资源通过一个同轴电缆、光纤电脑、双绞线、数字用户线(DSL)或以例如红外、无线和微波等无线方式传输的也被包含在所定义的电脑可读媒介中。所述的碟片(disk)和磁盘(disc)包括压缩磁盘、镭射盘、光盘、DVD、软盘和蓝光光盘，磁盘通常以磁性复制数据，而碟片通常以激光进行光学复制数据。上述的组合也可以包含在电脑可读媒介中。

以上所述的具体实施方式，对本发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施方式而已，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种叠前反演的薄层含油气性识别方法，其特征在于，所述方法包括：

（1）利用地震构造解释资料和测井数据，基于沉积模式建立初始弹性参数模型；

（2）基于Bayesian原理的确定性反演得到随机反演的初始模型和模型约束范围，包括：假设地震数据噪声和模型空间满足高斯分布，则反演似然函数和先验概率分布满足高斯分布，根据Bayesian原理，综合反演似然函数和先验分布函数得到后验概率分布函数，且后验概率分布函数满足高斯分布；对后验概率分布函数求梯度，并令梯度为零，得到线性方程组；利用共轭梯度法求解所述线性方程组得到最大后验概率解，即确定性叠前地震反演得到的弹性参数模型；

（3）采用岩石物理诊断与回归分析建立的纵波速度、横波速度和密度三参数相关的Monte Carlo模拟进行随机抽样获得测井尺度模型；

（4）利用Hashin-Shtrikman边界平均方法对随机抽样得到测井尺度模型进行计算获得地震尺度下弹性参数模型；

（5）基于地震尺度下弹性参数模型和完整的Zoeppritz非线性方程正演模拟角度域叠前地震记录，由正演模拟角度域叠前地震记录和实际地震记录计算价值函数，并计算弹性参数模型的接受概率，并按概率确定新的测井尺度弹性参数模型；

（6）重复迭代以上步骤（3）、（4）和（5），并通过所述价值函数和最大迭代次数确定最优的测井尺度弹性参数模型作为地震三参数叠前反演最终结果；

（7）采用纵波速度、横波速度和密度三参数三维空间交汇方法进行油藏薄层含油气性识别，包括：

将所述随机反演得到的地震三参数叠前反演最终结果，包括纵波速度、横波速度和密度，进行三维空间交汇分析，并将交汇结果与基于实际测井纵波速度、横波速度、密度数据和含水饱和度建立的含油气性识别图版叠合，根据交汇图叠合结果，确定储层岩性与含油气特征，以实现油藏薄层含油气性识别。

2.如权利要求1所述方法，其特征在于，

所述初始弹性参数模型包括：纵波速度模型、横波速度模型和密度模型。

3.如权利要求1所述方法，其特征在于，

所述测井尺度模型包括；纵波速度模型、横波速度模型和密度模型。

4.如权利要求1所述方法，其特征在于，

所述地震尺度下弹性参数模型包括：纵波速度模型、横波速度模型和密度模型。

Images