莫修文, 張強， 陸敬安

1 吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院, 長春　130026 2 廣州海洋地質(zhì)調(diào)查局, 廣州　510760

模擬退火法建立數(shù)字巖心的一種補充優(yōu)化方案

莫修文1, 張強1， 陸敬安2

1 吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院, 長春130026 2 廣州海洋地質(zhì)調(diào)查局, 廣州510760

摘要針對傳統(tǒng)的模擬退火建立數(shù)字巖心方法建模效率低、計算結(jié)果存在較多的孤立點等缺陷，設(shè)計了后續(xù)的補充搜索方案，定義了孔隙與骨架的邊界點選取準(zhǔn)則，提出了基于擇多算子算法選取的對象點與孔隙骨架邊界點交換作為新系統(tǒng)的生成方案.使用傳統(tǒng)模擬退火算法對隨機的初始模型優(yōu)化作為第一階段，用新系統(tǒng)生成方案繼續(xù)優(yōu)化作為第二階段，最后減少模型中孤立孔隙點或骨架點，輸出結(jié)果.試驗表明，新的補充優(yōu)化方案不僅可行而且效果良好.最終優(yōu)化模型與傳統(tǒng)模擬退火建模結(jié)果相比，無論是在孔隙形態(tài)還是在孔隙連通性上都與原始巖心較為接近，由反映巖石孔隙空間特征的變差函數(shù)可以看出補充搜索方案優(yōu)化結(jié)果優(yōu)于傳統(tǒng)模擬退火法建模結(jié)果.關(guān)鍵詞模擬退火； 數(shù)字巖心； 自相關(guān)函數(shù)； 線性路徑函數(shù)； 變差函數(shù)

1引言

目前，盡管巖石物理實驗可以描述多孔介質(zhì)中的微觀機理，但隨著研究的不斷深入，巖石物理實驗方法已不能完全滿足一些特殊儲層的研究需要，如天然氣水合物儲層巖心研究、裂縫發(fā)育的碳酸鹽巖儲層的巖心獲取問題、深部鉆探巖心物性分析等(姚軍和趙秀才,2010).隨著計算機應(yīng)用技術(shù)的發(fā)展，數(shù)值模擬方法逐漸成為研究巖石物理性質(zhì)的重要輔助手段，而數(shù)字巖心技術(shù)是數(shù)值模擬的前提，其建模的好壞直接影響巖石物理數(shù)值計算的結(jié)果.

數(shù)字巖心建模方法主要有兩種，第一種是物理實驗方法構(gòu)建數(shù)字巖心，即采用物理實驗方法對巖心進行掃描，用實際巖心掃描數(shù)據(jù)直接建立三維數(shù)字巖心(Lymberopoulos and Payatakes,1992).第二種是數(shù)值模擬方法如高斯場法(Joshi,1974; Quiblier,1984)、模擬退火法(Hazlett,1997)、過程法(Bakke and Oren,1997)、順序指示法(Keehm,2003)、多點地質(zhì)統(tǒng)計法(Okabe and Blunt,2005)、馬爾可夫鏈-蒙特卡洛法(Wu et al.,2004)等.其中，模擬退火算法是20世紀(jì)80年代初期發(fā)展起來的一種求解大規(guī)模組合優(yōu)化問題的隨機性方法，與高斯模擬法相比，模擬退火法的優(yōu)勢在于它能夠?qū)⒏嗟膸r石信息考慮進來，使所建模型與真實巖石在孔隙空間結(jié)構(gòu)上更接近(姚軍等,2007).過程法建模結(jié)果在孔隙連通性上優(yōu)于模擬退火法(Bakke and Oren,2003)，但過程法建模過程較為復(fù)雜，更適合于成巖過程簡單的巖石(聶昕,2014).多點地質(zhì)統(tǒng)計法所建的數(shù)字巖心孔隙連通性好，但計算速度較慢且只適合于各向同性介質(zhì)(張挺等,2009).可見，每種建模方法都有自身的優(yōu)點與局限，基于優(yōu)缺互補的原則，近期陸續(xù)發(fā)展了模擬退火法與其他建模方法相結(jié)合的混合數(shù)字巖心建模法，如將高斯模擬法的輸出作為模擬退火法的輸入，使得建模時間大大減少(Hidajat et al.,2002)；針對過程法的缺陷，劉學(xué)鋒等將過程法和模擬退火法相結(jié)合，提出了三維數(shù)字巖心重構(gòu)的混合法(劉學(xué)鋒,2010)等.

