多孔可变临界孔隙度模型及储层孔隙结构表征

张佳佳， 曾庆才， 印兴耀， 李红兵， 姜仁，黄家强， 张连群， 张广智， 谷一鹏

1 中国石油大学(华东)地球科学与技术学院， 山东青岛 266580 2 深层油气地质与勘探教育部重点实验室， 山东青岛 266580 3 中国石油勘探开发研究院， 北京 100083

0 引言

随着油气勘探开发程度的不断提高，碳酸盐岩、致密砂岩和页岩等储层已经成为油气增储上产的重点领域.这些储层孔隙类型多样，孔隙结构复杂，非均质性很强，属于典型的多重孔隙储层(含有多种孔隙类型，如孔隙、裂缝和溶洞等).多重孔隙储层预测的难点就是寻找高可采储量和高经济价值的优质储层，而岩石孔隙结构是影响多重孔隙储层油气储量和产能的重要因素(Anselmetti and Eberli，1993；Lucia，1995；Baechle et al., 2005；Li and Chen，2013).因此，岩石孔隙结构表征是多重孔隙储层预测和流体识别的关键.

实验室测量或者测井数据都有丰富的数据或直观的手段来表征储层孔隙结构，例如岩石薄片分析、成像测井等，但是实验室测量或者测井数据很难进行横向预测.而善于进行横向预测的地震数据只有纵波速度、横波速度以及密度三种弹性参数，无法直观的描述储层孔隙结构.更重要的是，多重孔隙储层的弹性参数与孔隙结构之间的关系非常复杂，例如相同孔隙度和流体成分、不同孔隙结构的碳酸盐岩之间的速度差异可以达到2500 m·s-1(Lønøy，2006；Sayers，2008；Baechle et al.，2005，2008，2009).利用地震数据表征储层孔隙结构存在两个核心问题：一个是如何描述储层孔隙结构，另一个是如何建立孔隙结构与弹性参数之间的关系.很多学者在这两方面做了尝试和探索.

岩石物理中多利用孔隙纵横比(等效椭圆的短轴与长轴之比)描述储层孔隙形状，因此最常见的方法是利用孔隙纵横比表征储层孔隙结构，再利用岩石物理模型建立弹性参数与孔隙纵横比之间的关系.其中一类方法是利用总的等效孔隙纵横比定性的描述储层孔隙结构，例如总的等效孔隙纵横比偏大，认为铸模孔等硬孔隙较多，反之偏小则认为裂缝等软孔隙较多.例如Sun等(2012)根据实际测量的纵波速度反演碳酸盐岩储层的平均孔隙纵横比，再对平均孔隙纵横比进行剖分表征储层孔隙结构.李宏兵等(2013)提出了利用Gassmann方程(Gassmann，1951)和DEM解析模型(Li and Zhang，2010，2011，2012)由纵波速度反演等效孔隙纵横比预测横波速度以及描述储层孔隙结构.

另一类方法是利用不同的孔隙纵横比来表征不同的孔隙类型(如裂缝的孔隙纵横比值偏小，为0.01左右，基质孔隙的孔隙纵横比值居中，往往在0.10左右，铸模孔的孔隙纵横比偏大，往往在0.60左右)，再由弹性参数等反演不同孔隙类型的体积含量，相比于总的等效纵横比方法来说更加定量化表征储层孔隙结构.例如Kumar和Han(2005)提出了利用DEM模型(Berryman，1980)和Gassmann方程估算岩石不同孔隙类型的孔隙纵横比以及它们各自体积含量的方法.Xu和Payne(2009)将适用于砂泥岩的Xu和White模型(1995)推广到碳酸盐岩，提出了一种适用于碳酸盐岩的Xu-Payne模型，然后建立岩石纵波速度与岩石不同孔隙类型的孔隙纵横比以及体积含量之间的关系图版.Zhao等(2013)参照Kumar和Han(2005)以及Xu和Payne(2009)的方法，提出了利用地震数据进行岩石物理反演表征岩石的孔隙结构的方法，以此推测碳酸盐岩储层的地质演化历史.

