一种孔裂隙介质中纵波速度与衰减的预测方法与流程

本发明属于地震岩石物理领域，特别涉及一种孔裂隙介质中纵波速度与衰减的预测方法。

背景技术：

地下岩石往往具有孔隙与裂隙并存的特征，裂隙的存在不仅仅影响岩石的弹性性质，同时控制孔裂隙间的孔隙流体流动，也被称为喷射流，(müller等，2010；carcione，2014)，尤其是这种局部流体流动(简称局域流)是造成孔裂隙介质中地震波频散和衰减的主要原因。而地震波的频散和衰减又蕴含着丰富的地下岩性，岩石孔隙结构与孔隙流体信息等(quintal等，2011；yao等，2015；khalid等，2016)，因此，开展孔裂隙介质中地震波传播特征的研究，有利于地震解释与地下储层的识别。

国内外众多学者通过等效介质理论(eshelby，1957；walsh，1965；mori和tanaka，1973；kuster和1974；berryman，1980；norris，1985；thomsen，1985；song等，2016)与岩石物理实验(fortin等，2007；adam等，2013；wang，2018；yin，2018)研究裂隙对岩石弹性性质的影响，其中，david和zimmerman(2011)计算了岩石中包含随机分布的椭球状裂隙的弹性模量。由于裂隙在岩石中所占含量较小(通常小于0.1％)，但对弹性性质的影响比较提出。考虑到岩石中裂隙的特征与压力变化紧紧相关，岩石的裂隙的分布特征(裂隙纵横比，裂隙密度，裂隙孔隙度，裂隙半径)可从与压力相关的弹性模量中进行提取(cheng和1979；tran等，2008；david等，2012)。尽管这些方法可估算岩石内部裂隙分布特征，但并没有考虑孔裂隙间的孔隙流体流动。

为解决孔裂隙间的孔隙流体流动，最经典的方法是利用流体力学研究孔裂隙间的孔隙流动流动，提出众多模型以解释该现象(mavko和nur，1975；murphy等，1986；gurevich等，2009；carcione和gurevich，2011)，但这些模型中并没有体现裂隙的特征(裂隙密度，裂隙纵横比与裂隙半径)并且计算过程复杂。唐晓明(2011,2012)将裂隙密度与裂隙纵横比引入以解释孔裂隙间的孔隙流体流动(以下简称为唐模型)，发现裂隙纵横比影响地震波传播的驰豫频率，而裂隙密度影响地震波频散和衰减的大小。上诉模型中预测的地震波速度在低频极限无法与gassmann理论预测结果重合(gassmann理论在孔隙介质理论中被普遍认为是岩石弹性模量的下限)，yao等(2015)通过将孔裂隙间的孔隙流体流动引入到流体模量中(以下简称为dfm模型)，已解决该问题，但该方法预测的地震波频散和衰减相对于唐模型预测的要略小。

综上所诉，现有的技术研究中仍存在一下问题：

(1)实验室室内岩石内部微观孔隙结构特征的研究需考虑孔裂隙间孔隙流体流动的影响；

(2)孔裂隙介质中地震波传播理论的建立需要充分考虑裂隙的特征(裂隙密度，裂隙纵横比以及裂隙长度半径)；

(3)地震波的预测结果低频极限需与gassmann理论预测结果重合。

技术实现要素：

发明目的：为了克服现有技术中存在的问题，本发明提供孔裂隙介质中纵波速度与衰减的预测方法，利用硬币型的嵌入体模拟孔裂隙间的孔隙流体流动，进而研究孔裂隙介质中地震波的传播特征。

技术方案：为实现上述目的，本发明提供一种孔裂隙介质中纵波速度与衰减的预测方法，包括如下步骤：

(1)获取随压力变化的超声波速度实验数据；

(2)建立符合孔裂隙储层微观孔隙结构特征的等效介质理论，结合所获得的实验数据，求取岩石内部的微观孔隙结构特征；

(3)推导孔裂隙介质的应力应变关系；

(4)推导孔裂隙介质的运动方程；

(5)通过将裂隙模拟为硬币型的嵌入体，以模拟孔裂隙间的孔隙流体流动，推导孔裂隙介质的局域流控制方程；

(6)利用平面波分析，预测孔裂隙介质的纵波速度和衰减；

(7)通过得到的孔裂隙介质的纵波速度和衰减,可预测地下储层的油气分布规律，进一步指导地下油气储层的勘探与开发。

进一步的，所述步骤(2)中求取岩石内部的微观孔隙结构特征的具体步骤如下：

