CN103984027A - 基于椭球体双重孔隙模型的岩石纵波速度预测方法 - Google Patents

Abstract

本发明提供一种基于椭球体双重孔隙模型的岩石纵波速度预测方法，包括以下步骤：通过获取测井数据，实验观测数据，得到岩石物理参数，并生成岩石的干骨架模型和孔隙流体模型；建立双重孔隙介质椭球斑块饱和模型，并计算孔隙介质的势能/动能，以及内嵌入体中流体的动能和耗散方程，推导出拉格朗日方程组，并求取纵横波速度；根据平面波分析方法，得到波动方程的频散关系，并得到纵波速度频散和衰减计算公式。本发明提高了岩石纵波速度的预测准确性。

Description

基于椭球体双重孔隙模型的岩石纵波速度预测方法

技术领域

本发明涉及一种地震岩石物理领域的速度频散和衰减的预测方法，特别涉及一种包含椭球形状非均匀掺杂体/流体饱和区域的孔隙介质的斑块模型的纵波速度预测方法。

背景技术

近年来工业高速发展，油气资源需求日益增长，寻找新的资源，勘探油区成为我们的首要目标。随着油气资源的大幅度开采，地下油呈现“整体高度分散，局部相对富集”的状态，给开采带来了很大的困难。这就对勘探技术提出了更高的要求，需要提高声波勘探的分辨率，为最佳井位部署提供依据。利用声波的方法来圈闭油气藏和描述储层特性，是目前我们广泛使用的方法之一，其中纵横波速度是我们的重点描述对象。

在石油勘探过程中，常常需要研究声波在复杂地质构造中的传播问题。岩石中孔隙流体的性质发生变化时，会改变纵、横波速度，地震波诱导的孔隙流体局部流动也会产生影响，并造成弹性波的速度频散和衰减。弹性波穿过饱含流体的多孔介质时，在孔隙流体中会诱发压力梯度，导致孔隙流体流动，直到孔隙压力平衡。同时孔隙流体的流动可能发生在不同的尺度下，因此利用岩石物理模型预测储层纵、横波速度，关键在于明确不同尺度下的流体速度场。大部分油气储层都可看为具有双重孔隙特征的介质。在水-油、水-气、气-油饱和模型中，都可以通过其弹性模量、粘性、密度等参数来计算动能，势能和耗散函数，得到岩石纵、横波速度。

然而，现有技术中岩石纵波速度的预测方法的所有模型都将孔隙介质中的斑块饱和单元体等同为一致的同心球体。但现实中的孔隙介质饱和斑块区域往往不会是理想球体，这导致我们建立的动能与耗散方程与实际的有所出入，不可避免地引起地震波速度预测误差，造成预测的不准确。

发明内容

综上所述，确有必要提供一种能够提高预测准确度的岩石纵波速度的预测方法。

一种基于椭球体双重孔隙模型的岩石纵波速度预测方法，包括以下步骤：步骤S10，通过获取测井数据，实验观测数据，得到岩石物理参数，并生成岩石的干骨架模型和孔隙流体模型；步骤S20，以岩石的干骨架模型和孔隙流体模型为基础，建立双重孔隙介质椭球斑块饱和模型，所述双重孔隙介质椭球斑块饱和模型包括三层具有不同孔隙度和不同流体饱和的椭球壳区域模型，每层球壳代表具有不同孔隙特征和不同流体特征的区域，并计算孔隙介质的势能/动能，以及内嵌入体中流体的动能和耗散方程，推导出拉格朗日方程组，并求取纵横波速度；步骤S30，根据平面波分析方法，得到波动方程的频散关系，并得到纵波速度频散和衰减计算公式。

与现有技术相比较，本发明通过引入椭球形斑块饱和模型，并且包含了固体骨架和流体的空间非均匀性，因此其速度预测更具灵活性和贴近实际的特点，可以解决较为复杂的多流体、多孔隙类型岩石纵波预测问题。

附图说明

图1为本发明第一实施例提供的椭球体斑块饱和模型的结构示意图。

图2为椭球体三层斑块孔隙度和流体分布示意图；三层椭球壳具有不同的孔隙度φ₁,φ₂,φ₃和不同流体f₁,f₂,f₃，不同区域之间交界面为S₁,S₂。

图3为椭球体流体饱和斑块主轴半径和角度示意图。a,b,c为椭球的主半径，为椭球面上任意一点对应的空间角度

图4为三层椭球体斑块饱和模型的孔隙介质纵波速度预测实施方式流程图。

图5为球体斑块饱和模型与其他三种方法对Fort Union砂岩纵波速度预测结果对比。纵波频率为5千赫兹，其中BR代表Biot-Rayleigh方法；Johnson代表Johnson2001年提出的方法，White代表修正后的Whtie方法；TLP代表球体斑块饱和方法；Sw表示含水饱和度，Vp代表纵波速度，单位是米/秒(m/s)。

图6为椭球形斑块部分饱和模型与其他三种Fort Union砂岩5千赫兹纵波速度预测结果对比；其中BR代表Biot-Rayleigh方法；Johnson代表Johnson的方法，White代表修正后的Whtie方法；代表椭球体斑块饱和方法，椭球主轴比例r_ba采用统一等效值0.005;代表椭球体斑块饱和方法，椭球主轴比例r_ba随含水饱和度变化；Sw表示含水饱和度，Vp代表纵波速度，单位是米/秒(m/s)。

图7为椭球部分饱和模型与其他两种纵波速度预测的北海砂岩波速对比（数据50赫兹，50千赫兹和500千赫兹）；其中，BGW和BGH代表低频和高频极限理论预测速度，Gassmann代表低频纵波的Gassmann理论预测速度，White代表修正后的White方法，Johnson代表Johnson提出的方法；TLP代表椭球体斑块饱和方法，椭球主轴比例r_ba采用统一等效值0.001；Vp代表纵波速度，单位是米/秒(m/s)。

具体实施方式

以下将结合附图详细说明本发明提供的岩石纵波速度的预测方法。

请参阅图1，本发明提供的岩石纵波速度的预测方法包括如下步骤：

步骤S10，通过获取测井数据，实验观测数据，得到包括流体饱和孔隙介质的渗透率、孔隙度的岩石数据及孔隙中的流体数据，并生成岩石的干骨架模型和孔隙流体模型；

步骤S20，建立双重孔隙介质椭球形斑块饱和模型，并获取孔隙介质的势能/动能，以及内嵌入体中流体的动能和耗散方程，推导出拉格朗日方程组，求取相应情况的岩石纵、横波速度；

