李周，滕奇志＊，張余強(qiáng)

（1.四川大學(xué)電子信息學(xué)院圖像信息研究所，成都610065；2.成都西圖科技公司，成都 610024）

0 引言

在石油地質(zhì)行業(yè)中，確定巖石中各礦物成分比例是薄片礦物鑒定的重要內(nèi)容[1]，對(duì)巖石顆粒進(jìn)行準(zhǔn)確分割則是確定礦物成分的前提。在巖石薄片圖像中，顆粒表面紋理復(fù)雜，形狀各異，顆粒間差異不明顯，導(dǎo)致巖石顆粒圖像自動(dòng)分割的難度較大。

不少學(xué)者對(duì)巖石顆粒圖像的自動(dòng)分割進(jìn)行了研究，如文獻(xiàn)[2]中B.Obara將巖石圖像轉(zhuǎn)換到不同的顏色空間，然后再進(jìn)行閾值分割。文獻(xiàn)[3]中Gorsevski等使用了基于元胞自動(dòng)機(jī)演化的邊緣檢測(cè)方法對(duì)巖石圖像進(jìn)行分割。雖然已有較多的顆粒分割算法研究，但是能夠?qū)嶋H應(yīng)用的卻比較少。因?yàn)樵谧詣?dòng)分割之后，仍需要進(jìn)行大量的人工交互，對(duì)巖石顆粒的分割結(jié)果進(jìn)行修正，自動(dòng)化程度不高。

文獻(xiàn)[4]中提到，巖石樣品中晶體存在消光特性，在不同正交偏光角度下，顆粒的亮度會(huì)發(fā)生變化，顆粒之間的邊界也能清晰地顯示出來(lái)，所以可以利用這種特性來(lái)輔助巖石顆粒的分割。以此為依據(jù)，本文通過(guò)對(duì)采集到的不同正交偏光角度下的巖石薄片序列圖的分析處理，最終實(shí)現(xiàn)顆粒的自動(dòng)分割。

1 偏光序列圖采集處理

采集巖石薄片圖像一般是通過(guò)在顯微鏡下利用投射單偏光來(lái)獲取圖像，但是由于單偏光下巖石薄片圖像中的顆粒缺乏清晰的輪廓，顆粒之間沒(méi)有良好的區(qū)分度，如圖 1（a）所示。

本文使用正交偏光系統(tǒng)下采集的序列圖進(jìn)行圖像分割。正交偏光系統(tǒng)是在單偏光系統(tǒng)的基礎(chǔ)上再加上一個(gè)單偏光鏡，上、下偏光鏡振動(dòng)方向相互垂直，將巖石薄片放在上、下偏光鏡之間的載物臺(tái)上，視域中將呈現(xiàn)一系列光學(xué)現(xiàn)象，如消光、干涉色等[6]。采用自主研制的偏光圖像采集裝置，每旋轉(zhuǎn)正交偏光鏡15°獲取一幅圖像，總共旋轉(zhuǎn)120°，一共采集9張序列圖像。圖1（b）（c）（d）為其中 0°，60°，90°角度的 3 張圖像，可以看出，正交偏光圖像中的巖石顆粒清晰，但是由于顆粒的消光現(xiàn)象，在單張正交偏光圖中并不能把所有顆粒都顯示出來(lái)。

2 顆粒分割

圖1 單偏光圖像和不同角度下的正交偏光圖像

利用一組正交偏光序列圖像進(jìn)行顆粒分割，比較直觀的做法是對(duì)每一幅偏光序列圖下的圖像進(jìn)行分割，然后將多幅圖像的分割結(jié)果進(jìn)行合并，得到最后的分割結(jié)果。但是通過(guò)實(shí)驗(yàn)發(fā)現(xiàn)，該做法存在這樣的問(wèn)題：由于巖石圖像本身的復(fù)雜性和多樣性，分割中存在許多誤分割，對(duì)多幅圖像的處理圖進(jìn)行疊加處理，容易將多幅圖像的錯(cuò)誤分割結(jié)果相疊加。因此本文未選擇上述方法，而采用先將序列圖像疊加到一起，再做分割處理的方法。

