許康，滕奇志

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

0 引言

數(shù)字巖心由于其獨(dú)特的優(yōu)勢(shì)，近些年來在石油和天然氣行業(yè)已經(jīng)得到了廣泛地運(yùn)用［1］。建立一個(gè)準(zhǔn)確客觀的三維數(shù)字巖心模型是應(yīng)用數(shù)字巖心技術(shù)進(jìn)行各種數(shù)值模擬實(shí)驗(yàn)的基礎(chǔ)［2］，本文所使用到的一種建模的方法是，先通過CT逐層掃描巖心得到序列圖，然后再使用相關(guān)軟件批處理序列圖得到處理后的二值圖像，最后再逐層疊加圖層即得到三維的巖心結(jié)構(gòu)。

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

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

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

1 算法原理

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

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

圖1 裂縫平移提取流程

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

1.1 預(yù)提取裂縫

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

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

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

因此在閾值分割之后必須對(duì)裂縫目標(biāo)進(jìn)行人為的修補(bǔ)和去噪。人工交互處理的過程主要是人為手動(dòng)將斷裂的裂縫連接好，以及將噪點(diǎn)、巖石、孔隙等非裂縫目標(biāo)置為背景色。

1.2 裂縫平移

由于本文的算法采取逐條提取裂縫的模式，因此可以只針對(duì)感興趣的裂縫目標(biāo)做分割提取，進(jìn)而排除一些細(xì)微的有干擾的裂縫項(xiàng)。

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

圖3 裂縫平移方向及距離計(jì)算

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

（1）在劃定的矩形區(qū)域內(nèi)將前一幀中粗線表示的裂縫目標(biāo)像素點(diǎn)提取成一個(gè)集合，便于后續(xù)的計(jì)算操作。先掃描矩形框內(nèi)的點(diǎn)尋找裂縫目標(biāo)，以第一個(gè)目標(biāo)色點(diǎn)做為起始點(diǎn)，執(zhí)行八鄰域的聚類操作。具體過程如下。

定義一個(gè)隊(duì)列Q，設(shè)起始點(diǎn)為p(x,y)，（x,y為其對(duì)應(yīng)的坐標(biāo)）。并將其放入隊(duì)列中，設(shè)待提取目標(biāo)點(diǎn)的集合為P。

1）將p存入隊(duì)列Q中。

2）取隊(duì)列Q中的首元素，設(shè)為x，將該點(diǎn)存入點(diǎn)集P。彈出首元素，同時(shí)把x的八鄰域點(diǎn)中為目標(biāo)色且未訪問過的點(diǎn)存入Q中。

3）當(dāng)Q不為空時(shí)循環(huán)執(zhí)行第二步，直到?jīng)]有點(diǎn)滿足條件此時(shí)聚類完成。

圖4 裂縫目標(biāo)的八領(lǐng)域聚類操作

執(zhí)行聚類操作后裂縫將以坐標(biāo)數(shù)據(jù)的形式存儲(chǔ)在集合P中，通過索引坐標(biāo)值可以獲得該裂縫點(diǎn)位置在巖心原始圖像層的具體灰度值與其在二值圖層的狀態(tài)，即目標(biāo)或背景。

計(jì)算該集合像素點(diǎn)的灰度平均值E（P）。其中灰度平均值的計(jì)算公式如下：

（2）在當(dāng)前矩形選框內(nèi)，以E（P）為閾值進(jìn)一步分割后一幀圖像，達(dá)到粗提取淺色部分裂縫的目標(biāo)，設(shè)后一幀圖像粗提取后的集合為B。

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

解得：

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

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

1.3 改變裂縫平移形態(tài)

由于裂縫存在一定的寬度，所以在線性擬合時(shí)生成的直線更趨近于裂縫區(qū)域的中軸線。因此在平移裂縫的過程中也會(huì)將裂縫移動(dòng)到下一幀裂縫目標(biāo)的中軸線附近?？紤]在中軸線吻合的情形下，裂縫的形態(tài)變化就基本只集中在其外層的部分區(qū)域，即二維點(diǎn)集的周邊區(qū)域。因此可以保持裂縫目標(biāo)中心區(qū)域不變化，只改變周邊區(qū)域就可以模擬出裂縫逐幀形態(tài)的變化過程。