除了利用孔隙纵横比描述孔隙结构之外，很多学者还提出了构建一个新参数表征储层孔隙结构，实际上这些构建的新的孔隙结构参数往往都是与孔隙纵横比有直接或者间接的关系.例如蒋炼等(2011)运用Eshelby(1957)椭球包体裂缝理论和Gassmann方程，推导了碳酸盐岩孔隙度与岩石孔隙纵横比以及岩石弹性参数之间的关系式，用来模拟非均质性较强的生物礁滩储层的孔隙结构，该关系式包含了孔隙纵横比的影响.Sun(2004)提出了一种孔隙结构参数用来描述岩石孔隙结构和孔隙大小对其弹性参数的影响，该孔隙结构参数与孔隙纵横比也存在间接的关系.

在岩石物理模型除了传统的经验模型、理论模型之外，还有一类综合经验模型和理论模型的启发式模型(Heuristic model)在行业中应用广泛(Zhang et al.，2020)，因为启发式模型综合了二者的优点.临界孔隙度模型(Nur，1992)是一种常见的经验模型，因为其具有一定的物理意义和地质含义，并且简单实用，所以在岩石物理建模过程经常使用.很多专家学者也对该模型进行了发展，Markov等(2012，2013)提出的GDEM(Generalized Differential Effective-Medium)模型，在传统的DEM模型的基础上，考虑了岩石的临界孔隙度对岩石弹性特征的影响.临界孔隙度可以用来建立岩石骨架与岩石基质之间的关系，那么能否利用岩石的临界孔隙度来描述储层孔隙结构，是需要解答的问题.

针对这一问题，本文利用临界孔隙度模型和等效介质理论的Kuster-Toksöz方程(Kuster and Toksöz，1974)推导了临界孔隙度与孔隙纵横比之间的关系，进而利用极化(形状)因子建立临界孔隙度与弹性参数之间的关系，构建了能够包含多种孔隙类型的临界孔隙度模型——多孔可变临界孔隙度模型.基于多孔可变临界孔隙度模型，本文提出了一种利用储层的弹性参数反演储层不同孔隙类型的体积含量的方法.

1 多孔可变临界孔隙度模型

1.1 临界孔隙度与孔隙纵横比之间的关系

等效介质理论或者经验模型通常用来描述复杂储层的孔隙结构与弹性参数之间的关系.Kuster和Toksöz(1974)基于长波长一阶散射理论，使得弹性波通过等效介质产生位移场和波经每个包含物(孔隙)散射引起的位移场相同，从而推导出包含N种孔隙类型的岩石的等效弹性模量的表达式.

(1)

(2)

Nur(1992)提出了临界孔隙度的概念，所谓临界孔隙度就是当岩石的孔隙度超过一定临界值时，岩石的组成矿物之间相互分散就不再是岩石了，这时候对应的孔隙度就是临界孔隙度，如图1所示，利用临界孔隙度建立了岩石骨架和岩石基质弹性模量之间的关系：

(3)

(4)

式中，Kdry和μdry分别是岩石骨架的体积模量和剪切模量，φ为岩石的孔隙度，φc为岩石的临界孔隙度.

图1 岩石临界孔隙度示意图(据Nur，1992)Fig.1 Schematic diagram of critical porosity of rock (from Nur, 1992)

Nur(1992)认为同一岩性岩石的临界孔隙度都是相同的，譬如说砂岩取40%，灰岩60%等，岩石的临界孔隙度大小取决于岩石内部结构，它对于颗粒岩石可能居中，对于裂缝岩石可能非常小，而对于泡沫状的岩石可能非常能够大.

岩石孔隙结构通常用孔隙纵横比来表征，即模拟孔隙的等效椭圆的短轴与长轴之比，如图2所示.

结合临界孔隙度模型和Kuster-Toksöz方程可以推导出岩石的临界孔隙度与岩石孔隙结构(孔隙纵横比)之间的关系，具体推导过程详见附录B.

(5)

(6)

式中，φcK为体积模量的临界孔隙度值，φcμ为剪切模量的临界孔隙度值，由式(5)和(6)可以看到岩石的临界孔隙度值是岩石孔隙结构(孔隙纵横比)以及岩石基质体积模量和剪切模量的函数.

假设纯净砂岩(岩石基质体积模量Km=38 GPa和剪切模量μm=44 GPa)含有单一孔隙类型的孔隙，孔隙纵横比分别为α=0.12，孔隙度为φ=0.30，那么岩石的临界孔隙度可以利用式(5)和(6)计算可得到：φcK=0.36，φcμ=0.45，如图3所示.

