地球物理学报  2021, Vol. 64 Issue (3): 1016-1033   PDF    
基于等效介质理论的孔隙纵横比分布反演
欧阳芳, 赵建国, 李智, 肖增佳, 贺艳晓, 赵皓, 任静     
中国石油大学(北京)油气资源与探测国家重点实验室, 北京 102249
摘要:孔隙纵横比是描述多孔岩石微观孔隙结构特征的重要参数,目前用于获取岩石完整孔隙纵横比分布的经典模型为David-Zimmerman (D-Z)孔隙结构模型,该模型假设岩石由固体矿物基质、一组纵横比相等的硬孔隙以及多组纵横比不等的微裂隙构成,并认为固体矿物基质和硬孔隙均不受压力影响,在此基础上,利用超声纵横波速度的压力依赖性反演岩石硬孔隙和各组微裂隙的孔隙纵横比及孔隙度.该方法的关键点在于以累积裂隙密度为桥梁,借助等效介质理论建立了岩石弹性模量和孔隙纵横比之间的内在联系.但在D-Z模型中,多重孔隙岩石累积裂隙密度的计算直接由单重孔隙裂隙密度公式实现,这种近似导致该模型在许多情况下难以获得良好的反演精度.为了完善经典D-Z模型,本文提出了一种基于虚拟降压的孔隙纵横比分布反演策略,通过多个假想降压过程实现累积裂隙密度的准确计算,并将基于DEM和MT的经典D-Z模型推广到KT和SCA中,结合四种等效介质理论建立了一套完整的反演流程.采用一系列砂岩和碳酸盐岩样品,测试了反演流程在实际岩芯孔隙纵横比提取中的应用效果,研究结果表明:与D-Z模型相比,本文方法的模拟结果与实际数据吻合更好,并同时适用于砂岩和碳酸盐岩;此外,通过分析四种等效介质理论的模拟结果发现,本文方法并不十分依赖于等效介质理论的选择,这些理论获得的孔隙结构参数随压力的变化趋势基本一致,数值上仅存在略微差异,且这种差异随着压力的增大逐渐消失.本文方法是经典D-Z孔隙结构模型的重要补充,对岩石孔隙结构表征、流体饱和岩石速度预测以及孔间喷射流效应的模拟具有十分重要的意义.
关键词: 岩石孔隙结构      孔隙纵横比      虚拟降压      等效介质理论     
Inversion of pore aspect ratio distribution based on effective medium theories
OUYANG Fang, ZHAO JianGuo, LI Zhi, XIAO ZengJia, HE YanXiao, ZHAO Hao, REN Jing     
State Key Laboratory of Petroleum Resource and Prospecting, China University of Petroleum(Beijing), Beijing 102249, China
Abstract: The pore aspect ratio is an important parameter to describe pore structure characteristics of porous rocks. At present,the most classical model used to obtain complete pore aspect ratio distribution of rocks is the David-Zimmerman (D-Z) pore structure model. The model assumes that the rock is composed of solid mineral matrix,a group of stiff pores having the same aspect ratio and a distribution of micro-cracks with different aspect ratios,in which the material matrix and stiff pores are pressure-independent. Based on the assumptions,the pore aspect ratio and porosity of stiff pores and micro-cracks are extracted from the pressure dependence of ultrasonic velocities. The key point of this method is that the cumulative crack density is taken as a bridge to establish the intrinsic relationship between rock elastic modulus and pore aspect ratio through effective medium theory. However,in the D-Z model,the cumulative fracture density of multi-porosity rock is directly calculated using the formula of crack density theoretically only applicable for rock with a single pore type,which makes it difficult to achieve good inversion accuracy in many cases. In order to improve the classical D-Z model,this paper proposes an inversion method for calculating pore aspect ratio distribution using thought experiment. In this method,the cumulative crack density is accurately computed through a fictitious unloading path,and the classical D-Z model based on DEM and MT is extended to the case of KT and SCA. As a result,a detailed inversion workflow combining with four effective medium theories is established. This workflow is applied in the aspect ratio distribution extraction for a set of sandstones and carbonate samples to illustrate the capability of the method. The results show that the proposed method is in better agreement with the actual data as compared to the D-Z model,and can be applied to both sandstone and carbonate samples. In addition,the proposed method is almost independent of the selection of the effective medium theories. The pore structure parameter obtained from the four theories changes with pressure in a similar way. Although there is a slight difference in these results,it gradually disappears with increasing pressure. Our method is an important supplement to the classical D-Z pore structure model. It is of great significance in the characterization of rock pore structure,prediction of fluid saturated rock velocity and simulation of pore-scale squirt flow.
Keywords: Pore structure model of rock    Pore aspect ratio    Fictitious unloading method    Effective medium theory    
0 引言

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

