CN104459778A - 一种基于双相介质固液解耦的叠前地震反演方法及系统 - Google Patents

Abstract

本发明提供一种基于双相介质固液解耦的叠前地震反演方法及系统，所述方法包括：采集目的储层的叠前地震资料、测井资料以及地质成果资料；构建固液解耦近似方程；根据所述的叠前地震资料、测井资料、地质成果资料以及所述的固液解耦近似方程提取子波；根据所述的叠前地震资料、测井资料、地质成果资料以及子波确定绝对弹性阻抗数据体；根据所述的绝对弹性阻抗数据体确定流体体积模量数据体；根据所述的流体体积模量数据体识别当前储层孔隙中流体的类型。得到的流体体积模量作为流体因子不仅具有更高的流体指示灵敏性，而且有效的提高了储层流体识别和储层特征描述的可靠性。

Description

一种基于双相介质固液解耦的叠前地震反演方法及系统

技术领域

本发明关于石油地球物理勘探领域，特别是关于叠前地震的反演技术，具体的讲是一种基于双相介质固液解耦的叠前地震反演方法及系统。

背景技术

在石油工业中，利用流体因子进行地下储层流体识别是业界常用的烃类描述与预测方式。随着全球地球物理勘探的重点从构造油气藏向岩性油气藏转变，早先在流体识别领域应用广泛的“亮点”及AVO属性分析技术已无法满足现阶段复杂岩性油气藏的储层流体识别需求。

针对上述问题，现有技术中许多学者提出了不同类型的流体因子，主要包括两类：

(1)通过岩石物理参数的权差运算来突出体现流体异常效应，如交会图角、流体因子角、泊松阻抗等，但是此类流体因子的岩石物理意义模糊且普适性较差；

(2)将对岩石孔隙流体敏感的弹性参数直接作为流体因子，如拉梅参数，泊松比，Russell流体因子以及Gassmann流体项f等，该类流体因子岩石物理意义明确且普适性较好，在现阶段储层流体识别中占据主要地位，但是在岩石固结程度较为复杂的地区应用该类流体因子时，受岩石固液两相组分的耦合效应影响，特别是受孔隙度影响，在实际储层的孔隙流体判识中会出现流体识别假象。另外，常规流体因子的计算多依赖弹性参数之间的间接代数运算，受地震反问题固有的“病态”问题影响，最终计算的流体因子无法避免会引入累计误差，降低了储层流体指示作用的可靠性。

因此，随着以复杂岩性因素为主导作用的岩性地层油气藏勘探程度的加深，急需研究一种能够准确反应孔隙流体类型的高灵敏流体因子类型以及更为可靠的流体因子直接反演方法。

发明内容

为了克服现有技术存在的上述问题，本发明提供了一种基于双相介质固液解耦的叠前地震反演方法及系统，以岩石固液双相效应解耦并反演完全体现了流体效应的流体体积模量参数K_f为出发点，利用叠前弹性阻抗反演实现了流体因子的直接提取，得到的流体体积模量作为流体因子不仅具有更高的流体指示灵敏性，而且有效的提高了储层流体识别和储层特征描述的可靠性。

本发明的目的之一是，提供一种基于双相介质固液解耦的叠前地震反演方法，包括：采集目的储层的叠前地震资料、测井资料以及地质成果资料；构建固液解耦近似方程；根据所述的叠前地震资料、测井资料、地质成果资料以及所述的固液解耦近似方程提取子波；根据所述的叠前地震资料、测井资料、地质成果资料以及子波确定绝对弹性阻抗数据体；根据所述的绝对弹性阻抗数据体确定流体体积模量数据体；根据所述的流体体积模量数据体识别当前储层孔隙中流体的类型。

本发明的目的之一是，提供了一种基于双相介质固液解耦的叠前地震反演系统，包括：资料采集装置，用于采集目的储层的叠前地震资料、测井资料以及地质成果资料；固液解耦近似方程构建装置，用于构建固液解耦近似方程；子波提取装置，用于根据所述的叠前地震资料、测井资料、地质成果资料以及所述的固液解耦近似方程提取子波；弹性阻抗数据体确定装置，用于根据所述的叠前地震资料、测井资料、地质成果资料以及子波确定绝对弹性阻抗数据体；流体体积模量确定装置，用于根据所述的绝对弹性阻抗数据体确定流体体积模量数据体；流体类型识别装置，用于根据所述的流体体积模量数据体识别当前储层孔隙中流体的类型。

本发明的有益效果在于，提供了一种基于双相介质固液解耦的叠前地震反演方法及系统，以岩石固液双相效应解耦并反演完全体现了流体效应的流体体积模量参数K_f为出发点，利用叠前弹性阻抗反演实现了流体因子的直接提取，得到的流体体积模量作为流体因子不仅具有更高的流体指示灵敏性，而且有效的提高了储层流体识别和储层特征描述的可靠性，常规的流体因子岩石物理意义模糊且普适性较差，而且受岩石双相特性的耦合效应影响，在表征岩石流体效应的时候易受岩石固体效应(特别是孔隙度参数)影响，造成流体识别假象，而流体体积模量K_f只与岩石孔隙流体有关，不受其他岩石固体效应影响，作为一项对储层孔隙流体类型敏感的弹性参数在流体识别领域的效果更佳，利用叠前弹性阻抗反演实现了流体因子的直接提取，可以避免引入累计误差，提高了储层流体指示作用的可靠性。

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

附图说明

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

图1是利用某工区实际数据计算得到的常用流体因子的流体指示系数比较图；

图2是本发明提供的基于固液解耦近似方程的流体因子叠前地震反演方法在具体实施例中的流程图；

图3是根据某工区实际地震数据进行基于固液解耦近似方程的流体因子叠前地震反演后得到的流体因子K_f的剖面示意图；

图4为本发明实施例提供的一种基于双相介质固液解耦的叠前地震反演方法的实施方式一的流程图；

图5为本发明实施例提供的一种基于双相介质固液解耦的叠前地震反演方法的实施方式二的流程图；

图6为本发明实施例提供的一种基于双相介质固液解耦的叠前地震反演方法的实施方式三的流程图；

图7为本发明实施例提供的一种基于双相介质固液解耦的叠前地震反演方法的实施方式四的流程图；

图8为图4中的步骤S103的具体流程图；

图9为图4中的步骤S104的具体流程图；

图10为本发明实施例提供的一种基于双相介质固液解耦的叠前地震反演系统的实施方式一的结构框图；

图11为本发明实施例提供的一种基于双相介质固液解耦的叠前地震反演系统的实施方式二的结构框图；

图12为本发明实施例提供的一种基于双相介质固液解耦的叠前地震反演系统的实施方式三的结构框图；

图13为本发明实施例提供的一种基于双相介质固液解耦的叠前地震反演系统的实施方式四的结构框图；

图14为本发明实施例提供的一种基于双相介质固液解耦的叠前地震反演系统中的子波提取装置300的具体结构框图；

图15为本发明实施例提供的一种基于双相介质固液解耦的叠前地震反演系统中的弹性阻抗数据体确定装置400的具体结构框图。

具体实施方式

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

在现阶段，储层流体识别中占据主要地位的流体因子是将对岩石孔隙流体敏感的弹性参数直接作为流体因子，在岩石固结程度较为复杂的地区应用该类流体因子时，受岩石固液两相组分的耦合效应影响，特别是受孔隙度影响，在储层的孔隙流体判识中会出现流体识别假象。另外，常规流体因子的计算多依赖弹性参数之间的间接代数运算，受地震反问题固有的“病态”问题影响，最终计算的流体因子无法避免会引入累计误差，降低了储层流体指示作用的可靠性。本发明以实现岩石固液双相效应解耦为出发点，提出将完全体现流体效应、不受孔隙度影响的流体因子，并利用叠前弹性阻抗反演实现了该流体因子的直接提取。图4为本发明提出的一种基于双相介质固液解耦的叠前地震反演方法的实施方式一的具体流程图，由图4可知，在实施方式一中，所述的方法包括：

S101：采集目的储层的叠前地震资料、测井资料以及地质成果资料。

叠前地震资料的采集在具体的实施方式中可通过如下方式：在勘探工区中，根据目的层的深度和地质特点，考虑动校正拉伸、干扰波、多次波等因素来设计宽方位角的三维观测系统，保证有足够的偏移距、方位角。经过激发、接收，得到可以满足AVO分析需求的宽方位大角度的叠前地震资料。

测井资料的采集在具体的实施方式中可通过如下方式：在勘探工区中进行全波列测井，得到测井资料，主要包含纵波速度、横波速度、密度等全波列测井曲线，孔隙度等解释成果曲线，测井层位，录井资料等。

S102：构建固液解耦近似方程。

AVO(Amplitude Versus Offset，振幅随偏移距的变化)近似公式Aki公式是进行弹性参数估算的基础，本发明基于多孔弹性介质理论以及Nur岩石物理模型，推导得到的流体体积模量分离的AVO近似公式即固液解耦近似方程，其为基于地震资料的流体体积模量参数估算提供了理论基础。

Russell等人从Aki-Richards近似公式出发，推导了突出储层流体特征的反射特征近似公式：

$R_{\mathrm{PP}}\left(\mathrm{\θ}\right)=\left[(1-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2})\frac{{\sec }^2\mathrm{\θ}}{4}\right]\frac{\mathrm{\Δf}}{f}+[\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{4{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}]\frac{\mathrm{\Δ\μ}}{\mathrm{\μ}}+[\frac{1}{2}-\frac{{\sec }^2\mathrm{\θ}}{4}]\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}---\left(1\right)$

Dehua Han等人提出了Gassmann流体项的经验公式：

f＝G(φ)K_f      (2)

其中， $G\left(\mathrm{\φ}\right)=\frac{{(1-K_n)}^2}{\mathrm{\φ}},K_n=\frac{K_{\mathrm{dry}}}{K_m}.$

将公式(2)代入公式(1)，并考虑到剪切模量不受孔隙流体的影响，利用干岩石剪切模量μ_dry替换μ，进行相应变换，可以得到：

