秦杰1，吳運(yùn)強(qiáng)2，滕奇志1，王子強(qiáng)2，何小海1

　　(1.四川大學(xué) 電子信息學(xué)院 圖像信息研究所，四川 成都 610065；2.中石油新疆油田分公司，新疆 克拉瑪依 834000 )

    摘要：為了提高穩(wěn)相分級(jí)模擬退火算法在三維重建上的精度，引入隨機(jī)線性路徑函數(shù)，與兩點(diǎn)相關(guān)函數(shù)同時(shí)約束重建過(guò)程，借助穩(wěn)相模擬退火思想，最終重建出三維圖像。針對(duì)大量巖心二維圖像進(jìn)行實(shí)驗(yàn)分析，結(jié)果表明：相較于單約束條件，改進(jìn)的重建算法能夠顯著提高重建精度，保存更多原始圖像的形態(tài)學(xué)特性。

　　關(guān)鍵詞：分級(jí)模擬退火；兩點(diǎn)相關(guān)函數(shù)；線性路徑函數(shù)；三維重建

0引言

　　巖心的三維結(jié)構(gòu)分析在石油地質(zhì)巖心研究中有很重要的意義，通過(guò)對(duì)巖心的三維結(jié)構(gòu)分析可以得出其微觀結(jié)構(gòu)特性進(jìn)而可以計(jì)算出其宏觀特性參數(shù)［12］。但是在多數(shù)情況下，獲取目標(biāo)真實(shí)三維結(jié)構(gòu)難度較大［3］，其單張的二維結(jié)構(gòu)圖像卻可以方便得到，所以在二維圖像的基礎(chǔ)上進(jìn)行三維重建有非常大的研究?jī)r(jià)值。

　　本文所使用的模擬退火重建算法［4］最早由Yeong和TORQUATO S教授提出，應(yīng)用于儲(chǔ)集層巖心三維重建領(lǐng)域，但該算法計(jì)算量大、時(shí)間效率低。ALEXANDER S K等人提出分級(jí)模擬退火算法［5］，依次重建整體和細(xì)節(jié)信息，顯著降低了計(jì)算量，有效拓展了模擬退火算法的應(yīng)用范圍。但是在分級(jí)重建中，每一級(jí)依然需要遍歷所有的像素點(diǎn)，存在重復(fù)計(jì)算問(wèn)題。由陳冬冬提出的穩(wěn)相分級(jí)模擬退火算法［6］，針對(duì)分級(jí)策略進(jìn)行改進(jìn)，根據(jù)每級(jí)重建結(jié)果確定穩(wěn)定數(shù)據(jù)，進(jìn)一步降低計(jì)算量。然而，穩(wěn)相分級(jí)模擬退火算法僅使用兩點(diǎn)相關(guān)函數(shù)作為約束，重建時(shí)并未考慮線性關(guān)系。

　　本文針對(duì)該問(wèn)題，引入隨機(jī)線性路徑函數(shù)，與兩點(diǎn)相關(guān)函數(shù)同時(shí)約束重建過(guò)程，利用分級(jí)策略提高時(shí)間效率。對(duì)巖心圖像進(jìn)行對(duì)比實(shí)驗(yàn)，結(jié)果表明相較于單約束條件，改進(jìn)算法顯著提高了重建精度，有效完成了巖心圖像的三維重建。

1穩(wěn)相分級(jí)模擬退火算法描述

　　穩(wěn)相分級(jí)模擬退火算法在原有算法的基礎(chǔ)上，通過(guò)增加新相位，將兩相重建轉(zhuǎn)換為三相重建，使用新添加的相位有效地避免了兩相分級(jí)退火的重復(fù)計(jì)算問(wèn)題，大大提升了收斂速度。

　　分級(jí)過(guò)后，從最高級(jí)圖像開(kāi)始重建，穩(wěn)相分級(jí)采用三相兩點(diǎn)概率作為約束函數(shù)，3個(gè)相位可能產(chǎn)生的9種概率函數(shù)為：P11、P12、P13、P21、P22、P23、P31、P32、P33。根據(jù)其中參數(shù)對(duì)稱(chēng)性與非線性相關(guān)函數(shù)［6］，選取3個(gè)獨(dú)立參數(shù)作為參考對(duì)象，計(jì)算原始圖像與隨機(jī)重建圖像兩點(diǎn)概率的能量誤差，并作為模擬退火時(shí)的能量，方程如下：