模擬退火法作為較早建立數(shù)字巖心的一類方法，它能包含反映巖石空間結(jié)構(gòu)特征的較多信息，與其他方法相比有不可替代的優(yōu)勢，但目前應(yīng)用中的問題主要有建模效率低下、建模結(jié)果孔隙連通性差、存在過多的孤立點等(趙秀才等,2007).本文主要針對傳統(tǒng)模擬退火法運行到后期出現(xiàn)的問題，設(shè)計了基于模擬退火算法的新系統(tǒng)生成方案，作為后續(xù)補充優(yōu)化算法，并取得很好效果.對于其他隨機建模方法，如多點地質(zhì)統(tǒng)計法、馬爾可夫鏈蒙特卡羅法等，均可用此后續(xù)補充優(yōu)化方法來做進一步的細致優(yōu)化,有較好的實用性和移植性.

2數(shù)字巖心的建模函數(shù)

2.1初始隨機模型計算函數(shù)

假設(shè)孔隙度為φ的巖石中只考慮孔隙與巖石骨架兩相系統(tǒng)，并用1來表示孔隙相，0代表巖石骨架.那么產(chǎn)生一個0～1之間的隨機數(shù)r，若r在(0,φ)之間，那么此點Z(r)便為孔隙點，否則為骨架點.

(1)

式中，Z(r)為模擬點的狀態(tài)；r為(0,1)之間的隨機數(shù)；φ為巖石孔隙度.

對模擬三維巖心x,y,z三個方向依次掃描賦值便可得到完全隨機的初始巖心模型.

2.2兩點概率函數(shù)

多相系統(tǒng)中第j相的兩點概率函數(shù)定義如下：

(2)

式中，r1,r2為系統(tǒng)中的相距一定距離r的任意兩點.

(3)

并具有以下性質(zhì)：

(4)

對只考慮孔隙與骨架的兩相系統(tǒng)而言，建模過程中通常以孔隙相為研究對象，故以S(r)來簡化表示.對于各項異性系統(tǒng)，為了盡可能多地獲取各個方向的介質(zhì)信息，可選取不同的方向來計算兩點概率函數(shù)值，但計算量也會增加.對于各項同性介質(zhì)，為了節(jié)省運算時間只需沿一個主軸方向計算即可.

2.3線性路徑函數(shù)

線性路徑函數(shù)L(i)(r1,r2)是描述巖石微觀孔隙結(jié)構(gòu)的重要函數(shù)，其定義如下：

(5)

式中，r1，r2為系統(tǒng)中相距為r的任意兩點；rx為連接r1，r2線段上的任意點.

由以上定義可以看出，線性路徑函數(shù)體現(xiàn)了系統(tǒng)的連通性信息，在各向同性介質(zhì)中，線性路徑函數(shù)的值取決于r1，r2這兩點之間的距離r，故表達式可以簡化為L(i)(r)，對于孔隙相為φ的系統(tǒng)有：L(j)(0)=φ.

3傳統(tǒng)的模擬退火建模方法

3.1模擬退火法的基本理論

模擬退火算法(SA)是模擬物理退火過程的一種最優(yōu)化算法(Metropolis et al.,1953).據(jù)統(tǒng)計熱力學(xué)研究表明,在某一特定溫度T下,物體原子能量的分布概率P滿足如下關(guān)系：

(6)

式中，E(α)為狀態(tài)α的能量；KB為Boltzman常數(shù)；Z(T)為標(biāo)準(zhǔn)化因子.

當(dāng)系統(tǒng)溫度經(jīng)過退火而緩慢降低后，將達到能量最低的平衡狀態(tài).SA算法就是模擬系統(tǒng)在退火過程中原子能量的概率分布來進行優(yōu)化計算的.定義在第(K+1) 次搜索時原子能量概率分布的接受準(zhǔn)則(Metropolis準(zhǔn)則)為

(7)

式中，T為控制參數(shù)溫度.

相應(yīng)狀態(tài)的接受準(zhǔn)則為

(8)

式中，RAN(0,1)為在(0,1)內(nèi)的隨機數(shù).

