基于細(xì)胞自動(dòng)機(jī)的碳酸鹽巖沉積模擬

趙華偉 段太忠 劉彥鋒 廉培慶 張文彪 趙磊

(中國(guó)石化石油勘探開發(fā)研究院， 北京 100083)

在油氣藏地質(zhì)建模過程中，基于變差函數(shù)的兩點(diǎn)地質(zhì)統(tǒng)計(jì)學(xué)建模方法是一種主流地質(zhì)建模方法[1]。近年來，基于訓(xùn)練圖像的多點(diǎn)地質(zhì)統(tǒng)計(jì)學(xué)建模方法也逐漸得到發(fā)展應(yīng)用[2-3]。這類基于統(tǒng)計(jì)學(xué)的建模方法最明顯的缺陷是，對(duì)井點(diǎn)數(shù)據(jù)的依賴性極強(qiáng)，在數(shù)據(jù)較少時(shí)難以獲得滿意的結(jié)果。基于沉積過程的建模方法(沉積模擬方法)，是通過模擬真實(shí)的沉積歷史來建立地質(zhì)模型，其結(jié)果更加符合地質(zhì)規(guī)律，有望發(fā)展成為新的主流地質(zhì)建模方法[4-5]。沉積模擬方法可模擬碎屑巖、碳酸鹽巖及其混合物的沉積過程。其基本模擬方法包括基于擴(kuò)散方程及水動(dòng)力學(xué)方程的確定性模擬方法，以及基于兩點(diǎn)隨機(jī)游走或細(xì)胞自動(dòng)機(jī)(Cellular Automaton，CA)的隨機(jī)模擬方法。

細(xì)胞自動(dòng)機(jī)是一種離散迭代演化模型。它是由無限或有限個(gè)相同的網(wǎng)格組成，每個(gè)網(wǎng)格處于一種有限狀態(tài)。網(wǎng)格在t時(shí)刻的狀態(tài)，由t-1時(shí)的有限集合的狀態(tài)來確定。在每次演化過程中，所有網(wǎng)格遵守同一規(guī)則，同時(shí)演進(jìn)[6-7]。基于細(xì)胞自動(dòng)機(jī)的建模方法應(yīng)用較為廣泛[8]，近年來開始應(yīng)用于濁積巖、碳酸鹽巖等儲(chǔ)層模擬[9-12]。下面將詳細(xì)描述該方法的碳酸鹽巖沉積模擬過程，分析其優(yōu)缺點(diǎn)，探討其發(fā)展趨勢(shì)和應(yīng)用前景。

1 細(xì)胞自動(dòng)機(jī)基本原理

一個(gè)細(xì)胞自動(dòng)機(jī)系統(tǒng)(A)通常包含4個(gè)元素，分別是細(xì)胞空間(Ld)、細(xì)胞狀態(tài)集合(S)、包括中心細(xì)胞在內(nèi)的鄰域內(nèi)細(xì)胞組合(N)，以及對(duì)應(yīng)的演化規(guī)則R。其關(guān)系如式(1)所示：

A=(Ld，S，N，R)

(1)

細(xì)胞在空間上的網(wǎng)格點(diǎn)集合就是細(xì)胞空間。目前的細(xì)胞自動(dòng)機(jī)研究，主要方向側(cè)重于一維細(xì)胞機(jī)和二維細(xì)胞機(jī)研究。一維細(xì)胞自動(dòng)機(jī)的細(xì)胞空間只有一種，而二維細(xì)胞自動(dòng)機(jī)的細(xì)胞空間通常可按三角形、正方形或六邊形這3種網(wǎng)格進(jìn)行排列。其中，二維細(xì)胞自動(dòng)機(jī)的細(xì)胞空間中，正方形網(wǎng)格的應(yīng)用最為廣泛，其邊界條件分為周期型、反射型、定值型等不同類型。

細(xì)胞的狀態(tài)可以是二進(jìn)制形式(通常用0或1表示)，也可以是整數(shù)形式的離散集。理論上，細(xì)胞自動(dòng)機(jī)的細(xì)胞只能有1個(gè)狀態(tài)變量，但實(shí)際應(yīng)用中其狀態(tài)變量往往有多個(gè)。

