基于数字岩心孔喉特征的等效孔隙纵横比有效性研究

谭开俊， 赵建国， 滕团余， 刘欣泽， 闫博鸿

1 中国石油勘探开发研究院西北分院， 兰州 730020 2 中国石油大学(北京)油气资源与探测国家重点实验室, 北京 102249 3 中国石油天然气股份有限公司玉门油田环庆分公司， 甘肃酒泉 735019

0 引言

岩石微观孔隙结构形态对岩石物性具有重要的影响，以定量化储层预测为目标，油气勘探中对岩石微观孔隙结构研究的要求日益提高.碳酸盐岩储层由于成岩后期强胶结、溶蚀、以及矿物充填等作用次生孔隙非常发育，这使其孔隙结构极其复杂(Bathurst, 1972; Lucia, 1983)，常见的有铸模孔、粒内溶孔、粒间溶孔、晶间孔、裂缝等，Anselmetti和Eberli(1993, 1997)发现这些复杂多变的孔隙结构对弹性性质的控制作用甚至比孔隙度还大；Zhang等(2015)研究发现在中国四川碳酸盐岩储层中，在相同孔隙度下由于多种孔隙结构的发育其纵波速度会出现约有1000 m·s-1的速度差异；Sayers(2008)也指出具有相同饱和度、孔隙度的数据点在解释量版上的分布并不集中.这导致我们在储层预测时很容易将孔隙结构对速度、弹性性质造成的影响归结在孔隙度、孔隙流体等其他因素上，使结果产生很大的偏差(李宏兵等，2019).另一方面，Anselmetti和Eberli(1999)、Eberli等(2003)还通过分析裂缝与孔洞的分布发现裂缝发育会显著增强渗透率这一规律，并试图用之来近似预测岩石渗透率.因此，孔隙结构特征的研究对碳酸盐岩等岩性的储层预测、含油性评价具有重要意义.

为了尽可能定量化地预测与表征孔隙结构，国内外许多学者从多种角度进行了大量的实验与理论研究.从岩石物理实验研究的角度出发，Anselmetti等(1998)通过图像数字分析技术(Digital Image Analysis，DIA)提取了17块白云岩铸体薄片的孔隙尺寸分布、孔隙形状(孔隙周长与孔隙面积的比值)等作为孔隙结构表征参数，随后又使用Wyllie时间平均方程计算速度与实测速度的差值为标准，定义了300多块碳酸盐岩岩心的孔隙结构变化情况，最后将孔隙结构用于渗透率的预测(Anselmetti and Eberli，1999；Eberli et al.，2003).Verwer等(2010)分析指出具有较高孔隙比表面的碳酸盐岩样品剪切模量会明显降低.目前使用高分辨率的扫描电镜还能够更进一步研究纳米级的微小孔隙，三维全岩心CT成像技术可以更准确地获得岩石内部骨架与孔隙的复杂几何形状(Coenen et al.，2004)，建立真实的岩石模型.此外随着实验手段的进步，Volokitin等(2001)将核磁共振技术应用于岩心的孔隙结构描述：首先通过离心过程获取岩心处于不同饱和情况时的横向弛豫时间T2谱，再由一定的标定方法与理论转换得到孔喉半径分布频率.还可以对岩心进行恒速或高压压汞，通过进汞与退汞压力变化曲线与相应理论依据能够较准确地推算孔隙、喉道、孔喉比大小及含量分布，实现孔隙结构定量化研究(Wang et al.，2019；王永诗等，2021).

