基于裂縫平移的巖心CT序列圖像裂縫批量提取

許康，滕奇志

（四川大學電子信息學院圖像信息研究所，成都 610065）

0 引言

數字巖心由于其獨特的優勢，近些年來在石油和天然氣行業已經得到了廣泛地運用［1］。建立一個準確客觀的三維數字巖心模型是應用數字巖心技術進行各種數值模擬實驗的基礎［2］，本文所使用到的一種建模的方法是，先通過CT逐層掃描巖心得到序列圖，然后再使用相關軟件批處理序列圖得到處理后的二值圖像，最后再逐層疊加圖層即得到三維的巖心結構。

巖心CT序列圖像的二維批量處理過程是構建巖心三維模型過程中的關鍵步驟，其目的就在于從復雜的巖心結構中提取出感興趣的目標，如巖心的孔隙顆粒和裂縫等，以便于三維模型的構建以及三維參數的計算。巖心當中的微觀裂縫對于研究巖心的連通性和滲流特性具有重要意義［3］［5］。研究某巖心裂縫的相關參數的一般步驟是，先對二維CT序列圖像進行巖心裂縫目標的提取，再通過一系列算法實現數字巖心裂縫三維結構的重建，進而計算裂縫的各類參數。然而目前CT成像的對比度往往不夠理想，且易受到噪聲影響，導致二維批處理提取裂縫困難，且一些細小的裂縫在提取過程中容易斷裂，較難做到完整的提取。

常見的數字巖心裂縫分割提取的辦法一般是先進行圖像增強，再通過閾值分割或區域生長等算法來提取裂縫。但是這類方法往往難以解決CT圖像對比度低和噪聲難以去除等問題，同時也會存在裂縫提取不完整、易斷裂等情況。如果要獲得好的分割效果，通常需要進行大量的人工修補工作。而在實際應用中，巖心CT序列圖的數量通常較大，在批處理過程中如果每一張都需要人為修補，將耗費大量的人力和時間。文獻［3］將CT序列看做一個連續變化的視頻幀序列，利用區域聚類和光流閾值分割進行幀間像素運動的預測來分割裂縫，然而其時間復雜度較高，所以在工程運用中存在局限性。文獻［4］采用相關性閾值分割，先對首幀序列圖計算閾值，再對之后的每張序列圖自動調整分割閾值，以獲得較好的分割效果，然而其主要應用于孔隙提取，對裂縫目標仍然難以做到完整提取。

本文的方法同樣考慮到CT序列幀間較強的相關性［7］，且經過進一步觀察，裂縫在序列圖逐幀遞進的過程中變化較緩慢，裂縫整體呈方向性移動趨勢，因此提出一種基于裂縫平移的巖心CT序列圖像裂縫提取的方法。以第一幀閾值分割預提取的裂縫為模板，將其人為手工補全并去噪，然后分別對圖像中的每條裂縫進行逐幀平移提取，同時根據序列圖中該裂縫的真實形態變化情況來微調平移后的裂縫形態。逐幀平移變化至最后一幀，從而實現序列圖裂縫的批量提取。

1 算法原理

巖心CT序列圖像二維批量處理的過程中存在兩個重要的圖層，其中一個是CT設備掃描巖心之后獲得的原始圖像層，另一個就是用于保存提取之后的孔隙、顆粒和裂縫等目標用于三維重建的二值圖像層，我們大部分的閾值分割和目標提取操作都在該層完成。

本文的算法是將序列圖中的裂縫目標以裂縫平移的方式逐條提取，具體流程如圖1所示。

圖1 裂縫平移提取流程

算法的具體原理如下：針對這種存在方向性移動趨勢的裂縫，通過計算出序列圖像每幀之間裂縫移動的方向以及距離，就可以將上一幀二值圖像層中提取到的裂縫目標移動到下一幀圖層的對應位置上，以實現裂縫目標的粗提取。且由于裂縫幀間差異較小，因此可以模擬圖像層裂縫的真實形態，使粗提取裂縫在保證不斷裂的前提下發生細微的變化。為使變化后的裂縫與裂縫真實情況盡可能相近，因此采取一種部分區域二值化的方法對裂縫周邊區域的像素點進行進一步的提取劃分，同時針對區域二值化過程中可能發生的細微裂縫斷裂情況，需要再結合Dijkstra算法以及區域生長算法進行裂縫的部分區域補全。再對CT圖像中的每條裂縫都執行本文的算法，即可實現CT序列圖像裂縫的批量提取。

1.1 預提取裂縫