滿足式(8)的狀態(tài)是第(k+1)次狀態(tài).這樣搜索既能向性能指標(biāo)優(yōu)化的方向迭代，又有一定概率接受性能指標(biāo)劣化的狀態(tài).開始時溫度較高，接受劣態(tài)的概率大，隨著迭代次數(shù)增加，系統(tǒng)得到改善，逐漸達到全局最優(yōu)解.

3.2模擬退火法數(shù)字巖心建模的實現(xiàn)

本文所建數(shù)字巖心模型是僅考慮骨架和孔隙的兩相系統(tǒng)，因流體的儲滲空間為孔隙，故以孔隙空間為研究對象計算各特征函數(shù)值.建模目標(biāo)函數(shù)設(shè)置為重構(gòu)系統(tǒng)與參考系統(tǒng)的自相關(guān)函數(shù)與線性路徑函數(shù)差值的平方和，某時刻系統(tǒng)的目標(biāo)函數(shù)可定義為

(9)

式中，αi,βi為不同約束函數(shù)的權(quán)重值；E為系統(tǒng)的總能量.

數(shù)字巖心建模步驟如下： (1) 對實際巖心CT圖像或電鏡掃描圖像進行二值化處理，統(tǒng)計建模基本資料φ、S(r)、L(r); (2) 隨機生成孔隙度為φ的三維數(shù)字巖心(0代表巖石骨架，1代表孔隙)，并計算初始模型目標(biāo)函數(shù)E，進入優(yōu)化進程; (3) 在孔隙和巖石骨架空間任意選取一孔隙點和骨架點，兩者互換位置生成新系統(tǒng)，并計算目標(biāo)函數(shù)值E′; (4) 新系統(tǒng)由Metropolis準(zhǔn)則來判定是否取代舊系統(tǒng); (5) 根據(jù)內(nèi)循環(huán)終止條件和外循環(huán)終止條件來決定是否終止循環(huán); (6) 最后輸出數(shù)字巖心模擬結(jié)果，結(jié)束進程.

4基于模擬退火法的補充優(yōu)化方案

上述傳統(tǒng)模擬退火法建立數(shù)字巖心，開始時對完全隨機的模型進行優(yōu)化，溫度最高，由Metropolis準(zhǔn)則可以看出，劣態(tài)系統(tǒng)的接受概率大，這將有效避免系統(tǒng)陷入局部最優(yōu)，導(dǎo)致循環(huán)終止.隨著溫度的降低，迭代次數(shù)的增加，系統(tǒng)孔隙結(jié)構(gòu)的形態(tài)逐步形成，但分布較為分散，連通性差，孤立的孔隙點或骨架點較多.如果再按傳統(tǒng)的隨機交換骨架點和孔隙點作為新系統(tǒng)的生成方案，就會導(dǎo)致原有已形成的孔隙結(jié)構(gòu)遭到破壞，孔隙結(jié)構(gòu)陷入形成-破壞-再形成的惡性循環(huán)，無形中增加了不必要的計算量.而孔隙中的骨架點和骨架中孤立的孔隙點都是不可能大量存在的，新系統(tǒng)生成方案選取的對象點主要以孤立孔隙點和骨架點為目標(biāo).故本文設(shè)計了兩個階段的優(yōu)化方案，第一階段用傳統(tǒng)的SA方法對初始的隨機模型做優(yōu)化.第二階段以擇多算子選取的對象點與孔隙骨架的邊界點交換生成新系統(tǒng)，作為補充優(yōu)化階段.

圖1　D3Q19模型網(wǎng)格結(jié)構(gòu)Fig.1　The D3Q19 lattice structure model

新系統(tǒng)生成方案的對象點選取標(biāo)準(zhǔn)是以D3Q19模型網(wǎng)格結(jié)構(gòu)(Michalewict et al.,1996)為基礎(chǔ)，借助二值模式下的19鄰域來定義擇多算子(劉學(xué)鋒等,2013；趙秀才等,2008)，見圖1，網(wǎng)絡(luò)結(jié)構(gòu)的中心體素0，若在其19鄰域中體素的狀態(tài)值與中心體素相反的個數(shù)大于等于n，則將該中心體素狀態(tài)判定為新系統(tǒng)生成方案選取的對象點.n的取值根據(jù)具體情況而定，若n取值過高，會導(dǎo)致符合條件的點越來越難選取，從而使程序在尋找孤立點上花費過多時間.若n的取值過低，則相當(dāng)于傳統(tǒng)的系統(tǒng)生成方案，導(dǎo)致優(yōu)化后模型與實際巖心孔隙形態(tài)相差較大.