$R_{\mathrm{PP}}\left(\mathrm{\θ}\right)=\left[(1-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2})\frac{{\sec }^2\mathrm{\θ}}{4}\right]\frac{\mathrm{\Δ}\left(G\left(\mathrm{\φ}\right)K_f\right)}{G\left(\mathrm{\φ}\right)K_f}+[\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{4{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}]\frac{\mathrm{\Δ\μ}}{\mathrm{\μ}}+[\frac{1}{2}-\frac{{\sec }^2\mathrm{\θ}}{4}]\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}---\left(3\right)$

Nur提出的临界孔隙度模型表达式如下式所示：

$\left\{\begin{array}{c}{K}_{\mathrm{dry}}=K_m(1-\frac{\mathrm{\φ}}{{\mathrm{\φ}}_c})\\ {\mathrm{\μ}}_{\mathrm{dry}}={\mathrm{\μ}}_m(1-\frac{\mathrm{\φ}}{{\mathrm{\φ}}_c})\end{array}\right.---\left(4\right)$

其中，φ_c表示临界孔隙度，K_dry表示干岩石的体积模量，μ_dry表示干岩石的剪切模量，K_m表示固体矿物基质的体积模量，μ_m表示矿物基质的剪切模量。

结合公式(4)对公式(3)进一步展开可以得到：

$\begin{array}{c}{R}_{\mathrm{PP}}\left(\mathrm{\θ}\right)=\left[(1-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2})\frac{{\sec }^2\mathrm{\θ}}{4}\right](\frac{\mathrm{\ΔG}\left(\mathrm{\φ}\right)}{G\left(\mathrm{\φ}\right)}+\frac{\mathrm{\Δ}{K}_f}{K_f})+[\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{4{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}](\frac{\mathrm{\Δ}{\mathrm{\μ}}_m}{{\mathrm{\μ}}_m}+\frac{\mathrm{\Δ}({\mathrm{\φ}}_c-\mathrm{\φ})}{{\mathrm{\φ}}_c-\mathrm{\φ}})\\ +[\frac{1}{2}-\frac{{\sec }^2\mathrm{\θ}}{4}]\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}\end{array}---\left(5\right)$

将代入G(φ)，进一步展开可以得到：

$\begin{array}{c}{R}_{\mathrm{PP}}\left(\mathrm{\θ}\right)=\left[(1-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2})\frac{{\sec }^2\mathrm{\θ}}{4}\right]+\frac{\mathrm{\Δ}{K}_f}{K_f}+[\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{4{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}]\frac{\mathrm{\Δ}{\mathrm{\μ}}_m}{{\mathrm{\μ}}_m}+[\frac{{\sec }^2\mathrm{\θ}}{4}-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{2{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}+\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}]\frac{\mathrm{\Δ\φ}}{\mathrm{\φ}}\\ +[\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{4{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^1}{\sin }^2\mathrm{\θ}]\frac{\mathrm{\Δ\φ}({\mathrm{\φ}}_c-\mathrm{\φ})}{\mathrm{\φ}({\mathrm{\φ}}_c-\mathrm{\φ})}+[\frac{1}{2}-\frac{{\sec }^2\mathrm{\θ}}{4}]\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}\end{array}---\left(6\right)$

令f_m＝φμ，并利用公式(4)，最终得到的固液解耦近似方程如下所示：

$\begin{array}{c}{R}_{\mathrm{PP}}\left(\mathrm{\θ}\right)=\left[(1-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2})\frac{{\sec }^2\mathrm{\θ}}{4}\right]\frac{\mathrm{\Δ}{K}_f}{K_f}+[\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{4{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}]\frac{\mathrm{\Δ}{f}_m}{f_m}+[\frac{1}{2}-\frac{{\sec }^2\mathrm{\θ}}{4}]\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}\\ +(\frac{{\sec }^2\mathrm{\θ}}{4}-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{2{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}+\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ})\frac{\mathrm{\Δ\φ}}{\mathrm{\φ}}\end{array}---\left(7\right)$

其中，f_m＝φμ为地下岩石的固体项，θ为入射角度，K_f为地下岩石的流体体积模量，f_m为地下岩石的固体项，φ为地下岩石的孔隙度，ρ为地下岩石的密度，μ为地下岩石的剪切模量，△K_f为界面两侧的流体体积模量的差值，△f_m为界面两侧固体项的差值，△ρ为界面两侧的密度的差值，△φ为界面两侧的孔隙度的差值，为干岩石纵横波速度比的平方，为饱和岩石纵横波速度比的平方。

S103：根据所述的叠前地震资料、测井资料、地质成果资料以及所述的固液解耦近似方程提取子波。振幅随偏移距的变化导致了子波随偏移距的空变特性，为消除不同角度道集间的相互影响，在此需要针对不同的入射角度对部分角度叠加数据体分别提取不同的子波。图8为步骤S103的具体流程图。

S104：根据所述的叠前地震资料、测井资料、地质成果资料以及子波确定绝对弹性阻抗数据体。综合叠前地震数据、测井数据、层位数据和子波，利用约束稀疏脉冲反演算法计算相对弹性阻抗数据体，然后通过补偿低频阻抗成分，最终得到绝对弹性阻抗数据体。图9为步骤S104的具体流程图。

S105：根据所述的绝对弹性阻抗数据体确定流体体积模量数据体。根据弹性阻抗与弹性参数之间的线性求解关系，利用四个角度的弹性阻抗数据体提取流体体积模量。

S106：根据所述的流体体积模量数据体识别当前储层孔隙中流体的类型。提取的流体体积模量K_f数据体与储层流体体积模量的理论值进行对比，可以识别储层孔隙中流体的类型。

图5为本发明实施例提供的一种基于双相介质固液解耦的叠前地震反演方法的实施方式二的流程图，由图5可知，在实施方式二中，该方法具体包括：

S201：采集目的储层的叠前地震资料、测井资料以及地质成果资料。

S202：对所述的叠前地震资料进行保幅处理；

S203：对经过保幅处理的叠前地震资料进行去噪处理；

S204：根据目的储层对经过去噪处理的叠前地震资料进行角度叠加，生成四个入射角度的角度叠加道集。

步骤S202至S204是反演之前需要对叠前地震资料进行的预处理。在保证地震资料品质的前提下根据目的储层进行相应部分角度叠加，最终生成四个角度的部分角度叠加道集。此处提及的四个角度，在具体的实施方式中，若叠前地震数据的角度范围大、大角度数据信噪比较高，可选择6度、14度、22度、30度；若叠前地震数据的角度范围不大、大角度数据信噪比较低，可选择5度、10度、15度、20度。

S205：构建固液解耦近似方程。

S206：根据所述的叠前地震资料、测井资料、地质成果资料以及所述的固液解耦近似方程提取子波。振幅随偏移距的变化导致了子波随偏移距的空变特性，为消除不同角度道集间的相互影响，在此需要针对不同的入射角度对部分角度叠加数据体分别提取不同的子波。图8为步骤S206的具体流程图。

S207：根据所述的叠前地震资料、测井资料、地质成果资料以及子波确定绝对弹性阻抗数据体。综合叠前地震数据、测井数据、层位数据和子波，利用约束稀疏脉冲反演算法计算相对弹性阻抗数据体，然后通过补偿低频阻抗成分，最终得到绝对弹性阻抗数据体。图9为步骤S207的具体流程图。

S208：根据所述的绝对弹性阻抗数据体确定流体体积模量数据体。根据弹性阻抗与弹性参数之间的线性求解关系，利用四个角度的弹性阻抗数据体提取流体体积模量。

S209：根据所述的流体体积模量数据体识别当前储层孔隙中流体的类型。提取的流体体积模量K_f与储层流体体积模量的理论值进行对比，可以识别储层孔隙中流体的类型。

图6为本发明实施例提供的一种基于双相介质固液解耦的叠前地震反演方法的实施方式三的流程图，由图6可知，在实施方式三中，该方法具体包括：

S301：采集目的储层的叠前地震资料、测井资料以及地质成果资料。

S302：对所述的叠前地震资料进行保幅处理；

S303：对经过保幅处理的叠前地震资料进行去噪处理；

S304：根据目的储层对经过去噪处理的叠前地震资料进行角度叠加，生成四个入射角度的角度叠加道集。

步骤S302至S304是反演之前需要对叠前地震资料进行的预处理。在保证地震资料品质的前提下根据目的储层进行相应部分角度叠加，最终生成四个角度的部分角度叠加道集。此处提及的四个角度，在具体的实施方式中，若叠前地震数据的角度范围大、大角度数据信噪比较高，可选择6度、14度、22度、30度；若叠前地震数据的角度范围不大、大角度数据信噪比较低，可选择5度、10度、15度、20度。

S305：对所述的测井资料进行环境校正；

S306：对经过环境校正后的测井资料进行奇异值消除；

S307：根据经过奇异值消除后的测井资料确定所述目的储层的干岩石纵横波速度比、饱和岩石纵横波速度比。具体的实施方式中，在勘探工区中进行全波列测井，得到的测井资料需要经过环境校正、奇异值消除等预处理后，结合岩石物理经验公式或者实验统计结果计算干岩石纵横波速度比γ_dry、饱和岩石纵横波速度比γ_sat的数值。

S308：根据所述的目的储层的干岩样纵横波速度比以及所述的测井资料确定与所述测井资料对应的流体体积模量曲线以及剪切模量曲线。在具体的实施方式中，经过换算得到相应的流体体积模量为K_f，剪切模量为μ曲线，后续用来计算井旁道的基于固液解耦近似方程的弹性阻抗伪测井曲线。

S309：构建固液解耦近似方程。

S310：根据所述的叠前地震资料、测井资料、地质成果资料以及所述的固液解耦近似方程提取子波。振幅随偏移距的变化导致了子波随偏移距的空变特性，为消除不同角度道集间的相互影响，在此需要针对不同的入射角度对部分角度叠加数据体分别提取不同的子波。图8为步骤S310的具体流程图。