先取序列圖的第一幅圖像進行裂縫目標的初步提取，采用常見的閾值分割方法提取即可，通過閾值算法計算出合適的閾值作為二值化圖片的界限，將裂縫目標在二值圖層中置為目標色，其余為背景色。具體情況如圖2 所示，深色部分即為提取到的裂縫目標，它疊加在原始圖像層的上方。

圖2 閾值分割后裂縫斷裂情況

由于裂縫相對背景對比度不高且存在噪點的干擾導致裂縫的提取效果不夠理想，具體表現為：①噪點和裂縫雜糅在一起難以去除。②一些裂縫在提取的過程中由于顏色較淺而出現了斷裂的情況，這會導致裂縫參數的計算出現較大誤差。

因此在閾值分割之后必須對裂縫目標進行人為的修補和去噪。人工交互處理的過程主要是人為手動將斷裂的裂縫連接好，以及將噪點、巖石、孔隙等非裂縫目標置為背景色。

1.2 裂縫平移

由于本文的算法采取逐條提取裂縫的模式，因此可以只針對感興趣的裂縫目標做分割提取，進而排除一些細微的有干擾的裂縫項。

下面是裂縫單條平移提取的具體過程。人工選擇覆蓋其中一條裂縫某段的部分矩形區域如圖3 所示，注意選擇的區域塊中只能包含一條裂縫否則將會出現干擾影響該裂縫的平移提取。選擇的部位盡量為裂縫目標中較直的部分以更好的擬合直線計算方向及距離。

圖3 裂縫平移方向及距離計算

圖3 中的粗線裂縫代表上一幀處理好的裂縫目標，細線裂縫表示裂縫需要平移到下一幀的目標位置，這是由第二幀裂縫的真實位置決定的。由于CT序列圖幀間較強的相關性，同時裂縫存在一定的厚度，所以在實踐中前后兩幀裂縫的位置偏移通常較小且存在重合的部分。

（1）在劃定的矩形區域內將前一幀中粗線表示的裂縫目標像素點提取成一個集合，便于后續的計算操作。先掃描矩形框內的點尋找裂縫目標，以第一個目標色點做為起始點，執行八鄰域的聚類操作。具體過程如下。

定義一個隊列Q，設起始點為p(x,y)，（x,y為其對應的坐標）。并將其放入隊列中，設待提取目標點的集合為P。

1）將p存入隊列Q中。

2）取隊列Q中的首元素，設為x，將該點存入點集P。彈出首元素，同時把x的八鄰域點中為目標色且未訪問過的點存入Q中。

3）當Q不為空時循環執行第二步，直到沒有點滿足條件此時聚類完成。

圖4 裂縫目標的八領域聚類操作

執行聚類操作后裂縫將以坐標數據的形式存儲在集合P中，通過索引坐標值可以獲得該裂縫點位置在巖心原始圖像層的具體灰度值與其在二值圖層的狀態，即目標或背景。

計算該集合像素點的灰度平均值E（P）。其中灰度平均值的計算公式如下：

（2）在當前矩形選框內，以E（P）為閾值進一步分割后一幀圖像，達到粗提取淺色部分裂縫的目標，設后一幀圖像粗提取后的集合為B。

對集合P中的像素點做線性回歸擬合成一條直線，該直線反映了集合P總的數據趨勢［6］。直線方程的最終表達式為hθ(x) =θ0+θ1x, 為了讓擬合成的直線與實際值誤差最小，求得其損失函數為：

解得：

（3）計算集合B到直線hθ(x) =θ0+θ1x在x、y方向的平均移動距離，即可將前一幀的裂縫平行移動相應的距離到下一幀。

每次平移之后的裂縫在位置上基本與前一幀裂縫的位置相一致，但是考慮到CT序列圖像中的真實裂縫在每幀之間的形態肯定會存在細微的差異，因此為了模擬這種形態上的差異，考慮采用一種區域二值化的算法進行局部閾值分割。

1.3 改變裂縫平移形態

由于裂縫存在一定的寬度，所以在線性擬合時生成的直線更趨近于裂縫區域的中軸線。因此在平移裂縫的過程中也會將裂縫移動到下一幀裂縫目標的中軸線附近。考慮在中軸線吻合的情形下，裂縫的形態變化就基本只集中在其外層的部分區域，即二維點集的周邊區域。因此可以保持裂縫目標中心區域不變化，只改變周邊區域就可以模擬出裂縫逐幀形態的變化過程。

