一种岩石微观孔隙结构测试方法与流程

本发明涉及一种岩石内部裂缝孔隙度的测试方法，尤其涉及一种利用岩石的弹性模量来测试其微观孔隙分布特征的表征测试方法。

背景技术：

开展岩石微观孔隙结构对岩石弹性模量的影响分析，这对刻画地下储层特征具有重要的指导意义(tranetal,2008)。众多学者基于实验室测量的超声波速度(动态弹性模量)和应力应变关系(静态弹性模量)来研究该影响(king,1966；nur和simmons,1969；coyner,1984；baudetal,2017)。通过研究发现动弹性模量在低有效压力(有效压力p＝pc-pp，其中，pc和pp分别是围压和孔隙压力)下变化比较明显，这是由于纵横比(裂缝的纵横比定义为裂缝宽度一半与长度一半的比)较小的裂缝随着压力增大逐渐闭合，进而导致岩石的弹性性质变化明显。从本质上讲，与压力有关的岩石弹性性质与微观孔隙结构密切相关。walsh(1965)推导了干燥岩石(孔隙压力pp为零)的体积模量(体积模量为弹性模量的一种，定义为静水压力与体应变的比)随围压pc变化的表达式，发现裂缝(软孔隙)对弹性模量的影响比粒间孔隙(硬孔隙)更加显著，并估算了岩石样品的裂缝孔隙度。chengand(1979)利用kt等效介质模型(kuster和1974)估算了不同种类岩石中微观孔隙的孔隙度分布特征，即具有不同孔隙纵横比特征的孔隙，对应不同的孔隙度。由于kt模型无法估算高孔隙度岩石的弹性模量，tran等(2008)通过引入dem等效介质模型(berryman,1992)改进了cheng和的方法。然而，最终的计算结果并不是唯一的，且结果的可靠性依赖于输入的先验模型。为解决这个问题，izumotani和onozuka(2013)利用快速模拟退火算法(一种基于概率的算法)估算了微观孔隙的孔隙度分布特征。

eberhart-phillips等(1989)基于线性和指数项的组合推导了地震波速度和有效压力之间的经验关系式。zimmerman(1991)利用该关系式近似计算了三种砂岩样品的可压缩系数，并估算了微观孔隙的孔隙度分布特征。基于zimmerman的方法，david和zimmerman(2012)利用mori-tanaka理论(mori和tanaka，1973；benveniste，1987)从干燥岩石的地震波速度中获得微观孔隙的孔隙度分布特征(以下简称dz模型)，并通过gassmann方程预测饱和岩石的纵横波速度随有效压力的变化趋势。这些方法虽易实现，但估算的结果通常与实际数据不吻合，研究发现基于动态弹性模量预测的结果往往偏小(cheng和1979；pervukhina等，2010)。

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

(1)动态弹性模量对压力的敏感性未充分考虑；

(2)动态弹性模量预测的累计裂隙密度、裂缝孔隙度远远低于岩石内部的真实情况。

技术实现要素：

发明目的：针对以上问题，本发明提出一种岩石微观孔隙结构测试方法，该方法基于静态弹性模量来提取岩石内部微观孔隙的孔隙结构分布特征，提高累计裂隙密度、裂缝孔隙度等孔隙结构的测试准确度。

技术方案：本发明所采用的技术方案是：一种岩石微观孔隙结构测试方法，包括以下步骤：

s1：基于mori-tanaka理论建立岩石弹性模量与微观孔隙结构的关系，该关系包括公式(1)、公式(2)所示的岩石弹性模量与岩石中硬孔隙的孔隙度φstiff的关系：

式中kstiff和gstiff分别为仅含硬孔隙岩石的等效体积模量与等效剪切模量，k0和g0分别为不含孔隙岩石的体积模量与剪切模量，式中假设岩石中的硬孔隙为椭球形，φstiff为岩石中硬孔隙的孔隙度，p、q分别为硬孔隙的形状因子，与椭球形硬孔隙的纵横比α有对应关系；

以及公式(3)、(4)表示的岩石弹性模量与岩石中累积裂隙密度γ的关系：

式中keff和geff分别为岩石的等效体积模量与等效剪切模量，vstiff是仅含硬孔隙岩石的泊松比vstiff＝(3kstiff-2gstiff)/(6kstiff+2gstiff)，γ是累积裂隙密度；

其中，公式(1)、公式(2)中的p、q定义为：

其中，v＝(3k0-2g0)/(6k0+2g0)，g仅与硬孔隙的纵横比α有关，其表达式为：

