基于粒子群優化的巖石薄片三維圖像重建

滕奇志 唐 棠 李征驥 何小海

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

1 引言

為了分析巖石的微觀結構參數，常用顯微成像方式獲得巖石薄片的 2維圖像并對其進行定量分析。隨著分析技術的發展，對巖石3維微觀結構的研究越來越引起重視。由于巖石的硬度和致密性，均勻連續切片基本不可行[1]。CT(Computed Tomography)圖像重建，能根據接收的信號重建樣本的2維剖面圖像，獲得其3維序列圖像，但樣本尺寸或分辨率受到局限，且設備比較昂貴，尚不能在實驗室普遍使用[2]。因此人們致力于采用數學方法構建3維模型。

相對于砂巖的斷層掃描圖像序列，其高分辨率2維顯微圖像的獲取十分便利。因此考慮通過提取2維圖像的統計特征，在3維空間構建與2維圖像具有相同特征、與真實3維結構具有相似特性的3維圖像。同一幅2維圖像的重建結果可能在視覺上有所差異，但它們都具有相同的統計特征，這是本文的重建算法與常規CT圖像重建的顯著區別。重建結果不僅便于研究人員直觀認識砂巖的微觀結構，更具有2維圖像缺乏的深度方向的連通性信息，為砂巖滲流特性的研究提供平臺。

國外有學者提出了一些基于2維圖像的3維重建方法，較典型的有模擬退火算法[3,4]和基于高斯隨機場的重建方法[5]。模擬退火(SA)算法可以引入多個統計特性，重建結果能準確反映2維圖像的統計特征，但計算量大，耗時長。針對重建效率問題國內外學者開展了一系列工作[6]，本文作者也作了較深入的研究[7,8]，但在實際工作中重建大規模的3維結構仍然比較困難。而高斯隨機場變換法利用Wiener-Khinchin定理，能在較短時間內重建與 2維圖像具有相同低階統計特性的3維結構。該方法最早通過對隨機高斯場進行線性變換實現，但在重建過程中，很難用數值模擬的方法確定線型變換過程中的系數，使得該方法在應用上受到了很大的限制。Liang等人[5]提出利用傅里葉變換(FFT)重建3維結構，不需要進行線性變換，但由于使用經驗公式估計重建參數，重建結果與2維圖像之間統計特征誤差較大，文獻[9]在其基礎上作了進一步工作，但仍未解決這一問題。為此，本文提出一種新的重建方法，以3維結構與2維圖像的自相關分布誤差最小為目標，用粒子群優化(PSO)算法確定重建參數，提高重建結果的準確性。與模擬退火算法相比，達到相同的重建精度，需要的時間大大減少。

2 重建原理

重建的目的是構建與2維圖像具有相同一階和二階統計特性的3維結構，該結構具有與真實結構表現相似的物理特性(主要指滲流特性)。重建的流程如圖1所示，下面分別介紹其中的關鍵環節。

圖1 重建流程

2.1 一階和二階統計

2維圖像和3維結構由孔隙和顆粒構成，稱為兩種“相”，其相函數I(u)表示為[4]

其中u是表示3維結構中任意一點位置的矢量。孔隙度ε是I(u)的一階統計，表示巖石中孔隙所占體積的比例[4]。

自相關函數C(r)是I(u)的二階統計，表示2維圖像或3維結構中距離為r的兩點同時屬于孔隙相的概率，也稱為兩點相關函數[2]。

將C(r)歸一化為RI(r)，使

本文假設2維圖像和3維結構都是各向同性的均勻介質，其自相關函數只與r的大小有關，與方向無關。因此 2維圖像的自相關函數可用RI(r)表示。自相關函數具有對稱性，對大小為N×N的2維圖像，只需計算0≤r≤N/2時的RI(r)。

2.2 自相關函數的擴展

重建的3維結構與2維圖像具有相同的自相關函數，可將 2維圖像的自相關函數RI(r)擴展為 3維自相關函數RY(i,j,k)[5]。

其中0≤i≤Nx/2, 0≤j≤Ny/2, 0≤k≤Nz/2,Nx,Ny,Nz為分別3維結構的長、寬、高，3維自相關函數的其余部分通過自相關函數的對稱性獲得。

2.3 通過傅里葉變換重建3維結構

Wiener-Khinchin定理指出，隨機過程f(t)的自相關函數R(r)與其功率譜|F(u)|2互為傅里葉變換對[5]。

如果已知隨機過程的自相關函數，可利用傅里葉變換得到具有相同自相關函數的隨機過程。已知3維結構的自相關函數RY(i,j,k)，可以通過傅里葉變換，計算 3維結構的功率譜|F(u,v,w)|2，進而通過設置適當的輻角 θ(u,v,w)，構建該3維結構的傅里葉變換[5]。

