基于自適應(yīng)H-極大值的粘連顆粒分割算法

何磊，王正勇，滕奇志，田剛

（1.四川大學(xué)電子信息學(xué)院，成都 610065；2.新疆油田公司行政事務(wù)中心，克拉瑪依 834000）

0 引言

隨著CT巖心掃描成像技術(shù)日益成熟，利用CT巖心序列圖重建出巖心三維結(jié)構(gòu)成為一種重要的數(shù)字巖心分析手段。由于巖心中巖石顆粒相互擠壓和粘連，在對巖心三維圖像進(jìn)行顆粒的粒度分析時(shí)，存在的粘連顆粒會(huì)影響顆粒的參數(shù)計(jì)算準(zhǔn)確性。目前針對粘連顆粒分割的問題，研究者們根據(jù)巖心三維體顆粒（如圖1所示）具有類球性差、輪廓信息不規(guī)則以及數(shù)據(jù)量大等特征，提出了不同的思路：

圖1 CT巖心序列圖

（1）基于腐蝕膨脹的粘連分割方法[1]。此方法對粘連顆粒進(jìn)行連續(xù)腐蝕操作，直至粘連顆粒分開，再進(jìn)行膨脹操作恢復(fù)顆粒。特點(diǎn)是原理簡單，運(yùn)算速度快，對粘連程度較小的顆粒能實(shí)現(xiàn)快速準(zhǔn)確分割，但是針對粘連程度較大的顆粒，存在腐蝕次數(shù)不易控制和目標(biāo)形變嚴(yán)重等問題。

（2）基于骨架的分割方法[2]。此方法通過統(tǒng)計(jì)垂直于顆粒骨架的截面的變化規(guī)律來進(jìn)行分割。但是骨架提取算法和截面變化規(guī)律統(tǒng)計(jì)的計(jì)算量大，并且容易受到噪聲的影響而出現(xiàn)一些冗余的骨架段，不適合形狀復(fù)雜的巖石顆粒分割。

（3）基于分水嶺的分割方法[3]。此方法利用二值圖的距離變換圖中的局部極大值作為種子點(diǎn)，實(shí)現(xiàn)分水嶺分割。但是因?yàn)閹r石顆粒的形狀復(fù)雜以及噪聲點(diǎn)，距離變換矩陣中會(huì)出現(xiàn)冗余的局部極大值點(diǎn)，會(huì)導(dǎo)致嚴(yán)重的過分割現(xiàn)象。本文借鑒分水嶺算法的思想，提出了一種三維粘連顆粒的分割算法：利用自適應(yīng)的H-極大值變換抑制冗余的局部極大值[4]，然后以巖石顆粒的最優(yōu)形狀因子作為目標(biāo)函數(shù)，利用基于三維顆粒粘連程度的合并算法進(jìn)行區(qū)域合并，從而得到最終的分割結(jié)果[5]。

1 巖石顆粒形狀度量指標(biāo)

顆粒的球度公式定義由Wadell[6]提出，其文中通過計(jì)算顆粒等效球的表面積和顆粒表面積（顆粒等效球?yàn)榕c顆粒體積相等的球），得到兩者的比值即為球度，球度公式如下：

其中Se為顆粒等效球的表面積，Sp為顆粒的表面積。公式（1）是作為經(jīng)典公式被廣泛使用，文獻(xiàn)[7-9]均采用此球度公式作為顆粒的形狀度量，其中文獻(xiàn)[9]將公式（1）變換成二維的圓度公式，分別計(jì)算三維顆粒在x、y和z方向投影對應(yīng)的三個(gè)圓度值，取三個(gè)圓度值的加權(quán)平均值作為顆粒的形狀度量。

另外，文獻(xiàn)[10]提出形狀因子計(jì)算公式：

其中，V表示顆粒體積，Sp表示顆粒表面積。Φ的取值范圍是（0,1]，顆粒形狀越接近于球，Φ的取值越接近于1，若顆粒形狀呈現(xiàn)面狀，則Φ的取值趨向于0。公式（1）形狀因子表達(dá)式與公式（2）球度表達(dá)式具有相近的數(shù)學(xué)內(nèi)涵，把顆粒體積V代入式（1）中得到式（3）：

在程序的計(jì)算過程中，式（1）需要計(jì)算顆粒等效球的表面積Se，Se的計(jì)算首先需要計(jì)算顆粒體積V，然后由球的體積與表面積之間的關(guān)系式推導(dǎo)出Se，而式（2）沒有推導(dǎo)Se的步驟,得到顆粒體積V和顆粒表面積Sp即可計(jì)算形狀因子。本文對于巖石顆粒的形狀度量采用公式（3）。