假设纯净碳酸盐岩(岩石基质体积模量Km=77 GPa和剪切模量μm=32 GPa)含有三种不同类型的孔隙，孔隙纵横比分别为α1=0.01，α2=0.12和α3=0.8，对应为裂缝、基质孔隙、铸模孔，各自的孔隙体积含量分别为x1=0.04，x2=0.2和x3=0.06，那么岩石的临界孔隙度可以利用式(5)和(6)计算可得到：φcK1=0.03，φcK2=0.36和φcK3=0.99，φcμ1=0.05，φcμ2=0.45和φcμ3=0.99，如图4所示.

图4 三种孔隙类型碳酸盐岩的临界孔隙度值Fig.4 Critical porosity values for carbonate with three pore types

1.2 多孔可变临界孔隙度模型

建立岩石临界孔隙度与孔隙结构之间的关系，再将其代入到Kuster-Toksöz理论模型中就可以获得包含多种孔隙类型的临界孔隙度模型表达式.如果考虑包含多种孔隙类型的岩石，各种孔隙类型的体积含量为xi，φcKi分别为各种孔隙的体积含量的临界孔隙度值，φcμi分别为各种孔隙的剪切含量的临界孔隙度值，那么就可以得到包含多种孔隙类型的临界孔隙度模型的体积模量和剪切模量表达式(具体推导过程见附录C)：

(7)

(8)

其中求和是对多种孔隙类型，N为岩石包含的孔隙类型的数目.

式(7)和(8)可以称为多孔可变临界孔隙度模型，新模型可以包含多种孔隙类型，孔隙类型可以用该孔隙类型对应的临界孔隙度表征.不同孔隙类型的临界孔隙度值是不同的，取决于孔隙纵横比以及岩石基质体积模量与剪切模量之比.值得注意的是，由于要满足临界孔隙度的物理意义，即岩石的孔隙度要小于临界孔隙度，所以对于多重孔隙储层来说，每种孔隙类型的体积含量要小于其对应的临界孔隙度，即xi<φci，因此计算过程中需要进行迭代求解，具体迭代过程同与Xu和White模型(Xu and White，1995)一样.

选取Han(1986)75个干燥砂岩中的10个纯净砂岩，利用单孔可变临界孔隙度模型(即式(7)和(8)中的孔隙类型N=1)进行岩石物理建模，其中岩石基质体积模量Km=38 GPa和剪切模量μm=44 GPa，见图5(a和b).从图中可以看，单孔可变临界孔隙度模型只能用单一的孔隙纵横比表征储层孔隙结构，对于储层孔隙结构变化则需要选取不同的孔隙纵横比值，才能够模拟得到较好的岩石弹性模量，并且只是定性描述孔隙结构，无法定量描述孔隙类型的体积含量等参数.

选取Anselmetti和Eberli(1993)77个碳酸盐岩样品，利用双孔可变临界孔隙度模型进行岩石物理建模，即式(7)和(8)中的孔隙类型N=2.这些碳酸盐岩岩石样品的孔隙类型包括致密岩石、粒内孔、粒间孔、铸模孔(>300 μm)、铸模孔(<300 μm)和微孔隙等六种，如图5所示.这里假设有三种孔隙类型在这些碳酸盐岩样品中，硬孔隙、基质孔隙和软孔隙，它们的孔隙纵横比分别为0.6、0.15和0.03，岩石基质体积模量Km=77 GPa和剪切模量μm=32 GPa.这里虽然划分为三种孔隙类型，但是仍然要利用双孔可变临界孔隙度模型，主要原因是目前实际岩石物理模拟过程大多数都假设岩石中包含两种孔隙类型，即同时包含硬孔隙和基质孔隙，或者同时包含基质孔隙和软孔隙(Kumar and Han，2005； Xu and Payne，2009).孔隙类型确定好以后，就需要模拟不同孔隙类型的体积含量的变化，图6显示了77个碳酸盐岩样品的纵波速度随孔隙度以及硬孔体积含量或软孔体积含量变化曲线，从图中可以看到碳酸盐岩的纵波速度随着孔隙度增加以及软孔含量增加或者硬孔含量减小而减小.这就表明利用双孔可变临界孔隙度模型同时包含两种孔隙类型(软孔隙与基质孔隙组合或者硬孔隙与基质孔隙组合)，改变不同孔隙类型的体积含量，就能够模拟得到准确的碳酸盐岩速度.相对于单孔可变临界孔隙度模型，多孔可变临界孔隙度模型可以定量描述储层孔隙类型的体积含量等.

