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

?

基于聚類(lèi)算法的巖石CT圖像分割及量化方法*

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

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

關(guān)鍵詞:聚類(lèi)分析;閾值分割;巖石CT圖像

0　引言

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

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

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

1　CT掃描原理

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

式中u為被測(cè)物體的X射線(xiàn)線(xiàn)性吸收系數(shù); um為水的X射線(xiàn)線(xiàn)性吸收系數(shù)。

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

圖1　CT掃描原理圖Fig．1 CT scanning principle diagram1射線(xiàn)源　2被測(cè)物　3接收器　4旋轉(zhuǎn)臺(tái)

2　聚類(lèi)算法分割CT圖像

2. 1聚類(lèi)分析原理

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

文中所使用的算法為全局K均值算法，它產(chǎn)生的結(jié)果不依賴(lài)于任何初始的參數(shù)設(shè)置，解決了由于隨機(jī)初始化聚類(lèi)中心導(dǎo)致的陷入局部最小的可能［7］。算法的基本思想是:先對(duì)當(dāng)前的每一類(lèi)求均值，然后按新生成的均值對(duì)對(duì)象進(jìn)行分類(lèi)(將對(duì)象歸入均值最近的類(lèi))，對(duì)新生成的類(lèi)再迭代執(zhí)行前面的步驟。具體步驟如下［8］

1)隨機(jī)選取K個(gè)對(duì)象作為初始聚類(lèi)的聚類(lèi)中心;

2)計(jì)算對(duì)象與各個(gè)類(lèi)的聚類(lèi)中心的距離，將對(duì)象劃分到距離其最近的一類(lèi);

3)重新計(jì)算每個(gè)新類(lèi)的均值，作為下一次聚類(lèi)的中心;

4)若類(lèi)的聚類(lèi)中心不再變化，則返回聚類(lèi)結(jié)果，否則轉(zhuǎn)步驟2)

2. 2聚類(lèi)算法對(duì)巖石CT圖像的分割

在巖石CT掃描圖像中，對(duì)應(yīng)不同密度的對(duì)象就有不同的CT值，表現(xiàn)為不同的灰度值［9］。在一幅巖石CT掃描圖像中，其灰度共有256個(gè)級(jí)別，0代表黑色，255代表白色，在0到255之間則代表著不同程度的灰色地帶。圖2為巖石下層掃描層面的CT檢測(cè)結(jié)果，圖中物質(zhì)的密度越大表現(xiàn)在圖像的亮度越低，反之亮度越高。巖體顆粒的密度最大，表現(xiàn)在圖像中為灰色或灰黑色，而裂紋、孔洞的密度較巖石顆粒小，表現(xiàn)在圖像中為灰白色。從圖2可看出圖像整體偏暗，相對(duì)于邊緣而言圖像中心較亮，巖石顆粒的邊界與裂紋和空隙的分界不明顯，圖像的對(duì)比度低。

圖2　巖石CT掃描原圖Fig．2 Rock CT scanning diagram

2. 2. 1相似性量度

將聚類(lèi)算法應(yīng)用于巖石CT圖像的裂紋識(shí)別中，首先要求有一個(gè)“相似性”的量度［10］，通常情況下采用絕對(duì)距離或歐氏距離作為相似性指標(biāo)，但該方法忽略了像素點(diǎn)與周?chē)袼攸c(diǎn)CT數(shù)之間的大小，也就是沒(méi)有考慮到巖石內(nèi)部裂紋幾何上的模糊性，因此無(wú)法將裂紋區(qū)分開(kāi)。文中提出一種新的距離指標(biāo)，該指標(biāo)是絕對(duì)距離的改進(jìn)，其具體描述如下

式中dij為第i行j列的像素點(diǎn)與聚類(lèi)中心之間的距離指標(biāo); f(t)為經(jīng)過(guò)t次循環(huán)后聚類(lèi)中心值;

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

2. 2. 2算法應(yīng)用

算法應(yīng)用具體分為以下幾個(gè)步驟:

第一步初始化聚類(lèi)中心。根據(jù)巖石CT圖像基本特征，認(rèn)為CT圖像中的像素有3種相對(duì)純粹的物質(zhì)，分別是巖石區(qū)、損傷區(qū)以及背景區(qū)，由于CT圖像的CT值與物質(zhì)的密度呈正比，因此可以判斷這種物質(zhì)的密度關(guān)系是

巖體顆粒＞損傷＞背景

圖3為原圖像圖2的灰度直方圖，該直方圖形成了3個(gè)峰，根據(jù)灰度值與密度的關(guān)系，可以初步判斷代表巖石顆粒的灰度范圍是0到100左右，代表?yè)p傷的灰度范圍是100到150左右，而背景則全是白色，即灰度值為255所代表的像素。于是初始聚類(lèi)中心設(shè)可設(shè)為其范圍的中心值

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