最后由傅里葉反變換得到自相關函數為RY(i,j,k)的3維結構Y(u)。輻角 θ(u,v,w)在[0,2π)間隨機取值，滿足一定的對稱規則，保證F(u,v,w)的傅里葉反變換Y(u)屬于實數域。

2.4 離散化

根據Wiener-Khinchin定理，由傅里葉反變換得的3維結構Y(u)，是服從正態分布的連續實數，其分布函數如式(8)所示。為了使重建結果的孔隙度與2維圖像一致，并采用與相函數一樣的方式表示重建結果中的孔隙和顆粒，Adler等人[10]采用式(8)得到離散Z(u)。

2.5 自相關函數的校正

式(9)中的離散化是非線性變換，Adler等人[10]給出了Y(u)和Z(u)的自相關函數RY(r)和RZ(r)之間的關系如式(10)，根據式(10)及Cm的定義可推導出式(11)及式(12)。

或

由式(11)和式(12)可知RZ(r)與RY(r)的關系僅與孔隙度ε有關，文獻[10]中給出了不同孔隙度下RZ(r)與RY(r)的關系圖示。

重建的目的是使離散結果Z(u)的自相關函數RZ(r)與2維圖像的自相關函數RI(r)相同，因此重建時不能直接擴展2維圖像的自相關函數RI(r)，而應根據式(10)中關系得到RY(r)。再按式(5)將RY(r)擴展為Y(u)的自相關函數RY(i,j,k)。

2.6 存在問題

基于傅里葉變換的重建算法解決了基于高斯隨機場重建方法中的瓶頸問題。但式(11)及式(12)中給出的是RY(r)與RZ(r)關系的近似，采用該方法估計RY(r)會造成重建3維結構與2維圖像兩點相關函數之間存在較大的誤差。導致重建結果的多方面參數特性與2維圖像吻合較差，無法真實地重現出3維物體的結構空間。

3 基于粒子群優化的3維重建

針對上面提到的傅里葉變換算法的不足，本文作者提出應用粒子群優化算法確定適當的RY(r)，使重建3維結構Z(u)與2維圖像的自相關分布RI(r)更接近。粒子群優化(Particle Swarm Optimization,PSO)算法是一種迭代尋優算法，由于其概念簡單，易于實現，近年來被廣泛用于解決各種優化問題[11,12]。

將粒子群優化算法引入高斯隨機場重建法的目的是找到連續高斯場Y(u)的適當的自相關分布RY(r)，使重建3維結構與2維圖像的自相關分布更接近。

3.1 粒子編碼

群體中粒子個數為n，每個粒子的位置對應于一組RY(r)，粒子搜索空間的維數D與RY(r)中r的范圍相同。對粒子i，其位置ax∈RD，速度v∈iiRD,i=1,2 ,…,n。ax(j)和v(j)是粒子i第j維的ii位置和速度，j=1,2,…,D。是粒子i搜索發現的最優位置。是整個種群中所有粒子發現的最優位置。按照兩點相關函數非負的限制，及歸一化后兩點相關函數的范圍，設置axmin=0,axmax=1，速度限制設為vmin=-1,vmax=1。

3.2 粒子群初始化

一般采用隨機初始化方法，使每個粒子的位置和速度在[axmin,axmax]和[vmin,vmax]之間均勻分布。本文作者在研究中發現，與隨機產生的粒子位置相比，用式(10)估計的粒子位置具有相對較低的適應值，因此以式(10)估計的結果作為初始種群中的一個粒子的位置，其余粒子的初始位置隨機產生。

3.3 適應度計算

根據兩點相關函數僅與距離r相關，與方向無關的性質，將axi映射到 3維空間，得到自相關函數Ri(x,y,z)[5]。

其中N為3維結構的長、寬、高。根據兩點相關函數的性質，3維自相關函數的其余部分通過自相關函數的對稱性獲得。根據Wiener-Khinchin定理，可根據2.3節和2.4節中的方法重建3維結構。

為了度量按照粒子i的位置axi重建的3維結構Zi(u)與2維圖像的統計特征的不同，適應度函數定義為Zi(u)的自相關函數Ri(r)與2維圖像兩點相關函數RI(r)之間的誤差為

適應度越小說明重建結果與2維圖像的自相關函數分布越相似。

3.4 更新粒子狀態

其中c1和c2分別為局部和全局加速因子，用于調整和的影響強度。r1和r2∈RD，是[0,1]之間均勻分布的隨機數。ω是慣性權重，ω值較大時，全局尋優能力強，局部尋優能力弱，ω較小則反之。如果f(axi) ＜f()，則更新局部最優位置=axi。如果f(axi) ＜f()，則更新全局最優位置=axi。

