曹 磊,張 達,李 寧,荊錫貴

(中國石油化工股份有限公司東北油氣分公司儲層預測團隊,吉林長春130062)

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

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

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

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

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

1 疊前同時反演面臨的問題

1.1 歸一化問題

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

1.2 低頻模型問題

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

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

圖1 橫波阻抗和縱波阻抗交會分析結果

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

1.3 物理特性識別

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

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

2 方法原理

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

2.1 馬爾可夫隨機場原理

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

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

2.2 馬爾可夫隨機場的地震反演

將地震反演看作一個統(tǒng)計問題,貝葉斯定理可以寫成:

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

(1)

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

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

(2)

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

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

L(Sreal|Z)?

(3)

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

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

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

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

(4)

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

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

p(Z|F)=

(5)

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

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

(6)

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

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

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

(7)

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

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

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

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

圖3 模型反演分析

4 實例分析

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

由密度概率分布圖(圖4)可以看出,密度能夠對儲層進行有效識別,因此選取密度參數(shù)作為敏感參數(shù)參與到后續(xù)反演之中。

圖4 密度概率分布

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

圖5 反演剖面對比

5 結論

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

從馬爾可夫隨機場反演結果來看:

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

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

3) 反演得到的密度體與井震標定后的測井曲線有較好的匹配關系,可以用于該地區(qū)火山巖機構內部儲層的物性預測。

該反演方法適合于火山巖機構內部儲層物性及含氣性預測,對其它類似火山巖儲層的地區(qū)具有推廣意義。