馬爾可夫隨機(jī)場(chǎng)反演在火山巖儲(chǔ)層預(yù)測(cè)中的應(yīng)用

曹 磊,張 達(dá),李 寧,荊錫貴

(中國(guó)石油化工股份有限公司東北油氣分公司儲(chǔ)層預(yù)測(cè)團(tuán)隊(duì),吉林長(zhǎng)春130062)

火山巖儲(chǔ)層地震勘探一直是地球物理界關(guān)注的重點(diǎn)之一,伴隨著油氣勘探領(lǐng)域的不斷拓展,火山巖儲(chǔ)層的勘探開(kāi)發(fā)逐漸受到重視。張新榮等[1]分析了火山巖儲(chǔ)層作為復(fù)雜特殊油氣儲(chǔ)層在其勘探開(kāi)發(fā)中面臨的諸多難題;宋宗平[2]和張芝銘等[3]應(yīng)用基于寬帶約束的模擬退火反演方法并結(jié)合地震數(shù)據(jù)、測(cè)井?dāng)?shù)據(jù)等獲取火山巖縱波速度、橫波速度、密度等彈性參數(shù),從而為地質(zhì)研究人員提供重要的解釋依據(jù),它已成為火山巖儲(chǔ)層勘探開(kāi)發(fā)的重要研究?jī)?nèi)容;陳永波等[4]通過(guò)古地貌恢復(fù)研究火山巖沉積;朱超等[5]通過(guò)疊后波阻抗反演預(yù)測(cè)火山巖優(yōu)質(zhì)儲(chǔ)層平面展布特征;楊曉光等[6]采用分頻屬性反演和基于神經(jīng)網(wǎng)絡(luò)的儲(chǔ)層參數(shù)反演有效預(yù)測(cè)了火山巖儲(chǔ)層分布;毛慶鑫[7]認(rèn)為火山巖結(jié)構(gòu)復(fù)雜、巖性橫向變化快、具有較強(qiáng)的各向異性特征,在地震剖面上表現(xiàn)為同相軸橫向連續(xù)性差、反射呈現(xiàn)低頻和強(qiáng)振幅特征,可見(jiàn)空白反射或雜亂反射結(jié)構(gòu)。

然而,基于疊后地震數(shù)據(jù)的地震屬性(相干、波形聚類(lèi)等)分析和疊后波阻抗反演等技術(shù)手段對(duì)火山巖機(jī)構(gòu)內(nèi)部期次劃分、巖性和物性的預(yù)測(cè)存在不足,越來(lái)越不能夠滿足火山巖勘探開(kāi)發(fā)的需求[8]。針對(duì)這一問(wèn)題,何輝等[9]應(yīng)用疊前方位各向異性法,通過(guò)優(yōu)選衰減起始頻率屬性,綜合預(yù)測(cè)了佳木河組火山巖儲(chǔ)層裂縫分布特征;張新培等[10]利用測(cè)井響應(yīng)和多井約束下的波阻抗反演方法實(shí)現(xiàn)了優(yōu)質(zhì)火山巖儲(chǔ)層的識(shí)別和預(yù)測(cè);印興耀等[11]認(rèn)為,目前應(yīng)用縱波速度、橫波速度和密度參數(shù)等結(jié)合井參與的低頻背景模型進(jìn)行疊前同時(shí)反演來(lái)識(shí)別火山巖機(jī)構(gòu)內(nèi)部?jī)?chǔ)層效果不佳的原因是,地質(zhì)模型與實(shí)際地質(zhì)情況差異較大而導(dǎo)致反演結(jié)果準(zhǔn)確率不高。如何弱化低頻背景模型的影響,使反演結(jié)果既符合地質(zhì)規(guī)律又有較高的準(zhǔn)確度是目前火山巖機(jī)構(gòu)內(nèi)部?jī)?chǔ)層預(yù)測(cè)的關(guān)鍵所在。

