戶瑞林，滕奇志，何小海，龔劍

（1.四川大學電子信息學院圖像信息研究所，成都 610065；2.成都西圖科技公司，成都 610024）

0 引言

近來，由于世界常規(guī)油氣產(chǎn)量的持續(xù)下降，致密氣、頁巖氣、致密油等非常規(guī)油氣資源越來越受到世界各國的關(guān)注[1]。非常規(guī)油氣儲集層以納米級孔喉為主，局部發(fā)育在微米-毫米級，研究儲集層的孔隙結(jié)構(gòu)對于非常規(guī)油氣的開采具有重要意義[2]。

聚焦離子束掃描電鏡(FIB-SEM)是研究非常規(guī)油氣儲集層結(jié)構(gòu)的一種新方法。FIB-SEM使用離子束對巖石進行連續(xù)剝蝕切割，同時利用電子束進行成像，能夠還原真實巖心三維孔隙結(jié)構(gòu)[2]。對巖心進行物相區(qū)分，提取巖心孔隙結(jié)構(gòu)是數(shù)字巖石分析的核心步驟，但在對巖石的FIB-SEM圖像進行孔隙提取時，存在著如下的難點：

（1）FIB-SEM成像時，由于觀察面與電子束不垂直，底部信號比上部弱，掃描圖由上至下會有變暗的效果，這在利用量化方式進行物相區(qū)分時會產(chǎn)生很大的干擾[4]；

（2）在對巖石成像時，由于每張圖片都是分別拍攝，當前幀相對上一幀圖像存在著一定的偏移，這在利用序列圖像層間相關(guān)性進行孔隙提取時，會產(chǎn)生極大的影響；

（3）巖石中除孔隙之外的其他結(jié)構(gòu)，如有機質(zhì)、黏土礦物等，在FIB-SEM圖像中也會呈現(xiàn)深色區(qū)域，這提高了區(qū)分孔隙和巖石的難度；

（4）由于SEM的分辨率較高，當巖石的觀察面出現(xiàn)較大孔隙時，SEM成像會將其內(nèi)部細節(jié)表現(xiàn)出來，再加上荷電的作用，孔隙內(nèi)部一般還伴隨著高亮的特點[4]，這意味著孔隙在圖像中所呈現(xiàn)的形式有所不同，增加了鑒別孔隙的難度。

圖1 致密碳酸巖FIB-SEM圖片

由于各種巖石結(jié)構(gòu)差別很大，F(xiàn)IB-SEM自身成像的特點又對孔隙提取造成干擾，目前尚沒有一種統(tǒng)一的提取巖石孔隙的方法[4]，在實踐中，經(jīng)常需要根據(jù)巖石樣品的特點，使用各種方式嘗試。Martin Salzer[5]等在2012年提出了一種兩階段提取多孔材料FIB-SEM圖像的方法，其利用FIB-SEM圖像孔隙出現(xiàn)的高光效果，提取出高光區(qū)域，然后進行反向傳播，但這是針對多孔材料的分割方法，要求孔隙外的其他部分不會出現(xiàn)明顯灰度的差異，無法對含有多種雜質(zhì)干擾的巖心提取孔隙。王羽[6]等在2016年對頁巖FIB-SEM圖像分別進行了邊緣檢測分割法、流域分割法和手動或自動閾值分割方法的實驗，但效果依賴人為調(diào)整，序列圖中孔隙提取效果表現(xiàn)不佳。

本文利用巖石FIB-SEM圖像具有孔隙存在較為明顯的邊緣，每張圖像變化差異不大的特點，針對孔隙個數(shù)較少的巖石，提出了一種基于主動輪廓的巖心序列圖像孔隙提取方法。我們先在單張圖像上初始化孔隙的邊緣控制點，利用Snake算法收斂到當前孔隙的邊緣，然后將該區(qū)域利用圖像形態(tài)學操作，將邊緣映射到下一張序列圖，繼續(xù)采用Snake算法，進行孔隙提取。由于巖石FIB-SEM圖像中含有黏土礦物等干擾，在各幀圖像進行Snake算法收斂時，部分會出現(xiàn)錯誤分割，本文針對該種情況提出了一種自動尋找最佳分割效果的方法，最終提取出序列圖像孔隙，提取效果良好。

1 預處理

由于FIB-SEM的成像特點，所生成的FIB-SEM圖像會產(chǎn)生位移。我們可以發(fā)現(xiàn)，相鄰兩張巖石FIBSEM圖像的真實觀察面變化很小，若圖像未產(chǎn)生位移，則兩張圖像的差值圖整體灰度應較小，故我們可將后一張圖像逐像素移動，尋找使得其與前一張圖像的差值圖整體灰度最小的位置，則此時應為配準后的圖片。

