高分辨率波阻抗貝葉斯序貫隨機反演

李祺鑫 羅亞能 張 生 張 璐 楊亞迪 黃捍東

(①中海油研究總院，北京 100028； ②東方地球物理公司物探技術研究中心，河北涿州 072751；③太原理工大學礦業工程學院，山西太原 030024； ④中國石油勘探開發研究院，北京 100083；⑤中國石油大學(北京)地球物理學院，北京 102249)

0 引言

地震波阻抗反演作為儲層描述的工具之一，在儲層預測方面具有重要作用[1-5]。地震波阻抗反演方法眾多，可以分為確定性反演與統計學反演兩大類[6-8]。確定性反演是一類成熟的方法，這類方法的共性是構建一個目標函數，求取使目標函數取得極小值的解作為反演結果[9-10]。確定性反演常見的方法包括模型法反演[11]和稀疏脈沖反演[12]，在特定的解約束條件下，只能進行單一的最優解估計，無法評價反演結果的不確定性。實際上，地震反演是一個多解性問題，相比于確定性反演只給出一個某種范數意義下的最優解，統計學反演則嘗試用概率語言完整地表達解空間及其每個解的可能性。因此，統計學反演給出的是一簇滿足約束條件的解，而非一個唯一解[13]。對于統計學反演輸出的多個結果，利用概率統計的手段評價反演結果的不確定性，是分析勘探風險的一個重要方面[14-15]，這是統計學反演的優勢之一。為了完全探索反演解空間的結構，統計學反演利用采樣求解的策略，反演結果的分辨率高于確定性反演，在薄層刻畫方面取得了很好的效果[16-18]，這是統計學反演的優勢之二。

作為統計學反演方法的一個分支，Haas等[19]首次將地質統計學應用于地震反演，將地震資料作為約束條件，控制測井資料序貫模擬結果的取舍，構成了傳統地質統計學反演的基本思想。隨著研究的深入，地質統計學由兩點地質統計學向多點地質統計學發展，地質統計學反演相應地也由兩點地質統計學反演轉向多點地質統計學反演[20-22]。González等[22]首次將多點地質統計學算法與地震反演結合反演河道砂波阻抗，探討了多點地質統計學反演的可行性與優越性。目前，多點地質統計學反演的難點在于如何合理建立訓練圖像以及提高計算效率，這兩方面制約著多點地質統計學反演的實際應用效果[23-24]。

目前地質統計學反演主流算法均為地質統計學隨機模擬方法與馬爾科夫鏈—蒙特卡洛方法的某種結合[25-28]，并沒有將地震反演、測井數據約束、隨機模擬方法統一在一個理論框架下。為此，本文首先在貝葉斯理論框架下，將地震反演理論、測井數據約束、地質統計信息統一在波阻抗求解這個問題中，推導了基于對數波阻抗的地震數據與測井數據聯合方程，然后通過序貫高斯隨機模擬算法對反演方程進行高效采樣求解，最后通過理論模型試算和實際資料應用驗證了方法的可行性和有效性。

1 方法原理

1.1 地震正演模型

假設地層界面兩側為各向同性、彈性介質，當地震波垂直入射時，反射系數r(t)和波阻抗z(t)的關系可表達為[29]

(1)

式中t為時間。在小反射系數假設條件下，將式(1)改寫為矩陣形式[30]

(2)

式中：r=[r1,r2,…,rN-1]T代表一道反射系數序列，下標N代表地震道的時間采樣點數；z=[z1,z2,…,zN]T代表聲波阻抗；而

(3)

代表差分矩陣。根據褶積模型,觀測地震數據dobs為

dobs=Wr+e=WDm+e=Gm+e

(4)

式中：W為N×(N-1)維子波褶積矩陣；e為誤差項；m=(1/2)lnz為波阻抗中間變量(模型參數)；而

G=WD

為N×N維正演算子矩陣。

1.2 貝葉斯空間結構約束后驗模型

地震反演的目的為從dobs中估計m，但由于地震資料的帶限特征，存在無數解滿足式(4)。傳統確定性反演直接對解空間添加約束條件[31]，如最小長度解、最光滑解、或者最稀疏解等，從而獲得單一的解估計。在此類方法中，正則化參數的選擇過于主觀，反演結果受人為因素影響大，而且不能評估反演結果的不確定性。