2.1 偏光序列圖像的疊加

不同的顆粒在正交偏光鏡下的明暗變化的趨勢(shì)不同，這種趨勢(shì)如果相同，則可以合理的認(rèn)為是一個(gè)顆粒的概率較大[3]，本文將偏光鏡下采集的不同序列圖像之間的變化差值累加后得到疊加圖像，如式（1）、式（2）所示：

2.2 標(biāo)記分水嶺分割

本文對(duì)疊加后的圖像采用標(biāo)記分水嶺方法進(jìn)行顆粒分割，分水嶺分割[7]是一類(lèi)目前使用最為廣泛的顆粒形態(tài)學(xué)分割方法。其中H-minima[8]是一種常見(jiàn)的標(biāo)記選擇方法，采用梯度圖像中的區(qū)域極小值作為種子點(diǎn)[9]，實(shí)現(xiàn)目標(biāo)分割。本文根據(jù)文獻(xiàn)[10]選取方向測(cè)度的極小點(diǎn)作為標(biāo)記種子點(diǎn)，將每個(gè)像素點(diǎn)的鄰域擴(kuò)充到7×7，有效抑制了因?yàn)樵肼暬蛘哳w粒形狀不規(guī)則等因素而產(chǎn)生的多個(gè)虛假種子點(diǎn)，在一定程度上解決了基于梯度標(biāo)記分水嶺分割方法產(chǎn)生的嚴(yán)重過(guò)分割現(xiàn)象。

如圖3為采用分水嶺分割后用隨機(jī)顏色填充顆粒，并將背景色置為黑色后的圖，由此可以清晰看出，基本所有的顆粒都被提取出來(lái)。但在部分區(qū)域，如圖3中的標(biāo)記區(qū)域，仍存在過(guò)分割的現(xiàn)象。

圖2 正交偏光下采集的序列圖像的疊加圖像

圖3 分水嶺分割后的圖像

3 過(guò)分割區(qū)域聚類(lèi)

圖4 過(guò)分割的粘連區(qū)域

由圖3中標(biāo)記可見(jiàn)，對(duì)于相毗鄰的顆粒，可能會(huì)出現(xiàn)將其中的顆粒過(guò)分割成為若干個(gè)小區(qū)域的現(xiàn)象。圖4是將這部分區(qū)域放大后以一定透明度附加到原本的顆粒上的圖像。如圖4（a）所示，由于右上角的顆粒被過(guò)分割，導(dǎo)致原本只有兩個(gè)顆粒的該區(qū)域被分割成多個(gè)區(qū)域。針對(duì)這種現(xiàn)象，本文首先對(duì)過(guò)分割的圖像區(qū)域進(jìn)行角度，面積，位置特征提取，然后采用“自底向上”的AGNES層次聚類(lèi)算法進(jìn)行聚類(lèi)，將該區(qū)域聚合為符合實(shí)際情況的分割效果。

3.1 圖像區(qū)域特征表示

以圖5中的兩個(gè)顆粒為例，標(biāo)記上方顆粒為Grain-U,下方顆粒為Grain-D,其中Grain-U由于過(guò)分割導(dǎo)致被分成若干小的區(qū)域。

圖5 被過(guò)分割的兩個(gè)顆粒

本文從以下幾個(gè)方面來(lái)描述區(qū)域的特征：

（1）區(qū)域的角度特征

由于巖石顆粒正交偏光的周期消光性，正交偏光鏡下旋轉(zhuǎn)360°，顆粒會(huì)出現(xiàn)4次消光現(xiàn)象，每90°為一個(gè)周期[11]。如圖 6（a）（b）（c）（d）：分別是在正交偏光0°，45°，60°，105°下的顆粒圖像，可見(jiàn)上下兩個(gè)顆粒在不同偏光角度下呈現(xiàn)出來(lái)的亮度是不一樣的，每個(gè)區(qū)域在一個(gè)周期內(nèi)都存在亮度的最小值和最大值，且不同顆粒的亮度峰值所在角度不同。圖6（e）是Grian-U中的區(qū)域亮度變化折線圖，由于共同的消光性，雖然被過(guò)分割為多個(gè)不同的區(qū)域，但其亮度變化的趨勢(shì)一樣，均在75°達(dá)到峰值；而未被過(guò)分割的Grian-D顆粒在60°時(shí)亮度達(dá)到峰值。