s2：采用实验室常规的通过应力应变关系测量体积模量的方法，测量当岩石中裂缝在有效压力的作用下闭合时对应的岩石等效静态体积模量并通过公式(1)计算得到硬孔隙的纵横比

s3：基于实验室测量各个有效压力下对应的等效体积模量keff或剪切模量geff，利用公式(3)或(4)计算，可得到各个有效压力p下对应的累积裂隙密度γp,st(α)；具体包括以下过程：

s31：基于实验室测量各个有效压力下对应的等效体积模量keff或剪切模量geff；

s32：通过s2步骤中获得的硬孔隙的纵横比利用公式(1)和公式(2)计算得到仅含硬孔隙岩石的等效体积模量kstiff与等效剪切模量gstiff；

s33：计算泊松比vstiff＝(3kstiff-2gstiff)/(6kstiff+2gstiff)；

s34：利用公式(3)或者公式(4)计算得到各个有效压力p下对应的累积裂隙密度γp,st(α)。

s4：利用裂隙密度随有效压力的变化呈指数衰减定律其中是在有效压力为零时的初始裂隙密度，为一个与压力p同数量级的压力常数，根据s3所得的各个有效压力p下所对应的累积裂隙密度γp,st(α)计算拟合该关系式中的参数，从而得到岩石的裂隙密度随有效压力的变化的拟合曲线或关系式；

s5：通过裂缝孔隙度的关系式建立裂缝孔隙度与纵横比之间的关系，从而获得岩石内部裂缝孔隙度的分布特征；具体包括以下过程：

s51：首先将实验室测量压力区间分成n等份，利用裂隙纵横比与压力之间的关系式，得到各压力下对应的裂隙纵横比；其中，所述裂隙纵横比与压力之间的关系式为：

式中是高有效压力下的等效静态杨氏模量，可由公式计算得到，泊松比由计算得到。

s52：其次通过s4得到压力与累计裂隙密度关系式，求取各压力下对应的累积裂隙密度；

s53：最后通过裂缝孔隙度的关系式建立裂缝孔隙度与纵横比之间的关系，从而获得岩石内部裂缝孔隙度的分布特征。

s6：基于s5所得到岩石内部裂缝孔隙度的分布特征，对所有的裂缝孔隙度与密度进行求和，计算得到岩石中累积裂缝孔隙度和累计裂隙密度曲线，则二者曲线上的渐进值分别为待测岩石内部的裂缝孔隙度与裂隙密度。

有益效果：本发明设计岩石微观孔隙结构提取方法，利用静态弹性模量替换动态弹性模量来估算微观孔隙的孔隙度分布特征。本发明充分考虑静态弹性模量对压力变化更敏感的特点，利用静态弹性模量替换动态弹性模量来测试岩石内部的微观孔隙结构分布特征，测试结果更加贴近岩石的真实情况，可实现岩石内部微观孔隙结构的准确刻画。

附图说明

图1是本发明基于静态弹性模量测试岩石微观孔隙结构方法的流程示意图；

图2是navajo和weber砂岩的应力-应变曲线图；

图3是navajo砂岩纵横波速度随有效压力的变化曲线；

图4是weber砂岩纵横波速度随有效压力的变化曲线；

图5是navajo和weber砂岩裂缝孔隙度的分布特征；

图6是navajo和weber砂岩的累积裂缝孔隙度分布特征；

图7是navajo和weber砂岩的累积裂隙密度分布特征；

图8是navajo和weber砂岩的硬孔隙、裂缝和总孔隙度随有效压力的变化曲线。

具体实施方式

下面结合附图和实施例对本发明的具体实施方式作进一步的详细说明。

本发明提供了一种岩石微观孔隙结构测试方法，如图1所示，具体包括以下步骤：

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

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

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

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

其中，keff、geff分别为岩石的等效体积与剪切模量，vstiff＝(3kstiff-2gstiff)/(6kstiff+2gstiff)是仅含硬孔隙岩石的泊松比，γ是裂隙密度，可由下式得到：

其中，a是裂缝的平均半径，n是嵌入在单元体(体积为v)中的裂缝总数，括号表示平均值。

步骤s2：考虑到walsh(1965)给出了孔隙的闭合压力关系式pclose＝πe0α0/(4(1-(v)²))，其中e0＝3k0(1-2v)是岩石颗粒的杨氏模量，α0是p＝0时的孔隙纵横比。则在典型砂岩(e0＝50gpa)中，纵横比为0.01的孔隙的闭合压力为500mpa，这远远超过在实验室测量中所使用的压力。因此，david和zimmerman(2012)认为硬孔隙的纵横比在0.01<α<1的范围内，软孔隙的纵横比小于0.01。