在岩石物理理论模型方面，众多学者通常采用等效介质理论将孔隙结构与弹性性质联系起来.经典的等效介质理论主要有K-T模型(Kuster and Toksöz，1974)、Mori-Tanaka(MT)模型(Benveniste，1987)和微分等效介质(Differential Equivalent Medium，DEM)模型(Berryman et al.，2002)等，这些模型均将孔隙理想化为纵横比不同的椭球体，在此基础上计算多孔岩石等效弹性模量.Xu和White(1995)考虑碎屑岩中的孔隙在砂岩组分与泥岩组分中具有不同的结构(孔隙纵横比)，以DEM、K-T模型为基础发展了著名的Xu-White模型，用以估算碎屑岩储层的横波速度.Kumar和Han(2005)引入了大、中、小三种纵横比的孔隙，给出了联合DEM模型与Wyllie时间平均方程估算三种孔隙的纵横比及体积分数的实施流程.Xu和Payne(2009)应用Kumar的工作思路并加入泥岩孔隙，对碳酸盐岩测井曲线进行了三种孔隙结构及含量的预测，其对孔隙结构的刻画通过实测横波速度与预测值的对比得到验证.Zhao等(2013)将Kumar的孔隙分类及计算方法应用于地震数据，并通过电阻率成像测井对孔隙结构进行了验证.Pan等(2015)、赵建国等(2021a，b)基于数字岩心也对地震孔隙结构预测方面进行了研究.而李宏兵等(2013)认为在研究孔隙结构对弹性性质的影响时，即使岩石中含有多种孔隙类型都可将其等效为具有单一纵横比的孔隙来表示，并通过具有多重孔DEM模型与单一孔隙类型DEM模型之间的对比证明了这种等效的合理性.除了上述将岩石中孔隙等效为一种或有限的几种结构外，David和Zimmerman(2012)考虑了岩石中含多种孔隙的情况，并通过变压力超声干燥测量速度数据与DEM理论反演得到了具有不同纵横比孔隙的分布，结果通过饱和测量数据进行验证.国内欧阳芳等(2021)也通过超声实验与多种等效介质理论对此方面进行了分析.然而这种方法很难应用于测井、地震等大尺度的情况.

总体上目前对孔隙结构的描述可以分成两大类，一种是从实验角度获取，包括薄片、CT扫描数字岩心、扫描电镜等直接成像手段，以及通过核磁共振、压汞等物理实验方法得到孔隙、喉道半径等微观几何形状参数；另外一种是将预先获取的弹性模量联合相应的等效介质理论，通过反演的思路计算等效孔隙纵横比或类似参数，相比其他孔隙结构参数而言基于等效介质理论孔隙纵横比更容易获得，并且能够应用于测井、地震等大尺度数据.然而，等效孔隙纵横比毕竟是从弹性性质出发通过理论间接计算得到的，能否确切地表征岩石真实的复杂孔隙结构特征，在于其是否能够与岩石中直接观测到的孔隙几何形状参数之间建立起定量的关系(Baechle et al.，2008).虽然Anselmetti等(1998)、Baechle等(2008)使用DIA技术能够提取铸体薄片的孔隙结构参数，可以与DEM模型中的孔隙纵横比建立一定的联系，但是铸体薄片毕竟属于二维图像，与三维岩心实际的孔隙空间相比有很大的差距，且不同位置的薄片结构变化或许很大.此外，使用薄片的孔隙结构参数提取缺乏明确的孔喉分割理论支撑，可能不够精准.因而如果能够定量化地建立孔隙纵横比与实际岩石三维图像的孔隙特征参数之间的联系，那么在利用等效介质理论表征孔隙结构时，其合理性与准确性就具有了强有力的支撑，CT扫描技术无疑是最合适的手段.

本文以CT扫描数字岩心为基础，从模拟的弹性性质和图像提取的孔隙网络几何特征两方面研究了孔隙结构的定量化表征方法，建立了等效孔隙纵横比与孔喉特征参数间的定量关系：我们首先分别对基于数字岩心的弹性模拟技术与孔隙网络提取方法进行介绍；然后，对砂岩与碳酸盐岩两类不同岩性岩心进行CT扫描，并对扫描图像进行图像处理以保证能精确重构三维数字岩心模型；接下来将数字岩心分割为子块，应用上述两种方法获得每个子块的等效弹性模量与多种孔喉几何特征参数；最后，构建等效介质理论与模拟弹性模量之间的损失函数方程，通过优化算法来计算等效孔隙纵横比，再分别与各子块的不同孔喉几何特征参数进行了交会分析， 最终发现平均喉道长度和等效孔隙纵横比有较为明显的线性拟合关系，表明它们对孔隙结构的表征具有良好的一致性.

1 基于数字岩心的弹性模拟与孔隙网络提取

对于物体弹性性质的数值模拟，我们使用的是Garboczi和Day(1995)提出的有限元模拟方法，它使用了静态线弹性方程的变分形式.在数值计算时首先给定一初始位移场，利用共轭梯度算法对线弹性方程离散后的有限元控制方程进行迭代求解，计算三维数字岩心数据体内任意一体素的位移量，由位移得到应力应变值并平均得到宏观量，再使用胡克定律即可计算岩石等效弹性模量.这种方法以弹性模量会影响物体中的应力应变分布为前提，根据最小位能原理，物体在受到应力产生应变后，会向着应变位能减小的趋势变化，最终在应变位能最小时达到平衡状态，利用此时物体的应力应变，计算得到的是物体处于小应变状态下的等效弹性模量，因此是一种静态的模拟.