馬爾可夫隨機(jī)場(chǎng)反演方法在儲(chǔ)層建模領(lǐng)域早有應(yīng)用,TARANTOLA[12]最早提出貝葉斯框架下的概率反演方法,基于馬爾可夫鏈蒙特卡洛(markov chain monte carlo,MCMC)統(tǒng)計(jì)方法獲得待反演參數(shù)的后驗(yàn)概率密度函數(shù),統(tǒng)計(jì)分析反演結(jié)果,獲取最優(yōu)的反演解,并對(duì)反演結(jié)果進(jìn)行不確定性分析;SPIKES[13]將多種蒙特卡洛抽樣方法引入地質(zhì)統(tǒng)計(jì)建模,聯(lián)合估算阻抗等彈性參數(shù)和孔隙度等儲(chǔ)層屬性;MARIT等[14]使用馬爾可夫隨機(jī)場(chǎng)建立巖性/流體分類(lèi)概率分布的先驗(yàn)?zāi)Ｐ?基于MCMC方法獲取彈性參數(shù)、孔隙度以及巖性/流體分類(lèi)后驗(yàn)概率分布的預(yù)測(cè)和模型參數(shù)的不確定性評(píng)估;MICHEL等[15]將分級(jí)貝葉斯模型引入MCMC方法;ODD等[16]使用MCMC隨機(jī)反演方法反演縱波阻抗;陳裕華[17]應(yīng)用縱波、轉(zhuǎn)換波振幅和旅行時(shí)數(shù)據(jù),反演孔隙度等儲(chǔ)層物性參數(shù);田玉昆[18]采用Huber-Markov隨機(jī)場(chǎng)約束邊緣保護(hù)方法,求解疊前縱橫波速度和密度三參數(shù)反演問(wèn)題?？梢?jiàn),在學(xué)者們的不斷深入研究下,馬爾可夫隨機(jī)場(chǎng)反演方法在最優(yōu)解、去除邊界效應(yīng)等方面不斷完善,但大部分的研究重點(diǎn)是改進(jìn)算法和模型試算上,只有少部分研究應(yīng)用到實(shí)際數(shù)據(jù),而且是碎屑巖和碳酸鹽巖數(shù)據(jù)中。

針對(duì)火山巖地層內(nèi)部地震反射波橫向變化快、呈現(xiàn)中弱反射和非均質(zhì)性強(qiáng)的特點(diǎn),目前仍沒(méi)有較好的反演方法。馮玉輝等[19]系統(tǒng)地進(jìn)行了火山巖地震響應(yīng)特征描述;代春萌等[20]進(jìn)行了火山巖相控反演;有研究者應(yīng)用體序域建模識(shí)別火山巖儲(chǔ)層,對(duì)于沒(méi)有井鉆遇的火山巖反演結(jié)果與體序域模型具有較強(qiáng)的相似性。

本文將馬爾可夫隨機(jī)場(chǎng)引入到火山巖地震疊前同時(shí)反演算法中,針對(duì)火山巖巖性多樣且測(cè)井曲線數(shù)值疊合嚴(yán)重的特點(diǎn),繞開(kāi)火山巖巖性的預(yù)測(cè),利用測(cè)井?dāng)?shù)據(jù)劃分巖相,對(duì)定義的每個(gè)巖相都計(jì)算出低頻背景模型,并將每個(gè)巖相應(yīng)用馬爾可夫隨機(jī)場(chǎng)求取先驗(yàn)分布函數(shù),從而有效避免了反演結(jié)果由于地震能量差異和頻率不足等引起的假象,提高了預(yù)測(cè)精度。

1 疊前同時(shí)反演面臨的問(wèn)題

1.1 歸一化問(wèn)題

在應(yīng)用疊前確定性同時(shí)反演方法進(jìn)行巖石物理交會(huì)分析時(shí),需要確定縱波阻抗和橫波阻抗的線性關(guān)系。圖1為典型的橫波阻抗與縱波阻抗交會(huì)圖。分析圖1可知,圖中線性關(guān)系需要3條線才能夠表述清楚,兩條通過(guò)傾向砂巖的藍(lán)色數(shù)據(jù)團(tuán)(含水和含烴砂巖)以及一條通過(guò)傾向泥巖的橙色數(shù)據(jù)團(tuán)。

1.2 低頻模型問(wèn)題