貝葉斯反演方法能夠綜合觀測數據中的多種先驗信息，通過概率統計理論獲得完備的后驗解空間，并且能夠表征反演結果的不確定性[9,13,32]。在貝葉斯反演中，反演結果由后驗概率密度函數σ(m)表示，其為先驗概率密度函數ρ(m)與地震似然函數L(m)的乘積，即

σ(m)=cρ(m)L(m)

(5)

式中c為歸一化常數。

假設m服從多維高斯分布，即m～NN(mprior,Cm)，則

(6)

式中：mprior=lnz0為N維列向量，表示波阻抗的先驗低頻背景模型，z0為初始波阻抗模型；c1為歸一化常數；Cm為N×N維模型先驗協方差矩陣

Cm(p,q)=C0[1-v(h)]

(7)

(8)

式中：c2為歸一化常數；Cd為N×N維地震數據協方差矩陣，通過測井合成記錄與實際井旁地震道資料估算得到。

(9)

(10)

1.3 井震聯合貝葉斯序貫隨機反演

當建立了后驗期望(式(9))和后驗協方差矩陣(式(10))表達式后，一般通過兩種傳統方法獲得反演結果： ①直接代入相應矩陣求解[34]； ②采用蒙特卡洛方法實現隨機采樣[35]。在第一種方法中，由于Cm和Cd非常龐大，需要很大的物理存儲空間存儲。如一個二維剖面包含200個時間采樣點、300道數據點，在地震反演中僅存儲Cm所需的物理存儲空間為(200×300)2×8B≈26.8GB，因此在實際三維地震資料反演中，對所有點進行一次性直接求解幾乎不太現實。在第二種方法中，盡管蒙特卡洛方法能較好地對高維數據降維處理，但是此類方法通常需要經過一系列的“拒絕—接受”重復迭代采樣過程，計算效率低[36]。

不同于以上兩種方法，地質統計學中序貫高斯模擬能夠充分利用求解問題的線性和高斯分布特性，是一種非迭代的高效采樣方法。此方法并非一次性求解式(9)和式(10)中龐大線性矩陣的所有點，而是在模型空間中隨機訪問單個點，并且采用“鄰域”思想將大型矩陣求解問題轉換為許多小型線性矩陣求解問題，適用于任意高維問題。

序貫模擬的基本思想為將已模擬點作為條件數據，同原始已知點一起，在給定鄰域內參與下一個點的計算，直至模擬完成隨機路徑中所有節點。借鑒這種思想，本文采用井震聯合的貝葉斯序貫隨機反演方法，為了求解方便，在鄰域內將觀測數據劃分為A、B兩類數據。A類數據為模型參數的直接測量值，包括測井數據和先前模擬的點，也稱為“硬數據”或者“點數據”； B類數據為模型參數的線性平均觀測值，即通過線性正演算子獲得的測量數據，也稱為“軟數據”或者“體數據”，這里指地震數據。聯合兩類觀測數據，式(4)可改寫為

(11)

dobs(A+B)=GA+BmA+B+eA+B

(12)

式中：dobsA、GA、mA和eA分別代表A類觀測數據、正演算子、模型參數和誤差項；dobsB、GB、mB和eB分別代表B類觀測數據、正演算子、模型參數和誤差項；dobs(A+B)、GA+B、mA+B和eA+B代表A、B兩類數據的聯合參數。這里需要注意的是，與地震正演算子GB不同，直接觀測的A類數據正演算子GA為簡單的單位矩陣，不依賴于任何物理規律。

Cd(A+B)]-1[dobs(A+B)-Gmprior(A+B)]

(13)

Cd(A+B)]-1GA+BCm(A+B)

(14)

式中：mprior(A+B)和dobs(A+B)均為(NA+NB)維列向量，NA、NB分別為鄰域內A、B類數據的采樣點數；Cm(A+B)和Cd(A+B)分別為(NA+NB)2維模型協方差矩陣和數據協方差矩陣，表達式為

(15)

式中：CmAA、CmBB分別為A、B類數據的模型協方差矩陣；CmAB為A、B類數據的互協方差矩陣；CdBB為B類數據的數據協方差矩陣。在Cd(A+B)中，假設A類數據為準確值，并假設A、B兩類數據的誤差不相關，因此A類數據的協方差矩陣和兩類數據的互協方差矩陣都為0。

為了專心準備金腰帶，鄒推掉了很多品牌代言。他現在每年要負擔200萬元～300萬元人民幣訓練費用，他靠比賽獎金維持收入，當然，這離一名真正的明星職業拳手依然相差甚遠。

圖1 算法流程

2 數值算例

