基于CT图像灰度分布的含杂质煤体三值化方法<sup>1)</sup>

引用本文

钟江城, 周宏伟, 任伟光, 王子辉, 陈吉. 基于CT图像灰度分布的含杂质煤体三值化方法¹⁾ . 力学与实践, 2018,40(2): 140-147
ZHONG Jiangcheng, ZHOU Hongwei, REN Weiguang, WANG Zihui, CHEN Ji. A THREE-VALUE-SEGMENTATION METHOD OF COAL CONTAINING INCLUSION BASED ON GRAY DISTRIBUTION OF COMPUTED TOMOGRAPHY IMAGE¹⁾ . MECHANICS IN ENGINEERING, 2018,40(2): 140-147. 复制到剪切板

Doi: 10.6052/1000-0879-18-003
Permissions

《力学与实践》编辑部所有

基于CT图像灰度分布的含杂质煤体三值化方法¹⁾

钟江城²⁾, 周宏伟³⁾, 任伟光, 王子辉, 陈吉

中国矿业大学(北京)力学与建筑工程学院力学系,北京 100083;

2) 钟江城,博士生. E-mail: zhongjc1988@126.com;3) 周宏伟,教授,主要从事岩石和渗流力学研究. E-mail: zhw@cumtb.edu.cn

收稿日期: 2018-01-02

资助项目: 1) 国家重点研发计划(2016YFC0600704) 和中国矿业大学(北京) 越崎杰出学者奖励计划(2017A03) 资助项目.

摘要

煤体CT图像阈值分割是其三维真实结构重构的前提. 基于单轴压缩状态下煤体的CT图像,提出了一种基于图像灰度分布的孔隙度计算方法, 结合灰度直方图,利用实测孔隙度反推法对煤体CT图像进行了裂隙--煤基质--煤杂质三值化,同时与最大类间方差法和最大熵法的三值化结果进行了对比分析. 结果表明, 最大熵法对于煤基质和煤杂质不能较好区分; 最大类间方差法对于图像微小区域灰度的突变(裂隙的产生)敏感性低,无法有效区分裂隙与煤基质; 本文改进的灰度直方图法避免了最大类间方差法和最大熵法的劣势,具有较好的图像三值化效果,且煤基质与煤杂质间的阈值基于实验结果得到,具有更高的可信度和准确度.

关键词: CT图像; 三值化; 孔隙度; 最大类间方差法; 最大熵法; 改进的灰度直方图法

中图分类号:TU45,TP391.41 文献标志码:A

A THREE-VALUE-SEGMENTATION METHOD OF COAL CONTAINING INCLUSION BASED ON GRAY DISTRIBUTION OF COMPUTED TOMOGRAPHY IMAGE¹⁾

ZHONG Jiangcheng²⁾, ZHOU Hongwei³⁾, REN Weiguang, WANG Zihui, CHEN Ji

School of Mechanics and Civil Engineering, China University of Mining and Technology (Beijing), Beijing 100083, China

Abstract

The threshold segmentation of CT image of coal is the prerequisite for the reconstruction of 3D real structure. A porosity calculating method was proposed based on grey level distribution of CT images obtained under uniaxial compression. Combined with the gray scale histogram, CT images were experienced three-value-segmentation of fracture, coal matrix and inclusion by inverse process from experimental porosity. The three-value-segmentation results obtained by Otsu method, maximum entropy method and the improved gray histogram method show that maximum entropy method can not get a good segmentation between coal inclusion and coal matrix, and Otsu method can not effectively split crack and coal matrix for its low sensitivity to dramatic variation for gray value in small area of the image. However, the improved gray histogram method eliminates the disadvantages of the two methods mentioned above. The threshold of coal matrix and inclusion is obtained with higher reliability and accuracy based on the experimental results, and a better three-value-segmentation of image is obtained by the improved grey histogram method.

Key words: computed tomography image; three-value-segmentation; porosity; Otsu method; maximum entropy method; improved gray histogram method

文章图片

煤体是一种成分复杂的沉积岩多孔介质, 内部结构主要由固体基质和孔隙裂隙两部分组成. 固体基质一般可以细分为煤基质和煤杂质, 对煤体物理力学性质起决定性作用. 孔隙度是表征岩石物性的重要参数, 直接影响岩石的物理力学性质和渗透特性^[1]. 通过CT技术对煤体真实结构进行三维重构, 并准确地建立数值几何模型已成为数值模拟发展的一种重要方法^[2].

众多学者基于CT实验展开了各种各样的研究, 例如单轴压缩过程中孔隙度演化特征研究^{[1, 3]}, 含裂隙煤体、沥青混合料、土石混合体等材料真实结构的三维重构研究^{[4, 5, 6, 7]}以及基于三维真实结构重构的数值模拟研究^{[2, 8, 9]}等. 在这些研究中, 由于研究对象及目的的差异, 并没有通用的CT图像阈值分割方法. 张引^[10]、Talab等^[11]通过最大类间方差法(Otsu法)分别对不同的车牌照图像和混凝土图像进行了二值化研究; Brink^[12]、黄春艳等^[13]运用几种基于熵的图像算法分别对灰度图像进行了阈值分割, 并对比分析了其二值化效果; 万成等^[6]、宫伟力等^[14]通过灰度直方图具有双峰特点对CT图像进行了二值化; 于庆磊等^[2]提出了环状分区与分割阈值自动识别相结合的CT图像分割算法, 即整个图像灰度直方图不具有双峰, 但各个环状分区的灰度直方图具有双峰; Glasbey^[15]提出了11种基于图像灰度直方图的阈值分割算法; 张青成等^[4]认为含裂隙图像二值化效果最佳的阈值为孔隙面积与灰度阈值间的关系曲线拐点对应的灰度值; Peng 等^[16]利用实测孔隙度反推法对9种不同岩石CT切片进行了孔隙--固体基质二值化, 并得到了较好的二值化效果. Zhou等^[17]结合实测孔隙度反推法和最大熵法对不同应力状态下煤样的CT图像进行了孔隙--固体基质二值化. 刘东菊^[18]和张琪^[19]专门针对图像阈值分割方法展开了详细研究; 吴一全等^[20]对近20年来图像阈值分割方法研究进展进行了总结. 以上大多针对图像进行二值化, 关于三值化的研究很少. 在这些阈值分割方法中, 灰度直方图双峰法是最简单实用的方法, 但是很多CT图像的灰度分布并不具有双峰特征; 理论算法中, 最大类间方差法和最大熵法应用比较广泛且阈值分割效果较好, 但是具有各自的最佳适用条件和范围; 而实测孔隙度反推法作为一种基于实验结果的非理论阈值分割方法, 具有较高的准确度和可信度, 但是当材料中含裂隙时, 很难通过实验测量其孔隙度来确定CT图像的分割阈值.