(2.1)基于mori-tanaka理论建立岩石弹性模量与微观孔隙结构的关系，则岩石介质的等效弹性模量表达为：

其中kstiff和gstiff分别为仅含硬孔隙岩石的等效体积与剪切模量，k0和g0分别为岩石颗粒的体积与剪切模量，φstiff为岩石中硬孔隙的孔隙度，p、q分别为硬孔隙的形状因子，与椭球形孔隙的纵横比α以及岩石颗粒的泊松比v有关，其中纵横比定义为裂缝宽度一半与长度一半的比，定义为：

其中，v＝(3k0-2g0)/(6k0+2g0),

以包含硬孔隙的介质为背景相，考虑裂隙对岩石弹性性质的影响，则岩石介质的等效弹性模量将表达为：

其中，keff、geff分别为岩石的等效体积与剪切模量，

vstiff＝(3kstiff-2gstiff)/(6kstiff+2gstiff)是仅含硬孔隙岩石的泊松比，γ是裂隙密度即嵌入在单元体中的裂隙总数；

(2.2)由于高有效压力下软孔隙即裂隙几乎完全关闭，岩石中仅有硬孔隙存在，因此，通过实验室所测量的高有效压力下岩石的纵横波速度估算岩石高压下的等效体积模量k^hp与剪切模量g^hp，则有：

其中上标“hp”为高有效压力，ρ为岩石的密度，即ρ＝(1-φ)ρs+φρf，ρs和ρf分别为颗粒密度和流体密度；基于上述kstiff的公式，利用最小二乘法求取最吻合的硬孔隙的纵横比

(2.3)基于上述keff和geff的公式，通过对各个有效压力下测量的纵横波速度进行最小二乘法计算各个压力下的累积裂隙密度γp(α)；

(2.4)基于所给出各有效压力p下的裂隙密度，建立裂隙密度随有效压力变化的定量关系：

其中，γⁱ是在有效压力为零时的初始裂隙密度，为一个与压力p同数量级的压力常数；

(2.5)建立初始状态的孔隙纵横比与有效压力之间的定量关系，表达为

其中，是高有效压力下的等效杨氏模量，定义为考虑到有效压力p逐渐增大时，具有不同纵横特征的裂隙随压差增量dp的变化量是相同的；因此，当dp足够小时，裂隙密度的减少量归因于纵横比小于的裂隙闭合造成的；故可以基于实验室测量的纵横波速度获得岩石内部裂隙孔隙度与密度的分布特征。

进一步的，所述步骤(3)中推导孔裂隙介质的应力应变关系的具体步骤如下：

孔裂隙介质的应力应变关系表达为：

其中，δij为克罗内克符号，i,j表示笛卡尔坐标系中的x,y,z三个方向，τij为总应力，pfm为孔隙流体压力，其中m＝1,2分别表示背景介质与嵌入体孔隙区域；为固体应变分量，体应变为流体形变增量，w^(m)＝φm(u^(m)-u),u和u分别为流体和固体位移分量；φ为总孔隙度，φ＝φ1+φ2,φ1和φ2分别为背景相介质和嵌入体的绝对孔隙度，φm＝vmφm0，vm和φm0为孔隙区域的体积分数与局部孔隙度；为局域流流动造成的体应变增量，μ＝g0，λc，α1，α2，m1，m2为孔裂隙介质中的弹性参数，具体表达式为：

其中kf为孔隙流体的体积模量，kb1和kb2分别为背景相与嵌入体的干骨架体积模量，kb为岩石干骨架的体积模量。

进一步的，所述步骤(4)中推导孔裂隙介质的运动方程的具体步骤如下：

孔裂隙介质中的运动方程可写为：

其中索引前面的“,”表示空间方向的偏导数，变量上的“.”表示在时间方向的偏导数，η为流体黏度，κ1和κ2分别为背景相与嵌入体的渗透率，参数m1、m2被定义为：

这里τ1与τ2分别为背景相与嵌入体的曲折度。

进一步的，所述步骤(5)中得出孔裂隙介质的局域流控制方程的具体步骤如下：

采用如下公式计算孔裂隙介质中的局域流控制方程：

其中ro为硬币状嵌入体的长度半径，局域流流动的特征长度。

进一步的，所述步骤(6)中利用平面波分析，预测孔裂隙介质的纵波速度和衰减的具体步骤如下：