S311：根据所述的叠前地震资料、测井资料、地质成果资料以及子波确定绝对弹性阻抗数据体。综合叠前地震数据、测井数据、层位数据和子波，利用约束稀疏脉冲反演算法计算相对弹性阻抗数据体，然后通过补偿低频阻抗成分，最终得到绝对弹性阻抗数据体。图9为步骤S311的具体流程图。

S312：根据所述的绝对弹性阻抗数据体确定流体体积模量数据体。根据弹性阻抗与弹性参数之间的线性求解关系，利用四个角度的弹性阻抗数据体提取流体体积模量。

S313：根据所述的流体体积模量数据体识别当前储层孔隙中流体的类型。提取的流体体积模量K_f与储层流体体积模量的理论值进行对比，可以识别储层孔隙中流体的类型。

图7为本发明实施例提供的一种基于双相介质固液解耦的叠前地震反演方法的实施方式四的流程图，由图7可知，在实施方式四中，该方法具体包括：

S401：采集目的储层的叠前地震资料、测井资料以及地质成果资料。

S402：对所述的叠前地震资料进行保幅处理；

S403：对经过保幅处理的叠前地震资料进行去噪处理；

S404：根据目的储层对经过去噪处理的叠前地震资料进行角度叠加，生成四个入射角度的角度叠加道集。

步骤S402至S404是反演之前需要对叠前地震资料进行的预处理。在保证地震资料品质的前提下根据目的储层进行相应部分角度叠加，最终生成四个角度的部分角度叠加道集。此处提及的四个角度，在具体的实施方式中，若叠前地震数据的角度范围大、大角度数据信噪比较高，可选择6度、14度、22度、30度；若叠前地震数据的角度范围不大、大角度数据信噪比较低，可选择5度、10度、15度、20度。

S405：对所述的测井资料进行环境校正；

S406：对经过环境校正后的测井资料进行奇异值消除；

S407：根据经过奇异值消除后的测井资料确定所述目的储层的干岩石纵横波速度比、饱和岩石纵横波速度比。具体的实施方式中，在勘探工区中进行全波列测井，得到的测井资料需要经过环境校正、奇异值消除等预处理后，结合岩石物理经验公式或者实验统计结果计算干岩石纵横波速度比γ_dry、饱和岩石纵横波速度比γ_sat的数值。

S408：根据所述的目的储层的干岩样纵横波速度比以及所述的测井资料确定与所述测井资料对应的流体体积模量曲线以及剪切模量曲线。在具体的实施方式中，经过换算得到相应的流体体积模量为K_f，剪切模量为μ曲线，后续用来计算井旁道的基于固液解耦近似方程的弹性阻抗伪测井曲线。

S409：根据所述的叠前地震资料、测井资料、地质成果资料对目的储层进行精细层位追踪，得到地震层位数据。

在具体的实施方式中，结合工区的地质背景，在保证目的层位符合工区地质沉积模式与层序格架的基础上，对研究的目的层进行精细地质解释，解释得到的地质成果资料不仅用来标识时深关系，还可结合弹性阻抗伪测井曲线建立低频阻抗数据体。

S410：构建固液解耦近似方程。

S411：根据所述的叠前地震资料、测井资料、地质成果资料以及所述的固液解耦近似方程提取子波。振幅随偏移距的变化导致了子波随偏移距的空变特性，为消除不同角度道集间的相互影响，在此需要针对不同的入射角度对部分角度叠加数据体分别提取不同的子波。图8为步骤S411的具体流程图。由图8可知，该步骤具体包括：

S501：从所述的测井资料中提取出地下岩石的密度以及孔隙度。在具体的实施方式中，地下岩石的密度通过ρ来表示，地下岩石的孔隙度通过φ来表示。

S502：根据所述的剪切模量曲线以及所述的孔隙度确定地下岩石的固体项，在具体的实施方式中，地下岩石的固体项通过f_m来表示，则f_m＝φμ。

S503：基于所述的固液解耦近似方程推导出弹性阻抗方程。

与普通弹性阻抗方程推导过程类似，固液解耦近似方程的反射系数可以表示为：

$R_{\mathrm{pp}}\left(\mathrm{\θ}\right)\≈\frac{1}{2}\frac{\mathrm{\ΔEI}}{\mathrm{EI}}\≈\frac{1}{2}\mathrm{\Δ}\ln \left(\mathrm{EI}\right)---\left(8\right)$

将公式(7)代入上式，得到：

$\begin{array}{c}\frac{1}{2}\mathrm{\Δ}\ln \left(\mathrm{EI}\right)=\left[(1-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2})\frac{{\sec }^2\mathrm{\θ}}{4}\right]\frac{\mathrm{\Δ}{K}_f}{K_f}+[\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{4{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}]\frac{\mathrm{\Δ}\left(f_m\right)}{f_m}+[\frac{1}{2}-\frac{{\sec }^2\mathrm{\θ}}{4}]\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}\\ +(\frac{{\sec }^2\mathrm{\θ}}{4}-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{2{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}+\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ})\frac{\mathrm{\Δ\φ}}{\mathrm{\φ}}\end{array}---\left(9\right)$

将弹性参数的相对变化量用对数形式表示，可以得到：

$\begin{array}{c}\mathrm{\Δ}\ln \left(\mathrm{EI}\right)=\left[(1-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2})\frac{{\sec }^2\mathrm{\θ}}{2}\right]\mathrm{\Δ}\ln \left(K_f\right)+[\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{2{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{4}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}]\mathrm{\Δ}\ln \left(f_m\right)+[1-\frac{{\sec }^2\mathrm{\θ}}{2}]\mathrm{\Δ}\ln \left(\mathrm{\ρ}\right)\\ +(\frac{{\sec }^2\mathrm{\θ}}{2}-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{2{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}+\frac{4}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ})\mathrm{\Δ}\ln \left(\mathrm{\φ}\right)\end{array}---\left(10\right)$

进一步变形，可以得到：

$\mathrm{\Δ}\ln \left(\mathrm{EI}\right)=\mathrm{\Δ}\ln \left\{{K_f}^{\left[(\frac{1}{2}-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{2{\mathrm{\γ}}_{\mathrm{sat}}^2}){\sec }^2\mathrm{\θ}\right]}{f_m}^{[\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{2{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{4}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}]}{\mathrm{\ρ}}^{[1-\frac{{\sec }^2\mathrm{\θ}}{2}]}{\mathrm{\φ}}^{(\frac{{\sec }^2\mathrm{\θ}}{2}-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}+\frac{4}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ})}\right\}---\left(11\right)$

两边取积分并将其指数化，消掉等式两边的微分项和对数项，进一步取积分常数为0，得到：

$\mathrm{EI}={K_f}^{\left[(\frac{1}{2}-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{2{\mathrm{\γ}}_{\mathrm{sat}}^2}){\sec }^2\mathrm{\θ}\right]}{f_m}^{[\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{2{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{4}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}]}{\mathrm{\ρ}}^{[1-\frac{{\sec }^2\mathrm{\θ}}{2}]}{\mathrm{\φ}}^{(\frac{{\sec }^2\mathrm{\θ}}{2}-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}+\frac{4}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ})}---\left(12\right)$

即：

$\mathrm{EI}\left(\mathrm{\θ}\right)={K_f}^a{f_m}^b{\mathrm{\φ}}^c{\mathrm{\ρ}}^d---\left(13\right)$

AVO近似公式是叠前反演的基础，本发明以包含流体体积模量K_f的固液解耦近似方程为基础，推导标准化的弹性阻抗方程，消除了不同角度之间存在的量纲不统一问题。通过叠前弹性阻抗反演得到相对阻抗数据体，经过低频分量补偿得到绝对弹性阻抗数据体，最终利用弹性参数提取直接得到反映孔隙流体效应的流体因子K_f，从而实现敏感流体因子(即流体体积模量K_f)的叠前地震反演方法。

对公式(13)进行标准化可以推导出弹性阻抗方程即为：

$\mathrm{EI}\left(\mathrm{\θ}\right)={\mathrm{EI}}_0{\left(\frac{K_f}{K_{f0}}\right)}^a{\left(\frac{f_m}{f_{m0}}\right)}^b{\left(\frac{\mathrm{\φ}}{{\mathrm{\φ}}_0}\right)}^c{\left(\frac{\mathrm{\ρ}}{{\mathrm{\ρ}}_0}\right)}^d---\left(14\right)$

其中， $a=(1-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2})\frac{{\sec }^2\mathrm{\θ}}{2},b=\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{2{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{4}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ},c=\frac{{\sec }^2\mathrm{\θ}}{2}-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}+\frac{4}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ},$ $d=1-\frac{{\sec }^2\mathrm{\θ}}{2};{\mathrm{EI}}_0={\left(K_f\right)}^{-\frac{1}{2}}{f}_m{\mathrm{\φ}}^{-\frac{3}{2}}{\mathrm{\ρ}}^{\frac{1}{2}}{\mathrm{\γ}}_{\mathrm{sat}}{({\mathrm{\γ}}_{\mathrm{sat}}^2-{\mathrm{\γ}}_{\mathrm{dry}}^2)}^{\frac{1}{2}},$ K_f、f_m、ρ和φ分别表示地下岩石的流体体积模量，固体项，密度和孔隙度，θ表示入射角度，表示干岩石纵横波速度比的平方，表示饱和岩石纵横波速度比的平方，K_f0、f_m0、ρ₀和φ₀分别为K_f、f_m、ρ和φ的平均值。

S504：根据所述的入射角度、地下岩石的密度、固体项、孔隙度、流体体积模量曲线以及弹性阻抗方程确定弹性阻抗伪测井曲线。在具体的实施方式中，对叠前地震资料预处理后得到四个入射角度，将每个入射角度、地下岩石的密度、固体项、孔隙度、流体体积模量曲线带入弹性阻抗方程中，皆能确定出一条弹性阻抗伪测井曲线，因此四个入射角度即能确定出四条弹性阻抗伪测井曲线。在反演中，弹性阻抗伪测井曲线不仅可以用来作为必备的先验约束条件，还可以以此结合地质解释层位建立低频模型来补偿地震波传播过程中损失的低频成分。