从图形几何学的角度，在利用数字岩心图像直接确定孔隙结构特征方面，我们使用了孔隙网络法.这主要是因为在理想情况下，我们期望岩心内部的孔隙能够通过联通、分割算法(夏良正和李久贤，1999)将每部分连在一起的孔隙空间单独标定出来，这样就可以很容易确定孔径、球度、表比面积等几何形状特征参数.然而大多数联通在一起的孔隙空间都是极度不规则的，很难定量评价形状特征，或者有些岩心中的孔隙空间是整体联通的，无法对其进行分割与评价，因此需要采用特定的方法将较为复杂的孔隙空间转化成一种形状比较规则的空间来等价表示.受到渗流性质模拟领域的启发，在研究中我们采用最大球方法提取了三维数字岩心的孔隙网络.这种方法首先将不规则的孔隙空间转化成了由大量内切球体充填的空间，对所有球体进行一定的划分后将孔隙空间中较宽的区域定义为孔隙，孔隙之间的窄区域通过圆柱型的喉道相连接，两者组合成为孔隙网络.这种网络能够近似等效地表示岩心的孔隙空间，因此其几何形状特征参数就能够表示孔隙结构特征(Dong and Blunt, 2009).

(1)

(2)

其中，ra、rb表示孔隙对应的最大球a、b的半径，rt表示喉道最大球的半径，k表示孔喉长度分割系数，喉道长度lt定义为喉道总长度lab减去两个孔隙长度la、lb后的值.以上就能够得到孔隙网络所有的几何特征参数. 从式(1)、(2)中可以发现，这种孔喉划分方法是在孔、喉对应的最大球半径与分割系数的约束下分别定义了孔隙与喉道的范围，在喉道长度范围内的所有最大球都归属于喉道，而在孔隙长度范围内的所有最大球则归属于孔隙.这种方法虽然能够简便的获得喉道与孔隙的长度特征，但是其中分割系数是一个人为控制的比例系数，Dong和Blunt (2009)在计算时给定为0.6，在本文中为了分析该参数对孔喉分割的影响，我们分别计算了不同系数值时的喉道长度参数并进行对比.

图1 孔隙-喉道划分示意图Fig.1 Pore-throat segmentation

为了更好地理解孔隙空间转换为孔隙网络后孔、喉位置对应的孔隙空间区域，为孔隙结构的定量分析提供基础，我们构建了由球体规则排列而成的岩石模型并提取孔隙网络，见图2.孔隙网络中8个较宽的区域被定义为孔隙(红色球体)，孔隙间的连接区域为喉道(蓝色圆柱体)，这有效地实现了不同孔型孔隙空间的定量判别.为进一步分析孔隙网络模型是否能够对岩石宏观整体的孔隙结构有比较定量化的表征，我们尝试提取了有明显孔隙结构特征岩心的孔隙网络.图3a中的数字岩心样品为典型的裂缝性孔隙结构，孔隙空间成扁平状延展，图3c中为孔洞型的数字岩心样品，其孔隙空间多为连续的类球状，延展小数量多，分布较为均匀；图3b与3d为对两者提取的孔隙网络模型，可以发现裂缝型样品对应的孔隙网络模型中，蓝色圆柱状喉道的长度大多相比孔洞型样品的喉道要明显更长一些，裂缝几乎都由喉道组成.而孔洞型样品的喉道较为短小，基本由孔隙与喉道连接而组成.统计喉道长度、半径频率分布谱，发现两种喉道相关参数在裂缝型样品中都要大于孔洞型样品，分布范围也更宽，但是其中喉道长度差别更明显，见图4、5，因此在下文中我们重点应用喉道长度对孔隙结构进行定量化分析.

2 数字岩心的图像预处理

为了分析与验证孔隙网络模型中的孔喉参数能否定量化描述岩石的孔隙结构，我们分别选择了孔隙结构较为复杂的碳酸盐岩样品和孔隙结构相对均一的砂岩样品，共同来进行孔隙结构定量化分析与对比.所有样品都为切割好的直径38 mm的圆柱体，如图6、7所示，使用本实验室的CT设备扫描成像后得到的图像分辨率可达每体素20.7678 μm.为了提取后续运算所需的孔隙结构信息，CT扫描原始图像需要经过对比度增强、滤波去噪、相边缘增强以及图像二值化分割几个主要的图像处理步骤，本文以图6中1号白云岩岩心为例来展示图像处理的过程及效果，见图8，处理仅对从原始图像中心剪裁的边长为1000体素的正方体图像进行.

