用孔隙、裂隙介质弹性波理论反演岩石孔隙分布特征

唐晓明， 王鹤鸣， 苏远大， 陈雪莲*

1 中国石油大学(华东)地球科学与技术学院， 青岛 266580 2 青岛海洋科学与技术试点国家实验室， 青岛 266580 3 中国石油大学(华东)深层油气重点实验室， 青岛 266580

0 引言

岩石中的孔隙分布特征对岩石的力学、声学和流体渗透性质有着十分重要的影响，是储层岩石声学关注的重点.用纵横比这一参数来度量，岩石孔隙分布的显著特征是包含纵横比～1的孔隙和纵横比≪1的裂隙；不同形态的裂隙分布可以用纵横比谱来很好地描述.多年来，确定岩石孔隙的纵横比谱一直是岩石物理学研究的一个方向；实验压力加载条件下的弹性波速测量为此研究提供了一条有效的途径.许多学者利用不同纵横比的裂隙对压力的不同响应来反演岩石的孔隙纵横比谱.Cheng(1978)和Cheng和Toksöz(1979)运用Kuster和Toksöz(1974)的裂缝模型(K-T模型)，结合孔隙体积随压力的变化理论(Toksöz et al.，1976)，提出了从实验流体饱和岩石加压测量的纵、横波速度数据反演孔隙纵横比谱的方法.David和Zimmerman(2012)提出了基于干燥岩石数据的类似方法.近年来又有很多学者在这方面做了大量的工作(Yan et al.，2014，2015；邓继新等，2015；Duan et al.，2018；Han et al.，2019；李闯等，2020).

本文利用唐晓明(2011)和Tang等(2012)的孔隙、裂隙介质的弹性波动理论来反演岩石孔隙的纵横比谱.对比于早期反演采用的K-T低频近似模型，新近的理论较好地描述了岩石中孔隙与裂隙相互作用产生的衰减和频散，更适宜于实验室超声频段测量的声波数据.事实上，唐晓明等(2013)应用该理论成功地模拟和反演了实验室数据，得到了岩石裂隙密度和纵横比随压力的变化曲线.不足的是，该反演模拟的是单一纵横比裂隙体系在压力作用下的变化，不能得到反映岩石孔隙形态分布的纵横比谱.为此本文将该理论进一步扩展，使之包括孔隙与多形态裂隙体系的相互作用；这样一来，岩石中裂隙随压力的变化就可以描述为不同形态(纵横比)裂隙在不同压力下的形变和闭合，对应的弹性波速变化便体现了不同纵横比裂隙的贡献.由此可见，孔隙与多形态裂隙体系的相互作用是反演孔隙形态分布的理论基础.

用扩展后的孔隙、裂隙介质的弹性波动理论替代K-T模型，对几套经典岩石样品(包括高孔隙度砂岩、致密砂岩和花岗岩)实验数据进行反演，得到这些岩石的孔隙纵横比谱.与岩芯的扫描电镜观测和分析结果对比，新理论的反演结果比原有结果更加符合实际观测结果，有效地提高了从实验室岩芯超声数据获取孔隙分布特征的精度和可靠性.

1 孔隙、裂隙弹性波理论对多裂隙岩石的推广

1.1 理论回顾

唐晓明(2011)和Tang等(2012)的孔隙、裂隙并存的双孔介质弹性波理论给出饱和岩石的体积模量、剪切模量的表达式为

(1)

(2)

其中Kd为干燥岩石的体积模量，α=1-Kd/Ks，β=(α-φ)/Ks+φ/Kf，Ks和Kf分别为岩石基质和孔隙流体的体积模量，φ为岩石的孔隙度，S(ω)为描述孔隙与裂隙相互作用的挤喷流函数，包含了裂隙密度和裂隙纵横比这两个描述裂隙的重要参数，K0、μ0分别为S(ω)=0时的饱和岩石体积模量与剪切模量.

针对硬币型的裂隙，唐晓明(2011)推导出的S(ω)表达式为

(3)

上述硬币模型中孔隙与裂隙的流体交换在硬币的边缘，但力学上裂隙在此是闭合的，作为模型的改进，Tang等(2012)提出了钹状的孔、裂隙模型，将流体交换放到硬币模型的中部，由此推导出的S(ω)表达式为