一條裂縫在巖心中的形態(tài)變化具體包括向外圍區(qū)域擴(kuò)張或者向內(nèi)部區(qū)域收縮（裂縫中途斷裂的情況暫不考慮），因此將裂縫目標(biāo)區(qū)域內(nèi)靠外一層的點(diǎn)集與裂縫區(qū)域最外圍一圈點(diǎn)集作為其大概率發(fā)生形態(tài)變化的區(qū)域。

假定裂縫的主要形態(tài)變化只集中在這些區(qū)域，如圖5所示。

圖5 裂縫周圍部分二值化區(qū)域

最靠近中心的一圈白色區(qū)域?yàn)榱芽p的中心部分，中間層的灰色部分為裂縫外圈的點(diǎn)集，最外圈一層為背景區(qū)域。因此我們保留中心部分的狀態(tài)，同時(shí)合并中間層與最外層，將其作為局部閾值分割的目標(biāo)區(qū)域，對(duì)移動(dòng)后的裂縫進(jìn)行調(diào)整優(yōu)化。

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

獲得目標(biāo)區(qū)域的具體做法如下：遍歷裂縫目標(biāo)，若其八鄰域內(nèi)的點(diǎn)存在不屬于裂縫目標(biāo)的點(diǎn)，則其為灰色區(qū)域的部分，同時(shí)其八鄰域中不屬于裂縫目標(biāo)的點(diǎn)即屬于白色區(qū)域。

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

設(shè)圖像的最大灰度級(jí)為L(zhǎng)-1，統(tǒng)計(jì)區(qū)域內(nèi)不同灰度級(jí)下像素點(diǎn)的個(gè)數(shù)n，區(qū)域內(nèi)像素點(diǎn)總數(shù)為MN=n0+n1+ …+nL-1，歸一化直方圖分量以及不同灰度等級(jí)k下的累積和分別為：

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

1.4 基于Dijkstra算法的裂縫自動(dòng)補(bǔ)全

當(dāng)裂縫寬度小于3 個(gè)像素時(shí)，在劃分裂縫區(qū)域時(shí)就會(huì)缺少上文所述的中心區(qū)域。在這種情形下，如果該部分目標(biāo)區(qū)域的原始灰度值又較高時(shí)，此時(shí)經(jīng)過局部閾值分割之后，裂縫就有可能會(huì)出現(xiàn)小部分?jǐn)嗔训那闆r。

針對(duì)這些小部分?jǐn)嗔烟?，需要采取一種基于Dijkstra最短路徑算法［10-12］的裂縫自動(dòng)補(bǔ)全方法。

設(shè)裂縫的兩處斷裂處的端點(diǎn)分別為start 和end，以這兩點(diǎn)為矩形對(duì)角線構(gòu)建一個(gè)矩形框，并對(duì)該矩形框構(gòu)造帶權(quán)圖G，表示如下：

其中，V代表矩形框中的像素點(diǎn)，即帶權(quán)圖中的頂點(diǎn)；E代表圖中的一個(gè)邊，圖中的每個(gè)頂點(diǎn)（除邊緣處）都有8 條邊與之相交；W代表帶權(quán)圖中每條邊的權(quán)值，具體對(duì)應(yīng)原始灰度圖像中頂點(diǎn)所代表的像素點(diǎn)的平均灰度值。

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

（1）設(shè)置斷裂裂縫的兩頭端點(diǎn)為start 和end，以兩點(diǎn)為對(duì)角線構(gòu)造矩形帶權(quán)圖G。

（2）創(chuàng)建兩個(gè)集合P，Q，其中P為已標(biāo)識(shí)點(diǎn)的最短路徑集合，初始為start 點(diǎn)，對(duì)應(yīng)的最短路徑d為0，Q為未標(biāo)識(shí)點(diǎn)的最短路徑集合，包含除start 點(diǎn)外的所有頂點(diǎn)，若與start 點(diǎn)相鄰則初始值為權(quán)值邊w，否則為無窮大。

（3）從Q集合中選取d最小的點(diǎn)，將該點(diǎn)設(shè)為x并劃為已標(biāo)識(shí)點(diǎn)，確認(rèn)點(diǎn)x的所有處于矩形框內(nèi)的鄰域點(diǎn)yi，使用已知的x至yi的權(quán)重邊w[x][yi]，重新計(jì)算yi到起始點(diǎn)的最短距離并更新P、Q集合，計(jì)算規(guī)則是：

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