在細(xì)胞自動(dòng)機(jī)中，演化規(guī)則是定義在局部作用范圍之內(nèi)的，即一個(gè)細(xì)胞在下一個(gè)時(shí)刻的狀態(tài)取決于其本身及其鄰域細(xì)胞的狀態(tài)。在一維細(xì)胞自動(dòng)機(jī)中，通常以半徑r來確定鄰域，距離某個(gè)細(xì)胞半徑內(nèi)的所有細(xì)胞，均被認(rèn)為是該細(xì)胞的鄰域。二維細(xì)胞自動(dòng)機(jī)的鄰域相對(duì)比較復(fù)雜，通常包括馮·諾依曼型、摩爾型、擴(kuò)展的摩爾型。如圖1所示，3種類型二維細(xì)胞機(jī)中心細(xì)胞(網(wǎng)絡(luò))周圍的陰影區(qū)域，即中心細(xì)胞的領(lǐng)域。

根據(jù)細(xì)胞當(dāng)前狀態(tài)及其鄰域狀態(tài)確定下一時(shí)刻該細(xì)胞狀態(tài)的動(dòng)力學(xué)函數(shù)，此過程稱為演化規(guī)則，此函數(shù)即狀態(tài)轉(zhuǎn)移函數(shù)。狀態(tài)轉(zhuǎn)移函數(shù)可表示為：

(2)

式中：R—— 狀態(tài)轉(zhuǎn)移函數(shù)，即演化規(guī)則；

演化規(guī)則是事先給出用于約束細(xì)胞自動(dòng)機(jī)狀態(tài)的條件集合。在實(shí)際應(yīng)用中，細(xì)胞自動(dòng)機(jī)模型設(shè)計(jì)是否成功，關(guān)鍵在于其演化規(guī)則是否合理，能否客觀反映系統(tǒng)的本質(zhì)特征。演化規(guī)則的設(shè)計(jì)是細(xì)胞自動(dòng)機(jī)模擬的核心環(huán)節(jié)。

圖1 二維細(xì)胞自動(dòng)機(jī)鄰域類型

2 基于細(xì)胞自動(dòng)機(jī)的碳酸鹽巖沉積模擬

當(dāng)前，細(xì)胞自動(dòng)機(jī)算法主要應(yīng)用在兩類儲(chǔ)層的沉積模擬中，分別是碳酸鹽巖和濁積巖儲(chǔ)層。其中，以Burgess、Teles等人的研究工作最為典型[11-12]。在此，以Burgess的研究工作為背景，介紹碳酸鹽巖的沉積模擬過程。

2.1 模擬過程

2.1.1 巖相平面演化

碳酸鹽巖的巖相平面演化遵循了細(xì)胞自動(dòng)機(jī)的基本規(guī)則，即某網(wǎng)格當(dāng)前時(shí)間步的巖相由上一時(shí)間步該網(wǎng)格及其鄰域范圍內(nèi)網(wǎng)格的巖相狀態(tài)所決定。巖相平面演化的鄰域范圍采用了擴(kuò)展的摩爾型鄰域，如圖2所示。巖相平面演化規(guī)則中，按照鄰域中與該巖相類型相同的鄰居數(shù)量n，分為式(3)、式(4)、式(5)所示不同條件：

4≤n≤10

(3)

n≤3或n≥11

(4)

6≤n≤10

(5)

當(dāng)某網(wǎng)格的狀態(tài)為非孔洞，且該網(wǎng)格及其鄰域內(nèi)網(wǎng)格中某巖相的數(shù)量滿足式(3)時(shí)，該網(wǎng)格的巖相類型保持不變。此條件的含義是，周圍的鄰居數(shù)量適宜該巖相繼續(xù)生存。

當(dāng)某網(wǎng)格的狀態(tài)為非孔洞，且該網(wǎng)格及其鄰域內(nèi)某巖相的數(shù)量滿足式(4)時(shí)，該巖相變?yōu)榭锥础４藯l件的含義是，周圍的同類巖相數(shù)量過少，或競(jìng)爭(zhēng)過于激烈，使得該巖相無法繼續(xù)存在。

圖2 巖相平面演化擴(kuò)展摩爾型鄰域

當(dāng)某網(wǎng)格的狀態(tài)為孔洞時(shí)，該網(wǎng)格及其鄰域內(nèi)某巖相的數(shù)量滿足式(5)時(shí)，該孔洞變?yōu)閹r相。此條件的含義是，周圍的環(huán)境適宜新巖相生成。

