巖石三維圖像重建算法分析

方瑩瑩，滕奇志，何小海，楊曉敏，李征驥

(1.四川大學(xué)電子信息學(xué)院，成都610064;2.四川大學(xué)圖像信息研究所，成都610064)

研究儲集層孔隙空間的幾何形狀和拓?fù)浣Y(jié)構(gòu)對石油地質(zhì)實(shí)驗(yàn)工作有重要的意義。其中一項(xiàng)最重要的應(yīng)用就是對儲集層中的油氣開采以及地質(zhì)科學(xué)的研究分析。而真實(shí)巖石三維微觀結(jié)構(gòu)的直接測量方法，現(xiàn)在均是通過X射線同步斷層顯微成像技術(shù)(XCT)、核磁共振成像技術(shù)等得到，但XCT設(shè)備價(jià)格昂貴，核磁共振成像得到的圖像分辨率不高，造成直接測量技術(shù)得不到廣泛應(yīng)用。故利用巖石樣本的二維薄片圖像進(jìn)行三維微觀結(jié)構(gòu)重建是目前國內(nèi)外研究的一個方向［1］。由于巖石的二維薄片圖像反映了其三維結(jié)構(gòu)固有的統(tǒng)計(jì)特性，利用薄片圖像重建的三維結(jié)構(gòu)，與真實(shí)三維結(jié)構(gòu)在統(tǒng)計(jì)學(xué)上等價(jià)，能反映真實(shí)樣本的特性，可用于儲集層的特性研究，如微觀孔隙結(jié)構(gòu)、滲流特性等［2］。

目前國內(nèi)外學(xué)者研究較多的三維圖像重建方法主要有:模擬退火重建算法(SA)、基于傅里葉變換的高斯域截?cái)嗨惴?FFT)、基于多點(diǎn)地質(zhì)統(tǒng)計(jì)重建算法(MPS)，以及基于這三種算法的一些改進(jìn)算法。本文先對三類重建算法進(jìn)行描述以及重建過程差別、重建時(shí)間進(jìn)行對比分析，最后用三類算法對同一幅CT圖進(jìn)行重建，并對重建結(jié)果進(jìn)行比較分析，得出綜合性結(jié)論。

1 重建算法

1.1 基于快速傅里葉變換的截?cái)喔咚褂蚍椒?/h3>

該算法基本思想是:利用Wiener－Khinchin定理，通過傅立葉變換生成具有與參考圖像相同二階統(tǒng)計(jì)分布的三維結(jié)構(gòu)。由于傅立葉變換的結(jié)果需要利用孔隙度對其截取，因此這一方法又被稱為基于快速傅里葉變換的截取高斯場重建方法。

基于FFT截?cái)喔咚褂蛩惴ǖ闹亟ㄟ^程為［3］: (1)首先計(jì)算二維參考圖像的歸一化兩點(diǎn)自相關(guān)函數(shù)和孔隙度，其中，二維參考圖像的歸一化兩點(diǎn)自相關(guān)函數(shù)是通過快速傅里葉變換得到;(2)將得到的歸一化自相關(guān)函數(shù)轉(zhuǎn)換為只與r相關(guān)的一維數(shù)組，并進(jìn)行校正，經(jīng)映射得到三維空間的自相關(guān)函數(shù)，此時(shí)三維空間自相關(guān)函數(shù)仍是與方向無關(guān)的一維數(shù)組;(3)根據(jù)對稱性，獲得三維空間三個方向的相關(guān)函數(shù);(4)根據(jù)Wiener－Khinchin定理，對三維自相關(guān)函數(shù)做快速傅里葉變換，得到實(shí)數(shù)域的三維空間高斯場;(4)依據(jù)二維參考圖像的孔隙度對該高斯場進(jìn)行截?cái)啵玫接蓛上嘟M成的(孔隙和顆粒)三維結(jié)構(gòu)。

文獻(xiàn)［2］提出了將粒子群優(yōu)化算法(PSO)引入FFT算法中，目的是在不斷的尋優(yōu)過程中得到與二維參考圖像最接近的三維空間的兩點(diǎn)相關(guān)分布，使得重建過程最優(yōu)。

1.2 基于模擬退火的重建算法

模擬退火算法的基本思想是［4］:通過不斷變換兩個不同相的點(diǎn)，尋找全局最優(yōu)的重建結(jié)果。