孔隙與骨架邊界點的判定：若中心體素點的狀態(tài)為孔隙，且周圍18個鄰域點的狀態(tài)均為孔隙，則中心體素點為完全孔隙點，空間位置用(i,j,k)表示.如圖2所示，用0～5表示6個不同的方向，0代表x軸正方向，3代表x軸負方向，同理1,4和2,5分別代表y軸和z軸的正，負方向.若O點為完全孔隙點，則生成0～5之間的一個隨機正整數(shù)t，例如t=5，那么計算按照z軸的負方向進行，即x、y坐標(biāo)值不變,令k=k-1，若此時點(i,j,k-1)的狀態(tài)仍為孔隙，則繼續(xù)判定(i,j,k-2)點的狀態(tài)，當(dāng)找到(i,j,k-n)點的狀態(tài)為骨架時循環(huán)結(jié)束，此時可以判定(i,j,k-n)點即為孔隙與骨架的邊界點.若計算超出所建模型的范圍且沒找到滿足條件的邊界點，則重新隨機生成一個0～5之間的整數(shù)，即選擇其他方向進行計算.

圖2　隨機方向選擇示意Fig.2　The choice of random direction diagram

選取新系統(tǒng)生成方案對象點與孔隙和骨架的邊界點交換，產(chǎn)生新系統(tǒng)并計算目標(biāo)函數(shù)E1,判斷是否滿足Metropolis準(zhǔn)則，若不滿足則重新產(chǎn)生新系統(tǒng)，若滿足則更新系統(tǒng).判斷是否滿足終止條件，若滿足則輸出優(yōu)化結(jié)果，循環(huán)結(jié)束.補充優(yōu)化算法流程圖見圖3.

圖3　SA補充優(yōu)化算法流程Fig.3　SA added optimization algorithm flow chart

5應(yīng)用實例與效果分析

5.1建模資料的提取

圖4　原始巖心二值化圖Fig.4　The two value map of original core

圖4為某砂巖巖心CT掃描二值化圖像，黑色代表巖石骨架，白色代表孔隙空間，大小為300×300像素點.提取實際巖心建模資料：孔隙度為24.93%.根據(jù)式(3)和式(5)計算數(shù)字巖心建模函數(shù)，并得到建模函數(shù)表達式，如圖5所示自相關(guān)函數(shù)與線性路徑函數(shù)均隨著r1,r2兩點之間距離的增大而減小，當(dāng)r=0，L(0)=S(0)=0.25=φ時,符合自相關(guān)函數(shù)和線性路徑函數(shù)的性質(zhì).

5.2補充優(yōu)化算法參數(shù)設(shè)置

考慮到重建系統(tǒng)與原始參考系統(tǒng)的自相關(guān)函數(shù)和線性路徑函數(shù)在自變量r較小時差異較大，并隨著r值的增大而減小，故設(shè)定αi,βi的值在r

以二維數(shù)字巖心建模為例.如圖6b所示，以(i,j)為中心點，若其周圍8個鄰域點與中心點狀態(tài)值相反的個數(shù)大于n，則中心點即為生成新系統(tǒng)所選取的對象點，由圖6a可以看出棕色的9點鄰域模板，若中心點的狀態(tài)為孔隙，其余的鄰域點中有1個點處于孔隙相，有7個點處于骨架相，且處于骨架相的點數(shù)遠大于孔隙相的點數(shù)，則可判定中心點狀態(tài)為孤立點.本文n的取值分別依次為5～8，n的每次取值對應(yīng)一次系統(tǒng)的遍歷掃描，直到滿足循環(huán)終止條件.系統(tǒng)是否更新取決于系統(tǒng)能量是否降低.對數(shù)字巖心模型遍歷掃描的次數(shù)用循環(huán)終止條來控制，當(dāng)系統(tǒng)能量降低到足夠小時可視為模型達到最優(yōu)(E< 10- 6).這樣第二階段的遍歷掃描運算取代第一階段的隨機搜索，程序在第一階段進行傳統(tǒng)SA隨機搜索優(yōu)化計算時效率低下，而第二階段的遍歷掃描運算則可加快優(yōu)化進程，提高建模效果.

5.3數(shù)字巖心建模結(jié)果