　　隨機(jī)交換重建圖像兩點(diǎn)，按照式（1）亦可計(jì)算出新?tīng)顟B(tài)與原始圖像的能量誤差Etp*，相減得到新老狀態(tài)模擬退火時(shí)的能量差ΔE=Etp*-Etp，若能量差小于或等于0則接受新?tīng)顟B(tài)，若大于0則按exp(-ΔET)概率接受新?tīng)顟B(tài)，否則拒絕新?tīng)顟B(tài)。

　　將得到的新?tīng)顟B(tài)按上述步驟循環(huán)迭代，將迭代的最大拒絕次數(shù)設(shè)為200 000，降溫系數(shù)設(shè)為0.99，迭代步數(shù)設(shè)為500，溫度T按0.99速率下降，最小溫度設(shè)為1.0e-38，按照經(jīng)典模擬退火算法［5］重建本級(jí)圖像。完成后，將重建好的圖像按穩(wěn)相規(guī)則［6］向上一級(jí)還原，繼續(xù)重建步驟，直至第0級(jí)重建完成。

2隨機(jī)線性路徑函數(shù)

　　2.1線性路徑函數(shù)概述

　　不同于兩點(diǎn)概率函數(shù)，線性路徑函數(shù)［79］描繪的是將一些隨機(jī)的標(biāo)量線段拋入?yún)⒖紙D像的像素點(diǎn)中，各個(gè)標(biāo)量線段落在特定相位上的概率。在兩點(diǎn)概率函數(shù)中，只是給出了線段首位兩點(diǎn)相位概率分布，而線性路徑函數(shù)彌補(bǔ)了兩點(diǎn)概率的不足，給出了整條線段的相位概率分布。線性路徑函數(shù)的數(shù)學(xué)表達(dá)式如下：

　　其中，Ln（r）表示將長(zhǎng)度為r的標(biāo)量線段隨機(jī)拋入樣本圖片中N次，其落在相位n的概率。ln（r）代表了長(zhǎng)度為r的標(biāo)量線段落在相位n的次數(shù)。顯而易見(jiàn)，線性路徑函數(shù)是一個(gè)關(guān)于r單調(diào)遞減函數(shù)，隨著線段r增長(zhǎng)，其線性路徑概率會(huì)逐漸降低?？紤]一個(gè)極限情況，如果r為0時(shí)，得到的線性路徑函數(shù)就變?yōu)榱舜讼辔徽颊麄€(gè)圖像的比重，即相位密度。如果r趨近于無(wú)窮時(shí)，其線性路徑函數(shù)最終結(jié)果為0，數(shù)學(xué)表達(dá)式如下：

　　Ln（0）=vn,Ln（∞）=0（3）

　　顯然，線性路徑函數(shù)表征了圖像目標(biāo)相的連通特征。與兩點(diǎn)概率函數(shù)不同，不同相位的線性路徑函數(shù)之間是相互獨(dú)立的，并沒(méi)有任何線性關(guān)系。

　　2.2隨機(jī)線性路徑函數(shù)

　　為簡(jiǎn)便計(jì)算，線性路徑函數(shù)一般只統(tǒng)計(jì)目標(biāo)相在固定X、Y方向上的連通特征，如圖1中的L1（r）和L2（r）。然而，這種方法導(dǎo)致線性路徑函數(shù)僅能表征固定方向的連通特征，缺乏對(duì)整體連通特征的表達(dá)，因此，本文提出隨機(jī)線性路徑函數(shù)。設(shè)想如圖2中L3（r）與L4（r）這種隨機(jī)傾斜線段，其方向并未落在X、Y方向上，卻一樣表現(xiàn)了整幅圖像的連通特性，但在以往算法中并沒(méi)有考慮入內(nèi)。本文提出隨機(jī)線性路徑函數(shù)中就包含了這些隨機(jī)方向上的線段。不同于以往線性路徑函數(shù)的計(jì)算，隨機(jī)方向上線性路徑的相位判斷需要特殊的方法，描述如下：假設(shè)需要判定的線段兩端點(diǎn)分別為（X1，Y1）和（X2，Y2），可以得出其方程為：

