煤岩CT图像的孔隙度和比表面积测量方法

金智敏 周宏伟 薛东杰

摘 要：为实现基于CT图像的煤岩孔隙度和比表面积快速测量，提出了一种基于MATLAB图像处理的测量方法。以平煤十二矿深部采集的煤岩为研究对象，应用Nano CT技术对样品进行了扫描，并分别借助Avizo和MATLAB软件计算孔隙度和比表面積，研究了Avizo三维建模和MATLAB图像处理2种方法的优缺点。发现2种方法计算出的结果误差仅在10%左右，且以复合Cotes公式的误差估计最小，孔隙度和比表面积的误差分别为9.392%和8.247%.结果表明：传统的基于Avizo三维建模的测量方法精确度较高，同时还能获取样品的三维孔隙网络模型和孔隙的形状、走向及孔径分布，但计算过程复杂，对于复数个试样必须进行多次计算，且对CT图像有较高的质量要求。基于MATLAB图像处理的测量方法则通过简单的数值积分实现了快速计算煤岩的孔隙度和比表面积，避免了三维建模，极大地降低了计算复杂度，且能够一次性处理复数个试样，较好地满足了实际工程需要。且对于CT图像具有更高的容错率，每次计算仅需其中数张扫描质量效果较好的CT图像即可。

关键词：煤岩;CT图像;孔隙度;比表面积;MATLAB;Avizo;

中图分类号：TP 391

文献标志码：A

文章编号：1672-9315（2020）01-0133-08

DOI：10.13800/j.cnki.xakjdxxb.2020.0118开放科学（资源服务）标识码（OSID）：

Measurement of porosity and specific surface

area of coal rock using CT images

JIN Zhi-min，ZHOU Hong-wei，XUE Dong-jie

（School of Mechanics and Civil Engineering，China University of Mining and Technology（Beijing），Beijing 100083，China）

Abstract：In order to realize the rapid measurement of porosity and specific surface area of coal using CT images，a measurement method based on MATLAB image processing is proposed.Taking the coal rock obtained from the deep part of Pingdingshan Coal Mine Group 12 as the research object，the sample was scanned by Nano-CT，and the porosity and specific surface area were calculated by Avizo and MATLAB software respectively.The advantages and disadvantages of Avizo 3D modeling and MATLAB image processing were studied.It is pointed that the errors calculated by the two methods are only about 10%，with the error done by the composite Cotes formulathe smallest.The errors of porosity and specific surface area are 9.392% and 8.247%，respectively.The results indicate that the traditional measurement method based on Avizo 3D modeling is more accurate，and the 3D pore network model and the shape，orientation and size distribution of pores can be obtained.However，the calculating process is comparatively complicated，and it is necessary to perform multiple calculations for a plurality of samples，with high quality requirements for CT images.The measurement method based on MATLAB image processing proposed in this paper realizes the rapid calculation of porosity and specific surface area of coal by simple numerical integration，which avoids 3D modeling，reduces the computational complexity greatly，and processes multiples at a time.The sample satisfies the engineering needs well.Moreover，for CT images，there is a higher fault tolerance rate，and only a few CT images with better scanning quality effects are needed for each calculation.Key words：coal rock;CT image;porosity;specific surface area;MATLAB;Avizo

0 引 言

煤岩是典型的多孔介质，力学上常表现出非均匀性、非连续性和各向异性。其内部结构可视为由固体基质和孔（裂）隙2部分组成，其中孔（裂）隙的微观结构特征表现为孔隙与喉道的几何形状、体积、分布及连通关系，与煤岩的力学性质、渗流能力密切相關[1]。在煤矿瓦斯治理等实际工程问题中，孔隙度和渗透率等参数都具有极其重要的意义。Kozeny-Carman（KC）方程是多孔介质渗流领域著名的半经验公式之一，即

式中 K为渗透率，μm2;为多孔介质的孔隙度;C为KC常数（一般取5）[2];S为比表面积，μm -1.由式（1）可知，多孔介质的渗透率可由其孔隙度和比表面积估算出。