S505：针对不同的入射角度，根据入射角度对应的角度叠加道集、弹性阻抗伪测井曲线、地震层位数据提取不同的子波。在具体的实施方式中，振幅随偏移距的变化导致了子波随偏移距的空变特性，为消除不同角度道集间的相互影响，在此需要针对不同的入射角度对部分角度叠加道集分别提取不同的子波。

由图7可知，在实施方式四中，该方法还包括：

S412：根据所述的叠前地震资料、测井资料、地质成果资料以及子波确定绝对弹性阻抗数据体。综合叠前地震数据、测井数据、层位数据和子波，利用约束稀疏脉冲反演算法计算相对弹性阻抗数据体，然后通过补偿低频阻抗成分，最终得到绝对弹性阻抗数据体。图9为步骤S412的具体流程图。由图9可知，该步骤具体包括：

S601：根据所述的叠前地震资料、测井资料、地震层位数据以及子波确定出相对弹性阻抗数据体。在具体的实施方式中，综合叠前地震资料、测井资料、地震层位数据和子波，利用约束稀疏脉冲反演算法计算相对弹性阻抗数据体。

S602：根据所述的弹性阻抗伪测井曲线以及地震层位数据构建低频阻抗数据体。在反演中，弹性阻抗伪测井曲线不仅可以用来作为必备的先验约束条件，还可以以此结合地震层位数据建立低频阻抗数据体来补偿地震波传播过程中损失的低频成分。

S603：根据所述的低频阻抗数据体对所述的弹性阻抗数据体进行补偿，得到绝对弹性阻抗数据体。即通过补偿低频阻抗成分，最终得到绝对弹性阻抗数据体。

由图7可知，在实施方式四中，该方法还包括：

S413：根据所述的绝对弹性阻抗数据体确定流体体积模量数据体。根据弹性阻抗与弹性参数之间的线性求解关系，利用四个角度的弹性阻抗数据体提取流体体积模量。

得到绝对弹性阻抗数据体之后需要基于阻抗与弹性参数之间的关系提取相应的弹性参数。弹性参数与阻抗之间是非线性指数关系，为了简化求解方法，对公式进行对数变换。考虑到入射角度相同时，各采样点所对应的系数是相同的，其不随时间发生变化，将纵横波速度与密度之前的系数分别用a(θ)、b(θ)、c(θ)和d(θ)表示，可以将公式(14)表示为：

考虑同一道地震数据的不同采样点(假设有n个采样点)，可以借助井旁道与尺度匹配之后的测井曲线构建如下矩阵：

$\left[\begin{array}{cccc}\ln \frac{K_f\left(t_1\right)}{K_{f0}}& \ln \frac{f_m\left(t_1\right)}{f_{m0}}& \ln \frac{\mathrm{\φ}\left(t_1\right)}{{\mathrm{\φ}}_0}& \ln \frac{\mathrm{\ρ}\left(t_1\right)}{{\mathrm{\ρ}}_0}\\ \ln \frac{K_f\left(t_2\right)}{K_{f0}}& \ln \frac{f_m\left(t_2\right)}{f_{m0}}& \ln \frac{\mathrm{\φ}\left(t_2\right)}{{\mathrm{\φ}}_0}& \ln \frac{\mathrm{\ρ}\left(t_2\right)}{{\mathrm{\ρ}}_0}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ \ln \frac{K_f\left(t_n\right)}{K_{f0}}& \ln \frac{f_m\left(t_n\right)}{f_{m0}}& \ln \frac{\mathrm{\φ}\left(t_n\right)}{{\mathrm{\φ}}_0}& \ln \frac{\mathrm{\ρ}\left(t_n\right)}{{\mathrm{\ρ}}_0}\end{array}\right]\left[\begin{array}{c}a\left(\mathrm{\θ}\right)\\ b\left(\mathrm{\θ}\right)\\ c\left(\mathrm{\θ}\right)\\ d\left(\mathrm{\θ}\right)\end{array}\right]=\left[\begin{array}{c}\ln \frac{\mathrm{EI}(t_1,\mathrm{\θ})}{{\mathrm{EI}}_0}\\ \ln \frac{\mathrm{EI}(t_2,\mathrm{\θ})}{{\mathrm{EI}}_0}\\ .\\ .\\ .\\ \ln \frac{\mathrm{EI}(t_n,\mathrm{\θ})}{{\mathrm{EI}}_0}\end{array}\right]---\left(16\right)$

求解以上矩阵，即可获得适用于研究工区的系数a(θ)、b(θ)、c(θ)和d(θ)的数值，针对不同的入射角度分别计算，便可得到四个入射角度(分别为θ₁、θ₂、θ₃和θ₄)对应的16个系数(即a(θ₁)、b(θ₁)、c(θ₁)、d(θ₁)；a(θ₂)、b(θ₂)、c(θ₂)、d(θ₂)；a(θ₃)、b(θ₃)、c(θ₃)、d(θ₃)；a(θ₄)、b(θ₄)、c(θ₄)、d(θ₄))。将16个系数分别代入方程组(15)，即可得到如下方程组：

通过求解以上方程组即可得到任意采样点处的流体体积模量数据体、固体项数据体、密度参数数据体以及孔隙度数据体。

S414：根据所述的流体体积模量数据体识别当前储层孔隙中流体的类型。提取的流体体积模量K_f数据体与储层流体体积模量的理论值进行对比，可以识别储层孔隙中流体的类型。

如上所述，即为本发明提供的一种基于双相介质固液解耦的叠前地震反演方法，该方法的实现过程需要首先进行叠前地震资料、测井资料以及地质成果资料的预处理，经过叠前道集的角度叠加，测井资料的预处理和地震层位精细解释之后，然后进行角度子波提取和约束稀疏脉冲反演以得到相对阻抗数据体，基于测井资料和地震层位数据构建的低频阻抗数据对相对数据进行补偿，最终得到绝对弹性阻抗数据体，通过弹性参数提取从中直接得到流体体积模量K_f数据体。常规流体因子岩石物理意义模糊且普适性较差，而且受岩石双相特性的耦合效应影响，在表征岩石流体效应的时候易受岩石固体效应(特别是孔隙度参数)影响，造成流体识别假象。流体体积模量K_f只与岩石孔隙流体有关，不受其他岩石固体效应影响，作为一项对储层孔隙流体类型敏感的弹性参数在流体识别领域的效果更佳。利用叠前弹性阻抗反演实现了流体因子的直接提取，可以避免引入累计误差，提高了储层流体指示作用的可靠性。

图10为本发明提出的一种基于双相介质固液解耦的叠前地震反演系统的实施方式一的具体结构框图，由图10可知，在实施方式一中，所述的系统包括：

资料采集装置100，用于采集目的储层的叠前地震资料、测井资料以及地质成果资料。

叠前地震资料的采集在具体的实施方式中可通过如下方式：在勘探工区中，根据目的层的深度和地质特点，考虑动校正拉伸、干扰波、多次波等因素来设计宽方位角的三维观测系统，保证有足够的偏移距、方位角。经过激发、接收，得到可以满足AVO分析需求的宽方位大角度的叠前地震资料。

测井资料的采集在具体的实施方式中可通过如下方式：在勘探工区中进行全波列测井，得到测井资料，主要包含纵波速度、横波速度、密度等全波列测井曲线，孔隙度等解释成果曲线，测井层位，录井资料等。

固液解耦近似方程构建装置200，用于构建固液解耦近似方程。

AVO(Amplitude Versus Offset，振幅随偏移距的变化)近似公式Aki公式是进行弹性参数估算的基础，本发明基于多孔弹性介质理论以及Nur岩石物理模型，推导得到的流体体积模量分离的AVO近似公式即固液解耦近似方程，其为基于地震资料的流体体积模量参数估算提供了理论基础。

Russell等人从Aki-Richards近似公式出发，推导了突出储层流体特征的反射特征近似公式：

$R_{\mathrm{PP}}\left(\mathrm{\θ}\right)=\left[(1-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2})\frac{{\sec }^2\mathrm{\θ}}{4}\right]\frac{\mathrm{\Δf}}{f}+[\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{4{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}]\frac{\mathrm{\Δ\μ}}{\mathrm{\μ}}+[\frac{1}{2}-\frac{{\sec }^2\mathrm{\θ}}{4}]\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}---\left(1\right)$

Dehua Han等人提出了Gassmann流体项的经验公式：

f＝G(φ)K_f       (2)

其中， $G\left(\mathrm{\φ}\right)=\frac{{(1-K_n)}^2}{\mathrm{\φ}},K_n=\frac{K_{\mathrm{dry}}}{K_m}.$

将公式(2)代入公式(1)，并考虑到剪切模量不受孔隙流体的影响，利用干岩石剪切模量μ_dry替换μ，进行相应变换，可以得到：

$R_{\mathrm{PP}}\left(\mathrm{\θ}\right)=\left[(1-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2})\frac{{\sec }^2\mathrm{\θ}}{4}\right]\frac{\mathrm{\Δ}\left(G\left(\mathrm{\φ}\right)K_f\right)}{G\left(\mathrm{\φ}\right)K_f}+[\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{4{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}]\frac{\mathrm{\Δ\μ}}{\mathrm{\μ}}+[\frac{1}{2}-\frac{{\sec }^2\mathrm{\θ}}{4}]\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}---\left(3\right)$

Nur提出的临界孔隙度模型，表达式如下式所示：

$\left\{\begin{array}{c}{K}_{\mathrm{dry}}=K_m(1-\frac{\mathrm{\φ}}{{\mathrm{\φ}}_c})\\ {\mathrm{\μ}}_{\mathrm{dry}}={\mathrm{\μ}}_m(1-\frac{\mathrm{\φ}}{{\mathrm{\φ}}_c})\end{array}\right.---\left(4\right)$

其中，φ_c表示临界孔隙度，K_dry表示干岩石的体积模量，μ_dry表示干岩石的剪切模量，K_m表示固体矿物基质的体积模量，μ_m表示矿物基质的剪切模量。

结合Nur公式(4)对公式(3)进一步展开可以得到：

$\begin{array}{c}{R}_{\mathrm{PP}}\left(\mathrm{\θ}\right)=\left[(1-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2})\frac{{\sec }^2\mathrm{\θ}}{4}\right](\frac{\mathrm{\ΔG}\left(\mathrm{\φ}\right)}{G\left(\mathrm{\φ}\right)}+\frac{\mathrm{\Δ}{K}_f}{K_f})+[\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{4{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}](\frac{\mathrm{\Δ}{\mathrm{\μ}}_m}{{\mathrm{\μ}}_m}+\frac{\mathrm{\Δ}({\mathrm{\φ}}_c-\mathrm{\φ})}{{\mathrm{\φ}}_c-\mathrm{\φ}})\\ +[\frac{1}{2}-\frac{{\sec }^2\mathrm{\θ}}{4}]\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}\end{array}---\left(5\right)$

