張 耀 滕奇志 卿粼波 吳曉紅

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

0 引 言

隨著世界常規(guī)油氣產(chǎn)量的持續(xù)下降與油氣勘探技術(shù)的發(fā)展，以頁巖氣、致密油、煤層氣等為代表的非常規(guī)油氣能源逐漸引起各國的重視[1]。與常規(guī)油氣資源不同，非常規(guī)油氣儲集層多以納米級孔喉為主，局部發(fā)育為微米-毫米級，因此微納米孔隙空間的結(jié)構(gòu)特征是影響儲集層的重要因素，準(zhǔn)確全面地表征微納米孔隙結(jié)構(gòu)已成為非常規(guī)油氣研究的重要內(nèi)容[1]。三維圖像重建技術(shù)可用于表征巖石內(nèi)部的微觀結(jié)構(gòu)，為研究非常規(guī)油氣提供了有力的技術(shù)支持。

聚焦離子束(FIB)與掃描電鏡(SEM)雙束系統(tǒng)，通過離子束的納米級表面切割和電子束的超高分辨率成像[2]，可達(dá)到精確重建微納米級孔隙的目的。但由于雙束系統(tǒng)的工作特性，其重建后的三維圖像縱向分辨率與橫向分辨率不一致，且大多表現(xiàn)為縱向分辨率偏低，這將對微納米孔隙后續(xù)研究產(chǎn)生較大影響，因此提高FIB-SEM三維圖像縱向分辨率顯得尤為重要。

提高三維圖像縱向分辨率的方法為層間插值，即通過向空間相鄰的兩幅圖像間插入新圖像的方式提高縱向分辨率。目前三維圖像層間插值方法主要有[3-4]：基于灰度插值、基于小波插值、基于配準(zhǔn)插值。

基于灰度的插值方法計(jì)算復(fù)雜度低且易于實(shí)現(xiàn)，如線性插值與三次樣條插值，但插值精度較低。

基于小波的插值方法通過對小波系數(shù)的操作來實(shí)現(xiàn)層間插值，生成的圖像質(zhì)量較高，但小波系數(shù)最多可表示一維點(diǎn)奇異信息，對圖像中重要的線段、輪廓邊界等二維信息的描述不能達(dá)到最優(yōu)[5-6]，對于圖像高頻信息的插值結(jié)果存在較大誤差[6]。Cunha A L D等[7]針對小波變換的局限性，提出了非下采樣輪廓波變換方法NSCT，可準(zhǔn)確捕獲圖像各方向的二維輪廓信息。

由于基于小波的插值方法的局限，王本峰等[8]將Curvelet變換與凸集投影算法相結(jié)合，較好地改善了 FIB-SEM序列圖層間距過大情況，但此算法僅可用于頁巖切片，通用性不佳。

基于配準(zhǔn)的插值方法[9-12]在醫(yī)學(xué)圖像領(lǐng)域應(yīng)用較為廣泛，其利用形態(tài)學(xué)相關(guān)知識，于插值前先對圖像進(jìn)行形態(tài)學(xué)分析[9]，以獲得待插值圖像像素點(diǎn)的形變信息。此類方法依賴于相鄰層間解剖學(xué)結(jié)構(gòu)的形狀變化及灰度變化為連續(xù)的假設(shè)。對于FIB-SEM圖像，由于FIB對巖心樣品縱深切割間隔為納米級，相鄰序列圖的特征結(jié)構(gòu)形狀與灰度的相對變化具有較強(qiáng)的連續(xù)性，因此基于配準(zhǔn)的插值方法適用于FIB-SEM圖像。