孔隙度和比表面积有多种测定方法，常用的有低压气体吸附法（LPGA）[3-4]、压汞法（MIP）[5]、核磁共振技术（NMR）[5]、聚焦离子束-扫描电镜技术（FIB-SEM）[6]、计算机断层扫描技术（CT）[6-7]、扫描电镜技术（SEM）[8]等。其中计算机断层扫描技术（Computed Tomography，CT）由于具备无损检测物体内部结构、组分及损伤情况的能力，近年来广泛应用于煤炭行业。目前有许多学者在这方面做了大量工作。Zhang等基于CT扫描图像利用Avizo软件对某无烟煤进行了三维重构并研究其孔隙度和渗透率的演化[9]。

Mayo等利用CT技术研究了不同煤样中Xe，Kr和CO2气体吸附量随时间的变化，并对气体扩散曲线和扩散模型进行了分析[10]。王刚等通过CT技术和三维重构技术建立了6种煤样的真实模型并进行模拟试验，对非达西渗流中各参数对渗流的影响进行了研究[11]。李伟等结合显微CT技术和多孔介质逾渗理论，对煤中孔隙进行了三维可视化表征，定量研究了不同煤体结构煤的孔隙连通性和渗透能力的差异[12]。

目前对于实验室尺度的煤样，主要采用显微CT技术表征其孔隙-裂隙结构。煤岩中瓦斯气体吸附、扩散和缓慢层流渗透的主要空间是微米级和纳米级孔隙，然而显微CT扫描精度一般最高只有几十微米，难以满足对这些孔隙的测量。Nano CT精度较显微CT高得多，但扫描样本尺度只有微米级，使用常规方法难以快速测定其孔隙度和比表面积。鉴于此，提出了一种基于MATLAB数字图像处理技术的孔隙度和比表面积测量方法。该方法基于少量的二维CT图像信息，通过简单的数值积分快速计算煤岩的孔隙度和比表面积，无须三维重构，为定量研究煤岩的孔隙结构提供了新的思路。并与基于Avizo三维建模的传统方法对比，结果表明，借助MATLAB处理二维CT图像来测量煤岩的孔隙度和比表面积，能较好地满足工程需要。

1

基于Avizo三维建模的孔隙度及比表面积测定

试验煤岩取自平煤十二矿深部，埋深约1 000 m.煤岩作为一种低渗岩体，其孔隙尺度一般在1 μm左右。试验前已用ACTIS300-320/225工业CT检测系统对煤岩进行成像，分辨率达30 μm时仍难见孔隙，因此改用Nano CT扫描煤岩的孔隙分布。

Nano CT是在Micro CT基础上发展出来的精度更高的计算机断层扫描技术。其原理与Micro CT基本一致，均是根据样品中不同部分对X射线的吸收和透过率差别，测量X射线在不同角度穿过样品时的衰减，并采用滤波反投影法（FBP）[13]、有序子集-最大期望法（OSEM）[14-15]、代数重建法（ART）[16]等算法进行图像重建，将光信号转变成数字信号，从而获得样品的断面CT图像[17]。CT图像的分辨率一般取决于射线等效束宽BW，即

式中 BW为射线等效束宽，μm;d为探测器孔径，μm;a为射线源焦点尺寸，μm，M=L/λ;M为放大倍数;L为射线源到探测器距离，m;λ为射线源到扫描中心距离[18]，m.射线等效束宽越小，CT图像的分辨率越高。

CT图像的锐利度又与半影H有关，即

式中 H为半影，μm.如图1所示，CT图像由本影和半影组成，其中本影指投影时没有任何光线到达的区域，半影指只有部分光线到达的区域。半影越小，CT图像的锐利度越高。由式（2）和式（3）可知，CT图像的分辨率受到射线源焦点尺寸、CCD分辨率（即探测器孔径d）及放大倍数的约束，而在射线源焦点尺寸和CCD分辨率一定的情况下，任意增大放大倍数以期提高CT图像的分辨率，可能会导致严重的半影模糊。因此，提高CT图像分辨率的主要途径就是缩小射线源焦点和提高CCD分辨率。