將此次初始聚類(lèi)中心選取方案命名為方案1．

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

第二步是進(jìn)行初始聚類(lèi)。首先計(jì)算各像素點(diǎn)與初始聚類(lèi)中心的距離，然后尋找與3個(gè)初始聚類(lèi)中心距離最小的那個(gè)聚類(lèi)，并將該像素分配給該類(lèi)，最后對(duì)3個(gè)聚類(lèi)中的所有像素值求平均值，作為下一次聚類(lèi)的聚類(lèi)中心。

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

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

由于聚類(lèi)算法的優(yōu)越性，編程時(shí)可在聚類(lèi)的過(guò)程中添加相應(yīng)的變量，分別對(duì)每次聚類(lèi)的終止聚類(lèi)中心、每個(gè)類(lèi)的像素個(gè)數(shù)、每類(lèi)像素個(gè)數(shù)在整幅圖像中所占的百分比以及完成聚類(lèi)所需迭代的次數(shù)進(jìn)行統(tǒng)計(jì)，統(tǒng)計(jì)結(jié)果見(jiàn)表1．此次聚類(lèi)迭代次數(shù)為10次，與終止聚類(lèi)中心的歐式距離為94，迭代終止時(shí)的聚類(lèi)中心為c(1) =28，c(2) =90，c(3)=255，背景在圖像中占21. 02%、損傷在圖像中占10. 18%，巖石在圖像中占68. 80%．

圖4　聚類(lèi)前后的對(duì)比圖Fig．4　Comparsion before and after clustering

表1　聚類(lèi)結(jié)果Tab．1 Clustering result

圖5　方案2，3的聚類(lèi)結(jié)果Fig．5　Clustering result based on method 2 and 3

2. 2. 3不同初始聚類(lèi)中心聚類(lèi)結(jié)果比較

為了驗(yàn)證初始聚類(lèi)中心的選取是否對(duì)聚類(lèi)結(jié)果造成影響，文中另外選取了2種初始聚類(lèi)中心的取值方案，分別命名為方案2，方案3．得到的聚類(lèi)結(jié)果，如圖5所示。觀(guān)察圖像可知，聚類(lèi)結(jié)果幾乎沒(méi)有差異，說(shuō)明該算法的穩(wěn)定性強(qiáng)。分析這3種方案的聚類(lèi)量化數(shù)據(jù)表，見(jiàn)表1，表2，表3，發(fā)現(xiàn)3種方案的終止聚類(lèi)中心沒(méi)有改變，都是28，90，255;聚類(lèi)的各類(lèi)所占的百分比也幾乎沒(méi)有變化;有變化的是聚類(lèi)完成所需的迭代次數(shù)，可以發(fā)現(xiàn)初始聚類(lèi)中心與實(shí)際聚類(lèi)中心距離越遠(yuǎn)，聚類(lèi)完成需要迭代的次數(shù)越多，反之亦然。因此，為了提高程序的運(yùn)行效率，建議初始聚類(lèi)中心的選值參照?qǐng)D像的直方圖進(jìn)行選取，如果圖像存在多峰的情況，則選取閾值為圖像的峰值;如果圖像為單峰圖，則選取閾值為整幅圖像灰度值的中間值，從而減少程序運(yùn)行迭代的次數(shù)。

表2　方案2聚類(lèi)結(jié)果Tab．2 Clustering result based on method 2

表3　方案3聚類(lèi)結(jié)果Tab．3 Clustering result based on method 3

3　與閾值分割法對(duì)比

閾值分割是一種基于區(qū)域的分割算法［11］?；驹硎遣捎貌煌奶卣鏖撝祵?duì)圖像像素點(diǎn)分進(jìn)行分類(lèi)。該方法的機(jī)理在于默認(rèn)相似灰度值對(duì)應(yīng)于同一類(lèi)對(duì)象，并以灰度直方圖中不同峰值加以區(qū)別。選取的閾值應(yīng)位于2個(gè)峰之間的谷，從而將各個(gè)峰分開(kāi)。

基于閾值分割算法的空隙及裂紋提取準(zhǔn)則如下

式中f(i，j)為第i行，第j列的灰度值;λ為空隙或裂紋的閾值。

可以說(shuō)，只要選取合適的閾值就能將巖石CT圖像進(jìn)行有效的分割，進(jìn)而識(shí)別出巖石內(nèi)部孔洞及裂紋的分布情況。但如何選取合適的閾值需要采取反復(fù)的實(shí)驗(yàn)方法，對(duì)閾值進(jìn)行調(diào)整。