因此可以用該峰值的角度來(lái)表征該區(qū)域的角度特征，如式（3）所示：

其中θLmax表示該區(qū)域亮度達(dá)到峰值時(shí)的角度，以π/2為底對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化操作，θ'為標(biāo)準(zhǔn)化之后的角度值。

（2）區(qū)域的面積特征

由于分水嶺的過(guò)分割現(xiàn)象，部分顆粒被分割成了許多小的區(qū)域。而這些小的區(qū)域的面積相對(duì)于其毗鄰顆粒的面積較小，具有相似性，因此面積可以作為聚類(lèi)的一個(gè)特征，如式（4）所示：

S'代表該區(qū)域的相對(duì)面積的大小，Ssum表示過(guò)分割區(qū)域的總面積，S為當(dāng)前區(qū)域的面積。

（3）區(qū)域的位置特征

由圖5（b）可以看出，Grian-U過(guò)分割區(qū)域的重心坐標(biāo)相對(duì)集中，而Grian-D的重心坐標(biāo)距離其他區(qū)域的相對(duì)較遠(yuǎn)。因此可以用區(qū)域的重心點(diǎn)坐標(biāo)來(lái)表征該區(qū)域所在的位置，如式（5）所示：

(X',Y')代表當(dāng)前區(qū)域重心坐標(biāo)的相對(duì)位置。其中W、H分別表示圖像的寬和高，(X,Y)為當(dāng)前區(qū)域的重心坐標(biāo)。

圖6 同一顆粒不同角度的明暗變化

根據(jù)公式（3）（4）（5）計(jì)算后的特征量值都被限定到[0,1]之間，從而組成一個(gè)四維的特征向量V={θ',S',X',Y'}.

3.2 過(guò)分割區(qū)域聚類(lèi)

本文采用“自底向上”AGNES層次聚類(lèi)算法，在不同層次對(duì)數(shù)據(jù)進(jìn)行劃分，從而形成樹(shù)形的聚類(lèi)結(jié)構(gòu)[12]。先將過(guò)分割區(qū)域內(nèi)所有的單個(gè)區(qū)域看成一個(gè)初始聚類(lèi)簇，然后找出距離最近的兩個(gè)簇進(jìn)行合并，該過(guò)程不斷重復(fù)，直至達(dá)到終止聚類(lèi)個(gè)數(shù)。本文計(jì)算簇間距離使用average linkage距離，簇間距離等于兩組對(duì)象之間的平均距離。給定聚類(lèi)簇Ci和Cj，其距離定義為：

其中d(Ci,Cj)表示不同簇Ci和Cj之間的距離，dist(Vp,Vq)是兩個(gè)區(qū)域p和q之間的歐氏加權(quán)距離，Vpi表示這個(gè)區(qū)域在i分量的特征值，ωi表示在i分量的權(quán)值。經(jīng)大量試驗(yàn)，區(qū)域的角度特征和面積特征在聚類(lèi)作用中比重較大，取亮度權(quán)重和面積權(quán)重ω1=ω2=5；坐標(biāo)位置權(quán)重ω3=ω4=3。

為了得到符合實(shí)際的聚類(lèi)個(gè)數(shù)，在使用層次聚類(lèi)的過(guò)程中，每一次兩個(gè)相近簇合并后，需要對(duì)該次合并的有效性進(jìn)行判斷，即衡量當(dāng)前聚類(lèi)的有效性。首先定義層次聚類(lèi)中的成本函數(shù)為：

其中，k是簇的個(gè)數(shù)，x表示簇中的數(shù)據(jù)，Ti是第i個(gè)簇中元素的集合；Ci表示第i個(gè)簇的中心，m為當(dāng)前簇中的元素個(gè)數(shù)。成本函數(shù)是各個(gè)簇畸變程度之和，若簇內(nèi)部的成員彼此間越緊湊則畸變程度越小，反之，若簇內(nèi)部的成員彼此間越分散則畸變程度越大。