试验仪器选用Xradia Ultra-XRM L200立体显微镜，该CT检测系统拥有65和15 μm 2种视场模式。基本工作原理如图2所示。同步辐射源发出的X射线经过单色器分离出单色X光后，再通过聚焦波带片聚焦，并用小孔过滤掉零级光和高级衍射光，这样在理论上就产生了一个焦点无穷小的射线源。通过小孔的X光照射到焦点附近的样本上，再通过物镜波带片和显微物镜放大成像，同时使用相衬环提高透明物体的清晰度，以提高CCD分辨率。

试验选用扫描直径为65 μm的大视场模式，测试电压8 kV.最终所测得的图像大小为1 024×1 024像素，每个像素点的边长约为0.063 μm.试验最终测得1014张二维CT层析图像，依次编号分别为1～1 014.为了使结果更加精确，截取了样品内部较清晰处的内接正方体以保证样品的代表性。选用正方体不仅便于统计孔隙度与孔隙结构参数，而且立体展示效果更佳。最终所截取的正方体边长为601像素，即只选取了其中601张601×601像素的图像进行分析。如图3所示，以编号702的图像为例，图中黄色边框内即为分析范围。

试验所得的图片均保存为tiff格式。所以在将图片导入三维图像数据分析软件Avizo前，需将得到的tiff图片进行格式变化。载入数据后，还需对二维CT层析图像进行图像处理。一个标准的岩心图像处理流程包括：灰度图像滤波处理、灰度图像的二值化、基于分水岭算法的图像分割、结果的分析以及渗流方向上的孔隙连通性判断等等[20]。在Avizo软件中，滤波降噪常用的命令模块包括Median filter，Non-local Means filter，Anisotropic Diffusion filter等。选用各向异性扩散滤波函数（Anisotropic Diffusion filter）对图像进行处理，并添加Interactive Thresholding命令模块进行图像分割，其中阈值参数为59.仍以编号702的图像为例，处理前后的图像如图4所示。图4（a）即为原始切片图像，图4（b）则为图4（a）经过Avizo滤波降噪、阈值分割和二值化等图像处理后得到的图像。其中蓝色部分为孔隙或裂缝，黑色部分为基质矿物。

通过Avizo重构的煤岩的三维数字模型如图5所示，该模型选用编号200～800的二维CT图像沿y轴排列并重构。煤岩的孔隙结构模型则如图6所示。Avizo软件主要基于分水岭算法识别孔隙边界并对各孔隙进行标记，其基本原理为将CT图像视为具有拓扑结构的地形图，用各点的灰度值表示该点的海拔高度，则每一个局部极小灰度值及其影响范围构成了一个盆地，而各盆地的交汇处则形成分水岭。图5和图6的三维重构模型中总计有6 360个孔隙。在Avizo中添加Label Analysis命令对各孔隙进行单独分析，通过统计各孔隙中的体素个数及边界上的像素个数来表征孔隙的体积及表面积，并按孔隙大小（即等效直径）分别对孔隙数量、表面积、体积的分布进行统计，结果如图7所示。最终计算得到孔隙的体积为3 600.286 μm 3，表面积为26 334.324 μm 2，所以该煤岩的孔隙度为6.633%，比表面积为0.485 μm -1.

2 基于MATLAB图像处理的孔隙度及比表面积测定

由二维CT图像，可以测出该图像所对应截面上的孔隙周长与面积。从图4可以看出，将二维CT图像沿y轴堆积即可得到煤岩的三维数字模型，各图像上的孔隙周长与孔隙面积均可视为轴向坐标y的函数。设试样整体体积为V，孔隙体积为Vp，孔隙内表面积为δ，垂直于截面方向的长度为l.再设截面上孔隙面积为Sp，孔隙周长为Pp.由此可得孔隙率为

式中 V，Vp分别为试样整体体积和孔隙体积，μm3;δ为孔隙内表面积，μm2;l为垂直于截面方向的长度，μm;Sp为再设截面上孔隙面积，μm2;Pp为孔隙周长，μm.

对于式（4）和式（5）中的积分部分，可以利用机械求积方法近似求解，将离散的CT图像数据连续化，即

