張嘉凡，張雪嬌，楊更社，劉　慧，張慧梅(.西安科技大學理學院，陜西西安70054; .西安科技大學建筑與土木工程學院，陜西西安70054)

?

基于聚類算法的巖石CT圖像分割及量化方法*

張嘉凡1，張雪嬌1，楊更社2，劉慧2，張慧梅1
(1.西安科技大學理學院，陜西西安710054; 2.西安科技大學建筑與土木工程學院，陜西西安710054)

摘要:為研究巖石CT圖像分割及量化方法，以識別巖石CT圖像中的巖石區(qū)、損傷區(qū)及背景區(qū)為目的，提出了一種聚類算法與數字圖像處理技術相結合的方法，該方法根據“物以類聚”的統(tǒng)計原理，按距離相近或相似程度對巖石CT圖像中的像素進行標定，從而實現圖像分割及量化。結果表明:該方法能夠準確地對巖石CT圖像中的不同區(qū)域進行分割并且實現了對損傷的量化表達;同時，對于結果不確定度影響的初始參數有完全的排異性，從而保證了結果的穩(wěn)定性;將該算法與閾值分割法進行比較，該算法可避免人為選擇閾值導致的誤差，從而保證結果的可靠性。

關鍵詞:聚類分析;閾值分割;巖石CT圖像

0　引言

隨著科學技術尤其是計算機技術的發(fā)展，CT掃描或探測技術的應用日趨成熟并逐漸成為各行業(yè)發(fā)展的熱點。以巖石損傷探測及量化工作為對象，CT技術目前已成為該領域最為先進的探傷及分析手段之一。毛林濤［1］等通過差值圖像的統(tǒng)計特征，分析了煤巖內部裂隙與宏觀變形之間的關系;楊更社［2］利用數字圖像增強技術實現了CT圖像的偽彩色增強，并根據灰度直方圖技術分析了巖石損傷規(guī)律;劉慧［3］通過遺傳算法與最大類間方差相結合，自動選取圖像的最優(yōu)閾值，完成了凍結巖石CT圖像的三值化分割;張青成［4］等基于圖像分割技術，對煤巖CT圖像的灰度級別設定不同的二值化閥值，得到不同閥值下的孔隙面積變化曲線，并提出以拐點處對應的閥值作為裂隙圖像二值化閥值時效果最佳;馬天壽，陳平［5］采用最大自動取閾值方法對頁巖水化CT圖像進行分割，得到了損傷變量與浸泡時間的關系。

事實上，圖像差值的可操性很低，它要求不同狀態(tài)的掃描位置不能發(fā)生改變，否則很容易造成圖像差值的結果毫無實際意義;偽彩色增強的方法雖然能很好地增強巖石裂紋識別的視覺效果，但無法做出定量評價;目前閾值分割法是最常用的圖像分割方法，對于不同類的物體灰度值或其他特征值相差很大時，它能很有效地對圖像進行分割，其關鍵技術是閾值的選取。最大類間方差(簡稱Otsu)法是常用的閾值圖像分割方法之一，但該技術不能良好解決由于圖像峰型復雜、灰度值差異不明顯所引起的圖像分割效率低甚至錯誤的問題。

要解決上述問題，主要涉及2種技術，其一為提高工業(yè)CT的分辨率和操作精度以提高圖像的質量，這依賴于工業(yè)CT掃描技術的發(fā)展與進步;其二為結合數字圖像處理技術與數學手段，科學合理地建立圖像分割指標。因此，文中將采用數字圖像處理技術與聚類分析相結合的方法，對巖石CT圖像中的巖石區(qū)、損傷區(qū)及背景區(qū)進行分割識別，并且量化損傷區(qū)，為進一步研究巖石內部裂隙的演化規(guī)律提供了信息依據。

1　CT掃描原理

如圖1所示，將被測物體放置于射線源與探測器之間，射線源發(fā)出射線并穿透被測物體，這必然會導致射線速度、強度、頻率等物理量數值的變化，而探測器則會檢測并記錄這些數據的變化，從而間接地反映被測物體內部的變化情況。就同一波長的X射線而言，密度越大的物質，其吸收X射線的能力越強。不同物質對同一波長X射線的吸收能力不同，物質密度越大，對X射線的吸收能力越強。根據這一物理原理，CT數與對應的巖石密度成正比，CT圖像每一像素點的值即為CT數，定義為［6］

式中u為被測物體的X射線線性吸收系數; um為水的X射線線性吸收系數。

由于巖石材料的組成和結構不均勻，導致各部位密度不均勻。而密度又與X射線吸收系數成正比，因此CT圖像就可看做是巖石掃描斷面密度分布圖。巖石內部一定區(qū)域內孔洞、裂紋的活動必然引起該區(qū)域密度的變化也可反應本區(qū)域孔洞、裂紋活動的集合效應。

