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

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

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

1 引 言

在石油地質(zhì)研究中，油氣主要儲(chǔ)存在儲(chǔ)層巖心的孔隙中，孔隙的微觀結(jié)構(gòu)直接影響著油氣的儲(chǔ)量及運(yùn)移能力[1]，因此，研究孔隙微觀結(jié)構(gòu)具有重要的意義.

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

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

2 算法詳解

2.1 模擬退火算法簡(jiǎn)介

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

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

算法以低分辨率三維巖心為基礎(chǔ)，低分辨率三維巖心中孔隙為硬數(shù)據(jù)，基于模擬退火算法，將高分辨率訓(xùn)練圖像中的小尺寸孔隙結(jié)構(gòu)直接融合重建在低分辨率三維巖心結(jié)構(gòu)中，實(shí)現(xiàn)低分辨率巖心到高分辨率三維巖心的融合重建.

2.2.1 巖心圖像預(yù)處理 我們先對(duì)不同分辨率巖心圖像進(jìn)行定標(biāo)，一個(gè)像素(體素)的長(zhǎng)度稱為點(diǎn)長(zhǎng)度，假如低分辨率巖心圖像的點(diǎn)長(zhǎng)度是高分辨率巖心圖像點(diǎn)長(zhǎng)度的N倍，則需要將低分辨率三維巖心通過插值放大以統(tǒng)一兩者點(diǎn)長(zhǎng)度[12]，低分辨率三維巖心中的每個(gè)體素最終在融合重建結(jié)果中占據(jù)N3個(gè)體素.

高低分辨率圖像的孔隙尺寸可能有重復(fù)，因此首先在低分辨率三維結(jié)構(gòu)中隨機(jī)抽取多張巖心圖像，統(tǒng)計(jì)其二維巖心圖像的平均孔隙分布情況，去除高分辨率二維巖心圖像中重復(fù)含有的大尺寸孔隙結(jié)構(gòu)，以此作為高分辨率訓(xùn)練圖像，指導(dǎo)低分辨率巖心的融合重建[13]過程.

2.2.2 模型建立 該算法模型的建立主要由解空間、目標(biāo)函數(shù)和初始解三部分組成[14].在本研究中我們采用了兩點(diǎn)相關(guān)函數(shù)作為目標(biāo)函數(shù)，算法開始迭代的起點(diǎn)為初始解，解空間是在低分辨率三維結(jié)構(gòu)的基礎(chǔ)上根據(jù)二維高分辨率訓(xùn)練圖像中小尺寸孔隙的比例所組成的三維結(jié)構(gòu).

(1)

式中，E的值越低說明低分辨率三維結(jié)構(gòu)中融合重建的小尺寸孔隙的相關(guān)函數(shù)分布和高分辨率二維訓(xùn)練圖像越接近.

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

(2)

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

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

圖1 基于模擬退火的多尺度巖心圖像融合重建實(shí)驗(yàn)

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

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

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

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

(4) 隨機(jī)選取巖石相和隨機(jī)布置的小孔隙相并進(jìn)行交換，產(chǎn)生新狀態(tài).

(5) 計(jì)算新狀態(tài)與舊狀態(tài)所對(duì)應(yīng)的目標(biāo)函數(shù)差△E，并依據(jù)Metropolis準(zhǔn)則判斷新狀態(tài)是否被接受[16]：若△E≤0，則接受新狀態(tài)作為新的當(dāng)前狀態(tài)；若△E>0,則以概率exp(-△E/T)接受新狀態(tài)作為當(dāng)前狀態(tài).

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

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