重建過程為:(1)計(jì)算二維參考圖像的歸一化兩點(diǎn)自相關(guān)函數(shù)(或者線性路徑等統(tǒng)計(jì)特性參數(shù))和孔隙度，并根據(jù)孔隙度隨機(jī)生成具有相同孔隙度的三維結(jié)構(gòu);(2)隨機(jī)尋找位于不同相位(孔隙相和巖石相)的兩點(diǎn)進(jìn)行交換，計(jì)算交換前后系統(tǒng)的能量;(3)根據(jù)Metropolis準(zhǔn)則判斷是否接受新生成的結(jié)構(gòu)，如果接受，則用新結(jié)構(gòu)代替原結(jié)構(gòu)，否則保留原結(jié)構(gòu);(4)不斷重復(fù)步驟(2)和步驟(3)，直至三維結(jié)構(gòu)的能量低于設(shè)定值，或者交換次數(shù)大于設(shè)定的門限值。

1.3 基于多點(diǎn)地質(zhì)統(tǒng)計(jì)的重建方法

MPS重建算法的基本思想是:使用參考圖像把先驗(yàn)?zāi)Ｐ兔鞔_而定量地引入到建模當(dāng)中，其最終模擬結(jié)果是由參考圖像中的概率信息決定［5－6］。

MPS重建過程［7］是:(1)利用設(shè)定模板掃描二維參考圖像，建立動態(tài)搜索樹;(2)進(jìn)行隨機(jī)布點(diǎn)得到初始結(jié)構(gòu);(3)對每個待模擬點(diǎn)，在搜索樹中搜取該待模擬點(diǎn)的條件數(shù)據(jù)概率，以此決定該待模擬點(diǎn)屬于孔隙相或者巖石相;(4)遍歷所有待模擬點(diǎn)，重復(fù)步驟(3)，直至所有待模擬點(diǎn)被準(zhǔn)確布點(diǎn)［8］。

2 三種重建方法分析

重建時(shí)間方面:(1)FFT算法在計(jì)算二維參考圖像兩點(diǎn)相關(guān)函數(shù)以及獲得三維空間高斯域時(shí)均采用了快速傅里葉變換，重建過程以秒為單位，重建時(shí)間最快;(2)PSO算法作為FFT算法優(yōu)化算法，基本過程和FFT算法一致，但由于增加了多次找尋與參考圖像最接近的自相關(guān)函數(shù)的步驟，所以其重建時(shí)間大幅增加，但重建時(shí)間仍少于SA和MPS算法;(3)MPS算法，重建過程中采用了搜索樹的優(yōu)化方法，將重建時(shí)間縮短，但重建過程中采用了三重模板進(jìn)行重建，導(dǎo)致重建時(shí)間增加;(4)由于SA算法是個尋找全局最優(yōu)的過程，在整個過程中，對每次交換前后的系統(tǒng)能量均需進(jìn)行計(jì)算比較，為了重建結(jié)果達(dá)到最優(yōu)，會在合適范圍內(nèi)選擇盡可能更低的溫度，以及更多的交換次數(shù)，故造成重建過程耗時(shí)較長。

重建的誤差分析:(1)FFT重建過程中，為了加快重建速度，首先通過傅里葉變換得到二維圖像的自相關(guān)函數(shù)，而其重建結(jié)果的三維自相關(guān)函數(shù)也是通過快速傅里葉變換得到的，但快速傅里葉變換對重建結(jié)構(gòu)造成影響，引入一定誤差;接著，由于重建得到的三維空間的函數(shù)是服從正太分布的連續(xù)實(shí)數(shù)，所以必須根據(jù)孔隙度對此連續(xù)實(shí)數(shù)進(jìn)行截?cái)啵拍艿玫椒弦蟮膬上嗳S結(jié)構(gòu)，這使得重建結(jié)果的兩相相對集中，導(dǎo)致結(jié)果出現(xiàn)誤差。(2)PSO算法也會進(jìn)行截?cái)啵蕰霈F(xiàn)結(jié)果偏差情況，但其三維結(jié)構(gòu)的兩點(diǎn)自相關(guān)函數(shù)是一個尋優(yōu)過程的最優(yōu)解，所以重建結(jié)果較FFT算法的形態(tài)以及連通性上更接近于真實(shí)結(jié)構(gòu)。(3)理論上，SA算法是一種全局最優(yōu)搜索算法，在約束函數(shù)的約束下，經(jīng)過多次不同相位的兩點(diǎn)交換、能量更新迭代后可以將初始的隨機(jī)孔隙空間結(jié)構(gòu)優(yōu)化到與真實(shí)三維結(jié)構(gòu)誤差最小的最優(yōu)狀態(tài)。但由于SA算法中，選取的約束函數(shù)只是描述兩點(diǎn)間關(guān)系的低階函數(shù)，比如兩點(diǎn)相關(guān)函數(shù)、線性路徑函數(shù)等，雖然較FFT算法可以采用更多的結(jié)構(gòu)信息，但并不能反映多孔介質(zhì)中多個點(diǎn)之間的關(guān)系情況，所以使得采用SA算法重建得到的三維結(jié)構(gòu)與原始多孔介質(zhì)在連通性上有一定的偏差，不能很好地描述三維孔隙結(jié)構(gòu)的長連通性。(4)MPS算法利用參考圖像中的先驗(yàn)?zāi)Ｐ瓦M(jìn)行重建，將參考圖像中多點(diǎn)間的相關(guān)性引入重建結(jié)構(gòu)，避免了FFT和SA僅利用低階統(tǒng)計(jì)信息而產(chǎn)生誤差的情況，能很好地恢復(fù)孔隙相(或巖性)的幾何形態(tài)(尤其是相邊界)，因此其重建結(jié)果優(yōu)于FFT和SA算法;但由于其重建結(jié)果需要參考圖像的更多信息，故對參考圖像的要求也更高。