圖1　CT掃描原理圖Fig．1 CT scanning principle diagram1射線源　2被測物　3接收器　4旋轉臺

2　聚類算法分割CT圖像

2. 1聚類分析原理

聚類分析(cluster analysis)，又叫集群分析，是一種“物以類聚”的統(tǒng)計描述方法，其實質是尋找一些能客觀反映研究對象之間親疏關系的統(tǒng)計量，然后根據這種統(tǒng)計量把研究對象按距離相近或者相似的原則分成若干類。聚類算法無需對樣本進行訓練，是一種無監(jiān)督的統(tǒng)計方法，通過迭代的方式對圖像進行分類。因此從某種意義上說，聚類是一種自我訓練的分類，避免了人為選擇閾值導致的誤差，是一種行之有效的方法。

文中所使用的算法為全局K均值算法，它產生的結果不依賴于任何初始的參數設置，解決了由于隨機初始化聚類中心導致的陷入局部最小的可能［7］。算法的基本思想是:先對當前的每一類求均值，然后按新生成的均值對對象進行分類(將對象歸入均值最近的類)，對新生成的類再迭代執(zhí)行前面的步驟。具體步驟如下［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數之間的大小，也就是沒有考慮到巖石內部裂紋幾何上的模糊性，因此無法將裂紋區(qū)分開。文中提出一種新的距離指標，該指標是絕對距離的改進，其具體描述如下

式中dij為第i行j列的像素點與聚類中心之間的距離指標; f(t)為經過t次循環(huán)后聚類中心值;

為以f(i，j)為中心的九宮格中所有像素點的均值作為該點后CT值。

2. 2. 2算法應用

算法應用具體分為以下幾個步驟:

第一步初始化聚類中心。根據巖石CT圖像基本特征，認為CT圖像中的像素有3種相對純粹的物質，分別是巖石區(qū)、損傷區(qū)以及背景區(qū)，由于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個聚類中的所有像素值求平均值，作為下一次聚類的聚類中心。

第三步是進行循環(huán)聚類，設定收斂條件:相鄰兩次的聚類中心的絕對差值小于0. 5．重復第二步直到達到收斂條件。

文中利用Matlab編程實現了聚類算法對巖石CT圖像的分割，分割結果如圖4所示?？梢钥吹?，聚類分割后，巖石區(qū)與損傷區(qū)黑白分明，邊界清晰，極為細小的孔隙也能很好的識別。

由于聚類算法的優(yōu)越性，編程時可在聚類的過程中添加相應的變量，分別對每次聚類的終止聚類中心、每個類的像素個數、每類像素個數在整幅圖像中所占的百分比以及完成聚類所需迭代的次數進行統(tǒng)計，統(tǒng)計結果見表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所示。觀察圖像可知，聚類結果幾乎沒有差異，說明該算法的穩(wěn)定性強。分析這3種方案的聚類量化數據表，見表1，表2，表3，發(fā)現3種方案的終止聚類中心沒有改變，都是28，90，255;聚類的各類所占的百分比也幾乎沒有變化;有變化的是聚類完成所需的迭代次數，可以發(fā)現初始聚類中心與實際聚類中心距離越遠，聚類完成需要迭代的次數越多，反之亦然。因此，為了提高程序的運行效率，建議初始聚類中心的選值參照圖像的直方圖進行選取，如果圖像存在多峰的情況，則選取閾值為圖像的峰值;如果圖像為單峰圖，則選取閾值為整幅圖像灰度值的中間值，從而減少程序運行迭代的次數。

表2　方案2聚類結果Tab．2 Clustering result based on method 2

表3　方案3聚類結果Tab．3 Clustering result based on method 3

3　與閾值分割法對比

閾值分割是一種基于區(qū)域的分割算法［11］?；驹硎遣捎貌煌奶卣鏖撝祵D像像素點分進行分類。該方法的機理在于默認相似灰度值對應于同一類對象，并以灰度直方圖中不同峰值加以區(qū)別。選取的閾值應位于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)聚類分析是一種基于統(tǒng)計學的分割算法，將其應用于巖石CT圖像的分割，能夠得到很好的分割效果，從而識別出巖石CT圖片中的損傷區(qū)、巖石區(qū)及背景區(qū)，并且實現了巖石圖像的量化表達，得到巖石CT圖像中巖石區(qū)占68. 8%，損傷區(qū)占10. 18%，剩下的背景區(qū)則占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］．石油勘探與開發(fā)，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］張朝宗．工業(yè)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］．北京:電子工業(yè)出版社，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