巖石鑄體薄片圖像的喉道分割算法

涂秉宇，熊淑華，滕奇志，田剛

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

0 引言

在儲集層巖石微觀結構中，由三個或三個以上巖石顆粒所包圍的空間被稱為孔腔，相鄰兩個孔腔之間的最窄連接部分稱為喉道，孔腔和連接它的部分喉道的總體稱為孔隙[1]。喉道是相鄰巖石孔隙之間最窄的部分，其大小和分布都是預測多孔介質的滲流性的重要因素[2]。因此對于喉道的研究是巖石孔隙結構分析的重要一環。

目前，微觀孔隙結構的分析過程是：利用鑄體薄片在光學顯微鏡下選取典型視域成像，再在獲取的圖像上根據鑄體薄片中孔隙的顏色特性將孔隙部分作為目標提取出來，然后對孔隙部分進行喉道的劃分。而喉道一般是處在目標區域凹陷或凸出最窄處，因此，對于喉道的分割可以通過對目標圖像拐點處的檢測、篩選與匹配來實現。最后再根據石油天然氣行業標準[1]進行參數計算，從而得出微觀孔隙結構的數據。

本文提出了一種基于拐點檢測和匹配的自動喉道分割方法，該方法通過檢測拐點信息，得到拐點的位置，再通過兩兩匹配拐點，形成拐點對，以此來確定喉道位置，從而實現對巖石孔隙連通區域的分割。

1 喉道分割算法

對于喉道分割給出以下示例，圖1是一幅孔隙原圖及其分割結果圖，（a）是孔隙原圖，（b）是孔隙提取結果圖，（c）是將孔隙部分疊加到原圖的效果圖，（d）是喉道分割結果圖，其中紅色分割線為喉道。由此可以看出，喉道分割的本質在圖像處理中的意義是找出孔腔與孔腔之間最窄的部分。由于喉道往往處在邊界的凹陷處，如果對目標邊界進行拐點檢測，再根據拐點的位置便能準確找出邊界凹陷處，基于此就能準確地進行喉道分割。所以本文提出了基于拐點的喉道分割算法，該算法首先通過八鄰域法獲得孔隙二值圖的邊界，引入邊界鏈碼代替曲率，通過鏈碼的變化找出孔隙邊界上的拐點。對找出的拐點進行篩選，過濾掉偽拐點，得出有效拐點。然后根據巖石孔隙結構的物理特性對有效拐點兩兩匹配，并處理拐點數為奇數的情況，從而形成喉道，最終得出喉道分割位置。

圖1 孔隙示例圖

1.1 獲得目標區域的邊界

孔隙提取后的效果圖一般為二值圖，而對于二值圖目標區域的邊界判定可以使用八鄰域連通標記法。所以本文對孔隙區域采用八鄰域法進行標記判別，具體步驟為：首先針對孔隙圖片的提取結果進行處理，其順序是從左上方的第一個像素點開始，從左到右，從上到下依次遍歷，并按照圖2所示掩膜的方法對孔隙區域進行標記，直至整幅圖搜索結束。最后根據八鄰域判別邊界，對標記的孔隙區域的所有像素點進行迭代，判斷其鄰域是否屬于該標記區域，若其不屬于，則斷定為該區域的邊界，否則判斷為該區域所屬的點。

圖2 八鄰域掩膜

1.2 邊界曲線的拐點檢測

（1）曲線的點曲率

邊界曲線進行拐點檢測需對曲線上的點進行曲率計算，然后根據曲率判斷拐點位置。

如圖3所示，假設函數y=f（x）代表的是一條曲線，P1和P2是曲線上任意兩點，曲線上P1點的曲率k定義為點沿著曲線的切線方向與水平X軸的夾角的變化[4-5]。

α1為點P1的切線夾角，α2為點P2的切線夾角，△l為曲線P1P2的弧長，△α為切線變化的角度。

圖3 曲率示意圖

（2）鏈碼

由曲率的定義可知，其計算需要求極限和乘除運算，相對較復雜[6]。本文引入計算更為簡便的鏈碼來表示邊界曲率，避免乘除法運算，可以有效地節省計算時間。

所謂鏈碼是用來反映像素點和其鄰近像素點方向的一種編碼[7-8]，其編碼方式如圖4，用當前像素點指向它的八鄰域方向來表示。

圖4 八連通域鏈碼

鏈碼可以表示邊界的斜率、曲率，即當前邊界點的切線方向。鏈碼差則是和曲率成正比的量[3]，其計算公式為：

式中D（xi）為邊緣曲線上第i個像素點的鏈碼值。但是由于只有0～7八個方向，細化程度差，導致精度缺失。文獻[9]算法根據鏈碼特性利用其計算得到的曲率進行局部平均從而得以提升算法效果。而文獻[10]引入了N點鏈碼差，利用進入和離開點Xi的N個鏈碼平均值之差，提升算法精確度。以上算法都只是針對八鄰域進行方向編碼，而八鄰域方向編碼將2π平面分為8個方向，每個方向精度的偏差達到π/4，這存在量化精度不足的問題。因此本文采用了K鄰域鏈碼中的16鄰域鏈碼方法[6]。16鄰域鏈碼在8鄰域鏈碼的基礎上向外擴展一圈，將2π平面分為16個方向，這使得16鄰域鏈碼相較于8鄰域鏈碼可以在不增加過多運算量的情況下提升方向精度和增加邊界平滑[10]，這對于邊界走勢復雜的孔隙區域而言具有實用性和可行性。

