基于疊前概率反演的致密砂巖甜點直接預測方法

趙晨，金鳳鳴，韓國猛，郭淑文，邢興，劉鴻洲

（中國石油大港油田公司，天津 300280）

0 引言

近年來，隨著油氣勘探的目標逐漸轉向非常規儲層，作為致密砂巖儲層預測中至關重要的地質目標——甜點的識別越來越受到重視［1］。甜點是非常規儲層中一種重要的地質和產能概念，用于描述致密砂巖儲層中相對高孔、高滲的區域［2］。因此，要實現甜點預測，首先要獲取砂體的空間展布特征，在砂體內部尋找甜點。然而，砂體及其內部甜點往往具有不同的敏感參數，因此通常需要進行反演參數之間的轉換，這往往增加了累積誤差，降低了反演結果的精度和效率。此外，致密砂巖內部甜點普遍厚度較小，傳統的確定性反演難以滿足預測精度的需求。因此，需要建立一種高分辨率且高效的致密砂巖甜點直接預測方法。

疊前地震反演技術是實現巖性和流體識別的一種重要手段［3-9］。作為一種常用的地震反演方法，概率反演方法利用隨機采樣獲得一系列關于后驗概率分布的樣本，充分利用測井數據所建立的先驗模型約束，具有較高的垂向分辨率。馬爾可夫鏈—蒙特卡洛（MCMC）算法［10］是一種較為常用的概率反演算法，它通過構建若干條與后驗概率分布相同的馬爾科夫鏈實現對貝葉斯公式的求解［11-12］。1970 年，Hastings［13］提出了最為常用的MCMC 算法形式，即Metropolis-Hastings（MH）算法，并且廣泛應用于地球物理反演領域。常規MH 算法收斂速度較慢，且穩定性較差。在MH 算法的基礎上，許多學者提出了一系列改進措施［14-19］。這些措施主要是通過改進搜索步長，充分利用多條馬爾科夫鏈之間的交互性，提高了算法的收斂性和穩定性。然而，MH 及其改進算法多依賴于參數先驗采樣空間的設置，對于一個較為復雜的參數采樣空間，MH 算法往往易陷入局部極值，影響了反演結果的精度與效率。

為此，本文利用疊前地震數據，提出了一種基于改進貝葉斯概率反演的致密砂巖甜點直接預測方法。以R 實際工區為靶區，首先，選取巖性敏感參數和甜點敏感參數，推導基于兩者的Zeoppritz 近似方程，并結合地震褶積模型，建立敏感參數與地震數據的直接關系，避免反演結果因參數轉換帶來的累積誤差；然后，在貝葉斯理論框架下，將統計混合高斯先驗信息［20］與平滑約束先驗相結合，構建反演參數的后驗概率表達式；最后，利用改進貝葉斯概率的反演方法實現巖性與甜點敏感參數的高精度預測。該方法建立反演參數之間的聯合概率分布，計算參數的條件概率分布，并將它作為隨機采樣過程中的約束條件，減少了反演參數的采樣空間，從而提高了反演的效率與精度。

1 方法原理

1.1 地震正演模型

地震正演模型是疊前概率反演的基礎，它主要由Zeoppritz 近似方程［21］和褶積公式組成。R 靶區巖性和甜點的敏感彈性參數分別為泊松比和縱波模量。為減少誤差的累積，從Aki-Richards近似式出發，推導了基于泊松比、縱波模量和密度的Zeoppritz 近似公式。

Aki-Richards近似式為

式中：R(θ)為PP 波反射系數；θ為入射角分別為地下縱波速度、橫波速度和密度的平均值；Δα、Δβ、Δρ分別代表地層分界面兩側的縱波速度、橫波速度和密度的變化量；γ為橫波速度與縱波速度之比。

式中：分別代表地下泊松比、縱波模量的平均值；Δσ、ΔM分別是地層分界面兩側的泊松比、縱波模量的變化量。

將式（2）、式（3）代入式（1），則可得到關于縱波模量M、泊松比σ及密度ρ的新的近似公式，即

給定子波W，地震數據可以利用褶積模型表征［22］，即

式中：d(t，θ)表征t時刻的振幅；τ為臨近點的時間。

利用式（4）、式（5），可以得到敏感參數與地震記錄之間的正演關系，這為貝葉斯概率反演中似然函數的構建奠定了基礎。

1.2 貝葉斯概率反演

在貝葉斯理論基礎上，貝葉斯概率反演根據先驗分布構建反演參數的采樣空間，利用似然函數對樣本進行優選，構建若干條樣本鏈以實現對后驗概率的表征。因此，需要首先根據貝葉斯定理構建反演參數的后驗概率分布表達式。

1.2.1 后驗概率分布構建

由于反演具有多解性，因此概率反演的主要目標是獲取反演參數的可能取值，即后驗概率分布。根據貝葉斯定理，后驗概率分布p(m|d)主要包含先驗概率分布p(m)及似然函數p(d|m)兩部分，可以表征為

式中：m=(M，σ，ρ)T是反演參數；d為觀測地震記錄。