本文提出了一种基于灰度分布的孔隙度反推法的图像双阈值分割方法, 成功地区分了煤体裂隙、煤基质和煤杂质, 并给出了基于煤体CT图像灰度分布的微观孔隙度计算方法. 此阈值分割方法为准确建立含杂质煤体三维真实结构的几何模型提供了一种可靠方法, 有利于提高数值模拟的准确度.

1 基于CT图像灰度分布的孔隙度计算

当一定能量的X射线射入具有一定密度的物体并穿出时, 由于反射、衍射以及物质吸收等作用导致X射线穿出时的能量相对于入射时的能量有所减弱, 而且物质密度不同导致X射线能量衰减量也不同. 根据Lambert--Beer’ s定律, X射线穿透物体前后的能量与物质密度满足如下关系^{[21, 22]}

$I = I_{0} \exp (- μ Δ x) (1)$

其中, $I_{0}$ 和 $I$ 分别为X射线入射和出射时的强度, $Δ x$ 为X射线穿透材料的厚度, $μ$ 为材料对X射线线性衰减系数.

X射线衰减系数 $μ$ 主要取决于康普顿散射和光电吸收, 可表达为^{[21, 22]}

$μ = ρ (a + b Z^{3.8} / E^{3.2}) (2)$

其中, $ρ$ 为材料密度, $a$ 为Klein--Nishina系数, $b$ 为常数, $Z$ 为等效原子序数, $E$ 为X射线粒子能量.

在医用CT中, 通过CT(Hounsfield)数来描述人体组织的衰减系数相对于水的衰减系数的大小, 其定义为^{[21, 22]}

$CT_{\rm number} = \dfrac{\mu - \mu _{\rm w} }{\mu _{\rm w} }\times 1\, 000 \quad(3)$

其中, $μ_{w}$ 为水对X射线的线衰减系数. 若采用医用CT数定义, 则比水轻的材料CT数为负值, 如空气的CT值为 $- 1 000$ , 使用起来很不方便. 为此可以将工业CT的CT数定义为^[23]

$CT_{\rm number} = \dfrac{\mu - \mu _{\rm w} }{\mu _{\rm w} }\times 1\, 000 + 1\, 000 \quad(4) $

当扫描样品直径为25 mm时, 一般工业 $μ$ CT的像素分辨率为几微米到几十微米, 但是一般致密岩石(煤、砂岩、盐岩等)的孔隙尺寸大部分分布在几纳米到几百纳米, 微米级别的大孔隙已经很少^{[16, 22]}. 因此, 岩石CT图像因为分辨率限制而很难准确识别孔隙, 只能识别一些大孔及裂隙. 鉴于此, 如果岩石的固体基质足够小, 岩石材料可以看成由固体基质和固体基质间孔隙组成, 并设基质的密度为常数 $ρ_{m}$ , 空气的密度为 $ρ_{a}$ , 则材料中任意像素大小处的体积密度为

$ρ = (1 - ϕ) ρ_{m} + ϕ ρ_{a} (5)$

其中, $ϕ$ 为材料任意像素大小中的孔隙度. 空气密度 $ρ_{a}$ 近似看成0, 则

$CT_{\rm number} = \dfrac{1\, 000\rho _{\rm m} \left( {a + bZ^{3.8} / E^{3.2}} \right)}{\mu _{\rm w} }(1 - \phi ) \quad(6) $

当材料基质以及X射线参数一定时, 式(6)等式右边中仅仅只有孔隙度 $ϕ$ 为变量, 说明CT数与孔隙度成负线性相关性.

计算机(Matlab)默认的图像存储位数为8位, 煤样CT切片图像为16位, 为了便于后期图像处理与计算, 在Matlab中将CT图像利用线性映射算法^[16]转换成对应的8位灰度图像, 其灰度值范围为0 $~$ 255, 物质密度越大对应的灰度值越大. 孔隙对应的像素灰度值最低为0, 图像中显示为黑色; 杂质集中区域灰度值较大, 显示为灰白色. Karacan等^[22]指出煤中主要存在三种基本的结构类型: 由基质孔隙度控制的有机煤基质区域, 由煤体孔隙度控制的有机煤基质区域和含少量孔隙的杂质(矿物)集中区域. 因此, 本文忽略杂质集中区域中相对较少的孔隙, 仅计算有机煤基质区域中的孔隙. 纯煤基质对应的像素灰度值设为 $T_{m} (0 < T_{m} < 255)$ . 根据式(6), 煤基质单个像素的灰度值与孔隙度成负线性相关性, 可以表示为 $T_{i} = - T_{m} ϕ + T_{m}$ . 则基于单张CT图像计算得到的孔隙度 $ϕ$ 为

$\phi \left( {T_{\rm m} } \right) = \dfrac{V_{\rm pore} }{V_{\rm coal} } = \dfrac{\sum_{T_{\rm i} = 0}^{T_{\rm m} } {n_{\rm i} s^2d\dfrac{T_{\rm m} - T_{\rm i} }{T_{\rm m} }} }{Ns^2d} = \\ \qquad \sum_{T_{\rm i} = 0}^{T_{\rm m} } {\dfrac{T_{\rm m} - T_{\rm i}}{T_{\rm m} }} p(T_{\rm i} ) \quad(7) $

其中, $p (T_{i}) = n_{i} / N$ , $n_{i}$ 为灰度值等于 $T_{i}$ 的像素个数, $N$ 为煤样切片图像的像素总个数, $s$ 为像素分辨率, $d$ 为扫描间距(垂直切片方向像素分辨率).

若每个煤样扫描 $M$ 张图片, 则整个煤样的孔隙度为

$ \varPhi (T_{\rm m} ) = \dfrac{1}{M}\sum_{j = 0}^M {\phi _j } \left( {T_{\rm m} } \right) \quad(8)$