将代入G(φ)，进一步展开可以得到：

$\begin{array}{c}{R}_{\mathrm{PP}}\left(\mathrm{\θ}\right)=\left[(1-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2})\frac{{\sec }^2\mathrm{\θ}}{4}\right]+\frac{\mathrm{\Δ}{K}_f}{K_f}+[\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{4{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}]\frac{\mathrm{\Δ}{\mathrm{\μ}}_m}{{\mathrm{\μ}}_m}+[\frac{{\sec }^2\mathrm{\θ}}{4}-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{2{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}+\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}]\frac{\mathrm{\Δ\φ}}{\mathrm{\φ}}\\ +[\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{4{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}]\frac{\mathrm{\Δ\φ}({\mathrm{\φ}}_c-\mathrm{\φ})}{\mathrm{\φ}({\mathrm{\φ}}_c-\mathrm{\φ})}+[\frac{1}{2}-\frac{{\sec }^2\mathrm{\θ}}{4}]\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}\end{array}---\left(6\right)$

令f_m＝φμ，并利用公式(4)，最终得到的固液解耦近似方程的近似公式如下所示：

$\begin{array}{c}{R}_{\mathrm{PP}}\left(\mathrm{\θ}\right)=\left[(1-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2})\frac{{\sec }^2\mathrm{\θ}}{4}\right]\frac{\mathrm{\Δ}{K}_f}{K_f}+[\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{4{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}]\frac{\mathrm{\Δ}\left(f_m\right)}{f_m}+[\frac{1}{2}-\frac{{\sec }^2\mathrm{\θ}}{4}]\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}\\ +(\frac{{\sec }^2\mathrm{\θ}}{4}-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{2{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}+\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ})\frac{\mathrm{\Δ\φ}}{\mathrm{\φ}}\end{array}---\left(7\right)$

其中，f_m＝φμ为地下岩石的固体项，θ为入射角度，K_f为地下岩石的流体体积模量，f_m为地下岩石的固体项，φ为地下岩石的孔隙度，ρ为地下岩石的密度，μ为地下岩石的剪切模量，△K_f为界面两侧的流体体积模量的差值，△f_m为界面两侧固体项的差值，△ρ为界面两侧的密度的差值，△φ为界面两侧的孔隙度的差值，为干岩石纵横波速度比的平方，为饱和岩石纵横波速度比的平方。

子波提取装置300，用于根据所述的叠前地震资料、测井资料、地质成果资料以及所述的固液解耦近似方程提取子波。振幅随偏移距的变化导致了子波随偏移距的空变特性，为消除不同角度道集间的相互影响，在此需要针对不同的入射角度对部分角度叠加数据体分别提取不同的子波。

弹性阻抗数据体确定装置400，用于根据所述的叠前地震资料、测井资料、地质成果资料以及子波确定绝对弹性阻抗数据体。综合叠前地震数据、测井数据、层位数据和子波，利用约束稀疏脉冲反演算法计算相对弹性阻抗数据体，然后通过补偿低频阻抗成分，最终得到绝对弹性阻抗数据体。

流体体积模量确定装置500，用于根据所述的绝对弹性阻抗数据体确定流体体积模量数据体。根据弹性阻抗与弹性参数之间的线性求解关系，利用四个角度的弹性阻抗数据体提取流体体积模量。

流体类型识别装置600，用于根据所述的流体体积模量数据体识别当前储层孔隙中流体的类型。提取的流体体积模量K_f数据体与储层流体体积模量的理论值进行对比，可以识别储层孔隙中流体的类型。

图11为本发明实施例提供的一种基于双相介质固液解耦的叠前地震反演系统的实施方式二的结构框图，由图11可知，在实施方式二中，该系统具体包括：

资料采集装置100，用于采集目的储层的叠前地震资料、测井资料以及地质成果资料。

叠前地震资料的采集在具体的实施方式中可通过如下方式：在勘探工区中，根据目的层的深度和地质特点，考虑动校正拉伸、干扰波、多次波等因素来设计宽方位角的三维观测系统，保证有足够的偏移距、方位角。经过激发、接收，得到可以满足AVO分析需求的宽方位大角度的叠前地震资料。

测井资料的采集在具体的实施方式中可通过如下方式：在勘探工区中进行全波列测井，得到测井资料。

保幅处理装置700，用于对所述的叠前地震资料进行保幅处理；

去噪处理装置800，用于对经过保幅处理的叠前地震资料进行去噪处理；

角度叠加装置900，用于根据目的储层对经过去噪处理的叠前地震资料进行角度叠加，生成四个入射角度的角度叠加道集。

保幅处理装置700至角度叠加装置900是反演之前需要对叠前地震资料进行的预处理。在保证地震资料品质的前提下根据目的储层进行相应部分角度叠加，最终生成四个角度的部分角度叠加道集。此处提及的四个角度，在具体的实施方式中，若叠前地震数据的角度范围大、大角度数据信噪比较高，可选择6度、14度、22度、30度；若叠前地震数据的角度范围不大、大角度数据信噪比较低，可选择5度、10度、15度、20度。

固液解耦近似方程构建装置200，用于构建固液解耦近似方程。

子波提取装置300，用于根据所述的叠前地震资料、测井资料、地质成果资料以及所述的固液解耦近似方程提取子波。振幅随偏移距的变化导致了子波随偏移距的空变特性，为消除不同角度道集间的相互影响，在此需要针对不同的入射角度对部分角度叠加数据体分别提取不同的子波。

弹性阻抗数据体确定装置400，用于根据所述的叠前地震资料、测井资料、地质成果资料以及子波确定绝对弹性阻抗数据体。综合叠前地震数据、测井数据、层位数据和子波，利用约束稀疏脉冲反演算法计算相对弹性阻抗数据体，然后通过补偿低频阻抗成分，最终得到绝对弹性阻抗数据体。

流体体积模量确定装置500，用于根据所述的绝对弹性阻抗数据体确定流体体积模量数据体。根据弹性阻抗与弹性参数之间的线性求解关系，利用四个角度的弹性阻抗数据体提取流体体积模量。

流体类型识别装置600，用于根据所述的流体体积模量数据体识别当前储层孔隙中流体的类型。提取的流体体积模量K_f与储层流体体积模量的理论值进行对比，可以识别储层孔隙中流体的类型。

图12为本发明实施例提供的一种基于双相介质固液解耦的叠前地震反演系统的实施方式三的结构框图，由图12可知，在实施方式三中，该系统具体包括：

资料采集装置100，用于采集目的储层的叠前地震资料、测井资料以及地质成果资料。

保幅处理装置700，用于对所述的叠前地震资料进行保幅处理；

去噪处理装置800，用于对经过保幅处理的叠前地震资料进行去噪处理；

角度叠加装置900，用于根据目的储层对经过去噪处理的叠前地震资料进行角度叠加，生成四个入射角度的角度叠加道集。

保幅处理装置700至角度叠加装置900是反演之前需要对叠前地震资料进行的预处理。在保证地震资料品质的前提下根据目的储层进行相应部分角度叠加，最终生成四个角度的部分角度叠加道集。此处提及的四个角度，在具体的实施方式中，若叠前地震数据的角度范围大、大角度数据信噪比较高，可选择6度、14度、22度、30度；若叠前地震数据的角度范围不大、大角度数据信噪比较低，可选择5度、10度、15度、20度。

环境校正装置1000，用于对所述的测井资料进行环境校正；

奇异值消除装置1100，用于对经过环境校正后的测井资料进行奇异值消除；

纵横速度比确定装置1200，用于根据经过奇异值消除后的测井资料确定所述目的储层的干岩石纵横波速度比、饱和岩石纵横波速度比。具体的实施方式中，在勘探工区中进行全波列测井，得到的测井资料需要经过环境校正、奇异值消除等预处理后，结合岩石物理经验公式或者实验统计结果计算干岩石纵横波速度比γ_dry、饱和岩石纵横波速度比γ_sat的数值。

剪切模量曲线确定装置1300，用于根据所述的目的储层的干岩样纵横波速度比以及所述的测井资料确定与所述测井资料对应的流体体积模量曲线以及剪切模量曲线。在具体的实施方式中，经过换算得到相应的流体体积模量为K_f，剪切模量为μ曲线，后续用来计算井旁道的基于固液解耦近似方程的弹性阻抗伪测井曲线。

固液解耦近似方程构建装置200，用于构建固液解耦近似方程。

子波提取装置300，用于根据所述的叠前地震资料、测井资料、地质成果资料以及所述的固液解耦近似方程提取子波。振幅随偏移距的变化导致了子波随偏移距的空变特性，为消除不同角度道集间的相互影响，在此需要针对不同的入射角度对部分角度叠加数据体分别提取不同的子波。

弹性阻抗数据体确定装置400，用于根据所述的叠前地震资料、测井资料、地质成果资料以及子波确定绝对弹性阻抗数据体。综合叠前地震数据、测井数据、层位数据和子波，利用约束稀疏脉冲反演算法计算相对弹性阻抗数据体，然后通过补偿低频阻抗成分，最终得到绝对弹性阻抗数据体。

流体体积模量确定装置500，用于根据所述的绝对弹性阻抗数据体确定流体体积模量数据体。根据弹性阻抗与弹性参数之间的线性求解关系，利用四个角度的弹性阻抗数据体提取流体体积模量。