(4)

孔、裂隙介质的弹性模量确定后，介质中的快纵波、慢纵波和横波的波数由下式计算(Tang et al.，2012)：

(5)

(6)

其中下标p和s分别代表纵波和横波，+和-分别代表快纵波和慢纵波，上述公式的符号表达式为

(7)

其中ρs和ρf分别为岩石固体基质和孔隙流体的密度，θ=iκ(ω)/(ηω)，其中κ(ω)为Johnson等(1987)推导出的动态渗透率，κ(ω)的具体表达式为

(8)

其中κ0为达西渗透率，τ为孔隙内流体的弯曲度.

由以上得到波数可进一步计算波的速度频散和衰减(Tang and Cheng，2004)：

v=ω/Re{k},

Q-1=2Im{k}/Re{k},

(9)

其中v与Q分别为速度与品质因子，Re{k}与Im{k}分别取k的实部与虚部.

1.2 对含多形态裂隙岩石的推广

现在考虑岩石中存在多种形态裂隙的情况，在Tang等(2012)的单一形态裂隙的孔裂隙介质中，孔隙流体压力本构方程为：

(10)

其中w=φ(U-u)，U和u分别为流体和固体的位移，q为单位体积岩石由裂隙挤喷到孔隙空间的流体体积.当M个形态(用纵横比衡量)不一的裂隙并存时，q为各个形态裂隙挤喷的贡献之和，(10)式的压力本构方程变为

(11)

其中下标m(m=1,2,…,M)代表第m种纵横比(形态)的裂隙，用Sm(ω)=φqm/p表示该种裂隙的挤喷函数，(11)式可化简为

(12)

其中第m种裂隙的Sm(ω)由(4)式给出，其具体表达式为

(13)

基于孔、裂隙的弹性模量公式(1)和(2)，多裂隙体系岩石的体积模量、剪切模量的计算公式变为

(14)

(15)

在孔裂隙理论计算公式(5)—(9)式中，将Knew与μnew代替原来的KTang与μTang，就可以对多形态裂隙体系进行模拟计算.

根据上述多形态裂隙岩石的推广理论，模拟了这种岩石介质中的波速频散和衰减特征.模拟中假定介质中分布有纵横比从0.00003至0.01的多组裂隙，对应的裂隙密度从0.045至0.012依次降低(见表1)；理论计算所需的其他岩石参数取值在表2中给出.作为对比，也计算了每一组裂隙单独存在时的结果.图1、图2和图3分别给出各种情况下的快纵波、横波和慢纵波的频散曲线(a)与衰减曲线(b)，包括了单一形态裂隙(彩线)与多组裂隙叠加(黑线)结果的对比；前者的刻度由右侧的纵坐标给出，后者的刻度由左侧的纵坐标给出.单一形态裂隙的频散和衰减曲线存在两个弛豫频率段，分别对应于(由其纵横比决定的)挤喷流和(由岩石孔隙决定的)Biot流动的弛豫频率，后一频率对所有裂隙组都是相同的；当多个形态的裂隙组共同作用时，波的频散与衰减是Biot流动与每个裂隙组产生的挤喷流叠加的结果.相对于单组裂隙模型，多组形态裂隙模拟的快纵波和横波的波速与衰减没有单裂隙那样明显的弛豫特征，在全频域的变化特征与实际数据的更加相符(Adam et al.，2006).因此，多形态裂隙模型常用来模拟实际岩石的实验室测量数据(Yan et al.，2014，2015；邓继新等，2015).

表1 孔裂隙介质的裂隙分布Table 1 Distribution of cracks in a cracked porous medium

表2 孔裂隙介质参数Table 2 Parameters of a cracked porous medium

图1 多组裂隙条件下孔裂隙介质快纵波频散曲线(a)与衰减曲线(b)Fig.1 Dispersion (a) and attenuation (b) of fast compressional wave in a cracked porous medium with multiple cracks

图2 多组裂隙条件下孔裂隙介质横波频散曲线(a)与衰减曲线(b)Fig.2 Dispersion (a) and attenuation (b) of shear wave in a cracked porous medium with multiple cracks