根據(jù)文獻(xiàn)[13]，可以通過(guò)肘部法則來(lái)估計(jì)聚類(lèi)數(shù)量。肘部法會(huì)把不同k值的成本函數(shù)值反映在曲線上。隨著k值的增大，每個(gè)簇中包含的樣本數(shù)會(huì)減少，樣本離其中心會(huì)更近，因此平均畸變程度會(huì)減小。但是，隨著k值繼續(xù)增大，平均畸變程度的改善效果會(huì)不斷減低。當(dāng)Jk減少緩慢時(shí)，就認(rèn)為即使進(jìn)一步增大聚類(lèi)數(shù)，畸變程度的改善效果也并不能增強(qiáng)，則此時(shí)的這個(gè)“肘點(diǎn)”對(duì)應(yīng)的k值就是最佳聚類(lèi)數(shù)目，如圖7（a）。上述出現(xiàn)明顯“肘部”的情況比較常見(jiàn)，但是Jk值的變化也可能比較平緩沒(méi)有明顯的拐點(diǎn)，如圖7（b），這時(shí)可以認(rèn)為該區(qū)域最終聚為一類(lèi)。

圖7 J隨k值的變化曲線

本文設(shè)定當(dāng)聚類(lèi)數(shù)小于設(shè)定數(shù)kmax時(shí)，即當(dāng)k取k=1，2，…，kmax，計(jì)算其成本函數(shù)。根據(jù)計(jì)算的結(jié)果得到Jk值變化對(duì)應(yīng)的k值，則此k值就為最終的聚類(lèi)個(gè)數(shù)，如果隨著聚類(lèi)數(shù)目的減小，Jk的變化幅度始終相當(dāng)，變化曲線上沒(méi)有明顯的拐點(diǎn),則把這個(gè)區(qū)域歸為一類(lèi)。

聚類(lèi)過(guò)程為：首先設(shè)定一個(gè)過(guò)分割區(qū)域內(nèi)的所有區(qū)域各為一簇，根據(jù)公式（6）計(jì)算所有簇間距離，將距離最近的兩個(gè)簇合并，然后更新整個(gè)簇的個(gè)數(shù)；再計(jì)算新的簇之間的距離，根據(jù)距離最近原則再合并兩個(gè)簇，如此反復(fù)。當(dāng)簇個(gè)數(shù)小于等于設(shè)定值kmax后，開(kāi)始計(jì)算其成本函數(shù) Jk，并保存每次更新的簇的結(jié)果，直到最終聚為一類(lèi)。

每合并一簇后計(jì)算每次變化的幅度：

計(jì)算ηk的平均值kavg，比較每次聚類(lèi)后當(dāng)前幅度值與平均值的大小，若滿足ηk≥δ×kavg，則判斷其為“肘部”，對(duì)應(yīng)的k值即為最佳聚類(lèi)個(gè)數(shù)，k值對(duì)應(yīng)的聚類(lèi)結(jié)果即為最終結(jié)果。如果均不滿足，即不存在明顯“肘部”，此時(shí)k取1，聚為一類(lèi)。經(jīng)過(guò)試驗(yàn)，當(dāng)設(shè)定δ=2時(shí)，聚類(lèi)結(jié)果較符合實(shí)際情況。

將同一類(lèi)區(qū)域填充為同一種顏色，根據(jù)上述算法，圖5的最終聚類(lèi)圖如圖8所示：

圖8 Grain-U和Grain-D最終聚類(lèi)圖

再分別對(duì)圖4中的過(guò)分割區(qū)域進(jìn)行聚類(lèi)，結(jié)果如圖9所示，可以看出，本文算法成功將過(guò)分區(qū)域聚類(lèi)為符合實(shí)際情況的結(jié)果。

4 實(shí)驗(yàn)結(jié)果與分析

本文對(duì)二組實(shí)際拍攝的偏光序列圖進(jìn)行顆粒分割，并與文獻(xiàn)[6]中的顆粒分割方法算法進(jìn)行對(duì)比。

