基于構造約束聯合概率反演的油藏參數表征方法

張 健 李景葉 王建花 陳小宏 李遠強 周春雷

(①中國石油大學(北京)海洋石油勘探國家工程實驗室，北京 102249；②中國石油大學(北京)油氣資源與探測國家重點實驗室，北京 102249；③西南交通大學地球科學與環境工程學院，四川成都 611756；④中海油研究總院有限責任公司，北京 100028；⑤中國石油勘探開發研究院西北分院，甘肅蘭州 730020)

0 引言

彈性參數反演、巖石物理反演和儲層巖相識別是油藏表征和評價的有效手段[1-3]。測井資料可提供井位附近油藏參數的直接或間接測量結果，而地震數據提供了獲取井間油藏參數的唯一指導信息。但由于地震頻帶、觀測誤差以及巖石物理建模誤差等因素的影響，地震反問題的求解仍然面臨較大挑戰。傳統方法通常利用井數據低頻信息和單道地震數據通過多步驟反演策略依次得到儲層彈性參數、儲層物性參數和儲層巖相等信息，忽略了井數據高頻信息和地震數據橫向信息以及各反演環節不確定性信息對反演結果的影響，降低了最終反演精度，增加了油藏評價和開發的風險。因此，如何充分利用現有資料獲取客觀、準確表征反演結果不確定性的信息，從而指導儲層評價尤為重要。

對于油藏參數預測及其不確定性評價，前人做了大量的研究工作?；诟怕史椒ǖ呢惾~斯反演框架通過引入先驗信息約束項，有效降低了地球物理反演的多解性，提高了反演精度，被廣泛用于求解地球物理反問題。貝葉斯反演的最終解由后驗概率分布定義，表示在給定觀測數據條件下待求模型參數正確的可能性。Buland等[4]基于貝葉斯反演框架，結合線性AVO方程和高斯分布假設，推導了彈性參數后驗概率分布解析表達式，在獲得彈性參數的同時，定量表征了反演結果的不確定性。Grana等[5]在Buland等[4]的研究基礎上，將貝葉斯線性反演的高斯假設推廣到混合高斯假設，得到了彈性參數的解析解。為提高反演結果的橫向連續性，Hamid等[6]基于貝葉斯反演框架，引入橫向約束項開展多道疊后阻抗反演。張廣智等[7]基于Fatti波阻抗近似式，利用蒙特卡洛—馬爾科夫鏈(MCMC)算法預測了彈性三參數，并提供了相關的不確定性信息。Aleardi等[8]針對某一特定研究區，通過對比線性反演方法與非線性反演方法的結果，驗證了前者求取儲層參數及其不確定性的有效性。黃捍東等[9]利用非線性隨機反演方法預測了陸相薄砂巖儲層。Larsen等[10]闡述了地震反演不確定性對儲層巖相識別的影響，引入馬爾科夫鏈先驗信息提高了巖相識別精度，并實現了巖相隨機模擬。李軍等[11]將馬爾科夫鏈模型從二維形式推廣到三維，獲得了更準確的儲層巖相模擬結果。劉興業等[12]基于貝葉斯框架，引入核貝葉斯算法預測了物性參數、巖相及其不確定性信息。袁成等[13-14]系統介紹了由地震數據識別儲層巖相各環節的不確定性評價方法。然而，上述方法均為多步驟反演，很難考慮各個環節的不確定性。Bosch等[15]以地震數據作為輸入，利用MCMC方法同時預測波阻抗和孔隙度，并定量分析了反演結果的不確定性。de Figueiredo等[16]提出基于巖石物理的貝葉斯反演方法，同時估計了波阻抗、孔隙度和儲層巖相及其不確定性信息。Grana[17]提出了一種基于模型參數為混合非參數分布統計假設的貝葉斯反演方法，定量表征了儲層巖相和儲層流體。李海山等[18]利用平面波分解技術提取地層信息，建立更可靠的波阻抗初始模型用于反演，獲得了橫向更連續的反演結果。曹丹平等[19]、周曉越等[20]嘗試多屬性聯合反演。