图3 多组裂隙条件下孔裂隙介质慢纵波频散曲线(a)与衰减曲线(b)Fig.3 Dispersion (a) and attenuation (b) of slow compressional wave in a cracked porous medium with multiple cracks

由以上模拟结果可见，单一形态裂隙挤喷流的弛豫频段随纵横比增加从低频向高频移动，当裂隙的纵横比增大至0.01，挤喷流特征频率与Biot流特征频率接近.因此可以把纵横比为0.01的孔隙作为裂隙与硬孔的分界，这与David和Zimmerman(2012)设定的软、硬孔隙边界一致；邓继新等(2015)也说明了当纵横比大于0.01时，令其闭合所需要的压力远大于实验室的正常压力加载范围，这为以下的反演裂隙参数搜索范围提供了依据.

2 压力对孔隙结构的影响

处于压力加载条件下的孔、裂隙岩石，其孔隙结构发生重大变化，岩石中裂隙变得扁平，裂隙空间随着压力的增大而减小，裂隙逐渐闭合；而压力的变化对硬孔的形态没有太大影响.本文采用的钹状孔裂隙模型(Tang et al.，2012)，是在硬币裂隙上加上球形孔隙，其力学形变主要由扁形裂隙决定，可以用Toksöz等(1976)的包体理论来描述.该理论给出了扁球形包体的体积变化率和压力的关系；对于纵横比为γ的扁球形裂隙，其体积变化率dc/c与有效压力Pe的关系为

(16)

其中K*为干岩石的静态有效模量，Ei(i=1,2,3,4)为纵横比γ与岩石基质模量K和μ的函数；对扁球形包体模型，Ei的表达式为

(17)

Cheng(1978)和Cheng和Toksöz(1979)建立了压力与孔隙结构的函数关系.对于岩石在压力加载过程中的第n个压力点Pn(n=0,1,…,N)，第m组形态裂隙的体积含量cn m(即裂隙孔隙度)可以由0有效压力下的裂隙孔隙度c0m表示：

(18)

其中γ0m为0有效压力下第m组形态裂隙的纵横比.当dc/c<-1时，裂隙闭合.

对于压力作用下的扁球形裂隙，其短轴的变化远大于长轴变化，在忽略后一变化的条件下，纵横比的变化率即为孔隙体积变化率：

(19)

因此，压力Pn下的裂隙纵横比γn m可以由0压下的γ0m表示：

(20)

结合(18)—(20)式，得到：

(21)

上式的意义是，任意压力下的裂隙纵横比谱(即孔隙度与纵横比的函数关系)可以由0压时的纵横比谱计算得到.

以上公式给出了压力对裂隙孔隙度和裂隙纵横比的影响，本文用裂隙密度作为主控参数，因此需要将裂隙孔隙度转化为裂隙密度.利用上述的扁球形包体的形变理论，并结合(21)式，裂隙密度与裂隙孔隙度的关系可表示为

(22)

从上式及其推导过程可以看出，压力变化时裂隙长度不变因而裂隙密度也保持不变(即:εn m=ε0m)，达到裂隙的闭合压力时裂隙的作用消失.通过压力对裂隙的闭合作用，对岩石加压可直接影响其弹性模量.在理论计算上，裂隙闭合(裂隙密度改变)对弹性模量的影响可以通过将(14)和(15)式中的体积模量Kd和剪切模量μ0作为裂隙密度和孔隙度的函数来实现：

Kd=Kd(ε,φ),μ0=μ0(ε,φ),

(23)

具体的函数形式由Biot相洽理论(Thomsen，1985)给出.对于纵横比不同的多组裂隙体系，计算上式所需的裂隙密度为各组裂隙的裂隙密度之和.对第n个压力点Pn，计算时应将已经达到(18)式的闭合条件(dc/c<-1)的裂隙剔除，这时的裂隙密度为

(24)

其中M′为尚未闭合的裂隙组数.

3 孔隙结构反演