公式(8)计算结果虽然相对偏小, 但其充分利用了像素灰度信息, 克服了CT图像分辨率低的问题. 因此, 从测量精度上看此公式计算得到的孔隙度更接近于实验测得的孔隙度. 此公式只要知道煤基质的灰度值 $T_{m}$ , 就能计算出煤样孔隙度 $Φ$ , $T_{m}$ 可以通过压汞实验测得的孔隙度进行反推得到, 且认为同一地质条件下 $T_{m}$ 基本保持不变, 近似看成常数.

2 煤体CT图像三值化方法

对于含杂质的裂隙煤体(图1(a)), 有时除了需要得到孔隙裂隙结构外, 还需要得到煤基质和煤杂质各自的分布(例如精确的数值建模等), 此时对于CT切片进行阈值分割时不仅仅是区分裂隙和固体基质(二值化), 还需要将固体基质中的煤基质和煤杂质分开(三值化). 三值化原理如下

$ T'\left( {x, y} \right) = \left\{ \!\!\begin{array}{ll} {T_1 '}, & T\left( {x, y} \right) < T_1 \\ {T_2 '}, & T_1 \leq T\left( {x, y} \right) \leq T_2 \\ {T_3 '}, & T\left( {x, y} \right) > T_2 \end{array} \right. \quad (9)$

	Figure Option View Download New Window
	图1 含杂质裂隙煤体原始灰度图像及去噪后图像

其中, $T (x, y)$ , $T' (x, y)$ 分别为原始图像和三值化后图像中各个像素点的灰度值, $T_{1}$ 和 $T_{2}$ 为三值化的灰度阈值, ${T'}_{1}$ , ${T'}_{2}$ , ${T'}_{3}$ 为三值化的三个灰度值, 对应像素区域分别代表裂隙、煤基质和煤杂质.

为了排除噪声干扰, 在进行阈值分割前对原始图像进行3× 3窗口的中值滤波. 计算时为了区分背景与裂隙, 将背景颜色改为白色, 去噪后的白色背景图像如图1(b).

本文利用基于图像灰度直方图的孔隙度反推法对含杂质裂隙煤体进行了三值化, 并与最大类间方差法和最大熵法进行了对比分析.

2.1 最大类间方差法(Otsu法)

最大类间方差法^[24]是1975年日本学者大津展之提出的, 其根据图像预分割的各个部分灰度值之间的总方差最大来确定最佳阈值. 总方差值越大, 说明各个部分间差别越大, 分割出的图像的各个部分越准确.

三值化后图像被分为3个区域, 设分割阈值为 $T_{1}$ 和 $T_{2}$ , 则各区域占整幅图像的面积比为

$\left.\begin{array}{1} P_1 = \sum_{T_{\rm i} = 0}^{T_1 } p(T_{\rm i} ) \\ P_2 = \sum_{T_{\rm i} = T_1 + 1}^{T_2 } p(T_{\rm i} ) \\ P_3 = \sum_{T_{\rm i} = T_2 + 1}^{255} p(T_{\rm i} ) \end{array} \right\} \quad(10) $

整幅图像的加权平均灰度为

$T_{\rm I} = \sum_{T_{\rm i} = 0}^{255} {T_{\rm i} p(T_{\rm i} )} \quad(11)$

每个区域的加权平均灰度为

$T_{\rm I1} = \sum_{T_{\rm i} = 0}^{T_1 } {T_{\rm i} p(T_{\rm i} )} / P_1 \quad (12a) \\ T_{\rm I2} = \sum_{T_{\rm i} = T_1 + 1}^{T_2 } {T_{\rm i}p(T_{\rm i} )} / P_2 \quad (12b) \\ T_{\rm I3} = \sum_{T_{\rm i} = T_2 + 1}^{255} {T_{\rm i} p(T_{\rm i} )} / P_3 \quad (12c)$

则, 各个区域间的总方差可以表示为

$σ^{2} (T_{1}, T_{2}) = P_{1} {(T_{I 1} - T_{I})}^{2} + P_{2} {(T_{I 2} - T_{I})}^{2} +$

$P_{3} {(T_{I 3} - T_{I})}^{2} (13)$

总方差 $σ^{2} (T_{1}, T_{2})$ 的值越大表示阈值选取得越好, 即最佳阈值为

$(T_{1}, T_{2}) = \underset{0 < T_{1} < T_{2} < 255}{argmax} {σ^{2} (T_{1}, T_{2})} (14)$

2.2 最大熵法

最大熵法^[25]是1985年由Kapur提出的, 其阈值选取准则是图像分割后的各个部分的总熵值最大, 即信息量最大. 这种方法的目的在于将图像的灰度直方图分成独立的类, 使得各类的总熵最大.

同样, 对图像进行三值化, 假设双阈值分割后图像分为3个区域, 阈值为 $T_{1}$ 和 $T_{2}$ , 那么图像分割后各部分的熵之和为

$ H(T_1 , T_2 ) = - \sum_{T_{\rm i} = 0}^{T_1 } {\dfrac{p(T_{\rm i} )}{P_1 }} \ln \dfrac{p(T_{\rm i} )}{P_1 } - \\ \qquad \sum_{T_{\rm i} = T_1 + 1}^{T_2 } {\dfrac{p(T_{\rm i} )}{P_2 }} \ln \dfrac{p(T_{\rm i} )}{P_2 } - \\ \qquad \sum_{T_{\rm i} = T_2 + 1}^{255} {\dfrac{p(T_{\rm i} )}{P_3 }} \ln \dfrac{p(T_{\rm i} )}{P_3 } \quad(15) $

最佳阈值的选取准则为

$(T_{1}, T_{2}) = \underset{0 < T_{1} < T_{2} < 255}{argmax} {H (T_{1}, T_{2})} (16)$

2.3 本文改进的灰度直方图法

