基于聚类算法的岩石CT图像分割及量化方法*

张嘉凡，张雪娇，杨更社，刘　慧，张慧梅(.西安科技大学理学院，陕西西安70054; .西安科技大学建筑与土木工程学院，陕西西安70054)



基于聚类算法的岩石CT图像分割及量化方法*

张嘉凡1，张雪娇1，杨更社2，刘慧2，张慧梅1
(1.西安科技大学理学院，陕西西安710054; 2.西安科技大学建筑与土木工程学院，陕西西安710054)

摘要:为研究岩石CT图像分割及量化方法，以识别岩石CT图像中的岩石区、损伤区及背景区为目的，提出了一种聚类算法与数字图像处理技术相结合的方法，该方法根据“物以类聚”的统计原理，按距离相近或相似程度对岩石CT图像中的像素进行标定，从而实现图像分割及量化。结果表明:该方法能够准确地对岩石CT图像中的不同区域进行分割并且实现了对损伤的量化表达;同时，对于结果不确定度影响的初始参数有完全的排异性，从而保证了结果的稳定性;将该算法与阈值分割法进行比较，该算法可避免人为选择阈值导致的误差，从而保证结果的可靠性。

关键词:聚类分析;阈值分割;岩石CT图像

0　引言

随着科学技术尤其是计算机技术的发展，CT扫描或探测技术的应用日趋成熟并逐渐成为各行业发展的热点。以岩石损伤探测及量化工作为对象，CT技术目前已成为该领域最为先进的探伤及分析手段之一。毛林涛［1］等通过差值图像的统计特征，分析了煤岩内部裂隙与宏观变形之间的关系;杨更社［2］利用数字图像增强技术实现了CT图像的伪彩色增强，并根据灰度直方图技术分析了岩石损伤规律;刘慧［3］通过遗传算法与最大类间方差相结合，自动选取图像的最优阈值，完成了冻结岩石CT图像的三值化分割;张青成［4］等基于图像分割技术，对煤岩CT图像的灰度级别设定不同的二值化阀值，得到不同阀值下的孔隙面积变化曲线，并提出以拐点处对应的阀值作为裂隙图像二值化阀值时效果最佳;马天寿，陈平［5］采用最大自动取阈值方法对页岩水化CT图像进行分割，得到了损伤变量与浸泡时间的关系。

事实上，图像差值的可操性很低，它要求不同状态的扫描位置不能发生改变，否则很容易造成图像差值的结果毫无实际意义;伪彩色增强的方法虽然能很好地增强岩石裂纹识别的视觉效果，但无法做出定量评价;目前阈值分割法是最常用的图像分割方法，对于不同类的物体灰度值或其他特征值相差很大时，它能很有效地对图像进行分割，其关键技术是阈值的选取。最大类间方差(简称Otsu)法是常用的阈值图像分割方法之一，但该技术不能良好解决由于图像峰型复杂、灰度值差异不明显所引起的图像分割效率低甚至错误的问题。

要解决上述问题，主要涉及2种技术，其一为提高工业CT的分辨率和操作精度以提高图像的质量，这依赖于工业CT扫描技术的发展与进步;其二为结合数字图像处理技术与数学手段，科学合理地建立图像分割指标。因此，文中将采用数字图像处理技术与聚类分析相结合的方法，对岩石CT图像中的岩石区、损伤区及背景区进行分割识别，并且量化损伤区，为进一步研究岩石内部裂隙的演化规律提供了信息依据。

1　CT扫描原理

如图1所示，将被测物体放置于射线源与探测器之间，射线源发出射线并穿透被测物体，这必然会导致射线速度、强度、频率等物理量数值的变化，而探测器则会检测并记录这些数据的变化，从而间接地反映被测物体内部的变化情况。就同一波长的X射线而言，密度越大的物质，其吸收X射线的能力越强。不同物质对同一波长X射线的吸收能力不同，物质密度越大，对X射线的吸收能力越强。根据这一物理原理，CT数与对应的岩石密度成正比，CT图像每一像素点的值即为CT数，定义为［6］

式中u为被测物体的X射线线性吸收系数; um为水的X射线线性吸收系数。

由于岩石材料的组成和结构不均匀，导致各部位密度不均匀。而密度又与X射线吸收系数成正比，因此CT图像就可看做是岩石扫描断面密度分布图。岩石内部一定区域内孔洞、裂纹的活动必然引起该区域密度的变化也可反应本区域孔洞、裂纹活动的集合效应。