基于前人的工作，本文提出基于構造約束聯合概率反演的油藏參數表征方法?；谪惾~斯反演框架提出油藏參數同時反演策略，利用最小二乘井數據插值將地質構造信息和井信息整合到同時反演框架，經推導得到油藏參數后驗概率分布解析表達式。最終，利用獲得的后驗概率分布信息定量評價和表征反演結果。所提方法通過同時反演策略客觀表征各環節的不確定性，構造約束保證反演結果橫向空間耦合性，最終反演結果為油藏表征和評價提供了更可靠的參考信息。實際資料應用和井測試驗證了方法的有效性。

1 方法理論

為了研究基于巖石物理分析統計的油藏參數關系，引入構造約束地質模型，構建油藏參數同時反演解析表達式，獲取構造約束的油藏參數(包括波阻抗、孔隙度和儲層巖相)聯合概率分布信息，并開展定量評價。理論方法包括巖石物理建模及分析、構造約束和井信息融合以及油藏參數聯合概率表征。

1.1 線性巖石物理建模及分析

巖石物理模型(rock physics model，RPM)是連接儲層物性參數和儲層彈性參數的橋梁，是聯合反演的基礎，也是參數敏感性分析以及參數優選的重要工具。文中主要針對疊后數據，其RPM(即波阻抗和儲層物性參數之間關系)為

Zp=fRPM(R)+ε

(1)

式中：fRPM(R)為線性或理論RPM，R為儲層物性參數；Zp為聲阻抗；ε為觀測誤差，代表建模結果與實際觀測值之間的失配程度，通過分析井數據確定。

為得到油藏參數聯合后驗概率解析表達式，需建立儲層物性參數與聲阻抗自然對數線性RPM，即

lnZp=a1+a2D+a3φ+a4Sw+a5Vsh

(2)

式中：D為深度；φ、Sw、Vsh分別為孔隙度、含水飽和度和泥質含量；a1～a5為利用井數據或者巖心數據得到的回歸系數。

對M區井數據(圖1)利用多元線性回歸方法得到a1～a5，則由式(2)得

lnZp=9.0542+8×10-5D-2.4415φ+

0.0092Sw-0.5464Vsh

(3)

式(3)表明：孔隙度在聲阻抗預測中起決定作用，其次是泥質含量和含水飽和度；由于有限的深度范圍，深度對阻抗預測的影響非常小。

結合M區地質與地球物理特征以及測井阻抗和物性參數之間的響應關系，選取基于Hertz-Mindlin接觸理論的RPM(表1)建立阻抗與物性參數的關系(圖2)，以證明線性RPM在M區的可行性[21]。Hertz-Mindlin模型已被證明可以較好地表征砂巖物性參數與彈性參數之間的關系。

表1 Hertz-Mindlin理論RPM參數

由圖1、圖2可見：兩種RPM阻抗預測結果與實測阻抗吻合，證明了RPM的有效性(圖1)；兩種RPM參數敏感性分析結果一致，證明線性RPM可以很好表征阻抗與物性參數的關系，且敏感性分析結論與式(3)一致，即孔隙度對阻抗預測起決定性作用，另外兩種物性參數對于阻抗預測作用較小(圖2)。因此，本文重點研究阻抗與孔隙度之間的關系。

圖1 M區井數據(a)φ；(b)Vsh；(c)Sw；(d)lnZp；(e)巖相

圖2 Hertz-Mindlin理論RPM(a)與多元線性回歸得到的線性RPM(b)的Zp與φ(左)、Vsh(中)、Sw(右)敏感性分析圖

將井數據φ—lnZp聯合向量記為m，則地下某一時間點(t時刻)的阻抗和孔隙度聯合分布(圖3a的彩色底圖，服從多元高斯分布)為