有时通过图像灰度直方图中某灰度值出现概率是否突变可以简单直接判断出图像阈值分割点, 比如当灰度直方图出现双峰间的明显波谷时, 则波谷处灰度值一般为阈值分割点. 三值化时若灰度直方图具有三峰, 三峰间的两个波谷则为两个阈值点. 一般天然岩石CT图像灰度直方图都不具有双峰或三峰, 例如图2为含杂质裂隙煤体的灰度直方图, 并没有双峰, 其中, 裂隙与固体基质间的阈值分割点 $T_{1}$ (突变点)可以较容易得到. 由于满足此阈值的灰度区间(图2圆形区域)范围很小, 且灰度阈值前后对应的 $p (T_{i})$ 差值应该很大, 因此在此区域范围内, 取 $p (T_{i})$ 相对增长速率最快的点对应的灰度值作为 $T_{1}$ 的值. 然而, 煤基质与煤杂质间的阈值分割点无法从灰度直方图直接找到. 为此, 在CT图像分辨率不够的情况下, 本文充分利用CT图像像素的灰度信息, 依据像素灰度值与物质平均密度线性正相关, 基于CT图像灰度近似计算了煤样孔隙度. 可见, 像素灰度值大小与其中的孔隙度线性负相关, 若一个像素空间里全部为纯煤基质, 则该像素的灰度值为煤基质的最大灰度值(灰度上限), 即该值可以作为煤基质和煤杂质的阈值分割点. 那么, 将实验测得的孔隙度值代入式(8), 反算得到煤基质和煤杂质间的阈值 $T_{m}$ , 此方法称作基于图像灰度直方图的孔隙度反推法. 即, $T_{2} = T_{m} = f^{- 1} (Φ_{0})$ , $f^{- 1} (Φ)$ 为式(8)的反函数, $Φ_{0}$ 为目标煤样的孔隙度.

	Figure Option View Download New Window
	图2 含杂质裂隙煤体的灰度直方图

3 算例与分析

3.1 取样

本文实验煤块取自平煤12K己15-17220采面, 埋深约900m, 煤层牌号为焦煤, 煤质松软, 煤体呈粉末状, 半光亮型, 属脆性煤体. 煤块经严密封装后, 送至室内试验室, 加工成直径为25 mm的圆柱形试样. 在同一块煤相同位置上钻取1个高为50 mm以及3个高为25 mm的圆柱形试件, 分别用于单轴加载的CT扫描实验和压汞实验.

3.2 CT扫描及压汞实验

取直径为25 mm、高为50 mm的圆柱形煤样进行单轴加载下的CT扫描实验. CT实验设备采用中国矿业大学(北京)煤炭资源与安全开采国家重点实验室与美国BIR公司合作研发的ACTIS300-320/225CT/DR高分辨率工业CT实时成像系统, 并研发了配套的单轴加载设备. 在单轴压缩条件下对煤样加载到不同应力水平, 进行CT扫描实验, 直至试件破坏. 煤样扫描载荷依次为0, 3, 4, 4.5, 5, 6, 7, 8, 9, 10 MPa. 每个载荷下扫描1001张图片, 扫描间距和切片厚度均为0.05 mm, 图像尺寸为750× 750像素, 视场为35 mm, 每个像素大小(分辨率)为 46× 46$\mu$m².

对3个圆柱形煤样M1、M2、M3进行压汞实验. 实验在中国石油勘探开发研究院北京院区进行, 实验仪器采用美国core lab CMS300和美国Auto Pore Ⅳ 9505压汞仪. 试验进行前先通过普通渗透仪测得煤样的孔隙度, 结果如表1所示. 压汞实验测得的三个煤样的孔隙尺寸分布如图3所示, 可以看出孔喉尺寸主要集中在0.008 $~$ 5 m 区域, 区域两端的孔隙比较多, 且大部分为纳米和微米级孔隙, 符合式(5)推导的假设. 然而, CT图像单个像素分辨率为46$\mu$m, 因此通过CT像素识别的孔隙分布是不准确的, 而通过式(8)计算的煤样孔隙度是可以精确到煤基质尺度的, 具有更高的精度, 与压汞法测得的孔隙度的精度更接近. 虽然压汞法测的是开孔孔隙度, 相对于CT测的全孔孔隙度值偏小, 但是CT计算的孔隙度式(8)中同样忽略了少量煤杂质中的孔隙, 结果也是偏小的. 因此, 定性分析可以知道采用压汞法测得的孔隙度代入式(8)反算得到的煤基质与煤杂质间的阈值 $T_{2}$ 的误差相对被缩小.

表1 压汞实验煤样的物理几何参数

煤样编号	直径/cm	高度/cm	质量/g	密度/(g $\cdot$ cm< sup> -3< /sup> )	孔隙度/ M1	2.485	2.449	15.763	1.360	5.770
M2	2.452	2.438	15.611	1.470	5.580
M3	2.461	2.414	14.825	1.590	5.770

表1 压汞实验煤样的物理几何参数

	Figure Option View Download New Window
	图3 煤样孔隙分布图

3.3 三值化结果与分析

为了简化计算, 等间距地取5张(第300, 400, 500, 600, 700张)加载前不含裂隙的CT扫描图片进行孔隙度计算. 5张图片孔隙度与临界灰度值 $T_{m}$ 间的关系曲线如图4所示. 取压汞实验得到的3个煤样的平均孔隙度作为整个煤块的孔隙度, 即 $Φ_{0} = 5.71 %$ , 代入式(8)解得 $T_{2} = T_{m} = 90$ .

	Figure Option View Download New Window
	图4 CT图像孔隙度随临界灰度值的变化曲线

为了说明上述简化计算带来的误差, 当 $T_{m} = 90$ 时, 通过式(7)计算了等间距的25张CT切片的孔隙度, 结果如图5所示. 通过式(8)计算得到的试件的孔隙度 $Φ = 5.69 %$ , 与压汞测得的真实孔隙度的相对误差为0.35说明这样简化计算是可行的.

	Figure Option View Download New Window
	图5 孔隙度随切片位置的变化曲线( $T_{m} = 90$ )

