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

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

0 引言

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

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

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

1 喉道分割算法

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

圖1 孔隙示例圖

1.1 獲得目標(biāo)區(qū)域的邊界

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

圖2 八鄰域掩膜

1.2 邊界曲線的拐點(diǎn)檢測

（1）曲線的點(diǎn)曲率

邊界曲線進(jìn)行拐點(diǎn)檢測需對曲線上的點(diǎn)進(jìn)行曲率計算，然后根據(jù)曲率判斷拐點(diǎn)位置。

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

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

圖3 曲率示意圖

（2）鏈碼

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

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

圖4 八連通域鏈碼

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

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

（3）拐點(diǎn)檢測算法

本文基于鏈碼的拐點(diǎn)檢測算法如下：

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

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

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

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

為了進(jìn)一步擴(kuò)大拐點(diǎn)和非拐點(diǎn)處曲率的差距，所以使用以下公式，得到最終曲率值：

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

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

孔隙連通區(qū)域內(nèi)部存在內(nèi)部孔洞時，內(nèi)部孔洞目標(biāo)的邊界為內(nèi)輪廓，孔隙連通區(qū)域邊界為外輪廓。對內(nèi)輪廓取其凸點(diǎn)，經(jīng)由以下篩選和匹配等步驟后與外輪廓上的凹點(diǎn)形成喉道分割線，用以進(jìn)行有效的喉道分割。

1.3 篩選算法

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

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

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

圖5 由于太“淺”不能作為拐點(diǎn)

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

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

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

1.4 匹配算法

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

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

（1）喉道線的有效凹點(diǎn)對應(yīng)該是位于孔隙連通區(qū)域兩側(cè)，如圖6所示。

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

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

圖7 最短距離匹配

A1點(diǎn)對應(yīng)候選的匹配點(diǎn)存在A2、B2兩個點(diǎn)，計算A1A2和A1B2的歐氏距離，取距離最短的A2進(jìn)行匹配，歐氏距離公式如下：

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

其中C為孔隙連通區(qū)域等效圓周長，D為兩匹配凹點(diǎn)間距離。

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

（5）對于孔隙連通區(qū)域內(nèi)部存在孔洞目標(biāo)的情況，應(yīng)該將內(nèi)輪廓上檢測到的凸點(diǎn)進(jìn)行篩選，再與外輪廓上的凹點(diǎn)根據(jù)以上四點(diǎn)匹配準(zhǔn)則進(jìn)行有效匹配，形成分割線，如圖8所示。

圖8 最短距離匹配

1.5 奇數(shù)拐點(diǎn)的處理方法

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

圖9 匹配非拐點(diǎn)

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

通過此法，我們可以處理拐點(diǎn)數(shù)為奇數(shù)的情況，也防止了欠分割的情況，使得基于拐點(diǎn)的喉道分割算法更具有實(shí)際應(yīng)用的價值。

2 測試與分析

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

圖10 實(shí)驗(yàn)結(jié)果圖

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

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

圖11 分割對比圖

3 結(jié)語

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