當(dāng)前基于配準(zhǔn)的插值算法均存在配準(zhǔn)過程中出現(xiàn)異常點(diǎn)的問題[12]，且由于配準(zhǔn)算法本身未考慮圖像全局相關(guān)性，因此插值后還會存在圖像高頻細(xì)節(jié)信息丟失現(xiàn)象[13-15]。為解決高頻細(xì)節(jié)丟失現(xiàn)象及運(yùn)算耗時(shí)較長的問題，文獻(xiàn)[15]提出了融合多分辨率分析與曲率配準(zhǔn)的層間插值算法，一定程度上提升了插值效率。但由于每層分辨率的配準(zhǔn)結(jié)果均基于上一分辨率的配準(zhǔn)結(jié)果，因此插值效果較大地依賴于最低級分辨率圖像形變場，算法魯棒性較低；且算法須對修改原圖并手動(dòng)調(diào)整分辨率級數(shù)，但即使在最優(yōu)分辨率級數(shù)下，圖像重要細(xì)節(jié)依然會丟失。

對于異常點(diǎn)問題，文獻(xiàn)[12]提出了基于曲面擬合的配準(zhǔn)插值算法，在一定程度上抑制了配準(zhǔn)過程中異常點(diǎn)的生成，但圖像高頻信息丟失現(xiàn)象依然存在，重建的圖像在一定程度上依然存在模糊現(xiàn)象。

針對當(dāng)前基于配準(zhǔn)的插值算法中存在的高頻細(xì)節(jié)丟失現(xiàn)象，本文提出一種針對FIB-SEM圖像的多策略層間配準(zhǔn)與插值算法。該算法首先對原始FIB-SEM序列圖進(jìn)行NSCT分解，在不破壞原圖信息的基礎(chǔ)上將原圖重要高頻特征與低頻概貌剝離，對二者采用不同的插值算法策略進(jìn)行針對性處理。結(jié)合像素點(diǎn)空間位置信息與NSCT系數(shù)，對NSCT低頻子圖采用基于配準(zhǔn)的層間插值算法進(jìn)行插值，對NSCT高頻方向子帶采用多項(xiàng)式擬合算法進(jìn)行插值，最終將新插值圖像的低頻子圖及高頻方向子帶利用逆NSCT生成新切片圖。該算法將NSCT變換與基于配準(zhǔn)的插值算法有機(jī)結(jié)合，并融入多項(xiàng)式擬合方法，充分利用連續(xù)多張序列圖的高頻信息，較好地抑制了配準(zhǔn)過程中出現(xiàn)異常點(diǎn)以及圖像高頻信息缺失的情況，提高了插值效率與生成圖片的質(zhì)量。

1 NSCT

NSCT繼承自Contourlet變換，解決了Contourlet變換中由于下采樣操作而導(dǎo)致的圖像Gibbs效應(yīng)以及低頻頻譜泄漏等問題[8]。NSCT內(nèi)部采用非下采樣塔式濾波器組NSPFB(Nonsubsampled Pyramid Filter Bank)及非下采樣方向?yàn)V波器組NSDPFB(Nonsubsampled Directional Filter Bank)來捕捉二維圖像幾何結(jié)構(gòu)與紋理細(xì)節(jié)，最終以NSCT系數(shù)來表征各方向的圖像輪廓信息。NSCT首先利用NSPFB對圖像進(jìn)行金字塔分解以便捕捉圖像奇異點(diǎn)信息，分解將獲得一個(gè)低頻圖和一個(gè)帶通圖；然后利用NSDFB將帶通圖沿一系列指定方向進(jìn)行分解，提取二維方向性高頻信息(即圖像輪廓信息)[8]。

每一層NSFDB分解將產(chǎn)生一個(gè)非下采樣低通圖b和b與原始圖的差值圖(帶通圖)，對b繼續(xù)降采樣分解可得下一層的低頻圖及帶通圖，如圖1(a)所示。LP分解后可獲得一系列帶通圖與一個(gè)低頻子圖，NSDFB將對每個(gè)帶通圖進(jìn)行處理以提取二維方向性特征，如圖1(b)所示。對于某一張帶通圖，NSDFB可將其分解為2的任意次冪個(gè)方向，即生成2的任意次冪個(gè)高頻方向子帶。

圖1 NSCT

2 結(jié)合NSCT的多策略層間配準(zhǔn)與插值

2.1 基于配準(zhǔn)的NSCT低頻系數(shù)插值方法