表2和表3分别为等间距5张不含裂隙CT图像和含裂隙CT图像在最大类间方差法、对数最大熵法以及本文改进灰度直方图3种方法下的三值化结果, 其对应单轴压缩载荷分别为0 MPa和5 MPa. 由于5张图片阈值结果很接近, 取加载前后第500张(A500)图片为特例来说明3种方法下CT图像的三值化效果. 图6和图7分别给出了不含裂隙煤体CT图像和含裂隙煤体CT图像的原始图像以及最大类间方差法、对数最大熵法和本文改进灰度直方图3种方法下的三值化图像.根据图6和图7并结合表2和表3数据可知, 最大类间方差法可以粗略分割出煤基质和煤杂质, 相对本文改进灰度直方图法的求解阈值偏小, 且无法分割出裂隙和固体基质; 当裂隙较为明显时, 对数最大熵法能够较好地区分出裂隙与固体基质, 但是固体基质中的煤基质和煤杂质无法完全区分; 本文改进的灰度直方图法既能很好区分出裂隙与固体基质, 也能很好区分出煤基质和煤杂质. 图6(b)和图7(b)以及表格中数据表明, 最大类间方差法的阈值分割点基本不受裂隙产生的影响, 即图片中微小区域的灰度突变其无法检测到, 而仅仅是对图片中主要成分组成进行检测区分, 这与文献[26]的观点“ Otsu方法适用于图像中被分割各部分的面积所占比例接近的情况, 否则, Otsu分割方法将失效” 是一致的. 对比图6(c)和图7(c)以及表格中数据可以发现, 对数最大熵法阈值分割点对裂隙的产生很敏感, 能够较好检测区分, 同时发现其能检测出煤杂质中密度相对较高的区域, 说明其主要是识别区分相对主要区域变化较大的微小区域, 即对灰度突变的区域(裂隙的产生, 高密度杂质等)敏感性高. 因此, 这两种较为成熟的理论阈值分割算法都有各自的主要优势和劣势, 应该根据研究重点和需要正确选取和应用. 分别对比图6和图7可知, 无论是否含裂隙, 本文改进灰度直方图法三值化效果均很好, 避免了最大类间方差法和最大熵法的劣势, 且煤基质与煤杂质间的阈值由实验实测孔隙度反推得到, 阈值求解具有更高的准确度.

表2 不含裂隙CT图像的三值化阈值分割点

双阈值分割方法	不含裂隙CT图像三值化阈值点 $(T_{1}, T_{2})$
	A300	A400	A500	A600	A700
最大类间方差法	(76, 86)	(77, 87)	(74, 88)	(74, 88)	(76, 86)
对数最大熵法	(107, 161)	(112, 171)	(109, 169)	(110, 167)	(107, 160)
改进灰度直方图法	(63, 90)	(64, 90)	(62, 90)	(62, 90)	(63, 90)

表2 不含裂隙CT图像的三值化阈值分割点

表3 含裂隙CT图像三值化阈值分割点

双阈值分割方法	含裂隙CT图像三值化阈值点 $(T_{1}, T_{2})$
	A300	A400	A500	A600	A700
最大类间方差法	(76, 86)	(77, 87)	(74, 88)	(74, 88)	(76, 86)
对数最大熵法	(59, 133)	(61, 125)	(57, 129)	(59, 115)	(62, 136)
改进灰度直方图法	(63, 90)	(64, 90)	(62, 90)	(62, 90)	(63, 90)

表3 含裂隙CT图像三值化阈值分割点

	Figure Option View Download New Window
	图6 不含裂隙煤体CT灰度图像(A500)及三值化后图像

	Figure Option View Download New Window
	图7 含裂隙煤体CT灰度图像(A500)及三值化后图像

基于本文三值化方法, 对各个载荷下煤样的1 000张CT图像进行三值化, 然后导入到Mimics中进行三维重构, 能够分别得到孔裂隙、煤基质和煤杂质的三维空间分布情况. 以5 MPa载荷下的CT图像为例, 图8给出了基于本文三值化方法的煤样(裂隙--煤基质--煤杂质)三维重构的真实几何模型. 将三维几何模型导入到数值软件中, 并根据数值模拟需求对不同组分赋予不同的物理力学参数或控制方程, 则可以实现煤体真实结构的数值计算与分析. 这对提高数值模拟的准确度具有重要意义.

	Figure Option View Download New Window
	图8 5 MPa载荷下煤样的三维真实结构

4 结论

(1) 在CT图像像素分辨率无法辨别孔隙的情况下, 充分利用图像单个像素的灰度信息, 根据灰度值与材料密度的线性关系提出了一种基于CT图像灰度分布的孔隙度计算方法, 该方法能够计算到微观尺度的孔隙. 只要确定煤基质的灰度值, 就能计算出任意一张CT图像的孔隙度, 进而计算整个煤样的孔隙度.

(2) 通过实验测得的孔隙度结合本文提出的孔隙度计算方法反推煤基质与煤杂质的阈值分割点, 同时通过CT图像灰度直方图得到裂隙与煤基质的阈值分割点, 以此对煤体CT图像成功进行了裂隙--煤基质--煤杂质三值化, 是一种实用而可靠的三值化方法.

(3) 最大类间方差法只能粗略区分出煤基质和煤杂质, 对灰度突变的区域(裂隙的产生)敏感性低, 即无法区分裂隙与固体基质; 而最大熵法仅能较好地区分出裂隙与煤基质, 却无法区分出煤基质与煤杂质. 相比而言, 本文改进的灰度直方图方法避免了以上两种方法的劣势, 煤基质与煤杂质间的阈值基于实验结果得到, 对于含杂质煤体CT图像, 具有更好的三值化效果, 且可信度和准确度更高.

(4)本文方法计算量较大, 寻找类似最大类间方差法和最大熵法的理论计算方法来准确对类似含裂隙含杂质煤体的CT图像进行三值化值得研究. 依据本文三值化方法对煤体真实结构进行三维重构建模并进行力学或渗流等数值模拟值得进一步研究和讨论.

The authors have declared that no competing interests exist.

参考文献

文献选项