式中 yk为求积节点;Ak为节点yk的权，仅与节点yk的选择相关。因此，只需数张二维CT层析图像，测出它们的孔隙周长与孔隙面积，即可通过数值积分的方法计算出试样整体的孔隙率与比表面积。

为了方便与前一节的对比，同样从编号200～800共601张图像中选取了编号为200，275，350，425，500，575，650，725，800的二维CT图像。由于原始图像位深度为32位，MATLAB中的imcrop函数无法识别，故需将其导入画图、ACDSee等图像编辑工具中转存为位深度24位的bmp格式后再导入MATLAB。图像在该过程中仅丢失了代表透明度的Alpha通道信息，对原始数据并无影响。再参照图3剪裁其中边长为601像素的正方形以便于分析，其中该正方形左上角像素在原图中的坐标为（240，200）。

CT图像中的噪声普遍存在。为了有效地抑制噪声的干扰，必须对原始图像进行滤波以突出图像中的有效信息。滤波主要可分为3类：线性滤波[21]、中值滤波[22-23]和自适应滤波[24-25]。由于消除噪音的同时还需保持图像细节清晰，且图像中需统计的细节如点、线、尖顶部等较多，因此选择自适应滤波对原始图像进行处理。自适应滤波是一种新型的信号处理方法，它是一种基于最小均方误差准则的最优估计，对高斯白噪声的去除效果尤其明显。它既能一定程度上克服线性滤波后图形细节模糊的问题，也能避免中值滤波导致的图像细节缺失对统计结果的影响。

对原始图像进行滤波后，可以获得更平滑的灰度直方图信息。以编号200的CT图像为例，其自适应滤波前后的灰度直方图变化如图8所示。从图8（a）可以看出，滤波前该灰度直方图有3个波峰，但其中2个分别在左端（暗部）和右端（亮部）产生了溢出。而滤波后该灰度直方图两端溢出消失，波形虽有一定的偏移，但波峰和波谷的偏差不大，如图8（b）所示。

当基于双峰直方图的谷底阈值法推广到三峰直方图时，可以得到：三值化时若灰度直方图具有3峰，3峰间的2个波谷则为2个阈值点[26]。因此选取了图8（b）中2个波谷作为阈值点。这样对CT图像进行三值化处理，不仅可以区分孔（裂）隙和固体基质，还能区分固体基质中的煤基质和煤杂质。第1個波谷为[30，80]区间（图8（b）中红色区间），在此区间中图像可近似为水平线，其误差不超过1%;第2个波谷则在225附近。因此，可将2个阈值点分别近似取为55和225，即灰度值在[0，55]的像素点代表孔（裂）隙，灰度值在[56，225]的像素点代表煤基质，灰度值在[226，255]的像素点代表煤杂质。

图9（a）和（b）分别给出了编号200的CT原始切片图像及其三值化结果。图9（b）中黑色区域为孔（裂）隙，灰色区域为煤基质，白色区域为煤杂质。其三值化后的灰度直方图则如图9（c）所示。但由于文中主要测量的是CT图像的孔隙周长与孔隙面积，因此只需对原始图像进行二值化处理即可。按照灰度特性将图像划分为背景和目标2部分，其中将孔（裂）隙视为目标，煤基质和煤杂质等固体基质视为背景。阈值分割点的灰度值取55，二值化后的图像则如图10所示。其中黑色像素点代表孔（裂）隙，白色像素点代表固体基质。图10黑白颜色对调后所有像素的像素值求和的结果即为编号200图像的孔隙面积，单位为像素。

利用for循环语句对图10中所有黑色像素点的4邻域像素值进行求和并映射到一个新的601×601矩阵中，该矩阵中的任一元素表示的是所对应像素点在孔隙边界上的长度，单位为像素边长。该矩阵所有元素之和即为图10中的孔隙周长。编号200，275，350，425，500，575，650，725，800的维CT图像的孔隙周长与孔隙面积见表1，其中每个像素点边长约为0.063 μm.