（3）拐點檢測算法

本文基于鏈碼的拐點檢測算法如下：

①對提取得到的邊緣輪廓曲線進行16鄰域鏈碼編碼。

②按照曲率定義，使用差分替代微分，則邊緣曲線第i點的切線角度變化的差分表示值θi，可通過公式（3）求得：

其中，C（16,i）表示第i個像素的16鄰域鏈碼值。

由于θi是拐角的差分表示，所以應該將其角度歸一化到[0,π]。而在16鄰域鏈碼編碼中，[0,π]的區間對應的鏈碼區間為[0,8]，所以點Pi處的曲率計算為：

為了進一步擴大拐點和非拐點處曲率的差距，所以使用以下公式，得到最終曲率值：

③根據曲率值得到曲率局部峰值的位置。即確定一個閾值 t，當 ei≥t>ei-1時曲率曲線為上升，當 ei>t≥ei+1時，曲線為下降，那么在此范圍內與波峰位置左右相對應的邊緣點即是待選拐點。其中t值應為邊界曲率值的平均值，計算公式如下：

④判斷凹凸性。以相距拐點前后相鄰n點的中點為判別條件，如果中點像素值屬于目標區域則為凸點，否則為凹點。

孔隙連通區域內部存在內部孔洞時，內部孔洞目標的邊界為內輪廓，孔隙連通區域邊界為外輪廓。對內輪廓取其凸點，經由以下篩選和匹配等步驟后與外輪廓上的凹點形成喉道分割線，用以進行有效的喉道分割。

1.3 篩選算法

由于孔隙區域邊界凹陷處呈塊狀，而凹陷塊中曲率的變化會導致檢測出多個拐點，同時檢測出的拐點中也可能有不屬于喉道位置的噪聲點，所以經由拐點檢測后得出的拐點僅是待選拐點。還應對待選拐點進行篩選，去除相鄰點和噪聲點的干擾。為了達到消除噪聲和冗余點的目的，本文采用的篩選算法如下：

（1）對于孔隙連通區域目標邊緣凹陷塊的情況，根據孔隙形狀特征，本文通過某個距離閾值來判斷待選凹點是否相鄰，在閾值范圍內的待選凹點屬于相鄰凹點。相鄰的待選凹點，將其看作一個凹點群，取其中曲率最大的點作為該組凹點的有效凹點，代表這組凹點。而對于孔隙連通區域內部的孔洞，邊緣檢測形成的內輪廓上會有凸起塊，通過拐點檢測后找到的相鄰凸點即為一個凸點群，取其中曲率最大的點作為有效凸點。因此，檢測到的前后兩個拐點之間的距離小于某個閾值的時候，我們認為這兩個拐點是相鄰的，應舍棄曲率較小的點。通過實驗發現，兩點間距離閾值一般設定為4～6個像素。

（2）由于孔隙邊界情況的復雜性，所以由第一步篩選之后所得的拐點，容易出現拐點比較“淺”，即如圖5所示的情況：通過第一步篩選算法后，得出A為有效拐點，但是A到直線BC的距離卻很短，即A到直線的距離值非常的小。根據喉道的定義可知，喉道為兩孔腔之間最短的地方，所以如圖5所示的點不應該計入有效拐點。

圖5 由于太“淺”不能作為拐點

因此本文提出確定一個有效的d值對此類拐點進行判別。由于全薄片圖像中孔隙區域的邊緣情況復雜，如果d值的選取過小，會導致一些原本不是匹配點的拐點被誤判；取值過大，則會導致漏掉一些真正的拐點。所以d的取值不能為固定值。本文算法采用了自適應d值，即使用所有孔隙目標根據如下公式計算得到平均面積大小：

最后根據該面積大小的等效圓半徑確定適中的d值。等效圓半徑r計算如式（8）：

（3）對于所有已得的拐點進行以上兩步篩選迭代，直至所有拐點篩選完畢。

1.4 匹配算法

喉道處在孔隙邊界兩端凹陷處，而要找到相應喉道則必須對邊界兩端凹陷處的拐點進行配對，如圖1（d）所示。因此通過篩選算法得出的有效拐點還需進行拐點對匹配算法，從而形成喉道線對孔隙連通區域進行有效的分割。

本文采用的匹配算法如下：

（1）喉道線的有效凹點對應該是位于孔隙連通區域兩側，如圖6所示。