压力加载下的岩石波速测量是岩石物理学研究的重要手段，已有大量的实验数据.我们采用实验室测量的不同压力下饱和岩石纵、横波速数据来反演岩石的裂隙密度与纵横比随压力的变化，用(22)式将裂隙密度转化为裂隙孔隙度，由孔隙度与纵横比之间的关系得到岩石孔隙纵横比谱.在反演中，将理论波速与实验波速数据的均方差作为目标函数，调节模型参数使目标函数达到最小值，目标函数E的具体表达式如下：

(25)

(26)

(26)式的目标函数只需反演0压下的裂隙纵横比值，再由此0压下的纵横比值及其随压力的变化关系计算各压力点下的纵横比值，此时需要反演的参数减少为2×M+3个.裂隙形态随压力变化关系的运用大大减少了反演未知参数的个数，同时也降低了反演的不确定性.

下面给出0压下裂隙纵横比及其谱元素量M的取值方法.一般来说，这种取值受到裂隙闭合、压力加载范围、测量压力点数N以及岩性的影响.为了达到数据拟合的最佳效果，采用Cheng(1978)建议的方法，令每个压力点下都有裂隙闭合.对于裂隙纵横比的取值，应满足如下条件：

(1)由于把纵横比为0.01的孔隙作为裂隙与硬孔的分界，在孔隙、裂隙介质的弹性波动理论中硬孔不易变形，且对硬孔纵横比没有限定，因此硬孔的纵横比可在0.01至1之间任取；

(2)裂隙纵横比的设置要考虑压力的影响，基于(20)式，可以设置0压下纵横比的取值范围，以保证每个压力点下都有裂隙闭合，且设置的裂隙纵横比谱的取值点应均匀分布；

(3)考虑岩性的影响，如致密岩石一般在高压段波速变化较小，此段的裂隙已接近闭合，故岩石中裂隙纵横比应该较小.

对于裂隙数量M的确定，应满足如下条件：

(1)考虑到每个压力点下都有裂隙闭合，M的数量通常应大于N；

(2)不同岩性的M有所不同，致密岩石(或高压段速度变化较小的岩石)的M个数较少，甚至可以小于N.

为了提高反演结果的准确性，将0压下的总孔隙度(即裂隙孔隙度和纵横比大于0.01的硬孔的孔隙度之和)作为约束条件：

(27)

在孔裂隙理论中，硬孔的孔隙度φp远大于裂隙的孔隙度，因此对于(27)式，可令0压下的硬孔孔隙度φp保持不变，反演时调节0压下的裂隙密度，使0压下的裂隙孔隙度与φp之和满足(27)式.

孔隙结构的反演可以看作是搜寻目标函数最小值的过程，为了降低反演结果的非唯一性，我们采取全局优化的遗传(GA)算法(Goldberg，1988)对(26)式进行求解.

4 实例分析

将上述反演方法应用于几组经典数据，分别是Kayenta砂岩，较为致密的Weber砂岩，以及Westerly花岗岩的实验室测量数据(Coyner，1984).上述岩样数据都有早期方法反演的裂隙纵横比谱.更为重要的是，这些岩样还有扫描电镜的分析结果，可以用来验证反演结果的正确性并比较本文结果和已有反演结果的异同.

4.1 Kayenta砂岩

首先讨论Kayenta砂岩岩芯测量数据(Coyner，1984)的反演结果.Kayenta砂岩为高孔隙度岩石，孔隙度为22.2%，骨架密度为2017 kg·m-3，渗透率取为63 mD(King and Marsden，2002)，声速测量的频率为1 MHz.这些实验数据提供了理论计算所需的岩石参数(参见表2)值.图4a给出了该砂岩样品在苯饱和条件下纵、横波速在1到70 MPa压力范围的测量数据(圆形标识符).反演得到的岩石基质的体积模量和剪切模量分别为31.00 GPa和24.04 GPa.图4a中的黑色实、虚线是用反演的岩石弹性模量和裂隙纵横比谱计算的理论饱和纵、横波速度，与实测数据吻合很好.图4b中红色部分是反演的0压下的裂隙纵横比谱，其分布在3.5×10-4至1×10-2的纵横比范围内.