為了驗證貝葉斯序貫隨機反演算法的有效性，利用無條件序貫高斯隨機模擬算法生成一個二維薄互層波阻抗模型，作為實際地下地層的參考模型，并且從二維模型中抽取一道作為測井數據，即A類數據(圖2)。

圖3為根據圖2生成的合成地震記錄與常數背景初始模型，其中合成地震記錄(圖3a)作為B類數據，常數背景模型(圖3b)在地質統計學模擬中作為全局期望，在地震反演中作為反演初始模型。以合成地震記錄(圖3a)中的測井數據(A類)和地震數據(B類)作為觀測數據，分別進行B類數據最小二乘法確定性反演、A類數據序貫高斯模擬及A、B兩類數據聯合貝葉斯序貫反演測試，分析貝葉斯序貫反演方法的優勢。圖4為地震數據(B類)最小二乘法反演結果，反演輸入數據只有地震數據dobsB和地震數據協方差矩陣Cd，此時的最小二乘解為

mLSI=m0+GT(GGT+Cd)-1(dobsB-Gm0)

(16)

由圖4可見，與實際波阻抗參考模型(圖2)相比，最小二乘反演結果的縱向分辨率與地震資料相當，只能反映參考模型中參數的變化趨勢，不能刻畫測井曲線的細節。事實上，最小二乘反演為一種過度光滑的估計方法，而且沒有添加先驗信息約束，反演結果僅依賴地震資料，因此反演結果分辨率低。

圖2 合成波阻抗參考模型

模型共301道,道間距為10m,縱向采樣點數為150，采樣間隔為1ms。對數波阻抗均值為16.15kg·m-2·s-1，方差為0.01(kg·m-2·s-1)2，水平變差函數變程為1000m(或100道)，垂直變差函數變程為10ms(或10個采樣點)，設置變差函數為球狀變差函數

圖3 根據圖2生成的合成地震記錄(a)與常數背景初始模型(b)

圖4 地震數據(B類)最小二乘法反演結果

圖5為測井數據(A類)序貫高斯模擬2次實現結果，反演輸入數據只有測井數據dobsA和模型協方差矩陣Cm。由圖5可見：模擬結果的分辨率與參考模型(圖2)相當，且在井點處的結果與測井資料完全吻合；由于缺乏地震數據約束，模擬結果的不確定性高(圖5箭頭處)，2次實現結果的差別非常大。因此在缺乏井資料時，傳統的地質統計學方法不適用于井間模型參數預測。

圖5 測井數據(A類)序貫高斯模擬2次實現結果

圖6為井震聯合(A+B類)數據貝葉斯序貫反演3次實現結果，輸入數據為模型協方差矩陣Cm和數據協方差矩陣Cd。由圖6可見，3次實現的分辨率與參考模型(圖2)相當，不僅在井點處與測井曲線完全吻合，且降低了實現結果的不確定性(圖6箭頭處)，3次實現結果整體特征相似，反映了地震數據約束的一致性；在小尺度條件下，3次反演結果間具有微小的差異，這種差異反映了不確定的地下介質非均質性。

圖7為不同結果合成地震記錄與實際地震記錄殘差。由圖可見：最小二乘反演結果(圖7b)和序貫高斯模擬(圖7c)的殘差均與實際噪聲(圖7a)相差較大，其中后者與實際噪聲的差異更大，而前者逐道反演的痕跡明顯，反映了沒有利用模型參數橫向相關性的缺點；貝葉斯序貫反演結果的殘差(圖7d)與實際地震噪聲(圖7a)在幅度和形態上較接近(圖7a和圖7d箭頭處)。上述結果表明，井震聯合數據貝葉斯序貫反演方法得到的模型參數精度最高。

圖6 井震聯合(A+B類)數據貝葉斯序貫反演3次實現結果

圖8為水平、垂直方向反演結果變差函數。由圖可見：由于A+B類數據貝葉斯序貫反演過程使用了模型的空間結構信息Cm，因此所得水平、垂直變差函數均與參考模型變差函數匹配較好，其中水平方向變程約為1000m(圖8a)，垂直方向變程約為10ms(圖8b)；最小二乘反演結果(B類)由于沒有利用模型參數的空間結構信息，反演結果的水平(圖8a)、垂直(圖8b)變差函數與參考模型變差函數差別較大。