針對當(dāng)前基于配準(zhǔn)的插值算法普遍存在的高頻細(xì)節(jié)丟失以及出現(xiàn)異常點(diǎn)的情況，本文算法將NSCT與配準(zhǔn)算法結(jié)合，對圖像低頻概貌與圖像高頻信息“解耦剝離”，以提升后續(xù)步驟中圖像配準(zhǔn)質(zhì)量及輪廓信息的插值精確度。令FIB-SEM切片圖的低頻子圖為αj0，NSCT第j層第k方向高頻子帶圖為βj,k，則原始FIB-SEM圖像素點(diǎn)的NSCT系數(shù)可表示為：

(1)

(2)

式中：k=?x/δx」-1;l=?y/δy」-1;u=x/δx-?x/δx」;v=y/δy?y/δy」;Bi表示第i階B樣條基函數(shù)；δx與δy分別為x方向與y方向的網(wǎng)格間距。如圖2所示，(a)與(c)為相鄰原始圖像，(b)為以(a)為浮動(dòng)圖像并以(c)為參考圖像獲得的形變場；(d)與(f)分別為(a)與(c)的NSCT低頻子圖，(e)為以(d)為浮動(dòng)圖像并以(c)為參考圖像獲得的形變場。由于?i,j決定非剛體形變程度，因此須對Φd進(jìn)行優(yōu)化?；趯IB-SEM圖像具有非剛性形變特性的考慮，本文采用Rueckert等[16]提出的針對非剛性形變的基于歸一化互信息測度與平滑函數(shù)的控制點(diǎn)自適應(yīng)移動(dòng)優(yōu)化方法對形變場進(jìn)行迭代優(yōu)化：

(3)

(4)

▽C=?C(Φd)/?Φd

(5)

Φd=Φd+μ▽C/‖▽C‖

(6)

式(4)中，H代表圖像熵或圖像聯(lián)合熵。式(3)中，D為圖像面積。根據(jù)修正后的Φd重新計(jì)算梯度向量直至其模小于給定閾值結(jié)束。

圖2 相鄰FIB-SEM序列圖配準(zhǔn)

經(jīng)過優(yōu)化后的形變場TΦd表征相鄰圖像中像素的位移變化，位移變化即為位移向量Kn。如圖3所示，(a)為根據(jù)圖2(e)形變場繪制的圖2(d)與圖2(f)的位移向量二維投影示意圖，為方便讀者觀察圖像配準(zhǔn)結(jié)果，圖中位移向量經(jīng)過同比例放大處理。圖3(b)為位移向量示意圖，其中虛直線即表示位移向量。根據(jù)待插值圖像在左右兩圖的空間相對位置參數(shù)η(0<η<1)以及每個(gè)像素的位移向量Kn(n為NSCT低頻子圖上像素的位置)，可確定待插值圖像的NSCT低頻系數(shù)值與位置：

iη=(1-η)id+ηid+1

(7)

(8)

(9)

式中：φ3代表三次樣條插值。

圖3 位移向量示意圖

2.2 NSCT高頻系數(shù)層間插值方法

由于NSCT將圖像低頻概貌與高頻細(xì)節(jié)剝離，對于NSCT高頻子帶圖βj,k，其表征圖像在某一方向的輪廓細(xì)節(jié)信息。由于FIB-SEM圖像對巖心內(nèi)部結(jié)構(gòu)或孔隙輪廓在縱深方向變化具有良好的連續(xù)性，出于對巖心結(jié)構(gòu)縱深連續(xù)性的考慮，本文根據(jù)Hermite多項(xiàng)式擬合插值算法原理，結(jié)合NSCT高頻系數(shù)特征，對待插值圖高頻部分βj,k進(jìn)行處理。由于多項(xiàng)式擬合插值方法存在龍格現(xiàn)象[17]，因此本文采用三次Hermite插值方法進(jìn)行高頻系數(shù)插值。

以兩節(jié)點(diǎn)為例，設(shè)插值節(jié)點(diǎn)集為{x0,x1,…,xn}，相應(yīng)的節(jié)點(diǎn)值集為{f0,f1,…,fn},設(shè)區(qū)間[xk-1,xk]位于[x0,xn]內(nèi)，且兩邊端點(diǎn)值分別為fk與fk-1，則可定義此區(qū)間上的三次Hermite插值函數(shù)Hk(x)如下：

