0 引言

岩石的弹性性质受孔隙度、孔隙结构、矿物成分、及流体性质等多种因素的综合影响，其中孔隙结构是影响岩石特性的重要因素.由于岩石孔隙空间分布十分复杂，孔隙结构难以定量表征，因此需要借助理想的几何形状对复杂的孔隙结构进行近似处理.目前，最常用的孔隙结构表征方法是将岩石中的孔隙包含物近似为椭球体，并采用椭球体的短长轴之比(即孔隙纵横比)作为描述孔隙结构特征的主要参数.为了研究椭球状孔隙包含物的弹性性质，Eshelby(1957)基于弹性理论和假想实验建立了椭球包含物弹性质及其几何形状与岩石基质弹性模量之间的内在联系，并给出了椭球包含物的弹性场表达式.在Eshelby理论的基础上，Wu(1966)与Dunn和Ledbetter(1995)进一步推导了椭球孔隙压缩系数P和剪切柔度Q的解析表达式，其中P和Q又称为孔隙几何因数(Mavko et al., 2009)，它们均为孔隙包含物弹性模量、孔隙纵横比、以及其所在背景基质弹性模量的函数，是描述椭球孔隙结构特征的重要参数.这一系列工作为基于椭球孔隙包含物假设的各种岩石物理模型的建立和发展奠定了重要理论基础.

等效介质理论是联系岩石微观孔隙结构与岩石等效弹性性质的重要桥梁.目前，运用较为广泛的等效介质理论包括：Kuster-Toksöz(KT)模型(Kuster and Toksöz，1974)、Mori-Tanaka(MT)模型(Mori and Tanaka, 1973；Benveniste，1987；Ferrari，1994)、自相容近似(Self-Consistent Approximation, SCA)模型(Budiansky，1965；Berryman，1980;Berryman and Wang, 1995；Wu，1966)和微分等效介质(Differential effective medium, DEM)模型(Zimmerman，1984；Norris，1985；Berryman et al., 2002；李宏兵和张佳佳，2014).KT模型由Kuster和Toksöz(1974)在长波长假设条件下基于散射理论推导得到，具有计算岩石等效弹性性质的显式表达式.与KT模型类似，MT模型也具有显式的解析理论公式，因此KT模型和MT模型的数值实现十分简单，且更便于实际应用.但已有研究表明，KT模型和MT模型仅适用于包含物稀疏分布的情形(Kuster and Toksöz，1974；Berryman and Bergr, 1996；Mavko et al., 2009)，并且在某些条件下，KT模型和MT模型的等效弹性模量可能会超出Hashin-Shtrikman界限(Berryman，1980；Norris，1989；Ferrari，1991).Berryman和Bergr(1996)指出，只有当岩石基质占岩石总体积的比例高于70~80%时，KT模型和MT模型才能得到比较合理的结果.因此，在运用KT和MT模型时需十分谨慎.与显式的KT模型和MT模型不同，SCA模型和DEM模型属于隐式耦合模型，这些模型均涉及一系列迭代或微分运算.尽管SCA模型和DEM模型的计算相对复杂，但它们考虑了包含物之间的弹性相互作用，从而可以将孔隙包含物推广到含量稍高的情况(Mavko et al., 2009).在上述等效介质模型中，岩石等效弹性模量均为孔隙纵横比和孔隙度的非线性复杂函数.由此可见，孔隙结构参数的获取对等效介质理论的有效运用十分重要.

一直以来，许多学者尝试采用各种方式从实际数据中提取孔隙纵横比，根据这些方法的反演结果，大致可以分为两类：等效纵横比反演法和纵横比分布反演法.等效纵横比反演法通常假设岩石仅含一组或两组孔隙，并利用实际数据(如岩芯数据、测井数据等)反演计算各组孔隙的等效纵横比，比如：Kumar和Han (2005) 假设岩石孔隙空间仅由两种孔隙构成，即背景孔和球形孔或背景孔和微裂隙，利用Wyllie时间平均方程和DEM模型，由含水饱和碳酸盐岩的实测速度反演计算两种孔隙的等效孔隙纵横比及体积含量；之后，Xu和Payne(2009)在此基础上提出了适用于碳酸盐岩的Xu-Payne模型；Zhao等(2013)进一步将该建模思路应用于测井和地震数据的孔隙纵横比反演；Fortin等(2007)与Adelinet等(2011)假设岩石由球形硬孔和单一纵横比的微裂隙构成，在各个压力下分别对砂岩和玄武岩进行了孔隙纵横比反演；De Paula等(2012)提出了一种由干燥岩石纵横波速度-压力数据获取中孔(即纵横比介于0.01和0.2之间的孔隙)和软孔(即纵横比小于0.01的孔隙)等效孔隙纵横比的反演方法；郭继亮等(2017)采用单重孔隙岩石物理模型对碳酸盐岩孔隙形态参数(或孔隙纵横比)和孔隙度进行了反演；李宏兵等(2019)基于单一孔隙纵横比的DEM模型建立了复杂孔隙储层的三维岩石物理量版.