根據大量的實驗發現，動態地調整ω能獲得比固定值更好的尋優結果，本文中令ω按照凹函數遞減：

其中maxω和minω分別為ω的最大值和最小值，ML和CL分別為最大迭代次數和當前迭代次數。

4 重建結果及分析

本文實驗所用硬件平臺配置為 CPU: P8600(2.4 G)，內存：DDR3-10664 G(2 G×2)。

本文的研究內容是針對只能獲取2維圖像時進行的3維重建，但為了說明重建結果，本文采用高分辨率CT掃描序列圖進行實驗，即隨機選取CT序列圖中的一幅圖像作為2維參考圖像進行3維重建。如圖2中，(a)為CT掃描序列圖，(b)是該樣本中一層切片的圖像，圖像尺寸為128×128, 2維面孔率為17.6697%，其中白色部分為孔隙。

圖2 巖心CT掃描圖

同時，將CT序列圖直接進行3維重構，并計算其孔隙度、兩點相關函數，與3種重建方法的結果進行對照。圖3為直接重構的真實3維圖像，N=128，孔隙度為19.2968%。圖4為真實3維結構圖像和選取的2維圖像的兩點相關函數對照，從圖中可以看出，隨機選取的2維圖像與真實3維結構在兩點相關函數特征上有一定的差異。

圖3 巖心真實3維結構圖

圖4 3維真實結構與參考圖像的兩點相關函數

下面分別用基于粒子群優化的重建算法、基于傅里葉變換的重建算法以及模擬退火算法重建大小為128×128×128的3維結構。然后分別從孔隙度、兩點相關函數以及計算時間上進行對比。

4.1 粒子群優化重建與傅里葉變換重建

粒子群優化重建算法的參數設置為：粒子數量n=10, ML=10,c1=c2=2, maxω=1, minω=0.5，重建時間為27 min，適應度為 5.9 × 10?2。重建結構的孔隙度為17.6697%，與2維參考圖像吻合度很好。

傅里葉變換重建采用文獻[5]的算法進行，重建時間為 35 s，適應度為 3.94 × 10?1，重建結構的孔隙度為17.4662%，與2維參考圖像相差0.2035%。兩種算法的重建結果見圖5。

圖5 3種重建算法效果對照

4.2 模擬退火重建

模擬退火重建方法是一種經典的重建算法，根據2維參考圖像，直接以孔隙的3維結構作為重建問題的解。為了與基于粒子群優化的重建算法比較，本文中模擬退火重建方法也以2維圖像的自相關函數作為重建的約束條件。重建的起始溫度設為使第1次交換接受的概率為0.5，重建結束的條件設為溫度低于 1.0 × 10?37。大量實驗結果表明，重建結果與溫度下降速度有關，每迭代5000次溫度下降1%的重建結果，其適應度與粒子群重建算法的適應度基本一致，重建結構的孔隙度為17.6697%，與2維參考圖像吻合很好，后面也就基于此結果進行對比。

4.3 重建結果比較

圖5為3種重建算法的直觀效果對照，將此圖與圖3的真實3維結構相比較，可以看出，本文提出的基于粒子群優化算法的結果與真實3維結構相似度最好。

圖6為3種重建算法與2維參考圖像的兩點相關函數進行對比的結果。

圖6 2維參考圖像與3種重建算法兩點相關函數對比

由實驗數據看出，對于兩點相關函數特征，模擬退火算法最逼近2維參考圖像，粒子群算法次之，傅里葉變換算法最差。但由于2維參考圖像的兩點相關函數與3維真實結構的兩點相關函數本身就存在差異，因此，再用真實3維結構的兩點相關函數與3種重建算法進行比較，從圖7可以看出，本文提出的基于粒子群優化算法與真實3維結構的兩點相關函數最逼近。

圖7 真實3維結構與重建算法結果的兩點相關函數對照

4.4 實驗結果分析

根據上述實驗，從直觀視覺效果看，基于粒子群優化的重建方法最接近真實樣本的3維結構。

從實驗數據看，在孔隙度、兩點相關函數及重建時間幾個方面的對比如下：

(1)孔隙度：基于粒子群優化的重建算法和模擬退火算法所得3維結構的孔隙度與原始2維圖像完全一致，基于傅里葉變換的重建算法結果的孔隙度與原始圖片有 0.2035%的誤差，3種重建結果與真實3維圖像的孔隙度均有差異，這是由于重建算法都是以 2維參考圖像的孔隙度作為約束條件造成的。