基于孔裂隙介质的应力应变关系，运动方程和局域流控制方程，最终得到孔裂隙介质中地震波的传播方程组；通过平面波分析，将位移场平面波解代入到所述方程组中，并将方程组转换到频率波数域，得到christoffel方程；

由复波数k生成纵波的速度预测公式和衰减其中，ω＝2πf，f为频率。

有益效果：本发明与现有技术相比具有以下优点：

本发明充分考虑了孔裂隙介质中裂隙的分布特征(裂隙密度、裂隙孔隙度、裂隙纵横比与裂隙半径)，将裂隙模拟为硬币型嵌入体，通过硬币型嵌入体的周期性振荡以模拟孔裂隙间的孔隙流体流动，建立孔裂隙介质中地震波传播的理论方程组，实现孔裂隙介质中纵波速度和衰减的预测，预测结果低频极限与gassmann理论预测结果完全重合。

附图说明

图1为本发明的总体流程图；

图2为实施例二中唐模型(双划线)、dfm模型(虚线)、gassmann理论(实线)以及本发明(点划线)预测的纵波速度和衰减曲线图；

图3为实施例三中地震频带(虚线)、测井频带(双划线)与超声波频带(实线)的纵波速度与衰减随裂隙半径的变化示意图；

图4为实施例三中纵波速度与衰减随不同裂隙孔隙度变化的示意图；图中带菱形的实线、带菱形的双划线、带菱形的虚线、双划线、虚线与点划线分别对应ε＝0.16，γ＝0.002；ε＝0.18，γ＝0.002；ε＝0.2，γ＝0.002；ε＝0.2，γ＝0.0018；ε＝0.2，γ＝0.005；ε＝0.2，γ＝0.0016；

图5为实施例四中干燥状态下样本1和样本2中纵横波速度与有效压力变化的示意图；实线为mori-tanaka理论拟合结果，圆圈表示实验室测量干燥岩石的纵波速度，菱形表示实验室测量干燥岩石的横波速度

图6为实施例四中样本1和样本2中裂隙孔隙度与密度随有效压力变化的示意图；带圆圈的实线为样本1，带菱形的实线为样本2；

图7为实施例四中利用gassmann理论，dfm模型与本发明预测样本1和样本2水饱和状态下的纵波速度随有效压力变化的示意图；圆圈为实验室测量数据，实线为本发明预测结果，带正方形的双划线为gassmann预测结果，带菱形的双划线为dfm模型预测结果。

具体实施方式

下面结合附图对本发明作更进一步的说明。

实施例一

一种孔裂隙介质中纵波速度与衰减的预测方法，包括如下步骤：

(1)获取随压力变化的超声波速度实验数据；

(2)建立符合孔裂隙储层微观孔隙结构特征的等效介质理论，结合所获得的实验数据，求取岩石内部的微观孔隙结构特征；

(3)推导孔裂隙介质的应力应变关系；

(4)推导孔裂隙介质的运动方程；

(5)通过将裂隙模拟为硬币型的嵌入体，以模拟孔裂隙间的孔隙流体流动，推导孔裂隙介质的局域流控制方程；

(6)利用平面波分析，预测孔裂隙介质的纵波速度和衰减；

(7)通过得到的孔裂隙介质的纵波速度和衰减,可预测地下储层的油气分布规律，进一步指导地下油气储层的勘探与开发。

进一步的，所述步骤(2)中求取岩石内部的微观孔隙结构特征的具体步骤如下：

(2.1)基于mori-tanaka理论建立岩石弹性模量与微观孔隙结构的关系，则岩石介质的等效弹性模量表达为：

其中kstiff和gstiff分别为仅含硬孔隙岩石的等效体积与剪切模量，k0和g0分别为岩石颗粒的体积与剪切模量，φstiff为岩石中硬孔隙的孔隙度，p、q分别为硬孔隙的形状因子，与椭球形孔隙的纵横比α以及岩石颗粒的泊松比v有关，其中纵横比定义为裂缝宽度一半与长度一半的比，定义为：

其中，v＝(3k0-2g0)/(6k0+2g0),

以包含硬孔隙的介质为背景相，考虑裂隙对岩石弹性性质的影响，则岩石介质的等效弹性模量将表达为：

其中，keff、geff分别为岩石的等效体积与剪切模量，

vstiff＝(3kstiff-2gstiff)/(6kstiff+2gstiff)是仅含硬孔隙岩石的泊松比，γ是裂隙密度即嵌入在单元体中的裂隙总数；