图5 利用单孔可变临界孔隙度模型模拟砂岩弹性模量(a) 体积模量; (b) 剪切模量.Fig.5 Simulating sandstone elastic modulus with single-porosity variable critical porosity model(a) Bulk modulus; (b) Shear modulus.

图6 利用双孔可变临界孔隙度模型模拟碳酸盐岩纵波速度Fig.6 Simulating P-wave velocity of carbonate rocks with the dual-porosity variable critical porosity model

2 基于多孔可变临界孔隙度模型的储层孔隙结构表征

多孔可变临界孔隙度模型可以包含多种孔隙类型，如孔隙、裂缝和溶洞等，不同的孔隙类型的孔隙结构不相同，临界孔隙度值也不相同.因此，可以利用多孔可变临界孔隙度模型进行孔隙结构表征.当储层岩石含有多种孔隙类型的时候，如碳酸盐岩，可以仿照Kumar和Han的思想，利用多孔可变临界孔隙度模型反演不同孔隙类型的体积含量.具体实现步骤为：

(1)选取岩石基质矿物组分的弹性模量和孔隙流体组成成分的弹性模量和密度，计算岩石基质的体积模量Km、剪切模量μm和纵波速度VPm和孔隙流体的体积模量Kfl和纵波速度VPfl：

(9)

(10)

(11)

(12)

(13)

ρfl=Soilρoil+Sgasρgas+Swaterρwater,

(14)

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

(7)采用非线性全局寻优算法修改步骤(3)给定的不同孔隙类型的体积含量，计算新的碳酸盐岩储层岩石骨架的弹性模量，重复步骤(3)—(6)，计算相应的误差，比较误差之间的大小，最后得到最优的不同孔隙类型的体积含量.

图7为利用多孔可变临界孔隙度模型对某碳酸盐岩储层不同孔隙类型体积含量反演的结果.这里假设有三种孔隙类型，硬孔、基质孔隙和软孔，它们的孔隙纵横比α分别为0.8、0.12和0.01.利用多孔可变临界孔隙度模型反演得到三种不同孔隙类型的体积含量，图7f为硬孔孔隙度，图7g为软孔孔隙度.

图8为利用多孔可变临界孔隙度模型对某致密砂岩储层软孔隙体积含量反演的结果.这里假设有两种孔隙类型，基质孔隙和软孔，它们的孔隙纵横比α分别为0.12和0.01.利用多孔可变临界孔隙度模型反演得到两种不同孔隙类型的体积含量，图8g为软孔隙的体积含量，可以看到在3775～3880 m软孔隙的体积含量占总的孔隙度的5%以上，指示了软孔隙含量较多，实际钻井表明该层段裂缝含量较多，渗透性很好，是非常好的“甜点”储层.

3 讨论

需要指出的是，本文的提出的多孔可变临界孔隙度模型对传统的临界孔隙度模型的发展和修正，建立了岩石临界孔隙度与储层孔隙结构(孔隙纵横比)之间的理论关系，是一种综合经验模型和理论模型的启发式模型(Heuristic model).多孔可变临界孔隙度模型是基于等效介质理论Kuster-Toksöz方程推导出来的，虽然Kuster-Toksöz方程要求孔隙彼此之间是孤立的，近似高频假设条件，但是多孔可变临界孔隙度模型计算过程中采用了类似Xu-White模型的计算过程，每一步计算得到的等效介质又被作为新的背景介质，因此这种迭代计算方法既考虑了岩石孔隙之间的相互作用，同时也可以一定程度上解决岩石物理中不同的尺度问题(频率问题).另外，针对岩石不同频率(超声波高频与地震低频)的弹性特征差异还可以利用临界孔隙度值差异来进行等效模拟.

图7 利用多孔可变临界孔隙度模型反演碳酸盐岩储层硬孔孔隙度和软孔孔隙度(a) 密度； (b) 泥质含量； (c) 孔隙度； (d) 含水饱和度； (e) 纵波速度； (f) 反演得到的硬孔孔隙度； (g) 反演得到的软孔孔隙度.Fig.7 Inversion of porosity of different pore types of carbonate reservoir based on multiple-porosity variable critical porosity model(a) Density； (b) Shale volume content； (c) Porosity； (d) Water saturation ；(e) P-wave velocity； (f) Inverted porosity of stiff pore； (g) Inverted porosity of soft pore.