(2)兩點相關函數：對于2維參考圖像，模擬退火算法的重建結果與2維參考圖像基本一致，粒子群優化算法的重建結構次之，而基于傅里葉變換的重建結果與2維參考圖像相差最大。而對于真實3維結構，粒子群優化算法的兩點相關函數與真實 3維結構最接近，好于其他兩種算法。

(3)重建時間： 基于傅里葉變換的算法耗時約35 s，基于粒子群優化的重建算法耗時約為27 min。模擬退火算法要達到與粒子群算法同樣的重建效果，需要的重建時間為74 min。粒子群算法比傅里葉變換算法耗時多，但換來了明顯更好的重建效果。

5 結束語

本文作者提出了一種新的基于2維圖像的砂巖3維結構重建算法，即利用粒子群優化算法確定 3維結構的自相關函數，再通過傅里葉變換進行重建。為了驗證算法的有效性，分別從孔隙度、兩點相關函數以及重建時間上對實驗結果進行了分析對比。從實驗結果可知，該方法解決了傅里葉變換重建算法使用經驗公式估計重建參數，重建結果與2維參考圖像之間統計特征誤差較大的問題，得到的重建結果在孔隙度和兩點相關函數方面，均與2維參考圖像及真實3維結構更接近。與經典的模擬退火重建算法相比，在重建效果相當的情況下，基于粒子群優化的重建方法重建速度提高很多，具有更高的效率。

本文提出的粒子群優化算法只采用了兩點相關函數作為特征函數進行重建，特征約束比較單一。另外，從實驗中可看出，選取的2維參考圖像在3維重建中起著重要的作用，因此，如何在重建中加入更多的特征條件以及如何從巖石樣品中選取最能代表其結構的2維切片圖像，是需要進一步研究的問題。

[1] ?ren P E and Bakke S. Process based reconstruction of sandstones and prediction of transport properties [J].Transport in Porous Media, 2002, 46(2): 311-343.

[2] Politis M G, Kikkinides E S, Kainourgiakis M E,et al.. A hybrid process-based and stochastic reconstruction method of porous media [J].Microporous and Mesoporous Materials,2008, 110(1): 92-99.

[3] Yeong C L Y and Torquato S. Reconstructing random media[J].Physical Review E, 1998, 57(1): 495-506.

[4] Yeong C L Y and Torquato S. Reconstructing random media.II. Three-dimensional media from two-dimensional cuts[J].Physical Review E, 1998, 58(1): 224-233.

[5] Liang Z R, Fernandes C P, Magnani F S,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.

[6] Liu Xue-feng, Sun Jian-meng, and Wang Hai-tao.Reconstruction of 3-D digital cores using a hybrid method[J].Applied Geophysics, 2009, 6(2): 105-112.

[7] 滕奇志, 唐棠, 何小海, 等. 基于交換-單親遺傳算法的砂巖三維顯微圖像重建[J]. 數據采集與處理, 2010, 25(3):364-368.Teng Q, Tang T, He X,et al.. Reconstruction of 3D microstructure of sandstone based on swap-aprtheno-genetic algorithm[J].Journal of Data Acquisition and Processing,2010, 25(3): 364-368 .

[8] Tang T, Teng Q, and He X,et al.. A pixel selection rule based on the number of different phase neighbors for the simulated annealing reconstruction of sandstone microstructure [J].Journal of Microscopy-Oxford, 2009,234(3): 262-268.

[9] 趙凱, 李強, 宣益民. 基于圖象處理和傅里葉變換的三維多孔介質重構方法[J]. 工程熱物理學報, 2008, 29(2): 287-290.Zhao K, Li Q, and Xuan Y. A reconstruction technique for three dimensional porous media using image analysis and Fourier transform [J].Journal of Engineering Thermophysics,2008, 29(2): 287-290.

[10] Adler P M, Jacquin C G, and Quiblier J A. Flow in simulated porous media [J].International Journal of Multiphase Flow, 1990, 16(4): 691-712.

[11] 陶新民, 徐晶, 楊立新, 等. 一種改進的粒子群和K均值混合聚類算法[J]. 電子與信息學報, 2010, 32(1): 92-97.Tao Xin-min, Xu Jing, Yang Li-biaoet al.. Improved cluster algorithm based on K-means and particle swarm optimization [J].Journal of Electronics&Information Technology, 2010, 32(1): 92-97.

[12] 陳榮元, 張飛艷, 張斌, 等. 基于數據同化和粒子群優化算法的遙感影像融合[J]. 電子與信息學報, 2009, 31(10):2509-2513.Chen Rong-yuan, Zhang Fei-yan, Zhang Bin,et al.. Remote sensing image fusion based on data assimilation and particle swarm optimization [J].Journal of Electronics&Information Technology, 2009, 31(10): 2509-2513.