尽管等效纵横比法实现简单且容易应用，但球状硬孔隙的假设通常难以满足实际情形；此外，微裂隙对压力的变化十分敏感，假设所有微裂隙仅具有一个“等效纵横比”而不考虑压力对孔隙形态的影响往往是不够的.与等效纵横比反演不同，纵横比分布反演法认为岩石的孔隙空间是由一系列纵横比不同的孔隙构成的，因而更符合实际条件.早在20世纪60年代，Walsh(1965)便指出，岩石纵横波速度的压力依赖性是由岩石中孔隙的闭合引起的；不久，Morlier(1971)和Toksöz等(1976)根据速度与压力之间的变化关系，提出了孔隙纵横比分布的反演思路；之后，Cheng和Toksöz(1979)、Zimmerman(1991)与David和Zimmerman(2012)进一步将孔隙纵横比分布反演法与等效介质理论结合起来.在这些孔隙纵横比反演方法中，David和Zimmerman(2012)的孔隙结构模型(简称D-Z模型)最为经典，该模型假设岩石由固体矿物基质、一组纵横比相等的硬孔隙以及多组纵横比不等的微裂隙构成，并认为固体矿物基质和硬孔隙均不受压力影响，然后以累积裂隙密度为桥梁，借助等效介质理论将岩石弹性模量和孔隙纵横比联系起来；并在此基础上，利用超声纵横波速度的压力依赖性反演岩石硬孔隙和各组微裂隙的孔隙纵横比及孔隙度.尽管D-Z模型能够获得岩石的完整孔隙纵横比分布，但在该模型中，多重孔隙岩石的累积裂隙密度是由单重孔隙条件下的裂隙密度公式近似计算得到的，而这种近似处理导致D-Z模型在许多情况下无法获得良好的反演精度.此外，D-Z模型仅考虑了DEM和MT理论，而在许多情况下，KT和SCA等其他等效介质理论的应用也十分广泛.

为此，本文提出一种基于虚拟降压的孔隙纵横比分布反演方法，将基于DEM和MT的经典D-Z模型推广到KT和SCA中，结合四种等效介质理论建立了一套完整的反演流程，该流程通过模拟假想降压过程实现累积裂隙密度的准确反演，从而有效克服经典D-Z模型的精度不足.本文首先简要回顾了四种等效介质模型的基本理论；然后，根据压力对孔隙形态的影响规律，推导了微裂隙孔隙纵横比和裂隙密度的计算公式；在此基础上，给出了基于虚拟降压的孔隙纵横比分布反演流程；最后，应用实际岩芯数据对该反演方法进行了测试和分析.

1 等效介质理论的基本公式

等效介质理论是联系微观孔隙结构与岩石弹性性质的重要工具.本文假设岩石孔隙形状均为扁圆椭球体状(见图 1，图中α表示孔隙纵横比，a和c分别为椭球孔的长轴和短轴)，并利用KT模型、MT模型、DEM模型和SCA模型四种不同的等效介质理论实现硬孔隙和微裂隙孔隙纵横比的反演计算.这里我们首先回顾四种等效介质理论的基本公式.

1.1 KT模型

基于散射理论，Kuster和Toksöz(1974)推导了多重孔隙岩石的Kuster-Toksoz表达式，其一种推广形式可以表示为

(1)

以及

(2)

式中，(K_KT^*, G_KT^*) 为岩石的KT等效体积模量和剪切模量；(K_m, G_m) 为岩石矿物基质的体积模量和剪切模量，通常由Voight-Reuss-Hill平均公式得到(Mavko et al., 2009)；(K_i, G_i) 为第i种孔隙包含物的体积模量和剪切模量，其中i=1, 2, …, N；在干燥条件下，岩石孔隙不充填任何流体，此时K_i=G_i=0；v_i表示第i种孔隙包含物在总孔隙度ϕ中所占的体积百分比；(P^mi, Q^mi) 为几何因数，其上标“mi”表示以岩石矿物基质m作为背景介质并在其中加入第i种孔隙包含物；对于椭球状孔隙，P和Q是背景介质弹性模量、孔隙包含物弹性模量及其孔隙纵横比的复杂函数，其具体表达式参见附录A.

对于单重孔隙的干燥岩石，即N=1且K_i=G_i=0，式(1)可以简化为

(3)

由式(1)和式(3)可见，KT模型的理论公式是显式且非耦合的，因此由上述公式可以直接计算出岩石的KT等效弹性模量(K_KT^*, G_KT^*).

1.2 MT模型

根据Mori-Tanaka理论(Mori and Tanaka, 1973)，多重孔隙岩石的等效弹性模量满足如下方程(Berryman and Berge, 1996)：

(4)

式中，(K_MT^*, G_MT^*) 为岩石的MT等效体积模量和剪切模量；(K_i, G_i) 为第i相材料的体积模量和剪切模量；f_i表示第i相材料在岩石中所占的体积百分比；几何因数(P^mi, Q^mi) 的上标“mi”表示以岩石的主相材料作为背景介质并在其中加入第i种材料.

如果岩石仅由体积百分比分别为1-ϕ和ϕ的矿物基质和单重孔隙构成，则式(4)可以退化为

(5)

式中，(K_m, G_m) 为岩石矿物基质的体积模量和剪切模量；(K_i, G_i)为孔隙包含物的体积模量和剪切模量.上述过程运用了关系P^mm=Q^mm=1，此式表明当孔隙包含物材料与主相材料相同时，几何因数为1.在干燥条件下，式(5)还可以进一步简化为

(6)

从式(4)-(6)可以看出，MT公式中的(P^mi, Q^mi) 与等效模量(K_MT^*, G_MT^*) 无关，因此，与KT模型类似，MT模型公式也是显式且非耦合的.

1.3 DEM模型

在多重孔隙情形下，Norris形式的DEM公式可以写成(Norris, 1985; 李宏兵和张佳佳, 2014)

(7)

上述微分方程满足初始条件：

(8)

