基于三维CT扫描图像的致密砂岩孔隙数据提取方法

孙 泽

(华北科技学院 安全工程学院，北京 东燕郊 065201)

0 引言

近年来，随着三维重构技术的发展，研究的不断深。致密砂岩作为一种非均匀性的多孔介质，其孔隙结构特征的研究对于深入理解致密砂岩岩石内部流体的储存和运移机制具有重要意义，也是决定其孔渗特征的重要因素[1]。基于现阶段国内外学者对致密砂岩微观孔隙结构的研究技术，比如铸体薄片、扫描电子显微镜(SEM)、压汞法(MICP)、核磁共振(NMR)、物性测试分析等，这些技术或方法都只能反映出微观孔隙的二维结构信息，难以获取三维分布及孔隙间连通性等信息[1]，为了得到测试结果，部分方法会破坏砂岩样品本身结构，如磨制铸体薄片，扫描电子显微镜(SEM)要喷金，压汞测试后岩样中会残留重金属有害物质等等，使采集的岩样无法进行其他实验。三维CT优点在于，它不仅能够在不损坏岩样的条件下，对复杂孔隙结构进行360°检测扫描，而且扫描过后的岩样仍可进行其他试验。美国学者Rosenberg等利用CT扫描技术建立了枫丹白露砂岩的三维数字岩心，Arns也通过三维CT扫描得到了孔隙度分别为8%、13%、15%和22%四块枫丹白露砂岩的三维数字岩心，扫描分辨率为5.68μm/像素[7-9]。Coenen等通过CT技术成功建立了分辨率小于1μm的三维数字岩心模型[10，11]，可见用CT扫描获取的岩心二维切片及三维立体图像是当前建立数字岩心三维重构可视化最精确最直接的办法，实现针对被测物体内部结构、密度分布和尺寸等信息的定性和定量分析，而不受被测对象内外结构复杂程度影响等特点[2-4]，本文利用Avizo数据可视化软件对CT扫描重建数据分析，可以形象生动地展现所研究问题的实质。

本次研究采用的是天津三英科技公司生产的nanoVoxel 3000 高分辨率CT扫描扫描仪，首先通过NanovoxelScan软件控制CT扫描仪对岩石样本进行三维透视扫描，再采用VoxelStudio Recon三维重构软件还原岩石样本，将重构的数字岩心添加入Dragenfly软件中进行数据的修整、裁切，以及相关信息的提取。基于上述准备工作完毕且保证数据体正确的前提下，利用Avizo计算和提取岩石样本的孔隙度、孔喉大小及连通方式、孔隙占比、骨架结构可视化等信息。实现微米级别孔隙的数据可视化，有助于对地下气液流动机理的认识。

1 实验样品

本文研究实验样品采自于山西省境内华晋吉宁矿区的石炭二叠纪太原组致密砂岩含水层组，岩样埋深450.18 m，平均厚度14.33 m。以浅灰色，灰色细粒砂岩为主，巨厚层状，含岩屑及少量云母，分选中等，次圆-棱角状，黏土胶结，坚硬，下部夹粉砂岩薄层。实验样品的物性分析及核磁共振测试的孔隙度、渗透率结果，见表1。

表1 岩样的孔隙度和渗透率实测结果

从表1中可以看出，两种方法在孔隙度的测试结果上，具有较大的差异(即物性孔隙度1.16%小于核磁孔隙度5.01%)，这是由于物性分析法(气测法)的特殊性导致，其结果仅能反映连通孔隙部分，而核磁共振则能反映出总孔隙(连通孔隙和孤立孔隙)的数据。因此，总孔隙为5.01%，连通孔隙为1.16%，推导出孤立孔隙为3.85%，渗透率范围在0.00183 mD～0.00210 mD之间。

2 三维CT扫描图像中岩心目标提取及重建

在传统的二维CT扫描图像中，选取样本时尽量选择形状规则的，表面平整度较高的样本，其目的是为了降低岩心周边的噪声点，便于后期降噪处理，而三英精密三维CT扫描仪，可以做到在不规则形状的岩心上进行三维CT扫描，其扫描后的数据，随即可运用三维重构软件将边缘噪声点以及硬化校正，实现岩心样本的三维可视化，如图1所示。