　　分別從橫坐標(biāo)與縱坐標(biāo)區(qū)間中開(kāi)始遍歷，以橫坐標(biāo)為例，設(shè)X3∈（X1，X2），將其帶入式（4）中得到Y(jié)3坐標(biāo)，若Y3為整型，則直接判斷（X3，Y3）點(diǎn)是否在目標(biāo)相位，若Y3為圖1隨機(jī)線性路徑函數(shù)示意圖浮點(diǎn)數(shù)，則將其分別向上向下取整得到minY3和maxY3，判斷（X3，minY3）與（X3，maxY3）是否同時(shí)在目標(biāo)相位，若是，則證明在上述線段橫坐標(biāo)為X3上的點(diǎn)在目標(biāo)相位上，繼續(xù)往下判斷，直至遍歷結(jié)束；若不是則證明此條線段并不同相，終止循環(huán)遍歷。縱坐標(biāo)的判斷與橫坐標(biāo)類(lèi)似，并且僅當(dāng)橫縱坐標(biāo)上所有點(diǎn)都判斷屬于同一目標(biāo)相位時(shí)才能確定此條線段落在目標(biāo)相位上。

　　2.3隨機(jī)線性路徑函數(shù)在穩(wěn)相重建中的運(yùn)用

　　一般巖心圖像僅包含孔隙與巖石兩相，因此本文算法僅統(tǒng)計(jì)此兩相的隨機(jī)線性路徑，由于穩(wěn)相策略而產(chǎn)生的附加相并不予以統(tǒng)計(jì)。利用線性路徑函數(shù)相互獨(dú)立特性，計(jì)算能量時(shí)也只需要考慮感興趣的兩個(gè)相位，其能量定義如下：

　　結(jié)合式（2），將兩種限制函數(shù)能量相加，可得到退火時(shí)的總能量：

　　E=Etp+Elp(6)

　　將得到的總能量替換之前的能量值，繼續(xù)按照上文描述的步驟重建圖像，這樣隨機(jī)線性路徑函數(shù)就與兩點(diǎn)概率函數(shù)同時(shí)作為限制函數(shù)加入到穩(wěn)相模擬退火重建算法中。

3重建結(jié)果分析

　　將本文算法應(yīng)用于儲(chǔ)集層巖心圖像重建中，進(jìn)行對(duì)比實(shí)驗(yàn)，驗(yàn)證本文算法的真實(shí)性與有效性。在處理過(guò)程中，白色相作為1相位，黑色相作為2相位，灰色相記為3相位。

　　3.1本文重建結(jié)果與原始圖像對(duì)比分析

　　給出尺寸為512×512的訓(xùn)練圖像，如圖2（a）所示。將圖像分為三級(jí)，利用本文重建算法重建后的結(jié)果如圖2（b）所示。圖3（a）和圖3（b）分別給出了重建圖和原始圖的兩點(diǎn)概率和線性路徑函數(shù)的對(duì)比，從圖中可以看出重建結(jié)果與原始圖像的兩點(diǎn)概率和線性路徑的統(tǒng)計(jì)信息非常吻合，證明了本文算法的真實(shí)性。但是因?yàn)閮牲c(diǎn)概率和線性路徑函數(shù)并不能包含原始圖像的所有信息特性，如邊緣信息等，所以在人眼觀測(cè)下重建圖與原始圖依然存在一定的差別。

　　3.2三維重建結(jié)果對(duì)比分析

　　本文針對(duì)尺寸為256×256的訓(xùn)練圖像（如圖4(a)所示），分別使用單約束條件下穩(wěn)相模擬退火算法與本文算法進(jìn)行三維重建，得到結(jié)果如圖4(b)和圖4(c)所示。為了更清晰地觀察三維結(jié)構(gòu)體，本文對(duì)上述三維結(jié)構(gòu)圖進(jìn)行切片顯示，如圖4(d)和圖4(e)所示。

　　在圖中很難直觀地看出兩重建結(jié)果與訓(xùn)練圖像之間的差異，但借助形態(tài)特性函數(shù)可以很好地分析圖像的結(jié)構(gòu)。

　　本文使用另一個(gè)圖像形態(tài)特性描述函數(shù)——兩點(diǎn)簇函數(shù)［10］來(lái)評(píng)價(jià)重建結(jié)果圖像。兩點(diǎn)簇函數(shù)C2(r)定義了圖像上兩個(gè)相距歐氏距離為r的相位點(diǎn)落在同一相位簇上的概率，表征了圖像一定的連接特性。

　　限定r的范圍在(10,40)內(nèi)，對(duì)比兩組重建結(jié)果與原圖兩點(diǎn)簇函數(shù)概率分布，如圖5所示。從圖中可以看出使用本文算法的重建結(jié)果較單約束條件下穩(wěn)相法重建結(jié)果更貼近原始圖像，這也證實(shí)了加入隨機(jī)線性路徑函數(shù)后重建的結(jié)果保存了更多原圖的連接特性，使得重建結(jié)果更加精確。