目前,建立低頻模型通常是將井的波阻抗曲線進(jìn)行低通濾波和內(nèi)插,這種方法僅適用于較平緩的地層,對(duì)于復(fù)雜的地下地質(zhì)情況該方法主要存在以下兩個(gè)方面的不足:

1) 對(duì)于高凈毛比值(有效儲(chǔ)層厚度與儲(chǔ)層總厚度的比值)的井和低凈毛比值的井,由于無(wú)法確定某個(gè)位置的準(zhǔn)確凈毛比,因此期望在兩者之間內(nèi)插獲得可靠的低頻背景模型不現(xiàn)實(shí);

圖1 橫波阻抗和縱波阻抗交會(huì)分析結(jié)果

2) 低頻模型的實(shí)現(xiàn)算法只包含針對(duì)平行于頂層、平行于底層、頂?shù)讓又g等比例模式內(nèi)插的3種算法,適用于平行地層、超覆和尖滅地層,但是對(duì)于碳酸鹽縫洞、火山巖內(nèi)幕、復(fù)雜斷塊等復(fù)雜地層問(wèn)題適用性較差,反演結(jié)果會(huì)出現(xiàn)較大誤差。

1.3 物理特性識(shí)別

圖2給出了縱波阻抗模型及其合成地震記錄與疊前同時(shí)反演的縱波阻抗結(jié)果。圖2中黑色曲線代表縱波阻抗真實(shí)模型(橫波阻抗和密度也是一樣,未顯示),紅色曲線表示縱波阻抗疊前同時(shí)反演結(jié)果,黑色虛線是計(jì)算縱波阻抗用到的低頻背景模型,右側(cè)是零入射角的合成地震記錄。分析圖2可知:合成地震記錄中同相軸最大振幅值位置對(duì)應(yīng)的是巖相變換的地方,疊前同時(shí)反演結(jié)果只能在這些巖相變換位置附近提供可靠的結(jié)果,遠(yuǎn)離這些巖相變換的位置,縱波阻抗結(jié)果趨向于低頻模型(如箭頭所示)。

圖2 縱波阻抗模型及其合成地震記錄和疊前同時(shí)反演的縱波阻抗結(jié)果

2 方法原理

在理想情況下,希望用一種不依賴于低頻背景模型的內(nèi)插方法來(lái)得到阻抗信息。采用迭代方法,先用地震數(shù)據(jù)反演出相(給定一個(gè)均勻的阻抗模型),然后由反演的相得到地震數(shù)據(jù)的阻抗信息,再在這些阻抗信息基礎(chǔ)上,由地震數(shù)據(jù)再反演出相,迭代至收斂。為了實(shí)現(xiàn)這一反演策略,引入馬爾可夫隨機(jī)場(chǎng)反演理論。

2.1 馬爾可夫隨機(jī)場(chǎng)原理

馬爾可夫隨機(jī)場(chǎng)(MRF)方法建立在馬爾可夫隨機(jī)模型(MR)和Bayes理論的基礎(chǔ)上,馬爾可夫隨機(jī)場(chǎng)模型提供了不確定性描述與先驗(yàn)知識(shí)聯(lián)系的紐帶,并利用觀測(cè)圖像根據(jù)統(tǒng)計(jì)決策和估計(jì)理論中的最優(yōu)準(zhǔn)則,確定分割問(wèn)題的目標(biāo)函數(shù),求解滿足這些條件或目標(biāo)函數(shù)的最大可能分布,從而將分割問(wèn)題轉(zhuǎn)化為最優(yōu)化問(wèn)題[21-22]。