得到了各CT图像的孔隙周长与孔隙面积，则可通过式（4）和式（5）计算煤样的孔隙度和比表面积。式（6）中常用的机械求积公式有复合梯形公式、复合Simpson公式、复合Cotes公式等，分别记为

最终各数值积分方法计算出的孔隙度、比表面积见表2.

若以Avizo三维建模测定的孔隙度及比表面积为精确值，复合梯形公式、复合Simpson公式、复合Cotes公式计算出的孔隙度、比表面积误差则见表3.从表3可以看出，各数值积分方法计算出的孔隙度和比表面积误差均是复化Cotes公式最小，复化Simpson公式次之，复化梯形公式最大。各误差均在10%左右，表明提出的测量方法能较好地估计煤样的孔隙率和比表面积。

3 结 论

1）传统的基于Avizo三维建模的测量方法精确度较高，同时还能获取样品的三维孔隙网络模型和孔隙的形状、走向及孔径分布，但计算过程复杂，对于复数个试样必须进行多次计算，且对CT图像有较高的质量要求。

2）提出的基于MATLAB图像处理的测量方法则通过简单的数值积分实现了快速计算煤岩的孔隙度和比表面积，避免了三维建模，极大地降低了计算复杂度，且能够一次性处理复数个试样，更加简便经济，较好地满足了工程需要。且对于CT图像具有更高的容错率，每次计算仅需其中数张扫描质量效果较好的CT图像即可。

参考文献（References）：

[1]Song L，Ning Z F，Duan L.Research on reservoir characteristics of Chang 7 tight oil based on nano-CT[J].Arabian Journal of Geosciences，2018，11（16）：472.

[2]徐 鹏，邱淑霞，姜舟婷，等.各向同性多孔介质中Kozeny-Carman常数的分形分析[J].重庆大学学报，2011，34（4）：78-82.XV Peng，QIU Shu-xia，JIANG Zhou-ting，et al.Fractal analysis of Kozeny-Carman constant in the homogenous porous media[J].Journal of Chongqing University，2011，34（4）：78-82.

[3]Zhang R，Liu S，Bahadur J，et al.Changes in pore structure of coal caused by coal-to-gas bioconversion[J].Scientific Reports，2017，7（1）：3840.

[4]Rodrigues C F，Sousa M J L D.The measurement of coal porosity with different gases[J].International Journal of Coal Geology，2002，48（3）：245-251.

[5]Li X，Kang Y，Haghighi M.Investigation of pore size distributions of coals with different structures by nuclear magnetic resonance（NMR）and mercury intrusion porosimetry（MIP）[J].Measurement，2018，116：122-128.

[6]Liu S，Sang S，Wang G，et al.FIB-SEM and X-ray CT characterization of interconnected pores in high-rank coal formed from regional metamorphism[J].Journal of Petroleum Science and Engineering，2017，148：21-31.

[7]Karacan C O，Okandan E.Adsorption and gas transport in coal microstructure：investigation and evaluation by quantitative X-ray CT imaging[J].Fuel，2001，80（4）：509-520.

[8]宫伟力，李 晨.煤岩结构多尺度各向异性特征的SEM图像分析[J].岩石力学与工程学报，2010，29（4）：2681-2689.GONG Wei-li，LI Chen.Multi-scale and anisotropic characterization of coal structure based on SEM image analysis[J].Chinese Journal of Rock Mechanics & Engineering，2010，29（4）：2681-2689.

[9]Zhang G，Ranjith P G，Perera M S A，et al.Characterization of coal porosity and permeability evolution by demineralisation using image processing techniques：A micro-computed tomography study[J].Journal of Natural Gas Science & Engineering，2018，56（8）：384-396.

[10]Mayo S，Josh M，Kasperczyk D，et al.Dynamic micro-CT study of gas uptake in coal using Xe，Kr and CO2[J].Fuel，2018，212：140-150.

[11]王 剛，杨鑫祥，张孝强，等.基于CT三维重建的煤层气非达西渗流数值模拟[J].煤炭学报，2016，41（4）：931-940.WANG Gang，YANG Xin-xiang，ZHANG Xiao-qiang，et al.Numerical simulation on non-Darcy seepage of CBM by means of 3D reconstruction based on computed tomography[J].Journal of China Coal Society，2016，41（4）：931-940.