图2 最大球法提取简单球体堆积模型的孔隙网络(a) 球体堆积模型骨架空间； (b) 球体堆积模型孔隙空间； (c) 球体堆积模型孔隙网络.Fig.2 Process of pore-network extracted by maximal ball method for a simple sphere packing model(a) Skeleton of sphere packing model； (b) Pore space of sphere packing model； (c) Pore-network of sphere packing model.

图3 不同孔隙类型数字岩心样品孔隙网络提取(a) 典型裂缝型样品数字岩心模型； (b) 裂缝型岩心孔隙网络模型； (c) 典型孔洞型样品数字岩心模型； (d) 孔洞型岩心孔隙网络模型.Fig.3 Extraction of pore-network for digital core of different pore structure(a) Digital core model of typical fracture core； (b) Pore-network model of fracture core； (c) Digital core model of typical porous core； (d) Pore-network model of porous core.

图4 裂缝、孔洞型样品喉道长度分布谱Fig.4 Histogram of throat length in fracture and porous cores

图5 裂缝、孔洞型样品喉道半径分布谱Fig.5 Histogram of throat radius in fracture and porous cores

对数字岩心的图像处理是进行弹性性质模拟及孔隙网络模型提取的关键步骤，直接决定着孔隙结构定量化表征的准确性.图6中第二排展示了7块碳酸盐岩样品CT扫描成像处理后的三维孔隙空间分布情况，可以发现该类样品具有比较复杂的孔隙结构类型、非均质性较强、孔隙度都很小，图7中第二排展示了砂岩岩心的孔隙三维分布，直观看来砂岩样品的孔隙结构明显比较单一且分布均匀，和碳酸盐岩样品的区别较大，孔隙度也远大于碳酸盐岩.

图6 碳酸盐岩样品Fig.6 Carbonate sample

图7 砂岩样品Fig.7 Sandstone sample

图8 CT扫描图像处理中间过程及结果(a) 原始扫描图像切片； (b) 对比度增强结果； (c) 滤波； (d) 二值化分割结果.Fig.8 Result and process of CT image processing(a) Original scanned image slice； (b) Contrast enhancement result； (c) Filtering for image； (d) Binarized image segmentation result.

在后续的模拟运算前，我们还需对整个数字岩心数据进行子块分割，这主要有两方面的原因，首先是受限于计算机的处理能力、内存和运算所需要的时间，一般较难处理10003个体素的数据量；另外单个或少量岩心计算结果较为单一，无法分析岩心性质的统计性规律.对此Dvorkin等(2011)认为如果将数字岩心内部子块数据集的模拟数据点组成趋势线，那么这将在多个尺度上都比较稳定，并且这样的趋势线不会随岩心取样空间位置和尺度的变化而产生剧烈波动， 因此我们需要将大正方体的岩心分割成大小合适的子块数据.但是分割后的岩心子块也应足够大，这样才能包含足够充分的孔隙信息，与真实岩心的孔隙结构特征也更相近.姜黎明等(2012)通过代表体积元(REV)方法测试发现，岩心块尺寸选在2003时通常能够包含较丰富的孔隙结构且计算较快，因此本文使用了该尺寸进行子块分割，之后就可以对子块岩心进行弹性性质模拟与孔隙网络提取.

在弹性性质模拟时还需给定岩心中骨架与孔隙部分各自的弹性模量， 那么根据等效介质理论的思想，岩石的等效弹性性质除孔隙结构、孔隙度外，还会受到骨架的矿物组成与岩石饱和流体情况两个因素的影响.其中碳酸盐岩的骨架矿物组成较为纯净，常为方解石或白云石.本文样品经XRD(X-ray diffraction，X射线衍射)矿物分析发现白云石占90%以上，矿物组成单一，因此在比较所有碳酸盐岩弹性模量时骨架矿物组成对各样品实际弹性性质的影响应很小，在此我们给定骨架部分为白云岩的弹性参数(体积模量为80 GPa，剪切模量为58 GPa，该模量通过对同井位中、孔隙度0.1%左右的一非常致密岩心测量得到)；与之相反砂岩骨架的矿物组成通常比较复杂，因此在对比砂岩弹性性质时矿物成分是重要的影响因素，但是考虑到我们的研究核心在于孔隙结构的定量化表征，孔隙结构对弹性性质的影响是我们主要考虑的因素，所以仍旧需要假定岩石骨架矿物成分单一，这样所有样品的模拟弹性模量才仅会受到孔隙结构的影响，而不会受不同样品矿物组成不同所影响.考虑到较纯净的砂岩(如枫丹白露砂岩)基本为石英，本文砂岩经XRD分析石英含量也基本占主导，因此给定砂岩的骨架部分为石英的弹性参数(体积模量为37 GPa，剪切模量为44 GPa；Mavko et al.，1998).同样为了排除孔隙流体对样品模拟弹性模量的影响，我们假设两种岩性样品的孔隙中都完全饱和空气.