由于高有效压力下软孔隙(裂缝)几乎完全关闭，岩石中仅有硬孔隙存在。因此，可通过实验室所测量的高有效压力下岩石的等效静态体积模量基于应力应变关系计算：

其中上标“hp”为高有效压力，下标“st”表示该值是基于静态弹性模量的计算值，δv为单元体体积v的变化量。最后，基于公式(1)，采用最小二乘法，即求取的体积模量与静态体积模量之间误差的平方为最小时，对应的最吻合的纵横比α，即为硬孔隙的纵横比

本步骤中所述的高有效压力，对于navajo砂岩来讲，其范围一般取为80～100mpa，weber砂岩则一般取为60～100mpa，对于任意一种待测试的岩石来说，可以通过裂缝孔隙度随压强的变化曲线，选择当裂缝孔隙度处于0％附近时对应的有效压力值作为本步骤中，也可以直接取实验测试时所能达到的压强的最大值，例如本实验中该值为100mpa，然后利用最终拟合得到的裂缝孔隙度随压强的变化曲线来验证该压力取值是否满足要求。

步骤s3：由于岩石体积模量的压力相关性与裂隙密度联系密切，当裂隙密度已知时，可以通过公式(5)与(6)获得岩石的弹性模量，反之也可基于实验室测量的弹性模量估算岩石的裂隙密度。因此，可通过对各个有效压力下测量的弹性模量进行最小二乘法计算各个压力下的累积裂隙密度γp,st(α)。具体包括以下过程：

s31：基于实验室测量各个有效压力下对应的等效体积模量keff或剪切模量geff；

s32：通过s2步骤中获得的硬孔隙的纵横比利用公式(1)和公式(2)计算得到仅含硬孔隙岩石的等效体积模量kstiff与等效剪切模量gstiff；

s33：计算泊松比vstiff＝(3kstiff-2gstiff)/(6kstiff+2gstiff)；

s34：利用公式(3)或者公式(4)计算，采用最小二乘法，当利用公式(3)或公式(4)计算的等效体积模量或剪切模量与实验测量值之间的误差平方最小时，得到对应的累积裂隙密度γp,st(α)。

步骤s4：国内外学者提出裂隙密度随有效压力的变化遵循指数衰减定律(如，shapiro，2003)：

其中，是在有效压力为零时的初始裂隙密度，为一个与压力p同数量级的压力常数。基于步骤s3中所给出各有效压力p下的裂隙密度，可得到公式(9)中的拟合参数。

步骤s5：考虑到岩石中的裂缝(软孔隙)纵横比及含量在一定范围内呈谱分布，因此，裂缝随压力增大闭合时，每个有效压力p未闭合的所有裂缝的最小初始纵横比为

其中，kst(γp)是有效压力p下的岩石等效静态体积模量，可由公式(5)计算得到。

将公式(9)代入到公式(10)中：

对公式(11)进行γⁱ到γp(α)之间的积分，可得到

其中，是高有效压力下的等效静态泊松比。

结合公式(9)和(12)，可以得到孔隙纵横比和压力p之间的关系式：

其中，是高有效压力下的等效静态杨氏模量，定义为泊松比由计算得到。

考虑到有效压力p逐渐增大时，具有不同纵横特征的裂缝随压差增量dp的变化量是相同的。因此，当压力的变化量即dp足够小时，裂隙密度的减少量可归因于纵横比小于的裂缝闭合造成的。故可以基于实验室测量的静态体积模量获得岩石内部裂缝孔隙度的分布特征。包括以下过程：

s51：首先将实验室测量的压力区间分成n等份，利用孔隙纵横比与压力之间的关系式，得到各压力下对应的纵横比；

s52：其次通过s4得到压力与累计裂隙密度关系式，求取各压力下对应的累积裂隙密度；

s53：最后通过裂缝孔隙度的关系式建立裂缝孔隙度与纵横比之间的关系，从而获得岩石内部裂缝孔隙度的分布特征。

步骤s6：基于砂岩实验室测量的应力应变数据，计算得到岩石中累积裂缝孔隙度和累计裂隙密度曲线，则二者曲线上的渐进值即分别为待测岩石内部的裂缝孔隙度与裂隙密度。

实施例

以navajo和weber砂岩为实施例来进一步说明，使用的应力应变与超声波速度实验数据，来自于coyner(1984)实验室测量的navajo和weber砂岩数据。实验测量条件为：干燥状态、有效压力可达100mpa。navajo砂岩的孔隙度为0.118，干燥岩石的颗粒密度为2316kg/m³，基质颗粒的体积与剪切模量分别为30gpa与33gpa。weber砂岩的孔隙度为0.095，干燥岩石的颗粒密度为2392kg/m³，基质颗粒的体积模量与剪切模量分别为42gpa与30gpa。