图4 Kayenta砂岩反演结果(a) 压力-波速数据(标识符)及其理论拟合(曲线)； (b) 反演得到的孔隙纵横比谱； (c) 本文(红色)与Cheng方法(黑色)反演结果(标识符)与扫描电镜孔隙纵横比谱(直方图)的对比； (d) 硬孔(红色)、裂隙(绿色)孔隙度和裂隙密度(蓝色)随压力的变化Fig.4 Kayenta sandstone inversion results(a) Laboratory pressure-velocity data (markers) and theoretical fitting (curves) using inversion results; (b) Inverted pore aspect ratio spectrum; (c) The new (red) and existing (black) inversion results (markers) versus SEM pore aspect ratio spectrum (histogram); (d) Stiff (red) and crack (green) porosity and crack density (blue) variation with differential pressure

Burns等(1990)采用Cheng(1978)方法反演了Kayenta砂岩的孔隙纵横比谱，同时还给出了利用扫描电镜(SEM)技术观测得到的岩石孔隙纵横比谱，二者分别由图4c中的黑色数据点和直方图表示.图中直方的个数为(不同分辨率的)电镜扫描次数，直方的宽度表示某一分辨率下看到的裂隙纵横比的范围，其高度为该范围内裂隙的归一化孔隙度(用总孔隙度归一).将两种方法的反演结果放在直方图上，可以看到二者在反映孔隙的1到10-1的纵横比范围内基本是吻合的；但在反映裂隙分布的10-2到10-3的纵横比范围，新的反演结果(红圆圈)比已有的结果(黑圆圈)明显要低，且与扫描电镜结果更加吻合.

图4d给出了反演得到的硬孔孔隙度、裂隙密度以及由其换算(参见(22)式)得到的裂隙孔隙度随压力变化的数据.硬孔和裂隙的孔隙度均随压力的增大而减小，前者(后者)的变化十分微小(明显)；裂隙密度随压力的变化最为显著，随压力的增加迅速降低，对应着图4a的波速在低压段的快速增加，其原因是图4c分布谱上的低纵横比裂隙在压力作用下的闭合，这一物理机理使得我们能够从波速-压力数据反演裂隙的纵横比分布谱.

4.2 Weber砂岩

接下来考察Weber砂岩的岩芯测量数据(Coyner，1984)和结果.该砂岩较为致密，孔隙度为9.5%，骨架密度为2392 kg·m-3；理论计算所需的其他参数与Kayenta砂岩类似.图5a给出了该砂岩样品在苯饱和条件下纵、横波速在1到100 MPa压力范围的测量数据(圆形标识符).反演得到的岩石基质的体积模量和剪切模量分别为28.00 GPa和26.90 GPa，图5a中的黑色实、虚线是用反演的岩石弹性模量和裂隙纵横比谱计算的理论饱和纵、横波速度，与实测数据吻合很好.图5b中红色部分是反演的0压下的裂隙纵横比谱，其分布在2×10-4至6×10-3的纵横比范围内.

图5 Weber砂岩反演结果(a) 压力-波速数据(标识符)及其理论拟合(曲线)； (b) 反演得到的孔隙纵横比谱； (c) 本文(红色)与Cheng方法(黑色)反演结果(标识符)与扫描电镜孔隙纵横比谱(直方图)的对比； (d) 硬孔(红色)、裂隙(绿色)孔隙度和裂隙密度(蓝色)随压力的变化.Fig.5 Weber sandstone inversion results(a) Laboratory pressure-velocity data (markers) and theoretical fitting (curves) using inversion results; (b)Inverted pore aspect ratio spectrum; (c) The new (red) and existing (black) inversion results (markers) versus SEM pore aspect ratio spectrum (histogram); (d) Stiff (red) and crack (green) porosity and crack density (blue) variation with differential pressure.

图5c是Weber砂岩在1至10-3纵横比范围内扫描电镜的直方图(Burns et al.，1990)与两种反演结果(标识符)的对比.与Kayenta砂岩类似，新的反演(红色)与已有的反演(黑色)结果在1到10-1的纵横比范围内基本吻合；但在反映裂隙分布的10-2到10-3的纵横比范围，新结果与扫描电镜结果的吻合明显优于已有的结果.图5d的孔隙度和裂隙密度随压力变化的规律与图4d类似，不复赘述.对比于高孔的Kayenta砂岩，致密的Weber砂岩的孔隙度要小很多，但二者的裂隙发育度(裂隙密度)和形态分布(实测和反演的裂隙纵横比谱)都较为相似；因此，二者的波速随压力变化的形态和变化量较为一致.