3 等效孔隙结构定量化表征对比

根据研究思路，我们首先联合模拟弹性模量与等效介质理论估算等效孔隙纵横比，以该参数表征数字岩心整体等效的孔隙结构特征.目前已经发展了多种常用的以孔隙纵横比来表征孔隙形状的等效介质模型，它们的假设条件、适用情况以及加入孔隙的方式都有所不同，考虑到计算效率本文使用了较为简单的显式MT模型，但是为了对比不同理论计算的等效孔隙纵横比是否会对孔隙结构表征产生显著的影响， 我们也给出了DEM模型的计算结果，见附录B.通过等效孔隙纵横比分布的对比表明，等效介质理论不同，对最后孔隙结构刻画造成的影响很小.便于计算的包含N种包裹物的MT模型表示为式(3)(Benveniste，1987)：

(3)

图9a展示了所有碳酸盐岩数字岩心模拟纵波速度与孔隙度交会特征，可以发现交会图中数据点分布总体比较分散，在孔隙度为10%左右时纵波速度的差异高达2500 m·s-1左右，这只能是不同的孔隙结构引起的；而在砂岩数字岩心模拟纵波速度与孔隙度交会图中数据点非常集中，呈明显的近线性分布，同一孔隙度下的纵波速度差距很小，见图9b.在求取等效孔隙纵横比时将岩心孔隙度代入式(3)中，计算的等效模量作为输出，以模拟弹性模量Kmodeling、μmodeling作为约束值作差，构建包含体积模量与剪切模量两个约束条件下的方程组(4)，通过优化算法求解就可以得到当前弹性模量与孔隙度条件下对应的等效孔隙纵横比.式(4)中OF表示目标函数.

(4)

图9中每个数据点的颜色代表该点等效孔隙纵横比，颜色越偏暖色值表示纵横比值越大，反之纵横比值越小.图中蓝色虚线表示理论计算的不同等效孔隙纵横比时纵波速度随孔隙度变化曲线.在碳酸盐岩样品中MT模型计算的等效纵横比分布在0.03～0.4范围内，分布区间比较宽(见频率直方图10中红色所示)，同一孔隙度下都表现为速度越快等效孔隙纵横比越大，总体上所有数据点可以根据颜色分为下部的冷色值与上部的暖色值两类；砂岩样品的等效纵横比分布在0.2～0.5之间，分布区间都要小于碳酸盐岩(见频率直方图10中蓝色所示)，更加集中.

图9 数字岩心纵波速度-孔隙度交会中等效孔隙纵横比分布对比(色标表示孔隙纵横比)(a) 碳酸盐岩； (b) 砂岩.Fig.9 Distribution of effective pore aspect ratio of digital core in cross-plot between P-velocity and porosity (color represent pore aspect ratio)(a) Carbonate； (b) Sandstone.

图10 砂岩与碳酸盐岩等效孔隙纵横比频率分布直方图Fig.10 Histogram of effective pore aspect ratio in sandstone and carbonate

