CN103412336A - 一种非均质油藏中岩石系统的纵波速度预测方法 - Google Patents

Abstract

本发明提供一种非均质油藏中岩石系统的纵波速度预测方法，包括：获取渗透率、孔隙率、泥质含量、矿物组分、矿物体积比率以及孔隙流体数据；建立岩石干骨架模型和嵌入体流体模型；根据嵌入体流体模型，计算嵌入体内部的流体速度；根据嵌入体内部的流体速度，计算双重孔隙介质的动能函数和耗散能函数；根据哈密顿原理和拉格朗日方程，利用双重孔隙介质的动能函数和耗散能函数，导出改进后的Biot-Rayleigh方程组；通过平面波分析及改进后的Biot-Rayleigh方程组，预测非均质油藏的岩石纵波速度。本发明充分考虑了岩石系统的非均质性以及嵌入体内部的流体速度场引起的流体动能和耗散能，可以预测“水包油”双重孔隙岩石系统的纵波速度；不仅弥补了原Biot-Rayleigh方程组的劣势，还保留了原Biot-Rayleigh方程组格式简洁的优点。

Description

一种非均质油藏中岩石系统的纵波速度预测方法

技术领域

本发明涉及岩石纵波预测技术领域，具体地，涉及一种非均质油藏中岩石系统的纵波速度预测方法。

背景技术

面对国际能源需求日益增长的严峻形势，油气资源勘探仍是各石油公司的主要工作目标。尽管随着几十年的开采，国内大部分油田已进入高含水后期甚至特高含水期，但剩余油产量仍非常可观，地下剩余油呈“整体高度分散、局部相对富集”的状态，因此在高含水期剩余油的挖掘过程中，开展储层物性参数的预测研究是十分必要的。而岩石物理作为油气勘探的基础，可以提供各种对储层识别及含油气性分析的敏感岩石物理参数，尤其包括连接岩石物理与地震勘探及测井桥梁的纵、横波速度。

目前人们已普遍认识到，不仅岩石中孔隙流体的性质发生变化时，会改变纵、横波速度，甚至波诱导的孔隙流体流动也会产生影响，并造成弹性波的速度频散和衰减，由于弹性波穿过饱含多相流体的多孔介质时，在各种孔隙流体之间会诱发压力梯度，导致孔隙流体流动，直到孔隙压力平衡，因此，利用岩石物理模型分析，可预测部分饱和储层的纵、横波速度，建立其与流体饱和情况之间的联系。

流固双相介质中最早的纵、横波速度预测方法是基于Gassmann和Biot提出的孔隙介质理论，该模型考虑的是波长尺度下的孔隙流体流动，由于忽略了岩石内部的微观非均质性，在超声波频带预测的速度和频散明显偏低，其预测结果在声波频带和地震频带更差，甚至与实验观测冲突。为了解释超声波频段的高衰减与高频散，Mavko和Nur考虑到孔隙规模的流体流动，提出了喷射流机制。考虑到喷射流机制和Biot理论同时存在的可能性，Dvorkin等建立了可同时处理宏观Biot流动机制和微观喷射流机制的BISQ理论。国内杨顶辉等将BISQ理论推广到了一般孔隙各向异性的情况，拓展了BISQ理论的应用范围，但其缺点在于引入了一些并不具有明确物理意义或并不易于直接物理实现的参数或系数（如特征喷射流长度），使得这些理论虽然在进行岩石波速现象的描述与预测中非常有效，但相关理论与核心参数的数学基础与物理内涵却难以得到实现或验证。基于Biot理论和BISQ理论，唐晓明导出了一个可以描述孔隙、裂隙并存的弹性波统一理论，该模型能够较好地解释声波频带的速度频散和衰减，但引入了两个不易直接测量的参数：裂隙密度和孔隙纵横比，并且不能描述地震频带的速度频散和衰减。

对于高含水储层，孔隙中的油水分布形式，较易形成“水包油”的形态，具有双重孔隙的特征，由于两者在弹性模量、密度以及粘滞性方面的差异性与水、气的差异性相比不是很大，在计算局域流引起的动能和耗散函数时，嵌入体内部的流体速度是不可被忽略的。因此，基于流体力学的基本原理，计算嵌入体内部的流体速度场，完善描述整个局域流的振荡过程，进而完善“一类骨架、两类流体”双重孔隙介质理论，不仅可以提供部分饱和不同孔隙流体岩石的纵、横波速度预测方法，而且可以预测“两类骨架、一类流体”岩石中的纵、横波速度，完善后的模型将非常利于储层地震特征正演。