需要說明的是，孔洞具體變?yōu)楹畏N巖相是與時(shí)間變化相關(guān)的。改變每次檢索的順序，可以使各相的轉(zhuǎn)化概率相等。在圖3所示隨機(jī)初始狀態(tài)下的巖相平面演化過程中，可以看出，隨著演化時(shí)間的延長(zhǎng)，同種巖相逐漸聚集，整個(gè)模擬區(qū)域的熵降低，系統(tǒng)逐漸進(jìn)入有序狀態(tài)。

2.1.2 巖相垂向生長(zhǎng)

將3類不同巖相的垂向生長(zhǎng)速率與水深關(guān)系進(jìn)行了對(duì)比，其關(guān)系曲線如圖4所示。其中，巖相1的微生物為淺水喜光微生物，在深水中生長(zhǎng)速率很低，但是在淺水中生長(zhǎng)速率極高；與巖相1的微生物相比，巖相2的微生物能夠在深水中生長(zhǎng)，但是在淺水中的最大生長(zhǎng)速率較低；巖相3的微生物生長(zhǎng)速率不受水深的影響，保持勻速的狀態(tài)。

圖3 巖相平面演化結(jié)果

圖4 3類巖相垂向生長(zhǎng)速率與水深關(guān)系曲線

Bosscher、Schlager在研究中給出了描述碳酸鹽巖生長(zhǎng)速率的函數(shù)，認(rèn)為碳酸鹽巖垂向生長(zhǎng)速率與水深存在雙曲正切關(guān)系[13]：

(6)

式中：g(w)—— 垂向生長(zhǎng)速率；

gm—— 最大生長(zhǎng)速率；

I0—— 表層光照強(qiáng)度；

Ic—— 飽和光照強(qiáng)度；

k—— 速率常數(shù)；

w——水深。

根據(jù)式(6)，可以進(jìn)一步描述碳酸鹽巖的累積速率：

υ(w)=υmtanh(kexp(dw(t)))

(7)

其中：υ(w)—— 累積速率；

υm—— 最大累積速率；

d—— 衰減常數(shù)。

如果進(jìn)一步考慮周圍巖相的影響，則某網(wǎng)格中巖相的真實(shí)垂向生長(zhǎng)速率a可以用式(8)所示方程組表示：

(8)

式中：nmin—— 該巖相在下一時(shí)間步保持狀態(tài)不變的最小網(wǎng)格數(shù)量；

nmax—— 該巖相在下一時(shí)間步保持狀態(tài)不變的最大網(wǎng)格數(shù)量；

noptimum—— 保證最大生長(zhǎng)速率的網(wǎng)格數(shù)量。

根據(jù)式(8)可知，隨著鄰域內(nèi)同巖相網(wǎng)格的數(shù)量逐漸增加，該網(wǎng)格中巖相的垂向生長(zhǎng)速率也不斷增大。當(dāng)達(dá)到最大生長(zhǎng)速率后，隨著同巖相網(wǎng)格數(shù)量進(jìn)一步增加，該網(wǎng)格中巖相的垂向生長(zhǎng)速率又逐漸降低。巖相的垂向變化，同時(shí)受縱向水深和平面狀態(tài)變化的影響。

2.1.3 沉積物搬運(yùn)及相對(duì)海平面變化

盡管大部分碳酸鹽巖臺(tái)地是由生物生長(zhǎng)而形成，但沉積物搬運(yùn)作用對(duì)于碳酸鹽巖臺(tái)地生長(zhǎng)的影響不能忽略。沉積物搬運(yùn)過程可以用“梯度算法”簡(jiǎn)化表示。如果某網(wǎng)格的海拔高度高于周圍的網(wǎng)格，且其周圍存在孔洞網(wǎng)格，則該網(wǎng)格產(chǎn)生的沉積物一部分留在該網(wǎng)格，剩余部分均勻地分配給周圍海拔高度較低的孔洞網(wǎng)格；反之，其產(chǎn)生的沉積物留在原地。在此，以圖5所示沉積物運(yùn)移過程為例進(jìn)行分析。當(dāng)沉積物不發(fā)生運(yùn)移時(shí)，中心網(wǎng)格厚度為0.5 m，周圍4個(gè)網(wǎng)格的海拔均高于該網(wǎng)格；當(dāng)沉積物發(fā)生運(yùn)移時(shí)，周圍4個(gè)網(wǎng)格的海拔較低，則沉積物均勻地分配給周圍的4個(gè)網(wǎng)格，按照上述規(guī)則進(jìn)一步搬運(yùn)，直至厚度小于0.5 cm。

