基于模擬退火的多尺度巖心三維圖像融合重建

孫本耀， 滕奇志， 馮俊羲， 李 洋

(四川大學電子信息學院， 成都 610065)

1 引 言

在石油地質研究中，油氣主要儲存在儲層巖心的孔隙中，孔隙的微觀結構直接影響著油氣的儲量及運移能力[1]，因此，研究孔隙微觀結構具有重要的意義.

近年來，采用計算機斷層掃描(CT)獲取巖心序列切片圖像并構建三維圖像的方法逐漸得到應用，但由于成像機理的限制導致CT掃描圖像在成像視域和分辨率之間存在矛盾.在傳統石油地質研究中，一般利用全直徑巖心或柱塞狀巖心開展物性實驗，因此也希望得到相應視域的三維結構圖像，但由于CT圖像分辨率的不足，導致無法完全表征微米級尺寸的小孔隙.為了獲取微米級及以下小尺寸孔隙結構的清晰圖像，只能將巖石切割成幾毫米甚至更小的樣本，這在一定程度上導致了樣本的代表性有所欠缺.另一方面，在巖心分析研究中，通過對巖心薄片的光學顯微鏡成像以及巖石樣本的掃描電鏡成像等，高分辨率二維巖心圖像的獲取相對容易和普遍，且基于二維巖心圖像的三維重建研究已比較成熟[2-5].因此，我們考慮，利用易獲取的高分辨率二維巖心圖像，借鑒三維重建算法，將高分辨率圖像的小孔隙信息與低分辨率CT圖像進行融合，以獲得較大視域的高分辨率圖像.

目前，有不少學者在多尺度巖心圖像融合方面開展了研究，Okabe和Blunt[6]使用了兩組不同尺度的圖像.將低分辨率3D X射線圖像作為基礎，并且隨機地重建2D高分辨率圖像中所示的亞微米級孔隙，結合其中可以觀察到的小規模和大規模信息作為最終的重建結果.Tahmase[7]等人使用來自兩種不同分辨率的2D圖像，分別建立了納米級和微米級模型并疊加，合成了同時具有納米級孔隙和微米級孔隙的3D巖心圖像.這種直接疊加的方法可能使巖心三維結構中部分小尺寸孔隙被大尺寸孔隙完全覆蓋，導致使巖心的孔隙度發生較大偏差.本研究提出了一種基于模擬退火算法將高分辨率二維巖心圖像與低分辨率三維巖心融合重建為高分辨率三維巖心的算法.

2 算法詳解

2.1 模擬退火算法簡介

模擬退火算法是1953年由Metropoli依據金屬冶煉中退火原理提出，1983年由Kirkpatric[8]將其成功應用于求解組合優化問題中.1997年Hazlett[9-10], Yeong和Torquat等人[11]將其應用到三維隨機重建中.模擬退火算法有效地避免了迭代過程中容易出現局部最優化的現象.并且當重建結構的目標函數收斂于二維圖像的目標函數時，重建結構中一定會呈現出目標函數所描述的特征.

2.2 基于模擬退火的多尺度巖心圖像融合重建

算法以低分辨率三維巖心為基礎，低分辨率三維巖心中孔隙為硬數據，基于模擬退火算法，將高分辨率訓練圖像中的小尺寸孔隙結構直接融合重建在低分辨率三維巖心結構中，實現低分辨率巖心到高分辨率三維巖心的融合重建.

2.2.1 巖心圖像預處理 我們先對不同分辨率巖心圖像進行定標，一個像素(體素)的長度稱為點長度，假如低分辨率巖心圖像的點長度是高分辨率巖心圖像點長度的N倍，則需要將低分辨率三維巖心通過插值放大以統一兩者點長度[12]，低分辨率三維巖心中的每個體素最終在融合重建結果中占據N3個體素.

高低分辨率圖像的孔隙尺寸可能有重復，因此首先在低分辨率三維結構中隨機抽取多張巖心圖像，統計其二維巖心圖像的平均孔隙分布情況，去除高分辨率二維巖心圖像中重復含有的大尺寸孔隙結構，以此作為高分辨率訓練圖像，指導低分辨率巖心的融合重建[13]過程.