流体类型识别装置600，用于根据所述的流体体积模量数据体识别当前储层孔隙中流体的类型。提取的流体体积模量K_f与储层流体体积模量的理论值进行对比，可以识别储层孔隙中流体的类型。

图13为本发明实施例提供的一种基于双相介质固液解耦的叠前地震反演系统的实施方式四的结构框图，由图13可知，在实施方式四中，该系统具体包括：

资料采集装置100，用于采集目的储层的叠前地震资料、测井资料以及地质成果资料。

保幅处理装置700，用于对所述的叠前地震资料进行保幅处理；

去噪处理装置800，用于对经过保幅处理的叠前地震资料进行去噪处理；

角度叠加装置900，用于根据目的储层对经过去噪处理的叠前地震资料进行角度叠加，生成四个入射角度的角度叠加道集。

保幅处理装置700至角度叠加装置900是反演之前需要对叠前地震资料进行的预处理。在保证地震资料品质的前提下根据目的储层进行相应部分角度叠加，最终生成四个角度的部分角度叠加道集。此处提及的四个角度，在具体的实施方式中，若叠前地震数据的角度范围大、大角度数据信噪比较高，可选择6度、14度、22度、30度；若叠前地震数据的角度范围不大、大角度数据信噪比较低，可选择5度、10度、15度、20度。

环境校正装置1000，用于对所述的测井资料进行环境校正；

奇异值消除装置1100，用于对经过环境校正后的测井资料进行奇异值消除；

纵横速度比确定装置1200，用于根据经过奇异值消除后的测井资料确定所述目的储层的干岩石纵横波速度比、饱和岩石纵横波速度比。具体的实施方式中，在勘探工区中进行全波列测井，得到的测井资料需要经过环境校正、奇异值消除等预处理后，结合岩石物理经验公式或者实验统计结果计算干岩石纵横波速度比γ_dry、饱和岩石纵横波速度比γ_sat的数值。

剪切模量曲线确定装置1300，用于根据所述的目的储层的干岩样纵横波速度比以及所述的测井资料确定与所述测井资料对应的流体体积模量曲线以及剪切模量曲线。在具体的实施方式中，经过换算得到相应的流体体积模量为K_f，剪切模量为μ曲线，后续用来计算井旁道的基于固液解耦近似方程的弹性阻抗伪测井曲线。

地质解释装置1400，用于根据所述的叠前地震资料、测井资料、地质成果资料对目的储层进行精细层位追踪，得到地震层位数据。

在具体的实施方式中，结合工区的地质背景，在保证目的层位符合工区地质沉积模式与层序格架的基础上，对研究的目的层进行精细地质解释，解释得到的地震层位数据不仅用来标识时深关系，还可结合弹性阻抗伪测井曲线建立低频阻抗数据体。

固液解耦近似方程构建装置200，用于构建固液解耦近似方程。

子波提取装置300，用于根据所述的叠前地震资料、测井资料、地质成果资料以及所述的固液解耦近似方程提取子波。振幅随偏移距的变化导致了子波随偏移距的空变特性，为消除不同角度道集间的相互影响，在此需要针对不同的入射角度对部分角度叠加数据体分别提取不同的子波。图14为本发明实施例提供的一种基于双相介质固液解耦的叠前地震反演系统中的子波提取装置300的具体结构框图，由图14可知，子波提取装置具体包括：

孔隙度提取模块301，用于从所述的测井资料中提取出地下岩石的密度以及孔隙度。在具体的实施方式中，地下岩石的密度通过ρ来表示，地下岩石的孔隙度通过φ来表示。

固体项确定模块302，用于根据所述的剪切模量曲线以及所述的孔隙度确定地下岩石的固体项，在具体的实施方式中，地下岩石的固体项通过f_m来表示，则f_m＝φμ。

弹性阻抗方程推导模块303，用于基于所述的固液解耦近似方程推导出弹性阻抗方程。

与普通弹性阻抗方程推导过程类似，固液解耦近似方程的反射系数可以表示为：

$R_{\mathrm{pp}}\left(\mathrm{\θ}\right)\≈\frac{1}{2}\frac{\mathrm{\ΔEI}}{\mathrm{EI}}\≈\frac{1}{2}\mathrm{\Δ}\ln \left(\mathrm{EI}\right)---\left(8\right)$

将公式(7)代入上式，得到：

$\begin{array}{c}\frac{1}{2}\mathrm{\Δ}\ln \left(\mathrm{EI}\right)=\left[(1-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2})\frac{{\sec }^2\mathrm{\θ}}{4}\right]\frac{\mathrm{\Δ}{K}_f}{K_f}+[\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{4{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}]\frac{\mathrm{\Δ}\left(f_m\right)}{f_m}+[\frac{1}{2}-\frac{{\sec }^2\mathrm{\θ}}{4}]\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}\\ +(\frac{{\sec }^2\mathrm{\θ}}{4}-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{2{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}+\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ})\frac{\mathrm{\Δ\φ}}{\mathrm{\φ}}\end{array}---\left(9\right)$

将弹性参数的相对变化量用对数形式表示，可以得到：

$\begin{array}{c}\mathrm{\Δ}\ln \left(\mathrm{EI}\right)=\left[(1-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2})\frac{{\sec }^2\mathrm{\θ}}{2}\right]\mathrm{\Δ}\ln \left(K_f\right)+[\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{2{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{4}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}]\mathrm{\Δ}\ln \left(f_m\right)+[1-\frac{{\sec }^2\mathrm{\θ}}{2}]\mathrm{\Δ}\ln \left(\mathrm{\ρ}\right)\\ +(\frac{{\sec }^2\mathrm{\θ}}{2}-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}+\frac{4}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ})\mathrm{\Δ}\ln \left(\mathrm{\φ}\right)\end{array}---\left(10\right)$

进一步变形，可以得到：

$\mathrm{\Δ}\ln \left(\mathrm{EI}\right)=\mathrm{\Δ}\ln \left\{{K_f}^{\left[(\frac{1}{2}-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{2{\mathrm{\γ}}_{\mathrm{sat}}^2}){\sec }^2\mathrm{\θ}\right]}{f_m}^{[\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{2{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{4}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}]}{\mathrm{\ρ}}^{[1-\frac{{\sec }^2\mathrm{\θ}}{2}]}{\mathrm{\φ}}^{(\frac{{\sec }^2\mathrm{\θ}}{2}-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}+\frac{4}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ})}\right\}---\left(11\right)$

两边取积分并将其指数化，消掉等式两边的微分项和对数项，进一步取积分常数为0，得到：

$\left(\mathrm{EI}\right)={K_f}^{\left[(\frac{1}{2}-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{2{\mathrm{\γ}}_{\mathrm{sat}}^2}){\sec }^2\mathrm{\θ}\right]}{f_m}^{[\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{2{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{4}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}]}{\mathrm{\ρ}}^{[1-\frac{{\sec }^2\mathrm{\θ}}{2}]}{\mathrm{\φ}}^{(\frac{{\sec }^2\mathrm{\θ}}{2}-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}+\frac{4}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ})}---\left(12\right)$

即：

$\mathrm{EI}\left(\mathrm{\θ}\right)={K_f}^a{f_m}^b{\mathrm{\φ}}^c{\mathrm{\ρ}}^d---\left(13\right)$

AVO近似公式是叠前反演的基础，本发明以包含流体体积模量K_f的固液解耦近似方程为基础，推导标准化的弹性阻抗方程，消除了不同角度之间存在的量纲不统一问题。通过叠前弹性阻抗反演得到相对阻抗数据体，经过低频分量补偿得到绝对弹性阻抗数据体，最终利用弹性参数提取直接得到反映孔隙流体效应的流体因子K_f，从而实现敏感流体因子(即流体体积模量K_f)的叠前地震反演方法。

对公式(13)进行标准化可以推导出的弹性阻抗方程即为：

$\mathrm{EI}\left(\mathrm{\θ}\right)={\mathrm{EI}}_0{\left(\frac{K_f}{K_{f0}}\right)}^a{\left(\frac{f_m}{f_{m0}}\right)}^b{\left(\frac{\mathrm{\φ}}{{\mathrm{\φ}}_0}\right)}^c{\left(\frac{\mathrm{\ρ}}{{\mathrm{\ρ}}_0}\right)}^d---\left(14\right)$

其中， $a=(1-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2})\frac{{\sec }^2\mathrm{\θ}}{2},b=\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{2{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{4}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ},c=\frac{{\sec }^2\mathrm{\θ}}{2}-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}+\frac{4}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ},$ $d=1-\frac{{\sec }^2\mathrm{\θ}}{2};{\mathrm{EI}}_0={\left(K_f\right)}^{-\frac{1}{2}}{f}_m{\mathrm{\φ}}^{-\frac{3}{2}}{\mathrm{\ρ}}^{\frac{1}{2}}{\mathrm{\γ}}_{\mathrm{sat}}{({\mathrm{\γ}}_{\mathrm{sat}}^2-{\mathrm{\γ}}_{\mathrm{dry}}^2)}^{\frac{1}{2}},$ K_f、f_m、ρ和φ分别表示地下岩石的流体体积模量，固体项，密度和孔隙度，θ表示入射角度，表示干岩石纵横波速度比的平方，表示饱和岩石纵横波速度比的平方，K_f0、f_m0、ρ₀和φ₀分别为K_f、f_m、ρ和φ的平均值。

弹性阻抗伪测井曲线确定模块304，用于根据所述的入射角度、地下岩石的密度、固体项、孔隙度、流体体积模量曲线以及弹性阻抗方程确定弹性阻抗伪测井曲线。在具体的实施方式中，对叠前地震资料预处理后得到四个入射角度，将每个入射角度、地下岩石的密度、固体项、孔隙度、流体体积模量曲线带入弹性阻抗方程中，皆能确定出一条弹性阻抗伪测井曲线，因此四个入射角度即能确定出四条弹性阻抗伪测井曲线。在反演中，弹性阻抗伪测井曲线不仅可以用来作为必备的先验约束条件，还可以以此结合地质解释层位建立低频模型来补偿地震波传播过程中损失的低频成分。