最常用的非均匀孔隙介质中的纵、横波速度预测方法包括：White等于1975年提出的空间周期排列的球状Patchy模型，但其缺点在于对嵌入体内气体的物理参数做了简化和忽略处理。Pride等提出双重孔隙介质理论，可以描述含气孔与含水孔之间的中观流，但其不足在于方程组复杂，难于应用到实际；方程中引入了过多的参数，不利于理论实验验证与实际证明；动力学方程组的推导过程中采用了统计热力学关系，使得方程的一些核心参数在计算过程中仅能通过频率域的关系式进行上、下限的估算，而无法做到更准确的地震响应的定量预测。基于以上考虑，巴晶、曹宏等于2012年针对含气饱和岩石，推导了格式尽可能简洁、物理参数尽可能少、各参量均具备物理可实现性的Biot-Rayleigh方程组，并证实了此类储层中弹性波在地震勘探频段内可能出现强烈的频散现象。但该模型在计算过程中，没有考虑嵌入体内部气体的速度，因此，局域流振荡的整个过程是不完备的，无法预测嵌入体内部的流体速度不可被忽略情形（如“水包油”）下岩石的纵、横波速度。

综上所述，现有技术的研究中主要存在以下问题：

（1）常规流体替换的方法是Biot-Gassmann方程，这种方程的假设前提是岩石内部的所有孔隙是完全均一的，因此这种假设不能满足非均匀复杂地层的工程应用。

（2）巴晶等于2012年提出的双重孔隙介质模型，忽略了嵌入体内的局域流速度场，不能完备描述局部流振动的物理过程，不可适用于嵌入体流体动能不可被忽略的情况，因此主要应用于非均质气藏岩石物理建模，而无法应用于非均质油藏的地震检测工程。

发明内容

本发明实施例的主要目的在于提供一种非均质油藏中岩石系统的纵波速度预测方法，以克服现有技术忽略嵌入体内部的流体速度场以及没有考虑岩石系统的非均质性问题。

为了实现上述目的，本发明实施例提供一种非均质油藏中岩石系统的纵波速度预测方法，包括：

根据测井数据、地质报告、录井数据和岩石实验观测获取渗透率、孔隙率、泥质含量、矿物组分、矿物体积比率，以及根据流体实验测量获取孔隙流体数据；

利用所述渗透率、孔隙率、泥质含量、矿物组分、矿物体积比率建立岩石干骨架模型，以及利用孔隙流体数据建立嵌入体流体模型；

根据所述嵌入体流体模型，计算嵌入体内部的流体速度；

根据所述嵌入体内部的流体速度，计算双重孔隙介质的动能函数和耗散能函数；

根据哈密顿原理和拉格朗日方程，利用所述双重孔隙介质的动能函数和耗散能函数，导出改进后的Biot-Rayleigh方程组；

通过平面波分析及所述改进后的Biot-Rayleigh方程组，预测非均质油藏的岩石纵波速度。

借助于上述技术方案，本发明充分考虑了嵌入体内部的流体速度场引起的流体动能和耗散能，具有可以预测“水包油”双重孔隙岩石系统的纵横波速度能力；充分考虑了岩石系统的非均质性，可以对非均质油藏进行流体替换以及AVO正演；不仅弥补了原Biot-Rayleigh方程组忽略嵌入体内流体物理特征的劣势，还保留了原先Biot-Rayleigh方程组格式简洁的优点，所有相关参数物理可测，具有较好的物理可实现性。

附图说明

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

图1是本发明实施例一提供的非均质油藏中岩石系统的纵波速度预测方法的流程示意图；

图2是本发明实施例二提供的频率为0.01～10000Hz范围内“水包油”双重孔隙介质改进前后的理论预测的纵波速度与衰减变化关系示意图

图3是本发明实施例二提供的饱和度为0～100%范围内“水包油”双重孔隙介质改进前后理论预测的纵波速度变化关系示意图；

图4是本发明实施例三提供的“水包油”双重孔隙介质实验数据与改进前后BR理论预测结果的对比示意图；

图5是本发明实施例四提供的某井深度927.125～992m段改进BR理论预测的速度曲线与实测曲线的对比示意图；

图6（a）、图6（b）、图6（c）、图6（d）是本发明实施例四提供的某井深度927.125～992m段预测纵横波速度的AVO分析结果与实测纵横波速度的AVO分析结果；

图7是本发明实施例四提供的某井深度933.375～992m段改进BR理论预测的速度曲线与实测曲线的对比示意图；

图8（a）、图8（b）、图8（c）、图8（d）是本发明实施例四提供的某井深度927.125～992m段预测纵横波速度的AVO分析结果与实测纵横波速度的AVO分析结果。

具体实施方式

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

实施例一

本实施例提供一种非均质油藏的岩石纵波速度预测方法，如图1所示，该方法包括:

步骤S100，根据测井数据、地质报告、录井数据和岩石实验观测获取渗透率、孔隙率、泥质含量、矿物组分、矿物体积比率，以及根据流体实验测量获取孔隙流体数据；

该步骤的目的是获取准确可信的岩石物理参数和流体参数，从而为预测非均质油藏的岩石纵波速度提供可靠的数据基础。

步骤S101，利用所述渗透率、孔隙率、泥质含量、矿物组分、矿物体积比率建立岩石干骨架模型，以及利用孔隙流体数据建立嵌入体流体模型；