圖6為閾值分別為λ1= 100，λ2= 200;λ1= 120，λ2=240的裂紋提取結(jié)果，可以看到閾值分割提取裂紋簡(jiǎn)單易行，但是要確定合理的閾值需要進(jìn)行大量的閾值實(shí)驗(yàn)，且算法存在較大的差異，不同的閾值對(duì)提取結(jié)果影響較大。

圖6　不同閾值分割的結(jié)果Fig．6　Different threshold segmentation results

4　結(jié)論

1)聚類(lèi)分析是一種基于統(tǒng)計(jì)學(xué)的分割算法，將其應(yīng)用于巖石CT圖像的分割，能夠得到很好的分割效果，從而識(shí)別出巖石CT圖片中的損傷區(qū)、巖石區(qū)及背景區(qū)，并且實(shí)現(xiàn)了巖石圖像的量化表達(dá)，得到巖石CT圖像中巖石區(qū)占68. 8%，損傷區(qū)占10. 18%，剩下的背景區(qū)則占21. 02%;

2)初始聚類(lèi)中心的選值會(huì)影響聚類(lèi)的迭代次數(shù)，但對(duì)聚類(lèi)的結(jié)果不會(huì)造成影響。文中使用3種初始聚類(lèi)中心的選取方案進(jìn)行聚類(lèi)，3種方案的終止聚類(lèi)中心都是(28，90，255)，與初始聚類(lèi)中心的距離分別是94，13，163，對(duì)應(yīng)所用的迭代次數(shù)是10次，2次，28次，由此可得，初始聚類(lèi)中心離終止聚類(lèi)中心越近，迭代次數(shù)越少，反之亦然;

3)閾值分割算法和聚類(lèi)分析算法均可應(yīng)用于巖石CT圖像。前者受閾值影響較大，確定合理的閾值往往需要大量的試驗(yàn);后者只要設(shè)定收斂值，即可得到很好的結(jié)果，且可以避免人為因素的影響，效率更高，算法更可靠。

參考文獻(xiàn)References

［1］毛林濤．煤樣力學(xué)特性與內(nèi)部裂隙演化關(guān)系CT實(shí)驗(yàn)研究［J］．遼寧工程技術(shù)大學(xué)學(xué)報(bào)，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圖像處理技術(shù)的巖石損傷特性研究［J］．煤炭學(xué)報(bào)，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圖像處理技術(shù)的凍結(jié)巖石細(xì)觀(guān)結(jié)構(gòu)及損傷力學(xué)特性研究［D］．西安:西安科技大學(xué)，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］張青成，王萬(wàn)富．煤巖CT圖像二值化閥值選取及三維重構(gòu)技術(shù)研究［J］．CT理論與應(yīng)用研究，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掃描技術(shù)研究頁(yè)巖水化細(xì)觀(guān)損傷特性［J］．石油勘探與開(kāi)發(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技術(shù)和原理［M］．北京:科學(xué)出版社，2009. ZHANG Chao-zhong．Industrial CT technique and principle［M］．Beijing: Science Press，2009.

［7］張?chǎng)我埃诰垲?lèi)分析的圖像分割方法研究［D］．大連:大連海事大學(xué)，2012. ZHANG Xin-ye．Research on image segmentation method based on clustering analysis［D］．Dalian: Maritime Affairs University of Dalian，2012.

［8］王文平．聚類(lèi)分析及其在圖像分割中的應(yīng)用［D］．濟(jì)南:山東師范大學(xué)，2007. WANG Wen-ping．Application of clustering analysis on image segmentation［D］．Jinan: Shandong Normal University，2007.

［9］陳超．MATLAB應(yīng)用實(shí)例精講——圖像處理與GUI設(shè)計(jì)篇［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］郝景宏，姜袁，梅世強(qiáng)．基于X射線(xiàn)CT的混凝土內(nèi)部裂紋識(shí)別研究［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－)，女，山西大同人，工學(xué)博士，教授，E-mail: zhanghuimei168@163．com

基金項(xiàng)目:國(guó)家自然科學(xué)基金項(xiàng)目(11172232，41272340) ;教育部新世紀(jì)人才支持計(jì)劃項(xiàng)目(NECT－12－1044) ;陜西省教育廳專(zhuān)項(xiàng)基金項(xiàng)目(11JK0542)

*收稿日期:2015－11－20責(zé)任編輯:李克永

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

文章編號(hào):1672－9315(2016) 02－0171－05

中圖分類(lèi)號(hào):TU 455

文獻(xiàn)標(biāo)志碼:A