圖3 孔隙度與阻抗對數交會圖(a)lnZp和φ聯合分布(彩色底圖)、線性RPM(黑色曲線)及φ—lnZp交會圖(綠色散點)疊合圖；(b)不同巖相的φ—lnZp交會圖(散點)及其先驗分布

p(mt)=p[ln(Zpt),φt]～N(mt;μ,Σ)

(4)

式中：下標t表示t時刻的相應變量(下同)；N表示均值為μ、協方差為Σ的多元高斯分布

(5)

其中Σ對角線上的元素為各個隨機變量的方差，非對角線上的元素為兩兩隨機變量之間的協方差。

將式(4)擴展到整個地震道，同時引入垂向空間耦合信息，有

p(m)=p[ln(Zp),φ]～N(m;μm,Σm)

(6)

其中

(7)

式中Γ代表垂直相關矩陣，可定義為[4]

(8)

式中：τ為時間變量矩陣；l為垂直方向耦合系數，通過統計井數據獲得，本研究中l=9。

1.2 構造約束和井信息融合

橫向連續性和分辨率是儲層參數表征的重要要求，而地質構造信息和井數據可分別提供橫向空間耦合關系和高頻細節信息。為進一步提高反演結果的橫向連續性和分辨率，建立儲層參數m和構造約束井插值結果f的關系，即

f=m+ef

(9)

式中ef代表誤差容忍參數，即m與f之間的相似性，可以利用距井距離的函數表示。

構造約束井插值過程可表示為[22-23]

Mf=w

(10)

式中：M為采樣算子；w為已知井數據。對式(10)利用最小二乘方法求解，有

(11)

(12)

1.3 油藏參數聯合概率表征

高精度儲層參數表征及其不確定性評價是儲層表征的另一重要需求。與傳統多步反演方法相比，同時反演方法可以更準確地表征儲層參數并量化其不確定性。為實現同時反演，聯合觀測地震數據和構造約束井插值結果建立如下方程

(13)

式中：d為地震響應；ed為觀測誤差；G=[WD0]為地震響應正演算子，其中W為子波矩陣，而為一階差分矩陣。

p(m|y)～N(m;μm|y,Σm|y)

(14)

其中

(15)

(16)

而

(17)

(18)

在某一巖相(對應m)服從高斯分布(圖3b)假設條件下，有

p(mt|κ)～N(mt;μκ,Σκ)

(19)

式中μκ和Σκ分別為巖相κ對應的儲層參數m的均值和協方差。

根據鏈式法則，得到

(20)

其中

p(κt|mt)∝p(mt|κt)p(κt)

(21)

p(mt|κt)、p(mt|yt)均服從高斯分布，有p(mt|κt)～N(m;μ1,Σ1)，p(mt|yt)～N(m;μ2,Σ2)。其中，μi、Σi(i=1,2)分別表示高斯分布中的均值和協方差矩陣。結合式(20)和式(21)，根據概率統計法則[24-25]，巖相后驗分布為(附錄A)

(22)

其中

(23)

而

(24)

(25)

假定式(22)積分項為1，同時引入馬爾科夫鏈先驗，則式(22)改寫為

p(κt|yt)∝∏p(κt|κt-1)cc

(26)

其中p(κ1|κ0)=p(κ1)表示巖相先驗分布，p(κt|κt-1)可利用馬爾科夫鏈轉移概率矩陣(通過統計井數據獲得)計算。

最后，可以得到油藏參數(阻抗、孔隙度、巖相)預測結果及其不確定性評價信息。

1.4 算法實現流程

(1)數據準備，包括井數據和地震數據；

(2)利用平面波分解濾波技術由地震數據提取地層傾角信息；

(3)結合提取的地層傾角信息和井數據建立構造約束與井融合的地下參數模型；

(4)統計、分析井數據，獲取巖相依賴阻抗和孔隙度聯合先驗分布、垂向耦合系數以及轉移概率矩陣等相關信息；

(5)最后，計算式(15)、式(16)、式(26)，分別得到阻抗、孔隙度、巖相及其后驗概率分布。