等效介质理论是联系岩石微观孔隙结构与岩石等效弹性性质的重要桥梁.目前,运用较为广泛的等效介质理论包括:Kuster-Toksöz(KT)模型(Kuster and Toksöz,1974)、Mori-Tanaka(MT)模型(Mori and Tanaka, 1973Benveniste,1987Ferrari,1994)、自相容近似(Self-Consistent Approximation, SCA)模型(Budiansky,1965Berryman,1980;Berryman and Wang, 1995Wu,1966)和微分等效介质(Differential effective medium, DEM)模型(Zimmerman,1984Norris,1985Berryman et al., 2002李宏兵和张佳佳,2014).KT模型由Kuster和Toksöz(1974)在长波长假设条件下基于散射理论推导得到,具有计算岩石等效弹性性质的显式表达式.与KT模型类似,MT模型也具有显式的解析理论公式,因此KT模型和MT模型的数值实现十分简单,且更便于实际应用.但已有研究表明,KT模型和MT模型仅适用于包含物稀疏分布的情形(Kuster and Toksöz,1974Berryman and Bergr, 1996Mavko et al., 2009),并且在某些条件下,KT模型和MT模型的等效弹性模量可能会超出Hashin-Shtrikman界限(Berryman,1980Norris,1989Ferrari,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,图中α表示孔隙纵横比,ac分别为椭球孔的长轴和短轴),并利用KT模型、MT模型、DEM模型和SCA模型四种不同的等效介质理论实现硬孔隙和微裂隙孔隙纵横比的反演计算.这里我们首先回顾四种等效介质理论的基本公式.

图 1 扁圆椭球状孔隙 Fig. 1 The oblate-ellipsoid pore
1.1 KT模型

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

(1)

以及

(2)

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

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

(3)

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

1.2 MT模型

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

(4)

式中,(KMT*, GMT*) 为岩石的MT等效体积模量和剪切模量;(Ki, Gi) 为第i相材料的体积模量和剪切模量;fi表示第i相材料在岩石中所占的体积百分比;几何因数(Pmi, Qmi) 的上标“mi”表示以岩石的主相材料作为背景介质并在其中加入第i种材料.

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

(5)

式中,(Km, Gm) 为岩石矿物基质的体积模量和剪切模量;(Ki, Gi)为孔隙包含物的体积模量和剪切模量.上述过程运用了关系Pmm=Qmm=1,此式表明当孔隙包含物材料与主相材料相同时,几何因数为1.在干燥条件下,式(5)还可以进一步简化为

(6)

从式(4)-(6)可以看出,MT公式中的(Pmi, Qmi) 与等效模量(KMT*, GMT*) 无关,因此,与KT模型类似,MT模型公式也是显式且非耦合的.

1.3 DEM模型

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

(7)

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

(8)

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

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

(9)

式中,a, b, S0i, S1i, S2iS3i的具体表达式详见附录B.

1.4 SCA模型

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

(10)

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

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

(11)

在干燥条件下,上式变为

(12)

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

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

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

(13)

图 2 受压后岩石内部硬孔隙和微裂隙的形态变化示意图 Fig. 2 Schematic diagram of stiff pores and micro-cracks in rock at different hydrostatic pressures