2 基于自適應(yīng)H-極大值的粘連顆粒分割算法

巖石顆粒具有邊緣復(fù)雜和形狀不規(guī)則的特征，傳統(tǒng)分水嶺變換能將粘連的顆粒分割開，但是容易產(chǎn)生如圖2所示的過分割現(xiàn)象。因?yàn)轭w粒的距離變換圖中普遍存在不止一個(gè)局部最大值的情況，直接選取每個(gè)局部最大值點(diǎn)作分水嶺變換的種子點(diǎn)會(huì)產(chǎn)生嚴(yán)重的過分割現(xiàn)象。本文通過使用H-極大值變換來抑制冗余的局部極大值點(diǎn)，同時(shí)加入第1節(jié)定義的巖石顆粒形狀度量指標(biāo)實(shí)現(xiàn)了自適應(yīng)的H-極大值變換,為了防止h值的增大導(dǎo)致正確的局部極大值點(diǎn)被過濾掉，H-極大值變換后的子分割區(qū)域集合還加入了區(qū)域合并算法。另外，由于經(jīng)過三維重建的巖石顆粒的表面以及鄰域空間會(huì)產(chǎn)生許多細(xì)紋以及噪聲點(diǎn)，在進(jìn)行H-極大值變換之前利用開重建運(yùn)算去除因三維重建過程中產(chǎn)生的噪聲點(diǎn)[11]。

圖2 分水嶺變換的過分割現(xiàn)象

2.1 自適應(yīng)H-極大值變換

分水嶺變換[12]的過程可以理解為降雨的過程，最后得到集水盆的集合，不同的集水盆對應(yīng)不同的目標(biāo)區(qū)域。因?yàn)橐粋€(gè)集水盆對應(yīng)于距離變換矩陣的一個(gè)局部極大值區(qū)域。所以過分割意味著有冗余的局部極大值點(diǎn)，要解決過分割就是要去除冗余的局部極大值點(diǎn)。

圖3 顆粒的距離變換圖

如圖3所示，黑色像素點(diǎn)為背景，不同顏色的像素點(diǎn)代表不同的距離值，距離值由目標(biāo)邊緣向目標(biāo)中心逐漸增大，局部極大值存在于目標(biāo)中心區(qū)域。上圖中我們期望的顆粒數(shù)是1，即只出現(xiàn)一個(gè)局部極大值，但圖中存在3個(gè)局部極大值區(qū)域1、4、5（對應(yīng)的距離值分別為17、19和20）。如果直接用分水嶺變換，會(huì)得到三個(gè)目標(biāo)，這顯然是錯(cuò)的。因?yàn)榫嚯x變換受目標(biāo)的局部噪聲的影響出現(xiàn)了多個(gè)局部極大值，如果有辦法過濾掉多余的局部極大值點(diǎn)，使得過濾后的局部極大值點(diǎn)的數(shù)目與我們期望的顆粒數(shù)目一致，就能得到正確的結(jié)果，達(dá)到抑制過分割的目的。

H-極大值變換[13]的作用是過濾掉距離值小于閾值h的冗余局部極大值點(diǎn)。h=0為起始點(diǎn)，通過以Δh=0.5逐步增大h值，逐步減少子分割區(qū)域的數(shù)目，找到能使當(dāng)前子分割集合的加權(quán)平均形狀因子最大的h值，從而實(shí)現(xiàn)自適應(yīng)的H-極大值變換。加權(quán)平均形狀因子公式和變換公式如下：

上述公式中Vtotal表示子分割區(qū)域的總體積，vi表示子分割顆粒的體積表示對應(yīng)h值變換下的加權(quán)平均形狀因子，maxpts表示子分割區(qū)域最大的局部極大值，hmax則表示子分割區(qū)域最終的H-極大值變換使用的h值。其中表示對標(biāo)記圖像I的測地膨脹，標(biāo)記圖像I是原圖像的距離變換矩陣中局部極大值點(diǎn)的集合。

圖4展示了不同h值對應(yīng)不同巖石顆粒子分割集合的加權(quán)平均形狀因子。圖4（b）是原始圖像（a）的分水嶺算法分割結(jié)果，可見一顆完整的顆粒被過分割成三顆顆粒。逐步增大h值，直到h值增加到5也不能抑制過分割的情況，說明冗余的局部極大值大于5。當(dāng)h值為5.5時(shí)，如圖4（c）所示，一個(gè)冗余局部極大值被過濾掉。當(dāng)h值為8時(shí)，圖4（a）中兩個(gè)冗余的局部極大值點(diǎn)均被過濾掉，得到子分割結(jié)果（d），此時(shí)加權(quán)平均形狀因子最大，所以h=8對應(yīng)的H-極大值變換得到的子分割集合作為候選結(jié)果。以此類推，對每個(gè)連通區(qū)域都可自適應(yīng)地選取最優(yōu)h值。