（4）重復(fù)迭代至Q集合為空時(shí)結(jié)束，此時(shí)P集合中end與起始點(diǎn)start的距離即最短路徑。

（5）經(jīng)過回溯找到start 點(diǎn)和end 點(diǎn)的最短路徑路線，算法結(jié)束。

此時(shí)將位于最短路徑路線上的點(diǎn)在圖形層中置為裂縫目標(biāo)，這樣裂縫的小部分?jǐn)嗔烟幘蛯?shí)現(xiàn)了自動(dòng)連接工作。

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

具體過程如下：統(tǒng)計(jì)單像素寬的裂縫連接處的所有像素點(diǎn)的平均灰度值G，再將裂縫骨架上的所有目標(biāo)點(diǎn)當(dāng)做區(qū)域生長(zhǎng)的初始種子點(diǎn)，并加入?yún)^(qū)域生長(zhǎng)點(diǎn)集合A中。掃描集合A中的每一個(gè)點(diǎn)，若該點(diǎn)的八鄰域范圍內(nèi)像素點(diǎn)的灰度值g（x，y）滿足如下的關(guān)系表達(dá)式：

其中T為人為設(shè)定的閾值，則將該點(diǎn)在二值圖層中置為裂縫目標(biāo)，并將該點(diǎn)加入到A中，繼續(xù)下一個(gè)點(diǎn)的判定。重復(fù)循環(huán)至A集合為空時(shí)，即實(shí)現(xiàn)了單像素裂縫的橫向區(qū)域生長(zhǎng)，此時(shí)的裂縫的小部分?jǐn)嗔烟幘蛯?shí)現(xiàn)了橫向與縱向的總體自動(dòng)補(bǔ)全。

2 實(shí)驗(yàn)結(jié)果分析

本文對(duì)存在裂縫平移趨勢(shì)的多組CT序列圖進(jìn)行了實(shí)驗(yàn)，展示其中的兩組序列圖二值分割的結(jié)果如圖6、圖7所示。

圖6 第一組CT序列圖分割結(jié)果對(duì)比（續(xù)）

圖6 第一組CT序列圖分割結(jié)果對(duì)比

圖7 第二組CT序列圖分割結(jié)果對(duì)比

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

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

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

對(duì)于第二組裂縫存在較多微小裂縫干擾的情況，本文算法先通過人工交互的方式修補(bǔ)連接第一幀斷裂裂縫，去除噪點(diǎn)并刪去干擾的微小裂縫項(xiàng)，簡(jiǎn)化裂縫結(jié)構(gòu)，隨后平移提取目標(biāo)裂縫項(xiàng)，成功實(shí)現(xiàn)了裂縫的針對(duì)性提取。因此采取本文的方法既可以選擇性地排除微小裂縫的干擾，同時(shí)還能夠保持了裂縫提取的完整性。

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

圖8 序列圖100幀斷裂裂縫放大處對(duì)比

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

圖9 本文方法三維重建模型結(jié)果

可以進(jìn)一步總結(jié)出使用本文方法存在的兩大優(yōu)勢(shì)：

（1）可以針對(duì)性地對(duì)感興趣的裂縫目標(biāo)進(jìn)行提取，進(jìn)一步排除細(xì)小裂縫和干擾裂縫的影響，使重建后的三維模型結(jié)構(gòu)較為清晰；

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

3 結(jié)語

本文依據(jù)巖心CT序列圖像中裂縫目標(biāo)在幀間存在較強(qiáng)的相關(guān)性這一特點(diǎn)，提出了一種基于裂縫平移的裂縫批量提取算法。在序列圖首幀進(jìn)行裂縫的預(yù)提取，再通過人工交互的方式補(bǔ)全裂縫，之后再逐條平移提取裂縫，同時(shí)還根據(jù)裂縫在圖像層的真實(shí)灰度值來改變平移裂縫的形態(tài)以適應(yīng)真實(shí)的裂縫變化，并實(shí)現(xiàn)了基于Dijkstra算法的小部分?jǐn)嗔蚜芽p的自動(dòng)補(bǔ)全，做到了針對(duì)性地批處理提取裂縫。

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

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