式中，(K_DEM^*, G_DEM^*) 为岩石的DEM等效体积模量和剪切模量；几何因数(P^*i, Q^*i) 的上标“*i”表示以DEM等效材料作为背景介质并在其中加入第i种孔隙包含物.当N=1时，式(7)退化为经典的DEM模型(Mavko et al., 2009).

在干燥条件下，上述耦合微分方程具有解析近似式，即(李宏兵和张佳佳, 2014)

(9)

式中，a, b, S_0i, S_1i, S_2i和S_3i的具体表达式详见附录B.

1.4 SCA模型

Berryman(1980; Berryman and Wang, 1995)给出了N相混合物的自相容近似(SCA或CPA)的一般形式：

(10)

式中，(K_SCA^*, G_SCA^*) 为岩石的SCA等效弹性模量；(K_i, G_i) 为第i相材料的弹性模量，其中i=1, 2, …, N；f_i表示第i相材料在岩石中所占的体积百分比.与前面3种等效介质模型不同，SCA模型可以由各相材料直接构造等效介质，而无需给定主相材料.

对于由矿物基质和单重孔隙构成的两相混合物，Wu(1966)给出了如下自相容模量的计算公式：

(11)

在干燥条件下，上式变为

(12)

这里需要注意的是，尽管式(10)和式(11)均为SCA模型公式，但两者来源却完全不同.通常，如果一个模型既能由准静态分析得到，也能通过散射近似获得，则称该模型为“物理可实现”，否则称“物理不可实现”.在上述公式中，Berryman提出的SCA模型为“物理可实现”模型，而Wu提出的SCA模型为“物理不可实现”模型.事实上，KT模型和MT模型均为“物理不可实现”模型，而DEM模型为“物理可实现”模型(Berryman and Bergr, 1996).虽然“物理不可实现”模型无法适用于所有类型的孔隙形状，但对于许多形状(如球形)，它们仍能给出十分准确的估计结果.

2 微裂隙孔隙纵横比及裂隙密度

已有研究表明，在不考虑孔隙弹性相互作用的情况下，孔隙的闭合压力与孔隙纵横比近似成正比，即p≈αE₀，其中E₀为岩石基质的杨氏模量(David and Zimmerman, 2012; 邓继新等, 2015).这意味着，若岩石为典型砂岩(E₀~50 GPa)，则需要施加500 MPa的压力才能使纵横比为0.01的孔隙完全闭合，而此闭合压力远高于实验测量压力.由此可见，在实验条件下，α≥0.01的孔隙几乎不受压力影响.为此，我们将最高实验压力下无法闭合的孔隙统称为硬孔.在本文中，我们假设岩石由一组纵横比为α_s的硬孔和多组纵横比不同的硬币形微裂隙构成，其中微裂隙的总孔隙度为ϕ_c，硬孔的孔隙度为ϕ_s=ϕ-ϕ_c，ϕ为岩石的总孔隙度.由于硬孔几乎不受实验压力的影响，故纵横波速度的压力依赖性完全可以归因于微裂隙的闭合.一般而言，纵横比越小的微裂隙越容易被压缩(Shapiro, 2003)，所以随着有效压力的增大，岩石中越“扁”的微裂隙越先闭合；当压力达到某一临界值P_c后，岩石内部的所有微裂隙全部闭合，此时岩石的纵横波速度几乎不再变化，且干燥岩石可以视为由等效矿物基质和单重硬孔构成的两相混合物(见图 2d)，而其弹性模量(K_D^h, G_D^h)可以利用如下经验关系拟合得到(Shapiro, 2003)：

(13)

式中，(K_D, G_D) 表示干燥岩石的弹性模量，由实测超声速度V_{P, D}(P)，V_{S, D}(P)和密度ρ_D计算得到；(K_Dⁱ, G_Dⁱ) 为零压力下干燥岩石的初始模量；(K_D^h, G_D^h) 为干燥岩石的高压模量；为拟合参数.如果实验压力足够高并能使岩石内部的所有微裂隙全部闭合，则可直接由V_{P, D}^h=V_{P, D}(P_max)和V_{P, D}^h=V_{S, D}(P_max) 获得高压速度，进而计算出高压模量(K_D^h, G_D^h).

尽管硬孔隙不引入任何压力效应，但硬孔隙的存在仍会影响岩石的弹性性质.因此，在后续微裂隙孔隙纵横比的反演计算中，我们将硬孔隙和矿物基质构成的两相混合物作为背景介质(即弹性模量为(K_D^h, G_D^h) 的材料)，通过在其中逐渐加入微裂隙以引入压力对岩石弹性性质的影响.

2.1 孔隙纵横比的理论公式

下面我们从微裂隙随有效压力的变化规律出发，推导微裂隙孔隙纵横比的计算公式.考虑压力P下某一个纵横比为α的干燥硬币形微裂隙(即近似于孔隙纵横比很小的扁椭球孔隙)，由压缩定律可知，其相对体积变化为

(14)

式中，C_c(P, α) 表示压力P下纵横比为α的微裂隙的压缩系数；V_c(P) 表示压力P下微裂隙的体积，dV_c(P) 为相应的体积变化量.在微裂隙的闭合过程中，通常假设其长轴a保持不变，即dV_c/V_c=dα/α. 因此，式(14)变为

(15)

当α很小(即α→0)时，有(Zimmerman, 1991)

(16)

即αC_c(P, α)与α无关.式中，K和ν分别为体积模量和泊松比.

将式(15)从0到P进行积分，得到