在本實(shí)驗(yàn)中，樣本1的CT掃描高分辨率三維巖心圖像如圖2(a)所示，分辨率為10 μm/px，尺寸為256 px×256 px×256 px.通過下采樣模擬得到與其對(duì)應(yīng)的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.通過下采樣模擬得到與其對(duì)應(yīng)的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倍融合重建實(shí)驗(yàn)的待訓(xùn)練圖像.圖3(b)是在圖2(c)所示樣本二高分辨率巖心圖像中獲取的一張高分辨率二維巖心圖像，分辨率為10 μm/px，尺寸為256 px×256 px，作為樣本二高低分辨率4倍融合重建實(shí)驗(yàn)的待訓(xùn)練圖像.由于低分辨率巖心圖像和高分辨率二維巖心圖像分辨率不同.因此需要將低分辨率三維巖心圖像進(jìn)行插值放大，以統(tǒng)一低分辨率三維巖心和高分辨率二維巖心圖像的點(diǎn)長(zhǎng)度.

(a)樣本一2倍融合重建實(shí)驗(yàn)待訓(xùn)練圖像. (b) 樣本二4倍融合重建實(shí)驗(yàn)待訓(xùn)練圖像

圖3 高分辨率待訓(xùn)練圖像

Fig.3 High resolution training image

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

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

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

圖4 高低分辨率巖心孔數(shù)目頻率對(duì)比

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

3.1 巖心視覺效果對(duì)比

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

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

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

Fig.6 High resolution core image of sample one

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

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

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

Fig.7 High resolution core image of sample two

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

3.2 巖心孔隙分布對(duì)比

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

(a)樣本一2倍融合重建實(shí)驗(yàn)孔隙分布對(duì)比 (b) 樣本二4倍融合重建實(shí)驗(yàn)孔隙分布對(duì)比

圖8 融合重建高分辨率巖心和真實(shí)高分辨率巖心孔隙分布對(duì)比圖

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

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

3.3 統(tǒng)計(jì)特征函數(shù)分布對(duì)比

多孔介質(zhì)的微觀結(jié)構(gòu)十分復(fù)雜，在巖心分析中，一般使用兩點(diǎn)相關(guān)函數(shù)和線性路徑函數(shù)等作為描述巖心微觀結(jié)構(gòu)的重要統(tǒng)計(jì)參數(shù)[18].在本驗(yàn)證實(shí)驗(yàn)中樣本一的2倍融合重建高分辨巖心的平均兩點(diǎn)相關(guān)函數(shù)(S2)和平均線性路徑函數(shù)(L)等特征函數(shù)和真實(shí)高分辨率巖心的對(duì)比如圖9.

(a)X方向兩點(diǎn)相關(guān)函數(shù) (b)Y方向兩點(diǎn)相關(guān)函數(shù) (c)Z方向兩點(diǎn)相關(guān)函數(shù)

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

圖9 樣本一統(tǒng)計(jì)參數(shù)對(duì)比

Fig.9 Comparison of statistical parameters of sample one

(a)X方向兩點(diǎn)相關(guān)函數(shù) (b)Y方向兩點(diǎn)相關(guān)函數(shù) (e)Y方向線性路徑函數(shù)

(c)Z方向兩點(diǎn)相關(guān)函數(shù) (d)X方向線性路徑函數(shù) (f)Z方向線性路徑函數(shù)

圖10 樣本二統(tǒng)計(jì)參數(shù)對(duì)比

Fig.10 Comparison of statistical parameters of sample two

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

由圖9和圖10可知，實(shí)驗(yàn)融合重建的高分辨率巖心和真實(shí)高分辨率三維巖心的兩點(diǎn)相關(guān)函數(shù)和線性路徑函數(shù)的分布十分接近.驗(yàn)證了算法的準(zhǔn)確性.

由實(shí)驗(yàn)及結(jié)果分析可以看到，本文提出的基于模擬退火的多尺度巖心融合重建算法可以有效的將高分辨率訓(xùn)練圖像的小孔隙信息與低分辨率CT圖像進(jìn)行融合，以重建高分辨率巖心，且融合重建結(jié)果有效，準(zhǔn)確.

4 結(jié) 論

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