图1　CT扫描原理图Fig．1 CT scanning principle diagram1射线源　2被测物　3接收器　4旋转台

2　聚类算法分割CT图像

2. 1聚类分析原理

聚类分析(cluster analysis)，又叫集群分析，是一种“物以类聚”的统计描述方法，其实质是寻找一些能客观反映研究对象之间亲疏关系的统计量，然后根据这种统计量把研究对象按距离相近或者相似的原则分成若干类。聚类算法无需对样本进行训练，是一种无监督的统计方法，通过迭代的方式对图像进行分类。因此从某种意义上说，聚类是一种自我训练的分类，避免了人为选择阈值导致的误差，是一种行之有效的方法。

文中所使用的算法为全局K均值算法，它产生的结果不依赖于任何初始的参数设置，解决了由于随机初始化聚类中心导致的陷入局部最小的可能［7］。算法的基本思想是:先对当前的每一类求均值，然后按新生成的均值对对象进行分类(将对象归入均值最近的类)，对新生成的类再迭代执行前面的步骤。具体步骤如下［8］

1)随机选取K个对象作为初始聚类的聚类中心;

2)计算对象与各个类的聚类中心的距离，将对象划分到距离其最近的一类;

3)重新计算每个新类的均值，作为下一次聚类的中心;

4)若类的聚类中心不再变化，则返回聚类结果，否则转步骤2)

2. 2聚类算法对岩石CT图像的分割

在岩石CT扫描图像中，对应不同密度的对象就有不同的CT值，表现为不同的灰度值［9］。在一幅岩石CT扫描图像中，其灰度共有256个级别，0代表黑色，255代表白色，在0到255之间则代表着不同程度的灰色地带。图2为岩石下层扫描层面的CT检测结果，图中物质的密度越大表现在图像的亮度越低，反之亮度越高。岩体颗粒的密度最大，表现在图像中为灰色或灰黑色，而裂纹、孔洞的密度较岩石颗粒小，表现在图像中为灰白色。从图2可看出图像整体偏暗，相对于边缘而言图像中心较亮，岩石颗粒的边界与裂纹和空隙的分界不明显，图像的对比度低。

图2　岩石CT扫描原图Fig．2 Rock CT scanning diagram

2. 2. 1相似性量度

将聚类算法应用于岩石CT图像的裂纹识别中，首先要求有一个“相似性”的量度［10］，通常情况下采用绝对距离或欧氏距离作为相似性指标，但该方法忽略了像素点与周围像素点CT数之间的大小，也就是没有考虑到岩石内部裂纹几何上的模糊性，因此无法将裂纹区分开。文中提出一种新的距离指标，该指标是绝对距离的改进，其具体描述如下

式中dij为第i行j列的像素点与聚类中心之间的距离指标; f(t)为经过t次循环后聚类中心值;

为以f(i，j)为中心的九宫格中所有像素点的均值作为该点后CT值。

2. 2. 2算法应用

算法应用具体分为以下几个步骤:

第一步初始化聚类中心。根据岩石CT图像基本特征，认为CT图像中的像素有3种相对纯粹的物质，分别是岩石区、损伤区以及背景区，由于CT图像的CT值与物质的密度呈正比，因此可以判断这种物质的密度关系是

岩体颗粒＞损伤＞背景

图3为原图像图2的灰度直方图，该直方图形成了3个峰，根据灰度值与密度的关系，可以初步判断代表岩石颗粒的灰度范围是0到100左右，代表损伤的灰度范围是100到150左右，而背景则全是白色，即灰度值为255所代表的像素。于是初始聚类中心设可设为其范围的中心值

c(1) =50，c(2) =130，c(3) =255．

将此次初始聚类中心选取方案命名为方案1．

图3　灰度直方图分析图Fig．3 Grey histogram analysis diagram

第二步是进行初始聚类。首先计算各像素点与初始聚类中心的距离，然后寻找与3个初始聚类中心距离最小的那个聚类，并将该像素分配给该类，最后对3个聚类中的所有像素值求平均值，作为下一次聚类的聚类中心。

第三步是进行循环聚类，设定收敛条件:相邻两次的聚类中心的绝对差值小于0. 5．重复第二步直到达到收敛条件。