2.2.2 模型建立 該算法模型的建立主要由解空間、目標函數和初始解三部分組成[14].在本研究中我們采用了兩點相關函數作為目標函數，算法開始迭代的起點為初始解，解空間是在低分辨率三維結構的基礎上根據二維高分辨率訓練圖像中小尺寸孔隙的比例所組成的三維結構.

(1)

式中，E的值越低說明低分辨率三維結構中融合重建的小尺寸孔隙的相關函數分布和高分辨率二維訓練圖像越接近.

2.2.3 狀態產生函數和接受函數 基于模擬退火的多尺度巖心圖像融合重建算法按隨機的方式產生新狀態：將低分辨率三維結構中的孔隙作為硬數據點不做任何改動，在剩余的三維結構中隨機選擇不同相的兩個像素，交換它們的相，得到新狀態，保證三維結構中各種成分所占比例不變.計算新狀態的能量E′和舊狀態的能量E的差值：ΔE=E′ -E，并按Metropolis準則接受新狀態[15].

(2)

其中，T是溫度，其初始值和變化速度由冷卻進度表決定.

2.2.4 實驗流程 為了驗證算法的有效性，需要對同一巖心分別獲取高低分辨率三維圖像，并對低分辨率三維巖心的融合重建結果與高分辨率三維巖心進行定量分析和視覺效果的比較.但在實際中，對同一巖心分別進行高低分辨率掃描并進行準確的配準很難實現.因此，在本驗證實驗中，利用已有的CT掃描高分辨率三維巖心圖像，對其進行下采樣模擬得到與其對應的低分辨率三維巖心圖像.在高分辨率圖像中去掉大尺寸孔隙，在留下小尺寸孔隙的序列圖像中選取二維巖心圖像作為待融合的訓練圖像.融合重建算法流程圖如圖1所示.圖1中具體流程步驟如下.

圖1 基于模擬退火的多尺度巖心圖像融合重建實驗

Fig.1 Multiscale core image fusion reconstruction experiment based on simulated annealing

(1) 在真實CT高分辨率巖心圖像中選取一張高分辨率二維巖心圖像.將高分辨率三維巖心圖像進行下采樣模擬得到低分辨率巖心圖像.

(2) 將低分辨率三維巖心圖像進行插值放大，使其與高分辨率二維巖心圖像的點長度相同.統計低分辨率三維巖心圖像二維切片的孔隙分布，記錄其最小孔隙尺寸.統計二維高分辨率巖心圖像的孔隙分布，只保留孔隙尺寸小于低分辨率圖像最小孔隙尺寸的孔隙，以此作為二維高分辨率訓練圖像.

(3) 將低分辨率三維巖心圖像中的孔隙相作為硬數據，在重建過程中不做改變.根據二維高分辨率訓練圖像的孔隙度在三維低分辨率巖心圖像的巖石相中隨機布點，將其設置為初始狀態，并設定初始溫度.

(4) 隨機選取巖石相和隨機布置的小孔隙相并進行交換，產生新狀態.

(5) 計算新狀態與舊狀態所對應的目標函數差△E，并依據Metropolis準則判斷新狀態是否被接受[16]：若△E≤0，則接受新狀態作為新的當前狀態；若△E>0,則以概率exp(-△E/T)接受新狀態作為當前狀態.

(6) 重復步驟(4)和步驟(5)，直至交換次數大于設定的門限值或三維結構的能量低于設定值[16],高分辨率三維巖心融合重建完成.

3 實驗結果及分析