二維數(shù)字巖心建模大小為300×300×300像素，將SA算法中隨機生成的數(shù)字巖心稱為初始巖心,生成的初始巖心孔隙度為24.72%，如圖7a所示孔隙點幾乎均勻分布, 無法反映巖心孔隙空間的任何特征.初始巖心經(jīng)傳統(tǒng)SA 算法處理一段時間后輸出的數(shù)字巖心稱第一階段建模結(jié)果, 傳統(tǒng)SA 算法在建模參考函數(shù)的約束下對初始巖心進行有效的優(yōu)化, 使系統(tǒng)能量達到最低，由圖7b可以看出模型孔隙分布較為紊亂、孔隙空間結(jié)構(gòu)較為獨立、連通性差、含有較多的孤立孔隙點.這是由兩方面原因造成的：第一，模擬退火建模使用的建模函數(shù)有限，只能反映有限的孔隙空間結(jié)構(gòu)特征; 第二，理論上只要運行足夠多的時間，建模巖心就會和真實巖心足夠的接近，但本文為了提高運行效率，只將傳統(tǒng)的SA算法運行一段時間的結(jié)果輸出.經(jīng)補充優(yōu)化算法輸出的數(shù)字巖心稱第二階段建模結(jié)果,如圖7c所示，經(jīng)補充優(yōu)化算法處理后的結(jié)果與實際巖心相比雖不能完全吻合，但從孔隙形態(tài)，孔隙連通性上都有很大改善.三維數(shù)字巖心建模大小為80×80×80像素，建?？紫抖葹?4.93%.由圖8所示，最終數(shù)字巖心模型較第一階段結(jié)果，空隙分布較為集中，連通性較好，孤立點明顯減少，建模效果明顯加強.

圖5　原始巖心的統(tǒng)計函數(shù)(a) 原始巖心自相關(guān)函數(shù); (b) 原始巖心線性路徑函數(shù).Fig.5　Statistical functions of original core (a) Auto-correlation function of original core; (b) Linear path function of original core.

圖6　對象點選取示意圖(a) 孔隙骨架分布； (b) 九點鄰域模板.Fig.6　Target selection diagram(a) The distribution of pore and matrix; (b) The nine point field template.

圖7　二維數(shù)字巖心建模結(jié)果對比(a) 完全隨機的模型; (b) 第一階段建模結(jié)果; (c) 第二階段建模結(jié)果.Fig.7　Comparison chart of 2D digital core modeling results (a) Completely random model; (b) Modeling results from the first stage; (c) Modeling results from the second stage.

圖8　三維數(shù)字巖心建模結(jié)果對比(藍色為孔隙，紅色為骨架)(a) 第一階段數(shù)字巖心模型; (b) 第二階段數(shù)字巖心模型.Fig.8　Comparison chart of 3D digital core modeling results (The blue color represents the pore,the red represents the skeleton) (a) Digital core model from the first stage; (b) Digital core model from the second stage.

變差函數(shù)是與圖像結(jié)構(gòu)相關(guān)，評價圖像性質(zhì)的重要函數(shù)，定義如下：

(10)

式中，γ(h)為變差函數(shù)值，ri為第i個點的坐標(biāo)；f(ri)和f(ri+h)分別為ri及ri+h這兩點處的觀測值，h為兩個觀測點之間的距離.N(h)為距離為h的兩點的數(shù)量，即(x，x+h)的個數(shù).一般來說，對于孔隙度相同的介質(zhì)，函數(shù)曲線的形態(tài)越平緩，介質(zhì)的孔隙連通性就越好.

為了評價數(shù)字巖心建模效果，分別對實際巖心、第一階段建模結(jié)果、第二階段建模結(jié)果的變差函數(shù)對比分析.如圖9所示，在X，Y方向上第一階段建模結(jié)果與原始巖心的變差函數(shù)有一定的偏差.而經(jīng)第二階段處理后模型顯著改善, 其變差函數(shù)與真實巖心的幾乎完全吻合，說明最終的模型更能反映真實巖心的孔隙空間結(jié)構(gòu)特征.

圖9　數(shù)字巖心建模變差函數(shù)對比(a) X方向變差函數(shù)圖; (b) Y方向變差函數(shù)圖.Fig.9　Comparison chart of variation function for digital core modeling(a) The variation function for X direction; (b) The variation function for Y direction.

6結(jié)論

