1. 引言
多孔岩石的孔隙结构发育情况对岩石的物理力学性质、工程地质特征和渗流特性有着强烈的影响 [1] [2] [3] ，一直以来都有很多学者致力于测量和表征多孔岩石的孔隙发育情况。孔隙率作为一个可以宏观表征多孔岩石孔隙体积占比的指标，不但被广泛认为可以用来体现多孔岩石孔隙的发育程度，还被认为会直接影响到它的力学性能 [4] [5] 。然而对于内部不可见、孔隙结构复杂、孤立孔隙广泛存在的多孔岩石，孔隙率指标并不易获得。
目前针对多孔岩石孔隙率的测量方法有很多种，包括传统的质量体积法 [6] 和饱水法 [7] 。前者通过岩石材料的密度和岩石的干质量来计算出岩石内部孔隙的体积，而后者则通过用液体水饱和多孔岩石并计算饱和用水的体积来获得孔隙部分的体积。随着技术手段的进步，一些新的测量技术陆续出现。比如侵入法，包括低温氮气吸附法 [8] 和水银压汞法 [9] ，两者都是通过向多孔岩石内部注入外部流体，通过判断流体的表现和侵入量来计算多孔岩石的孔隙率。还有基于CT扫描技术、核磁共振技术和扫描电镜技术的图像法，通过获得多孔岩石内部的图像像素信息计算代表孔隙部分像素/体素的数量从而获得孔隙率数值 [10] [11] [12] 。多孔类岩土内部的声波波速和岩土体的孔隙率也存在着一定的关系，可以用来判断甚至计算其孔隙率 [13] 。此外，一些多孔岩石孔隙率的计算模型也被许多学者提出 [14] ，这些模型对于某些类型的多孔岩石的孔隙率计算具有较高的准确度。
珊瑚礁灰岩是一种生物成因沉积灰岩，化学成分以CaCO3为主，矿物组分主要为文石和方解石 [15] ，在我国广泛分布于东南沿海地区及南海岛礁上 [16] 。它是一种由珊瑚虫遗骸及其分泌物经长期地质沉积而形成的独特多孔岩石 [17] ，这一岩石也继承了珊瑚虫生物骨骼发育、孔隙结构复杂的特点。为了明确珊瑚礁灰岩的孔隙结构特征与其宏观物理力学性质间的关系，需要对其孔隙率进行准确测量。现有多种多孔岩石孔隙率测量方法都有各自的适用范围，尚未见有关不同方法在珊瑚礁灰岩孔隙率测量中适用性的探讨。
本研究以海南地区某类珊瑚礁灰岩试样为研究对象，讨论三种常用的多孔岩石孔隙率测量方法，饱水法、水银压汞法和图像法对珊瑚礁灰岩这种孔隙结构发育而复杂多孔岩石孔隙率测量的适用性。
2. 珊瑚礁灰岩试样
所分析试样为从海南省某珊瑚礁灰岩地层工程现场获得的一类珊瑚礁灰岩，经过加工得到高度为10 cm，直径为5 cm的标准圆柱试样，如图1所示。试样整体呈灰白色，质地较脆，用手可以轻易从其表面将岩屑剥落下来。从外观来看，这一块珊瑚礁灰岩孔隙结构并不发育，但近距离观测可在其表面观察到大量复杂的孔隙特征，如图1右侧所示。因其孔隙结构较为致密，将其命名为CR-D (coral reef limestone with dense pore structure)。
Figure 1. A standard sample of CR-D and its partially enlarged picture
图1. 珊瑚礁灰岩标准试样及其局部放大图
经过烘干处理后，测量CR-D的干质量为244.9 g，计算干重度为12.22 kN/m3，质量较轻，孔隙结构发育，将其置入水中可以观察到大量的气泡从试样内部冒出。准确对其孔隙率进行测量是研究其孔隙结构特征的重要基础。
3. 三种孔隙率测量方法简介
对于珊瑚礁灰岩试样CR-D这类孔隙结构复杂发育的多孔岩石，并不是所有孔隙率测量方法都适用。首先质量体积法并不适合，因为应用质量体积法的前提是能够明确所测多孔岩石的矿物组成成分及每种成分的质量占比。但对于成因复杂、孔隙结构发育的珊瑚礁灰岩，很难准确测量试样每个部位的矿物成分以获取每种矿物的质量。而低温氮气吸附法则主要用于测量孔隙发育程度较差的岩石，如页岩的孔隙结构信息 [9] ，在较大的孔隙中，其或许可以测量出岩石内部结构的比表面积，但无法测量孔隙率。而一般的孔隙率计算模型因其特殊的适用范围也无法应用到珊瑚礁灰岩的孔隙率测量上。
因此，本文选择饱水法、水银压汞法(MIP)和基于CT扫描图像的阈值分割法三种通用性较强的方法来对珊瑚礁灰岩试样的孔隙率进行测量，其基本原理及计算结果简述如下。
3.1. 饱水法
饱水法是一种较为常规的孔隙率测量方法，它适用于大部分由不溶于水材料组成的多孔介质的孔隙率测量。当然，对于不同孔隙发育程度的多孔介质，测量的精度也不同，孔隙较不发育不连通的多孔介质将较难饱和，测量结果也将不准。该方法通过室内饱水试验，获得岩样的孔隙率。首先对试样进行烘干，称量其干质量md。然后将试样置入纯净水中进行饱和。考虑珊瑚礁灰岩试样孔隙结构的复杂性，整个饱和过程持续48 h以保证试样完全被水饱和，并测量饱和后的试样质量ms。可以通过下式计算试样的孔隙率p1 [7] 。
(1)
其中ρw为水的密度，在数值上等于1 g/cm3。V为CR-D试样的表观体积，为196.35 cm3。md和ms经测量分别为313.9 g和244.9 g。最后计算得孔隙率p1为0.3514。
3.2. 水银压汞法(MIP)
水银压汞试验MIP (Mercury Intrusion Porosimetry)基于非浸润液体只有在施加压力的情况下才会进入固体孔隙结构这一原理 [18] ，在不断的增压下，可被液体侵入孔隙的孔径与所用压力呈函数关系式。水银压汞法也是一种适用范围较广的多孔介质孔隙结构测量方法，但容易造成试样的破坏。对于圆柱形孔隙，孔径和压力之间的关系可按Washburn公式进行计算：
(2)
其中，pH为所施加的压力，对于MIP试验，即为水银的压力；γ为水银的表面张力，试验中取为0.48 N/m；θ为水银和固体材料之间的接触角，试验中取水银和珊瑚礁灰岩材料之间的接触角为140˚；r为圆柱形孔隙的半径。
上式主要用来测量被测固体的孔径分布规律，而并不能直接获得孔隙率。水银在低压的情况下无法充满试样内部的孔隙结构，而在高压的情况下则可以充满，如图2所示。故通过不断地提高水银侵入的压力，直到当压力连续提高多次后侵入试样内部的水银量维持不变时，可认定为水银已充满整个孔隙结构，即可通过下式计算被测物体的孔隙率p2。
(3)
式中，VP为被测物体孔隙部分的体积，VPH为在高压条件下水银所侵入的体积，Vt为被测物体的总体积。
Figure 2. Schematic diagram of porosity measurement by MIP
图2. 水银压汞法测量孔隙率原理图
由于水银压汞试验无法在CR-D这样大尺寸的试样上进行，因此从CR-D的不同部位选取三个边长为8 mm的微小的立方体试样，如图3所示，即Vt等于1.536 cm3。经过烘干处理后，三个微小试样进行水银压汞试验。
试验所用仪器为美国麦克仪器公司生产的AutoPoreⅣ9500全自动压汞仪，其最大进汞压力可达413.4 MPa，孔径测量范围为0.003~360 μm。在最高进汞压力下，被测试样中所有孔径大于0.003 μm的连通孔隙都将被水银填满。对于图3中的三块珊瑚礁灰岩试样，考虑其表观孔隙率较大，因此当其中0.003 μm孔径的孔隙都被水银填满后，完全可以认为试样的孔隙结构中已经充满了水银，则可通过这部分水银的体积计算试样的孔隙率。水银压汞试验结束后，计算这三块珊瑚礁灰岩小试样的平均孔隙率为0.4763。
Figure 3. The samples used for MIP measurement, which were taken from the same rock matrix as CR-D
图3. 水银压汞法所用试样(与CR-D取自同一块岩样)
3.3. 基于CT扫描图像的阈值分割方法
CT (Computed Tomogaphy)扫描技术作为一种可以在不破坏被扫描物体的前提下，准确探测到物体内部结构信息的技术，已经被广泛应用于包括医学、物理、土木工程等多个领域。它通过不断地向被扫物体发射X射线而获得物体内部的信息，当X射线穿过一种物质时，部分射线会由于光电效应和散射等而呈现出特定的射线衰减效应 [19] ，通过综合分析不同X射线发射角度时射线的衰减程度即可以判断出被测物体内部不同部位的物质信息，达到区分它们的目的。X射线强度衰减公式如式4所示。
(4)
其中I0是原始的X射线强度，I是透过物体的X射线强度，μ是X射线衰减系数，h是物体厚度，对于不同的物质，衰减系数是不同的。
利用CT扫描图像计算多孔岩石的孔隙率的基本思路如下：首先，对图片进行二值化处理，即要区分开岩石内部的基质部分和孔隙部分；然后，基于岩石的诸多二值化切片图像将整个试样在计算机中进行三维重构；最后，通过三维可视化软件或者体素计算算法将岩石中孔隙部分的占比计算出来，即孔隙率。CT扫描与三维重构过程如图4所示。
Figure 4. CT scanning procedure and three-dimensional reconstruction
图4. CT扫描与三维重构流程
本文CT扫描在德国依科视朗(YXLON)公司的YXLON FF85计算机断层扫描检测系统上利用螺旋扫描法进行的，最后扫描的图像精度达到了18.4174 μm/像素。在整个扫描过程中，扫描电压设置为209.6 KV，扫描电流为110.2 μA。三维重构主要在三维可视化软件Avizo中进行，将珊瑚礁灰岩CT扫描图像导入Avizo中并经过图像降噪、阈值分割等操作，即可对试样进行三维重构从而计算孔隙率。整个CR-D数字岩芯包含有约400亿个体素，这将会经常造成计算的中断。为保证计算稳定，提高计算效率，整个CR-D数字岩芯被沿着试样轴向方向分为十段，每段为一个直径5 cm，高度1 cm的圆盘，每个圆盘的孔隙率将被分别计算最后取平均即为整个CR-D的孔隙率。
在利用基于CT扫描技术的图像法计算多孔岩石的孔隙率时，扫描图像的二值化处理尤为重要。岩石的原始切片图像往往灰度值分布范围广、模糊不清，需要准确判断代表基质像素值和代表孔隙像素值的界限，即阈值。目前的图像阈值分割方法有很多种，常用的方法包括自适应阈值分割方法、最大类间方差法、最大熵法、矩保持法(moment-preserving)和灰度值双峰法等。
其中，灰度值双峰法是考虑某些图片的灰度直方图呈现出两个峰值分别代表着两种物质，则可将双峰间谷底处的灰度值定为阈值。但如图5所示，珊瑚礁灰岩CR-D的灰度值分布并不是“双峰型”，两个不同位置扫描图像的灰度直方图的双峰效应都不明显，谷底处像素值的像素个数仍很多，且两个切片的谷底处所对应的像素值也并不相同，所以这一方法并不适用。将其它的一些全局和局部的阈值分割方法在CR-D的切片图像上进行尝试，发现都呈现出较大程度肉眼可见的阈值过分割和欠分割的情况。本文采用自适应阈值分割方法、最大类间方差法和最大熵法来计算CR-D试样的孔隙率。
Figure 5. Gray histogram of CR-D’s CT scanned images. (a) Slice 25 mm from the bottom of the sample; (b) Slice 75 mm from the bottom of the sample
图5. CT扫描CR-D切片灰度值直方图。(a)距离试样底部25 mm处切片；(b)距离试样底部75 mm处切片
3.3.1. 自适应阈值分割方法
自适应阈值分割方法，又可被称为动态或者局部阈值分割方法，针对图像中具有不同阴影、对比度的位置，图像分割将不设定单一的界限，而是对图像按照坐标进行区域划分，对每一处都选择特定的阈值进行分割 [20] 。
自适应阈值分割方法二值化效果如图6所示，图片选取距CR-D试样底部5 cm处截面的一部分，其中蓝色部代表试样的基质部分。仅在200像素 × 200像素的视窗内即可清晰地观察到珊瑚礁灰岩试样CR-D较为发育且复杂的孔隙结构。根据这种阈值分割方法计算的CR-D孔隙率数值为0.4877。
Figure 6. Binarization effect of the adaptive threshold segmentation method (the position of the picture is 5 cm from the bottom of the sample, where the blue part is the rock matrix)
图6. 自适应阈值分割方法二值化效果(图片位置距试样底部5 cm，其中蓝色部分为基质)
3.3.2. 最大类间方差法
最大类间方差法，又称大津法，由日本学者大津在1979年提出 [21] ，这种方法主要基于图像的灰度值信息，将图像区分为背景和目标两个部分，以寻找一个灰度值使得这两部分的灰度值类间方差最大，这一灰度值即为阈值。这种方法由于计算简单，不受图像亮度和对比度的影响，因此错分概率较小 [22] 。
从图5可以看出，CR-D扫描图像的灰度值分布并不均匀，而是主要分布在30~110之间。因此，为了提高计算效率，保证阈值判断的准确性，最大类间方差法和最大熵法都将基于灰度值的最小值和最大值分别为30和110这一设定进行计算。
最大类间方差法二值化效果如图7所示，其中蓝色部代表试样的基质部分。根据这种阈值分割方法计算的CR-D孔隙率数值为0.5475。
Figure 7. Binarization effect of the OTSU method (the position of the picture is 5 cm from the bottom of the sample, where the blue part is the rock matrix)
图7. 最大类间方差法二值化效果(图片位置距试样底部5 cm，其中蓝色部分为基质)
3.3.3. 最大熵法
最大熵法的核心思想和最大类间方差法相似，但与最大类间方差法通过判断背景与目标之间的灰度值方差不同，最大熵法主要通过信息熵来衡量目标与背景之间的差距以确定阈值 [23] 。
最大熵法二值化效果如图8所示，其中蓝色部代表试样的基质部分。根据这种阈值分割方法计算的CR-D孔隙率数值为0.5553。
Figure 8. Binarization effect of the maximum entropy method (the position of the picture is 5 cm from the bottom of the sample, where the blue part is the rock matrix)
图8. 最大熵法二值化效果(图片位置距试样底部5 cm，其中蓝色部分为基质)
4. 结果与讨论
4.1. 三种方法对孔隙率测量结果影响分析
将上述三种方法测量到的珊瑚礁灰岩试样CR-D孔隙率结果进行归纳如表1所示。可以发现饱水法所测量的CR-D孔隙率数值仅为0.3514，明显低于其它两种方法所获得的孔隙率。而通过水银压汞法获取到的CR-D孔隙率数据则与图像法中利用自适应阈值分割法获取到的孔隙率相近，前者为0.4763，而后者为0.4877，相差不到0.015。
Table 1. Porosity of CR-D from three segmentation methods
表1. 三种方法测量CR-D孔隙率结果
实际上，珊瑚礁灰岩试样的连通性对水银压汞试验所获得孔隙率数据的准确性有着很大的影响，如果三个微小试样的孔隙连通性很好，则这一方法所获得的孔隙率也将比较接近准确值。因此，对整个珊瑚礁灰岩试样CR-D的孔隙结构连通性进行分析，如图9所示，图9右侧为其中的孤立孔隙部分，红色虚线区域内的一个孤立孔隙被放大展示。经计算，CR-D试样中不仅孤立孔隙的体积占比极小，孤立部分孔隙的体积也仅约占总孔隙体积的0.002%，而且这些孤立孔隙的分布并不集中而是较为离散。CR-D试样中心某个边长为500个像素立方体区域内的连通孔隙部分和孤立孔隙部分如图9所示，可以看出孤立孔隙部分所占的体积极小。由此可推断三个与CR-D来自于同一岩样的珊瑚礁灰岩微小试样中的孔隙也极为连通，且不存在体积较大的孤立孔隙。因此，通过水银压汞试验所获得的CR-D及其原始岩样的孔隙率p1是最为准确的。
饱水法测量CR-D孔隙率偏低是因为试样内部微小而复杂的孔隙结构。在大气压条件下，微小的孔隙结构由于毛细力的存在，水无法充满所有孔隙，错综复杂的孔隙结构更是增大了水渗透的难度，从而降低了实际试样的饱和程度，故饱水法孔隙率测量结果偏低。另一方面，在CR-D饱和操作完成后，需要将其从水中拿出秤其质量，在这个过程中会有一部分水由于试样孔隙过于发育而从试样表面流出，造成饱和质量测量值偏小，孔隙体积测量值也由此偏小。故孔隙结构发育的多孔类岩石，常压条件下的饱水法会对孔隙率的测量造成较大的误差。
Figure 9. Connectivity analysis of CR-D’s pore structure (blue parts are connected pores, purple parts are isolated pores)
图9. CR-D试样孔隙结构连通性分析(蓝色部分为连通孔隙，紫色部位为孤立孔隙)
4.2. 三种阈值分割方法对图像法测量结果影响讨论
采用三种不同阈值分割方法利用CT扫描图像测量的CR-D孔隙率数据中，自适应阈值法得到的孔隙率最为接近水银压汞法测得的试样的孔隙率，因此更加准确。与全局阈值分割方法不同，自适应阈值分割法会将图像分为多个区域，对不同的区域进行差异化的阈值分割，这种方法较适用于像素灰度值，即图像特征分布不均匀的图像。这表明对于珊瑚礁灰岩这一类孔隙结构特征极为发育的多孔类岩石，基于自适应阈值分割方法的图像法所测得的孔隙率更为准确。
而基于最大类间方差法所测得的CR-D孔隙率要明显偏高。这意味着在阈值分割过程中，很多岩石部分被划分为孔隙部分。因此，对于珊瑚礁灰岩，通过全局灰度值分布特征获得的阈值并不是适用于图像的每个部位的。最大熵法有着同样的问题。
此外，由于最大类间差法和最大熵法本质上都是通过区分目标灰度值和背景灰度值来确定阈值的，因此孔隙部分和基质部分的灰度值差异越明显，即灰度值分布图越呈现双峰特征，这两个方法的图像分割效果越好。但如图5所示，珊瑚礁灰岩CR-D的灰度值分布图的双峰特征并不十分明显，因此这两种阈值分割方法对于珊瑚礁灰岩的适用性均较差。
5. 结论
1) 珊瑚礁灰岩孔隙结构极度发育且连通性很强，对于本文所采用的珊瑚礁灰岩标准试样，孤立孔隙的体积仅占孔隙总体积的0.002%。
2) 对于珊瑚礁灰岩和其它类孔隙结构发育且连通的多孔类岩石，水银压汞法测量得到的孔隙率较准确，而饱水法由于孔隙中难以充满水且已经赋存在孔隙内部的水易流出而造成孔隙率的测量结果偏小。
3) 当采用图像法测量多孔类岩石孔隙率时，不同的阈值分割方法会对孔隙率计算结果产生很大的影响。珊瑚礁灰岩由于灰度直方分布图的双峰特征不明显，最大类间方差方法和最大熵法并不适用于其CT扫描图像的阈值分割，故二值化结不适用于测量孔隙率。
4) 对于孔隙结构复杂的多孔类岩石，当采用图像法计算其孔隙率时，局部阈值分割方法所获得二值化图像所计算出的孔隙率要比根据全局阈值分割方法计算出来的孔隙率更加准确。
基金项目
国家自然科学基金(52008307)。