(10)

式中：dk-1與dk為插值函數(shù)Hk(x)在[xk-1,xk]端點(diǎn)處的一階導(dǎo)數(shù)。對于Hk(x)，其由區(qū)間端點(diǎn)函數(shù)值和一階導(dǎo)數(shù)決定。對于區(qū)間[x0,xn]內(nèi)中間節(jié)點(diǎn)xi(i=1,2,…,n-1),其相應(yīng)導(dǎo)數(shù)可用左右相鄰區(qū)段的一階差商加權(quán)方式近似：

(11)

對于區(qū)間端點(diǎn)，無法獲得其雙側(cè)的一階差商，因此無法使用上式方法計(jì)算導(dǎo)數(shù)值。針對此情況，本文采取區(qū)間端點(diǎn)雙側(cè)一階差商相等的策略。

本文針對FIB-SEM圖像高頻系數(shù)插值節(jié)點(diǎn)數(shù)的選取問題，對大量FIB-SEM圖像進(jìn)行了實(shí)驗(yàn)驗(yàn)證對比分析。實(shí)驗(yàn)結(jié)果表明，針對FIB-SEM圖像，其高頻系數(shù)插值節(jié)點(diǎn)數(shù)取4為最佳(即考慮連續(xù)4張F(tuán)IB-SEM圖像的輪廓變化信息)，可充分利用巖心縱向結(jié)構(gòu)輪廓的連續(xù)變化信息。因此本文對NSCT高頻部分的處理均采用4節(jié)點(diǎn)插值(即“四圖策略”)。

結(jié)合上文所述，F(xiàn)IB-SEM圖像多策略層間配準(zhǔn)與插值算法可總結(jié)為以下幾個(gè)步驟：

1) 確定原始FIB-SEM序列圖中新層間圖的空間位置xη,與新層間圖關(guān)聯(lián)的4幅原始圖組成一個(gè)關(guān)聯(lián)組Fxη。

2) 對原始序列圖執(zhí)行NSCT變換，獲得低頻子圖與高頻方向子帶圖。

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

3.1 NSCT低頻子圖配準(zhǔn)效果驗(yàn)證實(shí)驗(yàn)

為驗(yàn)證本文算法低頻子圖配準(zhǔn)部分處理效果，對6組原始FIB-SEM序列圖(圖像均由中石油勘探開發(fā)研究院提供)，分別進(jìn)行原始圖像配準(zhǔn)實(shí)驗(yàn)和NSCT低頻子圖配準(zhǔn)實(shí)驗(yàn)。對于此兩種實(shí)驗(yàn)，其配準(zhǔn)方式均采用以后一張圖像為參考圖并以前一張圖為浮動(dòng)圖的方式進(jìn)行配準(zhǔn)，并將配準(zhǔn)優(yōu)化后的結(jié)果圖與參考圖作比較。實(shí)驗(yàn)結(jié)果采用業(yè)界廣泛認(rèn)可的均方差MSD(Mean Square Difference)與互相關(guān)系數(shù)CC(Cross Correlation)[18]：

(12)

式中：I1為原始切片圖，I2為層間插值后圖像，切片圖像大小為m×n，R(i,j)代表參考圖(i,j)處的像素灰度值，F(xiàn)(i,j)代表浮動(dòng)圖(i,j)處的像素灰度值。MSD值越小，CC值越趨近于1，則圖像配準(zhǔn)精度越高，效果越好。由表1可知，對原始FIB-SEM圖的NSCT低頻圖像配準(zhǔn)，其配準(zhǔn)精度與配準(zhǔn)效果要優(yōu)于直接對原始FIB-SEM圖像配準(zhǔn)。實(shí)驗(yàn)結(jié)果如表1所示。

表1 原始FIB-SEM圖配準(zhǔn)與NSCT低頻圖配準(zhǔn)比較

續(xù)表1

3.2 FIB-SEM圖像配準(zhǔn)與插值實(shí)驗(yàn)