模型測試分析結果表明，貝葉斯序貫隨機反演方法采用鄰域思想將高維數據降維處理，適用于任意的高維數據反演，同時能夠克服最小二乘法反演的光滑效應以及序貫高斯模擬不確定性大的問題，可獲得高分辨率且相對穩定的反演結果。

圖7 不同結果合成地震記錄與實際地震記錄殘差

圖8 水平(a)、垂直(b)方向反演結果變差函數

3 實際應用

為了進一步檢驗貝葉斯序貫隨機反演方法的實際應用效果，以中國東海盆地西湖凹陷N井區進行試驗，試驗數據來自一條過井三維地震剖面(圖9)，其中測井合成地震記錄與實際井旁道對應較好(圖9箭頭處)。圖10為W1井對數域彈性參數正態分布檢驗Q-Q圖。由圖可見，數據點基本位于正態分布參考線之上，表明貝葉斯序貫隨機反演方法對于研究區的適用性。

圖9 實際地震剖面與測井合成記錄標定

縱向時間采樣間隔為1 ms，采樣點數為620個，道間距為25 m，共400道，地震資料主頻為32 Hz。黑色曲線為井震標定后的測井合成地震記錄，測井合成記錄選用主頻為32 Hz的雷克子波

圖10 W1井對數域彈性參數正態分布檢驗Q-Q圖

圖11為模型法、稀疏脈沖法對圖9的反演結果，圖12為貝葉斯序貫隨機反演3次實現結果。由圖可見，3種反演方法得到的波阻抗剖面在整體上具有相似性，但模型法反演結果在縱向上更光滑(圖11a黑色箭頭處)，稀疏脈沖法反演結果則展示了一個更“塊狀”化的結果(圖11b黑色箭頭處)，貝葉斯序貫隨機反演方法則提供了多個與井曲線分辨率相當的反演結果(圖12)。

事實上，之所以產生上述反演結果是由采用的反演方法決定的。模型法反演直接對低頻初始模型進行迭代擾動，當正演模擬數據與觀測數據達到一定的匹配標準時停止迭代；稀疏脈沖反演是在反射系數稀疏準則的約束下，對地震資料反褶積獲得寬帶反射系數，然后再通過遞推等算法獲得絕對阻抗；貝葉斯序貫隨機反演將測井數據和已模擬點作為硬數據，在地質統計信息和地震數據的共同約束下，參與計算下一個點，因此測井信息在遍歷隨機路徑的過程中得到傳播，從而使反演結果具有高分辨率特征，且橫向變化特征清晰、自然。

圖11 模型法(a)、稀疏脈沖法(b)對圖9的反演結果

圖13為盲井反演結果、序貫隨機反演概率解。由圖可見：在恢復井曲線波阻抗的中低頻相對變化時，3種方法給出了等價結果，而序貫隨機反演方法則展示了更高頻的阻抗變化細節(圖13a)；對比概率分布與真實井曲線，能夠獲得在當前數據條件下反演結果的不確定性程度(圖13b)。模型法反演與稀疏脈沖反演這類確定性反演方法，由于只能給出一個唯一解，無法評價反演結果的不確定性與風險，而風險評估在勘探決策中具有重要作用，這是利用貝葉斯序貫隨機反演方法進行地質解釋的另一優勢。圖14為三維貝葉斯序貫隨機反演最大后驗解。由圖可見，沿目的層頂界切片明顯展示了兩條河道特征(圖中箭頭處)，與研究區三角洲前緣水下分流河道沉積背景相符合。

圖12 貝葉斯序貫隨機反演3次實現結果

圖13 盲井反演結果(a)、序貫隨機反演概率解(b)

圖14 三維貝葉斯序貫隨機反演最大后驗解

4 結束語

本文采用高分辨率貝葉斯波阻抗序貫隨機反演方法，在貝葉斯理論框架下整合地震數據、測井數據和地質空間結構先驗信息，并且將線性反演理論與克里金算法相聯系，采用序貫模擬方法高效隨機實現后驗概率分布。數值算例表明：與常規波阻抗最小二乘法反演相比，該方法分辨率高、反演結果受先驗統計學參數約束，且輸出的多個波阻抗反演實現可用于不確定性評價；與純粹基于測井數據的序貫高斯隨機模擬方法相比，該方法最大的優點就是考慮了地震數據約束，降低了反演結果的不確定性。針對實際數據，對比、分析了模型法反演、稀疏脈沖法反演、貝葉斯序貫隨機反演結果，表明貝葉斯序貫隨機反演方法在提高反演縱向分辨率與不確定性評價方面具有良好的應用前景。