(17)

式中，α⁰表示零压力下的微裂隙初始纵横比；第二项表示压力由0变化至P的过程中微裂隙纵横比的变化量.此外，硬币形微裂隙的纵横比通常较小(< 0.01)，因此第二项几乎是与α无关的(参见式(16))，这便意味着，在相同的压力变化区间内，所有微裂隙的纵横比变化量是相同的.

假设该微裂隙的闭合压力p等于P，即α(P)=0，则式(17)变为

(18)

式中，α⁰(p)表示闭合压力为p的微裂隙的初始孔隙纵横比，在后文中，闭合压力均用p表示，以与P相区分.上式表明，当压力由0增至该微裂隙的闭合压力p时，其纵横比变化量刚好等于其初始纵横比大小，即当施加的有效压力P等于微裂隙的闭合压力p时，微裂隙刚好完全闭合.注意，利用式(18)为计算单一微裂隙的孔隙纵横比公式，需要已知该微裂隙的压缩系数C_c，而此参数通常难以测量，因此，我们还需要进一步将式(18)与实验可测量的物理量(如岩石压缩系数)联系起来.

设岩石的总孔隙空间Ω_p由N组纵横比为α_i的微裂隙组成，其中i=1, 2, …, N，取C_p(P)为压力P下总孔隙空间的压缩系数，于是有

(19)

联立式(14)和式(19)，得到

(20)

式中，ϕ_c=Ω_p/Ω为微裂隙的总孔隙度，ϕ_i=V_c(α_i)/Ω为第i组微裂隙的孔隙度，其中Ω表示岩石体积.由于岩石中充满了各种各样孔隙纵横比的微裂隙，因此，式(20)中的求和符号可以转换为积分符号，且积分区间应包含压力P下保持开孔的所有微裂隙的孔隙纵横比.为此，不妨令

(21)

则式(20)变为

(22)

式中，α表示保持张开状态微裂隙的最小纵横比；c(α) 称为孔隙度分布函数，根据孔隙度与裂隙密度Γ之间的关系，即ϕ=4παΓ/3，c(α) 可以类似地表示为

(23)

式中，称为裂隙密度分布函数，其中表示纵横比大于α的所有微裂隙的累积裂隙密度(David and Zimmerman, 2012).联立式(22)-(23)，得到

(24)

当微裂隙纵横比较小时，αC_c(P, α)与α几乎无关，可以提到积分符号之外.于是，有

(25)

由此可见，累积裂隙密度表示的是压力P下所有开孔微裂隙的裂隙密度总和.

利用式(24)，式(18)可以重新写成

(26)

此外，ϕ_cC_p=1/K_D-1/K_D^h，于是有

(27)

显然，相较于式(18)，上式中的K_D和K_D^h均可以与实验数据联系起来.因此，一旦获得累积裂隙密度，我们便能利用式(27)反演获得各组微裂隙的初始纵横比α⁰(p₁), α⁰(p₂), …, α⁰(p_N)，其中p₁, p₂, …, p_N表示各组微裂隙对应的闭合压力.

2.2 裂隙密度的理论公式

由式(27)可见，微裂隙初始纵横比α⁰的计算主要涉及3个参数：干燥岩石的体积模量K_D、高压体积模量K_D^h以及累积裂隙密度，其中K_D^h可以直接由超声数据和拟合关系式(13)得到，而体积模量K_D和累积裂隙密度的计算则需借助等效介质理论.为了建立岩石弹性模量和累积裂隙密度之间的关系，这里主要考虑4种常用的等效介质模型，即KT模型、MT模型、DEM模型和SCA模型.

假设岩石仅包含一组纵横比为α=c/a的硬币形微裂隙，其中a和c分别为微裂隙的长轴和短轴，则该组微裂隙的裂隙密度Γ可以表示为

(28)

式中，Ω为岩石体积；L为岩石中微裂隙的数量.于是，微裂隙的总孔隙度ϕ_c可以写成

(29)

对于纵横比较小的硬币形微裂隙，其几何因数P和Q具有显式解析表达式(David and Zimmerman, 2011)：

(30)

(31)

式中，ν为岩石基质的泊松比，而岩石基质的选取与等效介质理论有关，比如，对于KT模型和MT模型，ν表示背景基质的泊松比，而对于DEM模型和SCA模型，ν则表示DEM或SCA等效介质的泊松比.另外，由式(29)-(31)可以发现，ϕP和ϕQ与α无关，且均为裂隙密度Γ的函数.由此可见，将ϕP和ϕQ引入等效介质模型中便可以将岩石弹性模量与裂隙密度Γ联系起来.

对于KT模型，我们可以将ϕP和ϕQ代入式(3)中得到裂隙密度的计算公式：

(32)

式中，K_b，G_b和ν_b分别表示背景基质的体积模量、剪切模量和泊松比.注意，这里的背景基质可以是等效矿物基质m，也可以是其它等效介质，其选择视具体问题而定.

对于MT模型，将ϕP和ϕQ代入式(6)中得到

(33)

考虑到微裂隙孔隙度ϕ_c≪1，我们在上述公式中运用了近似关系1-ϕ_c≈1.

对于DEM模型，David和Zimmerman(2011)推导了如下裂隙密度的解析近似式：

(34)

此近似式的精度略依赖于泊松比v_b，但在裂缝密度很高的情况下，其计算误差仍小于2%(David and Zimmerman, 2011).

对于SCA模型，将ϕP和ϕQ代入式(12)中有

(35)