步骤S30，根据平面波分析方法，得到波动方程的频散关系，并得到纵波速度频散和衰减计算公式。

在步骤S10中，首先，通过测井数据、实验观测数据获得可靠的岩石物理参数，为测量双重孔隙岩石中的纵、横波速度提供参数。

具体的，根据目标区域的地质报告、测井资料、岩芯切片等手段，得到矿物成分、矿物体积比率、渗透率、孔隙率、泥质含量参数；根据目标区域地层的温度、压力、矿化度等信息，对流体特征分析，确定流体的密度、粘性、弹性模量等参数。

所述岩石物理参数包括岩石的矿物成分、矿物体积比率、渗透率、孔隙率、泥质含量，以及孔隙中的流体数据包括流体的密度、粘性、弹性模量。所述岩石物理参数可由目标区域的地质报告、测井资料、岩芯切片、地层的温度、压力、矿化度等信息获得。然后再采用Voigt-Reuss-Hill平均模型计算岩石本身的等效弹性模量：

$M_{\mathrm{VRH}}=\frac{1}{2}(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{v}_i{M}_i+1/\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{v_i}{M_i});$

其中：

M_VRH：岩石本身的等效弹性模量；

v_i：第i种矿物的体积率；

M_i：第i种矿物的弹性模量；

n：岩石中矿物的总种类数；

采用Pride半经验性的公式来计算干骨架的等效体积模量和剪切模量:

$K_b=\frac{1-\mathrm{\φ}}{1+\mathrm{c\φ}}{K}_s,{\mathrm{\μ}}_b=\frac{1-\mathrm{\φ}}{1+c\text{'}\mathrm{\φ}}{\mathrm{\μ}}_s;$

K_b,μ_b:岩石干骨架的体积模量和剪切模量；

K_s,μ_s:岩石基质的体积模量和剪切模量；

φ:岩石的孔隙度；

c，c′:为经验性参数，与岩石的固结程度有关。

在步骤S20中，请一并参阅图2，以岩石的干骨架模型和孔隙流体模型为基础建立所述椭球板块饱和模型，所述椭球斑块饱和模型包含双重孔隙和双重流体，以模拟固体骨架非均匀性和流体非均匀性同时存在的情况。为了建立该模型的波动方程，首先需要计算椭球内嵌入体中的流体速度场。因为椭球三条主轴的不同，导致内部流体流动发生很大变化。

所述椭球体斑块模型的波动方程建立包括以下步骤：

步骤S21，获取椭球形嵌入体内部的流体速度特征，对双重孔隙介质中流体动能函数、耗散函数进行求取。

步骤S22，通过哈密顿原理和拉格朗日方程，建立双重孔隙介质的波动方程。

请一并参阅图3，根据三层具有不同孔隙度和不同流体饱和的椭球壳区域模型，每层球壳代表具有不同孔隙特征和不同流体特征的区域，设三层椭球壳具有不同的孔隙度φ₁,φ₂,φ₃和不同流体f₁,f₂,f₃，不同区域之间交界面为S₁,S₂。在不同区域之间存在流体的局部流动。由于模型中存在不同流体分界面和不同孔隙度固体骨架分界面，两种非均匀性相互关系较为复杂，假设存在内外两种流体F1和F2，内外两种固体Φ₁和Φ₂，则可以分别考虑两种情况：

情况1：流体分界面在固体分界面内部。这种情况表示流体F1完全包含在固体掺杂物Φ₁内部。此时，三层区域的孔隙特征参数满足φ₁=φ₂=Φ₁,φ₃=Φ₂，流体特征参数满足f₂=f₃=F₂,f₁=F₁。

情况2：流体分界面在固体分界面外部。这种情况表示流体F1超出了固体掺杂物Φ₁范围。此时，三层区域的孔隙特征参数满足φ₂=φ₃=Φ₂,φ₁=Φ₁，流体特征参数满足f₁=f₂=F₁,f₃=F₂。

从三层椭球体斑块部分饱和模型出发，根据流体饱和度和固体掺杂体含量确定斑块模型的岩石参数。这包括不同孔隙度的掺杂体半径、斑块椭球体半径比例关系。将确定后的岩石干骨架/固体有机质参数和流体参数带入双重孔隙双重流体波动方程。

所述双重孔隙介质中椭球体斑块的动能和耗散能量计算方法包括：

首先，计算椭球体斑块的动能，势能和耗散能量。

在讨论椭球体形状瞬时变化时，为了简化方程，将椭球体看为等比例变化，即有：

$\frac{a}{a_{m0}}=\frac{b}{b_{m0}}=\frac{c}{c_{m0}}=\frac{r}{r_{m0}};$

其中，a,b,c为椭球的主轴半径瞬时值，a₀,b₀,c₀为椭球的主轴半径初始值，为椭球面上某一点对应的空间角度r为该点对应的半径(见图3)，下标m=1,2,3表示三层球壳的不同区域。

所述三层椭球斑块饱和系统的动能函数可写为：

$T=\frac{1}{2}\underset{i=1}{\overset{3}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{00}{\stackrel{\·}{u}}_i^2+\underset{i=1}{\overset{3}{\mathrm{\Σ}}}\underset{m=1}{\overset{3}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{0m}{\stackrel{\·}{u}}_i{\stackrel{\·}{U}}_i^{\left(m\right)}+\frac{1}{2}\underset{i=1}{\overset{3}{\mathrm{\Σ}}}\underset{m=1}{\overset{3}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{mm}}{\stackrel{\·}{U}}_i^{{\left(m\right)}^2}+T_L;$

其中表示第m类孔隙区域内的流体空间位移分量，u₁,u₂,u₃表示固体位移分量，ρ_0m=ρ_m-ρ_mm，ρ_m=φ_mρ_fm, 这里ρ_fm表示三层斑块模型中第m区域的流体密度，a是与孔隙几何特征有关的系数。φ为岩石总的孔隙度，v_m分别表示分别表示三个区域所占体积比，φ_m分别表示三个区域内部的局部孔隙度。

T_L是双重孔隙之间的差异引起的局部流动动能：

$\begin{array}{c}{T}_L=\frac{1}{24\mathrm{\π}}{\mathrm{\ρ}}_{f_1}{a}_{10}^2{\mathrm{\φ}}_1{\mathrm{\φ}}_2^2+{\stackrel{\·}{\mathrm{\ζ}}}_1^{2}E+\frac{{\mathrm{\ρ}}_{f_2}}{18}\frac{{\mathrm{\φ}}_{10}}{{\mathrm{\φ}}_{20}}{a}_{10}^2{\mathrm{\φ}}_1{\mathrm{\φ}}_2^2(1+\frac{b_{10}^2}{a_{10}^2}+\frac{c_{10}^2}{a_{10}^2}){\stackrel{\·}{\mathrm{\ζ}}}_1^2+\\ \{-{\mathrm{\ρ}}_{f_2}+{\mathrm{\ρ}}_{f_3}\frac{{\mathrm{\φ}}_{20}}{{\mathrm{\φ}}_{30}}[1-{(\frac{{\mathrm{\φ}}_1}{{\mathrm{\φ}}_{10}}+\frac{{\mathrm{\φ}}_2}{{\mathrm{\φ}}_{20}})}^{1/3}\left]\right\}\frac{1}{18}(1+\frac{b_{10}^2}{a_{10}^2}+\frac{c_{10}^2}{a_{10}^2})a_{20}^2{\mathrm{\φ}}_{20}{\mathrm{\φ}}_3^2(\frac{{\mathrm{\φ}}_1}{{\mathrm{\φ}}_{10}}+\frac{{\mathrm{\φ}}_2}{{\mathrm{\φ}}_{20}}){\stackrel{\·}{\mathrm{\ζ}}}_2^2\end{array};;$

其中，ζ为局域流流动造成的体应变增量，表示随时间的变化率，a₁₀,b₁₀,c₁₀是椭球体初始主轴半径，E是与椭球体积分有关的项：

椭球形双重孔隙模型的势能可以表示为：

$\frac{1}{2}[(A+2N)I_1^2-4N{I}_2+2Q_1{I}_1({\mathrm{\ξ}}_1-{\mathrm{\φ}}_2{\mathrm{\ζ}}_1)+R_1{({\mathrm{\ξ}}_1-{\mathrm{\φ}}_2{\mathrm{\ζ}}_1)}^2+2Q_2{I}_1({\mathrm{\ξ}}_2+{\mathrm{\φ}}_1{\mathrm{\ζ}}_1-{\mathrm{\φ}}_3{\mathrm{\ζ}}_2)+;$

$R_2{({\mathrm{\ξ}}_2+{\mathrm{\φ}}_1{\mathrm{\ζ}}_1-{\mathrm{\φ}}_3{\mathrm{\ζ}}_2)}^2+2Q_3{I}_1({\mathrm{\ξ}}_3+{\mathrm{\φ}}_2{\mathrm{\ζ}}_2)+R_3{({\mathrm{\ξ}}_3+{\mathrm{\φ}}_2{\mathrm{\ζ}}_2)}^2\]$

其中为第一和第二应力不变量，ξ_m（m=1,2,3）表示流体的第一应变不变量。ζ_m表示流体在不同孔隙区域边界间流动引起的流体体积含量变化，N是固体骨架剪切模量，A是固体弹性模量，Q_m，m=1,2,3是固体、流体耦合弹性模量，R_m，m=1,2,3是流体弹性模量。

基于孔隙流体与固体骨架的摩擦耗散机制,双重孔隙介质的耗散函数具体形式可以表示为：

$D=\frac{1}{2}\underset{m=1}{\overset{3}{\mathrm{\Σ}}}{b}_m\underset{i=1}{\overset{3}{\mathrm{\Σ}}}{({\stackrel{\·}{U}}_i^{\left(m\right)}-{\stackrel{\·}{u}}_i)}^2{D}_L;$

其中是耗散系数，η_m,κ_m,v_m分别是第m层区域的流体粘性、固体骨架渗透率和第m层区域所占体积比。D_L为局域流振荡引起的耗散函数，其形式为：

$\begin{array}{c}{D}_L=\frac{1}{24\mathrm{\π}}\frac{{\mathrm{\η}}_1}{k_1}{a}_{10}^2{\mathrm{\φ}}_{10}{\mathrm{\φ}}_1{\mathrm{\φ}}_2^2{\stackrel{\·}{\mathrm{\ζ}}}_1^{2}E+\frac{1}{18}\frac{{\mathrm{\η}}_2}{k_2}{a}_{10}^2{\mathrm{\φ}}_{10}{\mathrm{\φ}}_1{\mathrm{\φ}}_2^2(1+\frac{b_{10}^2}{a_{10}^2}+\frac{c_{10}^2}{a_{10}^2}){\stackrel{\·}{\mathrm{\ζ}}}_1^2+\\ \frac{1}{18}\{-\frac{{\mathrm{\η}}_2}{k_2}+\frac{{\mathrm{\η}}_3}{k_3}[1-{(\frac{{\mathrm{\φ}}_1}{{\mathrm{\φ}}_{10}}+\frac{{\mathrm{\φ}}_2}{{\mathrm{\φ}}_{20}})}^{1/3}\left]\right\}(1+\frac{b_{10}^2}{a_{10}^2}+\frac{c_{10}^2}{a_{10}^2})a_{20}^2(\frac{{\mathrm{\φ}}_1}{{\mathrm{\φ}}_{10}}+\frac{{\mathrm{\φ}}_2}{{\mathrm{\φ}}_{20}}){\mathrm{\φ}}_{20}^2{\mathrm{\φ}}_3^2{\stackrel{\·}{\mathrm{\ζ}}}_2^2\end{array};$

其次，根据经典力学哈密顿原理和拉格朗日方程，得到改进后的Biot-Rayleigh方程组：

$\begin{array}{c}N{\▿}^{2}u+(A+N)\▿I_1+Q_1\▿({\mathrm{\ξ}}_1-{\mathrm{\φ}}_2{\mathrm{\ζ}}_1)+Q_2\▿({\mathrm{\ξ}}_2+{\mathrm{\φ}}_1{\mathrm{\ζ}}_1-{\mathrm{\φ}}_3{\mathrm{\ζ}}_2)+Q_3\▿({\mathrm{\ξ}}_3+{\mathrm{\φ}}_2{\mathrm{\ζ}}_2)=\\ {\mathrm{\ρ}}_{00}\stackrel{\·\·}{u}+\underset{m=1}{\overset{3}{\mathrm{\Σ}}}[{\mathrm{\ρ}}_{0m}{\stackrel{\·\·}{U}}^{\left(m\right)}+b_m(\stackrel{\·}{u}-{\stackrel{\·}{U}}^{\left(m\right)})]\end{array};;$

$Q_1\▿I_1+R_1\▿({\mathrm{\ξ}}_1-{\mathrm{\φ}}_2{\mathrm{\ζ}}_1)={\mathrm{\ρ}}_{01}\stackrel{\·\·}{u}+{\mathrm{\ρ}}_{11}{\stackrel{\·\·}{U}}^{\left(1\right)}+b_1({\stackrel{\·}{U}}^{\left(1\right)}-\stackrel{\·}{u});$

$Q_2\▿I_2+R_2\▿({\mathrm{\ξ}}_2{\mathrm{\φ}}_1{\mathrm{\ζ}}_1-{\mathrm{\φ}}_3{\mathrm{\ζ}}_2)={\mathrm{\ρ}}_{02}\stackrel{\·\·}{u}+{\mathrm{\ρ}}_{22}{\stackrel{\·\·}{U}}^{\left(2\right)}+b_2({\stackrel{\·}{U}}^{\left(2\right)}-\stackrel{\·}{u});$

$Q_3\▿I_1+R_3\▿({\mathrm{\ξ}}_3-{\mathrm{\φ}}_2{\mathrm{\ζ}}_2)={\mathrm{\ρ}}_{03}\stackrel{\·\·}{u}+{\mathrm{\ρ}}_{33}{\stackrel{\·\·}{U}}^{\left(3\right)}+b_3({\stackrel{\·}{U}}^{\left(3\right)}-\stackrel{\·}{u});$

${\mathrm{\φ}}_2{Q}_1{I}_1+{\mathrm{\φ}}_2{R}_1({\mathrm{\ξ}}_1-{\mathrm{\φ}}_2{\mathrm{\ζ}}_1)-{\mathrm{\φ}}_1{Q}_2{I}_1-{\mathrm{\φ}}_1{R}_2({\mathrm{\ξ}}_2+{\mathrm{\φ}}_1{\mathrm{\ζ}}_1-{\mathrm{\φ}}_3{\mathrm{\ζ}}_2)=-{\mathrm{\α}}_1{\stackrel{\·\·}{\mathrm{\ζ}}}_1-{\mathrm{\β}}_1{\stackrel{\·}{\mathrm{\ζ}}}_1;$

${\mathrm{\φ}}_3{Q}_2{I}_1+{\mathrm{\φ}}_3{R}_2({\mathrm{\ξ}}_2-{\mathrm{\φ}}_1{\mathrm{\ζ}}_1-{\mathrm{\φ}}_3{\mathrm{\ζ}}_2)-{\mathrm{\φ}}_2{Q}_3{I}_1-{\mathrm{\φ}}_2{R}_3({\mathrm{\ξ}}_3+{\mathrm{\φ}}_2{\mathrm{\ζ}}_2)=-{\mathrm{\α}}_2{\stackrel{\·\·}{\mathrm{\ζ}}}_2-{\mathrm{\β}}_2{\stackrel{\·}{\mathrm{\ζ}}}_2;$

其中：

${\mathrm{\α}}_1=\frac{{\mathrm{\ρ}}_{f_1}}{12\mathrm{\π}}{a}_{10}^2{\mathrm{\φ}}_1{\mathrm{\φ}}_2^{2}E+\frac{{\mathrm{\ρ}}_{f_2}}{9}\frac{{\mathrm{\φ}}_{10}}{{\mathrm{\φ}}_{20}}{\mathrm{\φ}}_1{\mathrm{\φ}}_2^2{a}_{10}^2(1+\frac{b_{10}^2}{a_{10}^2}+\frac{c_{10}^2}{a_{10}^2});$

${\mathrm{\β}}_1=\frac{1}{12\mathrm{\π}}\frac{{\mathrm{\η}}_1}{k_1}{a}_{10}^2{\mathrm{\φ}}_{10}{\mathrm{\φ}}_1{\mathrm{\φ}}_2^{2}E+\frac{1}{9}\frac{{\mathrm{\η}}_2}{k_2}{a}_{10}^2{\mathrm{\φ}}_{10}{\mathrm{\φ}}_1{\mathrm{\φ}}_2^2(1+\frac{b_{10}^2}{a_{10}^2}+\frac{c_{10}^2}{a_{10}^2});$

${\mathrm{\α}}_2=\frac{1}{9}\{-{\mathrm{\ρ}}_{f_2}+{\mathrm{\ρ}}_{f_3}\frac{{\mathrm{\φ}}_{20}}{{\mathrm{\φ}}_{30}}[1-{(\frac{{\mathrm{\φ}}_1}{{\mathrm{\φ}}_{10}}+\frac{{\mathrm{\φ}}_2}{{\mathrm{\φ}}_{20}})}^{1/3}\left]\right\}(1+\frac{b_{10}^2}{a_{10}^2}+\frac{c_{10}^2}{a_{10}^2})a_{20}^2{\mathrm{\φ}}_{20}{\mathrm{\φ}}_3^2(\frac{{\mathrm{\φ}}_1}{{\mathrm{\φ}}_{10}}+\frac{{\mathrm{\φ}}_2}{{\mathrm{\φ}}_{20}});$

${\mathrm{\β}}_2=\frac{1}{9}\{-\frac{{\mathrm{\η}}_2}{k_2}+\frac{{\mathrm{\η}}_3}{k_3}[1-{(\frac{{\mathrm{\φ}}_1}{{\mathrm{\φ}}_{10}}+\frac{{\mathrm{\φ}}_2}{{\mathrm{\φ}}_{20}})}^{1/3}\left]\right\}(1+\frac{b_{10}^2}{a_{10}^2}+\frac{c_{10}^2}{a_{10}^2})a_{20}^2(\frac{{\mathrm{\φ}}_1}{{\mathrm{\φ}}_{10}}+\frac{{\mathrm{\φ}}_2}{{\mathrm{\φ}}_{20}}){\mathrm{\φ}}_{20}^2{\mathrm{\φ}}_3^2.$

在步骤S30中，根据平面波分析，将e^i(ωt-kx)代入波动方程即改进后的Biot-Rayleigh方程组，其中ω为角频率，k为波数，得到：

$\left|\begin{array}{cccc}{a}_{11}{k}^2+b_{11}& a_{12}{k}^2+b_{12}& a_{13}{k}^2+b_{13}& a_{14}{k}^2+b_{14}\\ a_{21}{k}^2+b_{21}& a_{22}{k}^2+b_{22}& a_{23}{k}^2+b_{23}& a_{24}{k}^2+b_{24}\\ a_{31}{k}^2+b_{31}& a_{32}{k}^2+b_{32}& a_{33}{k}^2+b_{33}& a_{34}{k}^2+b_{34}\\ a_{41}{k}^2+b_{41}& a_{42}{k}^2+b_{42}& a_{43}{k}^2+b_{43}& a_{44}{k}^2+b_{44}\end{array}\right|=0;$

其中，a_ij,b_ij是方程系数。从中求取波数k和角频率ω的解，得到纵波波速和逆品质因子来描述衰减和耗散。纵波速度频散的预测公式为衰减预测公式Im和Re表示虚部和实部。将固体骨架和流体参数带入纵波速度的预测公式和衰减预测公式可得到纵波速度频散和衰减变化曲线。

进一步，根据实际情况通过改变椭球体主轴半径比例，可以预测不同注水和排干情况下对应的储层岩石纵波速度变化情况。根据观测数据和预测数据的拟合关系，还可以反演得到不同含水饱和度下斑块形状变化情况。

本发明通过引入椭球形斑块饱和模型，并且包含了固体骨架和流体的空间非均匀性，因此其速度预测更具灵活性和贴近实际的特点，可以解决较为复杂的多流体、多孔隙类型岩石纵波预测问题。该方法进一步拓展了以斑块模型为基础的孔隙介质速度预测方法，其有益效果主要体现在如下方面：

首先，通过建立更具有一般性的椭球形流体饱和模型，提高了模拟实际岩石流体饱和斑块的真实性。椭球体也可以通过调整三条主轴的比例关系，退化为Patchy模型和Brown的Penny模型。

其次，通过考虑双孔介质中嵌入体内部的流体速度场，使模型不再局限于嵌入体内为气体的情况。

再次，将岩石干骨架孔隙度非均匀性和流体非均匀性结合起来，与单一流体的双重孔隙模型相比更具有一般性。

以下结合具体实施例详细说明本发明。

实施例一椭球形斑块部分饱和模型预测低孔隙砂岩低频纵波的速度

本实施例采用1984年发表的低频低孔隙非饱和砂岩的波速观测数据(Murphy,Acoustic Measures of Partial Gas Saturation in Tight Sandstones,JOURNAL OF GEOPHYSICAL RESEARCH,1984)，与本发明预测的结果进行对比分析，同时还与其他三种球形斑块部分饱和模型的预测结果进行了对比。

Fort Union砂岩是河流杂砂岩，其中包含石英约从65至95％。这些岩石包含超过25％的不稳定的材料（即长石和岩石碎片），砾石含量比长石多，具有较多空隙或矿物碎石（通常是碳酸盐）。Fort Union砂岩晶粒直径在0.125至0.15毫米之间，岩石参数为：基质体积模量35GPa，骨架体积模量7.14Gpa，骨架剪切模量9.06Gpa，水体积模量2.25GPa，空气体积模量0.8MPa，水黏度0.001Pa*s，空气黏度0.00001Pa*s，基质平均密度2.65g/cm3，水密度0.997g/cm3，空气密度0.1g/cm3，孔隙度为0.085，渗透率0.5mD。

在椭球形斑块部分饱和模型中，选取斑块主半径平均尺寸a=1.2mm，另外两个椭球主轴半径分别为b=a·r_ba,c=a·r_ca，其中r_ba,r_ca是比例系数。通过设置三个主轴半径比例关系，可以模拟具有不同体积/面积比V/A的饱和情况。在模型中，低孔隙度掺杂体所占比例为25%，内部低孔掺杂的孔隙度是平均孔隙度(0.085)的50%，内部渗透率是平均渗透率（0.5mD）的1%。

对于Fort Union砂岩低频范围(5kHz)声波实验数据，分别采用了White-Dutta、Johnson、Biot-Rayleigh、球形三层饱和斑块模型和椭球形三层饱和斑块模型来预测纵波速度。得到如下结果（见图5）:(1)White-Dutta、Johnson、Biot-Rayleigh方法的速度在BGW和BGH理论预测范围内，(2)White-Dutta和Johnson的预测两者近似，Biot-Rayleigh方法更接近实验数值，采用椭球形三层饱和斑块模型的TLP方法预测的速度可以突破理论范围，更符合实验数据。

在三层饱和斑块模型基础上，本发明提出的椭球体斑块部分饱和模型可以通过调整三个主轴半径比例关系，考察局部流体饱和斑块形状对波速频散和衰减的影响。在其他参数保持不变的情况下，改变第二主轴与第一主轴半径比例系数r_ba，发现速度预测能够更符合实验观察(图6)，当r_ba=0.005时，纵波速度与实验观察值吻合的很好。

这里r_ba=0.005是不同水饱和情况下斑块形状的等效平均值。实际上，斑块形状会随含水饱和度发生改变。如果考虑到岩石样本注水和干燥过程的不同，在相同含水饱和度情况下斑块形状也会不一样。另外，如果在实验中岩石边界的封闭和开放条件不一样，也会影响斑块形状的变化。因此可以断定，含局部饱和的斑块形状肯定不会一成不变，更不会是规则的球形，这一点已经得到Y.Tserkovnyak和D.L.Johnson的研究证实(Dvorkin and others1995)。

基于本发明，可以利用实验观察纵波速度随饱和度的变化，反演出不同含水饱和度下斑块的形状。图6给出了经过实验观察速度拟合的椭球斑块部分饱和模型速度预测曲线。这条曲线上不同饱和度对应的每个速度点具有不同的主轴半径比例r_ba，而不是使用上面提到的等效平均值。这表明本发明提出的三层椭球斑块部分饱和模型能够更精准的刻画孔隙岩石非均性特征，较其他方法具有明显优势。

实施例二预测高渗透率北海砂岩的地震波、声波和超声波速度

本实施例采用北海砂岩的波速观测数据(Pride and others2004;White1975)，计算了不同频率范围的波速，与包括本发明在内的三种斑块部分饱和模型预测结果进行了对比分析。

该北海砂岩样本的矿物构成包含80%石英，15%长石，5%粘土。砂岩晶粒直径约在0.1左右，岩石参数为(Boruah and Chatterjee2010)：基质体积模量39.47GPa，骨架体积模量5.33Gpa，骨架剪切模量3.54Gpa，盐水体积模量2.48GPa，空气体积模量0.01MPa，水黏度0.0011Pa·s，空气黏度1.81x10^-5Pa·s，基质平均密度2.63g/cm3，盐水密度1.06g/cm3，空气密度1.2x10^-3g/cm3，孔隙度为0.35，渗透率8.7D。

不同水饱和度下，斑块的尺寸会发生变化(Dvorkin and others1995)，斑块尺寸和孔隙中的慢波(扩散波)波长之间的相对关系，是影响频散和衰减的重要因素。在椭球形斑块部分饱和模型中，我们的选取斑块主半径平均尺寸a为慢波波长的1/5，另外两个椭球主轴半径分别为b=a·r_ba,c=a·r_ca，其中r_ba,r_ca是比例系数。通过设置三个主轴半径比例关系，可以模拟具有不同体积/面积比V/A的饱和情况。在模型中，外部低孔隙度掺杂体所占比例为80%，其孔隙度是平均孔隙度（0.35）的65%，密度为平均密度(2.63g/cm3)的85%。

图7中给出了岩石样本纵波速度的实验数据和模型预测结果对比。其中五组实验数据覆盖了不同的频率范围：(1)5-50Ｈz;(2)75-200Hz;(3)0.3-1kHz;(4)1-2.5kHz;(5)超声波（ultrasonic）。对于低频范围(50Hz)地震波数据，还利用Gassmann模型进行了预测。

对于不同频率范围的实验数据，分别采用了White、Johnson和椭球形三层部分饱和斑块模型来预测纵波速度，并做了对比分析。得到如下结果:

(1)White、Johnson方法的速度在BGW和BGH理论预测范围内(图7(a));

(2)Whitea和Johnson的预测两者近似(图7(a));

(3)采用椭球形三层饱和斑块模型的TLP方法预测的速度可以突破理论范围(图7(b))，在其他参数保持不变的情况下，改变第二主轴与第一主轴半径比例系数r_ba，当r_ba=0.001时，我们发现速度预测能够更接近实验测量数据(图7)。

数据对比分析表明，本发明提出的三层椭球斑块部分饱和模型，具有斑块几何形状可变的特点，能够更精准的刻画孔隙岩石非均性特征，较其他方法具有明显优势。

另外，本领域技术人员还可在本发明精神内作其它变化，当然这些依据本发明精神所作的变化，都应包含在本发明所要求保护的范围内。

Claims (10)

1.一种基于椭球体双重孔隙模型的岩石纵波速度预测方法，包括以下步骤：步骤S10，通过获取测井数据和实验观测数据，得到岩石物理参数，并生成岩石的干骨架模型和孔隙流体模型；

步骤S20，以岩石的干骨架模型和孔隙流体模型为基础，建立双重孔隙介质椭球斑块饱和模型，所述双重孔隙介质椭球斑块饱和模型包括三层具有不同孔隙度和不同流体饱和的椭球壳区域模型，每层球壳代表具有不同孔隙特征和不同流体特征的区域，并计算孔隙介质的势能/动能，以及内嵌入体中流体的动能和耗散方程，推导出拉格朗日方程组，并求取纵横波速度；

步骤S30，根据平面波分析方法，获得波动方程的频散关系，并得到纵波速度频散和衰减计算公式。

2.如权利要求1所述的基于椭球体双重孔隙模型的岩石纵波速度预测方法，其特征在于，所述岩石物理参数包括岩石的矿物成分、矿物体积比率、渗透率、孔隙率、泥质含量，以及孔隙中的流体数据包括流体的密度、粘性、弹性模量。

3.如权利要求2所述的基于椭球体双重孔隙模型的岩石纵波速度预测方法，其特征在于，岩石物理参数由目标区域的地质报告、测井资料、岩芯切片、地层的温度、压力、矿化度信息获得。

4.如权利要求3所述的基于椭球体双重孔隙模型的岩石纵波速度预测方法，其特征在于，所述岩石本身的等效弹性模量采用Voigt-Reuss-Hill平均模型计算：

$M_{\mathrm{VRH}}=\frac{1}{2}(\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{v}_i{M}_i+1/\underset{i=1}{\overset{n}{\mathrm{\Σ}}}\frac{v_i}{M_i});$

其中：

M_VRH：岩石本身的等效弹性模量；

v_i：第i种矿物的体积率；

M_i：第i种矿物的弹性模量；

n：岩石中矿物的总种类数；

采用Pride半经验性的公式来计算干骨架的等效体积模量和剪切模量:

$K_b=\frac{1-\mathrm{\φ}}{1+\mathrm{c\φ}}{K}_s,{\mathrm{\μ}}_b=\frac{1-\mathrm{\φ}}{1+c\text{'}\mathrm{\φ}}{\mathrm{\μ}}_s;$

K_b,μ_b:岩石干骨架的体积模量和剪切模量；

K_s,μ_s:岩石基质的体积模量和剪切模量；

φ:岩石的孔隙度；

c，c′:为经验性参数，与岩石的固结程度有关。

5.如权利要求1所述的基于椭球体双重孔隙模型的岩石纵波速度预测方法，其特征在于，所述椭球斑块饱和模型包含双重孔隙和双重流体，以模拟固体骨架非均匀性和流体非均匀性同时存在的情况。

6.如权利要求5所述的基于椭球体双重孔隙模型的岩石纵波速度预测方法，其特征在于，所述椭球体斑块模型的波动方程建立包括以下步骤：

步骤S21，获取椭球形嵌入体内部的流体速度特征，对双重孔隙介质中流体动能函数、耗散函数进行求取。

步骤S22，通过哈密顿原理和拉格朗日方程，建立双重孔隙介质的波动方程。

7.如权利要求6所述的基于椭球体双重孔隙模型的岩石纵波速度预测方法，其特征在于，所述椭球斑块饱和系统为一三层椭球体斑块饱和模型，其动能函数为：

$T=\frac{1}{2}\underset{i=1}{\overset{3}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{00}{\stackrel{\·}{u}}_i^2+\underset{i=1}{\overset{3}{\mathrm{\Σ}}}\underset{m=1}{\overset{3}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{0m}{\stackrel{\·}{u}}_i{\stackrel{\·}{U}}_i^{\left(m\right)}+\frac{1}{2}\underset{i=1}{\overset{3}{\mathrm{\Σ}}}\underset{m=1}{\overset{3}{\mathrm{\Σ}}}{\mathrm{\ρ}}_{\mathrm{mm}}{\stackrel{\·}{U}}_i^{{\left(m\right)}^2}+T_L;$

其中m=1,2,3表示三层球壳的不同区域，表示第m类孔隙区域内的流体空间位移分量，u₁,u₂,u₃表示固体位移分量， ${\mathrm{\ρ}}_{00}={\mathrm{\ρ}}_0-\underset{m=1}{\overset{3}{\mathrm{\Σ}}}({\mathrm{\ρ}}_m-{\mathrm{\ρ}}_{\mathrm{mm}}),$ ρ_0m=ρ_m-ρ_mm，ρ_m=φ_mρ_fm, ${\mathrm{\ρ}}_{\mathrm{mm}}={\stackrel{\‾}{\mathrm{\ρ}}}_m{\mathrm{\φ}}_m^2,$ ${\stackrel{\‾}{\mathrm{\ρ}}}_m={\mathrm{\ρ}}_{\mathrm{fm}}{a}^2;$ ρ_fm表示三层斑块模型中第m区域的流体密度，a是与孔隙几何特征有关的系数；φ为岩石总的孔隙度，v_m分别表示分别表示三个区域所占体积比，φ_m分别表示三个区域内部的局部孔隙度；

T_L是双重孔隙之间的差异引起的局部流动动能：

$\begin{array}{c}{T}_L=\frac{1}{24\mathrm{\π}}{\mathrm{\ρ}}_{f_1}{a}_{10}^2{\mathrm{\φ}}_1{\mathrm{\φ}}_2^2+{\stackrel{\·}{\mathrm{\ζ}}}_1^{2}E+\frac{{\mathrm{\ρ}}_{f_2}}{18}\frac{{\mathrm{\φ}}_{10}}{{\mathrm{\φ}}_{20}}{a}_{10}^2{\mathrm{\φ}}_1{\mathrm{\φ}}_2^2(1+\frac{b_{10}^2}{a_{10}^2}+\frac{c_{10}^2}{a_{10}^2}){\stackrel{\·}{\mathrm{\ζ}}}_1^2+\\ \{-{\mathrm{\ρ}}_{f_2}+{\mathrm{\ρ}}_{f_3}\frac{{\mathrm{\φ}}_{20}}{{\mathrm{\φ}}_{30}}[1-{(\frac{{\mathrm{\φ}}_1}{{\mathrm{\φ}}_{10}}+\frac{{\mathrm{\φ}}_2}{{\mathrm{\φ}}_{20}})}^{1/3}\left]\right\}\frac{1}{18}(1+\frac{b_{10}^2}{a_{10}^2}+\frac{c_{10}^2}{a_{10}^2})a_{20}^2{\mathrm{\φ}}_{20}{\mathrm{\φ}}_3^2(\frac{{\mathrm{\φ}}_1}{{\mathrm{\φ}}_{10}}+\frac{{\mathrm{\φ}}_2}{{\mathrm{\φ}}_{20}}){\stackrel{\·}{\mathrm{\ζ}}}_2^2\end{array};;$

其中，ζ为局域流流动造成的体应变增量，表示随时间的变化率，a₁₀,b₁₀,c₁₀是椭球体初始主轴半径，E是与椭球体积分有关的项：

8.如权利要求6所述的基于椭球体双重孔隙模型的岩石纵波速度预测方法，其特征在于，所述双重孔隙介质的耗散函数为：

$D=\frac{1}{2}\underset{m=1}{\overset{3}{\mathrm{\Σ}}}{b}_m\underset{i=1}{\overset{3}{\mathrm{\Σ}}}{({\stackrel{\·}{U}}_i^{\left(m\right)}-{\stackrel{\·}{u}}_i)}^2{D}_L;$

其中是耗散系数，η_m,κ_m,v_m分别是第m层区域的流体粘性、固体骨架渗透率和第m层区域所占体积比；D_L为局域流振荡引起的耗散函数，其形式为：

$\begin{array}{c}{D}_L=\frac{1}{24\mathrm{\π}}\frac{{\mathrm{\η}}_1}{k_1}{a}_{10}^2{\mathrm{\φ}}_{10}{\mathrm{\φ}}_1{\mathrm{\φ}}_2^2{\stackrel{\·}{\mathrm{\ζ}}}_1^{2}E+\frac{1}{18}\frac{{\mathrm{\η}}_2}{k_2}{a}_{10}^2{\mathrm{\φ}}_{10}{\mathrm{\φ}}_1{\mathrm{\φ}}_2^2(1+\frac{b_{10}^2}{a_{10}^2}+\frac{c_{10}^2}{a_{10}^2}){\stackrel{\·}{\mathrm{\ζ}}}_1^2+\\ \frac{1}{18}\{-\frac{{\mathrm{\η}}_2}{k_2}+\frac{{\mathrm{\η}}_3}{k_3}[1-{(\frac{{\mathrm{\φ}}_1}{{\mathrm{\φ}}_{10}}+\frac{{\mathrm{\φ}}_2}{{\mathrm{\φ}}_{20}})}^{1/3}\left]\right\}(1+\frac{b_{10}^2}{a_{10}^2}+\frac{c_{10}^2}{a_{10}^2})a_{20}^2(\frac{{\mathrm{\φ}}_1}{{\mathrm{\φ}}_{10}}+\frac{{\mathrm{\φ}}_2}{{\mathrm{\φ}}_{20}}){\mathrm{\φ}}_{20}^2{\mathrm{\φ}}_3^2{\stackrel{\·}{\mathrm{\ζ}}}_2^2\end{array}.$

9.如权利要求8所述的基于椭球体双重空隙模型的岩石纵波速度的预测方法，其特征在于，所述波动方程即Biot-Rayleigh方程组为：

$\begin{array}{c}N{\▿}^{2}u+(A+N)\▿I_1+Q_1\▿({\mathrm{\ξ}}_1-{\mathrm{\φ}}_2{\mathrm{\ζ}}_1)+Q_2\▿({\mathrm{\ξ}}_2+{\mathrm{\φ}}_1{\mathrm{\ζ}}_1-{\mathrm{\φ}}_3{\mathrm{\ζ}}_2)+Q_3\▿({\mathrm{\ξ}}_3+{\mathrm{\φ}}_2{\mathrm{\ζ}}_2)=\\ {\mathrm{\ρ}}_{00}\stackrel{\·\·}{u}+\underset{m=1}{\overset{3}{\mathrm{\Σ}}}[{\mathrm{\ρ}}_{0m}{\stackrel{\·\·}{U}}^{\left(m\right)}+b_m(\stackrel{\·}{u}-{\stackrel{\·}{U}}^{\left(m\right)})]\end{array};;$

$Q_1\▿I_1+R_1\▿({\mathrm{\ξ}}_1-{\mathrm{\φ}}_2{\mathrm{\ζ}}_1)={\mathrm{\ρ}}_{01}\stackrel{\·\·}{u}+{\mathrm{\ρ}}_{11}{\stackrel{\·\·}{U}}^{\left(1\right)}+b_1({\stackrel{\·}{U}}^{\left(1\right)}-\stackrel{\·}{u});$ $Q_2\▿I_2+R_2\▿({\mathrm{\ξ}}_2{\mathrm{\φ}}_1{\mathrm{\ζ}}_1-{\mathrm{\φ}}_3{\mathrm{\ζ}}_2)={\mathrm{\ρ}}_{02}\stackrel{\·\·}{u}+{\mathrm{\ρ}}_{22}{\stackrel{\·\·}{U}}^{\left(2\right)}+b_2({\stackrel{\·}{U}}^{\left(2\right)}-\stackrel{\·}{u});$

$Q_3\▿I_1+R_3\▿({\mathrm{\ξ}}_3-{\mathrm{\φ}}_2{\mathrm{\ζ}}_2)={\mathrm{\ρ}}_{03}\stackrel{\·\·}{u}+{\mathrm{\ρ}}_{33}{\stackrel{\·\·}{U}}^{\left(3\right)}+b_3({\stackrel{\·}{U}}^{\left(3\right)}-\stackrel{\·}{u});$ ${\mathrm{\φ}}_2{Q}_1{I}_1+{\mathrm{\φ}}_2{R}_1({\mathrm{\ξ}}_1-{\mathrm{\φ}}_2{\mathrm{\ζ}}_1)-{\mathrm{\φ}}_1{Q}_2{I}_1-{\mathrm{\φ}}_1{R}_2({\mathrm{\ξ}}_2+{\mathrm{\φ}}_1{\mathrm{\ζ}}_1-{\mathrm{\φ}}_3{\mathrm{\ζ}}_2)=-{\mathrm{\α}}_1{\stackrel{\·\·}{\mathrm{\ζ}}}_1-{\mathrm{\β}}_1{\stackrel{\·}{\mathrm{\ζ}}}_1;$

${\mathrm{\φ}}_3{Q}_2{I}_1+{\mathrm{\φ}}_3{R}_2({\mathrm{\ξ}}_2-{\mathrm{\φ}}_1{\mathrm{\ζ}}_1-{\mathrm{\φ}}_3{\mathrm{\ζ}}_2)-{\mathrm{\φ}}_2{Q}_3{I}_1-{\mathrm{\φ}}_2{R}_3({\mathrm{\ξ}}_3+{\mathrm{\φ}}_2{\mathrm{\ζ}}_2)=-{\mathrm{\α}}_2{\stackrel{\·\·}{\mathrm{\ζ}}}_2-{\mathrm{\β}}_2{\stackrel{\·}{\mathrm{\ζ}}}_2;$

其中：

${\mathrm{\α}}_1=\frac{{\mathrm{\ρ}}_{f_1}}{12\mathrm{\π}}{a}_{10}^2{\mathrm{\φ}}_1{\mathrm{\φ}}_2^{2}E+\frac{{\mathrm{\ρ}}_{f_2}}{9}\frac{{\mathrm{\φ}}_{10}}{{\mathrm{\φ}}_{20}}{\mathrm{\φ}}_1{\mathrm{\φ}}_2^2{a}_{10}^2(1+\frac{b_{10}^2}{a_{10}^2}+\frac{c_{10}^2}{a_{10}^2});$

${\mathrm{\β}}_1=\frac{1}{12\mathrm{\π}}\frac{{\mathrm{\η}}_1}{k_1}{a}_{10}^2{\mathrm{\φ}}_{10}{\mathrm{\φ}}_1{\mathrm{\φ}}_2^{2}E+\frac{1}{9}\frac{{\mathrm{\η}}_2}{k_2}{a}_{10}^2{\mathrm{\φ}}_{10}{\mathrm{\φ}}_1{\mathrm{\φ}}_2^2(1+\frac{b_{10}^2}{a_{10}^2}+\frac{c_{10}^2}{a_{10}^2});$

${\mathrm{\α}}_2=\frac{1}{9}\{-{\mathrm{\ρ}}_{f_2}+{\mathrm{\ρ}}_{f_3}\frac{{\mathrm{\φ}}_{20}}{{\mathrm{\φ}}_{30}}[1-{(\frac{{\mathrm{\φ}}_1}{{\mathrm{\φ}}_{10}}+\frac{{\mathrm{\φ}}_2}{{\mathrm{\φ}}_{20}})}^{1/3}\left]\right\}(1+\frac{b_{10}^2}{a_{10}^2}+\frac{c_{10}^2}{a_{10}^2})a_{20}^2{\mathrm{\φ}}_{20}{\mathrm{\φ}}_3^2(\frac{{\mathrm{\φ}}_1}{{\mathrm{\φ}}_{10}}+\frac{{\mathrm{\φ}}_2}{{\mathrm{\φ}}_{20}});$

${\mathrm{\β}}_2=\frac{1}{9}\{-\frac{{\mathrm{\η}}_2}{k_2}+\frac{{\mathrm{\η}}_3}{k_3}[1-{(\frac{{\mathrm{\φ}}_1}{{\mathrm{\φ}}_{10}}+\frac{{\mathrm{\φ}}_2}{{\mathrm{\φ}}_{20}})}^{1/3}\left]\right\}(1+\frac{b_{10}^2}{a_{10}^2}+\frac{c_{10}^2}{a_{10}^2})a_{20}^2(\frac{{\mathrm{\φ}}_1}{{\mathrm{\φ}}_{10}}+\frac{{\mathrm{\φ}}_2}{{\mathrm{\φ}}_{20}}){\mathrm{\φ}}_{20}^2{\mathrm{\φ}}_3^2.$

10.如权利要求9所述的基于椭球体双重孔隙模型的岩石纵波速度预测方法，其特征在于，通过将e^i(ωt-kx)代入波动方程获得纵波速度频散及纵波速度频散和衰减，所述纵波速度频散的预测公式为衰减预测公式Im和Re表示虚部和实部，k为波数和ω为角频率。