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

莫修文, 張強(qiáng)， 陸敬安

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

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

莫修文1, 張強(qiáng)1， 陸敬安2

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

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

1引言

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

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

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

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

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

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

(1)

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

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

2.2兩點(diǎn)概率函數(shù)

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

(2)

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

(3)

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

(4)

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

2.3線性路徑函數(shù)

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

(5)

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

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

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

3.1模擬退火法的基本理論

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

(6)

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

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

(7)

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

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

(8)

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

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

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

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

(9)

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

數(shù)字巖心建模步驟如下： (1) 對(duì)實(shí)際巖心CT圖像或電鏡掃描圖像進(jìn)行二值化處理，統(tǒng)計(jì)建?；举Y料φ、S(r)、L(r); (2) 隨機(jī)生成孔隙度為φ的三維數(shù)字巖心(0代表巖石骨架，1代表孔隙)，并計(jì)算初始模型目標(biāo)函數(shù)E，進(jìn)入優(yōu)化進(jìn)程; (3) 在孔隙和巖石骨架空間任意選取一孔隙點(diǎn)和骨架點(diǎn)，兩者互換位置生成新系統(tǒng)，并計(jì)算目標(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é)束進(jìn)程.

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

上述傳統(tǒng)模擬退火法建立數(shù)字巖心，開始時(shí)對(duì)完全隨機(jī)的模型進(jìn)行優(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)逐步形成，但分布較為分散，連通性差，孤立的孔隙點(diǎn)或骨架點(diǎn)較多.如果再按傳統(tǒng)的隨機(jī)交換骨架點(diǎn)和孔隙點(diǎn)作為新系統(tǒng)的生成方案，就會(huì)導(dǎo)致原有已形成的孔隙結(jié)構(gòu)遭到破壞，孔隙結(jié)構(gòu)陷入形成-破壞-再形成的惡性循環(huán)，無形中增加了不必要的計(jì)算量.而孔隙中的骨架點(diǎn)和骨架中孤立的孔隙點(diǎn)都是不可能大量存在的，新系統(tǒng)生成方案選取的對(duì)象點(diǎn)主要以孤立孔隙點(diǎn)和骨架點(diǎn)為目標(biāo).故本文設(shè)計(jì)了兩個(gè)階段的優(yōu)化方案，第一階段用傳統(tǒng)的SA方法對(duì)初始的隨機(jī)模型做優(yōu)化.第二階段以擇多算子選取的對(duì)象點(diǎn)與孔隙骨架的邊界點(diǎn)交換生成新系統(tǒng)，作為補(bǔ)充優(yōu)化階段.

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

新系統(tǒng)生成方案的對(duì)象點(diǎn)選取標(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)值與中心體素相反的個(gè)數(shù)大于等于n，則將該中心體素狀態(tài)判定為新系統(tǒng)生成方案選取的對(duì)象點(diǎn).n的取值根據(jù)具體情況而定，若n取值過高，會(huì)導(dǎo)致符合條件的點(diǎn)越來越難選取，從而使程序在尋找孤立點(diǎn)上花費(fèi)過多時(shí)間.若n的取值過低，則相當(dāng)于傳統(tǒng)的系統(tǒng)生成方案，導(dǎo)致優(yōu)化后模型與實(shí)際巖心孔隙形態(tài)相差較大.

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

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

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

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

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

5.1建模資料的提取

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

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

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

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

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

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

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

圖5　原始巖心的統(tǒng)計(jì)函數(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　對(duì)象點(diǎn)選取示意圖(a) 孔隙骨架分布； (b) 九點(diǎn)鄰域模板.Fig.6　Target selection diagram(a) The distribution of pore and matrix; (b) The nine point field template.

圖7　二維數(shù)字巖心建模結(jié)果對(duì)比(a) 完全隨機(jī)的模型; (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é)果對(duì)比(藍(lán)色為孔隙，紅色為骨架)(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)，評(píng)價(jià)圖像性質(zhì)的重要函數(shù)，定義如下：

(10)

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

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

圖9　數(shù)字巖心建模變差函數(shù)對(duì)比(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) 本文定義了孔隙骨架邊界點(diǎn)的選取標(biāo)準(zhǔn)，并以擇多算子選取的對(duì)象點(diǎn)與邊界點(diǎn)交換作為新系統(tǒng)的生成方案，設(shè)計(jì)了用遍歷掃描搜索代替?zhèn)鹘y(tǒng)隨機(jī)搜索的后續(xù)補(bǔ)充優(yōu)化方案，取得良好的效果.

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

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

致謝在本文研究和寫作過程中，得到了吉林大學(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.

附中文參考文獻(xiàn)

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

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

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

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

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

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

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

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

趙秀才, 姚軍, 衣艷靜等. 2008. 基于擇多算子的隨機(jī)搜索法建立數(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

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

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

doi:10.6038/cjg20160526 中圖分類號(hào)P631

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

莫修文, 張強(qiáng)， 陸敬安. 2016. 模擬退火法建立數(shù)字巖心的一種補(bǔ)充優(yōu)化方案. 地球物理學(xué)報(bào)，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.