这里需要注意的是，上述等效介质理论的裂隙密度公式均是在假设岩石仅含一组微裂隙的情况下推导得到的，因此理论上仅适用于单重孔隙岩石.然而，在前人的诸多研究中，这些公式被直接用来近似计算多重孔隙岩石的累积裂隙密度(Zimmerman, 1991; David and Zimmerman, 2012; 邓继新等, 2015).尽管这种做法可以很方便地获得累积裂隙密度，但在许多情况下该方法并不能很好地拟合实验数据，见“模型算例”部分.为此，我们提出一种基于虚拟降压的反演策略来实现微裂隙累积裂隙密度和孔隙纵横比的计算.

3 基于虚拟降压的反演策略 3.1 微裂隙的压力响应描述

我们将实际干燥岩石视为矿物基质(K_m, G_m)、硬孔ϕ_s和微裂隙ϕ_c的多相混合物，其中硬孔和微裂隙的方向均是随机的(如图 2a).由于岩石矿物基质与硬孔几乎不受实验压力的影响，因此可将两者的混合物(如图 2d)作为新的等效基质处理.如此，硬孔的贡献全部包含在新的等效基质中，而微裂隙所在的背景介质变为由硬孔和矿物基质构成的新材料，其弹性模量为(K_D^h, G_D^h). 在这种情况下，岩石仅由弹性模量为(K_D^h, G_D^h) 的等效基质与总孔隙度为ϕ_c的N组微裂隙构成(如图 3a所示).下面我们将在此基础上运用假想降压过程反演微裂隙的累积裂隙密度和孔隙纵横比分布.

设岩石中所有的微裂隙在压力P_c下刚好全部闭合，且压力区间[0, P_c] 由相等的N个子区间和N+1个压力点构成，其中压力点包括0, ΔP, 2ΔP, …, NΔP，此处ΔP=P_c/N表示压力增量.假设岩石中N组微裂隙的闭合压力刚好为p₁=ΔP, p₂=2ΔP, …, p_N=NΔP，且ΔP足够小时，各压力子区间内的孔隙纵横比可以认为是恒定的.当有效压力P_e∈[p_i, p_i+1) 时，岩石中初始孔隙纵横比小于α⁰(p_i+1) 的所有微裂隙全部闭合，而α⁰>α⁰(p_i) 的微裂隙仍处于张开状态.由式(16)-(17)可知，当压力增加ΔP时，所有微裂隙的纵横比变化量是相同的，这意味着，在压力P_e下岩石中仍处于开孔状态的各组微裂隙的孔隙纵横比均应相对于其初始值减小α⁰(p_i).此时，各组开孔微裂隙的孔隙纵横比应为

(36)

图 3展示了不同压力下岩石中各组微裂隙的开孔情况：当实验压力为0时，所有微裂隙均处于初始状态；当实验压力增至p₁，初始纵横比小于或等于α⁰(p₁)的微裂隙闭合，而其它微裂隙的纵横比均减小α⁰(p₁)；当实验压力继续增至p₂，满足α⁰≤α⁰(p₂)的微裂隙全部闭合，而仍处于开孔状态微裂隙的纵横比相对于其初始值减小α⁰(p₂)；以此类推下去，当实验压力增至压力P_c=NΔP时，岩石中的所有微裂隙全部闭合.由此可见，一旦获得各组微裂隙的初始纵横比(即微裂隙的初始纵横比分布)，利用上述过程便能求出不同压力下岩石的微裂隙纵横比分布.

3.2 孔隙纵横比反演流程

为了获得微裂隙的初始纵横比，首先需要求解微裂隙的累积裂隙密度(见式(27))，这里我们借助假想降压过程进行计算.注意，这里的假想降压过程是指增压逆过程，即对应于图 3中由右往左、由下至上的过程，而并不等同于实际降压过程.事实上，吴春燕等(2016)的实验研究表明，在实际降压过程中岩石孔隙空间通常难以恢复到原始状态.

基于虚拟降压的反演流程主要包括三部分，即微裂隙的累积裂隙密度反演、微裂隙孔隙纵横比分布计算以及硬孔纵横比反演，工作流程见图 4.

(1) 累积裂隙密度反演

根据假想降压思路，我们将降压过程分为N个变化状态，即p_N→p_N-1, …, p₂→p₁, p₁→0，其中p_N等于最高实验压力P_c. 在每个变化状态，干燥岩石均可以视为单重孔隙结构.这意味着，我们可以直接运用第2.2节推导的单重孔隙岩石的等效介质理论公式(32)-(35)计算微裂隙的裂隙密度，具体过程如下：

① 当p_N→p_N-1时，我们分别将压力p_N和p_N-1下的弹性模量作为岩石背景基质和岩石等效介质的弹性参数，即用高压模量(K_D^h, G_D^h)和压力p_N-1下的(K_D, G_D) 替换式(32)-(35)中(K_b, G_b)和(K^*, G^*)；然后，利用非线性最小二乘反演算法，得到初始纵横比为α⁰(p_N) 的微裂隙的裂隙密度Γ(α_N⁰)，其中α_N⁰=α⁰(p_N)；

② 当p_i→p_i-1(i=N-1, …, 2, 1)时，分别用压力p_i和p_i-1下的(K_D, G_D) 替换式(32)-(35)中的基质模量(K_b, G_b)和等效模量(K^*, G^*)，并结合非线性最小二乘反演算法计算初始纵横比为α⁰(p_i) 的微裂隙的裂隙密度Γ(α_i⁰)；