p(m)反映了利用測井數據統計的敏感參數概率分布特征。在實際地層中，不同巖性的彈性參數先驗信息是不同的。假設不同巖相下反演參數滿足不同的單一峰值高斯分布，則統計先驗概率分布可以分解為多個高斯概率分布的疊加。因此，統計先驗概率分布ps(m)滿足混合高斯分布。混合高斯分布主要通過不同高斯分量的均值、方差以及權重系數進行表征，可表示為

為提高概率反演的穩定性，通常需要引入先驗平滑約束［23］。假設平滑約束先驗信息pl(m)與實際模型之差滿足高斯分布，則可以表征為

式中：N為樣點數；m0為平滑約束先驗模型；C0是實際模型與平滑約束先驗模型之差的協方差矩陣。

p(d|m)表征反演參數所合成的地震記錄與實際地震記錄的相似程度，通常利用高斯分布進行表征，即

式中：Ce是地震噪聲的協方差矩陣；f表征正演算子，即地震正演模型。

綜上所述，后驗概率分布可以表示為

式中：c為常數；

J(m)表征反演結果與平滑約束模型、合成地震記錄與實際地震記錄之間的相似程度；gk(m)表征反演結果在第k個高斯分量下所滿足的高斯分布特征。

1.2.2 Metropolis-Hastings（MH）算法

要實現敏感參數的反演，需要求解反演參數的后驗概率分布。基于MH 算法的概率反演方法可以避免貝葉斯公式中復雜的邊緣積分計算問題，構建若干條平穩分布與后驗概率分布特征相同的馬爾科夫鏈可實現對貝葉斯公式的求解，馬爾科夫鏈達到平穩狀態下的一系列樣本，即可表征反演結果的后驗概率分布特征。

MH 算法主要分為隨機采樣與樣本優選兩個部分［24］：①根據參數先驗分布構建反演參數的采樣空間，通過對參數采樣空間的抽樣，獲得來自參數先驗分布的一系列樣本；②利用觀測數據與地震正演模型構建反演參數的似然函數，并利用似然函數與接受準則對馬爾科夫鏈樣本進行優選，使馬爾科夫鏈逐漸收斂到平穩狀態，即后驗概率分布。

其中，接受概率可以表示為

式中：π(?)為平穩分布，即后驗概率分布p(m|d)；q(m′，m*)為模型m′轉移到模型m*的建議分布。

當q(m′，m*) 為對稱分布時，式（14）可以改寫為［25］

平穩分布π(m)即反演參數的后驗概率分布。將式（11）代入式（15），并簡化常數項，即可得到概率反演的接受概率表達式

J(?)為反演的目標函數，表征實際地震記錄與合成地震記錄、反演結果與先驗約束模型之間的總體相似程度，可用來反映反演結果的收斂情況。

1.2.3 基于條件概率分布約束的MH 算法

傳統的MH 方法一般是將測井資料所統計的全局累計概率分布作為反演參數的采樣區間，通常采樣區間范圍較大而降低了反演結果的精度和效率。因此，本文充分考慮反演參數之間的相關性，通過建立反演參數的條件概率分布來減小采樣空間。

針對不同的巖性（即統計先驗概率分布中不同的高斯分量），構建泊松比的條件概率分布

根據上一次的迭代結果，計算當前巖相的估算結果；再根據估算結果分配每個樣點的條件概率分布，即采樣空間。通過對泊松比的條件概率分布進行抽樣，獲得泊松比模型σ*。

在此基礎上，可以用參數之間的聯合概率密度函數計算縱波模量與密度的條件概率分布，即

式中[σ*-ψσ，σ*+ψσ]表示當σ=σ*時泊松比的建議采樣區間。通過對縱波模量及密度的條件概率分布抽樣，即可獲得新的縱波模量M*及密度模型ρ*。

相比傳統 MH 反演，基于條件概率分布約束的MH 反演不是對全局概率分布函數進行采樣，而是通過對條件概率分布進行抽樣以獲得新的樣本，從而減小了反演參數的采樣空間，提高了算法的收斂性與穩定性。

1.3 巖相與甜點概率估計

在反演結果的基礎上，根據不同巖相和甜點類型的彈性參數特征，利用貝葉斯分類方法估算巖相和甜點等離散相態的后驗概率分布，進而得到巖相和甜點的預測結果。

給定巖性敏感參數泊松比預測結果～σ，砂體的后驗概率分布可以表征為［26］

式中：p(Isand) 為砂體的先驗認知，即p(Isand)=Ssand/(Ssand+Sshale)，其中Ssand、Sshale分別代表訓練數據中砂巖、泥巖的累計發生次數，兩者均可通過統計測井數據和解釋結果獲得；p(|Isand)為砂巖所對應彈性敏感參數的概率。

式中A表征甜點類型數目。

2 模型試算及實際數據應用

2.1 一維數據反演測試

為了驗證本文甜點預測方法的可行性，利用真實測井數據構建一維模型（圖1）進行疊前反演測試。給定35 Hz Ricker子波，利用正演模型計算合成角道集，并加入一定的隨機噪聲（SNR=10 dB），得到觀測地震記錄；再對一維模型進行平滑處理，得到反演初始模型；然后分別應用本文所提出的優化概率反演方法與傳統概率反演方法獲得敏感參數的反演結果及甜點估算結果。為提高反演的穩定性，同時運行3 條馬爾可夫鏈。