子波提取模块305，用于针对不同的入射角度，根据入射角度对应的角度叠加道集、弹性阻抗伪测井曲线、地震层位数据提取不同的子波。在具体的实施方式中，振幅随偏移距的变化导致了子波随偏移距的空变特性，为消除不同角度道集间的相互影响，在此需要针对不同的入射角度对部分角度叠加道集分别提取不同的子波。

由图13可知，在实施方式四中，该系统还包括：

弹性阻抗数据体确定装置400，用于根据所述的叠前地震资料、测井资料、地质成果资料以及子波确定绝对弹性阻抗数据体。综合叠前地震数据、测井数据、层位数据和子波，利用约束稀疏脉冲反演算法计算相对弹性阻抗数据体，然后通过补偿低频阻抗成分，最终得到绝对弹性阻抗数据体。图15为本发明实施例提供的一种基于双相介质固液解耦的叠前地震反演系统中的弹性阻抗数据体确定装置400的具体结构框图，由图15可知，弹性阻抗数据体确定装置具体包括：

相对弹性阻抗数据体确定模块401，用于根据所述的叠前地震资料、测井资料、地质成果资料以及子波确定出相对弹性阻抗数据体。在具体的实施方式中，综合叠前地震资料、测井资料、地震层位数据和子波，利用约束稀疏脉冲反演算法计算相对弹性阻抗数据体。

低频阻抗数据体构建模块402，用于根据所述的弹性阻抗伪测井曲线以及地震层位数据构建低频阻抗数据体。在反演中，弹性阻抗伪测井曲线不仅可以用来作为必备的先验约束条件，还可以以此结合地震层位数据建立低频阻抗数据体来补偿地震波传播过程中损失的低频成分。

补偿模块403，用于根据所述的低频阻抗数据体对所述的弹性阻抗数据体进行补偿，得到绝对弹性阻抗数据体。即通过补偿低频阻抗成分，最终得到绝对弹性阻抗数据体。

由图13可知，在实施方式四中，该系统还包括：

流体体积模量确定装置500，用于根据所述的绝对弹性阻抗数据体确定流体体积模量数据体。根据弹性阻抗与弹性参数之间的线性求解关系，利用四个角度的弹性阻抗数据体提取流体体积模量。

得到绝对弹性阻抗数据体之后需要基于阻抗与弹性参数之间的关系提取相应的弹性参数。弹性参数与阻抗之间是非线性指数关系，为了简化求解方法，对公式进行对数变换。考虑到入射角度相同时，各采样点所对应的系数是相同的，其不随时间发生变化，将纵横波速度与密度之前的系数分别用a(θ)、b(θ)、c(θ)和d(θ)表示，可以将公式(14)表示为：

考虑同一道地震数据的不同采样点(假设有n个采样点)，可以借助井旁道与尺度匹配之后的测井曲线构建如下矩阵：

$\left[\begin{array}{cccc}\ln \frac{K_f\left(t_1\right)}{K_{f0}}& \ln \frac{f_m\left(t_1\right)}{f_{m0}}& \ln \frac{\mathrm{\φ}\left(t_1\right)}{{\mathrm{\φ}}_0}& \ln \frac{\mathrm{\ρ}\left(t_1\right)}{{\mathrm{\ρ}}_0}\\ \ln \frac{K_f\left(t_2\right)}{K_{f0}}& \ln \frac{f_m\left(t_2\right)}{f_{m0}}& \ln \frac{\mathrm{\φ}\left(t_2\right)}{{\mathrm{\φ}}_0}& \ln \frac{\mathrm{\ρ}\left(t_2\right)}{{\mathrm{\ρ}}_0}\\ .& .& .& .\\ .& .& .& .\\ .& .& .& .\\ \ln \frac{K_f\left(t_n\right)}{K_{f0}}& \ln \frac{f_m\left(t_n\right)}{f_{m0}}& \ln \frac{\mathrm{\φ}\left(t_n\right)}{{\mathrm{\φ}}_0}& \ln \frac{\mathrm{\ρ}\left(t_n\right)}{{\mathrm{\ρ}}_0}\end{array}\right]\left[\begin{array}{c}a\left(\mathrm{\θ}\right)\\ b\left(\mathrm{\θ}\right)\\ c\left(\mathrm{\θ}\right)\\ d\left(\mathrm{\θ}\right)\end{array}\right]=\left[\begin{array}{c}\ln \frac{\mathrm{EI}(t_1,\mathrm{\θ})}{{\mathrm{EI}}_0}\\ \ln \frac{\mathrm{EI}(t_2,\mathrm{\θ})}{{\mathrm{EI}}_0}\\ .\\ .\\ .\\ \ln \frac{\mathrm{EI}(t_n,\mathrm{\θ})}{{\mathrm{EI}}_0}\end{array}\right]---\left(16\right)$

求解以上矩阵，即可获得适用于研究工区的系数a(θ)、b(θ)、c(θ)和d(θ)的数值，针对不同的入射角度分别计算，便可得到四个入射角度(分别为θ₁、θ₂、θ₃和θ₄)对应的16个系数(即a(θ₁)、b(θ₁)、c(θ₁)、d(θ₁)；a(θ₂)、b(θ₂)、c(θ₂)、d(θ₂)；a(θ₃)、b(θ₃)、c(θ₃)、d(θ₃)；a(θ₄)、b(θ₄)、c(θ₄)、d(θ₄))。将16个系数分别代入方程组(15)，即可得到如下方程组：

通过求解以上方程组即可得到任意采样点处的流体体积模量、固体项、密度参数以及孔隙度参数。

流体类型识别装置600，用于根据所述的流体体积模量数据体识别当前储层孔隙中流体的类型。提取的流体体积模量K_f与储层流体体积模量的理论值进行对比，可以识别储层孔隙中流体的类型。

如上所述，即为本发明提供的一种基于双相介质固液解耦的叠前地震反演系统，该系统的实现过程需要首先进行叠前地震资料、测井资料以及地质成果资料的预处理，经过叠前道集的角度叠加，测井资料的预处理和地质层位精细解释之后，然后进行角度子波提取和约束稀疏脉冲反演以得到相对阻抗数据体，基于测井资料和地震层位数据构建的低频阻抗数据对相对数据进行补偿，最终得到绝对弹性阻抗数据体，通过弹性参数提取从中直接得到流体体积模量K_f数据体。常规流体因子岩石物理意义模糊且普适性较差，而且受岩石双相特性的耦合效应影响，在表征岩石流体效应的时候易受岩石固体效应(特别是孔隙度参数)影响，造成流体识别假象。流体体积模量K_f只与岩石孔隙流体有关，不受其他岩石固体效应影响，作为一项对储层孔隙流体类型敏感的弹性参数在流体识别领域的效果更佳。利用叠前弹性阻抗反演实现了流体因子的直接提取，可以避免引入累计误差，提高了储层流体指示作用的可靠性。

下面结合具体的实施例，详细介绍本发明的技术方案。

图1是利用某工区实际数据计算得到的常用流体因子的流体指示系数比较图。为了说明流体体积模量在储层流体识别中的优越性，在此利用某工区的实际数据，利用基于固液解耦近似方程的流体因子叠前地震反演方法得到流体体积模量参数，将其与其他几种常用流体因子(图1中，1为μρ,2为I_p,3为4为f，5为K_f)进行流体敏感性比较，如表1所示为流体指示敏感性统计结果，通过比较流体指示系数(如图1)发现，流体体积模量K_f对流体的区分能力最强，说明K_f作为一项对储层孔隙流体类型敏感的弹性参数在流体识别领域的效果更佳，分析其原因是K_f只与岩石孔隙流体有关，不受其他岩石固体效应影响，而常规流体因子受岩石双相特性的耦合效应影响，在表征岩石流体效应的时候易受岩石固体效应(特别是孔隙度参数)影响，造成流体识别假象。

表1

图2是本实施例种基于固液解耦近似方程的流体因子叠前地震反演流程图。

图3是根据上述方法，对某工区实际地震数据进行了基于固液解耦近似方程的流体因子叠前地震反演方法，得到流体体积模量数据体K_f剖面，以此作为流体因子进行储层流体识别。根据反演得到的流体体积模量对储层含流体类型进行了识别。图中测井曲线是测井解释柱状图，其中三角形代表含气层，菱形代表含水层，从图中可以看到流体体积模量可以较好的体现储层的含流体信息，含气层主要表现为浅色区域，含水层则表现为较深色区域，其与实际钻井资料基本吻合，说明该方法可以较好的消除储层孔隙度引起的流体识别假象，在储层流体精细识别上具有一定的优越性。

综上所述，本发明提出的一种基于双相介质固液解耦的叠前地震反演方法及系统，在石油工业中，利用地下岩层中传播的地震波特点进行储层预测是描述和预测地下构造形态和含油气性的主要技术手段。相比较叠后地震资料，叠前地震资料中包含有更为丰富的动力学与运动学信息，在岩石物理理论的指导下，利用叠前地震反演技术可以从中提取对储层孔隙流体类型敏感的弹性参数，以此作为流体因子进行相应的储层特征描述与流体类型判识。本发明以实现岩石固液双相效应解耦并反演完全体现流体效应的流体体积模量参数K_f为出发点，利用叠前弹性阻抗反演实现了流体因子的直接提取。本发明计算的流体因子在具有更高的流体指示灵敏性的同时，有效的提高了储层流体识别和储层特征描述的可靠性。

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

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

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

Claims (16)

1.一种基于双相介质固液解耦的叠前地震反演方法，其特征是，所述的方法具体包括：

采集目的储层的叠前地震资料、测井资料以及地质成果资料；

构建固液解耦近似方程；

根据所述的叠前地震资料、测井资料、地质成果资料以及所述的固液解耦近似方程提取子波；

根据所述的叠前地震资料、测井资料、地质成果资料以及子波确定绝对弹性阻抗数据体；

根据所述的绝对弹性阻抗数据体确定流体体积模量数据体；

根据所述的流体体积模量数据体识别当前储层孔隙中流体的类型。

2.根据权利要求1所述的方法，其特征是，所述的方法在构建固液解耦近似方程之前还包括：

对所述的叠前地震资料进行保幅处理；

