煤巖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.