在本實驗中，樣本1的CT掃描高分辨率三維巖心圖像如圖2(a)所示，分辨率為10 μm/px，尺寸為256 px×256 px×256 px.通過下采樣模擬得到與其對應的CT掃描低分辨率三維巖心如圖2(b)所示.圖2(b)分辨率為20 μm/px，尺寸128 px×128 px×128 px，分辨率與樣本1高分辨率三維巖心相差2倍.樣本2的CT掃描高分辨率三維巖心圖像如圖2(c)所示，分辨率為10 μm/px，尺寸為256 px×256 px×256 px.通過下采樣模擬得到與其對應的CT掃描低分辨率三維巖心如圖2(d)所示.圖2(d)分辨率為40 μm/px，尺寸64 px×64 px×64 px，分辨率與樣本2高分辨率三維巖心相差4倍.

(a)樣本一高分辨率圖像 (b) 樣本一2倍低分辨率圖像 (c) 樣本二高分辨率圖像 (d) 樣本二4倍低分辨率圖像

圖2 巖心三維圖像

Fig.2 3D core image

圖3(a)是在圖2(a)所示樣本一高分辨率巖心圖像中獲取的一張高分辨率二維巖心圖像，分辨率為10 μm/px，尺寸為256 px×256 px，作為樣本一高低分辨率2倍融合重建實驗的待訓練圖像.圖3(b)是在圖2(c)所示樣本二高分辨率巖心圖像中獲取的一張高分辨率二維巖心圖像，分辨率為10 μm/px，尺寸為256 px×256 px，作為樣本二高低分辨率4倍融合重建實驗的待訓練圖像.由于低分辨率巖心圖像和高分辨率二維巖心圖像分辨率不同.因此需要將低分辨率三維巖心圖像進行插值放大，以統一低分辨率三維巖心和高分辨率二維巖心圖像的點長度.

(a)樣本一2倍融合重建實驗待訓練圖像. (b) 樣本二4倍融合重建實驗待訓練圖像

圖3 高分辨率待訓練圖像

Fig.3 High resolution training image

樣本一低分辨率三維巖心圖像二維切片的孔隙分布和高分辨率二維巖心圖像的孔隙分布如圖4(a)所示.根據兩者的孔隙分布對比，在低分辨率巖心二維切片圖像中直徑小于20 μm的孔隙不存在，而在高分辨率巖心二維圖像中，直徑小于20 μm的孔隙數目占總孔隙數目的28.23%.因此，只保留高分辨率二維巖心圖像中直徑小于20 μm的孔隙，并以此作為二維高分辨率訓練圖像.高分辨率訓練圖像如圖5(a)所示.

同理，根據如圖4(b)所示的樣本二中低分辨率三維巖心圖像二維切片的孔隙分布和高分辨率二維巖心圖像的孔隙分布，只保留高分辨率巖心圖像中直徑小于45 μm的孔隙，并以此作為二維高分辨率訓練圖像.高分辨率訓練圖像如圖5(b)所示.

(a)樣本一高低分辨率巖心孔數目頻率對比 (b) 樣本二高低分辨率巖心孔數目頻率對比

圖4 高低分辨率巖心孔數目頻率對比

Fig.4 Frequency comparison of pore number of high and low resolution core images

3.1 巖心視覺效果對比

圖6(a)為樣本一根據本研究提出的算法將2倍低分辨率三維巖心融合重建得到的高分辨率三維巖心圖像，圖6(b)和圖6(c)為融合重建高分辨率巖心以及原始高分辨率三維巖心在相同位置的二維巖心切片圖像.

(a) 2倍融合重建高分辨率巖心 (b) 2倍融合重建高分辨率巖心圖像 (c) CT掃描高分辨率巖心圖像

圖6 樣本一高分辨率巖心圖像

Fig.6 High resolution core image of sample one

圖7(a)為樣本二根據本研究提出的算法將4倍低分辨率三維巖心融合重建得到的高分辨率三維巖心圖像，圖7(b)和圖7(c)為融合重建高分辨率巖心以及原始高分辨率三維巖心在相同位置的二維巖心切片圖像.

(a) 4倍融合重建高分辨率巖心 (b) 4倍融合重建高分辨率巖心圖像 (c) CT掃描高分辨率巖心圖像

圖7 樣本二高分辨率巖心圖像

Fig.7 High resolution core image of sample two