[1]	王宇, 李晓, 阙介民等. 基于CT图像灰度水平的孔隙率计算及应用. 水利学报, 2015, 46(3): 357-365 [本文引用:2]
[2]	于庆磊, 杨天鸿, 唐世斌等. 基于CT 的准脆性材料三维结构重建及应用研究. 工程力学, 2015, 32(11): 51-62 [本文引用:3]
[3]	李建胜, 王东, 康天合. 基于显微CT试验的岩石孔隙结构算法研究. 岩土工程学报, 2010, 32(11): 1703-1708 [本文引用:1]
[4]	张青成, 王万富, 左建民等. 煤岩CT图像二值化阀值选取及三维重构技术研究. CT理论与应用研究, 2014, 23(1): 45-51 [本文引用:2]
[5]	李果, 张茹, 徐晓炼等. 三轴压缩煤岩三维裂隙CT图像重构及体分形维研究. 岩土力学, 2015, 36(6): 1633-1642 [本文引用:1]
[6]	万成, 张肖宁, 贺玲凤等. 基于真实细观尺度的沥青混合料三维重构算法. 中南大学学报(自然科学版), 2012, 43(7): 2813-2820 [本文引用:2]
[7]	孙华飞, 鞠杨, 行明旭等. 基于CT图像的土石混合体破裂——损伤的三维识别与分析. 煤炭学报, 2014, 39(3): 452-459 [本文引用:1]
[8]	王刚, 杨鑫祥, 张孝强等. 基于CT三维重建的煤层气非达西渗流数值模拟. 煤炭学报, 2016, 41(4): 931-940 [本文引用:1]
[9]	王聪. 基于CT技术的沥青混合料力学行为及三维重构研究. [硕士论文]. 大连: 大连海事大学, 2012 [本文引用:1]
[10]	张引. 基于空间分布的最大类间方差牌照图像二值化算法. 浙江大学学报(工学版), 2001, 35(3): 272-275 [本文引用:1]
[11]	Talab AMA, Huang Z, Xi F, et al. Detection crack in image using Otsu method and multiple filtering in image processing techniques. Optik-International Journal for Light and Electron Optics, 2016, 127(3): 1030-1033 [本文引用:1]
[12]	Brink AD. Thresholding of digital images using two-dimensional entropies. Pattern Recognition, 1992, 25(8): 803-808 [本文引用:1]
[13]	黄春艳, 杨国胜, 侯艳丽. 基于熵的图像二值化方法比较研究. 河南大学学报(自然科学版), 2005, 35(2): 76-78 [本文引用:1]
[14]	宫伟力, 安里千, 赵海燕等. 基于图像描述的煤岩裂隙CT图像多尺度特征. 岩土力学, 2010, 31(2): 371-376 [本文引用:1]
[15]	Glasbey CA. An analysis of histogram-based thresholding algorithms. CVGIP: Graphical Models and Image Processing, 1993, 55(6): 532-537 [本文引用:1]
[16]	Peng RD, Yang YC, Ju Y. Computation of fractal dimension of rock pores based on gray CT images. Chinese Science Bulletin, 2011, 56(26): 2256-2266 [本文引用:3]
[17]	Zhou HW, Zhong JC, Ren WG, et al. Characterization of pore-fracture networks and their evolution at various measurement scales in coal samples using X-ray µCT and a fractal method. International Journal of Coal Geology, 2018, 189: 35-49 [本文引用:1]
[18]	刘东菊. 基于阈值的图像分割算法的研究. [博士论文]. 北京: 北京交通大学, 2009 [本文引用:1]
[19]	张琪. 结合边缘检测的图像二值化算法. [硕士论文]. 长春: 吉林大学, 2011 [本文引用:1]
[20]	吴一全, 孟天亮, 吴诗婳. 图像阈值分割方法研究进展20年(1994—2014). 数据采集与处理, 2015, 30(1): 1-23 [本文引用:1]
[21]	Hsieh J. Computed Tomography: Principles, Design, Artifacts, and Recent Advances. Bellingham: SPIE Press, 2003 [本文引用:3]
[22]	Karacan CO, Okand an E. Adsorption and gas transport in coal microstructure: investigation and evaluation by quantitative X-ray CT imaging. Fuel, 2001, 80(4): 509-520 [本文引用:5]
[23]	Van GM, Swennen R, Wevers M. Towards 3-D petrography: application of micro-focus computer tomography in geological science. Computer & Geoscience, 2001, 27(9): 1091-1099 [本文引用:1]
[24]	Otsu N. A threshold selection method from gray-level histograms. IEEE Transactions on Systems, Man, and Cybernetics, 1979, 9(1): 62-66 [本文引用:1]
[25]	Kapur JN, Sahoo PK, Wong AKC. A new method for gray–level picture thresholding using the entropy of the histogram. Computer Vision, Graphics, and Image Processing, 1985, 29(3): 273-285 [本文引用:1]
[26]	周铭. 低信噪比红外小目标图像的分割方法. 激光与红外, 2004, 34(3): 225-228 [本文引用:1]

2015

0.0

王宇, 李晓, 阙介民等. 基于CT图像灰度水平的孔隙率计算及应用. 水利学报, 2015, 46(3): 357-365

摘　要：基于 CT图像的材料孔隙率分析是评价岩土多孔介质的重要手段。当前孔隙率的计算主要针对 CT单张图像采用阈值分割法计算，然而，由于阈值分割法在确定孔洞与固体基质灰度门槛值时的不确定性，即便针对相同的 CT图像由于不同的阈值分割方法而使计算结果相差较大。为了克服阈值分割法确定材料孔隙率的缺陷，探索一种基于 CT图像确定岩土介质孔隙的新方法。首先，借助于数字地面高程基本思路，提出基于 CT图像灰度水平的孔隙率计算方法，从图像的 CT数出发将 CT数分成两类，一类属于孔洞，另一类属于基质，对孔隙计算相关公式进行推导。其次，将灰度水平法应用于土石混合体材料的 CT单张图像孔隙率计算当中，所得计算结果并与两种阈值分割法进行对比研究，证明了灰度水平方法的优势与可行性。最后，在得到试样孔隙率的基础上，从试样细观孔隙演化角度对试样的变形破坏机理进行了分析。

... 孔隙度是表征岩石物性的重要参数, 直接影响岩石的物理力学性质和渗透特性^[1] ...

... 众多学者基于CT实验展开了各种各样的研究, 例如单轴压缩过程中孔隙度演化特征研究^[1,3], 含裂隙煤体、沥青混合料、土石混合体等材料真实结构的三维重构研究^[4,5,6,7]以及基于三维真实结构重构的数值模拟研究^[2,8,9]等 ...

2015

0.0

于庆磊, 杨天鸿, 唐世斌等. 基于CT 的准脆性材料三维结构重建及应用研究. 工程力学, 2015, 32(11): 51-62