2.2 合并規(guī)則

合并規(guī)則是作為h值的繼續(xù)迭代的約束條件，防止h值的增大把正確的局部極大值點(diǎn)也過濾掉，從而導(dǎo)致欠分割現(xiàn)象。由于巖心三維顆粒的類球性較差，粘連復(fù)雜，為此本文提出了基于粘連程度合并規(guī)則，這樣合并規(guī)則的適應(yīng)性會(huì)更強(qiáng)。粘連程度定義如下：

如圖5所示，展示了兩個(gè)粘連顆粒分水嶺變換的距離變換圖，b和c表示兩個(gè)粘連顆粒的局部最大值，即為分水嶺變換選取的種子點(diǎn)，a表示粘連處的局部最大值。當(dāng)兩個(gè)顆粒的粘連程度大于設(shè)定的閾值，就合并這兩個(gè)顆粒。如果粘連閾值過大，會(huì)導(dǎo)致部分過分割的顆粒無法合并。如果閾值過小，會(huì)讓本來不該合并的兩個(gè)顆粒合并。本文通過多次實(shí)驗(yàn)對比分割效果，選取分割效果最好的粘連閾值，本文實(shí)驗(yàn)結(jié)果對應(yīng)的粘連閾值為0.8。為了便于觀察且說明問題，采用程序隨機(jī)生成的球型數(shù)據(jù)為實(shí)驗(yàn)對象，如圖6所示，展示了不同粘連閾值對應(yīng)的不同的合并效果。

圖4 巖石顆粒的h值以及形狀因子

圖5 粘連程度的定義

圖6

2.3 算法步驟

定義圖像經(jīng)過開重建濾波處理后的N個(gè)連通分量集合I={Ij|j∈{1,2,…,N}}。對每一個(gè)Ij先后進(jìn)行距離變換和分水嶺分割，得到初始分割結(jié)果。定義Ij={Sn(Ij)|n∈{1,2,…,m}}為每個(gè)連通分量的初始分割結(jié)果，分割后的子區(qū)域?yàn)閙個(gè)，需要對m>1的連通分量進(jìn)一步處理，m=1的連通分量為最終分割結(jié)果。進(jìn)一步處理的算法如下：

（1）初始化連通分量的序號j=1；

（2）取連通分量Ij，及其初始分割結(jié)果Sn(Ij)，n∈{1,2,…,m},子分割數(shù)目是m，初始化閾值h=0，最優(yōu)平均形狀因子opt_shape=0；

（3）當(dāng)mj(h)>1且mj()h

（4）若opt_shape

（5）令h=h+Δh，進(jìn)行閾值為h的H-極大值變換，跳轉(zhuǎn)至 c）；

（7）j=j+1，若j

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

本文算法應(yīng)用在由巖心CT序列圖重建得到的三維模型。如圖7所示，由兩組CT序列圖提取出巖石顆粒分別得到圖7（a）和（b）所示的巖石顆粒三維模型，然后對三維模型中的巖石顆粒進(jìn)行分割，圖7（a）和（b）對應(yīng)的分割結(jié)果分別為圖8和圖9。圖8（a）為傳統(tǒng)分水嶺算法分割后的三維模型，圖8（b）為本文算法分割后的結(jié)果，兩圖均為偽彩色圖像，不同顏色的三維體目標(biāo)代表分割后的巖石顆粒。圖8中標(biāo)識(shí)的1、2、3和4顯示了兩種算法效果的差別，由圖可見本文算法有效抑制了傳統(tǒng)分水嶺的過分割現(xiàn)象。對于形狀復(fù)雜且不規(guī)則的顆粒圖像，在最優(yōu)加權(quán)平均形狀因子和合并規(guī)則約束下，H-極大值變換能保證每一個(gè)分割后的顆粒都接近于真實(shí)情況。

圖7

圖8

圖9

4 結(jié)語

針對三維模型中存在巖石顆粒粘連的問題，本文提出了一種基于自適應(yīng)H-極大值的粘連顆粒分割算法，對實(shí)際工程中遇到的三維粘連顆粒進(jìn)行了分割。實(shí)驗(yàn)結(jié)果表明，本文提出的算法具有分割準(zhǔn)確的優(yōu)點(diǎn)，能夠解決實(shí)際工程中三維粘連顆粒的分割問題。