為驗(yàn)證本文算法的正確性和有效性，將本文算法與線性插值法、penny插值法[11]及文獻(xiàn)[12]提出的USFBR算法作比較。線性插值法是被廣泛使用的基于場景插值的方法，penny插值法則為被廣泛認(rèn)可的基于配準(zhǔn)的層間插值算法，而USFBR算法則為近期在層間插值算法研究領(lǐng)域中一種優(yōu)秀的基于配準(zhǔn)的層間插值算法。本文共進(jìn)行6組實(shí)驗(yàn)，實(shí)驗(yàn)中涉及的圖像均為FIB-SEM圖像。實(shí)驗(yàn)結(jié)果如表2所示。

表2 切片圖MSD與PSNR比較

樣本圖像大小均為512×512。為保證實(shí)驗(yàn)數(shù)據(jù)結(jié)果真實(shí)性，輸入均為FIB-SEM原始掃描圖像。實(shí)驗(yàn)結(jié)果評判標(biāo)準(zhǔn)參數(shù)采用業(yè)界廣泛認(rèn)可的均方差MSD與峰值信噪比PSNR：

(13)

6組實(shí)驗(yàn)均選取五張連續(xù)原始圖像的最前兩張與最后兩張來重建中間第三張，并將重建圖與原始圖進(jìn)行比較。為使讀者更加清晰地觀察比較4種方法的插值結(jié)果的區(qū)別，本文在圖5中給出了圖4白色矩形框標(biāo)識區(qū)域的放大圖。

圖4 2組實(shí)驗(yàn)的結(jié)果圖

為客觀比較4種算法，實(shí)驗(yàn)計(jì)算了4種算法生成的各新層間圖的MSD與PSNR，如表2所示。根據(jù)表2結(jié)果可知，本文算法在客觀評價(jià)標(biāo)準(zhǔn)上優(yōu)于其他對比算法，PSNR相較于USFBR算法平均提升0.48 dB以上。

限于篇幅，圖5列出了6組實(shí)驗(yàn)中2組實(shí)驗(yàn)的結(jié)果示意圖：(a)-(e)為頁巖A樣本序列中某一張圖片的重建結(jié)果；(f)-(i)為致密碳酸鹽巖B樣本序列某一張圖的重建結(jié)果(為便于讀者觀察各算法結(jié)果的差異，此樣本重建結(jié)果的放大圖像經(jīng)過亮度增強(qiáng))。根據(jù)圖5所示結(jié)果可以看出，線性插值生成的圖像不僅會出現(xiàn)圖像信息明顯丟失的現(xiàn)象(尤以圖5(b)最為顯著)，而且圖像存在拖影的視覺現(xiàn)象和模糊現(xiàn)象；penny算法相較于線性插值算法，拖影現(xiàn)象明顯減少，孔隙準(zhǔn)確率有所提升，但模糊現(xiàn)象更為嚴(yán)重；USFBR算法相較于penny算法在圖片質(zhì)量上有了進(jìn)一步提升，塊狀異常區(qū)域減少，但模糊現(xiàn)象仍存在；本文算法在USFBR算法基礎(chǔ)上有了進(jìn)一步提升，減小了亮度差異較大區(qū)域與原圖的差異，減少了模糊現(xiàn)象。因此，本文算法是此4種算法中最優(yōu)的。

4 結(jié) 語

本文針對FIB-SEM圖像，提出多策略層間配準(zhǔn)與插值算法。該算法利用NSCT可有效捕獲二維圖像輪廓信息的特點(diǎn)，將圖像高頻輪廓細(xì)節(jié)與低頻概貌進(jìn)行分離。低頻概貌部分結(jié)合基于配準(zhǔn)的插值算法，提升低頻部分配準(zhǔn)與插值的準(zhǔn)確度；高頻輪廓細(xì)節(jié)部分利用融入端點(diǎn)優(yōu)化策略的分段三次Hermite多項(xiàng)式，提升新層間圖像輪廓結(jié)構(gòu)的精確度與質(zhì)量。實(shí)驗(yàn)結(jié)果顯示，本文算法的實(shí)驗(yàn)結(jié)果優(yōu)于各對比算法的實(shí)驗(yàn)結(jié)果與其他對比算法對比，PSNR平均提升0.48 dB以上。但本文算法須對圖像進(jìn)行多層結(jié)構(gòu)分解，并于分解過程中引入非下采樣技術(shù)以減少圖像信息的損失，算法內(nèi)部涉及迭代運(yùn)算，若對原始FIB-SEM序列圖采用串行運(yùn)算方式進(jìn)行處理，則運(yùn)算時(shí)耗大于線性插值，因此本文算法在速度上需進(jìn)一步優(yōu)化。