对经过保幅处理的叠前地震资料进行去噪处理；

根据目的储层对经过去噪处理的叠前地震资料进行角度叠加，生成四个入射角度的角度叠加道集。

3.根据权利要求2所述的方法，其特征是，所述的方法在构建固液解耦近似方程之前还包括：

对所述的测井资料进行环境校正；

对经过环境校正后的测井资料进行奇异值消除；

根据经过奇异值消除后的测井资料确定所述目的储层的干岩石纵横波速度比、饱和岩石纵横波速度比；

根据所述的目的储层的干岩样纵横波速度比以及所述的测井资料确定与所述测井资料对应的流体体积模量曲线以及剪切模量曲线。

4.根据权利要求3所述的方法，其特征是，所述的方法在构建固液解耦近似方程之前还包括：

根据所述的叠前地震资料、测井资料、地质成果资料对目的储层进行精细层位追踪，得到地震层位数据。

5.根据权利要求1所述的方法，其特征是，所述的固液解耦近似方程为：

$\begin{array}{c}{R}_{\mathrm{PP}}\left(\mathrm{\θ}\right)=\left[(1-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2})\frac{{\sec }^2\mathrm{\θ}}{4}\right]\frac{{\mathrm{\ΔK}}_f}{K_f}+[\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{4\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}]\frac{\mathrm{\Δ}\left(f_m\right)}{f_m}[\frac{1}{2}-\frac{{\sec }^2\mathrm{\θ}}{4}]\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}\\ +(\frac{{\sec }^2\mathrm{\θ}}{4}-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{2{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}+\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ})\frac{\mathrm{\Δ\φ}}{\mathrm{\φ}}\end{array}$

其中，f_m＝φμ为地下岩石的固体项，θ为入射角度，K_f为地下岩石的流体体积模量，f_m为地下岩石的固体项，φ为地下岩石的孔隙度，ρ为地下岩石的密度，μ为地下岩石的剪切模量，△K_f为界面两侧的流体体积模量的差值，△f_m为界面两侧固体项的差值，△ρ为界面两侧的密度的差值，△φ为界面两侧的孔隙度的差值，为干岩石纵横波速度比的平方，为饱和岩石纵横波速度比的平方。

6.根据权利要求4所述的方法，其特征是，根根据所述的叠前地震资料、测井资料、地质成果资料以及固液解耦近似方程提取子波具体包括：

从所述的测井资料中提取出地下岩石的密度以及孔隙度；

根据所述的剪切模量曲线以及所述的孔隙度确定地下岩石的固体项；

基于所述的固液解耦近似方程推导出弹性阻抗方程；

根据所述的入射角度、地下岩石的密度、固体项、孔隙度、流体体积模量曲线以及弹性阻抗方程确定弹性阻抗伪测井曲线；

针对不同的入射角度，根据入射角度对应的角度叠加道集、弹性阻抗伪测井曲线、地震层位数据提取不同的子波。

7.根据权利要求6所述的方法，其特征是，所述的弹性阻抗方程为：

$\mathrm{EI}\left(\mathrm{\θ}\right)={\mathrm{EI}}_0{\left(\frac{K_f}{K_{f0}}\right)}^a{\left(\frac{f_m}{f_{m0}}\right)}^b{\left(\frac{\mathrm{\φ}}{{\mathrm{\φ}}_0}\right)}^c{\left(\frac{\mathrm{\ρ}}{{\mathrm{\ρ}}_0}\right)}^d$

其中， $a=(1-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2})\frac{{\sec }^2\mathrm{\θ}}{2},b=\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{2\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{4}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ};$

$c=\frac{{\sec }^2\mathrm{\θ}}{2}-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}+\frac{4}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ},d=1-\frac{{\sec }^2\mathrm{\θ}}{2};$

K_f0、f_m0、ρ₀和φ₀分别为K_f、f_m、ρ和φ的平均值。

8.根据权利要求6所述的方法，其特征是，根据所述的叠前地震资料、测井资料、地质成果资料以及子波确定绝对弹性阻抗数据体具体包括：

根据所述的叠前地震资料、测井资料、地震层位数据以及子波确定出相对弹性阻抗数据体；

根据所述的弹性阻抗伪测井曲线以及地震层位数据构建低频阻抗数据体；

根据所述的低频阻抗数据体对所述的弹性阻抗数据体进行补偿，得到绝对弹性阻抗数据体。

9.一种基于双相介质固液解耦的叠前地震反演系统，其特征是，所述的系统具体包括：

资料采集装置，用于采集目的储层的叠前地震资料、测井资料以及地质成果资料；

固液解耦近似方程构建装置，用于构建固液解耦近似方程；

子波提取装置，用于根据所述的叠前地震资料、测井资料、地质成果资料以及所述的固液解耦近似方程提取子波；

弹性阻抗数据体确定装置，用于根据所述的叠前地震资料、测井资料、地质成果资料以及子波确定绝对弹性阻抗数据体；

流体体积模量确定装置，用于根据所述的绝对弹性阻抗数据体确定流体体积模量数据体；

流体类型识别装置，用于根据所述的流体体积模量数据体识别当前储层孔隙中流体的类型。

10.根据权利要求9所述的系统，其特征是，所述的系统还包括：

保幅处理装置，用于对所述的叠前地震资料进行保幅处理；

去噪处理装置，用于对经过保幅处理的叠前地震资料进行去噪处理；

角度叠加装置，用于根据目的储层对经过去噪处理的叠前地震资料进行角度叠加，生成四个入射角度的角度叠加道集。

11.根据权利要求10所述的系统，其特征是，所述的系统还包括：

环境校正装置，用于对所述的测井资料进行环境校正；

奇异值消除装置，用于对经过环境校正后的测井资料进行奇异值消除；

纵横速度比确定装置，用于根据经过奇异值消除后的测井资料确定所述目的储层的干岩石纵横波速度比、饱和岩石纵横波速度比；

剪切模量曲线确定装置，用于根据所述的目的储层的干岩样纵横波速度比以及所述的测井资料确定与所述测井资料对应的流体体积模量曲线以及剪切模量曲线。

12.根据权利要求11所述的系统，其特征是，所述的系统还包括：

地质解释装置，用于根据所述的叠前地震资料、测井资料、地质成果资料对目的储层进行精细层位追踪，得到地震层位数据。

13.根据权利要求9所述的系统，其特征是，所述的固液解耦近似方程为：

$\begin{array}{c}{R}_{\mathrm{PP}}\left(\mathrm{\θ}\right)=\left[(1-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2})\frac{{\sec }^2\mathrm{\θ}}{4}\right]\frac{{\mathrm{\ΔK}}_f}{K_f}+[\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{4\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ}]\frac{\mathrm{\Δ}\left(f_m\right)}{f_m}[\frac{1}{2}-\frac{{\sec }^2\mathrm{\θ}}{4}]\frac{\mathrm{\Δ\ρ}}{\mathrm{\ρ}}\\ +(\frac{{\sec }^2\mathrm{\θ}}{4}-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{2{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}+\frac{2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ})\frac{\mathrm{\Δ\φ}}{\mathrm{\φ}}\end{array}$

其中，f_m＝φμ为地下岩石的固体项，θ为入射角度，K_f为地下岩石的流体体积模量，f_m为地下岩石的固体项，φ为地下岩石的孔隙度，ρ为地下岩石的密度，μ为地下岩石的剪切模量，△K_f为界面两侧的流体体积模量的差值，△f_m为界面两侧固体项的差值，△ρ为界面两侧的密度的差值，△φ为界面两侧的孔隙度的差值，为干岩石纵横波速度比的平方，为饱和岩石纵横波速度比的平方。

14.根据权利要求12所述的系统，其特征是，所述的子波提取装置具体包括：

孔隙度提取模块，用于从所述的测井资料中提取出地下岩石的密度以及孔隙度；

固体项确定模块，用于根据所述的剪切模量曲线以及所述的孔隙度确定地下岩石的固体项；

弹性阻抗方程推导模块，用于基于所述的固液解耦近似方程推导出弹性阻抗方程；

弹性阻抗伪测井曲线确定模块，用于根据所述的入射角度、地下岩石的密度、固体项、孔隙度、流体体积模量曲线以及弹性阻抗方程确定弹性阻抗伪测井曲线；

子波提取模块，用于针对不同的入射角度，根据入射角度对应的角度叠加道集、弹性阻抗伪测井曲线、地震层位数据提取不同的子波。

15.根据权利要求14所述的系统，其特征是，所述的弹性阻抗方程为：

$\mathrm{EI}\left(\mathrm{\θ}\right)={\mathrm{EI}}_0{\left(\frac{K_f}{K_{f0}}\right)}^a{\left(\frac{f_m}{f_{m0}}\right)}^b{\left(\frac{\mathrm{\φ}}{{\mathrm{\φ}}_0}\right)}^c{\left(\frac{\mathrm{\ρ}}{{\mathrm{\ρ}}_0}\right)}^d$

其中， $a=(1-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2})\frac{{\sec }^2\mathrm{\θ}}{2},b=\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{2\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}-\frac{4}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ};$

$c=\frac{{\sec }^2\mathrm{\θ}}{2}-\frac{{\mathrm{\γ}}_{\mathrm{dry}}^2}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sec }^2\mathrm{\θ}+\frac{4}{{\mathrm{\γ}}_{\mathrm{sat}}^2}{\sin }^2\mathrm{\θ},d=1-\frac{{\sec }^2\mathrm{\θ}}{2};$

K_f0、f_m0、ρ₀和φ₀分别为K_f、f_m、ρ和φ的平均值。

16.根据权利要求14所述的系统，其特征是，所述的弹性阻抗数据体确定装置具体包括：

相对弹性阻抗数据体确定模块，用于根据所述的叠前地震资料、测井资料、地震层位数据以及子波确定出相对弹性阻抗数据体；

低频阻抗数据体构建模块，用于根据所述的弹性阻抗伪测井曲线以及地震层位数据构建低频阻抗数据体；

补偿模块，用于根据所述的低频阻抗数据体对所述的弹性阻抗数据体进行补偿，得到绝对弹性阻抗数据体。