文中利用Matlab编程实现了聚类算法对岩石CT图像的分割，分割结果如图4所示。可以看到，聚类分割后，岩石区与损伤区黑白分明，边界清晰，极为细小的孔隙也能很好的识别。

由于聚类算法的优越性，编程时可在聚类的过程中添加相应的变量，分别对每次聚类的终止聚类中心、每个类的像素个数、每类像素个数在整幅图像中所占的百分比以及完成聚类所需迭代的次数进行统计，统计结果见表1．此次聚类迭代次数为10次，与终止聚类中心的欧式距离为94，迭代终止时的聚类中心为c(1) =28，c(2) =90，c(3)=255，背景在图像中占21. 02%、损伤在图像中占10. 18%，岩石在图像中占68. 80%．

图4　聚类前后的对比图Fig．4　Comparsion before and after clustering

表1　聚类结果Tab．1 Clustering result

图5　方案2，3的聚类结果Fig．5　Clustering result based on method 2 and 3

2. 2. 3不同初始聚类中心聚类结果比较

为了验证初始聚类中心的选取是否对聚类结果造成影响，文中另外选取了2种初始聚类中心的取值方案，分别命名为方案2，方案3．得到的聚类结果，如图5所示。观察图像可知，聚类结果几乎没有差异，说明该算法的稳定性强。分析这3种方案的聚类量化数据表，见表1，表2，表3，发现3种方案的终止聚类中心没有改变，都是28，90，255;聚类的各类所占的百分比也几乎没有变化;有变化的是聚类完成所需的迭代次数，可以发现初始聚类中心与实际聚类中心距离越远，聚类完成需要迭代的次数越多，反之亦然。因此，为了提高程序的运行效率，建议初始聚类中心的选值参照图像的直方图进行选取，如果图像存在多峰的情况，则选取阈值为图像的峰值;如果图像为单峰图，则选取阈值为整幅图像灰度值的中间值，从而减少程序运行迭代的次数。

表2　方案2聚类结果Tab．2 Clustering result based on method 2

表3　方案3聚类结果Tab．3 Clustering result based on method 3

3　与阈值分割法对比

阈值分割是一种基于区域的分割算法［11］。基本原理是采用不同的特征阈值对图像像素点分进行分类。该方法的机理在于默认相似灰度值对应于同一类对象，并以灰度直方图中不同峰值加以区别。选取的阈值应位于2个峰之间的谷，从而将各个峰分开。

基于阈值分割算法的空隙及裂纹提取准则如下

式中f(i，j)为第i行，第j列的灰度值;λ为空隙或裂纹的阈值。

可以说，只要选取合适的阈值就能将岩石CT图像进行有效的分割，进而识别出岩石内部孔洞及裂纹的分布情况。但如何选取合适的阈值需要采取反复的实验方法，对阈值进行调整。

图6为阈值分别为λ1= 100，λ2= 200;λ1= 120，λ2=240的裂纹提取结果，可以看到阈值分割提取裂纹简单易行，但是要确定合理的阈值需要进行大量的阈值实验，且算法存在较大的差异，不同的阈值对提取结果影响较大。

图6　不同阈值分割的结果Fig．6　Different threshold segmentation results

4　结论

1)聚类分析是一种基于统计学的分割算法，将其应用于岩石CT图像的分割，能够得到很好的分割效果，从而识别出岩石CT图片中的损伤区、岩石区及背景区，并且实现了岩石图像的量化表达，得到岩石CT图像中岩石区占68. 8%，损伤区占10. 18%，剩下的背景区则占21. 02%;

2)初始聚类中心的选值会影响聚类的迭代次数，但对聚类的结果不会造成影响。文中使用3种初始聚类中心的选取方案进行聚类，3种方案的终止聚类中心都是(28，90，255)，与初始聚类中心的距离分别是94，13，163，对应所用的迭代次数是10次，2次，28次，由此可得，初始聚类中心离终止聚类中心越近，迭代次数越少，反之亦然;

3)阈值分割算法和聚类分析算法均可应用于岩石CT图像。前者受阈值影响较大，确定合理的阈值往往需要大量的试验;后者只要设定收敛值，即可得到很好的结果，且可以避免人为因素的影响，效率更高，算法更可靠。

参考文献References

［1］毛林涛．煤样力学特性与内部裂隙演化关系CT实验研究［J］．辽宁工程技术大学学报，2010，29(3) : 408 －411. MAO Lin-tao．Experimental study on relationship between coal mechanical characteristics and interior crack evolution using CT technique［J］．Journal of Liaoning Technical University，2010，29(3) : 408－411.