摘　要：准脆性材料的结构及其细观组分的力学性质决定了材料的力学性质及破裂机制,在数值模型中尽可能准确的考虑材料的真实结构已成为了数值计算与数值模拟发展的一种趋势。该文以混凝土为研究对象,借助于先进的CT技术获取混凝土内部结构切片图像,利用数字图像处理技术实现了基于CT图像的混凝土材料结构的识别和表征,并针对CT图像具有颜色亮度不均并呈环状分布的特点,提出了环状分区与分割阈值自动识别相结合的CT图像分割算法;在此基础上,建立了基于位图矢量化理论的三维实体材料结构模型的重建方法和三维网格化材料结构模型重建方法,并将三维网格化方法建立的材料结构模型与三维岩石破裂过程分析系统RFPA3D结合进行了初步应用,对混凝土单轴压缩破裂过程进行了数值模拟。通过数值试验与物理实验结果对比发现,考虑混凝土骨料真实分布的数值实验结果与物理实验具有一定可比性,数值试验结果从得到的力学参数和破裂模式方面比较接近于物理实验结果,为深入研究混凝土、岩石、复合材料等力学特征提供了一种可行的研究思路。

... 通过CT技术对煤体真实结构进行三维重构,并准确地建立数值几何模型已成为数值模拟发展的一种重要方法^[2] ...

... 于庆磊等^[2]提出了环状分区与分割阈值自动识别相结合的CT图像分割算法, 即整个图像灰度直方图不具有双峰, 但各个环状分区的灰度直方图具有双峰 ...

2010

0.0

李建胜, 王东, 康天合. 基于显微CT试验的岩石孔隙结构算法研究. 岩土工程学报, 2010, 32(11): 1703-1708

2014

0.0

张青成, 王万富, 左建民等. 煤岩CT图像二值化阀值选取及三维重构技术研究. CT理论与应用研究, 2014, 23(1): 45-51

运用工业CT实验系统(简称为DT-3)对煤样完成单轴加载实验,并对煤样进行一系列断层扫描,获取到不同截面处的破坏CT图像。基于图像分割技术,对图像的灰度级别设定不同的二值化阀值,得出不同阀值下的孔隙面积变化曲线图;并提出以拐点处对应的阀值作为裂隙图像二值化阀值时效果最佳。由上述研究结果结合灰度线性加权平均插值方法,编程实现了煤岩CT图像的三维重构,为后续研究煤岩裂隙三维表征提供了科学依据。

... 张青成等^[4]认为含裂隙图像二值化效果最佳的阈值为孔隙面积与灰度阈值间的关系曲线拐点对应的灰度值 ...

2015

0.0

李果, 张茹, 徐晓炼等. 三轴压缩煤岩三维裂隙CT图像重构及体分形维研究. 岩土力学, 2015, 36(6): 1633-1642

通过对煤岩试件进行常规三轴力学试验,并对破坏后的煤岩试件进行CT扫描,之后将处理过的CT图片导入Mimics10.01进行三维重构,得到了煤岩试件内部裂隙空间的分布情况。另外,利用自编的Matlab程序对煤岩试件的体分形维进行了初步分析。结果表明:煤岩试件在三轴加载条件下主要发生压剪破坏;且围压不同,煤岩试件表现出了特定的力学规律;利用Mimics10.01重构得到的煤岩三维裂隙空间分布情况显示,随着围压的增大,煤岩试件剪切破坏面的连通率和体密度有增大的趋势;体分形维计算的结果显示,煤岩内部破坏越严重,体分形维值越小。

2012

0.0

万成, 张肖宁, 贺玲凤等. 基于真实细观尺度的沥青混合料三维重构算法. 中南大学学报(自然科学版), 2012, 43(7): 2813-2820

... 万成等^[6]、宫伟力等^[14]通过灰度直方图具有双峰特点对CT图像进行了二值化 ...

2014

0.0

孙华飞, 鞠杨, 行明旭等. 基于CT图像的土石混合体破裂——损伤的三维识别与分析. 煤炭学报, 2014, 39(3): 452-459

Soil- rock mixture(SRM) is a unique type of discontinuous and inhomogeneous geomaterials.It is very difficult to analyze the deformation and failure of SRM using the experimental techniques and analytical theories of conventional soil mechanics.A novel 3D identification method was proposed to identify the fracture and damage of SRM based on CT and image processing techniques.Using the self- developed identification algorithm and computer program，we extracted and analyzed the mesoscopic characteristics of interior fractures of SRM with various rock contents subject to uniaxial compressive loads.A 3D model was reconstructed to represent the fracture and damage of SRM，based on which the damage evolution and characteristics of SRM under uniaxial compression and the effects of rock contents were determined.It is shown that the properties and evolution of fracture and damage of SRM can be characterized visually and quantitatively through the proposed method.The damage evolution of SRM and the effects of rock closely relate to the rock content.There is a critical ratio of rock content.

土石混合体是一种特殊的非连续非均匀地质材料，传统的岩土力学试验方法与分析理论难以直接用于分析其变形破坏行为。发展了一种基于CT扫描和图像处理技术的内部破裂与损伤的三维识别方法，利用自行开发的裂纹识别提取程序提取和分析了单轴压缩下不同含石率土石混合体变形破坏时的内部细观结构变化特征，建立了土石混合体破裂-损伤的三维重构模型，定量地分析了单轴压缩条件下土石混合体的内部损伤演化规律以及含石量的影响和作用。研究表明：该方法可以直观定量地表征土石混合体变形破坏的内部损伤特征与演化规律。石块对土石混合体损伤发展的影响与含石量密切相关，存在一个临界含石率。

2016

0.0

王刚, 杨鑫祥, 张孝强等. 基于CT三维重建的煤层气非达西渗流数值模拟. 煤炭学报, 2016, 41(4): 931-940

2012

0.0

2001

0.0

张引. 基于空间分布的最大类间方差牌照图像二值化算法. 浙江大学学报(工学版), 2001, 35(3): 272-275

车辆牌照识别VLPR(Vehicle License Plate Recognition)是智能交通系统ITS（IntelligentTransportation System）的重要部分，有着巨大的社会和经济效益。车辆牌照图像二值化方法的研究是VLPR中的关键技术，通常都采用经典的Bernsen算法和Otsu算法。但由于光照不均、摄像机畸变、曝光不足、动态范围太窄和车辆牌照被污染等原因，车辆牌照图像的质量往往不佳，常存在严重伪影和字符边缘模糊，极大地影响了牌照图像二值化效果，Bernsen算法和Otsu算法也不能很好地克服上述问题。为此，提出了一种新的牌照图像二值化算法CASDA(Cluster Algorithm Basedon Spatial Distributi on Analysis)，把牌照图像空间分布特性与最大类间方差的统计特性有机结合，大大加快了最大类间方差模式分类速度。与Bernsen、Otsu、LEVBB等算法比较，CASDA具有速度快，能消除不均匀光照引起的伪影，极少出现笔划断裂等优点，二值化效果好。