(2.2)由于高有效压力下软孔隙即裂隙几乎完全关闭，岩石中仅有硬孔隙存在，因此，通过实验室所测量的高有效压力下岩石的纵横波速度估算岩石高压下的等效体积模量k^hp与剪切模量g^hp，则有：

其中上标“hp”为高有效压力，ρ为岩石的密度，即ρ＝(1-φ)ρs+φρf，ρs和ρf分别为颗粒密度和流体密度；基于上述kstiff的公式，利用最小二乘法求取最吻合的硬孔隙的纵横比

(2.3)基于上述keff和geff的公式，通过对各个有效压力下测量的纵横波速度进行最小二乘法计算各个压力下的累积裂隙密度γp(α)；

(2.4)基于所给出各有效压力p下的裂隙密度，建立裂隙密度随有效压力变化的定量关系：

其中，γⁱ是在有效压力为零时的初始裂隙密度，为一个与压力p同数量级的压力常数；

(2.5)建立初始状态的孔隙纵横比与有效压力之间的定量关系，表达为

其中，是高有效压力下的等效杨氏模量，定义为考虑到有效压力p逐渐增大时，具有不同纵横特征的裂隙随压差增量dp的变化量是相同的；因此，当dp足够小时，裂隙密度的减少量归因于纵横比小于的裂隙闭合造成的；故可以基于实验室测量的纵横波速度获得岩石内部裂隙孔隙度与密度的分布特征。

进一步的，所述步骤(3)中推导孔裂隙介质的应力应变关系的具体步骤如下：

孔裂隙介质的应力应变关系表达为：

其中，δij为克罗内克符号，i,j表示笛卡尔坐标系中的x,y,z三个方向，τij为总应力，pfm为孔隙流体压力，其中m＝1,2分别表示背景介质与嵌入体孔隙区域；为固体应变分量，体应变为流体形变增量，w^(m)＝φm(u^(m)-u),u和u分别为流体和固体位移分量；φ为总孔隙度，φ＝φ1+φ2,φ1和φ2分别为背景相介质和嵌入体的绝对孔隙度，φm＝vmφm0，vm和φm0为孔隙区域的体积分数与局部孔隙度；为局域流流动造成的体应变增量，μ＝g0，λc，α1，α2，m1，m2为孔裂隙介质中的弹性参数，具体表达式为：

其中kf为孔隙流体的体积模量，kb1和kb2分别为背景相与嵌入体的干骨架体积模量，kb为岩石干骨架的体积模量。

进一步的，所述步骤(4)中推导孔裂隙介质的运动方程的具体步骤如下：

孔裂隙介质中的运动方程可写为：

其中索引前面的“,”表示空间方向的偏导数，变量上的“.”表示在时间方向的偏导数，η为流体黏度，κ1和κ2分别为背景相与嵌入体的渗透率，参数m1、m2被定义为：

这里τ1与τ2分别为背景相与嵌入体的曲折度。

进一步的，所述步骤(5)中得出孔裂隙介质的局域流控制方程的具体步骤如下：

采用如下公式计算孔裂隙介质中的局域流控制方程：

其中ro为硬币状嵌入体的长度半径，局域流流动的特征长度。

进一步的，所述步骤(6)中利用平面波分析，预测孔裂隙介质的纵波速度和衰减的具体步骤如下：

基于孔裂隙介质的应力应变关系，运动方程和局域流控制方程，最终得到孔裂隙介质中地震波的传播方程组；通过平面波分析，将位移场平面波解代入到所述方程组中，并将方程组转换到频率波数域，得到christoffel方程；

由复波数k生成纵波的速度预测公式和衰减其中，ω＝2πf，f为频率。

实施例二

本实施例设计一个孔裂隙介质模型，将唐模型、dfm模型、gassmann理论与本发明的预测结果进行对比分析，如图2所示。基本的岩石物理参数为：岩石颗粒的体积模量与剪切模量分别为37.9和32.6gpa，岩石颗粒的密度为2650；骨架的体积模量4.71gpa，剪切模量4.99gpa；水的体积模量2.25gpa，密度为1000kg/m³，黏度0.001pa·s；背景相与嵌入体的局部孔隙度φ10和φ20分别为0.25与0.32，渗透率分别为0.1d与100d，骨架模量分别为7.6与0.096gpa，裂隙密度为0.2，裂隙纵横比为0.002，裂隙半径为0.0053m。从图中可发现本说明预测的纵波速度在低频与gassmann预测结果完全吻合，并能解释喷射流造成的纵波频散和衰减。

实施例三