[12]李 伟，要惠芳，刘鸿福，等.基于显微CT的不同煤体结构煤三维孔隙精细表征[J].煤炭学报，2014，39（6）：1127-1132.LI Wei，YAO Hui-fang，LIU Hong-fu，et al.Advanced characterization of three-dimensional pores in coals with different coal-body structure by Micro-CT[J].Journal of China Coal Society，2014，39（6）：1127-1132.

[13]Pan X C，Xia D，Zou Y，et al.A unified analysis of FBP-based algorithms in helical cone-beam and circular cone-and fan-beam scans[J].Physics and Medicine &Biology，2004，49（18）：4349-4369.[14]Liu X，Comtat C，Michel C，et al.Comparison of 3-D reconstruction with 3D-OSEM and with FORE+OSEM for PET[J].IEEE Transactions on Medical Imaging，2001，20（8）：804-814.

[15]David S，Burion S，Tepe A，et al.Experimental validation of an OSEM-type iterative reconstruction algorithm for inverse geometry computed tomography[C]//Medical Imaging 2012：Physics of Medical Imaging. International Society for Optics and Photonics，2012，8313（6）：125.[16]Herman G T，Meyer L B.Algebraic reconstruction techniques can be made computationally efficient[positron emission tomography application][J].IEEE Transactions on Medical Imaging，1993，12（3）：600.

[17]Peyrin F，Dong P，Pacureanu A，et al.Micro-and nano-CT for the study of bone ultrastructure[J].Current Osteoporosis Reports，2014，12（4）：465-474.

[18]劉郁纪.X射线工业CT物理设计及图像重建[D].兰州：兰州大学，2010.LIU Yu-ji.The physical design and image reconstruction of X-ray CT[D].Lanzhou：Lanzhou University，2010.[19]李 光，罗守华，顾 宁.Nano CT成像进展[J].科学通报，2013，58（7）：501-509.LI Guang，LUO Shou-hua，GU Ning.Research progress of Nano CT imaging[J].Chinese Science Bulletin，2013，58（7）：501-509.

[20]陶 鹏.基于数字岩心的低渗储层微观渗流机理研究[D].成都：西南石油大学，2017.TAO Peng.Mechanism of micro seepage in low-permeability reservoirs based on digital core[D].Chengdu：Southwest Petroleum University，2017.

[21]Flamant J，Chainais P，Bihan N L.A complete framework for linear filtering of bivariate signals[J].IEEE Transactions on Signal Processing，2018，66（17）：4541-4552.

[22]孙宏琦，施维颖，巨永锋.利用中值滤波进行图像处理[J].长安大学学报（自然科学版），2003，23（2）：104-106.SUN Hong-qi，SHI Wei-ying，JU Yong-feng.Image processing with medium value filter[J].Journal of Changan University（Natural Science Edition），2003，23（2）：104-106.

[23]Huang T，Yang G，Tang G.A fast two-dimensional median filtering algorithm[J].IEEE Transactions on Acoustics，Speech and Signal Processing，1979，27（1）：13-18.

[24]张旭明，徐滨士，董世运.用于图像处理的自适应中值滤波[J].计算机辅助设计与图形学学报，2005，17（2）：295-299.ZHANG Xu-ming，XV Bin-shi，DONG Shi-yun.Adaptive median filtering for image processing[J].Journal of Computer-Aided Design & Computer Graphics，2005，17（2）：295-299.

[25]Althahab，Jumaah A Q.A new robust adaptive algorithm based adaptive filtering for noise cancellation[J].Analog Integrated Circuits and Signal Processing，2018，94（2）：217-231.

[26]钟江城，周宏伟，任伟光，等.基于CT图像灰度分布的含杂质煤体三值化方法[J].力学与实践，2018，40（2）：140-147.ZHONG Jiang-cheng，ZHOU Hong-wei，REN Wei-guang，et al.A three-value-segmentation method of coal containing inclusion based on gray distribution of computed tomography image[J].Mechanics in Engineering，2018，40（2）：140-147.