(1) 本文定義了孔隙骨架邊界點的選取標(biāo)準(zhǔn)，并以擇多算子選取的對象點與邊界點交換作為新系統(tǒng)的生成方案，設(shè)計了用遍歷掃描搜索代替?zhèn)鹘y(tǒng)隨機搜索的后續(xù)補充優(yōu)化方案，取得良好的效果.

(2) 基于擇多算子選取孤立點參數(shù)設(shè)置十分重要，n值依次為5～8，對n為不同的取值，依次進行優(yōu)化計算.與傳統(tǒng)的SA建立數(shù)字巖心方法相比優(yōu)化了建模效率，從建模結(jié)果上看孔隙空間形態(tài)和孔隙連通性上都有較大改善.

(3) 對于其他的隨機法建立數(shù)字巖心方法，如多點統(tǒng)計法、過程法等，可用本文的設(shè)計的基于模擬退火的后續(xù)優(yōu)化方案來進行更細致的優(yōu)化，具有很好的實用性.

致謝在本文研究和寫作過程中，得到了吉林大學(xué)李舟波教授的關(guān)心和建議，在此表示誠摯謝意.

References

Bakke S, ?ren P E. 1997. 3-D pore-scale modelling of sandstones and flow simulations in the pore networks.SPEJournal, 2(2): 136-149.

Bakke, S, Oren, P. E. 2003. Reconstruction of Berea sandstone and pore-scale modelling of wettability effects.JournalofPetroleumScienceandEngineering, 39: 177-199.

Hazlett R D. 1997. Statistical characterization and stochastic modeling of pore networks in relation to fluid flow.MathematicalGeology, 29(6): 801-822.Hidajat I, Rastogi A, Singh M. Transport properties of porous media from thin section.SPEJournal. 2002, 7(1): 40-48.

Joshi M. A class three-dimensional modeling technique for studying porous media [Ph. D Thesis]. Kansas: University of Kansas, 1974.Keehm Y. Computational Rock Physics: Transport Properties in Porous Media and Applications. Stanford: Stanford University, 2003.Liu X F. 2010. Numerical simulation of elastic and electrical properties of rock based on digital core[Ph. D. thesis] (in Chinese). Shandong: China University of Petroleum (Huadong). Liu X F, Zhang W W, Sun J M. 2013. Methods of constructing 3-D digital cores: A review.ProgressinGeophysics(in Chinese), 28(6): 3066-3072, doi: 10.6038/pg20130630.

Lymberopoulos D P, Payatakes A C. Derivation of topological, geometrical, and correlational properties of porous media from pore-chart analysis of serial section data.JournalofColloidandInterfaceScience, 1992, 150(1): 61-80.

Metropolis N, Rosenbluth A W, Rosenbluth M N. 1953. Equation of state calculations by fast computing machines.JournalofChemicalPhysics, 21(6): 1087-1092.

Michalewicz Z. 1995. A survey of constraint handling techniques in evolutionary computation methods.∥ Proceedings of the 4th Annual Conference on Evolutionary Programming. 135-155. Michalewicz Z, Dasgupta D, Le Riche R G, et al. 1996. Evolutionary algorithms for constrained Engineering problems.ComputersandIndustrialEngineering, 30(4): 851-870. Nie X. 2014. Modeling and numerical simulation of digital core of shale gas reservoir [Ph. D. thesis] (in Chinese). Beijing: China University of Geosciences (Beijing). Okabe H, Blunt M J. 2005. Pore space reconstruction using multiple-point statistics.JournalofPetroleumScienceandEngineering, 46(1-2): 121-137.

?ren P E, Bakke S. 2002. Process based reconstruction of sandstones and prediction of transport properties.TransportinPorousMedia, 46(2-3): 311-343.

Quiblier J A. 1984. A new three-dimensional modeling technique for studying porous media.JournalofColloidandInterfaceScience, 98(1): 84-102.

Wu K, Nunan N, Crawford J W, et al. 2004. An efficient Markov Chain model for the simulation of heterogeneous soil structure.SoilScienceSocietyofAmerica. 68:346-351.Yao J, Zhao X C. 2010. Theoretical Simulation, Digital Cores and Pore Level Seepage (in Chinese). Beijing: Petroleum Industry Press. Zhang T, Li D L, Lu D T, et al. 2010. Research on the reconstruction method of porous media using multiple-point geostatistics.ScienceChinaPhysics,Mechanics&Astronomy, 53(1): 122-134.