... 张引^[10]、Talab等^[11]通过最大类间方差法(Otsu法)分别对不同的车牌照图像和混凝土图像进行了二值化研究 ...

2016

0.0

... 张引^[10]、Talab等^[11]通过最大类间方差法(Otsu法)分别对不同的车牌照图像和混凝土图像进行了二值化研究 ...

1992

0.0

... Brink^[12]、黄春艳等^[13]运用几种基于熵的图像算法分别对灰度图像进行了阈值分割,并对比分析了其二值化效果 ...

2005

0.0

黄春艳, 杨国胜, 侯艳丽. 基于熵的图像二值化方法比较研究. 河南大学学报(自然科学版), 2005, 35(2): 76-78

对基于对数和指数图像熵的最大熵法和最小交叉熵法的图像二值化算法进行了分析. 通过仿真,对算法进行了比较研究,阐明了算法的局限性和比较适合的应用场合.

... Brink^[12]、黄春艳等^[13]运用几种基于熵的图像算法分别对灰度图像进行了阈值分割,并对比分析了其二值化效果 ...

2010

0.0

宫伟力, 安里千, 赵海燕等. 基于图像描述的煤岩裂隙CT图像多尺度特征. 岩土力学, 2010, 31(2): 371-376

摘　要：基于煤岩孔隙系统多尺度结构特征对深入认识多尺度流体运移机制的重要性，提出了基于图像描述的煤岩CT图像孔隙结构的多尺度精细描述方法。采用了图像的多策略分割技术提取目标，利用Freeman链码对目标的边界进行表达，研究了由形态学、统计矩、链码、计盒维数构造目标之间的关系、目标占有区域与边界的图像描绘子、以及分形描绘子；综合运用上述方法对煤岩CT图像中的大尺度宏观裂纹目标、小尺度细观裂隙目标进行了识别。结果表明，宏观裂纹可由灰度闽值法实现目标提取；小尺度细观裂隙需采用较复杂的分割策略，如基于索贝尔梯度算子的分水岭变换：进一步应用链码表达、图像描绘予和分形描绘子，实现了煤岩孔隙结构在欧氏空间与分形空间的多尺度精确描述。

... 万成等^[6]、宫伟力等^[14]通过灰度直方图具有双峰特点对CT图像进行了二值化 ...

1993

0.0

... Glasbey^[15]提出了11种基于图像灰度直方图的阈值分割算法 ...

2011

0.0

... Peng 等^[16]利用实测孔隙度反推法对9种不同岩石CT切片进行了孔隙--固体基质二值化, 并得到了较好的二值化效果 ...

... 当扫描样品直径为25 mm时,一般工业 μCT的像素分辨率为几微米到几十微米, 但是一般致密岩石(煤、砂岩、盐岩等)的孔隙尺寸大部分分布在几纳米到几百纳米, 微米级别的大孔隙已经很少^[16,22] ...

... 计算机(Matlab)默认的图像存储位数为8位,煤样CT切片图像为16位,为了便于后期图像处理与计算,在Matlab中将CT图像利用线性映射算法^[16]转换成对应的8位灰度图像,其灰度值范围为0 ~255, 物质密度越大对应的灰度值越大 ...

2018

0.0

... Zhou等^[17]结合实测孔隙度反推法和最大熵法对不同应力状态下煤样的CT图像进行了孔隙--固体基质二值化 ...

2009

0.0

... 刘东菊^[18]和张琪^[19]专门针对图像阈值分割方法展开了详细研究 ...

2011

0.0

... 刘东菊^[18]和张琪^[19]专门针对图像阈值分割方法展开了详细研究 ...

2015

0.0

... 吴一全等^[20]对近20年来图像阈值分割方法研究进展进行了总结 ...

2003

0.0

... s定律, X射线穿透物体前后的能量与物质密度满足如下关系^[21,22] ...

... X射线衰减系数 μ主要取决于康普顿散射和光电吸收, 可表达为^[21,22] ...

... 在医用CT中, 通过CT(Hounsfield)数来描述人体组织的衰减系数相对于水的衰减系数的大小, 其定义为^[21,22] ...

2001

0.0

... s定律, X射线穿透物体前后的能量与物质密度满足如下关系^[21,22] ...

... X射线衰减系数 μ主要取决于康普顿散射和光电吸收, 可表达为^[21,22] ...

... 在医用CT中, 通过CT(Hounsfield)数来描述人体组织的衰减系数相对于水的衰减系数的大小, 其定义为^[21,22] ...

... Karacan等^[22]指出煤中主要存在三种基本的结构类型: 由基质孔隙度控制的有机煤基质区域, 由煤体孔隙度控制的有机煤基质区域和含少量孔隙的杂质(矿物)集中区域 ...

2001

0.0

... 为此可以将工业CT的CT数定义为^[23] ...

1979

0.0

... 1 最大类间方差法(Otsu法)最大类间方差法^[24]是1975年日本学者大津展之提出的, 其根据图像预分割的各个部分灰度值之间的总方差最大来确定最佳阈值 ...

1985

0.0

... 2 最大熵法最大熵法^[25]是1985年由Kapur提出的, 其阈值选取准则是图像分割后的各个部分的总熵值最大, 即信息量最大 ...

2004

0.0

周铭. 低信噪比红外小目标图像的分割方法. 激光与红外, 2004, 34(3): 225-228

针对复杂背景条件下的红外图像分割问题,将遗传算法与人类视觉感知原理相结合,探讨了一种新的多阈值图像分割方法即基于遗传算法的图像阈值分割方法,该算法引入了一个自动判别且时空可变的目标背景条件和调整最佳分割区域的步骤,提高了分割算法的质量及鲁棒性,克服了传统阈值方法在图像分割中的局限性.通过计算机仿真实验,验证了该方法的有效性.

... 图6(b)和图7(b)以及表格中数据表明, 最大类间方差法的阈值分割点基本不受裂隙产生的影响, 即图片中微小区域的灰度突变其无法检测到, 而仅仅是对图片中主要成分组成进行检测区分, 这与文献[26]的观点#cod#x0201c ...