式中,(KD, GD) 表示干燥岩石的弹性模量,由实测超声速度VP, D(P),VS, D(P)和密度ρD计算得到;(KDi, GDi) 为零压力下干燥岩石的初始模量;(KDh, GDh) 为干燥岩石的高压模量;为拟合参数.如果实验压力足够高并能使岩石内部的所有微裂隙全部闭合,则可直接由VP, Dh=VP, D(Pmax)和VP, Dh=VS, D(Pmax) 获得高压速度,进而计算出高压模量(KDh, GDh).

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

2.1 孔隙纵横比的理论公式

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

(14)

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

(15)

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

(16)

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

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

(17)

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

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

(18)

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

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

(19)

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

(20)

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

(21)

则式(20)变为

(22)

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

(23)

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

(24)

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

(25)

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

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

(26)

此外,ϕcCp=1/KD-1/KDh,于是有

(27)

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

2.2 裂隙密度的理论公式

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

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

(28)

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

(29)

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

(30)

(31)

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

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

(32)

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

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

(33)

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

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

(34)

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

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

(35)

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

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

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

图 3 增压过程中岩石内各组微裂隙的形态变化示意图 Fig. 3 The change of each group of micro-cracks in rock during pressure increasing

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

(36)

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

3.2 孔隙纵横比反演流程

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

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

图 4 基于虚拟降压的孔隙纵横比反演流程图 Fig. 4 Inversion workflow of aspect ratio distribution based on a thought unloading method

(1) 累积裂隙密度反演

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

① 当pNpN-1时,我们分别将压力pNpN-1下的弹性模量作为岩石背景基质和岩石等效介质的弹性参数,即用高压模量(KDh, GDh)和压力pN-1下的(KD, GD) 替换式(32)-(35)中(Kb, Gb)和(K*, G*);然后,利用非线性最小二乘反演算法,得到初始纵横比为α0(pN) 的微裂隙的裂隙密度Γ(αN0),其中αN0=α0(pN);

② 当pipi-1(i=N-1, …, 2, 1)时,分别用压力pipi-1下的(KD, GD) 替换式(32)-(35)中的基质模量(Kb, Gb)和等效模量(K*, G*),并结合非线性最小二乘反演算法计算初始纵横比为α0(pi) 的微裂隙的裂隙密度Γ(αi0);

③ 依此类推,直至获得N组微裂隙的裂隙密度Γ(α10), Γ(α20), …, Γ(αN0);

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

(37)

式中,表示α0αk0的所有微裂隙的裂隙密度之和,也可以写成闭合压力pk的函数,即,表示压力pk下所有开孔微裂隙的裂隙密度之和.

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

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

(38)

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

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

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

(39)

(40)

(41)

(42)

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

(43)

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

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

(44)

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

(3) 硬孔纵横比反演

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

图 5 硬孔纵横比反演流程图 Fig. 5 Inversion workflow for aspect ratio of the stiff pores

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

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

(45)

式中,fi为各矿物组分的体积含量;Mi表示各矿物组分的体积模量或剪切模量.

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

(46)

式中,(VP, EMTh, VS, EMTh) 表示由等效介质理论模拟的高压速度;ε表示可容许的最小相对误差.

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).

表 1 砂岩和碳酸盐岩样品的基本参数 Table 1 Physical properties of the sandstone and carbonate samples
图 6 砂岩样品(a-b)和碳酸盐岩样品(c-f)的孔隙结构特征 Fig. 6 Pore structure characteristics of the sandstone (a-b) and carbonate (c-f) samples

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

图 7 干燥岩芯样品的纵横波速度-压力曲线 (a) 纯净砂岩样品S1; (b) 泥质砂岩样品S2;(c) 孔洞型致密碳酸盐岩样品C1;(d) 裂缝型致密碳酸盐岩样品C2. Fig. 7 Velocity versus confining pressure for dry rock samples (a) Pure sandstone sample S1; (b) Shaley sandstone sample S2; (c) Tight carbonate sample C1 with inter-granular pores; (d) Tight carbonate sample C2 with cracks.

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