三维CT扫描是建立三维数字岩心最常用的方法，但在实际应用过程中，也存在一些问题[5,6]：同一块岩样，由于分辨率的选定不同，会导致建立的孔隙网络模型存在明显差异，若选取的过低，就难以辨识岩样中的微细孔隙连通通道，计算机会自动忽略这部分孔隙通道，进而得到的孔隙度数据会远小于实际的孔隙度；若反之，分辨率过高，扫描出的数据体较大，这对计算机的硬件要求较高[1]。因此本次致密砂岩岩样，根据其实测孔隙度和渗透率数据，如表1所示。选取分辨率为2.98μm/像素，进行三维数字岩心模型的重建工作。在制作岩样及扫描岩样的过程中，应当注意以下几点：

图1 岩心样本及其CT图像分析

图2 三英精密CT扫描仪内部扫描装置工作原理(射线源、样品台及探测器)

(1) X射线源：射线源的位置应当调节适中，在没有触碰到样品的前提下，尽量将射线源靠近样品台，以求最大限度的提升目标扫描区的分辨率，为后续的数字岩心的三维建模提供一个较好的载体。若扫描过程中未能做到这一点，扫描结果极有可能由于分辨率过低而导致重建数据不可用。

(2) 样品台：样品需要制作成多大的体积，是根据研究目标体的致密程度来决定的，若样品过于致密(即孔隙度极低或射线源极不易穿透)的情况下，样品体积需要磨制的尽可能小，与此同时，样品台距X射线源的距离应当尽可能靠近。在固定样品时，应当尽可能地避免岩石样品在扫描过程中出现偏移、滑动等现象的发生，因为偏移和滑动会导致三维重建数据体出现图层间的不对接，最终成像模糊，孔隙反映呈“块状”分布，失去了三维成像可视化的意义。

(3) 探测器：探测器的位置关系到样品在计算机屏幕中的成像，当X射线源距离样品台足够近时，需将探测器的位置进行左右调整，直至样品整体出现在电脑屏幕的视域正中央，且样品整体占视域面积的三分之二时，为最佳状态。

(4) 电压和电流：在调节电压和电流的过程中，应当时刻关注样品区及非样品区的灰度值，切忌过高和过低。一般而言，样品区与无样品区灰度的比值在30%～50%为最佳，如图3、4所示。

图3 电脑控制CT扫描仪参数界面

图4 电脑控制CT扫描仪参数界面

(5) 重建：结束扫描后，将数据体完整地加载到重建软件中，在加载过程中，需要准确无误地记录下数据体X轴、Y轴、Z轴方向以及数据体分辨率的大小等相关信息，该信息也将在孔隙度分析软件Avizo中应用。

(6) 图像裁切：在重建好的数据体中进行裁切时，应确保非样品区完全裁切，保留下的样品成像部分需完整且符合研究要求。

(7) 阈值范围：在对裁切完成的岩心样本进行阈值分割时，切记时刻关注岩样每一层的孔隙度是否选中，确保最大限度的降低人为干预的影响。

3 基于三维CT图像的孔隙提取方法

通过上述操作，得到岩心的重建模型后，将重建数据导入Avizo图像处理软件中，根据所要完成或达到的效果，点选相应的命令。主要操作步骤如下：

(1) 提取总孔隙；

① 将原数据预处理，包括裁切和过滤(如图5、6所示)；

② 再进行阈值分割，即将孔隙和骨架提取为两种颜色以示区分(如图7所示)；

③ 导出的孔隙部分，计算孔隙体积占总体积的百分比；

(2) 提取连通孔隙，展示连通孔隙与孤立孔隙的位置关系(如图8所示)。

提取连通孔隙及孤立孔隙，是为了后期可以更好地将两类孔隙分类量化分析。提取连通孔隙，可以得到连通孔隙的球棒模型，即最大球算法模型如图9所示。其意义在于可更为直观的反映两个孔隙间的相互位置关系及连通情况，球状图形为孔隙，管状图形为连通通道(在储层的研究中，该部分被称为喉道；在水害及瓦斯防治研究中，该部分被称为渗流通道)，若球状图形连接的管状图形越多，则说明其配位数越高，即连通程度越强，渗透率就会越大；而提取出的孤立孔隙部分，可以在进行孔隙量化分析后，找出孤立孔隙中直径较大的孔隙，该部分孔隙，对孔隙研究者而言，其意义和价值并不高，而对裂隙研究者而言，应当高度重视，因为后续在岩样做压力—渗透试验时，该部分孤立孔隙会先于其他小孔隙发生破碎而产生裂隙或裂理，对岩样渗透率的贡献也远大于岩样本身连通孔隙的渗透数值。