本文實驗圖為致密碳酸巖的FIB-SEM圖像，我們在圖像的中間設定一個大小為500×500像素的模板窗口，取當前幀圖像在模板窗口位置的部分為I，令為將下一幀沿x方向移動m個像素，沿y方向移動n個像素后模板內(nèi)的圖像，令：

Gm,n表示兩幀的差值圖像的所有像素灰度絕對值之和，當Gm,n最小時，我們認為此時的I'm,n為與I配準的圖像，一般而言FIB-SEM圖像的偏移不會超過50個像素，故我們可以取m∈[- 50,50],n∈[- 50,50]。

我們?nèi)∨錅是昂笮蛄袌D的一個x-z切面觀察，可以發(fā)現(xiàn)，原始圖像的切面圖像連續(xù)性很差，這就是由于圖像間產(chǎn)生了位移導致的結(jié)果，而配準后的圖像明顯比配準前的圖像更加平滑，連續(xù)性更好，這說明了配準的圖像更加接近真實切面的結(jié)構(gòu)。

圖2 配準前后的x-z面圖像

由于后續(xù)操作是需要利用邊緣來提取孔隙，我們需要對配準后的圖像進行邊緣增強操作。FIB-SEM圖像存在一些小噪點，我們先進行中值濾波，然后再求圖像的梯度圖像，簡單起見，我們?nèi)√荻葓D像：

其中g(shù)x，gy為圖像分別在x，y方向上的微分，將梯度圖像疊加到原圖，這樣就得到了增強后的圖像。

圖3 增強前后的圖像

2 FIB-SEM孔隙提取

2.1 Snake模型

Snake模型是Kass[7]在1987年提出的一種主動輪廓線的圖像分割算法，其定義了一條受到內(nèi)、外部約束力的輪廓曲線，通過輪廓曲線的變形，讓其能量極小化，最終使得曲線逼近目標的邊緣。Snake模型由一組控制點：

組成，其中x(s)和y(s)是控制點的坐標位置，s是以傅里葉變換形式描述邊界的自變量。輪廓曲線的能量函數(shù)定義為：

一個作家創(chuàng)作觀的形成，是一個長期積淀的過程，他的生活經(jīng)歷對他的個性、氣質(zhì)、思維方式的形成和發(fā)展起著決定性的作用。“特別是那些印象深刻的經(jīng)驗往往給藝術(shù)家的一生涂上一種特殊的基調(diào)和底色，并在相當程度上決定著藝術(shù)家對于創(chuàng)作題材的選擇和作品的情緒情感或情節(jié)?！盵8]

在該定義中，Eint為內(nèi)部約束能量，Eimage為圖像力，Econ為外部約束力，而內(nèi)部約束能量可以表示為：

在極限的情況下，彈性能量會把輪廓線壓縮成一個點，彎曲能量會把輪廓線變?yōu)橐粋€光滑的曲線或直線[8]，圖像力會驅(qū)使各控制點到其附近梯度最大的位置。通過移動各控制點的位置使總能量最小化的過程中，輪廓線就會逐漸逼近目標的邊緣。Kass提出通過求解歐拉方程來使得總能量最小化，但該方法需要計算高階導數(shù)，計算量較大，對于FIB-SEM序列圖像的處理速度較慢。故本文采用了貪婪算法，每次移動一個控制點，將移動控制點到使得能量函數(shù)最小的位置，然后再移動下一個控制點，如此循環(huán)移動，最終可使得輪廓線逼近孔隙邊緣。

2.2 Snake模型提取單張圖像孔隙

圖4 主動輪廓提取孔隙

由于各類巖心的圖像相差很大，孔隙結(jié)構(gòu)也大有不同，甚至部分孔隙人眼也無法識別，所以第一張種子圖像的輪廓我們由人為給定，如圖4，在孔隙周圍設定控制點，然后自動插值，增加控制點個數(shù)，以盡可能逼近孔隙輪廓。然后Snake模型進行收縮，最終逼近孔隙真實輪廓。

2.3 提取序列圖像的孔隙

由于FIB-SEM圖像每幀間差異變化不大，而Snake算法會自動收縮至邊緣，故我們可以考慮將當前提取完成的區(qū)域進行膨脹，然后將其輪廓映射到下一幀圖像作為其初始輪廓。設我們當前幀所提取的孔隙區(qū)域為M，膨脹結(jié)構(gòu)元為B，下一幀圖像的初始輪廓為G'，G'所圍成的區(qū)域為M'，定義⊕為膨脹操作，則M'=M⊕B,G'=boundary(M')，其中boundary為取區(qū)域邊界輪廓。

在多數(shù)情況下該方法能夠獲得較好的結(jié)果，但膨脹使用固定的結(jié)構(gòu)元，存在提取結(jié)果不理想的情況：當黏土礦物等雜質(zhì)靠近孔隙邊緣時，其也具有強邊緣，若初始化區(qū)域膨脹過大，Snake模型在收縮時會受到雜質(zhì)干擾，容易收斂到雜質(zhì)的邊緣上；當孔隙變化較為劇烈時，若初始化區(qū)域膨脹過小，其輪廓覆蓋不夠大，則Snake模型容易收縮到孔隙內(nèi)部，而出現(xiàn)錯誤分割。