(47)

表 2 基于不同等效介质理论反演的硬孔隙纵横比 Table 2 Inversion results of the characteristic aspect ratio for the stiff pore obtained from different effective medium theories

表 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中可以看出,累积裂隙密度随孔隙纵横比的增大而减小,这是因为表示的是孔隙纵横比大于α的所有开孔微裂隙的裂隙密度之和,因此α越小意味着岩石中纵横比大于α的微裂隙越多,故累积裂隙密度越高.此外,由图可见,随着有效压力的增大,累积裂隙密度逐渐降低,曲线向低孔隙纵横比方向移动,且覆盖的孔隙纵横比范围变窄,其原因在于:压力增大后,岩石中纵横比较小的微裂隙先闭合,而纵横比较大的微裂隙逐渐被压缩,从而导致开孔微裂隙的孔隙纵横比减小且分布范围变窄.从图中还可以看到,在低有效压力下,由不同等效介质理论计算的 趋势基本一致,但数值上存在略微差异,且这种差异随着有效压力的增大逐渐消失.总体而言,等效介质理论的选择对孔隙纵横比的反演结果影响较小.

图 8 基于4种等效介质理论反演的不同有效压力下砂岩(a-b)和碳酸盐岩(c-d)样品的累积裂隙密度分布 Fig. 8 Cumulative crack density distribution as a function of the aspect ratio calculated by different effective medium theories at different pressures for sandstone (a-b) and carbonate samples (c-d)

图 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)的微裂隙.

图 9 基于4种等效介质理论反演的不同有效压力下砂岩(a-b)和碳酸盐岩(c-d)样品的微裂隙孔隙度 Fig. 9 Crack porosity as a function of the aspect ratio calculated by different effective medium theories at different pressures for sandstone (a-b) and carbonate samples (c-d)

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

图 10 基于不同等效介质理论的模型反演结果与实际测量结果对比 Fig. 10 Comparison of the measured velocities and the results inverted using different effective medium theories