一條裂縫在巖心中的形態變化具體包括向外圍區域擴張或者向內部區域收縮（裂縫中途斷裂的情況暫不考慮），因此將裂縫目標區域內靠外一層的點集與裂縫區域最外圍一圈點集作為其大概率發生形態變化的區域。

假定裂縫的主要形態變化只集中在這些區域，如圖5所示。

圖5 裂縫周圍部分二值化區域

最靠近中心的一圈白色區域為裂縫的中心部分，中間層的灰色部分為裂縫外圈的點集，最外圈一層為背景區域。因此我們保留中心部分的狀態，同時合并中間層與最外層，將其作為局部閾值分割的目標區域，對移動后的裂縫進行調整優化。

模擬的灰度圖如圖5，其中灰色與黑色代表裂縫區域。黑色為其中心區域，灰色為裂縫最外圈的部分，白色為背景中緊貼裂縫的一圈區域。

獲得目標區域的具體做法如下：遍歷裂縫目標，若其八鄰域內的點存在不屬于裂縫目標的點，則其為灰色區域的部分，同時其八鄰域中不屬于裂縫目標的點即屬于白色區域。

采用最大類間方差法（OSTU 法）計算目標區域的分割閾值［7］該算法能得到類間分離的最佳分割閾值，具體算法的流程如下：

設圖像的最大灰度級為L-1，統計區域內不同灰度級下像素點的個數n，區域內像素點總數為MN=n0+n1+ …+nL-1，歸一化直方圖分量以及不同灰度等級k下的累積和分別為：

求出使得σ2B(k)取得最大值的閾值k*，此時的k*即為最大類間方差分割閾值。通過該分割閾值對目標區域進行局部二值化分割，可以實現裂縫平移目標的形態變化，具體表現為：將平移而來但是實際不屬于該幀裂縫的像素點置為背景，同時將該幀圖像中裂縫擴展的像素點置為目標。

1.4 基于Dijkstra算法的裂縫自動補全

當裂縫寬度小于3 個像素時，在劃分裂縫區域時就會缺少上文所述的中心區域。在這種情形下，如果該部分目標區域的原始灰度值又較高時，此時經過局部閾值分割之后，裂縫就有可能會出現小部分斷裂的情況。

針對這些小部分斷裂處，需要采取一種基于Dijkstra最短路徑算法［10-12］的裂縫自動補全方法。

設裂縫的兩處斷裂處的端點分別為start 和end，以這兩點為矩形對角線構建一個矩形框，并對該矩形框構造帶權圖G，表示如下：

其中，V代表矩形框中的像素點，即帶權圖中的頂點；E代表圖中的一個邊，圖中的每個頂點（除邊緣處）都有8 條邊與之相交；W代表帶權圖中每條邊的權值，具體對應原始灰度圖像中頂點所代表的像素點的平均灰度值。

Dijkstra最短路徑算法的具體流程如下：

（1）設置斷裂裂縫的兩頭端點為start 和end，以兩點為對角線構造矩形帶權圖G。

（2）創建兩個集合P，Q，其中P為已標識點的最短路徑集合，初始為start 點，對應的最短路徑d為0，Q為未標識點的最短路徑集合，包含除start 點外的所有頂點，若與start 點相鄰則初始值為權值邊w，否則為無窮大。

（3）從Q集合中選取d最小的點，將該點設為x并劃為已標識點，確認點x的所有處于矩形框內的鄰域點yi，使用已知的x至yi的權重邊w[x][yi]，重新計算yi到起始點的最短距離并更新P、Q集合，計算規則是：

d[yi]= min(d[yi],d[x]+w[x][yi])

（4）重復迭代至Q集合為空時結束，此時P集合中end與起始點start的距離即最短路徑。

（5）經過回溯找到start 點和end 點的最短路徑路線，算法結束。

此時將位于最短路徑路線上的點在圖形層中置為裂縫目標，這樣裂縫的小部分斷裂處就實現了自動連接工作。

但是Dijkstra算法只能將斷裂部分以單像素寬的路徑連接成裂縫骨架，連接處與真實裂縫目標仍存在一定的差異，因此需要再通過一種區域生長算法將單像素寬的裂縫骨架橫向擴展，使其與真實的裂縫厚度相一致。

具體過程如下：統計單像素寬的裂縫連接處的所有像素點的平均灰度值G，再將裂縫骨架上的所有目標點當做區域生長的初始種子點，并加入區域生長點集合A中。掃描集合A中的每一個點，若該點的八鄰域范圍內像素點的灰度值g（x，y）滿足如下的關系表達式：