由于微米CT三维重建模型数据体较大，运用Avizo进行最大球算法(即球棒模型)，会使计算机运算变得缓慢甚至出现无法应用的现象，因此可采用SypiCore软件对建模数据体进行快速的最大球算法计算，其结果如下图10所示。连通孔隙呈现出局部层状分布，上下层之间处于弱连通或不连通状态，且没有两端贯通孔的出现，判断整体渗透率较弱。

图5 三维CT扫描岩心重建模型

图6 图像裁切及过滤后的骨架结构

图7 二维图像阈值分割

(3) 提取孤立孔隙；

① 利用总孔隙及连通孔隙，导出孤立孔隙部分，并对孤立孔隙进行量化分析，将孤立孔隙按照不同孔径大小划分出来。

② 筛选目标孔径的孤立孔隙部分，筛选目标为直径大于30μm的孤立孔隙。

图8 提取出的连通孔隙及连通孔隙-孤立孔隙的位置关系

图9 连通孔隙及最大球算法(球棒模型)

图10 整个微米CT重建的球棒模型

在提取出来的孤立孔隙中，找到大于30μm的孤立孔隙分布规律，该部分是对渗透率无贡献度或贡献程度极弱的孔隙，除非有外力致使岩样中心因破碎而产生裂隙，来增加连通程度的可能性，否则，它们是不参与导水、导气，甚至是储存油气任务的，它们是闭合的孔隙。

③ 将筛选出的孔径范围求和，导出目标孔径的孔隙度。

(4) 分析总孔隙、孤立孔隙及连通孔隙的孔隙度规律。

4 孔隙度计算结果

1) 孔隙度占比情况

基于上述模型运算操作过程，可以很快得出各类孔隙体积占比情况，见表2。模型切块大小为300×300×300个像素点，每1格代表1个像素点，即体像素为27000000，表2中的孔隙体积为相应单元体像素相加。

表2 模型数据与实测数据对比表

从表2中的信息不难看出：

① 总孔隙体积占切样体积比约为5.34%，与核磁共振的实测孔隙度5.01%较吻合，连通孔隙体积占切样体积比约为0.94%，与物性分析的实测孔隙度1.16%较吻合，说明模型预测效果较好。

② 总孔隙体积中，孤立孔隙体积处于主导地位，约占总孔隙体积的82.4%，模型计算的渗透率0.00197 mD与实际测量的渗透率0.00183 mD～0.00210 mD吻合度较高。

③ 在孤立孔隙中，直径大于30 μm的孔隙体积占的份额非常大，约为2.8%，根据孤立孔隙度约为4.4%，推算出直径小于30 μm的孔隙体积约占孤立孔隙的1.6%，因此直接大于30 μm的孤立孔隙是裂隙研究者们应当着重注意的部分，在给岩样施加压力时，这部分孔隙是产生裂隙的首要位置。

2) 孔隙度计算与核磁共振测试的数据对比

为验证整个数据体模型的有效性及准确性，将模型运算结果同核磁共振测试结果进行比对，如图11所示。

孔隙度计算结果与孔隙度测试结果大体上是接近的。由于X-微米CT扫描仪分辨率的影响，造成模型计算结果中缺失最小的孔分区间，而离心核磁共振是从最小的孔分区间开始记录的，在这一方面，对模型计算结果是一种补充和完善。

5 结论

(1) 图像处理软件Avizo可快速地将CT扫描重建后的岩心进行数字化操作，打破二维图像计算的局限性，使孔隙的大小、孔喉大小及连通程度、孔隙和基质骨架结构的占比等数据信息能够定量化和可视化，实现有针对性地围绕目标区域展开研究。如果我们将连通孔隙单独作为目标区域，它将帮助我们快速地将目标孔隙从模型整体中分离出来，排除其他区域对目标区域的影响。该技术不但缩小了研究区域，还提高了研究精度。

(2) 三维 CT扫描具有无损化、便捷化、快速化以及可视化等优势。该技术起初应用在生命科学领域，在岩心数字化的应用中，还处在起步阶段。本文的研究对象是吉宁矿石炭二叠纪太原组致密砂岩，模型计算结果同岩样实测数据较为吻合，数字岩心真实地反映了孔隙的结构表征。这将会拓宽煤系地层砂岩在微观层面上研究的深度与广度。

图11 孔隙度计算及测试结果对比图