本实施例设计一个孔裂隙介质模型，以分析裂隙特征对纵波传播特征的影响。基本的岩石物理参数为：岩石颗粒的体积模量与剪切模量分别为37.9gpa和32.6gpa，岩石颗粒的密度为2650kg/m³；骨架的体积模量4.71gpa，骨架的剪切模量4.99gpa；水的体积模量2.25gpa，密度为1000kg/m³，黏度0.001pa·s；背景相与嵌入体的局部孔隙度φ10和φ20分别为0.25与0.32，渗透率分别为0.1d与100d，骨架模量分别为7.6与0.096gpa，图3为裂隙密度为0.2，裂隙纵横比为0.002，裂隙半径为10-⁶～0.1m，地震频带(10hz)、测井频带(1000hz)与超声波频带(1mhz)纵波速度与衰减随裂隙半径变化的示意图。从图3中可发现，裂隙孔隙度(裂隙孔隙度φ2＝4/3πεγ)保持不变时，纵波速度随着裂隙半径的增大而增大，并趋于一定值。另外，喷射流控制的衰减曲线随着裂隙半径的增大向低频方向移动。图4为裂隙半径为0.01m，不同裂隙密度与纵横比组合所预测纵波速度和衰减变化示意图：ε＝0.16，γ＝0.002；ε＝0.18，γ＝0.002；ε＝0.2，γ＝0.002；ε＝0.2，γ＝0.0018；ε＝0.2，γ＝0.005；ε＝0.2，γ＝0.0016。从图4中可发现，当裂隙纵横比保持不变时，裂隙密度越大，喷射流造成的纵波频散和衰减的幅值越大；相反，当裂隙密度保持不变时，裂隙纵横比越大，喷射流造成的纵波频散和衰减的幅值越小。另外，当裂隙密度增大或纵横比减小时，喷射流所控制的纵波速度和衰减曲线向低频方向移动。同样可以发现裂隙密度主要影响纵波频散和衰减的大小，而裂隙纵横比主要影响地震波传播的驰豫频率。

实施例四

本实施例对四川盆地须家河组致密砂岩样本进行实验数据与理论预测数据的对比：

两块样本的主要成分是石英，样本1与2的颗粒体积模量皆为39gpa，剪切模量分别为为33gpa与35gpa，密度分别为2672kg/m³与2660kg/m³，背景相的局部孔隙度分别为0.0626与0.1391，背景相的渗透率分别为0.046md与1.37md，裂隙半径分别为0.000005m与0.000035m。实验中所采用的频率为1mhz，温度为80℃，有效压力(有效压力p＝pc-pp，其中，pc和pp分别是围压和孔隙压力)变化范围5-35mpa。实验中所采用的流体为水，水的体积模量2.41gpa，水的黏度0.000938pa·s，水的密度1010kg/m³。图5为实验室测量的两块致密砂岩样本干燥条件下纵横波速度随压力变化的示意图，图中实线为基于mori-tanaka理论的拟合结果。按照步骤1与2，可以得到样本1与2中裂隙孔隙度与密度分布特征，求取不同有效压力下的裂隙孔隙度与密度，如图5所示，可发现随着压力的增大，裂缝孔隙度与密度逐渐减小，趋于0。考虑到岩石内部的微观孔隙结构干燥与饱水状态下是相同的，因此，本发明基于从干燥纵横波速度数据中提取的裂隙孔隙度与密度分布特征，结合步骤3-5所建立的孔裂隙介质地震波传播方程组，通过步骤6预测饱水状态下样本1与2的纵波速度，通过与实测速度对比发现，本发明的预测结果与实验数据吻合较好。

综上所述，本发明实施例提供的一种孔裂隙介质中纵波速度与衰减的预测方法具有以下有益效果：

(1)本发明考虑了岩石内部的微观孔隙结构，将裂隙特征(裂隙密度，裂隙纵横比与裂隙半径)引入到孔裂隙介质地震波的传播方程组中，如实施例二中，对裂隙密度、裂隙纵横比以及裂隙半径对纵波传播特征的影响进行了分析，本发明说明，裂隙对孔裂隙介质中地震波传播特征的影响比较突出；

(2)采用硬币型嵌入体周期性振荡来刻画孔裂隙间的孔隙流体流动，物理过程明确，计算方法简单易实现

(3)本发明预测的低频速度与gassmann低频极限重合，并可分析喷射流造成的纵波频散和衰减，进而指导地震勘探，查明地下储层的流体信息。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。