3 重構(gòu)實(shí)驗(yàn)結(jié)果分析

在本文中，使用三類重建算法，均是基于二維參考圖片的統(tǒng)計(jì)特性，重建出與參考圖像具有相同統(tǒng)計(jì)特性的三維圖像。為了驗(yàn)證三類重建算法的正確性，我們選擇CT掃描真實(shí)巖心獲得的二維序列圖像中的一幅，作為重建算法的參考圖片，并以序列圖構(gòu)建真實(shí)巖心的三維結(jié)構(gòu)作為標(biāo)準(zhǔn)，同時(shí)，為了定量分析重建效果，采用線性路徑、局部孔隙度和局部滲透性三個分析函數(shù)，對三類算法重建結(jié)果進(jìn)行分析比較。

圖1是選取適合的參考圖像(其中，白色部分表示孔隙)，圖2(a)～圖2(e)分別為CT序列圖、FFT、PSO、SA以及MPS構(gòu)建的1283像素大小的三維孔隙空間結(jié)構(gòu)圖(其中，圖形中顯示部分表示孔隙)。視覺形態(tài)上，發(fā)現(xiàn)SA和PSO算法重建結(jié)果相較于FFT和MPS更接近于真實(shí)三維結(jié)構(gòu)，但是SA算法重建結(jié)果中存在許多孤立的顆粒(即孔隙)，造成連通性被破壞。同時(shí)，發(fā)現(xiàn)FFT算法重建三維結(jié)構(gòu)尺寸明顯大于真實(shí)結(jié)構(gòu)，并且孔隙分布不均勻，較集中。而MPS算法出現(xiàn)很多模糊的邊緣，需要進(jìn)行多次平滑。該結(jié)果與理論分析相一致。

圖1 砂巖二維參考圖像Fig.1 Image of a thin section from sandstone

3.1 數(shù)字巖心線性路徑分析

線性路徑函數(shù)L(u，u+r)指長為r的線段［u，u+r］完全位于孔隙中的概率，是多孔介質(zhì)結(jié)構(gòu)的形態(tài)描述函數(shù)之一，其中包含一定的連通性信息。

圖2 不同方法構(gòu)建的1283像素大小的三維空隙空間結(jié)構(gòu)圖Fig.2 The pore space in a 1283 voxels subregion of samples CT，F(xiàn)FT，PSO，SA and MPS

圖3為CT序列圖和FFT+PSO、SA、和MPS三類算法重建結(jié)果線性路徑函數(shù)的分布曲線。可以看出PSO和SA重建結(jié)果的曲線與CT序列圖基本吻合，MPS和FFT重建結(jié)果與CT序列圖差距較大，其中MPS曲線位于CT序列圖曲線下方，而FFT位于其上方，再次說明MPS重建結(jié)果中孔隙的線性長度要小于CT序列圖，而FFT則大于CT序列圖。

圖3 原始CT序列，F(xiàn)FT，PSO，SA以及MPS線性路徑曲線Fig.3 Linear path function for samples of CT，F(xiàn)FT，PSO，SA and MPS

3.2 數(shù)字巖心均質(zhì)性分析

定義K(O，L)表示一個以O(shè)為多孔介質(zhì)內(nèi)部點(diǎn)為中心、L為邊長的立方體的測量單元，在此測量單元中計(jì)算三維結(jié)構(gòu)局部特征，如局部孔隙度、局部滲透性等。局部孔隙度分布函數(shù):

式中:m為測量單元的個數(shù);δ(x)為狄拉克分布函數(shù)。當(dāng)測量單元長度一定時(shí)，該函數(shù)曲線為孔隙度為φ的測量單元占所有測量單元的比例，它能反應(yīng)三維結(jié)構(gòu)中孔隙度的分布情況，分布越集中，曲線開口越小，三維結(jié)構(gòu)均質(zhì)性越好。