图2为navajo和weber砂岩的应力-应变曲线图，利用zimmerman(1991)提出的应力拟合关系式进行拟合。针对navajo砂岩，可得到δv/v＝0.0489p+0.0007(1-e^-p/0.014)，其中压力p的单位为gpa，拟合优度r²＝0.9998。而对于weber砂岩，则拟合曲线关系式为δv/v＝0.0563p+0.0059(1-e^-p/0.0111)，其拟合优度r²＝0.9999。实验数据显示：由于裂缝的存在，两种砂岩的应力-应变曲线都展现出非线性特征。

图3与4分别为navajo和weber砂岩纵横波速度随有效压力的变化曲线。可发现有效压力较小时，随着裂缝的闭合，速度变化剧烈，该趋势与众多学者的研究发现一致(如carcione等，2003；agersborg等，2008；han等，2011；asefet等人，2013；yurikov等，2017)。正如步骤s1与s2所诉，利用高有效压力下的实验室测量数据估算硬孔隙的纵横比。静态弹性模量所估算的navajo和weber砂岩的硬孔隙纵横比分别为0.18与0.08，而动态弹性模量所估算硬孔隙纵横比分别为0.28与0.11。

按照步骤s3与s4可得到压力与裂隙密度的拟合关系。对于navajo砂岩，基于静态弹性模量得到的拟合关系式为γp,st＝0.8842e^-p/0.0089，拟合优度为r²＝0.9509；而基于动态弹性模量得到的关系式为γp,dy＝0.2007e^-p/0.0237，拟合优度为r²＝0.9894。针对weber砂岩，基于静态模量得到的关系式为γp,st＝5.1001e^-p/0.0089，拟合优度为r²＝0.9887；基于动态模量得到的关系式为γp,dy＝1.4596e^-p/0.0154，拟合优度为r²＝0.9951。

按照步骤s5可得到navajo和weber砂岩裂缝孔隙度的分布特征。采用现有技术，基于实验室测量的超声波速度数据，基于动态弹性模量来计算岩石中孔隙纵横比的分布特征、累计裂隙密度和裂缝孔隙度，并与步骤s5的计算结果进行对比分析。如图5所示，图中实线与双划线分别表示基于静态与动态弹性模量的计算结果。如图可知，基于静态体积模量估算出的裂缝孔隙度比动态弹性模量的预测高。基于静态弹性模量估算得到navajo砂岩中裂缝孔隙度的峰值为0.00035％，对应的纵横比为0.00023；weber砂岩中分别为0.0023％和0.00026。同样的，基于动态弹性模量估算得到的navajo砂岩裂缝孔隙度的峰值和纵横比分别是7×10^-5％和0.0006，weber砂岩分别为6×10^-4％和0.00045。

按照步骤s6可得到navajo和weber砂岩的累积裂缝孔隙度和累计裂隙密度分布特征。根据dz模型，当累积裂缝孔隙度或累计裂隙密度达到渐近值时，二者的渐进值分别为岩石内部的裂缝孔隙度与裂隙密度。采用现有技术，基于实验室测量的超声波速度数据，基于动态弹性模量来计算岩石中累积裂缝孔隙度和累计裂隙密度，并与步骤s6的计算结果进行对比分析，如图6和7所示。如图所示，基于静态体积模量估算的navajo砂岩的累积裂缝孔隙度为0.084％，累积裂隙密度为0.88，而基于动态体积模量的裂缝孔隙度为0.047％，裂隙密度为0.2。对于weber砂岩，基于静态体积模量估算的累积裂缝孔隙度为0.55％，累积裂隙密度为5.1，而基于动态体积模量的裂缝孔隙度为0.25％，裂隙密度为1.5。结果表明，基于静态体积模量估算的累积裂缝的孔隙度和密度皆高于动态体积模量得到的值。

为说明基于静态弹性模量测试得到的岩石孔隙结构数据更符合岩石内部真实情况，将测试所得的岩石孔隙度随压力变化的曲线与实验测量值进行对比，并与基于动态弹性模量测试得到的相应的数据进行比较，如图8所示。图中硬孔隙的孔隙度是基于高有效应力下数据的线性拟合关系估算得到(pervukhina等，2010)。对于navajo砂岩，基于静态弹性模量数据所估算总孔隙度的相关系数为0.9933，而动态弹性模量所估算的相关系数0.9849。另外，对于weber砂岩，相关系数分别为0.9820和0.9038。结果表明，基于静态体积模量估算的裂缝孔隙度与实验数据更吻合。