4.3 Westerly花岗岩

图6给出了Westerly花岗岩的岩芯数据(Coyner，1984)和扫描电镜(Hadley，1976)及数据反演结果.该岩石的孔隙度极低，仅为0.8%，骨架密度为2641 kg·m-3.图6a给出了岩样在苯饱和条件下纵、横波速在2.5到100 MPa压力范围的测量数据(圆形标识符)，反演得到的岩石基质的体积模量和剪切模量分别为57.60和32.30 GPa.图6a中的黑色实、虚线是用反演的岩石弹性模量和裂隙纵横比谱计算的理论饱和纵、横波速度，与实测数据吻合.图6b中红色部分是反演的0压下的裂隙纵横比谱，其分布在9.5×10-5至1×10-2的纵横比范围内.

图6 Westerly花岗岩反演结果(a) 压力-波速数据(标识符)及其理论拟合(曲线)； (b) 反演得到的孔隙纵横比谱； (c) 本文(红色)与Cheng方法(黑色)反演结果(标识符)与扫描电镜孔隙纵横比谱(直方图)的对比; (d) 硬孔(红色)、裂隙(绿色)孔隙度和裂隙密度(蓝色)随压力的变化.Fig.6 Westerly granite inversion results(a) Laboratory pressure-velocity data (markers) and theoretical fitting (curves) using inversion results; (b)Inverted pore aspect ratio spectrum; (c) The new (red) and existing (black) inversion results (markers) versus SEM pore aspect ratio spectrum (histogram); (d) Stiff (red) and crack (green) porosity and crack density (blue) variation with differential pressure.

Cheng和Toksöz(1979)反演了Westerly花岗岩的孔隙纵横比谱，Hadley(1976)利用扫描电镜技术给出了Westerly花岗岩在纵横比小于1×10-2的孔隙纵横比谱，二者分别由图6c中的黑色数据点和直方图表示，直方图的纵坐标为Hadley用纵横比小于1×10-2的裂隙总孔隙度归一化后的孔隙度，以线性刻度画出.该图同时也显示了新的反演结果(红色数据点)与已有结果的对比.与前述的砂岩结果很不相同的是，无论是实测还是反演的纵横比谱都显示了纵横比约为10-3附近的峰值；并且，新结果与实测数据的吻合大大优于之前的结果.该图的结果表明该致密花岗岩中发育的微裂隙的优势纵横比约为10-3.低纵横比值的裂隙在压力加载下很容易闭合，对应于图6d中岩样的裂隙密度在低压下的陡然下降；低纵横比值的裂隙闭合后，剩余的裂隙随压力缓慢闭合，岩石的波速变化也随之减缓，如图6a所示.

以上的几个不同岩性的反演实例表明：用孔隙与多形态裂隙相互作用的孔、裂隙模型替代低频近似的K-T裂缝模型，可以更好地模拟实验室超声频段的压力-弹性波速数据，由此反演得到的裂隙纵横比谱比原有结果更加符合实际观测结果，有效地提高了从实验室岩芯超声测量获取孔隙分布特征的精度和可靠性.

5 结论

本文将孔隙、裂隙介质弹性波理论扩展到含多形态裂隙岩石的情况，利用多裂隙体系中多个弛豫频率的裂隙挤喷效应的叠加，可以很好地模拟与解释弹性波的频散和衰减在全频域中的变化特征.运用该扩展理论描述岩石中不同形态裂隙在不同压力下的闭合对波速的影响，可以从压力-波速数据反演表征岩石孔、裂隙分布形态的孔隙纵横比谱.与基于K-T低频近似裂缝模型的反演方法相比，新的反演方法的精度和可靠性有较大的提高，结果与岩芯的扫描电镜分析结果符合更好.对多种岩石的应用结果表明：该反演方法对多种岩性的孔隙结构有很好的适用性，为岩石孔隙结构特征的实验分析提供了有效的途径和方法.