圖5 沉積物運(yùn)移過程示意圖

在模擬過程中，假設(shè)整個(gè)模擬區(qū)域以恒定不變的速率發(fā)生沉降，海平面的周期性變化采用正弦曲線表示。當(dāng)相對(duì)海平面發(fā)生變化，導(dǎo)致某一網(wǎng)格出露水面時(shí)，巖相生長(zhǎng)即中止；當(dāng)某網(wǎng)格再次沉入水下時(shí)，巖相生長(zhǎng)則恢復(fù)。值得注意的是，相對(duì)海平面的變化可能導(dǎo)致鄰域內(nèi)同巖相的網(wǎng)格數(shù)量發(fā)生變化，進(jìn)而影響巖相的平面演化和垂向生長(zhǎng)。這一點(diǎn)與真實(shí)的碳酸鹽巖臺(tái)地生長(zhǎng)過程相類似。

2.2 模擬結(jié)果

模擬區(qū)域的平面大小為50×50網(wǎng)格，單個(gè)網(wǎng)格的長(zhǎng)度為1 km。網(wǎng)格有4種狀態(tài)，即巖相1、巖相2、巖相3和孔洞。在初始時(shí)刻，整個(gè)模擬區(qū)域的巖相為隨機(jī)分布。單個(gè)時(shí)間步長(zhǎng)為1 000 a，模擬1 000個(gè)時(shí)間步，總時(shí)間步長(zhǎng)為1 Ma。分別針對(duì)3種情形下的沉積過程進(jìn)行模擬：(1) 海平面不變，無沉積物運(yùn)移；(2) 海平面變化，無沉積物運(yùn)移；(3) 海平面不變，有沉積物運(yùn)移。沉積模擬結(jié)果如圖6所示。

在第1種情形下，不同巖相的平面演化和垂向生長(zhǎng)形成了復(fù)雜的巖相結(jié)構(gòu)。此時(shí)，雖然模擬結(jié)束時(shí)每個(gè)網(wǎng)格的臺(tái)地高度基本相同，但整個(gè)過程中仍存在高度不同的情況。

在第2種情形下，海平面發(fā)生周期性變化，碳酸鹽巖臺(tái)地在相當(dāng)長(zhǎng)的時(shí)間內(nèi)出露水面，導(dǎo)致很長(zhǎng)時(shí)間的空白沉積期。與此同時(shí)，周期性的海平面升降也影響著巖相的平面演化，使得孔洞大量出現(xiàn)。2種機(jī)制共同作用的結(jié)果是，當(dāng)模擬結(jié)束時(shí)，碳酸鹽巖臺(tái)地的高度出現(xiàn)顯著差別。

在第3種情形下，沉積物搬運(yùn)作用以巖相4和巖相5來表示。由其剖面圖和柱狀圖可知，搬運(yùn)的沉積物具有一定厚度，不過其厚度明顯小于原地生長(zhǎng)的巖相。與第1種情形相類似，模擬結(jié)束時(shí)碳酸鹽巖臺(tái)地的高度基本一致。

3 細(xì)胞自動(dòng)機(jī)方法的優(yōu)缺點(diǎn)及應(yīng)用前景

3.1 優(yōu)缺點(diǎn)

與基于定量方程的確定性沉積模擬方法相比，該方法的優(yōu)點(diǎn)是所采用的規(guī)則十分簡(jiǎn)單，能夠保證極高的計(jì)算效率。與此同時(shí)，通過多次迭代演化，得到復(fù)雜結(jié)果，在一定程度上能夠重現(xiàn)復(fù)雜的地質(zhì)過程。

