技术领域
[0001] 本
发明涉及一种油气储层预测过程中的叠前地震资料反演方法，特别是基于流体体积压缩系数的叠前地震反演方法。
背景技术
[0002] 储
层流体识别是油气藏勘探与储层评价的重要环节。在人类油气勘探活动已经经历了两个多世纪的今天，我们面临的是日益复杂的地下条件，隐蔽性强的剩余油气资源分布，识别和描述难度较大的岩性
地层油气藏以及高额的勘探开发成本。与此同时，随着地震采集、处理以及解释技术的不断进步与发展，以地质知识和油气富集规律为指导，基于地震资料蕴含的振幅、
频率等信息实现储层含流体识别已成为可能。以地震资料为主体进行储层流体识别可以在一定程度上增加勘探成功率，提高生产效率以及减少勘探开发成本，是现阶段油气勘探领域的研究热点之一。
[0003] 为了降低勘探开发
风险，提高勘探准确率成为地球物理工作者近年来着重研究的问题，因此，基于地震资料的储层流体识别技术成为研究的热点。Smith和Gidlow(1987)提出可以利用叠前数据通过不同加权函数进行
叠加得到流体因子和伪泊松比剖面来预测岩性和流体，并首次提出了流体因子这个概念，该技术促进了用于流体异常识别的AVO截距-梯度交会技术的发展(Verm和Hilterman，1995；Castagna，1998)。Goodway等人(1997)提出了lambda-mu-rho技术，利用拉梅弹性参数进行储层流体识别。Hilterman(2001)在Goodway和Hedlin等人研究成果的
基础上，对基于AVO的流体识别技术进行了总结。Batzle(2001)
对流体指示因子进行了比较，提出对碎屑岩来说拉梅参数组合属性对流体类型最为敏感，且特别强调在实际应用中要根据区域特性对流体因子进行敏感性选择。George(2003)根据叠前AVO分析，提出了流体因子
角度和交会图角度的概念，通过模型试算和实际应用发现这两种属性对储层流体类型有较强的识别能
力。宁忠华等人(2006)等人在总结分析前人方法的基础上，提出了高灵敏度流体因子的概念。Mark等人(2006)提出了泊松阻抗的概念。在考虑多孔饱和弹性介质的前提下，Russell等人(2003，2006)总结了前人的观点，利用Biot-Gassmann方程对饱和流体条件下的的纵
波速度方程进行了改写，提出用作为流体指示因子，并且指出流体项可以直接作为一项流体因子参与流体检测和储层预测。Hilterman(2009)结合实例应用重点研究了非固结
砂岩的流体因子敏感性，并指出决定非固结砂岩储层流体识别敏感性的关键是纵横波速度的选取。
[0004] 现阶段流体因子的计算方式是基于弹性参数的间接组合运算，此类流体识别的
质量主要取决于两方面，一是基础弹性参数反演是否可靠；二是构建的流体因子对孔隙流体类型是否敏感。现阶段常用流体因子在反演过程中存在反演
精度差的
缺陷。考虑到叠前地震反演是提取基础弹性参数的主要手段，通过改善反演方法可以较好的提高弹性参数可靠性；另一方面，基于间接组合的流体因子不可避免的会造成累计误差，而流体因子直接反演可以较好的解决这个问题，基于双相介质
岩石物理理论构建敏感流体因子，通过研究其与地震动反射特征的内在联系，利用叠前地震资料即可实现流体因子直接反演，从而提高了流体因子的指示敏感性与其估算的可靠性。
发明内容
[0005] 本发明的目的在于，针对
现有技术的不足，提供一种基于流体体积压缩系数的叠前地震反演方法。
[0006] 本发明采用的技术方案如下。
[0007] 基于流体体积压缩系数的叠前地震反演方法，该基于流体体积压缩系数的叠前地震反演方法包括：基于孔隙介质理论，利用理论和经验岩石物理模型，建立流体体积压缩系数与其他弹性参数以及物性参数的联系；在这个关系的基础上，推导基于流体体积压缩系数的反射特征方程；利用反射特征方程建立地震数据与流体体积压缩系数之间的联系，进而实现对流体体积压缩系数的反演。
[0008] 进一步，基于流体体积压缩系数的叠前地震反演方法，其特征在于，包括如下步骤：
[0009] 步骤1：基于孔隙介质理论，利用理论和经验岩石物理模型，建立流体体积压缩系数与其他弹性参数以及物性参数的联系；其他弹性参数包括Gassmann流体因子、
剪切模量、
密度，所述物性参数包括含
水饱和度、孔隙度、泥质含量；
[0010] 步骤2：在平面波假设下，根据步骤1中得到的流体体积压缩系数与其他弹性参数或者物性参数的关系，推导基于流体体积压缩系数的AVO(
地震波振幅随着偏移距的变化)地震反射特征近似方程以及弹性阻抗方程；
[0011] 步骤3，利用步骤2中得到的弹性阻抗方程，建立地震资料-贝叶斯弹性阻抗反演-流体体积压缩系数反演流程，在贝叶斯理论
框架下实现对流体体积压缩系数的直接反演，形成了较为完整、稳定的基于流体体积压缩系数的叠前地震反演方法。
[0012] 进一步，反演完成后，还包括利用反演得到的流体体积压缩系数进行储层流体识别的步骤。
[0013] 进一步，储层流体识别时，还包括基于岩石物理理论以及
测井资料，对步骤1中的流体体积压缩系数进行敏感性分析，验证流体体积压缩系数流体识别精度的步骤。
[0014] 进一步，所述敏感性分析的法方法是：
[0015] 结合测井资料，将流体体积压缩系数与常规的基于单相介质理论的流体因子类型作流体敏感性对比；定义某一流体因子敏感性为：
[0016]
[0017] 其中，meangas,meanwater分别是井上目标层段该流体因子所对应的气层和水层的平均值，stdgas则是标准差。
[0018] 进一步，所述流体因子为纵横波速度、纵横波阻抗、泊松比、
拉梅常数*密度、剪切模量*密度、Gassmann流体因子、泊松比中的一种或数种组合。
[0019] 进一步，在步骤1中，所述孔隙介质理论为双相介质理论。
[0020] 一定
温度下，单位压强的体积相对缩小率
[0021]
[0022] 而双相介质的压缩系数C＝φCf+(1-φ)Cs，其中Cf＝SoCo+SwCw+SgCg，Cg,Co和Cw分别表示气、油和水的去流体体积压缩系数；Sg、So和Sw分别表示气，油和水的饱和度，且Sg+Sw+So＝1。
[0023] Batzle和Han通过研究发现，影响Gassmann流体项f取值的主要因素是孔隙流体体积模量与固体骨架孔隙度；Dehua Han和Batzle等人通过对碎屑岩进行岩石物理统计，研究了Biot-Gassmann理论中的孔隙流体与岩石骨架的固体效应(孔隙度，矿物模量等)对岩石模量信息的耦合作用，提出了突出去岩石骨架流体因子Cf的Gassmann流体项的经验公式：
[0024] f＝G(φ)/Cf (2)。
[0025] 其中， 其中增益函数G(φ)表示岩石骨架矿物与孔隙度的综合作用。
[0026] 结合岩石物理实验得到的经验关系式，利用测井资料可以对流体体积压缩系数Cf进行估算；因此，采用地球物理方法从地震资料中提取出流体体积压缩系数Cf，将流体体积压缩系数Cf作为一项流体因子参与流体识别，即可实现固体骨架与流体弹性效应的解耦，从而有效的提高了储层孔隙流体识别的可靠性。
[0027] 进一步，在步骤2中，利用流体体积压缩系数与Gassmann流体项的关系，推导基于流体体积压缩系数的AVO地震反射特征近似方程以及弹性阻抗方程。
[0028] 包含Gassmann流体项的反射系数近似公式如下所示：
[0029]
[0030] 其中，f,μ和ρ分别表示界面两侧介质的Gassmann流体项，剪切模量和密度的平均值；Δf,Δμ和Δρ则分别表示界面两侧的Gassmann流体项，剪切模量和密度的差值。
[0031] 将公式(2)代入Russell近似公式(4)，进行相应变换，可以得到：
[0032]
[0033] 考虑到剪切模量不受孔隙流体的影响，在此利用干岩石剪切模量μdry替换μ，公式(5)可进一步化简为：
[0034]
[0035] Nur通过大量研究指出，对于小于临界孔隙度的岩石来说，其干岩石的体积模量和剪切模量，可以用与临界孔隙度φc有关的线性函数表示，即临界孔隙度模型，表达式如下式所示：
[0036]
[0037] 其中，φc表示临界孔隙度，Kdry表示干岩石的体积模量，μdry表示干岩石的剪切模量，Km表示固体矿物基质的体积模量，μm表示矿物基质的剪切模量；以Nur模型为纽带对公式(6)进一步展开可以得到：
[0038]
[0039] 将 代入G(φ)，进一步展开可以得到：
[0040]
[0041] 假设 则有Fporo＝φcφμ；公式可以进一步化简为：
[0042]
[0043] 若设定fm＝φμ，则最终得到流体体积压缩系数近似公式如下所示：
[0044]
[0045] 其中，fm＝φμ，称为为固体刚性参数；
[0046] 借鉴Connolly推导弹性阻抗的思想，用弹性阻抗表示反射系数，得到：
[0047]
[0048] 将公式(12)代入(11)，得到：
[0049]
[0050] 将弹性参数的相对变化量用对数形式表示，可以得到：
[0051]
[0053]
[0054] 即
[0055]
[0056] 对上式两边去积分并将其指数化，消掉等式两边的微分项和对数项，进一步去积分常数为0，得到：
[0057] EI(θ)＝(Cf)a(θ)(fm)b(θ)(ρ)c(θ)(φ)d(θ) (17)。
[0058]
[0059] 其中，
[0060] 与常规弹性阻抗公式类似，公式(17)也存在数值量纲随角度变化的问题；在此引入四个参考常数，即A0，Cf0，fm0，ρ0以及φ0；将公式(17)进行标准化处理，可以得到标准化的基于去岩石骨架流体因子的弹性阻抗方程；
[0061] EI(θ)＝A0(Cf/Cf0)a(θ)(fm/fm0)b(θ)(ρ/ρ0)c(θ)(φ/φ0)d(θ) (18)。
[0062] Cf0，fm0，ρ0以及φ0分别定义为Cf0，fm，ρ以及φ的平均值；A0为标准化因子，具体表达式为：
[0063]
[0064] 进一步，在步骤3中，利用贝叶斯反演框架，将反映地震信息的似然函数与反映待反演参数的先验地质约束相结合，通过求解最大后验概率密度函数的方式建立反演目标函数，常规贝叶斯AVO反演对正演方程进行去相关处理，而没有在先验约束中考虑参数之间的相关特性，为了进一步提高反演质量，在对正演方程进行去相关处理的基础上，采用四变量柯西分布作为先验正则约束对贝叶斯AVO四参数反演进行改进，模型试算与实际应用表明该方法在一定程度上提高了四参数反演的可靠性。
[0065] 将公式(11)按照入射角度的不同表示为矩阵形式为：
[0066]
[0067] 其中，ai(i＝1,2…,m)，bi(i＝1,2…,m)，ci(i＝1,2…,m)和di(i＝1,2…,m)分别表示第i个入射角度的相应系数；将其推广到具有m个入射角度，n个界面的情况，并且将矩阵进行
块化处理，可以得到：
[0068]
[0069] 其中，Ri(i＝1,2…,m)表示第i个入射角度的反射系数向量，由n个元素组成；Ai(i＝1,2…,m)，Bi(i＝1,2…,m)，Ci(i＝1,2…,m)和Di(i＝1,2…,m)分别表示第i个入射角度对应的正演系数矩阵，分别是n×n维的斜对角矩阵； Rρ和Rφ则分别表示去岩石骨架流体因子、固体刚性参数、密度以及孔隙度相对变化率向量，分别由n个元素组成。
[0070] 基于地震记录符合褶积模型的假设，引入子波矩阵W，则公式(21)进一步变为：
[0071]
[0072] 其中，di(i＝1,2…,m)表示为第i个入射角度的地震数据组成的列向量，都包含n个元素；
[0073] 依据测井数据的样本统计的方法生成协方差矩阵，即：
[0074]
[0075] 其中， N表示样本的个数；
[0076] 经过去相关变换之后的待反演参数的协方差矩阵Cx的非对角线元素为零，这说明变换之后的参数变为相互独立的，有利于提高参数反演的可靠性。
[0077] 基于贝叶斯反演框架，通过求解最大后验概率密度函数构建反演目标函数，具体到该叠前AVO反演问题，后验概率密度函数可以表示为：
[0078] P(R|d,I)＝const0×P(d|R,I)P(R|I) (24)。
[0079] 其中，P(d|R,I)为似然函数，P(R|I)为先验分布函数，d表示随入射角度变化的叠前地震数据，I表示基本的地质信息，R表示待反演的模型参数，const0是概率归一化常数；
[0080] 由于基于贝叶斯理论的反演思想最终只是关心后验概率密度函数的形状，因此const0可以被忽略，则公式(24)可以进一步化简为：
[0081] P(R|d,I)＝P(d|R,I)P(R|I) (25)。
[0082] 最大后验概率密度函数的数值即两个随机函数的乘积最大值，而后验分布的宽度则对待估计参数的不确定性进行了表征；似然函数主要用于表示观测地震数据与待反演参数之间的关系，在此借助正演方程得到正演记录，通过研究观测数据与正演记录之间误差(即噪声)的特征来构建似然函数；假设地震噪声服从高斯分布，且不同的测量条件的噪声之间满足相互独立条件，将似然函数表示为：
[0083]
[0085] 先验函数主要表示待反演参数的统计特征，较为常见的先验分布有高斯分布和柯西分布；通常情况下使用的柯西先验约束是假设不同界面不同弹性参数分布特征相同，且相互独立，如对正演方程进行去相关处理时所论述的一样，实际情况中纵横波阻抗相对变化率与密度参数相对变化率之间存在一定的统计相关特性，因此，采用单变量的柯西先验正则约束，而忽略参数的相关特性会影响弹性参数的反演质量；在此假设不同界面参数分布符合独立特性，采用四变量柯西分布描述纵横波阻抗相对变化率与密度相对变化率的分布特征，从而充分考虑了四参数之间的相关特性。
[0086] 先验函数可表示如下：
[0087]
[0088] 其中， Cx是协方差矩阵，Di是4×4n的矩阵，其组成元素的取值定义如下：
[0089]
[0090] 因此，将似然函数与先验函数代入后验概率密度函数，可以得到：
[0091]
[0092] 省略常数，求解最大后验概率，可以得到反演目标函数，具体形式如下：
[0093] (G'TG'+2Q)R'＝G'Td (30)。
[0094] 其中， 公式左边第一项主要用来约束正演记录与实际地震记录之间的相近程度，第二项则是四变量柯西正则约束项，主要用来约束反演参数的稀疏程度，然后采用
迭代重加权最小二乘
算法对反演方程进行目标寻优。
[0095] 流体因子的计算方式是基于弹性参数的间接组合运算，此类流体识别的质量主要取决于两方面，一是基础弹性参数反演是否可靠；二是构建的流体因子对孔隙流体类型是否敏感。考虑到叠前地震反演是提取基础弹性参数的主要手段，通过改善反演方法可以较好的提高弹性参数可靠性；另一方面，基于间接组合的流体因子不可避免的会造成累计误差，而流体因子直接反演可以较好的解决这个问题，基于双相介质岩石物理理论构建敏感流体因子，通过研究其与地震动反射特征的内在联系，利用叠前地震资料即可实现流体因子直接反演，从而提高了流体因子的指示敏感性与其估算的可靠性。气体、液体以及固体都具有不同程度的压缩性，压缩系数定义为在一定温度下，压强增加一个单位体积的相对缩小率，可以定量地表征流体与固体的压缩性，一般情形下，气体、液体、固体的压缩系数存在如下关系：气体>液体>固体。本发明首先从构建孔隙流体敏感参数出发，考虑孔隙流体(气体、液体)与岩石骨架的压缩系数之间的差异，将流体体积压缩系数作为流体因子，进而研究流体体积压缩系数的反演方法，具有反演精度高、储层流体识别效果好的优点，对提高储层流体识别精度具有实际意义。
附图说明
[0096] 图1为拉梅常数*密度随孔隙度和含水饱和度的变化趋势图。
[0097] 图2为图1中虚线范围内(含水饱和度为10％-40％)拉梅常数*密度随孔隙度和含水饱和度的变化示意图。
[0098] 图3为流体体积压缩系数Cf随孔隙度和含水饱和度的变化趋势图。
[0099] 图4为图3中虚线范围内(含水饱和度为10％-40％)流体体积压缩系数Cf随孔隙度和含水饱和度的变化示意图。
[0100] 图5为不同流体因子敏感性分析图；其中，1-纵波速度；2-纵波阻抗；3-拉梅参数*密度；4-剪切模量*密度；5-流体因子；6-泊松比；7-横波阻抗；8:-横波速度；9-流体体积压缩系数。
[0101] 图6为模型一的反射系数计算结果.
[0102] 图7为模型二的反射系数计算结果。
[0103] 图8为本发明的所用的叠后地震剖面。
[0104] 图9是与图7对应的流体体积压缩系数反演结果剖面。
[0105] 图10为本发明的基于流体体积压缩系数的叠前地震反演方法详细技术路线图。
具体实施方式
[0107] 下面，结合附图和
实施例对本发明作进一步说明。
[0108] 实施例1。基于流体体积压缩系数的叠前地震反演方法，其特征在于，该基于流体体积压缩系数的叠前地震反演方法包括：
[0109] 步骤1，基于孔隙介质理论，利用理论和经验岩石物理模型，建立流体体积压缩系数与其他弹性参数以及物性参数的联系。
[0110] 步骤2，推导基于流体体积压缩系数的AVO近似方程以及弹性阻抗方程。
[0111] 步骤3，在贝叶斯理论框架下实现对流体体积压缩系数的直接反演，形成了较为完整、稳定的基于流体体积压缩系数的叠前地震反演方法。
[0112] 步骤4，将方法运用于实际资料。
[0113] 在步骤1中，该技术首先需要根据孔隙介质理论，利用理论和经验岩石物理模型，建立流体体积压缩系数与其他弹性参数以及物性参数的联系。
[0114] 一定温度下，单位压强的体积相对缩小率
[0115]
[0116] 双相介质理论充分考虑了介质的岩石骨架结构和孔隙流体性质以及局部特性与整体效应的关系，将含流体储层表述为固体相和流体相的
复合体，且分别考虑了固体和流体以及二者相互耦合对地震波传播的影响。而双相介质的压缩系数C＝φCf+(1-φ)Cs，其中Cf＝SoCo+SwCw+SgCg，Cg,Co和Cw分别表示气、油和水的去流体体积压缩系数；Sg、So和Sw分别表示气，油和水的饱和度，且Sg+Sw+So＝1。
[0117] Batzle和Han通过研究发现，影响Gassmann流体项f取值的主要因素是孔隙流体体积模量与固体骨架孔隙度。Dehua Han和Batzle等人通过对碎屑岩进行岩石物理统计，研究了Biot-Gassmann理论中的孔隙流体与岩石骨架的固体效应(孔隙度，矿物模量等)对岩石模量信息的耦合作用，提出了突出去岩石骨架流体因子Cf的Gassmann流体项的经验公式：
[0118] f＝G(φ)/Cf (2)。
[0119] 其中， 其中增益函数G(φ)表示岩石骨架矿物与孔隙度的综合作用。
[0120] 结合岩石物理实验得到的经验关系式，利用测井资料可以对流体体积压缩系数Cf进行估算。图1为不同孔隙度与含水饱和度情况下的Cf的变化趋势，从图中我们可以看到流体体积压缩系数Cf与含水饱和度成完全线性变化趋势，且完全不受孔隙度的影响。而且随着含气饱和度的增加，流体体积压缩系数Cf变化特别明显，表明流体体积压缩系数Cf对于流体敏感性很强，有利于提高流体识别的精度。因此，如果我们采用特定的地球物理方法从地震资料中提取出流体体积压缩系数Cf，将流体体积压缩系数Cf作为一项流体因子参与流体识别，即可实现固体骨架与流体弹性效应的解耦，从而有效的提高了储层孔隙流体识别的可靠性。
[0121] 结合测井资料，将去岩石骨架流体因子与常规的基于单相介质理论的流体因子类型(纵横波阻抗、泊松比，拉梅参数以及λρ等)和Gassmann流体项作流体敏感性对比。定义某一流体因子敏感性为：
[0122]
[0123] 其中，meangas,meanwater分别是井上目标层段该流体因子所对应的气层和水层的平均值，stdgas则是标准差。结合某工区的井资料数据，对常规9种流体因子做敏感性比较分析，得到图5,可以看到，流体体积压缩系数对储层流体具有最高的敏感性。在利用测井资料分析敏感性的时候不需要使用AVO近似方程，一般的做法是在得到敏感性最高的流体因子后再使用(或者推导)对应的AVO近似方程进行流体识别等后续工作。
[0124] 岩石的敏感性分析并不是当流体体积压缩系数被反演出来后再被使用的，是在之前使用的。岩石的敏感性分析指导开展流体体积压缩系数的构建，旨在为所构建的流体体积压缩系数进行实际可行性的验证。图5中所写的流体体积压缩系数指的是被反演出的流体体积压缩系数(在时间角度上与岩石敏感性分析时使用的流体体积压缩系数不一样)，岩石敏感性分析可以看到流体体积压缩系数的敏感性高，因此被反演出的流体体积压缩系数可以作为流体识别的指示因子。
[0125] 在步骤2中，利用流体体积压缩系数与Gassmann流体项的关系，推导基于流体体积压缩系数的AVO地震反射特征近似方程以及弹性阻抗方程。
[0126] 基于Biot-Gassmann理论，Russell等人对饱含流体多孔介质的AVO理论进行了研究，在2006年的研究报告中首次提出了包含Gassmann流体项的反射系数近似公式，并且在2011年在《Geophysics》正式发表论文发对其进行了讨论，Russell近似公式如下所示：
[0127]
[0128] 其中，f,μ和ρ分别表示界面两侧介质的Gassmann流体项，剪切模量和密度的平均值；Δf,Δμ和Δρ则分别表示界面两侧的Gassmann流体项，剪切模量和密度的差值。Russell等人通过模型试算表明该近似公式在入射角度小于50度的范围内满足近似精度要求。
[0129] 将公式(2)代入Russell近似公式(4)，进行相应变换，可以得到：
[0130]
[0131] 考虑到剪切模量不受孔隙流体的影响，在此利用干岩石剪切模量μdry替换μ，公式(5)可进一步化简为：
[0132]
[0133] Nur通过大量研究指出，对于小于临界孔隙度的岩石来说，其干岩石的体积模量和剪切模量，可以用与临界孔隙度φc有关的线性函数表示，即临界孔隙度模型，表达式如下式所示：
[0134]
[0135] 其中，φc表示临界孔隙度，Kdry表示干岩石的体积模量，μdry表示干岩石的剪切模量，Km表示固体矿物基质的体积模量，μm表示矿物基质的剪切模量。以Nur模型为纽带对公式(6)进一步展开可以得到：
[0136]
[0137] 将 代入G(φ)，进一步展开可以得到：
[0138]
[0139] 假设 则有Fporo＝φcφμ。公式可以进一步化简为：
[0140]
[0141] 若设定fm＝φμ，则最终得到流体体积压缩系数近似公式如下所示：
[0142]
[0143] 其中，fm＝φμ，称为为固体刚性参数。为了验证该公式精度，在此结合某碎屑岩工区的实际数据，利用流体替代的方法设计了二种两层砂岩模型对其进行测试。模型参数如表2-1所示，其中临界孔隙度φc均为27％。模型一为上层砂岩含气，下层砂岩含水，两层孔隙度相同；模型二的两层砂岩所含流体与孔隙度均不同，上层砂岩含气，孔隙度为20％，下层砂岩含水，孔隙度为25％。分别采用精确的Zoeppritz方程、Aki-Richard近似公式以及去岩石骨架AVO近似公式计算两种模型界面的反射系数，两个模型的计算结果分别如错误！未找到引用源。6、图7所示。
[0144] 表1两种砂岩模型参数
[0145]
[0146] 从图6、图7可以发现，当界面两侧孔隙度没有变化，仅孔隙流体不同时，去岩石骨架AVO近似的精度与Aki-Richard近似吻合较好，且均能够与Zoeppritz精确结果有较好的匹配；从图3可以看到，当界面两侧孔隙度和孔隙流体都存在差异的时候，去岩石骨架AVO近似仍然与Aki-Richard近似几乎一致，但两者同样随着入射角度的增加偏离精确Zoeppritz方程计算结果。分析其原因是，去岩石骨架AVO近似是从Aki-Richard近似推导而来，因此不难理解两者近似精度的一致性。此外，两者会在孔隙度变化的情况下产生较大误差的主要原因是：AVO线性近似公式的建立是基于界面两侧具有微小弹性扰动的假设前提，而孔隙度的变化对岩石模量信息的影响较大，进而造成了大角度入射情况下的较大误差。但是，基于分析结果我们可以看到，去岩石骨架AVO近似在小角度入射情况下是可以满足反射系数的近似精度的，考虑到实际应用中使用的角度道集一般不会超过30度，因此，可以认为在误差允许的范围内去岩石骨架AVO近似是完全满足反射系数精度要求的。
[0147] 借鉴Connolly推导弹性阻抗的思想，用弹性阻抗表示反射系数，得到：
[0148]
[0149] 将公式(12)代入(11)，得到：
[0150]
[0151] 将弹性参数的相对变化量用对数形式表示，可以得到：
[0152]
[0153] 进一步变形，可以得到：
[0154]
[0155] 即
[0156]
[0157] 对上式两边去积分并将其指数化，消掉等式两边的微分项和对数项，进一步去积分常数为0，得到：
[0158] EI(θ)＝(Cf)a(θ)(fm)b(θ)(ρ)c(θ)(φ)d(θ) (17)。
[0159]
[0160] 其中，
[0161] 与常规弹性阻抗公式类似，公式(17)也存在数值量纲随角度变化的问题。在此引入四个参考常数，即A0，Cf0，fm0，ρ0以及φ0。将公式(17)进行标准化处理，可以得到标准化的基于去岩石骨架流体因子的弹性阻抗方程。
[0162] EI(θ)＝A0(Cf/Cf0)a(θ)(fm/fm0)b(θ)(ρ/ρ0)c(θ)(φ/φ0)d(θ) (18)。
[0163] Cf0，fm0，ρ0以及φ0分别定义为Cf0，fm，ρ以及φ的平均值。A0为标准化因子，具体表达式为：
[0164]
[0165] 在步骤3中，利用贝叶斯反演框架，将反映地震信息的似然函数与反映待反演参数的先验地质约束相结合，通过求解最大后验概率密度函数的方式建立反演目标函数，常规贝叶斯AVO反演对正演方程进行去相关处理，而没有在先验约束中考虑参数之间的相关特性，为了进一步提高反演质量，本实施例在对正演方程进行去相关处理的基础上，采用四变量柯西分布作为先验正则约束对贝叶斯AVO四参数反演进行改进，模型试算与实际应用表明该方法在一定程度上提高了四参数反演的可靠性。贝叶斯反演是一种反演算法，图10所提到的黑箭头下面的框指的是贝叶斯反演需要输入的一些数据。
[0166] 反射系数线性近似公式是建立弹性参数与地震数据之间关系的
桥梁，不同的反射系数线性近似公式所包含的弹性参数种类以及信息量存在差异。本次研究采用去岩石骨架AVO反射系数线性近似公式建立正演模型。将公式(11)按照入射角度的不同表示为矩阵形式为：
[0167]
[0168] 其中，ai(i＝1,2…,m)，bi(i＝1,2…,m)，ci(i＝1,2…,m)和di(i＝1,2…,m)分别表示第i个入射角度的相应系数。将其推广到具有m个入射角度，n个界面的情况，并且将矩阵进行块化处理，可以得到：
[0169]
[0170] 其中，Ri(i＝1,2…,m)表示第i个入射角度的反射系数向量，由n个元素组成；Ai(i＝1,2…,m)，Bi(i＝1,2…,m)，Ci(i＝1,2…,m)和Di(i＝1,2…,m)分别表示第i个入射角度对应的正演系数矩阵，分别是n×n维的斜对角矩阵； Rρ和Rφ则分别表示去岩石骨架流体因子、固体刚性参数、密度以及孔隙度相对变化率向量，分别由n个元素组成。
[0171] 基于地震记录符合褶积模型的假设，引入子波矩阵W，则公式(21)进一步变为：
[0172]
[0173] 其中，di(i＝1,2…,m)表示为第i个入射角度的地震数据组成的列向量，都包含n个元素。
[0174] 为了提高反演问题的求解
稳定性，一些学者在对正演问题进行数学表征的过程中引入协方差矩阵对待反演参数进行去相关处理。依据测井数据的样本统计的方法生成协方差矩阵，即：
[0175]
[0176] 其中， N表示样本的个数。
[0177] 经过去相关变换之后的待反演参数的协方差矩阵Cx的非对角线元素为零，这说明变换之后的参数变为相互独立的，有利于提高参数反演的可靠性。
[0178] 基于贝叶斯反演框架，通过求解最大后验概率密度函数构建反演目标函数，具体到该叠前AVO反演问题，后验概率密度函数可以表示为：
[0179] P(R|d,I)＝const0×P(d|R,I)P(R|I) (24)。
[0180] 其中，P(d|R,I)为似然函数，P(R|I)为先验分布函数，d表示随入射角度变化的叠前地震数据，I表示基本的地质信息，R表示待反演的模型参数，const0是概率归一化常数。
[0181] 由于基于贝叶斯理论的反演思想最终只是关心后验概率密度函数的形状，因此const0可以被忽略，则公式(24)可以进一步化简为：
[0182] P(R|d,I)＝P(d|R,I)P(R|I) (25)。
[0183] 最大后验概率密度函数的数值即两个随机函数的乘积最大值，而后验分布的宽度则对待估计参数的不确定性进行了表征。似然函数主要用于表示观测地震数据与待反演参数之间的关系，在此借助正演方程得到正演记录，通过研究观测数据与正演记录之间误差(即噪声)的特征来构建似然函数。假设地震噪声服从高斯分布，且不同的测量条件的噪声之间满足相互独立条件，将似然函数表示为：
[0184]
[0185] 其中，σm表示噪声信号的标准方差。
[0186] 先验函数主要表示待反演参数的统计特征，较为常见的先验分布有高斯分布和柯西分布。通常情况下使用的柯西先验约束是假设不同界面不同弹性参数分布特征相同，且相互独立，如对正演方程进行去相关处理时所论述的一样，实际情况中纵横波阻抗相对变化率与密度参数相对变化率之间存在一定的统计相关特性，因此，采用单变量的柯西先验正则约束，而忽略参数的相关特性会影响弹性参数的反演质量。在此假设不同界面参数分布符合独立特性，采用四变量柯西分布描述纵横波阻抗相对变化率与密度相对变化率的分布特征，从而充分考虑了四参数之间的相关特性。
[0187] 先验函数可表示如下：
[0188]
[0189] 其中， Cx是协方差矩阵，Di是4×4n的矩阵，其组成元素的取值定义如下：
[0190]
[0191] 因此，将似然函数与先验函数代入后验概率密度函数，可以得到：
[0192]
[0193] 省略常数，求解最大后验概率，可以得到反演目标函数，具体形式如下：
[0194] (G'TG'+2Q)R'＝G'Td (30)。
[0195] 其中， 公式左边第一项主要用来约束正演记录与实际地震记录之间的相近程度，第二项则是四变量柯西正则约束项，主要用来约束反演参数的稀疏程度，然后采用迭代重加权最小二乘算法(Iterative Reweighted Least Squares,IRLS)对反演方程进行目标寻优。
[0196] 气体、液体以及固体都具有不同程度的压缩性，压缩系数定义为在一定温度下，压强增加一个单位体积的相对缩小率，可以定量地表征流体与固体的压缩性，一般情形下，气体、液体、固体的压缩系数存在如下关系：气体液体固体。本实施例首先从构建孔隙流体敏感参数出发，考虑孔隙流体(气体、液体)与岩石骨架的压缩系数之间的差异，将流体体积压缩系数作为流体因子，进而研究流体体积压缩系数的反演方法，对提高储层流体识别精度具有实际意义。
[0197] 本实施例首先基于岩石物理理论以及测井资料，对流体体积压缩系数进行敏感性分析，用于验证流体体积压缩系数流体识别精度。然后基于孔隙弹性理论，利用理论和经验岩石物理模型，充分考虑岩石模量、孔隙流体及孔隙大小对储层流体影响，推导基于流体体积压缩系数的平面波反射特征近似方程，并分析其适用条件和敏感性”对应图10中的从左边到右边的“双相介质理论”→“流体体积压缩系数”→“AVO近似方程”→“弹性阻抗方程”，双相介质理论是孔隙介质介质理论中的一种，平面波反射特征近似方程就是推导的AVO近似方程，因为目前使用的近似都是基于平面波假设下的，而在推导新方程过程中首先要验证的就是方程的适用条件和敏感性，如果这一点不满足，新的方程也就没有意义。
[0198] 实施例2。图10为本发明基于流体体积压缩系数的叠前地震反演方法的的技术路线图。本实施例的主要流程：
[0199] 基于孔隙弹性理论的流体体积压缩系数构建；
[0200] 推导基于流体体积压缩系数的地震反射特征方程；
[0201] 贝叶斯框架下流体体积压缩系数叠前地震反演方法；
[0202] 实际资料应用。
[0203] 步骤101：基于孔隙弹性理论的流体体积压缩系数构建
[0204] 一定温度下，单位压强的体积相对缩小率
[0205]
[0206] Batzle和Han通过研究发现，影响Gassmann流体项f取值的主要因素是孔隙流体体积模量与固体骨架孔隙度。Dehua Han和Batzle等人通过对碎屑岩进行岩石物理统计，研究了Biot-Gassmann理论中的孔隙流体与岩石骨架的固体效应(孔隙度，矿物模量等)对岩石模量信息的耦合作用，提出了突出流体体积压缩系数Cf的Gassmann流体项的经验公式：
[0207] f＝G(φ)/Cf (2)。
[0208] 其中， 其中增益函数G(φ)表示岩石骨架矿物与孔隙度的综合作用。
[0209] 结合测井资料，将流体体积压缩系数与常规的基于单相介质理论的流体因子类型作流体敏感性对比，如图5所示，可以看到在该工区流体体积压缩系数具有很高的敏感性。
[0210] 步骤102：推导基于流体体积压缩系数的地震反射特征方程
[0211] 基于Biot-Gassmann理论，Russell等人对饱含流体多孔介质的AVO理论进行了研究，在2006年的研究报告中首次提出了包含Gassmann流体项的反射系数近似公式，并且在2011年在《Geophysics》正式发表论文发对其进行了讨论，Russell近似公式如下所示：
[0212]
[0213] 根据相关岩石物理模型以及计算方法，最终得到流体体积压缩系数近似公式如下所示：
[0214]
[0215] 借鉴Connolly推导弹性阻抗的思想，用弹性阻抗表示反射系数，得到：
[0216]
[0217] 将公式(12)代入(11)，得到：
[0218]
[0219] 将弹性参数的相对变化量用对数形式表示，可以得到：
[0220]
[0221] 进一步变形，可以得到：
[0222]
[0223] 即
[0224]
[0225] 对上式两边去积分并将其指数化，消掉等式两边的微分项和对数项，进一步去积分常数为0，得到：
[0226] EI(θ)＝(Cf)a(θ)(fm)b(θ)(ρ)c(θ)(φ)d(θ) (17)。
[0227]
[0228] 其中，
[0229] 与常规弹性阻抗公式类似，公式(17)也存在数值量纲随角度变化的问题。在此引入四个参考常数，即A0，Cf0，fm0，ρ0以及φ0。将公式(17)进行标准化处理，可以得到标准化的基于流体体积压缩系数的弹性阻抗方程。
[0230] EI(θ)＝A0(Cf/Cf0)a(θ)(fm/fm0)b(θ)(ρ/ρ0)c(θ)(φ/φ0)d(θ) (18)。
[0231] Cf0，fm0，ρ0以及φ0分别定义为Cf0，fm，ρ以及φ的平均值。A0为标准化因子，具体表达式为：
[0232]
[0233] 通过不同模型可以测试(19)式的精确程度，如图6、图7所示，采集的地震数据偏移距/入射角在一定范围内，一般最大入射角小于35度。图6、7中粗虚线为精确方程反射系数曲线，描述了理论上的AVO特征，较细虚线代表AKI近似方程求取的反射系数曲线，描述了线性近似下理论AVO特征，实线为本实施例推导的新近似方程得到的反射系数曲线。在入射角小于30°范围内，从图6可以看到3条曲线几乎保持一致，图7中虽然3条曲线存在差异，而且其差异的主要原因是模型二物性差异明显，因此可以认为误差允许的范围内，因此推到的新方程在一定角度范围内具有较高的精确度，因此可以将其运用到实际生产中。
[0234] 步骤103：贝叶斯框架下流体体积压缩系数叠前地震反演方法
[0235] 最后，反射系数线性近似公式是建立弹性参数与地震数据之间关系的桥梁，不同的反射系数线性近似公式所包含的弹性参数种类以及信息量存在差异。本发明采用流体体积压缩系数AVO反射系数线性近似公式建立正演模型。将公式(11)按照入射角度的不同表示为矩阵形式为：
[0236]
[0237] 其中，ai(i＝1,2…,m)，bi(i＝1,2…,m)，ci(i＝1,2…,m)和di(i＝1,2…,m)分别表示第i个入射角度的相应系数。将其推广到具有m个入射角度，n个界面的情况，并且将矩阵进行块化处理，可以得到：
[0238]
[0239] 其中，Ri(i＝1,2…,m)表示第i个入射角度的反射系数向量，由n个元素组成；Ai(i＝1,2…,m)，Bi(i＝1,2…,m)，Ci(i＝1,2…,m)和Di(i＝1,2…,m)分别表示第i个入射角度对应的正演系数矩阵，分别是n×n维的斜对角矩阵； Rρ和Rφ则分别表示去岩石骨架流体因子、固体刚性参数、密度以及孔隙度相对变化率向量，分别由n个元素组成。
[0240] 基于地震记录符合褶积模型的假设，引入子波矩阵W，则公式(21)进一步变为：
[0241]
[0242] 其中，di(i＝1,2…,m)表示为第i个入射角度的地震数据组成的列向量，都包含n个元素。
[0243] 步骤104:实际资料应用
[0244] 在贝叶斯理论框架下，利用弹性阻抗反演方法得到流体体积压缩系数弹性阻抗反演结果，利用(22)式建立弹性阻抗与流体体积压缩系数的直接关系，即可以得到流体体积压缩系数最终结果。在具体的求解过程中，基于贝叶斯理论框架，将似然函数与先验函数代入后验概率密度函数，可以得到：
[0245]
[0246] 省略常数，求解最大后验概率，可以得到反演目标函数，具体形式如下：
[0247] (G'TG'+2Q)R'＝G'Td (30)。
[0248] 其中， 公式左边第一项主要用来约束正演记录与实际地震记录之间的相近程度，第二项则是四变量柯西正则约束项，主要用来约束反演参数的稀疏程度，然后采用迭代重加权最小二乘算法(Iterative Reweighted Least Squares,IRLS)对反演方程进行目标寻优。
[0249] 利用得到的弹性阻抗结果提取流体体积压缩系数，如图8、图9所示，A、B、C为3口井，其中有一口井为水井，其他两口为气井。通过本发明的方法得到的流体体积压缩系数反演剖面上很好地将气水分开，提高了流体识别的精度。
[0250] 本实施例一种基于流体体积压缩系数的叠前地震反演方法。气体、液体以及固体都具有不同程度的压缩性，体积压缩系数定义为在一定温度下，压强增加一个单位体积的相对缩小率，可以定量地表征流体与固体的压缩性，一般情形下，气体、液体、固体的压缩系数存在如下关系：气体>液体>固体，因此从构建孔隙流体敏感参数出发，考虑孔隙流体(气体、液体)与岩石骨架的压缩系数之间的差异，将流体体积压缩系数作为流体因子，对提高储层流体识别精度具有实际意义。首先基于岩石物理理论以及测井资料，对流体体积压缩系数进行敏感性分析，用于验证流体体积压缩系数流体识别精度。然后基于孔隙弹性理论，利用理论和经验岩石物理模型，充分考虑岩石模量、孔隙流体及孔隙大小对储层流体影响，推导基于流体体积压缩系数的平面波反射特征近似方程，并分析其适用条件和敏感性，以便于在进行反演时提供满足条件的输入数据以及选取合适的反演参数。本实施例中AVO近似方程是在平面波假设的基础上推导出来的，精确的平面波公式难以在实际中运用，所以需要在一定假设下推导近似方程，因此可以说AVO近似方程是平面波方程的一种特殊形式。最后推导基于流体体积压缩系数的弹性阻抗方程，在贝叶斯理论框架下建立起地震-弹性阻抗-流体体积压缩系数反演流程，形成了基于流体体积压缩系数的叠前地震反演方法。
[0251] 最后说明的是，以上优选实施例仅用以说明本发明的技术方案而非限制，尽管通过上述优选实施例已经对本发明进行了详细的描述，但本领域技术人员应当理解，可以在形式上和细节上对其作出各种各样的改变，而不偏离本发明
权利要求书所限定的范围。