在使用岩心孔隙网络模型提取的孔喉参数对孔隙结构进行表征时，考虑到图4、5中的分析，并且喉道此时更能表现孔隙空间中窄区域的特征，从而影响孔隙形状，我们提取了每块数字岩心所有喉道的平均长度、长度排序后的中间长度、经过喉道体积加权平均的喉道长度和长度分布频率直方图峰值4个喉道长度相关的统计性参数作为代表每块岩心的等效喉道长度，并与等效孔隙纵横比(基于MT模型计算)进行了交会分析.在图11中孔喉长度分割系数为0.6时可以明显发现碳酸盐岩的不同喉道长度参数相比砂岩都具有更宽的分布范围，这种特征和等效孔隙纵横比的分布特征是一致的；并且碳酸盐岩各喉道长度参数与等效孔隙纵横比交会都出现了近似线性趋势，其中喉道平均长度的相关系数明显更大、相关性更强，数据分布更集中，表现为纵横比增大对应喉道长度减小的负相关性.这可以解释为纵横比越小的孔隙越呈扁平状分布，延展更长，因此在提取的孔隙网络中会有更长的喉道长度，而纵横比越大的孔隙越接近球体，对应的喉道会缩短.图12将两类岩心单独分析发现砂岩平均喉道长度与等效孔隙纵横比交会无明显相关性，而碳酸盐岩的交会关系相比图11c中两类岩心时基本不变，因此岩性并不会影响两参数的相关性，图11c中的交会关系对两种岩性都是适用的.随后为了分析孔喉分割系数的影响，我们使用0.4与0.8两个不同系数重新计算了数字岩心喉道平均长度，并与等效孔隙纵横比进行了交会分析，见图13.其中虽然长度值的范围发生了变化，但是所有样品数据点的相对分布规律并没有明显的改变，这表明分割系数对我们的孔隙结构表征研究敏感度较小.根据以上交会分析我们选择了平均喉道长度作为孔隙结构的表征参数，同样使用纵波速度与孔隙度做交会，每个数据点颜色表示平均喉道长度值(由于喉道长度非常小，为了便于比较，我们将所有样品的平均喉道长度以其中的最大值为基准进行了归一化)，颜色越偏冷色值代表喉道越短，反之喉道越长.孔隙结构表征结果见图14所示，蓝色虚线为MT模型计算的不同等效孔隙纵横比时纵波速度随孔隙度变化曲线，作为以喉道长度来定量表征孔隙结构的参考.我们能够明显看出图14中数据点颜色分布与图9中相似，碳酸盐岩能够大致分为喉道较短的冷色点与较长的暖色点两类，砂岩则都表现为比较单一的喉道长度.

图11 孔喉长度分割系数0.6时各喉道长度参数-等效孔隙纵横比交会(a) 加权喉道长度-等效孔隙纵横比交会； (b) 喉道长度频率直方图峰值-等效孔隙纵横比交会； (c) 喉道平均长度-等效孔隙纵横比交会； (d) 喉道长度中值-等效孔隙纵横比交会.Fig.11 Cross-plot between effective pore aspect ratio and parameters which described throat length when pore-throat segmentation coefficient is 0.6(a) Effective pore aspect ratio and weighted throat length； (b) Effective pore aspect ratio and peak value of throat length histogram； (c) Effective pore aspect ratio and average value of throat length； (d) Effective pore aspect ratio and mid-value of throat length.

图12 不同岩心中孔喉长度分割系数0.6时喉道平均长度-等效孔隙纵横比交会(a) 碳酸盐岩； (b) 砂岩.Fig.12 Cross-plot between effective pore aspect ratio and average value of throat length when pore-throat segmentation coefficient is 0.6 in different core(a) Carbonate； (b) Sandstone.

图13 平均喉道长度-等效孔隙纵横比交会(a) 孔喉长度分割系数0.4； (b) 孔喉长度分割系数0.8.Fig.13 Cross-plot between effective pore aspect ratio and average value of throat length(a) Pore-throat segmentation coefficient is 0.4；(b) Pore-throat segmentation coefficient is 0.8.

图14 纵波速度-孔隙度交会图中数字岩心喉道平均长度分布(色标表示归一化喉道平均长度)(a) 碳酸盐岩； (b) 砂岩.Fig.14 Distribution of average length of throat from digital core in cross-plot between P-velocity and porosity (color represent normalized average throat length)(a) Carbonate； (b) Sandstone.

4 讨论与误差分析

经过对以上两种孔隙结构表征思路结果的对比与分析，证明虽然等效孔隙纵横比是通过弹性性质这一间接属性来定量化描述孔隙结构的，但是其对孔隙结构的表征效果与直接从岩心图像中提取的喉道长度的表征效果是相近的.然而目前我们的岩心样品中等效孔隙纵横比与平均喉道长度的交会关系并不是非常理想，因此需要对两参数可能产生的误差进行分析.这之中等效孔隙纵横比由等效介质理论与弹性模量计算得到，其误差应来自于两个方面：一个是弹性模拟的结果是否相对准确，另一个则是等效介质理论(MT模型)能否准确描述弹性模量和等效孔隙纵横比之间的关系.对此我们建立已知孔隙纵横比的模型介质，通过对比其模拟模量与等效介质理论计算模量间的差距来评估误差.首先使用位置随机分布、半径固定的球体作为孔隙，建立理论上孔隙纵横比为1的多孔介质；再使用4参数随机生长法(Quartet Structure Generation Set，QSGS；Wang et al.，2007)，通过控制生长参数生成长轴与短轴比不同的椭球体作为孔隙(同一个模型介质中每个椭球体长短轴都相同)，建立不同孔隙纵横比的多孔介质，具体的长、短轴长度可以通过对图像进行测量得到.以上多孔介质不考虑孔隙间的连通性，各个孔隙互不相交，见图15.