为了对比本文方法与现有方法,我们采用基于虚拟降压的反演方法(见第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反演方法的精度不足.

图 11 本文方法和D-Z方法基于DEM理论反演的干燥岩石弹性模量对比 (a) 纯净砂岩样品S1; (b) 泥质砂岩样品S2;(c) 孔洞型致密碳酸盐岩样品C1;(d) 裂缝型致密碳酸盐岩样品C2. Fig. 11 Elastic Moduli of dry rock calculated with DEM by using our method and D-Z method (a) Pure sandstone sample S1; (b) Shaley sandstone sample S2; (c) Tight carbonate sample C1 with inter-granular pores; (d) Tight carbonate sample C2 with cracks.
图 12 本文方法和D-Z模型反演结果对比 Fig. 12 Comparison of our method and D-Z method
5 结论

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

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

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

式中,

其中:

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

在DEM的解析近似式中,a, b, S0i, S1i, S2iS3i的具体表达式分别为

式中,

References
Adelinet M, Dorbath C, Le Ravalec M, et al. 2011. Deriving microstructure and fluid state within the Icelandic crust from the inversion of tomography data. Geophysical Research Letters, 38(3): L03305. DOI:10.1029/2010GL046304
Benveniste Y. 1987. A new approach to the application of Mori-Tanaka's theory in composite materials. Mechanics of Materials, 6(2): 147-157. DOI:10.1016/0167-6636(87)90005-6
Berryman J G. 1980. Long-wavelength propagation in composite elastic media I. Spherical inclusions. The Journal of the Acoustical Society of America, 68(6): 1809-1819. DOI:10.1121/1.385171
Berryman J G, Wang H F. 1995. The elastic coefficients of double-porosity models for fluid transport in jointed rock. Journal of Geophysical Research: Solid Earth, 100(B12): 24611-24627. DOI:10.1029/95JB02161
Berryman J G, Berge P A. 1996. Critique of two explicit schemes for estimating elastic properties of multiphase composites. Mechanics of Materials, 22(2): 149-164. DOI:10.1016/0167-6636(95)00035-6
Berryman J G, Pride S R, Wang H F. 2002. A differential scheme for elastic properties of rocks with dry or saturated cracks. Geophysical Journal International, 151(2): 597-611. DOI:10.1046/j.1365-246X.2002.01801.x
Budiansky B. 1965. On the elastic moduli of some heterogeneous materials. Journal of the Mechanics and Physics of Solids, 13(4): 223-227. DOI:10.1016/0022-5096(65)90011-6
Cheng C H, Toksöz M N. 1979. Inversion of seismic velocities for the pore aspect ratio spectrum of a rock. Journal of Geophysical Research: Solid Earth, 84(B13): 7533-7543. DOI:10.1029/JB084iB13p07533
David E C, Zimmerman R W. 2011. Elastic moduli of solids containing spheroidal pores. International Journal of Engineering Science, 49(7): 544-560. DOI:10.1016/j.ijengsci.2011.02.001
David E C, Zimmerman R W. 2012. Pore structure model for elastic wave velocities in fluid-saturated sandstones. Journal of Geophysical Research: Solid Earth, 117(B7): B07210. DOI:10.1029/2012JB009195
David E C. 2012. The effect of stress, pore fluid and pore structure on elastic wave velocities in sandstones[Ph. D. thesis]. London: Imperial College London.
De Paula O B, Pervukhina M, Makarynska D, et al. 2012. Modeling squirt dispersion and attenuation in fluid-saturated rocks using pressure dependency of dry ultrasonic velocities. Geophysics, 77(3): WA157-WA168. DOI:10.1190/geo2011-0253.1
Deng J X, Zhou H, Wang H, et al. 2015. The influence of pore structure in reservoir sandstone on dispersion properties of elastic waves. Chinese Journal of Geophysics (in Chinese), 58(9): 3389-3400. DOI:10.6038/cjg20150931
Dunn M L, Ledbetter H. 1995. Poisson's ratio of porous and microcracked solids: Theory and application to oxide superconductors. Journal of Materials Research, 10(11): 2715-2722. DOI:10.1557/JMR.1995.2715
Eshelby J D. 1957. The determination of the elastic field of an ellipsoidal inclusion and related problems. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 241(1226): 376-396. DOI:10.1098/rspa.1957.0133
Ferrari M. 1991. Asymmetry and the high concentration limit of the Mori-Tanaka effective medium theory. Mechanics of Materials, 11(3): 251-256. DOI:10.1016/0167-6636(91)90006-L
Ferrari M. 1994. Composite homogenization via the equivalent poly-inclusion approach. Composites Engineering, 4(1): 37-45. DOI:10.1016/0961-9526(94)90005-1
Fortin J, Guéguen Y, Schubnel A. 2007. Effects of pore collapse and grain crushing on ultrasonic velocities and Vp/Vs. Journal of Geophysical Research: Solid Earth, 112(B8): B08207. DOI:10.1029/2005JB004005
Guo J L, Li H B, Zhang Y, et al. 2017. Inverting carbonate reservoir porosity affected by pore aspect ratio. Progress in Geophysics (in Chinese), 32(1): 146-151. DOI:10.6038/pg20170120
Kumar M, Han D H. 2005. Pore shape effect on elastic properties of carbonate rocks. //75th Ann. Internat Mtg., Soc. Expi. Geophys. . Expanded Abstracts, 1477-1480, doi: 10.1190/1.2147969.
Kuster G T, Toksöz M N. 1974. Velocity and attenuation of seismic waves in two-phase media: Part I. Theoretical formulations. Geophysics, 39(5): 587-606. DOI:10.1190/1.1440450
Li H B, Zhang J J. 2014. A differential effective medium model of multiple-porosity rock and its analytical approximations for dry rock. Chinese Journal of Geophysics (in Chinese), 57(10): 3422-3430. DOI:10.6038/cjg20141028
Li H B, Zhang J J, Cai S J, et al. 2019. 3D rock physics template for reservoirs with complex pore structure. Chinese Journal of Geophysics (in Chinese), 62(7): 2711-2723. DOI:10.6038/cjg2019K0672
Mavko G, Mukerji T, Dvorkin J. 2009. The Rock Physics Handbook: Tools for Seismic Analysis of Porous Media. 2nd ed. New York: Cambridge University Press.
Mori T, Tanaka K. 1973. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metallurgica, 21(5): 571-574. DOI:10.1016/0001-6160(73)90064-3
Morlier P. 1971. Description de l'état de fissuration d'une roche à partir d'essais non-destructifs simples. Rock Mechanics, 3(3): 125-138. DOI:10.1007/BF01238439
Norris A N. 1985. A differential scheme for the effective moduli of composites. Mechanics of Materials, 4(1): 1-16. DOI:10.1016/0167-6636(85)90002-x
Norris A N. 1989. An examination of the Mori-Tanaka effective medium approximation for multiphase composites. Journal of Applied Mechanics, 56(1): 83-88. DOI:10.1115/1.3176070
Shapiro S A. 2003. Elastic piezosensitivity of porous and fractured rocks. Geophysics, 68(2): 482-486. DOI:10.1190/1.1567215
Toksöz M N, Cheng C H, Timur A. 1976. Velocities of seismic waves in porous rocks. Geophysics, 41(4): 621-645. DOI:10.1190/1.1440639
Walsh J B. 1965. The effect of cracks on the compressibility of rock. Journal of Geophysical Research, 70(2): 381-389. DOI:10.1029/JZ070i002p00381
Wu C Y, Zhang J G, Sun Z G, et al. 2016. An experimental study on sandstone reservoir parameters under variable pressures in Dongying sag. Geological Science and Technology Information (in Chinese), 35(1): 101-106.
Wu T T. 1966. The effect of inclusion shape on the elastic moduli of a two-phase material. International Journal of Solids and Structures, 2(1): 1-8. DOI:10.1016/0020-7683(66)90002-3
Xu S Y, Payne M A. 2009. Modeling elastic properties in carbonate rocks. The Leading Edge, 28(1): 66-74. DOI:10.1190/1.3064148
Zhao L X, Nasser M, Han D H. 2013. Quantitative geophysical pore-type characterization and its geological implication in carbonate reservoirs. Geophysical Prospecting, 61(4): 827-841. DOI:10.1111/1365-2478.12043
Zimmerman R W. 1984. Elastic moduli of a solid with spherical pores: New self-consistent method. International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, 21(6): 339-343. DOI:10.1016/0148-9062(84)90366-8
Zimmerman R W. 1991. Compressibility of Sandstones. New York: Elsevier.
邓继新, 周浩, 王欢, 等. 2015. 基于储层砂岩微观孔隙结构特征的弹性波频散响应分析. 地球物理学报, 58(9): 3389-3400. DOI:10.6038/cjg20150931
郭继亮, 李宏兵, 张研, 等. 2017. 受孔隙形态影响的碳酸盐岩孔隙度反演. 地球物理学进展, 32(1): 146-151. DOI:10.6038/pg20170120
李宏兵, 张佳佳. 2014. 多重孔岩石微分等效介质模型及其干燥情形下的解析近似式. 地球物理学报, 57(10): 3422-3430. DOI:10.6038/cjg20141028
李宏兵, 张佳佳, 蔡生娟, 等. 2019. 复杂孔隙储层三维岩石物理模版. 地球物理学报, 62(7): 2711-2723. DOI:10.6038/cjg2019K0672
吴春燕, 张金功, 孙志刚, 等. 2016. 东营凹陷增压减压条件下砂岩储层物性参数的变化规律实验研究. 地质科技情报, 35(1): 101-106.