③ 依此类推，直至获得N组微裂隙的裂隙密度Γ(α₁⁰), Γ(α₂⁰), …, Γ(α_N⁰)；

④ 利用各组微裂隙的裂隙密度计算累积裂隙密度：

(37)

式中，表示α⁰≥α_k⁰的所有微裂隙的裂隙密度之和，也可以写成闭合压力p_k的函数，即，表示压力p_k下所有开孔微裂隙的裂隙密度之和.

在上述过程中，我们自然而然地认为实验压力点是等间隔分布的.然而，在很多情况下，实际数据并不满足这些要求，此时我们可以利用关系式(13)对实验数据点进行拟合，然后利用扩展后的数据完成上述过程.

另外，我们还可以结合如下公式将累计裂隙密度扩展到更高的压力：

(38)

式中，为初始累积裂隙密度；为与压力量纲相同的参数.尽管上述公式并不具有任何物理含义，但是可以为累积裂隙密度提供很好的拟合关系.

(2) 孔隙纵横比分布函数计算

将已获得的累积裂隙密度代入式(36)中，便能计算出微裂隙的初始纵横比分布α⁰(p)，结合式(36)，还可以进一步给出任意非零有效压力下的孔隙纵横比分布α(p).除微裂隙的孔隙纵横比分布之外，根据如下关系还可以获得微裂隙的初始裂隙密度分布函数γ(α⁰)、初始孔隙度分布函数c(α⁰)、初始累积孔隙度分布函数C(α⁰)以及初始孔隙度ϕ_i(α⁰)：

(39)

(40)

(41)

(42)

利用上述初始分布函数，可以进一步计算任意非零有效压力下的分布函数γ(α)、c(α)、C(α)以及ϕ_i(α). 根据第4.1节知，当有效压力P_e∈[p_i, p_i+1) 时，岩石中仅存在初始纵横比大于α⁰(p_i) 的微裂隙，而α⁰≤α_i⁰=α⁰(p_i) 的微裂隙已全部闭合.由于闭合的微裂隙对裂隙密度没有任何贡献，因此有效压力P_e下的裂隙密度分布函数γ(α)应为

(43)

由此可见，对于α⁰＞α⁰(p_i) 的微裂隙，其压力P_e下的裂隙密度分布函数与初始分布函数γ(α⁰)相同，其原因在于：裂隙密度只与长轴a有关(见式(28))，而a不随压力变化，这意味着在微裂隙没有完全闭合的情况下，其裂隙密度不随外部压力而改变，而只有当微裂隙完全闭合后，该微裂隙对岩石裂隙密度的贡献方降为零.

一般而言，压力对岩石体积Ω的影响可以忽略不计，因此由微裂隙孔隙度的定义ϕ_i=4πa³α/(3Ω) 可知，孔隙度ϕ_i与纵横比α的比值等于常量.于是，在有效压力P_e下，ϕ_i(α) 可由初始孔隙纵横比α⁰和初始孔隙度分布函数ϕ_i(α⁰)表示为

(44)

此外，联立式(41)、式(42)和式(44)还可以求出有效压力P_e下的孔隙度分布函数c(α) 以及累积孔隙度分布函数C(α).

(3) 硬孔纵横比反演

硬孔纵横比反演是基于高压状态下的实验数据实现的.在高压条件下，干燥岩石可以视为由等效矿物基质和单重硬孔构成的两相混合物(见图 2d)，此时利用等效介质理论便可以从高压数据中提取出硬孔纵横比.图 5给出了硬孔纵横比的反演流程，具体步骤如下：

① 获取岩芯样品的基本信息，包括样品的总孔隙度ϕ、矿物组分及其含量、干燥岩石的密度ρ_D以及各有效压力下的超声纵横波速度V_{P, D}(P)和V_{S, D}(P)，并根据反演得到的微裂隙的总孔隙度ϕ_c，计算硬孔的孔隙度ϕ_s=ϕ-ϕ_c. 对于ϕ_c≪ϕ_s的情况，也可以由ϕ_s≈ϕ来估算硬孔的孔隙度；

② 利用Voigt-Reuss-Hill平均公式(简称V-R-H平均)计算等效矿物基质模量(K_m, G_m)，即(Mavko et al., 2009)

(45)

式中，f_i为各矿物组分的体积含量；M_i表示各矿物组分的体积模量或剪切模量.

③ 基于单重孔隙干燥岩石的等效介质理论公式(式(3)、式(6)、式(9)或式(12))，结合高压速度(V_{P, D}^h, V_{S, D}^h) 反演计算硬孔纵横比α_s.本文采用非线性最小二乘方法在孔隙纵横比区间[0.01, 1]内寻找最优的α_s，使其同时满足纵波速度和横波速度的相对误差最小，即

(46)

式中，(V_{P, EMT}^h, V_{S, EMT}^h) 表示由等效介质理论模拟的高压速度；ε表示可容许的最小相对误差.

4 实际应用分析

为了分析上述反演流程的可靠性和适用性，我们分别选取了具有代表性的2块砂岩样品S1和S2与2块致密碳酸盐岩样品C1和C2，并采用基于虚拟降压的反演方法对其孔隙结构参数进行了提取.表 1列出了这些岩芯样品的基本参数和数据来源，其中样品S1为纯石英砂岩，孔隙度为4%，岩芯数据来自于文献David (2012)；样品S2为泥质砂岩，孔隙度为23.9%，主要矿物成分为石英(66%)、长石(17%)、粘土(15%)以及极少量的黄铁矿和菱铁矿，其孔隙结构特征见图 6(a-b)；样品C1和C2为纯净的致密碳酸盐岩，孔隙度分别为3%和0.6%，其矿物基质均由100%的方解石构成，两块碳酸盐岩样品均为致密灰岩，其中样品C1主要发育晶间孔和粒内孔，而样品C2具有较强的非均质特征，孔隙以局部发育的平直微裂缝为主，孔隙结构特征见图 6(c-f).