該方法較明顯的缺點(diǎn)是，初始時(shí)刻的巖相分布對(duì)最終結(jié)果影響顯著。如圖7所示，在初始時(shí)刻，2種情形下整個(gè)模擬區(qū)域僅有1個(gè)網(wǎng)格的巖相不同。模擬持續(xù)一段時(shí)間后，2種情形下平面模擬結(jié)果開始顯現(xiàn)出不同。當(dāng)模擬結(jié)束時(shí)，2種情形下的平面模擬結(jié)果出現(xiàn)了顯著差別，同時(shí)垂向的巖相分布也有著明顯的不同。因此，選擇合適的初始巖相分布有助于獲得有價(jià)值的模擬結(jié)果。

平面的演化規(guī)則決定了巖相的平面變化。在此，針對(duì)保持巖相不變需要的鄰域中相同巖相的最小數(shù)量和最大數(shù)量作了定義，如式(3)所示。同時(shí)，針對(duì)保持巖相最大生長(zhǎng)速率的同巖相數(shù)量作了定義，如式(8)所示。這些定義均具有很強(qiáng)的主觀性，所定義數(shù)量是否符合真實(shí)的碳酸鹽巖臺(tái)地生長(zhǎng)過程，還有待進(jìn)一步研究。除此之外，對(duì)于沉積物的搬運(yùn)及相對(duì)海平面變化都作了簡(jiǎn)化處理，其過程有待深入討論。

3.2 應(yīng)用前景

就結(jié)果的可靠性來說，目前基于細(xì)胞自動(dòng)機(jī)的沉積模擬遠(yuǎn)不如基于定量方程的確定性沉積模擬。但是，如果能夠避免上述缺點(diǎn)，該方法就具有良好的應(yīng)用前景。具體思路是：(1) 用沉積初始時(shí)刻的巖相古地理作為模擬的初始條件；(2) 設(shè)計(jì)合適的具有地質(zhì)意義的平面演化規(guī)則；(3) 進(jìn)一步優(yōu)化沉積物搬運(yùn)的算法。

圖6 3種情形下的沉積模擬結(jié)果

圖7 巖相分布模擬對(duì)比

運(yùn)用基于定量方程的沉積模擬方法時(shí)，需要設(shè)置大量的參數(shù)，以精細(xì)刻畫各種沉積物搬運(yùn)和沉積過程，如沉積物沉降/抬升、沉積物供給變化。大量的參數(shù)難以全部準(zhǔn)確獲取，因此基于定量方程的沉積模擬方法存在一定的不確定性。此時(shí)，基于細(xì)胞自動(dòng)機(jī)的沉積模擬方法可以有效模擬不確定性沉積過程。我們可以運(yùn)用確定性沉積模擬方法來模擬沉積亞相和較明確微相的變化，而對(duì)于部分不明確的沉積微相，則采用基于細(xì)胞自動(dòng)機(jī)的沉積模擬方法來進(jìn)行刻畫。對(duì)于溶蝕、壓溶等控制方程不明確的后期成巖過程，也可以運(yùn)用基于細(xì)胞自動(dòng)機(jī)的沉積模擬方法來進(jìn)行模擬。將確定性模擬方法和細(xì)胞自動(dòng)機(jī)模擬方法結(jié)合起來，有望獲得比較滿意的結(jié)果。

4 結(jié) 語

基于細(xì)胞自動(dòng)機(jī)的沉積模擬方法是一種隨機(jī)模擬方法，可由鄰域范圍內(nèi)的網(wǎng)格狀態(tài)決定網(wǎng)格的演化方向。與確定性沉積模擬方法相比，該方法的最大優(yōu)點(diǎn)是，能夠通過簡(jiǎn)單的演化規(guī)律模擬復(fù)雜的沉積過程。基于細(xì)胞自動(dòng)機(jī)的沉積模擬方法可同時(shí)開展平面演化和垂向生長(zhǎng)，并將沉積物搬運(yùn)作用和相對(duì)海平面變化的影響也考慮了進(jìn)去。該方法涵蓋了碳酸鹽巖沉積的所有因素，基本上能夠反映碳酸鹽巖的沉積過程。該方法的主要優(yōu)點(diǎn)是，規(guī)則簡(jiǎn)單，計(jì)算效率高；主要缺點(diǎn)是，初始狀態(tài)對(duì)結(jié)果的影響顯著，且演化規(guī)則的地質(zhì)意義不明確，仍有待進(jìn)一步研究。該方法可以作為確定性沉積模擬的一種有效補(bǔ)充。