图8 利用多孔可变临界孔隙度模型反演致密砂岩软孔隙的体积含量(a) 密度; (b) 泥质含量; (c) 孔隙度; (d) 含水饱和度; (e) 纵波速度; (f) 横波速度; (g) 反演得到的软孔体积含量.Fig.8 Inversion of volume content of soft pore of tight sandstone based on multiple-porosity variable critical porosity model(a) Density; (b) Shale volume content; (c) Porosity; (d) Water saturation; (e) P-wave velocity; (f) S-wave velocity; (g) Inverted volume content of soft pore.

4 结论

本文结合经验模型临界孔隙度模型和等效介质理论Kuster-Toksöz方程，推导了岩石的临界孔隙度与岩石孔隙结构的关系，并且推导了能够包含多种孔隙类型的临界孔隙度模型，并将之定义为多孔可变临界孔隙度模型.基于多孔可变临界孔隙度模型可以利用储层的弹性参数反演不同孔隙类型的体积含量.实验室测量数据和实际测井数据表明，多孔可变临界孔隙度模型能够适用于多重孔隙储层岩石物理建模和孔隙结构表征.

致谢感谢匿名评阅人对本文提出的宝贵修改意见.

附录A 极化因子P和Q表达式

具有任意孔隙纵横比的椭球状包含物的P和Q系数可写为(Kuster和Toksöz，1974)：

(A1)

(A2)

其中的张量Tijkl将均匀远场应变场与椭球包含物的应变相联系(Wu，1966)，Beryman(1980)给出了计算P和Q所需要的相关的标量：

Tiijj=3F1/F2,

(A3)

(A4)

其中

+B(3-4R)+(A/2)(A+3B)(3-4R)

×[f+θ-R(f-θ+2θ2)],

F4=1+(A/4)[f+3θ-R(f-θ)],

F6=1+A[1+f-R(f+θ)]+B(1-θ)(3-4R)

F7=2+(A/4)[3f+9θ-R(3f+5θ)]

+Bθ(3-4R)

F8=A[1-2R+(f/2)(R-1)+(θ/2)(5R-3)]

+B(1-θ)(3-4R)

F9=A[(R-1)f-Rθ)]+Bθ(3-4R),

其中A，B和R为

A=μi/μm-1

且

R=(1-2νm)/2(1-νm)

函数θ为

以上分别针对扁长和扁圆椭球体，且

注意：对扁圆椭球体α<1，对扁长椭球体α>1.

附录B 临界孔隙度与孔隙纵横比之间的关系

(B1)

当岩石包含单一孔隙类型的时候，则孔隙的体积含量就等于岩石的孔隙度，即x=φ，那么式(B1)变为

(B2)

因此，当φ=φcK时，φcK定义为体积模量的临界孔隙度值，可以得到：

(B3)

式(B3)可改写为

(B4)

这就是体积模量的临界孔隙度值φcK与岩石孔隙结构α之间的关系表达式.

(B5)

当岩石包含单一孔隙类型时，孔隙的体积含量就等于岩石的孔隙度，即x=φ，那么式(B5)变为

(B6)

因此，当φ=φcμ时，φcμ定义为剪切模量的临界孔隙度值，可以得到：

(B7)

式(B7)可以改写为

(B8)

这就是剪切模量的临界孔隙度值φcμ与岩石孔隙结构α之间的关系表达式.

附录C 多孔可变临界孔隙度模型

如果考虑包含多种孔隙类型的岩石，各种孔隙类型的体积含量为xi，φcKi分别为各种孔隙的体积含量的临界孔隙度值，那么有

(C1)

因此

(C2)

(C3)

再将式(C3)代入式(B1)中，就可以得到包含多种孔隙类型的临界孔隙度模型的体积模量表达式

(C4)

其中求和是对多种孔隙类型，N为岩石包含的孔隙类型的数目.

同样的，如果考虑包含多种孔隙类型的岩石，各种孔隙类型的体积含量为xi，φcμi分别为各种孔隙的剪切含量的临界孔隙度值，那么有

(C5)

因此

(C6)

(C7)

再将式(C7)代入式(B5)中，就可以得到包含多种孔隙类型的临界孔隙度模型的剪切模量表达式

(C8)

其中求和是对多种孔隙类型，N为岩石包含的孔隙类型的数目.