基于自適應H-極大值的粘連顆粒分割算法

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

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

0 引言

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

圖1 CT巖心序列圖

（1）基于腐蝕膨脹的粘連分割方法[1]。此方法對粘連顆粒進行連續腐蝕操作，直至粘連顆粒分開，再進行膨脹操作恢復顆粒。特點是原理簡單，運算速度快，對粘連程度較小的顆粒能實現快速準確分割，但是針對粘連程度較大的顆粒，存在腐蝕次數不易控制和目標形變嚴重等問題。

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

（3）基于分水嶺的分割方法[3]。此方法利用二值圖的距離變換圖中的局部極大值作為種子點，實現分水嶺分割。但是因為巖石顆粒的形狀復雜以及噪聲點，距離變換矩陣中會出現冗余的局部極大值點，會導致嚴重的過分割現象。本文借鑒分水嶺算法的思想，提出了一種三維粘連顆粒的分割算法：利用自適應的H-極大值變換抑制冗余的局部極大值[4]，然后以巖石顆粒的最優形狀因子作為目標函數，利用基于三維顆粒粘連程度的合并算法進行區域合并，從而得到最終的分割結果[5]。

1 巖石顆粒形狀度量指標

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

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

另外，文獻[10]提出形狀因子計算公式：

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

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

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

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

圖2 分水嶺變換的過分割現象

2.1 自適應H-極大值變換

分水嶺變換[12]的過程可以理解為降雨的過程，最后得到集水盆的集合，不同的集水盆對應不同的目標區域。因為一個集水盆對應于距離變換矩陣的一個局部極大值區域。所以過分割意味著有冗余的局部極大值點，要解決過分割就是要去除冗余的局部極大值點。

圖3 顆粒的距離變換圖

如圖3所示，黑色像素點為背景，不同顏色的像素點代表不同的距離值，距離值由目標邊緣向目標中心逐漸增大，局部極大值存在于目標中心區域。上圖中我們期望的顆粒數是1，即只出現一個局部極大值，但圖中存在3個局部極大值區域1、4、5（對應的距離值分別為17、19和20）。如果直接用分水嶺變換，會得到三個目標，這顯然是錯的。因為距離變換受目標的局部噪聲的影響出現了多個局部極大值，如果有辦法過濾掉多余的局部極大值點，使得過濾后的局部極大值點的數目與我們期望的顆粒數目一致，就能得到正確的結果，達到抑制過分割的目的。

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

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

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

2.2 合并規則

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

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

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

圖5 粘連程度的定義

圖6

2.3 算法步驟

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

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

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

（3）當mj(h)>1且mj()h

（4）若opt_shape

（5）令h=h+Δh，進行閾值為h的H-極大值變換，跳轉至 c）；

（7）j=j+1，若j

3 實驗結果對比與分析

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

圖7

圖8

圖9

4 結語

針對三維模型中存在巖石顆粒粘連的問題，本文提出了一種基于自適應H-極大值的粘連顆粒分割算法，對實際工程中遇到的三維粘連顆粒進行了分割。實驗結果表明，本文提出的算法具有分割準確的優點，能夠解決實際工程中三維粘連顆粒的分割問題。