2 應用實例

為了驗證所提方法的有效性，利用M區資料進行試算。M區存在濁積砂巖油藏，在儲層段有含油砂巖、含水砂巖和泥巖等三種巖相。圖4為疊后地震剖面。利用3口井(A、B、D)數據直接參與反演，而另1口井(C)作為盲井驗證反演結果的可靠性。首先利用平面波分解濾波技術獲取地層傾角數據(圖5)，結合地層傾角信息和3口井數據，利用式(12)得到構造約束井插值結果(圖6)。

圖4 疊后地震剖面黑色曲線之間代表主要油藏目的層；A，B，C，D為井位

圖5 地層傾角

圖6 構造約束井插值結果(a)阻抗；(b)孔隙度；(c)阻抗方差；(d)孔隙度方差方差代表井數據對最終反演結果的影響，是距井的距離的函數，在井附近井約束權重較大，遠離井約束權重逐漸減小

將觀測地震數據和構造約束井插值結果代入式(13)，利用式(15)和式(16)得到阻抗及孔隙度反演結果(圖7)。可見：①本文方法反演結果(圖7a、圖7b)的橫向連續性較無約束反演結果(圖7e、圖7f)更好；②除約束井(A，B，D)附近以外，兩種方法獲得的阻抗方差(圖7c、圖7g)相近；③由于本文方法同時考慮了阻抗反演與孔隙度反演之間的耦合關系，除約束井附近以外，孔隙度方差(圖7d)大于無約束反演結果(圖7h)，更好地反映了反演結果的不確定性；④井數據既能提供低頻信息，又能提供高分辨率信息，本文方法同時引入井約束，在井附近可獲得更高分辨率反演結果，對應概率在井位處的方差接近0，即反演結果與井數據一致(圖7c、圖7d)。

將得到的巖相依賴孔隙度與阻抗聯合高斯分布均值和方差(圖3b)、阻抗與孔隙度及其方差(圖7)代入式(26)，可以得到M區儲層巖相概率，其中轉移概率矩陣T由統計已知測井數據得到圖8為巖相預測結果及其不確定性。由圖可見：本文方法識別的巖相(圖8a)較傳統多步法(圖8b)更可靠，通過引入構造約束信息提高了巖相的連續性；本文方法考慮了不確定性在各環節的傳遞，得到的巖相后驗概率(圖8c)較傳統多步法(圖8d)更準確，客觀表征了不確定性，為油藏表征、評價提供了重要參考依據。為進一步驗證本文方法的有效性，圖9、圖10分別展示了條件井(B)、盲井(C)測試結果。由圖可見：①在條件井井位處，本文方法考慮了已知井信息，阻抗以及孔隙度反演結果與實測井數據吻合很好，不確定性小(圖9a、圖9b)；多步法由于未考慮已知井信息和不確定性傳遞，波阻抗反演結果的分辨率較低(圖9a)，且孔隙度預測結果偏高(圖9b箭頭處)。②儲層巖相識別結果(圖9c～圖9e、圖10c～圖10e)也驗證了本文方法的有效性。圖11為圖9的巖相識別歸一化混淆矩陣，可見其對角線元素數值之和較大，定量驗證了本文方法的有效性。由盲井測試結果(圖10、圖12)可見：①本文方法的波阻抗反演結果與多步法相似，且與實測曲線基本吻合(圖10a)。②由于多步法存在較大累積誤差，孔隙度反演結果與實測井曲線偏差較大(圖10b)；本文方法同時考慮了波阻抗反演與孔隙度反演之間的耦合關系，孔隙度反演結果與井數據匹配更好(圖10b)。③多步法的混淆矩陣油砂預測值(圖12b，0.8)高于本文方法(圖12a，0.71429)，這是由于多步法在井頂端孔隙度預測值偏大(圖10)導致錯誤油砂預測結果所致，但本文方法的混淆矩陣對角線元素數值之和(圖12a，2.28756)大于多步法(圖12b，2.19781)，因此本文方法的巖相預測精度高于多步法。由圖7～圖12可知，本文方法的油藏參數預測性能較傳統方法更好，可以為油藏表征、評價提供客觀、可靠的參考信息。