圖1 一維數據模型

圖2 是一維數據模型所統計的不同巖相及砂巖內部不同甜點類型的敏感參數分布直方圖。由圖2可見，泊松比對巖性較為敏感，縱波模量對砂體內部的甜點較為敏感，且二者的參數先驗概率分布滿足高斯混合先驗分布。

圖2 巖相（左）與甜點（右）敏感參數分布直方圖

圖3 是利用模型數據構建的泊松比與縱波模量的聯合概率分布及縱波模量條件分布計算結果，可以發現條件概率分布的引入可以有效減小參數的采樣空間。

圖3 泊松比與縱波模量的聯合分布特征（左）和縱波模量條件分布計算結果（右）

圖4 是不同方法的敏感參數反演結果及巖相、甜點估算結果，圖5為統計的反演結果相對誤差直方圖。由圖4 和圖5 可以發現，相較于傳統方法，本文方法的多個反演結果構成了更加穩定的后驗概率分布特征，且相對誤差略低，具有更高的精度與穩定性。

圖4 基于不同反演方法的反演結果（SNR=10 dB）

圖5 不同方法反演結果與實際模型的相對誤差直方圖

圖6 為不同方法3 條鏈目標函數值隨迭代次數的變化。由圖可見，相比傳統概率反演方法，基于條件概率分布約束的反演方法的目標函數值在馬爾科夫鏈達到平穩狀態時更小，且收斂速度更快，前者需要25000 次迭代，傳統概率反演方法需要85000 次迭代。

圖6 不同方法目標函數值隨迭代次數的變化

圖7 是反演結果與實際模型參數分布直方圖的對比，可以發現反演結果的參數分布與實際模型參數分布具有較高的相似性，這說明反演結果滿足高斯混合先驗分布假設。

圖7 模型與反演結果反演參數統計直方圖對比

模型測試表明本文方法可以有效識別巖相及砂體內部甜點。

2.2 實際數據測試

為了驗證本文方法在實際應用中的可行性，利用中國北部R 區塊實際數據進行測試。R 區塊位于斜坡構造帶，發育水下分流河道、河口壩砂體，其中高能、厚層河道砂體粒度粗、分選好，儲層物性好，是致密砂巖甜點發育的主要區域；砂體內部，甜點具備相對較高的孔隙度與低彈性模量特征，因此疊前地震反演是實現該儲層甜點預測的有效手段。由于疊前地震資料品質的限制，為了提高反演的穩定性，將疊前角度道集按2°～10°、10°～18°、18°～26°分別進行角度疊加，三個部分疊加數據的中心角度分別為6°、14°、22°。利用井數據提取地震子波，并對井插值數據進行平滑處理得到初始模型。

圖8 為疊后與角度部分疊加地震剖面，可見研究層段地震分辨率較低，將會制約甜點的預測精度。

圖8 疊后與角度部分疊加地震剖面

圖9 為測井數據統計的不同巖相與砂體內部不同類型甜點泊松比、縱波模量的分布直方圖，可以確定巖性與甜點敏感參數分別為泊松比與縱波模量。

圖10與圖11分別為彈性參數反演剖面與甜點估算結果。與測井曲線對比可以發現，反演結果與井曲線吻合較好（圖10），且甜點估算結果符合測井解釋結果（圖11）。

圖10 本文方法疊前反演結果

圖11 本文方法甜點概率計算結果（上）與甜點的估算結果（下）

圖12 為沿層的甜點概率平面屬性。由圖可見，甜點主要發育于區塊西部與南部，這與地質認識相吻合，進一步驗證了本文方法反演結果的合理性。

圖12 甜點概率沿層平均屬性圖

圖13 為井旁道反演結果與實際井數據的對比，可以發現井旁道反演結果與井曲線吻合較好，且合成地震記錄與實際地震記錄相匹配，這說明了本文反演方法的有效性。

圖13 井旁道反演結果與實際井數據的對比

圖14 為井旁道反演結果與測井數據分別統計的參數分布直方圖。由圖可見，本文方法反演結果的參數分布特征與統計參數先驗一致，這說明了本文方法反演結果統計特征與先驗假設一致。

圖14 測井數據與井旁道反演結果統計直方圖對比

綜上所述，本文基于疊前概率反演的甜點直接估算方法可以有效應用于實際工區。

3 結束語

本文在Zoeppritz近似方程的基礎上，推導了基于致密砂巖巖性與甜點敏感參數的反射系數方程，并將統計混合高斯先驗與平滑約束先驗相結合，建立反演參數的后驗概率分布表達式，通過改進MH 算法實現對后驗概率分布的求解，從而建立基于疊前概率反演的甜點直接估算方法。在反演過程中，通過引入條件概率分布約束，減小反演參數的采樣空間，不僅加快了算法的收斂，還提高了反演的穩定性。模型測試與實際數據應用均表明，基于疊前概率反演的甜點直接估算方法可以有效地實現甜點預測。