MRF是基于圖像中像素的灰度信息,利用像素之間的互相作用與空間的關(guān)系,實(shí)現(xiàn)圖像分割[23]。MRF模型由于自身的特點(diǎn),既能反映圖像的隨機(jī)性,又能反映圖像的空間結(jié)構(gòu),從而能夠有效描述特定圖像的空間性質(zhì),保證了目標(biāo)函數(shù)有唯一的最小點(diǎn);在求解MRF描述的不確定性問(wèn)題時(shí),將圖像先驗(yàn)知識(shí)轉(zhuǎn)化為先驗(yàn)概率分布,基于貝葉斯理論,采用先驗(yàn)分布函數(shù)和似然函數(shù)來(lái)描述最大后驗(yàn)估計(jì)的圖像分布,同時(shí)保留了邊緣結(jié)構(gòu),便于邊緣識(shí)別。應(yīng)用數(shù)學(xué)公式推導(dǎo)出的參數(shù)具有明確的物理含義,而常規(guī)線性模型中的參數(shù)一般只作為被擬合了的參數(shù)而出現(xiàn),較少具有真實(shí)含義。

2.2 馬爾可夫隨機(jī)場(chǎng)的地震反演

將地震反演看作一個(gè)統(tǒng)計(jì)問(wèn)題,貝葉斯定理可以寫(xiě)成:

π(Z|Sreal)≈L(Sreal|Z)p(Z)

(1)

式中:Z是波阻抗;Sreal是實(shí)際地震數(shù)據(jù);π(Z|Sreal)是后驗(yàn)分布函數(shù);L(Sreal|Z)是似然函數(shù);p(Z)是先驗(yàn)分布。

將分布表示為多正態(tài)形式,通常比較合理,并更容易進(jìn)行數(shù)學(xué)處理。先驗(yàn)分布P(Z)可以從縱波阻抗的線性深度趨勢(shì)、縱波阻抗與橫波阻抗、縱波阻抗與密度之間的交會(huì)結(jié)果得到,并通過(guò)井?dāng)?shù)據(jù)進(jìn)行不確定性評(píng)估。這些趨勢(shì)擬合可以表示為一個(gè)多正態(tài)先驗(yàn)分布P(Z),形式為:

(2)

式中:CP是描述阻抗方差和相關(guān)性的協(xié)方差矩陣;Z0是波阻抗的初始模型。

似然函數(shù)L(Sreal|Z)可以表示為:

L(Sreal|Z)?

(3)

式中:F(Z)是從阻抗Z導(dǎo)出合成地震記錄的函數(shù);Cd是表示地震噪聲的協(xié)方差矩陣。

理想情況下,為了準(zhǔn)確描述地震反演問(wèn)題的物理特性,應(yīng)對(duì)每個(gè)巖相和阻抗進(jìn)行反演,并且每個(gè)相的巖石物理模型和低頻模型都應(yīng)該使用,這樣做在數(shù)學(xué)上的要求很高,同時(shí)也不能用標(biāo)準(zhǔn)的優(yōu)化公式來(lái)解決,但可以使用迭代方法解決這一問(wèn)題:

首先給定一個(gè)均勻的阻抗模型從地震數(shù)據(jù)中反演出相,然后從這些相中反演出阻抗信息,再在阻抗信息基礎(chǔ)上從地震數(shù)據(jù)中再反演出相,迭代至收斂[18-19]。在這種情況下,公式(1)可以寫(xiě)成:

π(Z,F|Sreal)≈L(Sreal|Z)p(Z|F)p(F)

(4)

式中:π(Z,F|Sreal)是后驗(yàn)分布函數(shù);p(Z|F)為巖相的先驗(yàn)分布;p(F)為巖相的概率密度。

先驗(yàn)分布p(Z|F)與(2)式非常相似,唯一的區(qū)別在于先驗(yàn)均值Z0和協(xié)方差矩陣CP取決于:

p(Z|F)=

(5)

對(duì)于巖相的概率密度p(F),應(yīng)用離散馬爾可夫隨機(jī)場(chǎng)求取[24-25]。簡(jiǎn)單地說(shuō),任何三維地震網(wǎng)格都是由許多像素組成,定義一系列像素對(duì),在反演時(shí)窗內(nèi)的每個(gè)像素能直接與相鄰像素相連接,因此每個(gè)像素屬于6個(gè)像素對(duì)(邊緣、角落除外)。然后,整個(gè)網(wǎng)格的配置F的概率由所有像素對(duì)上的勢(shì)能之和定義,即所謂的吉布斯分布:

P(F)≈exp[-∑Vp(Fp)]

(6)