　　3.3誤差分析

　　為了更加直觀地觀察兩種重建結(jié)果在兩點(diǎn)簇函數(shù)上與原圖的差異，本文采用相對(duì)誤差公式［11］分析重建結(jié)果在兩點(diǎn)簇函數(shù)上存在的誤差，如式（7）所示。

　　其中，ε-為重建結(jié)果與原圖的相對(duì)誤差，r代表點(diǎn)對(duì)的距離，C(r)與C-(r)分別代表了原圖和重建結(jié)果的兩點(diǎn)簇函數(shù)。

　　針對(duì)兩種重建方法，本文對(duì)10組尺寸為256×256的不同訓(xùn)練圖像進(jìn)行對(duì)比分析，得出重建結(jié)果在兩點(diǎn)簇函數(shù)上相對(duì)誤差的平均結(jié)果，相關(guān)實(shí)驗(yàn)結(jié)果如表1所示。 

　　從表1中的實(shí)驗(yàn)結(jié)果可以看出，兩種算法重建結(jié)果隨著距離增大，相對(duì)誤差總體呈上升趨勢(shì)，本文算法因采用了雙約束條件，相較于單約束條件重建結(jié)果平均相對(duì)誤差更小。相對(duì)誤差的產(chǎn)生是因?yàn)殡S機(jī)線性路徑函數(shù)與兩點(diǎn)相關(guān)函數(shù)都不能完全表征原始圖像連接特性，在距離增大時(shí)，誤差產(chǎn)生的概率越高。

4結(jié)論

　　本文主要討論了隨機(jī)方向上的線性路徑函數(shù)，給出了相應(yīng)計(jì)算方法。進(jìn)而介紹了如何將隨機(jī)線性路徑函數(shù)加入到穩(wěn)相分級(jí)模擬退火法中，并將改進(jìn)后與改進(jìn)前的重建結(jié)果與原圖在兩點(diǎn)概率函數(shù)、線性路徑函數(shù)以及兩點(diǎn)簇函數(shù)上做相應(yīng)的對(duì)比，論證了隨機(jī)線性路徑函數(shù)在重建方法中的真實(shí)性與有效性。

參考文獻(xiàn)

　?。?］ LATIEF F D E, BISWAL B, FAUZI U, et al. Continuum reconstruction of the pore scale microstructure for Fontainebleau sandstone［J］. Physica A, 2010, 389：16071618.

　?。?］ ALRAOUSH R. Microstructure characterization of granular materials［J］. Physica A Statistical Mechanics & its Application, 2007, 377(2)：545558.

　?。?］ JONES K W, FENG H, TOMOV S.Characterization of methane hydrate host sediments using synchrotroncomputed microtomography (CMT)［J］. Journal of Petroleum Science & Engineeriug, 2007, 56(13)：136145.

　　［4］ RINTOUL M D, TORQUATO S. Reconstruction of the structure of dispersions［J］. Journal of Colloid & Interface Science, 1997, 186(2)：186,467476.

　?。?］ ALEXANDER S K, FIEGUTH P, IOANNIDIS M A, et al. Hierarchical annealing for synthesis of binary images［J］. Mathematical Geosciences, 2009, 41(4)：357378.

　?。?］ Chen Dongdong, Teng Qizhi, He Xiaohai, et al. Stablephase method for hierarchical annealing in the reconstruction of porous media images［J］. Physical Review E, 2014, 89：013305.

　　［7］ Chen Dongdong, He Xiaohai, Teng Qizhi, et al. Reconstruction of multiphase microstructure based on statistical descriptors［J］. Physica A, 2014, 415：240250.

　?。?］ DEROSSI A, SEVERINI C, De PILLI T. Reconstruction of food microstructure via statistical correlation functions. The use of linealpath distribution functions［J］. Food Engineering, 2014, 142(6)：916.

　?。?］ DEROSSI A, De PILLI T, SEVERINI C. Use of linealpath distribution function as statistical descriptor of the crumb structure of bread［J］. Food Biophysics, 2013, 8(4)：223232.

　?。?0］ TORQUATO S, BEASLEY J D, CHIEW Y C. Twopoint cluster function for continuum percolation［J］. Chemical Physics, 1988, 88：65406546.

　?。?1］ 袁浩，滕奇志，卿粼波，等．三維模型兩點(diǎn)連通概率函數(shù)的快速計(jì)算.計(jì)算機(jī)工程與應(yīng)用［J］．計(jì)算機(jī)工程與應(yīng)用， 2015，51（3）：181185.