［2］杨更社．基于CT图像处理技术的岩石损伤特性研究［J］．煤炭学报，2007，32(5) : 463－468. YANG Geng-she．Study on the rock damage characteristics based on the technique of CT image processing［J］．Journal of China Coal Society，2010，29(3) : 408－411.

［3］刘慧．基于CT图像处理技术的冻结岩石细观结构及损伤力学特性研究［D］．西安:西安科技大学，2013. LIU Hui．Study on meso-structure and damage mechanical characteristics of frozen rock based on CT image processing［D］．Xi’an: Xi’an University of Science and Technology，2013.

［4］张青成，王万富．煤岩CT图像二值化阀值选取及三维重构技术研究［J］．CT理论与应用研究，2014，23 (1) : 45－50. ZHANG Qing-cheng，WANG Wan-fu．CT image binarization threshold selection of coal and 3D reconstruction technology research［J］．CT Theory and Applications，2014，23(1) : 45－50.

［5］马天寿，陈平．基于CT扫描技术研究页岩水化细观损伤特性［J］．石油勘探与开发，2014，41(2) : 227 －233. MA Tian-shou，CHEN Ping．Study of meso-damage characteristics of shale hydration based on CT scanning technology［J］．Petroleum Exploration and Development，2014，41(2) : 227－233.

［6］张朝宗．工业CT技术和原理［M］．北京:科学出版社，2009. ZHANG Chao-zhong．Industrial CT technique and principle［M］．Beijing: Science Press，2009.

［7］张鑫野．基于聚类分析的图像分割方法研究［D］．大连:大连海事大学，2012. ZHANG Xin-ye．Research on image segmentation method based on clustering analysis［D］．Dalian: Maritime Affairs University of Dalian，2012.

［8］王文平．聚类分析及其在图像分割中的应用［D］．济南:山东师范大学，2007. WANG Wen-ping．Application of clustering analysis on image segmentation［D］．Jinan: Shandong Normal University，2007.

［9］陈超．MATLAB应用实例精讲——图像处理与GUI设计篇［M］．北京:电子工业出版社，2011. CHEN Chao．MATLAB application examples-Image processing and GUI design［M］．Beijing: Electronic Industry Press，2011.

［10］LICHARD A Jornson，DEAN W Wekeen．Practical pluralistic statistical analysis［M］．Beijing: Tsinghua University Press，2008.

［11］郝景宏，姜袁，梅世强．基于X射线CT的混凝土内部裂纹识别研究［J］．混凝土，2010(10) : 44－47. HAO Jing-hong，JIANG Yuan，MEI Shi-qiang．Identification of crack in concrete interior based on X－ray CT ［J］．Concert，2010(10) : 44－47.

A method of rock CT image segmentation and quantification based on clustering algorithm

ZHANG Jia-fan1，ZHANG Xue-jiao1，YANG Geng-she2，LIU Hui2，ZHANG Hui-mei1
(1. College of Science，Xi’an University of Science and Technology，Xi’an 710054，China; 2. College of Civil and Architectural Engineering，Xi’an University of Science and Technology，Xi’an 710054，China)

Abstract:In order to study rock CT image segmentation and quantification，with the purpose of identifying the rock area，damage area and background area in rock CT image，a method combining clustering algorithm with digital image processing techniques was given．The application results show that: The method can accurately segment the different areas of rock CT image and quantify its percentage; Meanwhile，the initial parameters that can influence the result uncertainly have complementarily mutual exclusiveness completely，ensuring the stability of the result; Compared the algorithm with threshold segmentation method，the algorithm can avoid the error caused by human choice threshold，which ensured the reliability of the results．

Key words:clustering analysis; threshold segmentation; rock CT damage

通讯作者:张慧梅(1967－)，女，山西大同人，工学博士，教授，E-mail: zhanghuimei168@163．com

基金项目:国家自然科学基金项目(11172232，41272340) ;教育部新世纪人才支持计划项目(NECT－12－1044) ;陕西省教育厅专项基金项目(11JK0542)

*收稿日期:2015－11－20责任编辑:李克永

DOI:10．13800/j．cnki．xakjdxxb．2016．0204

文章编号:1672－9315(2016) 02－0171－05

中图分类号:TU 455

文献标志码:A