图 7给出了不同实验压力下干燥砂岩样品S1，S2与碳酸盐岩样品C1，C2的纵横波速度超声实验测量结果.由图 7可见，砂岩样品和碳酸盐岩样品对压力的弹性响应存在显著差异.当有效压力由0 MPa增至20 MPa时，砂岩样品的纵横波速度迅速增大；但随着有效压力的继续增大，其纵横波速度的增长趋势逐渐减缓，达到最大值后趋于平稳.这表明砂岩样品中存在大量小纵横比的微裂隙(如粒间孔)以及少量纵横比较大的微裂隙，由于小纵横比微裂隙的闭合压力较小且数量极多，所以当压力略微增加时大量的微裂隙闭合，从而使岩石的纵横波速度迅速增加.从图 7a-b可以推测出，砂岩样品中的大部分微裂隙在有效压力到达20 MPa时便已全部闭合，而仅剩余少数纵横比较大的微孔处于开孔状态；当压力继续增大并超过剩余微孔的闭合压力后，这些纵横比较大的孔隙才相继闭合，但由于其数量较少，岩石速度随压力的增加趋势也逐渐减缓.对于致密碳酸盐岩样品而言，纵横波速度随压力的变化则相对平缓(见图 7c-d).这意味着，与小纵横比微裂隙占主导的砂岩样品相比，碳酸盐岩样品中各种微裂隙的数量差异相对较小.因此，在有效压力较低的范围内，其纵横波速度随压力的增加趋势也相对平缓.

此外，图 7还给出了由经验关系式(13)得到的实验数据拟合曲线(图中虚线)，其中黑色箭头指示的是岩芯样品的高压速度(V_{P, D}^c, V_{S, D}^c)，该值由拟合参数(K_D^c, G_D^c) 计算得到.泥质砂岩样品S2的等效矿物基质模量通过对其所有矿物组分进行V-R-H平均(式(45))计算得到，模量大小分别为K_m=36.5 GPa和G_m=28.1 GPa.将岩石的高压速度和矿物基质模量代入单重孔隙干燥岩石的等效介质理论公式(式(3)、式(6)、式(9)或式(12))中，我们便能得到由不同等效介质理论反演的硬孔纵横比.表 2列出了基于KT、MT、DEM和SCA模型计算的硬孔纵横比α_s及其误差ε，其中ε定义为纵波速度相对误差和横波速度相对误差的平均值，即

(47)

由表 2可见，对于纯净砂岩样品S1和非均质程度较低的孔洞型碳酸盐岩样品C1，KT、MT、DEM和SCA模型均能取得很好的反演效果，其反演误差均小于0.5%，且这些模型估计的硬孔隙纵横比十分接近(样品S1：~0.4，样品C1：~0.18).对于泥质砂岩样品S2，四种模型的反演精度略低，其中DEM模型的效果最佳，其反演误差为2.05%，KT模型的效果最差，其反演误差为4.31%.由此可见，泥质及其他掺杂矿物对砂岩硬孔纵横比的反演精度影响较大，其原因在于：岩石弹性模量的主要贡献来源于矿物成分，而对于包含有多种矿物组分的非纯净砂岩，V-R-H平均仅能提供一个粗略的等效矿物基质模量估计值，因此存在一定的误差.另外，对于裂缝型碳酸盐岩样品C2，四种等效介质理论的反演结果均不理想，最大反演误差接近5%.由于C2为纯净碳酸盐岩样品，矿物组分的影响几乎可以忽略，因此，该误差主要来源于碳酸盐岩的复杂孔隙结构分布.由图 6f可见，样品C2存在局部发育的裂缝，具有较强的非均质性，而本文采用的反演算法理论上更适合非均质性较弱的岩石.

利用第3.2节的孔隙纵横比反演流程，我们可以进一步提取岩芯样品的微裂隙纵横比分布.图 8展示了4块样品在有效压力0 MPa、10 MPa、30 MPa和50 MPa下的累积裂隙密度分布. 从图 8中可以看出，累积裂隙密度随孔隙纵横比的增大而减小，这是因为表示的是孔隙纵横比大于α的所有开孔微裂隙的裂隙密度之和，因此α越小意味着岩石中纵横比大于α的微裂隙越多，故累积裂隙密度越高.此外，由图可见，随着有效压力的增大，累积裂隙密度逐渐降低，曲线向低孔隙纵横比方向移动，且覆盖的孔隙纵横比范围变窄，其原因在于：压力增大后，岩石中纵横比较小的微裂隙先闭合，而纵横比较大的微裂隙逐渐被压缩，从而导致开孔微裂隙的孔隙纵横比减小且分布范围变窄.从图中还可以看到，在低有效压力下，由不同等效介质理论计算的 趋势基本一致，但数值上存在略微差异，且这种差异随着有效压力的增大逐渐消失.总体而言，等效介质理论的选择对孔隙纵横比的反演结果影响较小.