兩個有效凹點待匹配的示意圖，其中的E1和E2為有效凹點，D1、F1和 D2、F2分別是 E1點和 E2點的等距前繼點與后繼點。將兩個待匹配凹點的連線向兩端進行延伸，如果E1一側的延伸線位于E1D1向量與E1F1向量形成的銳角區域內，E2一側的延伸線位于E2D2向量與E2F2向量形成的銳角區域內，且E2位于E1D1和E1F1的反向延長線所覆蓋的邊界區域內。只有上述條件均滿足的情況下，滿足匹配條件，該組有效凹點是相對應匹配的，否則該組有效凹點不對應匹配。

（2）通過判別步驟1進行匹配后，如果出現當前待匹配的有效凹點對應多個待匹配點的情況。根據孔隙喉道的定義，喉道是孔腔與孔腔之間最短的部分，所以需引入距離判別進一步判別正確的凹點匹配對，即只要分割線穿越目標區域的歐氏距離短即可滿足此條件，如圖7所示。

圖7 最短距離匹配

A1點對應候選的匹配點存在A2、B2兩個點，計算A1A2和A1B2的歐氏距離，取距離最短的A2進行匹配，歐氏距離公式如下：

（3）由于喉道是最窄處，而根據圖1中喉道的寬度和其所屬孔隙連通區域的等效圓直徑比較可知，分割線線長必定滿足不大于某個D值。即D應滿足公式（10）：

其中C為孔隙連通區域等效圓周長，D為兩匹配凹點間距離。

（4）由于分割線必須位于孔隙內部，所以喉道分割線上的像素點應該全部屬于孔隙連通區域。因此對匹配點對連線上的點進行抽樣分析，如果存在像素點屬于背景區域則舍棄該喉道，否則該連線可以作為喉道。

（5）對于孔隙連通區域內部存在孔洞目標的情況，應該將內輪廓上檢測到的凸點進行篩選，再與外輪廓上的凹點根據以上四點匹配準則進行有效匹配，形成分割線，如圖8所示。

圖8 最短距離匹配

1.5 奇數拐點的處理方法

上文的匹配算法是在一個孔隙連通區域中恰好篩選出了偶數個有效拐點的情況下進行一一匹配的。但是當遇到奇數個拐點的情況時，會在孔隙局部區域留下未做匹配的有效拐點。因此，根據喉道定義，本文通過選取對角邊緣延長線所覆蓋的邊界區域進行最短距離判斷，找到距離有效拐點最近的邊界點與此拐點進行匹配，如圖9所示。

圖9 匹配非拐點

反向延長線范圍可以找出邊界上距離A、C點最近的B、D點，同時通過上文的匹配算法可以有效的判斷AB、CD是否能夠成為分割線。

通過此法，我們可以處理拐點數為奇數的情況，也防止了欠分割的情況，使得基于拐點的喉道分割算法更具有實際應用的價值。

2 測試與分析

為了驗證本文算法的可行性，我們采用實際生產中的鑄體孔隙圖片，對不同情況下的孔隙進行實驗。圖10分別給出了不同情況下的孔隙原圖及其對應的分割效果圖，左邊是原圖，右邊是對應的分割效果圖。

圖10 實驗結果圖

其中（a）（b）為孔隙中有一個孔洞情況的示例。圖（c）（d）為孔隙中有多個孔洞情況的示例。圖（e）（f）為孔隙中無孔洞但是區域形狀復雜的情況示例。從圖11中可看出，本文算法能適應各種情況的喉道分割。

同時本文針對基于形態學流域分割目標方法和腐蝕膨脹分割目標方法的結果進行實驗比對。圖11給出了不同算法的分割結果對比圖，其中（a）為原圖。（b）為基于流域算法的分割結果，其能有效地分割出一部分喉道，但是由于其形成分水嶺的過程中受到孔隙部分局部極值影響，從而導致過分割現象。（c）（d）（e）分別為基于腐蝕膨脹算法17、20、25次后的分割結果圖，圖（c）在17次的腐蝕膨脹后會出現欠分割現象。圖（d）進行了20次腐蝕后會出現小顆粒被完全腐蝕掉，存在欠分割的現象。圖（e）腐蝕膨脹25次導致粘連性較大的區域得不到分割，而不是喉道區域的位置被分割。（f）為本文算法分割結果圖，可以看到孔隙中的喉道被全部分割出，沒有過分割和欠分割的現象。因為基于拐點的喉道分割算法從喉道定義出發，在分割喉道時只考慮目標的邊界信息，避免了孔隙區域局部極值和算法循環次數不同帶來極大誤差的影響，因此本文算法能有效地進行喉道分割。

圖11 分割對比圖

3 結語

本文提出了一種基于拐點的喉道分割算法，該算法主要通過檢測出圖像中孔隙連通區域的邊界，然后對邊界進行拐點檢測，并通過篩選算法確定有效拐點。最后進行有效拐點對的匹配，從而形成喉道分割線。通過理論分析和實驗表明，該算法具有以下特性：①算法具有很強的魯棒性，可以分割不同形態下的孔隙目標，例如內部存在孔洞情況。②算法準確度高，無論在外輪廓的凹陷處還是內輪廓的凸出處，均能準確地找出喉道位置。③算法計算量小，通過引入鏈碼表示曲率，可以有效減少計算量。