參 考 文 獻(xiàn)

[1] 賈承造,鄭民,張永峰.中國非常規(guī)油氣資源與勘探開發(fā)前景[J].石油勘探與開發(fā),2012,39(2):129-136.

[2] 孫亮,王曉琦,金旭,等.微納米孔隙空間三維表征與連通性定量分析[J].石油勘探與開發(fā),2016,43(3):490-498.

[3] Ale? N,Olivier S,Oscar A,et al.Constrained reverse diffusion for thick slice Interpolation of 3D volumetric MRI images [J].Computerized Medical Imaging & Graphics,2012,36(2):130-138.

[4] Jafari-Khouzani K.MRI Upsampling Using Feature-Based Nonlocal Means Approach[J].IEEE Transactions on Medical Imaging,2014,33(10):1969-1985.

[5] Chang S G,Cvetkovic Z,Vetterli M.Locally adaptive wavelet-based image interpolation[J].IEEE Transactions on Image Processing,2006,15(6):1471-1485.

[6] 武士想,尚鵬,王立功,等.小波-Lagrange方法進(jìn)行醫(yī)學(xué)圖像層間插值[J].中國圖象圖形學(xué)報(bào),2016,21(1):78-85.

[7] Cunha A L D,Zhou J,Do M N.The nonsubsampled contourlet Transform:Theory,Design,and Applications[J].IEEE Transactions on Image Processing,2006,15(10):3089-3101.

[8] 王本鋒,李景葉,陳小宏,等.基于Curvelet變換與POCS方法的三維數(shù)字巖心重建[J].地球物理學(xué)報(bào),2015(6):2069-2078.

[9] Grevera G J,Udupa J K.An objective comparison of 3D image interpolation methods[J].IEEE Transactions on Medical Imaging,1998,17(4):642-52.

[10] Herman G T,Zheng J,Bucholtz C A.Shape-Based Interpolation[J].IEEE Computer Graphics & Applications,1992,12(3):69-79.

[11] Penney G P,Schnabel J A,Rueckert D,et al.Registration-based interpolation[J].IEEE Transactions on Medical Imaging,2004,23(7):922-6.

[12] 張杰,何小海,卿粼波,等.基于曲面擬合及配準(zhǔn)的醫(yī)學(xué)圖像層間插值[J].四川大學(xué)學(xué)報(bào)(工程科學(xué)版),2016(S1):142-149.

[13] Frakes D H,Dasi L P,Pekkan K,et al.A new method for registration-based medical image interpolation[J].IEEE Transactions on Medical Imaging,2008,27(3):370-377.

[14] Laursen L F,Hansen M S.Registration-based interpolation realtime volume visualization[C]//Spring Conference on Computer Graphics [s.n.],2013:15-21.

[15] 王新征,雄洙,于靖.結(jié)合多分辨率修正曲率配準(zhǔn)的層間插值[J].光學(xué)精密工程,2016,24(5):1224-1231.

[16] Rueckert D,Sonoda L I,Hayes C,et al.Nonrigid registration using free-form deformations: application to breast MR images[J].IEEE Transactions on Medical Imaging,1999,18(8):712-21.

[17] Epperson J F.On the Runge Example[J].American Mathematical Monthly,1987,94(4):329-341.

[18] 高廣珠,李忠武,余理富,等.歸一化互相關(guān)系數(shù)在圖像序列目標(biāo)檢測中的應(yīng)用[J].計(jì)算機(jī)工程與科學(xué),2005(3):38-40.