圖5 孔隙提取較壞情況

本文提出一種自動尋找孔隙提取最優(yōu)效果的方法，其可以動態(tài)調(diào)整膨脹結(jié)構(gòu)元，使得孔隙提取區(qū)域達到最佳。思路如下，取相鄰兩幀的同一孔隙的分割結(jié)果，其差異區(qū)域為T，T的面積為Ta，由于兩幀之間的大部分孔隙變化差異不大，若Ta越小，我們越傾向于該次分割效果較好。同時，我們?nèi)∠噜弮蓭牟钪祱D像為Tsub，其包含兩幀間的變化信息，由于兩幀之間的變化集中在邊緣部分，取Tsub內(nèi)的T區(qū)域的灰度值之和為S，若S越大，則本次分割變化部分越接近邊緣，我們傾向于本次分割效果越好。由于S和Ta的量綱不同，我們可以將其歸一化后作差得R，若R值越大，則我們傾向于該次分割效果越好。算法具體步驟如下：

（1）設當前幀圖像為I，下一幀為I'，兩幀差值為為圓盤形結(jié)構(gòu)元，n為圓盤半徑，其范圍可根據(jù)具體情況指定，一般可取n∈[ ]1,10。使用Bn對當前提取的孔隙區(qū)域M分別膨脹，得

（4）取Isub內(nèi)的Tn區(qū)域內(nèi)的所有灰度值之和為STn，分別取 ATn，STn在 n∈[ ]1,10的最大值和最小值Amax，Amin，Smax，Smin，歸一化得:

圖6展示了利用上述算法的過程，圖6中，圖(e)為當前幀的提取區(qū)域，圖(f)為下一幀采用結(jié)構(gòu)元半徑為5的提取區(qū)域，其R值為0.18，圖(g)為下一幀采用結(jié)構(gòu)元半徑為7的提取區(qū)域，其R值為0.05，由于R值越大，效果越好，故圖(f)的提取孔隙效果更優(yōu)。視覺上看，圖(f)確實較圖(g)更優(yōu)，其原因是，圖(g)與圖(f)的差異部分對圖(g)而言，其所增加的A遠大于S，所以圖(g)所提取的區(qū)域的R值較小。

3 實驗結(jié)果與分析

為驗證本文所提出的算法的正確性，將本文的算法與閾值分割，Martin Salzer[5]文中所提出的反向傳播算法作比較。閾值分割是常用的圖像分割方法，其原理是設定分割區(qū)間，通過灰度值的差異來進行圖像分割，Martin Salzer[5]文中所提出的反向傳播算法是針對FIBSEM圖像提出的一種特定算法，在分割多孔材料中表現(xiàn)良好。本文以致密碳酸巖作為樣本，分別使用三種算法進行試驗，可以發(fā)現(xiàn)，由于FIB-SEM存在由上至下的灰度變暗的現(xiàn)象，而且孔隙內(nèi)部灰度不均勻，使用閾值分割，容易受外部黏土礦物的灰度影響，表現(xiàn)不好。而反向傳播算法，由于黏土礦物的邊緣也會出現(xiàn)類似孔隙內(nèi)部灰度突變，而誤將其也認為孔隙，提取效果也比較差，而利用本文所提出的算法，觀察孔隙提取結(jié)果，能夠明顯看出是三種算法中表現(xiàn)最好的。我們也可以取其中一個孔隙，人工識別標注該孔隙區(qū)域，與各算法所提取的區(qū)域相異或，這樣可以看出不同的部分。觀察圖7(s)、(t)、(v)，可以發(fā)現(xiàn)本算法所提取的區(qū)域與人工識別區(qū)域差別較小，提取效果良好。

圖6 自動調(diào)整孔隙提取效果

圖7 不同算法的效果對比

4 結(jié)語

本文提出了一種基于主動輪廓的巖心FIB-SEM序列圖像孔隙提取方法，在初始化首張圖像的孔隙輪廓后，利用FIB-SEM圖像的層間相關(guān)性，通過主動輪廓算法自動提取孔隙，同時也給出了一種自動尋找序列圖像中孔隙最優(yōu)提取效果的算法。本文方法優(yōu)于其他對比算法，能夠準確地提取出孔隙區(qū)域，在目前尚沒有一種統(tǒng)一的提取巖石孔隙的方法的情況下，本文所提出的方法對巖心FIB-SEM圖像具有一定的通用性。但本文方法需要人為給出孔隙的初始輪廓，還不能夠完全自動化提取孔隙，對于孔隙較多的巖石而言，工作量還比較大，因此可在這些方面進一步完善。