圖4為CT序列圖和FFT+PSO、SA、和MPS三類算法重建結(jié)果的局部孔隙度的分布曲線，取局部測量單元的邊長L=450μm。從圖4中可以看出MPS、SA、PSO曲線的峰值位置對應(yīng)的孔隙度非常接近，是由于峰值位置對應(yīng)的孔隙度為全局孔隙度，與理論一致。FFT曲線與其他四條曲線差異較大，因?yàn)長=450μm時(shí)，某些測量單完全位于FFT重建三維結(jié)構(gòu)的巖石顆粒中，導(dǎo)致測量單元整體的孔隙度較小，再次說明FFT重建結(jié)果中巖石和孔隙尺寸均偏大。同樣，SA曲線和MPS曲線的開口最小，其均質(zhì)性優(yōu)于PSO和FFT重建結(jié)果，而PSO均質(zhì)性優(yōu)于FFT。

圖4 局部孔隙度分布曲線Fig.4 Lacal porosity distributions for sam ples CT，F(xiàn)FT，PSO，SA and M PS

3.3 局部滲透性結(jié)果分析

局部滲透性［9］是衡量某一孔隙度測量單元的連通性能函數(shù)，首先定義滲透特征函數(shù)為:

式中:α的取值為x、y、z、3、c，取3表示三個方向，c表示三個方向任意一個方向。當(dāng)某個α方向上滲透特征函數(shù)為1時(shí)，表示該方向上，測量單元是連通的。

α方向上的局部滲流概率定義為［9］:

由式(3)可以看出，局部滲流概率表示孔隙度為φ的測量單元中，在α方向上連通的測量單元所占的比例。可以發(fā)現(xiàn)當(dāng)特征長度L一定時(shí)，局部滲流概率曲線會隨著局部孔隙度的增大而增大;滲流性越好的三維結(jié)構(gòu)曲線，變化速度越快。

圖5為CT序列圖和三類算法重建結(jié)果的局部滲透性分布曲線，此時(shí)選取XYZ三個方向均連通的情況。對比曲線，可以發(fā)現(xiàn)當(dāng)L=450μm時(shí)，MPS曲線最接近于真實(shí)三維結(jié)構(gòu)，PSO次之，這兩種重建結(jié)果均能較好的反映真實(shí)結(jié)構(gòu)的連通特性;FFT三個方向的連通性最差。

圖5 原始CT序列，F(xiàn)FT，PSO，SA以及MPS局部滲透性曲線Fig.5 Fraction percolation cells versus L for samples of CT，F(xiàn)FT，PSO，SA and MPS

4 結(jié)束語

本文中，實(shí)現(xiàn)FFT(包含優(yōu)化算法PSO)、SA、MPS三類重建算法對同一張真實(shí)的二維CT圖片進(jìn)行三維重建，通過分析比較三類重建結(jié)果與真實(shí)巖心的外觀形態(tài)、低階統(tǒng)計(jì)信息——線性路徑、局部信息——局部孔隙度以及局部滲流概率，得出三類重建算法中MPS和PSO算法，較SA和FFT算法更加接近于真實(shí)巖心的孔隙巖石特征以及滲流特性，更能有效地描述巖石的微觀特性。

［1］Kwiecien M J，Macdonald IF，Dullien F A L.Three－dimensional reconstruction of porous media from serial section data［J］.Journal of Microscopy，2011，159(3): 343－359.

［2］唐棠.基于二維圖像的砂巖結(jié)構(gòu)三.維重建及特性分析［M］.成都:四川大學(xué)，2009.

Tang Tang.Three－dimensional reconstruction and characterization of sandstonemicrostructure from two－dimensional images［M］.Chendu:Sichuan University，2009.

［3］Liang ZR，F(xiàn)ernandes C P，Magnani FS，et al.A reconstruction technique for three－dimensional porous media using image analysis and Fourier transforms［J］.Journal of Petroleum Science and Engineering，1998，21(3－4): 273－283.

［4］De Castro Martins.Image reconstruction using interval simulated annealing in electrical impedance tomography［J］.IEEE Engineering in Medicine and Biology Society， 2012，59(7):1861－1870.

［5］張挺.基于多點(diǎn)地質(zhì)統(tǒng)計(jì)的多孔介質(zhì)重構(gòu)方法及實(shí)現(xiàn)［M］.合肥:中國科學(xué)技術(shù)大學(xué)，2009.

［6］BurcArpat G，Jef Caers.Conditional simulation with patterns［J］.Mathematical Geology，2007，39:2177－203.

［7］HiroshiOkabe，Martin JBlunt.Pore spac reconstruction usingmultiple point statistics［J］.Journal of Petroleum Sisticscience and Engineering，2005，46(1－2):121－137.

［8］黃豐.多孔介質(zhì)模型的三維重構(gòu)研究［M］.合肥:中國科學(xué)技術(shù)大學(xué)，2009.

［9］Jing Hu，Piet Stroeven.Local porosity analysis of pore structure in cement paste［J］.Cement and ConcreteResearch.2005，35(2):233－242.