该步骤中，建立岩石干骨架模型主要是确定岩石干骨架的等效弹性模量和岩石干骨架密度，建立嵌入体流体模型主要是确定嵌入体的等效弹性模量和嵌入体流体密度；

步骤S102，根据所述嵌入体流体模型，计算嵌入体内部的流体速度；

步骤S103，根据所述嵌入体内部的流体速度，计算双重孔隙介质的动能函数和耗散能函数；

步骤S104，根据哈密顿原理和拉格朗日方程，利用所述双重孔隙介质的动能函数和耗散能函数，导出改进后的Biot-Rayleigh方程组；

步骤S105，通过平面波分析及所述改进后的Biot-Rayleigh方程组，预测非均质油藏的岩石纵波速度。

具体的，步骤S101中建立岩石干骨架模型具体包括：计算岩石干骨架的密度和等效弹性模量。

其中，岩石干骨架密度可通过计算岩石中各种矿物组分密度的算术平均值获得；确定岩石干骨架的等效弹性模量的具体步骤如下：

步骤A1，采用Voigt-Reuss-Hill平均模型，计算岩石骨架基质的等效弹性模量，Voigt-Reuss-Hill平均模型的表达式如下：

$M_{\mathrm{VRH}}=\frac{1}{2}(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{f}_i{M}_i+\frac{1}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{f_i}{M_i}})$         （公式1）

公式1中，M_VRH为岩石骨架基质的等效弹性模量，包括等效体积模量K_s和等效剪切模量μ_s；f_i、M_i分别为第i种矿物组分的体积率和弹性模量；N为岩石中矿物组分的种类数量；其中，M_i为岩石骨架基质的体积模量或剪切模量。

步骤A2，采用Pride半经验性的公式计算岩石干骨架的等效体积模量和等效剪切模量：