顆粒分割結(jié)果如圖10所示，其中：（a）（d）分別為兩組偏光序列圖進(jìn)行差分疊加后的圖像；（b）（e）分別為通過(guò)文獻(xiàn)[6]中SRM算法分割后的效果圖；（c）（f）分別為本文方法對(duì)顆粒的分割效果圖。文獻(xiàn)[6]中先對(duì)每一張圖片進(jìn)行分割后再合并，提取到的顆粒比較多，但是過(guò)分割現(xiàn)象嚴(yán)重。本文中先將序列圖合并后再進(jìn)行分割，并且對(duì)過(guò)分割區(qū)域進(jìn)行聚類(lèi)合并，過(guò)分割現(xiàn)象得到很好的抑制，使得結(jié)果較于文獻(xiàn)[6]的方法更符合實(shí)際需求。

5 結(jié)語(yǔ)

本文討論了一種正交偏光圖像序列圖的顆粒分割算法，首先采集不同偏光角度下的序列圖，對(duì)序列圖進(jìn)行疊加融合后再基于分水嶺算法分割將目標(biāo)顆粒提??；然后對(duì)過(guò)分割區(qū)域提取角度、面積、位置等特征；最后通過(guò)AGNES聚類(lèi)算法將過(guò)分割的顆粒聚類(lèi)到一起，從而完成顆粒的分割提取。

圖9 過(guò)分區(qū)域聚類(lèi)后效果圖

圖10 文獻(xiàn)[3]算法以及本文算法對(duì)比效果圖

實(shí)驗(yàn)結(jié)果顯示，本文算法通過(guò)先分離后聚合的方式，成功將過(guò)分割的區(qū)域融合。但對(duì)于沒(méi)有在分水嶺階段中沒(méi)有提取到的顆粒，本文算法不能將其分割出來(lái)，可在這方面繼續(xù)加以研究。

[1]徐耀鑒.巖石學(xué)[M].北京:地質(zhì)出版社,2007.

[2]Obara B.A New Algorithm Using Image Colour System Transformation for Rock Grain Segmentation[J].Mineralogy and Petrology,2007,91（3-4）:271-285.

[3]Gorsevski P V,Onasch C M,Farver J R,et al.Detecting Grain Boundaries in Deformed Rocks Using a Cellular Automata Approach[J].Computers&Geosciences,2012,42（3）:136-142.

[4]楊宗瑞.基于偏光序列圖像的巖石顆粒分割及礦物分[D].四川大學(xué),2014

[5]Gorsevski P V,Onasch C M,Farver J R,et al.Detecting Grain Boundaries in Deformed Rocks Using a Cellular Automata Approach[J].Computers&Geosciences,2012,42（3）:136-142.

[6]吳擁,蘇桂芬,滕奇志,等.巖石薄片正交偏光圖像的顆粒分割方法[J].科學(xué)技術(shù)與工程,2013,13（31）:9201-9206.

[7]Beucher S,Meyer F.Segmentation:The Watershed Transformation.Mathematical Morphology in Image Processing[J].Optical Engineering,1993,34（5）:433-481.

[8]Chanho J,Changick K.Segmenting Clustered Nuclei Using H-minima Transform-based Marker Extraction and Contour Parameterization[J].IEEE Transactions on Bio-medical Engineering,2010,57（10）:2600-4.

[9]Li B,Pan M,Wu Z.An Improved Segmentation of High Spatial Resolution Remote Sensing Image Using Marker-based Watershed Algorithm[C].International Conference on Geoinformatics.IEEE,2012:1-5.

[10]楊海軍,梁德群,畢勝.基于圖象方向性信息測(cè)度的圖象象素分類(lèi)[J].中國(guó)圖象圖形學(xué)報(bào),2001,6（5）:429-433.

[11]Fueten F,Goodchild J S.Quartz C-axes Orientation Determination Using The Rotating Polarizer Microscope[J].Journal of Structural Geology,2001,23（6-7）:895-902.

[12]周志華.機(jī)器學(xué)習(xí)[M].清華大學(xué)出版社,2016.

[13]Bholowalia P,Kumar A.EBK-Means:A Clustering Technique based on Elbow Method and K-means in WSN[J].International Journal of Computer Applications,2014.