其中,Vp表示每組像素對(duì)上的相系列Fp的“勢(shì)能”。

對(duì)于每個(gè)像素團(tuán),可以將Vp寫(xiě)為Vp=β·I(Fcentre,Fneighbour),其中,Fcentre為中心相,Fneighbour為臨近相。如果Fcentre和Fneighbour相同,則離散指標(biāo)函數(shù)I為0,如果它們不相同則為1,并且β為一個(gè)正連續(xù)參數(shù)。所以,公式(6)可改寫(xiě)為:

P(F)≈exp[-∑βI(Fcentre,Fneighbour]

(7)

因此,只要兩個(gè)相鄰相不同,就會(huì)減少概率P(F)。由于水平方向(主測(cè)線和聯(lián)絡(luò)線方向)和垂直方向的連續(xù)性不同,因此,需要對(duì)其賦予不同的β值。然而,由于地質(zhì)體很少是完全水平的,使用年代地層的年齡概念來(lái)確定的相同年齡的鄰點(diǎn),有可能在時(shí)間坐標(biāo)上并不是鄰近的像素。馬爾可夫隨機(jī)場(chǎng)可以阻止不合理的巖相組合,即可以確保不會(huì)在含氣砂巖的頂部得到含水砂巖。與常規(guī)疊前同時(shí)反演方法的不同之處在于增加了每個(gè)相的深度趨勢(shì)和交會(huì)圖進(jìn)行約束,而沒(méi)有使用變差函數(shù)。

3 模型數(shù)據(jù)分析

基于馬爾可夫隨機(jī)場(chǎng)的反演方法是貝葉斯同時(shí)反演方法的一種擴(kuò)展。兩種方法之間的最大區(qū)別在于:在貝葉斯形式的同時(shí)反演中,縱波阻抗與深度是一個(gè)線性趨勢(shì);而在馬爾可夫隨機(jī)場(chǎng)的反演中可以反演兩個(gè)或者多個(gè)巖相,所以可以指定兩個(gè)或更多這樣的趨勢(shì)。這意味著在反演時(shí)窗內(nèi),阻抗結(jié)果與每個(gè)相趨勢(shì)之間都可以建立聯(lián)系。

為了驗(yàn)證基于馬爾可夫隨機(jī)場(chǎng)反演算法的有效性,建立了一個(gè)如圖3a所示砂巖侵入泥巖的楔狀模型。模型頻寬0～60Hz,并應(yīng)用25Hz雷克子波進(jìn)行反演,為了更好地進(jìn)行對(duì)比分析,只考慮縱波速度的變化。圖3b給出了當(dāng)入射角為0時(shí)的地震正演結(jié)果。對(duì)正演結(jié)果進(jìn)行常規(guī)疊前同時(shí)反演和基于馬爾可夫隨機(jī)場(chǎng)反演,兩種反演結(jié)果分別如圖3c和圖3d所示。分析可知:在常規(guī)疊前同時(shí)反演中使用了一個(gè)低頻背景模型,其縱波速度值在泥巖縱波速度值和砂巖縱波速度值之間;基于馬爾可夫隨機(jī)場(chǎng)的反演使用了兩個(gè)低頻趨勢(shì),其縱波速度值分別為泥巖縱波速度值和砂巖縱波速度值,從反演結(jié)果對(duì)比來(lái)看,基于馬爾可夫隨機(jī)場(chǎng)的反演結(jié)果與模型更加吻合,預(yù)測(cè)準(zhǔn)確率更高。

圖3 模型反演分析

4 實(shí)例分析

將同時(shí)反演方法和基于馬爾可夫隨機(jī)場(chǎng)的反演方法應(yīng)用于長(zhǎng)嶺龍鳳山地區(qū)中基性火山巖油氣儲(chǔ)層的預(yù)測(cè)。由測(cè)井?dāng)?shù)據(jù)分析可知,研究區(qū)目的層系儲(chǔ)層以凝灰?guī)r為主,因此定義了3個(gè)相來(lái)進(jìn)行反演:將孔隙度大于5%且泥質(zhì)含量小于40%的凝灰?guī)r對(duì)應(yīng)的波阻抗定義為凝灰?guī)r高孔隙相;將孔隙度小于5%且泥質(zhì)含量小于40%的凝灰?guī)r對(duì)應(yīng)的波阻抗定義為凝灰?guī)r干層相;將泥質(zhì)含量大于40%的巖性對(duì)應(yīng)的波阻抗定義為泥巖相。