$\left\{\begin{array}{c}{K}_b=\frac{(1-\mathrm{\φ})K_s}{1+\mathrm{c\φ}}\\ {\mathrm{\μ}}_b=\frac{(1-\mathrm{\φ}){\mathrm{\μ}}_s}{1+c^{\′}\mathrm{\φ}}\end{array}\right.$         （公式2）

公式2中，K_b、μ_b分别为岩石干骨架的等效体积模量和等效剪切模量；K_s、μ_s分别为岩石骨架基质的等效体积模量和等效剪切模量；φ为岩石的孔隙度；c、c′分别为经验性参数，与岩石的固结程度有关。

具体的，步骤S101中建立嵌入体流体模型具体包括：利用孔隙流体数据，计算嵌入体流体的等效弹性模量和嵌入体流体的密度。

实际执行步骤S101时，若嵌入体流体的等效弹性模量与流体密度不能通过孔隙流体数据计算得到，则可引用本领域中已知的流体弹性模量和流体密度。

为了提高预测“水包油”双重孔隙岩石中纵波速度的能力，本实施例充分考虑了嵌入体内部的流体速度场引起的流体动能和耗散能。

步骤S102根据所述嵌入体流体模型，计算嵌入体内部的流体速度，具体包括如下过程：

将嵌入体等效为球体（以下简称球形嵌入体），假设球形嵌入体内部的流体是可压缩的，由于嵌入体震动会导致流体流动且球形嵌入体的形变微小，因此球形嵌入体满足流体质量守恒原理，则球坐标系下球形嵌入体的流体质量守恒方程可表示为：

$\frac{\mathrm{d\ρ}}{\mathrm{dt}}+\mathrm{\ρ}\left(\frac{1}{r^2}\frac{\∂}{\∂r}\left(r^2{\stackrel{\·}{r}}_{\mathrm{in}}\right)\right)=0$         （公式3）

公式3中，ρ为球形嵌入体内的流体密度；t表示时间；r为球形嵌入体受地震波激励后的半径；

为球形嵌入体内流体的径向速度；

经整理后，公式3变形为：

$\frac{r^2}{\mathrm{\ρ}}\frac{\mathrm{d\ρ}}{\mathrm{dt}}=-\frac{\∂}{\∂r}\left(r^2{U}_R\right)$         （公式4）

公式4中，

变化很小，可视为常数，同时对公式4等号两边的r积分，可得：

$U_R=-\frac{1}{3\mathrm{\ρ}}\left(\frac{\mathrm{d\ρ}}{\mathrm{dt}}\right)r$         （公式5）

设球形嵌入体的初始半径为R₀，在球形嵌入体变形后的时刻t₀（即当t＝t₀时），球形嵌入体受地震波激励后的动态半径为R，

是半径为R流体球表面的流体速度，C是球形嵌入体的初始质量，o(ε)为质量变化的高阶无穷小量，由于球形嵌入体满足流体质量守恒，则有ρR³＝C+o(ε)，故可以得到：

$\frac{\mathrm{d\ρ}}{\mathrm{dt}}=-\frac{3C{R}^2}{R^6}\stackrel{\·}{R}$         （公式6）

由公式6可得球形嵌入体内部（即r＜R）的流体速度为：

${\stackrel{\·}{r}}_{\mathrm{in}}=\frac{\stackrel{\·}{R}r}{R}$         （公式7）

步骤S103中，根据所述嵌入体内部的流体速度，计算双重孔隙介质的动能函数，具体过程如下：

双重孔隙介质的动能函数可写为：

$T=\frac{1}{2}{\mathrm{\ρ}}_{00}({\stackrel{\·}{u}}_1^2+{\stackrel{\·}{u}}_2^2+{\stackrel{\·}{u}}_3^2)+{\mathrm{\ρ}}_{01}({\stackrel{\·}{u}}_1{\stackrel{\·}{U}}_1^{\left(1\right)}+{\stackrel{\·}{u}}_2{\stackrel{\·}{U}}_2^{\left(1\right)}+{\stackrel{\·}{u}}_3{\stackrel{\·}{U}}_3^{\left(1\right)})$

$+{\mathrm{\ρ}}_{02}({\stackrel{\·}{u}}_1{\stackrel{\·}{U}}_1^{\left(1\right)}+{\stackrel{\·}{u}}_2{\stackrel{\·}{U}}_2^{\left(2\right)}+{\stackrel{\·}{u}}_3{\stackrel{\·}{U}}_3^{\left(2\right)})+\frac{1}{2}{\mathrm{\ρ}}_{11}({\stackrel{\·}{U}}_1^{\left(1\right),2}+{\stackrel{\·}{U}}_2^{\left(1\right),2}+{\stackrel{\·}{U}}_3^{\left(1\right),2})$        （公式8）

$+\frac{1}{2}{\mathrm{\ρ}}_{22}({\stackrel{\·}{U}}_1^{\left(2\right),2}+{\stackrel{\·}{U}}_2^{\left(2\right),2}+{\stackrel{\·}{U}}_3^{\left(2\right),2})+T_{\mathrm{LFF}}$

公式8中，T为双重孔隙介质的动能函数；T_LFF为局域流振荡的动能函数；u₁、u₂、u₃表示固体在直角坐标系三个坐标方向上的位移分量；U₁、U₂、U₃表示流体在直角坐标系三个坐标方向上的位移分量；上标(1)、(2)分别表示两类孔隙；ρ₀₀、ρ₀₁、ρ₀₂、ρ₁₁、ρ₂₂为密度参数；ρ_s为颗粒体积模量；

分别为背景相的密度和嵌入体内部的流体密度；φ为岩石的孔隙度；ν₁、ν₂分别为两种流体所占据的骨架占总骨架的体积比率；φ₁₀与φ₂₀分别为两个区域内部的局部孔隙度。

ρ₀₀、ρ₀₁、ρ₀₂、ρ₁₁、ρ₂₂是五个密度参数，其形式为：

$\left\{\begin{array}{c}{\mathrm{\ρ}}_{00}+{\mathrm{\ρ}}_{01}+{\mathrm{\ρ}}_{02}=(1-\mathrm{\φ}){\mathrm{\ρ}}_s\\ 2{\mathrm{\ρ}}_{01}=({\mathrm{\φ}}_1-v_1){\mathrm{\ρ}}_f^1\\ 2{\mathrm{\ρ}}_{02}=({\mathrm{\φ}}_2-v_2){\mathrm{\ρ}}_f^2\\ 2{\mathrm{\ρ}}_{11}=({\mathrm{\φ}}_1+v_1){\mathrm{\ρ}}_f^1\\ 2{\mathrm{\ρ}}_{22}=({\mathrm{\φ}}_2+v_2){\mathrm{\ρ}}_f^2\\ \mathrm{\φ}={\mathrm{\φ}}_1+{\mathrm{\φ}}_2\\ {\mathrm{\φ}}_1={\mathrm{\ν}}_1{\mathrm{\φ}}_{10}\\ {\mathrm{\φ}}_2={\mathrm{\ν}}_2{\mathrm{\φ}}_{20}\end{array}\right.$        （公式9）

公式9中，φ₁₀与φ₂₀分别表示两个区域内部的局部孔隙度；φ₁为嵌入体外孔隙占整个岩石的绝对孔隙度；φ₂为嵌入体内孔隙占整个岩石的绝对孔隙度；φ为岩石的孔隙度。

根据嵌入体的球体模型，在球坐标系下对嵌入体内部的流体速度进行积分，计算可得到局域流振荡的动能函数，具体形式为：

$T_{\mathrm{LFF}}=\frac{{\mathrm{\φ}}_2}{6}(\frac{{\mathrm{\ρ}}_f}{5}+\frac{{\mathrm{\φ}}_{20}}{{\mathrm{\φ}}_{10}}{\mathrm{\ρ}}_f){\mathrm{\φ}}_1^2{R}_0^2{\left(\stackrel{\·}{\mathrm{\ζ}}\right)}^2$         （公式10）

公式10中，ζ为局域流流动造成的体应变增量。

步骤S104中，根据所述嵌入体内部的流体速度，计算双重孔隙介质的耗散能函数，具体过程如下：

基于孔隙流体与固体骨架的摩擦耗散机制，双重孔隙介质的耗散函数具体形式：

$\left\{\begin{array}{c}D=\frac{1}{2}{b}_1(\stackrel{\·}{u}-{\stackrel{\·}{U}}^{\left(1\right)})(\stackrel{\·}{u}-{\stackrel{\·}{U}}^{\left(1\right)})+\frac{1}{2}{b}_2(\stackrel{\·}{u}-{\stackrel{\·}{U}}^{\left(1\right)})(\stackrel{\·}{u}-{\stackrel{\·}{U}}^{\left(1\right)})+D_{\mathrm{LFF}}\\ \stackrel{\·}{u}=(u_1,u_2,u_3)\\ {\stackrel{\·}{U}}^{\left(1\right)}=(U_1^{\left(1\right)},U_2^{\left(1\right)},U_3^{\left(1\right)})\\ {\stackrel{\·}{U}}^{\left(2\right)}=(U_1^{\left(2\right)},U_2^{\left(2\right)},U_3^{\left(2\right)})\end{array}\right.$         （公式11）

公式11中，D为双重孔隙介质的耗散能函数；D_LFF为局域流振荡引起的耗散函数；b₁、b₂分别为嵌入体外、内的耗散函数；κ₁为背景相中的渗透率；η₁为背景相中的流体黏度；η₂为嵌入体中的流体黏度；

b₁、b₂的形式为：

$\left\{\begin{array}{c}{b}_1=v_1{\mathrm{\φ}}_{10}^2\left(\frac{{\mathrm{\η}}_1}{{\mathrm{\κ}}_1}\right)={\mathrm{\φ}}_1{\mathrm{\φ}}_{10}\left(\frac{{\mathrm{\η}}_1}{{\mathrm{\κ}}_1}\right)\\ b_2=v_2{\mathrm{\φ}}_{20}^2\left(\frac{{\mathrm{\η}}_2}{{\mathrm{\κ}}_1}\right)={\mathrm{\φ}}_2{\mathrm{\φ}}_{20}\left(\frac{{\mathrm{\η}}_2}{{\mathrm{\κ}}_1}\right)\end{array}\right.$         （公式12）

根据嵌入体的球体模型，在球坐标系下对嵌入体内部的流体速度进行积分，可计算得到局域流振荡的耗散能函数D_LFF，D_LFF的具体形式为：

$D_{\mathrm{LFF}}=(2\mathrm{\π}{\mathrm{\φ}}_{2o}^2\left(\frac{{\mathrm{\η}}_1}{{\mathrm{\κ}}_1}\right)+\frac{2\mathrm{\π}}{5}{\mathrm{\φ}}_{2o}^2\left(\frac{{\mathrm{\η}}_2}{{\mathrm{\κ}}_1}\right))\frac{{\mathrm{\φ}}_2{\mathrm{\φ}}_1^2{R}_0^2{\left(\stackrel{\·}{\mathrm{\ζ}}\right)}^2}{12\mathrm{\π}{\mathrm{\φ}}_{20}}$         （公式13）

步骤S105中，根据哈密顿原理和拉格朗日方程，利用所述岩石干骨架模型、嵌入体流体模型、双重孔隙介质的动能函数和耗散能函数，导出改进后的Biot-Rayleigh方程组，该Biot-Rayleigh方程组具体为如下：

    （公式14）

公式14中，e、ξ分别为固相、液相体应变；ξ₁、ξ₂分别为嵌入体外部和内部液相体应变；

分别为局域流流动造成的体应变增量的速度和加速度；ρ_out、ρ_in分别为背景相、嵌入体内流体的密度；A、N、Q₁、R₁、Q₂、R₂分别为双孔隙介质中的Biot弹性参数；

分别为背景相和嵌入体中的流体体积模量。

步骤S106中，通过平面波分析及所述改进后的Biot-Rayleigh方程组，预测非均质油藏的岩石纵波速度，具体包括：

通过平面波分析，将位移场平面波解代入到所述改进的Biot-Rayleigh方程组中，并将所述改进后的Biot-Rayleigh方程组转换到频率波数域，得到Christoffel方程；

位移场平面波解包括：

$\left\{\begin{array}{c}{u}_j=A_j{e}^{i[k(x{l}_1+y{l}_2+z{l}_3)-\mathrm{\ωt}]}\\ U_j={\stackrel{\‾}{A}}_j{e}^{i[k(x{l}_1+y{l}_2+z{l}_3)-\mathrm{\ωt}]}\\ j=x,y,z\end{array}\right.$         （公式15）

公式15中ω是角频率，k是波数，

Christoffel方程具体为：

$\left|\begin{array}{ccc}{a}_{11}{k}^2+b_{11}& a_{12}{k}^2+b_{12}& a_{13}{k}^2+b_{13}\\ a_{21}{k}^2+b_{11}& a_{22}{k}^2+b_{12}& a_{23}{k}^2+b_{13}\\ a_{31}{k}^2+b_{11}& a_{32}{k}^2+b_{12}& a_{33}{k}^2+b_{13}\end{array}\right|=0$         （公式16）

公式16中，各参数如下：

再利用Christoffel方程的行列式等于零，得到岩石纵波的频散关系；

最后由复波数

生成纵波速度的速度预测公式

和逆品质因子 $Q_i^{-1}=\frac{2\mathrm{Im}\left(\widetilde{k_i}\right)}{\mathrm{Re}\left(\widetilde{k_i}\right)}.$

本实施例充分考虑了嵌入体内部的流体速度场引起的流体动能和耗散能，具有可以预测“水包油”双重孔隙岩石系统的纵波速度能力；充分考虑了岩石系统的非均质性，可以对非均质油藏进行流体替换以及AVO正演；不仅弥补了原Biot-Rayleigh方程组忽略嵌入体内流体物理特征的劣势，还保留了原先Biot-Rayleigh方程组格式简洁的优点，所有相关参数物理可测，具有较好的物理可实现性。

实施例二

本实施例对全频带“油包水”双重孔隙介质中纵波的速度与衰减进行定量预测：

本实施例设计一个油水非饱和岩石物理模型，基本的岩石物理参数为：基质的的体积模量35Gpa，骨架的体积模量7Gpa，骨架的剪切模量9Gpa，基质的平均密度2650Kg/m³，油的体积模量0.6Gpa，其密度为900Kg/m³，水的体积模量2.25Gpa，其密度为990Kg/m³，孔隙度为0.15，渗透率为0.1×10^-12m²，水的黏度0.001Pa·s，油的黏度0.006Pa·s，平均气泡尺寸0.25m，图2及图3分别为频率为0.01～10000Hz，饱和度为0～100%时改进前后理论预测的纵波速度与衰减变化示意图。

实施例三

本实施例对阿克苏市苏盖特布拉克地区油水饱和岩石样品实验数据与理论预测数据进行对比：

该样品主要成分为白云石,胶结程度很高，渗透率为0.174md，孔隙度为0.0547，设置的平均嵌入体尺寸为0.035mm，其主要物性参数分别为：

(i)岩石骨架的体积模量76.2Gpa，剪切模量30.67Gpa；

(ii)岩石颗粒的体积模94.9Gpa，剪切模量45Gpa，密度2870Kg/m³；

(iii)实验中所采用的流体为油和水，水的体积模量2.25Gpa，油的体积模量0.6Gpa，水的黏度0.001Pa·s，油的黏度0.006Pa·s，水的密度990Kg/m³，油的密度900Kg/m³；

(iv)实验中所采用的频率为800KHz，图4为实验数据与改进前后BR理论预测结果的对比示意图。

改进后的BR理论预测结果随饱和度的变化关系与实验数据吻合较好。

实施例四

本实施例根据测井数据预测井中地层纵波、横波速度：

使用测井数据中的孔隙度、含水饱和度和泥质含量，预测了某两口井深度段分别为927.125～992m和933.375～992m的纵、横波速度，目的岩层是油、水饱和状态，岩石骨架矿物采用石英和粘土。

岩石参数为：石英的体积模量38Gpa，剪切模量36Gpa，密度为2.65g/cm³，粘土的体积模量20Gpa，剪切模量6Gpa，密度为2.58g/cm³，基质的体积模量和剪切模量由Voigt-Reuss-Hill计算得到，基质的平均密度按Voigt公式计算，骨架体积模量和剪切模量由孔隙度和基质体积模量、剪切模量计算得到（Pride公式），水的体积模量2.4Gpa，其黏度为0.001Pa·s，其密度为1.01g/cm³，油的体积模量1.4Gpa，其黏度为0.02Pa·s，其密度为0.88g/cm³。选取的平均嵌入体尺寸9mm，声波测井所采用的频率为10kHz，图5、图7为改进BR理论预测的速度曲线与实测曲线的对比，图6(a)～(d)与图8(a)～(d)分别是预测纵横波速度的AVO分析结果与实测纵横波速度的AVO分析结果，预测与实测地震响应规律基本一致。

综上所述，本发明实施例提供的非均质油藏中岩石系统的纵、横波速度预测方法具有以下有益效果：

（1）本发明考虑了嵌入体内部的流体动能、耗散能，能够真实反映并准确预测地震波在储层中的频散和衰减，由于油与气相比，嵌入体流体的密度与体积模量是不可忽略的，因此本发明适用于非均质油藏的地震检测工程；

（2）不仅弥补了Biot-Rayleigh方程对嵌入体内的流体性质欠考虑的劣势，还保留了原先Biot-Rayleigh方程组的优点，因此本发明更具有一般性特征，例如在实施例三中，在对阿克苏市苏盖特布拉克地区油水非饱和岩石样品纵波速度的预测中，通过与Biot-Rayleigh方法对比，本发明“水包油”双孔算例的速度预测更加接近实验数据，也证明“水包油”模型更符合实际岩石内部的油、水分布的客观情况；

（3）采用改进后的Biot-Rayleigh方程组进行流体替换，用于估算“水包油”双重孔隙介质中的纵波速度尤其具有优势，通过估算的纵波速度进行AVO正演，进而指导地震勘探，查明地下目标区域内的流体饱和情况。

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

Claims (7)

1.一种非均质油藏中岩石系统的纵波速度预测方法，其特征在于，包括：

根据测井数据、地质报告、录井数据和岩石实验观测获取渗透率、孔隙率、泥质含量、矿物组分、矿物体积比率，以及根据流体实验测量获取孔隙流体数据；

利用所述渗透率、孔隙率、泥质含量、矿物组分、矿物体积比率建立岩石干骨架模型，以及利用孔隙流体数据建立嵌入体流体模型；

根据所述嵌入体流体模型，计算嵌入体内部的流体速度；

根据所述嵌入体内部的流体速度，计算双重孔隙介质的动能函数和耗散能函数；

根据哈密顿原理和拉格朗日方程，利用所述双重孔隙介质的动能函数和耗散能函数，导出改进后的Biot-Rayleigh方程组；

通过平面波分析及所述改进后的Biot-Rayleigh方程组，预测非均质油藏的岩石纵波速度。

2.根据权利要求1所述的方法，其特征在于，建立岩石干骨架模型具体包括：通过计算岩石中各种矿物组分密度的算术平均值获得岩石干骨架的密度；以及，通过如下过程计算岩石干骨架的等效弹性模量：

采用如下公式计算岩石骨架基质的等效弹性模量：

$M_{\mathrm{VRH}}=\frac{1}{2}(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{f}_i{M}_i+\frac{1}{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\frac{f_i}{M_i}})$

其中，M_VRH为岩石骨架基质的等效弹性模量，包括等效体积模量K_s和等效剪切模量μ_s；f_i、M_i分别为第i种矿物组分的体积率和弹性模量；N为岩石中矿物组分的种类数量；M_i为岩石骨架基质的体积模量或剪切模量；

采用如下公式计算岩石干骨架的等效体积模量和等效剪切模量：

$\left\{\begin{array}{c}{K}_b=\frac{(1-\mathrm{\φ})K_s}{1+\mathrm{c\φ}}\\ {\mathrm{\μ}}_b=\frac{(1-\mathrm{\φ}){\mathrm{\μ}}_s}{1+c^{\′}\mathrm{\φ}}\end{array}\right.$

其中，K_b、μ_b分别为岩石干骨架的等效体积模量和等效剪切模量；K_s、μ_s分别为岩石骨架基质的等效体积模量和等效剪切模量；φ为岩石的孔隙度；c、c′分别为经验性参数。

3.根据权利要求2所述的方法，其特征在于，采用如下公式计算嵌入体内部的流体速度：

$\left\{\begin{array}{c}\frac{\mathrm{d\ρ}}{\mathrm{dt}}+\mathrm{\ρ}\left(\frac{1}{r^2}\frac{\∂}{\∂r}\left(r^2{\stackrel{\·}{r}}_{\mathrm{in}}\right)\right)=0\\ \frac{r^2}{\mathrm{\ρ}}\frac{\mathrm{d\ρ}}{\mathrm{dt}}=-\frac{\∂}{\∂r}\left(r^2{U}^R\right)\\ U_R=-\frac{1}{3\mathrm{\ρ}}\left(\frac{\mathrm{d\ρ}}{\mathrm{dt}}\right)r\\ \mathrm{\ρ}{R}^3=C+o\left(\mathrm{\ϵ}\right)\\ \frac{\mathrm{d\ρ}}{\mathrm{dt}}=-\frac{3C{R}^2}{R^6}\stackrel{\·}{R}\\ {\stackrel{\·}{r}}_{\mathrm{in}}=\frac{\stackrel{\·}{R}r}{R}\end{array}\right.$

其中，ρ为嵌入体内的流体密度；t为时间；r为嵌入体受地震波激励后的半径，且r＜R；为嵌入体内流体的径向速度；R₀为嵌入体的初始半径；R为t＝t₀时刻嵌入体受地震波激励后的的动态半径；

是动态半径为R的嵌入体表面的流体速度，C为嵌入体的初始质量，o(ε)为质量变化的高阶无穷小量。

4.根据权利要求3所述的方法，其特征在于，采用如下公式计算双重孔隙介质的动能函数：

其中，T为双重孔隙介质的动能函数；T_LFF为局域流振荡的动能函数；u₁、u₂、u₃表示固体在直角坐标系三个方向上的位移分量；U₁、U₂、U₃表示流体在直角坐标系三个方向上的位移分量；上标(1)、(2)分别表示两类孔隙；ρ₀₀、ρ₀₁、ρ₀₂、ρ₁₁、ρ₂₂为密度参数；ρ_s为颗粒体积模量；

分别为背景相的密度和嵌入体内部的流体密度；ν₁、ν₂分别为两种流体所占据的骨架占总骨架的体积比率；φ₁₀与φ₂₀分别为两个区域内部的局部孔隙度；φ₁为嵌入体外孔隙占整个岩石的绝对孔隙度；φ₂为嵌入体内孔隙占整个岩石的绝对孔隙度；ζ为局域流流动造成的体应变增量。

5.根据权利要求4所述的方法，其特征在于，采用如下公式计算双重孔隙介质的耗散能函数：

$\left\{\begin{array}{c}D=\frac{1}{2}{b}_1(\stackrel{\·}{u}-{\stackrel{\·}{U}}^{\left(1\right)})(\stackrel{\·}{u}-{\stackrel{\·}{U}}^{\left(1\right)})+\frac{1}{2}{b}_2(\stackrel{\·}{u}-{\stackrel{\·}{U}}^{\left(1\right)})(\stackrel{\·}{u}-{\stackrel{\·}{U}}^{\left(1\right)})+D_{\mathrm{LFF}}\\ \stackrel{\·}{u}=(u_1,u_2,u_3)\\ {\stackrel{\·}{U}}^{\left(1\right)}=(U_1^{\left(1\right)},U_2^{\left(1\right)},U_3^{\left(1\right)})\\ {\stackrel{\·}{U}}^{\left(2\right)}=(U_1^{\left(2\right)},U_2^{\left(2\right)},U_3^{\left(2\right)})\\ b_1=v_1{\mathrm{\φ}}_{10}^2\left(\frac{{\mathrm{\η}}_1}{{\mathrm{\κ}}_1}\right)={\mathrm{\φ}}_1{\mathrm{\φ}}_{10}\left(\frac{{\mathrm{\η}}_1}{{\mathrm{\κ}}_1}\right)\\ b_2=v_2{\mathrm{\φ}}_{20}^2\left(\frac{{\mathrm{\η}}_2}{{\mathrm{\κ}}_1}\right)={\mathrm{\φ}}_2{\mathrm{\φ}}_{20}\left(\frac{{\mathrm{\η}}_2}{{\mathrm{\κ}}_1}\right)\\ D_{\mathrm{LFF}}=(2\mathrm{\π}{\mathrm{\φ}}_{2o}^2\left(\frac{{\mathrm{\η}}_1}{{\mathrm{\κ}}_1}\right)+\frac{2\mathrm{\π}}{5}{\mathrm{\φ}}_{2o}^2\left(\frac{{\mathrm{\η}}_2}{{\mathrm{\κ}}_1}\right))\frac{{\mathrm{\φ}}_2{\mathrm{\φ}}_1^2{R}_0^2{\left(\stackrel{\·}{\mathrm{\ζ}}\right)}^2}{12\mathrm{\π}{\mathrm{\φ}}_{20}}\end{array}\right.$

其中，D为双重孔隙介质的耗散能函数；D_LFF为局域流振荡引起的耗散函数；b₁、b₂分别为嵌入体外、内的耗散函数；κ₁为背景相中的渗透率；η₁为背景相中的流体黏度；η₂为嵌入体中的流体黏度。

6.根据权利要求5所述的方法，其特征在于，所述改进后的Biot-Rayleigh方程组为：

其中，e、ξ分别为固相、液相体应变；ξ₁、ξ₂分别为嵌入体外部和内部液相体应变；

分别为局域流流动造成的体应变增量的速度和加速度；ρ_out、ρ_in分别为背景相、嵌入体内流体的密度；A、N、Q₁、R₁、Q₂、R₂分别为双孔隙介质中的Biot弹性参数；

分别为背景相和嵌入体中的流体体积模量。

7.根据权利要求6所述的方法，其特征在于，通过平面波分析及所述改进后的Biot-Rayleigh方程组，预测非均质油藏的岩石纵波速度，具体包括：

通过平面波分析，将位移场平面波解代入到所述改进的Biot-Rayleigh方程组中，并将所述改进后的Biot-Rayleigh方程组转换到频率波数域，得到Christoffel方程；

利用Christoffel方程的行列式等于零，得到岩石纵波的频散关系；

由复波数

生成纵波速度的速度预测公式和逆品质因子

其中，所述位移场平面波解包括：

$\left\{\begin{array}{c}{u}_j=A_j{e}^{i[k(x{l}_1+y{l}_2+z{l}_3)-\mathrm{\ωt}]}\\ U_j={\stackrel{\‾}{A}}_j{e}^{i[k(x{l}_1+y{l}_2+z{l}_3)-\mathrm{\ωt}]}\\ j=x,y,z\end{array}\right.$

所述Christoffel方程为：

$\left|\begin{array}{ccc}{a}_{11}{k}^2+b_{11}& a_{12}{k}^2+b_{12}& a_{13}{k}^2+b_{13}\\ a_{21}{k}^2+b_{11}& a_{22}{k}^2+b_{12}& a_{23}{k}^2+b_{13}\\ a_{31}{k}^2+b_{11}& a_{32}{k}^2+b_{12}& a_{33}{k}^2+b_{13}\end{array}\right|=0$

所述Christoffel方程中各参数为：

其中，ω是角频率，k是波数，

Images