秦杰1，吳運強2，滕奇志1，王子強2，何小海1

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

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

　　關鍵詞：分級模擬退火；兩點相關函數；線性路徑函數；三維重建

0引言

　　巖心的三維結構分析在石油地質巖心研究中有很重要的意義，通過對巖心的三維結構分析可以得出其微觀結構特性進而可以計算出其宏觀特性參數［12］。但是在多數情況下，獲取目標真實三維結構難度較大［3］，其單張的二維結構圖像卻可以方便得到，所以在二維圖像的基礎上進行三維重建有非常大的研究價值。

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

　　本文針對該問題，引入隨機線性路徑函數，與兩點相關函數同時約束重建過程，利用分級策略提高時間效率。對巖心圖像進行對比實驗，結果表明相較于單約束條件，改進算法顯著提高了重建精度，有效完成了巖心圖像的三維重建。

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

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

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

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

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

2隨機線性路徑函數

　　2.1線性路徑函數概述

　　不同于兩點概率函數，線性路徑函數［79］描繪的是將一些隨機的標量線段拋入參考圖像的像素點中，各個標量線段落在特定相位上的概率。在兩點概率函數中，只是給出了線段首位兩點相位概率分布，而線性路徑函數彌補了兩點概率的不足，給出了整條線段的相位概率分布。線性路徑函數的數學表達式如下：

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

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

　　顯然，線性路徑函數表征了圖像目標相的連通特征。與兩點概率函數不同，不同相位的線性路徑函數之間是相互獨立的，并沒有任何線性關系。

　　2.2隨機線性路徑函數

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

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

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

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

　　結合式（2），將兩種限制函數能量相加，可得到退火時的總能量：

　　E=Etp+Elp(6)

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

3重建結果分析

　　將本文算法應用于儲集層巖心圖像重建中，進行對比實驗，驗證本文算法的真實性與有效性。在處理過程中，白色相作為1相位，黑色相作為2相位，灰色相記為3相位。

　　3.1本文重建結果與原始圖像對比分析

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

　　3.2三維重建結果對比分析

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

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

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

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

　　3.3誤差分析

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

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

　　針對兩種重建方法，本文對10組尺寸為256×256的不同訓練圖像進行對比分析，得出重建結果在兩點簇函數上相對誤差的平均結果，相關實驗結果如表1所示。 

　　從表1中的實驗結果可以看出，兩種算法重建結果隨著距離增大，相對誤差總體呈上升趨勢，本文算法因采用了雙約束條件，相較于單約束條件重建結果平均相對誤差更小。相對誤差的產生是因為隨機線性路徑函數與兩點相關函數都不能完全表征原始圖像連接特性，在距離增大時，誤差產生的概率越高。

4結論

　　本文主要討論了隨機方向上的線性路徑函數，給出了相應計算方法。進而介紹了如何將隨機線性路徑函數加入到穩(wě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.

　　［2］ 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.

　　［6］ 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.

　?。?］ 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］ 袁浩，滕奇志，卿粼波，等．三維模型兩點連通概率函數的快速計算.計算機工程與應用［J］．計算機工程與應用， 2015，51（3）：181185.