图 9给出了不同有效压力下岩石微裂隙的孔隙度曲线ϕ_i(α).由图 9可见，孔隙度ϕ_i(α)随孔隙纵横比和有效压力的变化规律特征与累积裂隙密度函数 类似，即随着有效压力的增大，微裂隙的孔隙纵横比分布区间变窄，曲线向低纵横比移动，且四种模型的反演结果逐渐趋于相同.从图中可以看出，在砂岩样品S1中，纵横比为0.2×10^-3的微裂隙占主导，其孔隙度为0.85×10^-3%，微裂隙的孔隙纵横比主要集中在区间0~0.7×10^-3；在砂岩样品S2中，纵横比为0.65×10^-3的微裂隙占主导，其孔隙度约为8×10^-3%，微裂隙的孔隙纵横比主要集中在区间0~2×10^-3；当有效压力为50 MPa时，两块砂岩样品中的微裂隙几乎完全闭合，微裂隙孔隙度趋近于零.由图 9c-d可见，在有效压力50 MPa下，碳酸盐岩样品C2中的微裂隙也基本完全闭合，而样品C1中则还存在少量的开孔孔隙.这个现象与岩芯实际情况是完全相符的，因为样品C1为孔洞型碳酸盐岩，而样品C2为裂缝型碳酸盐岩，裂缝的孔隙纵横比更小更容易闭合，故样品C2比样品C1更容易被压缩.此外，相较于裂缝型碳酸盐岩样品C2，孔洞型碳酸盐岩样品C1具有更宽的纵横比分布范围(0~1.7×10^-3)，并且其孔隙度曲线与横坐标轴之间并不是完全闭合的，这意味着该岩芯中还存在更大纵横比(＞1.7×10^-3)的微裂隙.

下面我们从数据统计的角度分析四种等效介质理论的反演效果.我们以31块干燥泥质砂岩和40块干燥碳酸盐岩的超声纵横波速度作为输入数据，其中砂岩样品的主要矿物成分为石英、长石、方解石、白云石和粘土，孔隙度分布范围为3.5%~30%；碳酸盐岩样品主要由方解石构成，其孔隙度分布范围为0.3%~3%.图 10给出了四种理论的模拟结果及超声实验数据.从图中可以看出，无论是泥质砂岩还是碳酸盐岩，四种等效介质理论的模拟速度与实测速度均具有很好的一致性.结合图 8和图 9可见，这些理论的反演精度和反演结果均十分接近.由此可见，本文提出的孔隙结构反演方法并不十分依赖于等效介质理论的选择，四种等效介质理论均能作为反演孔隙结构参数的有效工具.但是，在计算效率方面，由于KT和MT模型具有显式的理论表达式，其计算速度比隐式的DEM和SCA模型更快.

为了对比本文方法与现有方法，我们采用基于虚拟降压的反演方法(见第4.2节)和D-Z方法(Zimmerman, 1999; David and Zimmerman, 2012; 邓继新等, 2015)计算了不同有效压力下干燥岩石的弹性模量，以此评估两种方法的反演精度.本文方法与D-Z方法的主要区别在于，D-Z方法直接采用单重孔隙岩石的裂隙密度公式(式(33)和式(34))近似反演多重孔隙岩石的累积裂隙密度；而本文方法则借助假想降压过程实现多重孔隙岩石的累积裂隙密度计算.以DEM模型为例，图 11和图 12分别从单块岩芯与统计分析的角度对比了两种方法的反演结果.由图 11可见，对于纯净砂岩样品S1，本文方法和D-Z方法均能取得很好的反演效果；但对于泥质砂岩和碳酸盐岩，由D-Z方法计算的体积模量与实测数据存在较大误差，且计算误差随着有效压力的减小而增大，这表明D-Z模型所采用的近似做法(即直接利用单重孔隙公式(33)-(34)近似计算多重孔隙岩石的累积裂隙密度)存在一定局限性.另一方面，从图 11-12中可以看出，本文方法相较于D-Z方法具有更高的精度，无论是砂岩还是碳酸盐岩，由本文方法计算得到的体积模量和剪切模量与实际数据吻合更好.由此可见，本文方法比D-Z方法更具优势，能够有效克服经典D-Z反演方法的精度不足.

5 结论

本文对经典D-Z方法进行了改进和完善，并将其推广至KT、MT、DEM和SCA四种等效介质模型的情形，建立了一套完整的孔隙纵横比反演流程.以一系列砂岩样品和碳酸盐岩样品为例，分析了本文方法在实际岩芯孔隙纵横比提取中的应用效果，数值分析结果表明：基于本文方法从干燥岩芯速度-压力数据中提取的孔隙纵横比分布、累积裂隙密度分布和微裂隙孔隙度分布可以很好地解释实际岩芯中微裂隙的压力响应规律；此外，本文方法并不十分依赖于等效介质理论的选择，即四种理论均能取得十分好的反演精度，且反演的孔隙结构参数差异很小；与经典D-Z方法相比，本文方法有效改善了D-Z方法中利用单重孔隙公式近似计算多重孔隙岩石累积裂隙密度所引起的精度不足.在实际应用方面，本文所建立的反演流程可以为岩石物理实验和理论研究提供重要的孔隙结构参数提取工具；另外，结合深度和压力的对应关系，此方法也有望应用于测井数据的岩石孔隙结构参数反演中.

附录A 椭球状孔隙的几何因数

对于纵横比小于1的椭球状孔隙，其几何因数P和Q可以表示为(Mavko et al., 2009)

式中，

其中：

附录B DEM解析近似式的系数

在DEM的解析近似式中，a, b, S_0i, S_1i, S_2i和S_3i的具体表达式分别为

式中，