戶瑞林,滕奇志,何小海,龔劍
(1.四川大學電子信息學院圖像信息研究所,成都 610065;2.成都西圖科技公司,成都 610024)
近來,由于世界常規(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)錯誤分割,本文針對該種情況提出了一種自動尋找最佳分割效果的方法,最終提取出序列圖像孔隙,提取效果良好。
由于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 增強前后的圖像
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)移動,最終可使得輪廓線逼近孔隙邊緣。
圖4 主動輪廓提取孔隙
由于各類巖心的圖像相差很大,孔隙結(jié)構(gòu)也大有不同,甚至部分孔隙人眼也無法識別,所以第一張種子圖像的輪廓我們由人為給定,如圖4,在孔隙周圍設定控制點,然后自動插值,增加控制點個數(shù),以盡可能逼近孔隙輪廓。然后Snake模型進行收縮,最終逼近孔隙真實輪廓。
由于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值較小。
為驗證本文所提出的算法的正確性,將本文的算法與閾值分割,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 不同算法的效果對比
本文提出了一種基于主動輪廓的巖心FIB-SEM序列圖像孔隙提取方法,在初始化首張圖像的孔隙輪廓后,利用FIB-SEM圖像的層間相關(guān)性,通過主動輪廓算法自動提取孔隙,同時也給出了一種自動尋找序列圖像中孔隙最優(yōu)提取效果的算法。本文方法優(yōu)于其他對比算法,能夠準確地提取出孔隙區(qū)域,在目前尚沒有一種統(tǒng)一的提取巖石孔隙的方法的情況下,本文所提出的方法對巖心FIB-SEM圖像具有一定的通用性。但本文方法需要人為給出孔隙的初始輪廓,還不能夠完全自動化提取孔隙,對于孔隙較多的巖石而言,工作量還比較大,因此可在這些方面進一步完善。