基于馬爾柯夫地層序列分析的分布條件模擬方法及其應(yīng)用

0數(shù)值模擬方法receive模擬是計(jì)算區(qū)域地下水資源的基本方法。在數(shù)值建模過程中，應(yīng)注意非均質(zhì)和水文結(jié)構(gòu)的問題。水文地質(zhì)參數(shù)空間分布規(guī)律的認(rèn)識或刻畫程度決定了模擬結(jié)果的準(zhǔn)確性與可靠性。作為影響儲水介質(zhì)非均質(zhì)性的重要因素,巖性的空間分布規(guī)律是解決儲水介質(zhì)非均質(zhì)性的重要途徑,特別是在松散巖類為主的沉積平原地區(qū)。在國外,利用條件模擬方法進(jìn)行巖相空間分布規(guī)律研究已有較廣泛的發(fā)展,國內(nèi)目前仍處于起步階段。條件模擬在地質(zhì)統(tǒng)計(jì)學(xué)中占有十分重要的位置,在許多行業(yè),尤其是石油儲層建模研究中得到廣泛應(yīng)用。條件模擬方法在保持變量的空間自相關(guān)性函數(shù)不變的情況下使觀測點(diǎn)處的模擬值等于實(shí)測值,進(jìn)而實(shí)現(xiàn)三維空間的模擬。在建模過程中經(jīng)常采用實(shí)測數(shù)據(jù)對空間變量進(jìn)行結(jié)構(gòu)化分析,建立典型參數(shù)(如滲透率)的變差函數(shù),或采用Markov鏈進(jìn)行地層序列分析,利用轉(zhuǎn)移概率矩陣替代傳統(tǒng)變差函數(shù),利用隨機(jī)模擬方法建立二維或三維石油儲層模型。常用的隨機(jī)模擬方法有序貫?zāi)M、分形和轉(zhuǎn)向帶模擬、模擬退火、概率場模擬、混合方法等。本文將簡要介紹條件模擬方法基本原理,以華北平原為研究區(qū)域,模擬區(qū)內(nèi)松散巖類的空間分布,并對模擬結(jié)果進(jìn)行分析。1條件模擬的基本原則1.1馬爾柯夫轉(zhuǎn)移率矩陣法在地質(zhì)學(xué)理論中,研究區(qū)全部巖性所構(gòu)成的集合稱為巖相空間或相空間(Sk,k=1,2,…K),在地層剖面上巖相的變化過程具有馬爾柯夫性質(zhì)。垂向馬爾柯夫轉(zhuǎn)移概率矩陣采用矩陣指數(shù)形式表示:Τ(hΦ)=exp(RΦhΦ)(1)其中:Τ(hΦ)=[t11(hΦ)?t1k(hΦ)??tk1(hΦ)?tkk(hΦ)]；RΦ=[r11,Φ?r1k,Φ??rk1,Φ?rkk,Φ]式中:T表示轉(zhuǎn)移概率矩陣,tjk(hΦ)代表在Φ方向上、單位長度內(nèi),由狀態(tài)j轉(zhuǎn)移到狀態(tài)k的概率;R表示轉(zhuǎn)移率矩陣;hΦ表示在Φ方向上的采樣間隔,rjk,Φ代表在Φ方向上、單位長度內(nèi),由狀態(tài)j轉(zhuǎn)移到狀態(tài)k的變化率。對轉(zhuǎn)移概率矩陣進(jìn)行特征值分析,利用下式得到空間連續(xù)狀態(tài)離散的轉(zhuǎn)移概率矩陣中各元素。tjk(h)=pk+k-1∑i=1zjk,i+1exp(λi+1h)(2)式中:pk表示狀態(tài)k的極限概率,λi和Zi分別表示轉(zhuǎn)移率矩陣R的特征值(即矩陣的譜)和特征值對應(yīng)的譜元素矩陣。利用垂向鉆孔剖面資料,運(yùn)用馬爾柯夫轉(zhuǎn)移概率進(jìn)行地層序列分析,建立垂直方向轉(zhuǎn)移概率矩陣。根據(jù)Walther沉積相序定律,由垂向?qū)有蛴?jì)算得到的轉(zhuǎn)換概率矩陣將對側(cè)向變換性質(zhì)有較強(qiáng)預(yù)測能力,即側(cè)向巖相表現(xiàn)出與垂向巖相相同的疊置規(guī)律。由此,對于空間任意方向(Φ)的轉(zhuǎn)移率rjk,Φ可用下式求得:|rjk,Φ|=√(hxhΦrjk,x)2+(hyhΦrjk,y)2+(hzhΦrjk,z)2(3)式中:hΦ表示任意方向的延遲向量,hΦ=√h2x+h2y+h2z;hx、hy、hz分別表示x、y、z三個(gè)方向的延遲長度;rjk,x表示x方向轉(zhuǎn)移率矩陣的矩陣元素,其它相同;rjk,Φ表示任意方向轉(zhuǎn)移率矩陣的矩陣元素。1.2轉(zhuǎn)移概率最優(yōu)在對地質(zhì)變量進(jìn)行指示化的基礎(chǔ)上,區(qū)內(nèi)各網(wǎng)格的估計(jì)概率由協(xié)同克里格的方法進(jìn)行估計(jì),協(xié)同克里格方法是以全部已知變量(包括已知的原始數(shù)據(jù)及前次已模擬的數(shù)據(jù))為條件(n)的條件概率,下面公式表示的是在x0節(jié)點(diǎn)處的概率估值:Ρ{k(x0)|ij(xα);α=1,2,?,Ν;(4)j=1,2,??Κ}=[ik(x0)]coik=Ν∑α=1Κ∑j=1ij(xα)ωjk,α式中:ωjk,α表示權(quán)系數(shù),可由基于轉(zhuǎn)移概率的協(xié)同克里格方程組求得。對于研究區(qū)D,利用其中的鉆孔資料建立狀態(tài)離散、空間連續(xù)的馬爾柯夫模型T(h),建立協(xié)同克里格方程組:[Τ(x1-x1)?Τ(xΝ-x1)??Τ(x1-xΝ)?Τ(xΝ-xΝ)][W1?WΝ]=[Τ(x0-x1)?Τ(x0-xΝ)]式中:Wα=[ω11,α?ω1Κ,α??ωΚ1,α?ωΚΚ,α],方程組是由地質(zhì)統(tǒng)計(jì)學(xué)中的協(xié)同指示克里格技術(shù)建立的,具體推導(dǎo)過程略。1.3目標(biāo)函數(shù)產(chǎn)生的概率序貫?zāi)M是屬于地質(zhì)統(tǒng)計(jì)學(xué)中的條件模擬,是一種模擬思想,沒有現(xiàn)成的公式。它不但把鄰域內(nèi)所有已知原始數(shù)據(jù)作為后續(xù)模擬的條件數(shù)據(jù),而且還把已模擬實(shí)現(xiàn)的數(shù)據(jù)也作為后續(xù)模擬的條件數(shù)據(jù)。序貫指示模擬開始時(shí)利用了一個(gè)初始的模擬結(jié)構(gòu),即隨機(jī)地給定一個(gè)模擬值(確定的巖相代碼),將利用協(xié)同指示克里格方法所得到的概率估計(jì)值與根據(jù)初始結(jié)構(gòu)所測得的轉(zhuǎn)移概率值相比較,滿足目標(biāo)函數(shù)條件的假設(shè)被接受,排斥不符合目標(biāo)函數(shù)條件的假設(shè)并且重新模擬。目標(biāo)函數(shù)最小化是模擬值是否被采用的判斷標(biāo)準(zhǔn):min{Ο=Μ∑l=1Κ∑j=1Κ∑k=1[tjk(hl)Μeas-tjk(hl)mod]2}(5)式中:O代表一個(gè)目標(biāo)函數(shù);hl代表確定的延遲向量,l=1,2,…,M;tjk(hl)Meas、tjk(hl)Mod分別表示公式計(jì)算的與模擬的轉(zhuǎn)移概率值。2巖相轉(zhuǎn)移概率選取華北平原(面積約14×104km2)作為研究對象。研究區(qū)內(nèi)共收集到827個(gè)鉆孔資料,其中,70%的鉆孔深度在50m以上,最深的鉆孔已揭露第三系,但僅研究第四系沉積規(guī)律,因此不采用第三系巖性,最深模擬深度為400m。通過鉆孔巖性、粒度大小及水理性質(zhì),將區(qū)內(nèi)第四系巖性大致分為4類,建立研究區(qū)內(nèi)的巖相空間。以巖性顆粒粒度的大小為序,各類巖性分別用代號zcs、fxs、ft、fn表示。zcs表示巖性顆粒粒度最大的巖相類別,包括中砂、粗砂、砂礫石等;fxs表示巖性粒度次小,包括砂、中細(xì)砂、細(xì)砂、粉砂;ft包括黃土、粉士、耕作土、粉土與粉砂互層(原始剖面資料中,粉土所占的比例大于粉砂所占的比例)、黃土狀土;fn包括黏土、粉質(zhì)黏土等細(xì)顆粒物質(zhì)。首先,利用定間距(Δh=1m)對鉆孔巖性進(jìn)行采樣,建立垂直方向上空間離散、狀態(tài)離散的轉(zhuǎn)移概率矩陣T(Δh=1),fxs,zcs,ft,fn在整個(gè)所占的比例分別為0.2826,0.1342,0.1566,0.4266。fxszcsftfnΤ(Δh=1)=fxszcsftfn[6.18160.08530.27820.63650.11605.60710.27770.60630.42080.15607.13230.42320.63710.31410.04888.0641]為了便于表述,轉(zhuǎn)移概率矩陣T(Δh)對角線元素值分別表示巖相在z方向的平均長度。采用同樣方法,可得到Δh=2、3、4、…時(shí)的轉(zhuǎn)移概率矩陣及空間連續(xù)狀態(tài)離散的轉(zhuǎn)移概率矩陣,取每個(gè)相應(yīng)轉(zhuǎn)移概率值與對應(yīng)的間隔Δh,得到實(shí)測的轉(zhuǎn)移概率(MeasuredData)曲線(圖1黑色虛線所示)與理論上的馬爾柯夫轉(zhuǎn)移概率(MarkovChain)曲線(圖1淺色曲線所示),橫軸方向代表延遲長度,縱軸方向代表轉(zhuǎn)移概率。研究區(qū)共有4×4=16條擬合曲線,現(xiàn)只選擇了其中的4條曲線,圖1反映了研究區(qū)內(nèi)實(shí)測轉(zhuǎn)移概率曲線與理論轉(zhuǎn)移概率曲線的擬合情況。根據(jù)沃爾索相序定律,由垂向?qū)有蛴?jì)算得到的轉(zhuǎn)換概率矩陣對側(cè)向巖層性質(zhì)有較強(qiáng)的預(yù)測能力,即側(cè)向巖相表現(xiàn)出與垂向巖相相同的疊置規(guī)律,所不同的是,巖相側(cè)向(x、y)的平均長度大于垂向的平均長度。根據(jù)已有實(shí)測剖面,通過測量,合理給定初始巖相側(cè)向長度,然后經(jīng)過模型的多次調(diào)試,確定巖相在水平方面上的平均長度如表1所示。根據(jù)序貫指示模型建立步驟,在研究區(qū)D建立起協(xié)同指示克里格模型,采用平面4km×4km,垂向1m進(jìn)行剖分。模擬時(shí)首先利用一個(gè)初始的模擬結(jié)構(gòu),即隨機(jī)地給定一個(gè)模擬值(確定的巖相代碼),將利用協(xié)同指示克里格方法所得到的概率估計(jì)值與根據(jù)初始結(jié)構(gòu)所測得的轉(zhuǎn)移概率值相比,滿足目標(biāo)函數(shù)條件的假設(shè)被接受,并排斥不符合目標(biāo)函數(shù)條件的假設(shè)并且重新模擬,得到的模擬結(jié)果如圖2、圖3所示。條件模擬結(jié)果顯示,整個(gè)華北平原的沉積結(jié)構(gòu)具有一定的旋回性,自上而下第四系粒度由細(xì)變粗又變細(xì),構(gòu)成一個(gè)較完整的沉積旋回。從圖3由山前到濱海剖面可以看出,華北平原地區(qū)巖性顆粒呈現(xiàn)由粗到細(xì)的總體變化規(guī)律。在山前地帶,絕大部分為顆粒粒徑較大的zcs,到中部平原過渡為細(xì)砂,并且細(xì)顆粒越來越發(fā)育,靠近濱海地區(qū),細(xì)顆粒物質(zhì)含量達(dá)到最大。在濱海平原,大部分巖性由fxs組成,沉積韻律不強(qiáng),模擬結(jié)構(gòu)所表現(xiàn)出的沉積規(guī)律與利用地質(zhì)資料分析得到的華北平原沉積規(guī)律基本一致。模擬結(jié)果與相應(yīng)地段利用鉆孔資料手工繪制的地層剖面相比,兩者有一致的規(guī)律,但不完全一樣。筆者認(rèn)為產(chǎn)生不一樣的原因是鉆孔資料的欠缺、手工繪制的地層剖面存在人為性等。另外,在局部有鉆孔存在的區(qū)域,模擬結(jié)果與鉆孔巖性完全一致。3分區(qū)條件模型法(1)從以上建模過程可以看出,該方法以轉(zhuǎn)移概率矩陣的譜分析函數(shù)來代替?zhèn)鹘y(tǒng)地質(zhì)統(tǒng)計(jì)學(xué)中的變差函數(shù),然后運(yùn)用了比較成熟的協(xié)同克里格估值技術(shù)與序貫?zāi)M的思想建立地質(zhì)結(jié)構(gòu)模型,充分考慮了地層巖性的沉積序列的規(guī)律性。(2)與傳統(tǒng)手工繪制水文地質(zhì)剖面的方法相比,該方法能較客觀地反映出整個(gè)研究區(qū)的沉積規(guī)律,對于研究沉積區(qū)域的地質(zhì)及水文地質(zhì)方面的問題具有較高的利用價(jià)值,特別是對于研究地下水?dāng)?shù)值模擬中儲水介質(zhì)的非均質(zhì)性問題,具有重要的意義。(3)指示克里