由圖6和圖7可知，在低分辨率三維巖心和高分辨率訓練圖像分辨率相差2倍和4倍的情況下，融合重建的高分辨率三維巖心在三維外觀形態和二維外觀形態上與真實高分辨率巖心相似，不僅能夠較好的重現訓練圖像中的小尺寸孔隙形態特征及分布，并且可以很好的再現真實高分辨率巖心中的孔隙分布情況.

3.2 巖心孔隙分布對比

巖心孔隙的大小和分布狀態對于油層儲集與滲流特性有重要影響[17].由于真實巖心孔隙的形狀極其不規則，在實際研究中，一般利用孔隙分布的等效球直徑來反映孔隙的尺寸大小情況.為了進一步驗證算法的準確性和穩定性，我們又在圖2(a) 樣本一和圖2(c) 樣本二所示高分辨率巖心圖像中分別隨機挑取了4張高分辨率二維巖心切片圖像，并依照上文所述步驟進行實驗.樣本一和樣本二的真實高分辨率巖心孔隙的等效球直徑以及多組融合重建的高分辨率巖心的平均等效球直徑分布如圖8示.

(a)樣本一2倍融合重建實驗孔隙分布對比 (b) 樣本二4倍融合重建實驗孔隙分布對比

圖8 融合重建高分辨率巖心和真實高分辨率巖心孔隙分布對比圖

Fig.8 Comparison of pore distribution between low resolution core images fusion reconstructed high resolution core images and real high resolution core images

如圖8所示，樣本一在低分辨率三維巖心和高分辨率訓練圖像分辨率相差2倍的情況下，融合重建巖心和高分辨率巖心在孔數目頻率上存在少許誤差，樣本二在低分辨率三維巖心和高分辨率訓練圖像分辨率相差4倍的情況下曲線更為接近，誤差更小.

3.3 統計特征函數分布對比

多孔介質的微觀結構十分復雜，在巖心分析中，一般使用兩點相關函數和線性路徑函數等作為描述巖心微觀結構的重要統計參數[18].在本驗證實驗中樣本一的2倍融合重建高分辨巖心的平均兩點相關函數(S2)和平均線性路徑函數(L)等特征函數和真實高分辨率巖心的對比如圖9.

(a)X方向兩點相關函數 (b)Y方向兩點相關函數 (c)Z方向兩點相關函數

(d)X方向線性路徑函數 (e)Y方向線性路徑函數 (f)Z方向線性路徑函數

圖9 樣本一統計參數對比

Fig.9 Comparison of statistical parameters of sample one

(a)X方向兩點相關函數 (b)Y方向兩點相關函數 (e)Y方向線性路徑函數

(c)Z方向兩點相關函數 (d)X方向線性路徑函數 (f)Z方向線性路徑函數

圖10 樣本二統計參數對比

Fig.10 Comparison of statistical parameters of sample two

樣本二的4倍融合重建高分辨巖心的平均兩點相關函數(S2)和平均線性路徑函數(L)等特征函數和真實高分辨率巖心的對比如圖10所示.

由圖9和圖10可知，實驗融合重建的高分辨率巖心和真實高分辨率三維巖心的兩點相關函數和線性路徑函數的分布十分接近.驗證了算法的準確性.

由實驗及結果分析可以看到，本文提出的基于模擬退火的多尺度巖心融合重建算法可以有效的將高分辨率訓練圖像的小孔隙信息與低分辨率CT圖像進行融合，以重建高分辨率巖心，且融合重建結果有效，準確.

4 結 論

本文借鑒三維重建的思想，提出了一種基于模擬退火的多尺度巖心融合重建算法.為了方便實驗的開展和驗證算法的有效性，采用下采樣的方式來得到低分辨率三維結構，并將去除掉大尺寸孔隙的高分辨率巖心圖像的二維切片作為待融合的訓練圖像，通過模擬退火的方法將訓練圖像中的小尺寸孔隙融合重建在低分辨率巖心中，以重建高分辨率巖心.通過定量分析以及視覺效果的比較，本算法可以很好的將高分辨率信息融合重建在低分辨率巖心中，驗證了算法的有效性.