圖8 巖相預測結果及其不確定性(a)巖相(本文方法)；(b)巖相(多步法)；(c)巖相后驗概率(本文方法)；(d)巖相后驗概率(多步法)

圖9 條件井(B)測試結果(a)阻抗曲線與本文方法反演結果不確定性信息(底圖)疊合圖；(b)孔隙度曲線；(c)參考巖相；(d)本文方法預測巖相；(e)多步法預測巖相

圖10 盲井(C)測試結果(a)阻抗曲線與本文方法反演結果不確定性信息(底圖)疊合圖；(b)孔隙度曲線；(c)參考巖相；(d)本文方法預測巖相；(e)多步法預測巖相

圖11 圖9的巖相識別歸一化混淆矩陣(a)本文方法；(b)多步法

圖12 圖10的巖相識別歸一化混淆矩陣(a)本文方法；(b)多步法

(27)

3 結論與建議

本文通過分析實際資料提出一種基于構造約束聯合概率反演的油藏參數表征方法。新方法基于貝葉斯理論，利用阻抗和孔隙度的聯合先驗分布和巖相依賴先驗分布，引入線性巖石物理模型，得到油藏參數(阻抗、孔隙度、巖相)后驗概率分布解析表達式，客觀、準確地表征了不確定性傳遞。同時，將地質構造信息和井信息集成到反演過程，提高了反演結果的橫向連續性和分辨率。為驗證新方法的有效性，對M區實際數據集通過條件井和盲井測試，對比、分析了無構造約束多步方法與新方法的反演結果，后者基于線性化模型且服從高斯分布假設，獲得了較好的反演效果，為油藏表征、評價提供了有利依據。針對疊前地震資料，在油藏參數與地震響應之間滿足線性關系假設時，新方法可以實現疊前同時反演油藏參數。

然而，在一些較復雜地區線性表達式常常會降低反演精度，因此新方法很難直接求解非線性反演問題。非線性反演方法可以在一定程度上避免線性反演方法的一些假設和限制，從而提高油藏參數預測精度。

目前求解非線性問題主要有以下策略：①將復雜問題進行線性化近似，然后直接應用本文方法求解，但線性化不可避免地會降低反演精度；②引入全局尋優算法求解非線性反問題，如MCMC采樣算法，但該類算法計算效率較低，限制了實際應用；③在有足夠井數據情況下，利用基于數據驅動類方法訓練已知井數據，獲得油藏參數與地震數據之間的匹配函數式，從而求解非線性反問題。

另外，新方法中模型參數服從高斯分布假設可以較容易地推廣到混合高斯分布；針對一些更復雜的分布，可以利用上述求解非線性問題的策略②或③求解相應問題。

附錄A 后驗分布

根據貝葉斯定理，巖相κt后驗概率密度分布為

(A-1)

式中p(κt)為巖相先驗概率密度分布，積分項內p(mt|κt)、p(mt|yt)均服從高斯分布。

如果參數m的高斯分布記為N(m;μ,Σ)，即p(mt|κt)～N(m;μ1,Σ1)，p(mt|yt)～N(m;μ2,Σ2)，有

N(m;μ1,Σ1)×N(m;μ2,Σ2)=ccN(m;μc,Σc)

(A-2)

式(A-2)等式左邊表示式(A-1)中積分項內高斯分布乘積，式(A-2)等式右邊可表示為

(A-3)

式中d為數據采樣長度。

將式(A-3)展開并進行整理，得

ccN(m;μc,Σc)=

(A-4)

ccN(m;μc,Σc)=

(A-5)

其中

(A-6)

式(A-6)可簡記為

(A-7)

其中

(A-8)

(A-9)