其中T為人為設定的閾值，則將該點在二值圖層中置為裂縫目標，并將該點加入到A中，繼續下一個點的判定。重復循環至A集合為空時，即實現了單像素裂縫的橫向區域生長，此時的裂縫的小部分斷裂處就實現了橫向與縱向的總體自動補全。

2 實驗結果分析

本文對存在裂縫平移趨勢的多組CT序列圖進行了實驗，展示其中的兩組序列圖二值分割的結果如圖6、圖7所示。

圖6 第一組CT序列圖分割結果對比（續）

圖6 第一組CT序列圖分割結果對比

圖7 第二組CT序列圖分割結果對比

將本文算法與固定閾值分割算法與文獻［4］所提及的相關性閾值分割算法進行對比。其中的交互式閾值分割算法即通過人為交互的方式手動設定二值分割圖片的閾值，以實現肉眼可見的最佳分割效果，之后對CT序列圖的每一幀都采取第一幀的固定閾值來進行分割提取。相關性閾值分割則采取首幀固定閾值，之后再對相鄰幀的實際圖像改變分割的閾值，使分割的效果更佳。

從對比圖可以明顯看出，在CT序列圖中不論是采用固定閾值分割法還是相關性分割方法，其提取到的裂縫都會存在較多斷裂處，且裂縫雜糅、形態復雜。

對于第一組序列圖來說，其裂縫本身就較為完整，但是仍然存在部分細微裂縫的情況，本文算法相比固定閾值分割與相關性二值分割來說有較好的提取效果。對于該序列圖中長裂縫目標的局部斷裂處，本文算法成功地實現了完整提取，并延續到之后的每一幀提取上。

對于第二組裂縫存在較多微小裂縫干擾的情況，本文算法先通過人工交互的方式修補連接第一幀斷裂裂縫，去除噪點并刪去干擾的微小裂縫項，簡化裂縫結構，隨后平移提取目標裂縫項，成功實現了裂縫的針對性提取。因此采取本文的方法既可以選擇性地排除微小裂縫的干擾，同時還能夠保持了裂縫提取的完整性。

第二組提取結果的裂縫斷裂處放大圖如圖8所示，可以看到在序列圖100 幀時作為目標的裂縫仍然提取成功且沒有發生斷裂，同時排除了一些較淺或不明顯的其他裂縫的干擾。

圖8 序列圖100幀斷裂裂縫放大處對比

為了驗證CT 序列圖裂縫整體提取的正確性，將采用本文算法提取出的裂縫目標生成三維模型［9］［13］，如圖9所示。

圖9 本文方法三維重建模型結果

可以進一步總結出使用本文方法存在的兩大優勢：

（1）可以針對性地對感興趣的裂縫目標進行提取，進一步排除細小裂縫和干擾裂縫的影響，使重建后的三維模型結構較為清晰；

（2）維持了二維裂縫目標的連續性，使得在批處理提取裂縫不產生局部斷裂，最終重建的三維裂縫模型就會更加地密集且完整。

3 結語

本文依據巖心CT序列圖像中裂縫目標在幀間存在較強的相關性這一特點，提出了一種基于裂縫平移的裂縫批量提取算法。在序列圖首幀進行裂縫的預提取，再通過人工交互的方式補全裂縫，之后再逐條平移提取裂縫，同時還根據裂縫在圖像層的真實灰度值來改變平移裂縫的形態以適應真實的裂縫變化，并實現了基于Dijkstra算法的小部分斷裂裂縫的自動補全，做到了針對性地批處理提取裂縫。

通過對比二維提取結果和三維重建結果，證明了本文算法批量提取裂縫的優越性，在有效去除干擾裂縫及各類噪點的同時還能夠提高裂縫批量提取的完整度，在大批量巖心CT序列圖像的裂縫提取過程中有較好的應用。

本文的方法僅適用于部分裂縫幀間相關性較強的巖心CT圖像。如果單條裂縫在序列圖遞進的過程中分裂成了兩條，或者層數較多裂縫發生改變較大時，采用本文方法就難以做到一次性全部正確分割，需要考慮以第一張裂縫分裂的序列圖為節點前后分別平移提取分割裂縫，這增多了人工交互環節。另外，本文算法默認裂縫最外圍一圈為其大概率發生形態變化的區域，然而當裂縫幀間相差較多時或者裂縫擴展范圍較大時，本文算法就沒有做到很好地還原真實裂縫形態，在這種情形下就應該適當增加局部二值化區域的大小，實現裂縫目標在寬度上的完整提取。