Zhao X C, Yao J, Yi Y J, et al. 2008. A new method of constructing digital core utilizing stochastic search algorithm based on majority operator.RockandSoilMechanics(in Chinese), 29(5): 1339-1344.

附中文參考文獻

劉學(xué)鋒. 2010. 基于數(shù)字巖心的巖石聲電特性微觀數(shù)值模擬研究

[博士學(xué)位論文]. 東營: 中國石油大學(xué)(華東).

劉學(xué)鋒, 張偉偉, 孫建孟. 2013. 三維數(shù)字巖心建模方法綜述. 地球物理學(xué)進展, 28(6): 3066-3072, doi: 10.6038/pg20130630.

聶昕. 2014. 頁巖氣儲層巖石數(shù)字巖心建模及導(dǎo)電性數(shù)值模擬研究[博士學(xué)位論文]. 北京: 中國地質(zhì)大學(xué)(北京).

王晨晨, 姚軍, 楊永飛等. 2013. 碳酸鹽巖雙孔隙數(shù)字巖心結(jié)構(gòu)特征分析. 中國石油大學(xué)學(xué)報(自然科學(xué)版), 37(2): 71-74.

姚軍, 趙秀才, 衣艷靜等. 2005. 數(shù)字巖心技術(shù)現(xiàn)狀及展望. 油氣地質(zhì)與采收率, 12(6): 52-54.

姚軍, 趙秀才. 2010. 數(shù)字巖心及孔隙級滲流模擬理論. 北京: 石油工業(yè)出版社.

張挺, 李道倫, 盧德唐等. 2009. 基于多點地質(zhì)統(tǒng)計法的多孔介質(zhì)重構(gòu)研究. 中國科學(xué)G輯, 39(9): 1348-1360.

趙秀才, 姚軍, 衣艷靜等. 2008. 基于擇多算子的隨機搜索法建立數(shù)字巖心的新技術(shù). 巖土力學(xué), 29(5): 1339-1344.

(本文編輯汪海英)

A complement optimization scheme to establish the digital core model based on the simulated annealing method

MO Xiu-Wen1, ZHANG Qiang1, LU Jing-An2

1CollegeofGeo-ExplorationofScienceandTechnology,JilinUniversity,Changchun130026,China2GuangzhouMarineGeologicalSurvey,Guangzhou510760,China

AbstractThe purpose of this work is to overcome the disadvantages of low efficiency and much isolated points in building the digital core model by the traditional simulated annealing methods. We design the subsequent search scheme, defines the criterion to select the boundary point of pore and the matrix, and proposes the generation solutions of the new system by exchanging the isolated pore points and boundary points of pore and matrix based on the majority operator method. Using traditional simulated annealing algorithm to optimize the random initial model as the first stage, and using the optimizing plan generated by the new system as the second stage, we can finally reduce the isolated points of pore or matrix in the output model. A preliminary test shows the new optimization solution is feasible and very effective. For instance the results from the final optimization model are relatively closer to the original cores compared with that from traditional simulated annealing modeling, both in pore shape and connectivity. The variation function, which characterizes the rock pore space, also demonstrates that the optimization results of the complement search scheme are superior to that of the traditional simulated annealing methods.KeywordsSimulated annealing methods； Digital core; Auto-correlation function；Linear path function； Variation function

基金項目國家高技術(shù)研究發(fā)展計劃(2013AA092501)，國家自然科學(xué)基金項目(40874057)資助.

作者簡介莫修文，男，1970年生，教授，1998年畢業(yè)于長春科技大學(xué)地球物理系，獲工學(xué)博士學(xué)位，現(xiàn)從事地球物理測井教學(xué)與科研工作. E-mail:moxw@jlu.edu.cn

doi:10.6038/cjg20160526 中圖分類號P631

收稿日期2015-02-05，2016-03-17收修定稿

莫修文, 張強， 陸敬安. 2016. 模擬退火法建立數(shù)字巖心的一種補充優(yōu)化方案. 地球物理學(xué)報，59(5):1831-1838,doi:10.6038/cjg20160526.

Mo X W, Zhang Q, Lu J A. 2016. A complement optimization scheme to establish the digital core model based on the simulated annealing method.ChineseJ.Geophys. (in Chinese),59(5):1831-1838,doi:10.6038/cjg20160526.