模拟使用石英作为骨架，孔隙中充填空气.图16中不同标记点表示模型介质的模拟模量，粉色虚线表示MT模型计算的弹性模量随孔隙度变化曲线，两者所用孔隙纵横比互相对应.对比发现两模量相差较小，这表明对于简单的、孔隙结构单一且分布均匀的多孔介质，弹性模拟结果是可靠的，等效介质理论也基本能准确计算弹性模量.但是在孔隙结构复杂时等效介质理论的误差很难评价，这是因为很难像图15中的简单模型一样通过测量孔隙长短轴定量化提取复杂孔隙的等效纵横比，而只能由弹性模量反演计算一个理论值.因此裂缝型孔隙定向排列引起的各向异性、孔隙位置与区域的非均匀分布都可能使岩石不满足等效介质理论的假设，从而使等效孔隙纵横比无法反映真实的孔隙结构特征，导致图11c中的拟合误差.另外等效孔隙纵横比也无法反映孔隙之间的联通特性，忽略了孔隙结构对渗流性质的影响.

喉道长度参数的误差主要来源于从孔隙空间提取孔隙网络的过程.通过图2对一球体颗粒堆积形成的简单孔隙空间提取孔隙网络后发现，虽然最大球法对孔隙空间中局部宽区域的孔隙与局部窄区域的喉道有明确的定义方法，但是对两者中间剩余空间的归属还是比较模糊的，这也就导致了确定的喉道长度与孔隙长度可能并不准确.在上文中我们分析了最大球法分割孔喉时分割系数对喉道长度的影响，发现该参数基本只会影响喉道长度值的整体范围，对所用样品喉道长度分布的相对趋势与规律性没有影响.但是孔隙网络最开始是以研究渗流特性为目的而展开的，为了简便而将孔隙空间转化为孔隙网络，并认为孔隙网络能够近似代替实际的孔隙空间，这种近似可能不会对渗流特性产生明显影响，但是对原本孔隙结构的表征肯定存在误差；此外网络模型也相当于将原本联通在一起的孔隙空间分割开来，即使分割过程有最大球法作为理论依据，分割的位置也可能会对实际孔隙结构的刻画产生影响，这些因素都可能导致同一喉道长度下出现多个等效孔隙纵横比，使两者间的关系出现误差.因此在分析等效孔隙纵横比定量化表征复杂孔隙结构的误差方面，我们对弹性模拟、图像形状的定量化表征、孔隙网络提取及孔喉分割还需要更多的研究和尝试，对于孔隙结构复杂、非均质性强的碳酸盐岩也需要更多的研究样品，才能获得更加可靠的等效孔隙纵横比与平均喉道长度间的经验关系.

图15 不同单一孔隙纵横比模型(a) 球状孔隙模型，孔隙纵横比为1； (b) 椭球状孔隙模型， 控制孔隙纵横比定量变化.Fig.15 Model in different single pore aspect ratio(a) Sphere pore model, pore aspect ratio is 1； (b) Ellipsoid pore model, pore aspect ratio is designed quantitatively.

图16 不同单一孔隙纵横比模型模拟弹性模量与等效介质理论计算模量对比(图中不同符号点代表 不同单一孔隙纵横比模型的模拟模量，虚线为MT模型计算模量，数字为对应孔隙纵横比)(a) 体积模量； (b) 剪切模量.Fig.16 Comparison of simulated elastic modulus with calculated modulus by effective medium theory for different single pore aspect ratio model (Different symbolic points represent the simulated modulus of different single pore aspect ratio models, the dotted line represents calculated modulus by MT model, and the number represents the corresponding pore aspect ratio)(a) Bulk modulus； (b) Shear modulus.

5 结论与展望