由密度概率分布圖(圖4)可以看出,密度能夠?qū)?chǔ)層進(jìn)行有效識(shí)別,因此選取密度參數(shù)作為敏感參數(shù)參與到后續(xù)反演之中。

圖4 密度概率分布

研究區(qū)目的層段反演結(jié)果如圖5所示。其中,圖5a 為常規(guī)疊前同時(shí)反演得到的密度數(shù)據(jù)體;圖5b為基于馬爾可夫隨機(jī)場(chǎng)反演得到的巖相體與地震剖面疊合顯示結(jié)果,圖5b中紅色代表凝灰?guī)r高孔隙相,黃色代表凝灰?guī)r干層相,灰色代表泥巖相;圖5c為基于馬爾可夫隨機(jī)場(chǎng)反演得到的密度數(shù)據(jù)體。圖5a中反演結(jié)果在井眼處與密度曲線對(duì)應(yīng)較好,但同一火山巖機(jī)構(gòu)之間的兩口井橫向連續(xù)性和整體期次差異較大;由于劃分了3個(gè)巖相,進(jìn)行了3個(gè)低頻模型約束,圖5b 中反演的巖相體在橫向展布上與地震資料更加吻合,在縱向油氣水配置關(guān)系上更符合石油地質(zhì)理論和實(shí)際地質(zhì)情況;圖5c中的結(jié)果不僅在井眼處與密度曲線對(duì)應(yīng)較好,而且同一火山巖機(jī)構(gòu)之間的兩口井橫向連續(xù)性更好,能夠較好地識(shí)別火山巖期次,同時(shí)也能夠?qū)?chǔ)層物性進(jìn)行較好的預(yù)測(cè),可指導(dǎo)下一步井位部署。本文方法在松遼盆地龍鳳山地區(qū)火山巖的勘探開(kāi)發(fā)中發(fā)揮了重要作用,依據(jù)預(yù)測(cè)結(jié)果部署了3口直井均獲得工業(yè)氣流,4口水平井均獲得了高產(chǎn)氣流,以預(yù)測(cè)結(jié)果為依據(jù)在該地區(qū)提交火山巖天然氣探明儲(chǔ)量102×108m3。

圖5 反演剖面對(duì)比

5 結(jié)論

本文將馬爾可夫隨機(jī)場(chǎng)引入到火山巖地震疊前反演算法中,其優(yōu)點(diǎn)在于:采用火山巖多條測(cè)井?dāng)?shù)據(jù)交會(huì)劃分巖相,且測(cè)井?dāng)?shù)據(jù)既不直接參與低頻模型的建立也不直接參與反演,該方法對(duì)定義的每個(gè)巖相都計(jì)算出低頻背景趨勢(shì),使反演結(jié)果有效避免了由于地震能量差異、頻率不足等引起的假象。

從馬爾可夫隨機(jī)場(chǎng)反演結(jié)果來(lái)看:

1) 馬爾可夫隨機(jī)場(chǎng)反演方法不依賴傳統(tǒng)的低頻模型,反演結(jié)果更遵從于地震數(shù)據(jù);

2) 反演得到的巖相體符合石油地質(zhì)學(xué)規(guī)律,對(duì)地質(zhì)研究具有很大幫助;

3) 反演得到的密度體與井震標(biāo)定后的測(cè)井曲線有較好的匹配關(guān)系,可以用于該地區(qū)火山巖機(jī)構(gòu)內(nèi)部?jī)?chǔ)層的物性預(yù)測(cè)。

該反演方法適合于火山巖機(jī)構(gòu)內(nèi)部?jī)?chǔ)層物性及含氣性預(yù)測(cè),對(duì)其它類(lèi)似火山巖儲(chǔ)層的地區(qū)具有推廣意義。