本文从两种不同角度出发获取岩石孔隙结构特征的表征与描述结果，并对比了两种方法对孔隙结构进行表征的相似性.其中，基于弹性性质并结合等效介质理论计算等效孔隙纵横比的方法是目前最常用的孔隙结构描述方法，在地震以及测井数据的孔隙结构预测上已有比较多的应用，该方法并不是直接从孔隙空间的几何特征中求得孔隙结构参数，而是通过等效介质理论将这种几何特征间接与岩石弹性性质联系了起来.本研究中，借助于数字岩心可以提取储层岩石的孔隙网络模型，并由此出发提取具有明确地质意义的喉道长度，作为表征孔隙结构的几何参数.我们的研究发现，弹性性质联合等效介质理论估算的等效孔隙纵横比参数与直接从数字岩心提取的喉道长度这一几何参数，二者具有很好的相关性，喉道长度长对应于等效孔隙纵横比小的扁平裂缝型孔隙，反之对应于等效孔隙纵横比大的球形孔洞型孔隙.本研究在一定程度上强有力地证明，弹性性质联合等效介质理论这种间接反演等效孔隙纵横比的思路，对形状复杂的孔隙结构进行定量化表征具有准确性与合理性.相比在地震或测井孔隙结构预测时，仅使用成像测井、少量岩心薄片或横波预测误差等来验证等效孔隙纵横比对孔隙结构描述的效果，本文提供了更加定量化、更加直接的证据.此外，由于等效孔隙纵横比参数与喉道长度之间良好的相关性，而喉道特征通常是描述储层渗流特性的重要参数，这使得将基于测井数据或地震数据预测的等效孔隙纵横比属性体转化为更具开发价值的渗流属性体成为可能.

附录A

对于孔隙纵横比α小于1的极化因子，可以表示为(Mavko et al.，1998)：

式中，(Km,Gm)为骨架的弹性模量；(Ki,Gi)为孔隙包含物的弹性模量；α为孔隙纵横比.

附录B

包含N种孔隙的DEM模型可以表示为微分方程式(B1)(Berryman et al.，2002)：

(B1)

图B1与图B2中每个数据点的颜色代表该点孔隙纵横比，颜色越偏暖色值表示纵横比值越大，反之纵横比值越小，图中蓝色虚线表示理论计算的不同孔隙纵横比时纵波速度随孔隙度变化曲线.对比发现不同理论计算的结果有所区别，DEM模型计算的孔隙纵横比较MT模型偏大一些，但是不同岩性孔隙纵横比的分布特征是一致的.在碳酸盐岩样品中DEM模型计算的孔隙纵横比分布在0.04～0.5范围内，MT模型计算的纵横比分布在0.03～0.4范围内，不同模型的孔隙纵横比分布区间都比较宽(见频率直方图B3a、b中红色所示)，同一孔隙度下都表现为速度越快孔隙纵横比越大，总体上所有数据点可以根据颜色分为下部的冷色值与上部的暖色值两类；砂岩样品DEM模型计算的纵横比在0.3～0.5之间，MT模型则在0.2～0.5之间，不同模型孔隙纵横比的分布区间都要小于碳酸盐岩的分布区间且纵横比整体较大(见频率直方图B3a、b中蓝色所示)，从颜色上所有的数据点都为暖色，没有明显的冷色值.以上对比表明即使使用的等效介质理论不同，也不会对孔隙纵横比的计算及最后对孔隙结构的表征产生显著影响.

图B1 纵波速度-孔隙度交会图中碳酸盐岩数字岩心等效孔隙纵横比分布(色标表示孔隙纵横比)(a) 基于DEM模型； (b) 基于MT模型.Fig.B1 Distribution of effective pore aspect ratio of carbonate digital core in cross-plot between P-velocity and porosity (color represent pore aspect ratio)(a) Calculated by DEM model； (b) Calculated by MT model.

图B2 纵波速度-孔隙度交会图中砂岩数字岩心等效孔隙纵横比分布(色标表示孔隙纵横比)(a) 基于DEM模型； (b) 基于MT模型.Fig.B2 Distribution of effective pore aspect ratio of sandstone digital core in cross-plot between P-velocity and porosity (color represent pore aspect ratio)(a) Calculated by DEM model； (b) Calculated by MT model.

图B3 砂岩与碳酸盐岩等效孔隙纵横比频率分布直方图(a) MT模型计算； (b) DEM模型计算.Fig.B3 Histogram of effective pore aspect ratio in sandstone and carbonate(a) Calculated by MT model； (b) Calculated by DEM model.