確定性反演協同約束的疊后隨機地震反演方法

張豐麒 劉俊州 劉蘭鋒 時 磊 韓 磊

(①頁巖油氣富集機理與有效開發國家重點實驗室，北京 100083；②中國石化彈性波理論與探測技術重點實驗室，北京 100083；③中國石化石油勘探開發研究院，北京 102206)

0 引言

高分辨率地震反演技術對于預測薄儲層、薄互層油氣藏具有至關重要的意義。確定性反演技術發展日趨成熟，其代表性方法包括遞推反演、基于模型的反演、稀疏脈沖反演以及有色反演等。雖然確定性反演更忠實于地震資料本身，然而這也決定了反演結果難以突破地震資料的極限分辨率，方法本身受限于地震資料頻帶。隨機地震反演結合隨機模擬技術與地震反演理論，可以產生高分辨率反演結果。Joumel等[1]、Hass等[2]首先提出隨機地震反演方法。Dubrule[3]利用地震數據約束地質建模。De-beye等[4]、Sams等[5]將模擬退火和蒙特卡洛—馬爾科夫鏈(MCMC)算法引入隨機反演，提高了計算效率。Mosegaard等[6]將貝葉斯理念引入地球物理反演問題，修正了MCMC算法系列的Metropolises—Hastings(M-H)算法，提出了更高效的擴展M-H算法，該算法只需評估概率轉移前、后兩個狀態的似然函數值。Hansen等[7]結合序貫模擬和Gibbs采樣提出了序貫Gibbs采樣，根據地質統計學從復雜先驗分布中抽樣；實驗證明，該方法配合擴展M-H算法可快速收斂于目標分布，提高了隨機反演的計算效率。張繁昌等[8]將此隨機反演思想用于疊后波阻抗反演，并在擴展M-H算法中引入退火因子，進一步提高了計算效率。張廣智等[9]將經典M-H算法用于疊前地震反演，但未考慮子波效應，屬于帶限反演。王保麗等[10]引入快速傅里葉變換—滑動平均(FFT—MA)譜模擬與逐步變形算法(GDM)，形成了一種快速疊后隨機反演算法。劉興業等[11]聯合多點地質統計學反演與序貫高斯模擬隨機反演巖相與物性參數。趙晨等[12]聯合直接序貫協模擬與擴展M-H算法，構建了基于全局迭代反演策略的疊前隨機反演，并進一步引入平滑約束和二階差分橫向約束，以提高反演結果的橫向連續性。本文在前人的工作[6-8]基礎上，提出了一種更適合實際資料的高分辨率疊后隨機反演算法。該算法以序貫Gibbs擾動模擬、擴展M-H算法為核心，進一步引入層序地層網格，自適應融入地質統計學、構造以及沉積模式等信息，整個反演過程無需計算地層的局部傾角或進行復雜的坐標轉換。針對隨機地震反演結果的高頻成分仍具有較大的不確定性，造成不同實現展示的儲層特征差異較大，本文結合序貫Gibbs擾動模擬與同位協同克里金，通過引入確定性反演結果的協同約束，進一步限定隨機反演候選解空間，從而降低隨機地震反演高頻成分的不確定性。本文首先結合地質統計學、貝葉斯原理以及反演理論嚴謹、詳實地推導了算法；隨后在實際資料試算中，通過對比確定性反演結果以及隨機地震反演的四個不同實現，驗證了算法的有效性。

1 層序地層網格

傳統反演技術通常采用均勻地震網格計算(圖1b)，網格的橫向間距為CDP間距，垂向間距為地震采樣率。本文的隨機反演基于層序地層網格(圖1a)，網格的橫向間距為CDP間距，但是垂向間距并不等于地震采樣率。

基于構造解釋層位和沉積模式建立層序地層網格，通過結合不同的沉積模式(平行頂、底或平行于頂及平行于底等)在層段中建立層序地層，每個層序地層與CDP間隔線相交便剖分出層序地層網格(圖1a)。利用層序地層網格進行隨機反演，自適應融入了地層構造和沉積模式信息，因此只需沿著網格的水平方向利用水平變差約束隨機反演結果的橫向連續性；若采用均勻地震網格，則需要進一步計算地層的局部傾角或進行復雜的坐標變換，從而增加反演算法的難度以及降低計算效率。

圖1 層序地層網格(平行頂、底沉積模式)(a)與均勻地震網格(b)示意圖

2 先驗分布

模型參數的先驗分布作為貝葉斯反演框架的重要組成部分，具有限定反演解空間以及減少反演過程的不確定性和多解性等作用。在確定性反演中，先驗分布轉化為反演目標函數的約束項，可有效改善反演的病態性和不適定性，根據先驗分布是否服從“長尾分布”，決定反演結果是“塊化”還是“光滑”[13]；在隨機地震反演中，先驗分布還可有效限定解的搜索空間，如果利用先驗分布產生候選解，則能有效加速反演的收斂過程[7]。

通過在先驗分布中引入空間約束改善單道隨機反演的橫向不連續性問題。假設工區內待反演的地震道數為K，K道波阻抗的聯合概率服從多變量高斯分布，記為

(1)

式中：m=[1m2m…Km]為由K道波阻抗組成的列向量，每道波阻抗jm=[jm1jm2…jmN](j=1，2，…，K為道號)包含N個元素(N表示層序地層網格采樣點數)，因此m的維數為K×N；μm為期望，通過分層段統計工區內的測井數據或者建立反演低頻模型得到；Cm為協方差矩陣，為K×N行、K×N列的對稱矩陣，直接影響反演結果的垂向分辨率和橫向連續性。

Cm的對角線元素表示波阻抗的方差，決定了反演結果的動態范圍，通過分層段統計測井數據獲取；Cm的非對角線元素表示同一道、不同采樣點，不同道、相同采樣點以及不同道、不同采樣點之間的協方差。其中同一道、不同采樣點之間的協方差與垂向變差有關，通過分層段統計測井數據獲取，決定了反演結果的垂向分辨率；不同道、相同采樣點之間的協方差與水平變差有關，通過統計地震數據沿層切片獲取，決定了反演結果的橫向連續性；根據變差函數的幾何各向異性，計算水平變差和垂向變差得到不同道、不同采樣點之間的協方差[14]。

3 似然函數

在貝葉斯反演理論中，利用似然函數表征合成地震數據與實測地震數據之間的相似性，進而建立反演參數與地震數據之間的聯系。

經過前期去除相干噪聲(異常噪聲、線性干擾、面波以及多次波等)處理，殘留噪聲僅剩隨機噪聲，因此利用高斯分布描述疊后地震反演的似然函數。由于不同采樣點、不同道之間的隨機噪聲獨立不相干，因此似然函數l(d|m)為

l(d|m)∝

(2)

由于反演波阻抗基于層序地層網格采樣，地震數據基于均勻網格采樣，因此須在G(jm)中引入重采樣矩陣R(附錄A)，將波阻抗在均勻網格重采樣，并結合近似公式計算反射系數r

r≈?ln(m)/2

(3)

最后將反射系數與子波褶積，便得到js，則

(4)

式中：W為子波褶積矩陣；D為一階差分矩陣；ln(·)表示對向量中每一個元素取自然對數；Rj為第j道重采樣矩陣。可見，雖然反演波阻抗基于層序地層網格采樣，但仍在地震網格進行褶積正演。

式(2)為無測井數據約束的似然函數，當有Kw口井的測井數據參與約束時，則似然函數為

l(d,mw|m)∝

(5)

(6)

4 后驗分布的降維分解

根據貝葉斯定理，整合反演參數的先驗分布p(m)與似然函數l(d,mw|m)，可以得到反演參數的后驗分布

p(m|d,mw)∝p(m)·l(d,mw|m)

(7)

將式(1)與式(5)代入式(7)，得

(8)

實際上，確定性反演結果等價于反演參數的最大后驗概率解；隨機反演結果則可以看作是對后驗概率分布的一次隨機抽樣[15]。式(8)表明，后驗分布是復雜的高維分布，無法采用直接抽樣的方式從后驗分布中獲取一次隨機反演實現。本文結合序貫模擬的思想，將后驗分布降維分解

p(m|O)=p(1m|O)p(2m|1m,O)…

p(Km|1m,2m,…,K-1m,O)

(9)

式中O為觀測數據的集合，包含地震數據和參與約束的測井數據，即O=d∪mw。因此jO代表第j道觀測數據，如果第j道是井旁道，則jO=jd∪jmw，反之jO=jd。進一步，假設第j道的反演參數僅與第j道觀測數據相關，而與第i道(i≠j)觀測數據獨立不相關，則式(9)簡化為

p(m|O)=p(1m|1O)p(2m|1m,2O)…

p(Km|1m,2m,…,K-1m,KO)

(10)

由式(10)可見，對后驗分布p(m|O)的抽樣，可以借助序貫模擬的思想實現，即先對p(1m|1O)抽樣，再將抽樣結果1m作為條件數據，參與對p(2m|1m,2O)的抽樣，得到2m的抽樣結果，以此類推，直到完成最后一道Km的抽樣[17-21]。

利用序貫模擬的思想，將復雜的高維(維度為K×N)后驗分布p(m|O)化簡為K個低維(維度為N)后驗分布p(jm|1m,2m,…,j-1m,jO)的乘積。然而p(jm|1m,2m,…,j-1m,jO)仍然不便于直接抽樣。為此，進一步結合貝葉斯定理，將p(jm|1m,2m,…,j-1m,jO)轉化為

p(jm|1m,2m,…,j-1m,jO)=

(11)

由于jO為已知，因此p(jO)為一個不影響后驗概率分布形狀的常數，則

p(jm|1m,2m,…,j-1m,jO)∝

p(jm|1m,2m,…,j-1m)·l(jO|jm)

(12)

式中l(jO|jm)為第j道的似然函數。結合式(2)和式(5)可知，當該道為非井旁道時

l(jO|jm)∝

(13)

當該道為井旁道時

(14)

針對式(12)的抽樣，利用序貫Gibbs擾動模擬從先驗分布p(jm|1m,2m,…,j-1m)產生候選解，并配合擴展M-H算法進行概率轉移，通過多次迭代，便可收斂于后驗分布p(jm|1m,2m，…，j-1m,jO)。

5 序貫Gibbs同位協同擾動模擬與擴展M-H概率轉移

序貫Gibbs擾動模擬結合序貫模擬的思想與Gibbs采樣，從復雜先驗分布中抽樣。該算法配合擴展M-H算法形成一種改進的MCMC算法，進而可以從復雜的后驗分布中抽樣。Hansen等[7]指出，序貫Gibbs擾動模擬可以從更逼近真實情況的復雜先驗分布中抽樣，因此反演結果也更逼近真實解，并且可較大幅度地提高隨機反演效率。

擴展M-H算法由經典M-H算法發展而來，屬于MCMC算法的一個分支。MCMC算法利用建議分布產生候選解，并以一定概率接受候選解，經過多次迭代，便可以收斂于目標分布[6]。當建議分布為后驗分布對應的全條件概率分布時，MCMC算法轉化為Gibbs采樣；當建議分布為均勻分布時，MCMC算法轉化為經典M-H算法；當建議分布為先驗分布對應的條件概率分布時，MCMC算法轉化為擴展M-H算法。

相對于經典M-H算法，擴展M-H算法的候選解來自先驗分布而非廣泛的均勻分布，因此有效地限定了解空間的搜索范圍，其反演結果也可以更快速地收斂于真實解[6]。另外，與經典M-H算法相比，擴展M-H算法在計算轉移概率時，無需評估轉移前、后兩個狀態的后驗概率值，只需評估似然函數值，由于省略了評估先驗分布概率值，因此提高了計算效率。

實際上，隨機反演的高頻成分(60Hz～Nyquist頻率)僅受變差函數和測井數據約束，并不受地震數據約束。當工區內參與約束的測井數據較少時，反演結果的高頻成分仍然有較大的不確定性，從而造成隨機反演的不同實現展現的儲層特征差異較大，因此對后續的儲層精細預測帶來一定困擾。為此，本文結合同位協同克里金修改式(12)，提出確定性反演協同約束的隨機反演算法，利用確定性反演結果約束隨機反演的高頻成分。通過在先驗分布中引入確定性反演結果作為條件數據，則式(12)修正為

p(jm|1m,2m,…,j-1m,jO,jξ)∝

p(jm|1m,2m,…,j-1m,jξ)·l(jO|jm)

(15)

式中：jξ為第j道波阻抗確定性反演結果；p(jm|1m,2m,…,j-1m,jO,jξ)為修正后驗分布；p(jm|1m,2m,…,j-1m,jξ)為修正先驗分布。

式(15)表明，當前道的波阻抗隨機反演結果不僅與其他道的隨機反演結果相關，還與當前道的確定性反演結果相關。

由于從先驗分布p(jm|1m,2m,…,j-1m,jξ)抽樣產生候選解，需要結合序貫Gibbs擾動模擬和同位協同簡單克里金(序貫Gibbs同位協同擾動模擬)。以后驗分布p(jm|1m,2m,…,j-1m,jO,jξ)為目標分布，結合序貫Gibbs同位協同擾動模擬與擴展M-H算法構建的MCMC算法流程如下。

(2)從第j道N個采樣點中隨機選擇一點i(i=1,2,…,N)；

(3)利用同位協同簡單克里金計算第i個采樣點的條件均值和條件方差(附錄B)，進而構建第i個采樣點的局部條件概率密度函數p(jmi|jm-i,1mi,2mi,…，j-1mi,jξi)

式中：jξi為第j道的第i個采樣點的波阻抗確定性反演結果；jmi為第j道的第i個采樣點的波阻抗，其他以此類推；jm-i為第j道的第i個采樣點以外的波阻抗。

(5)計算接受概率

(6)從均勻分布U(0,1)產生一個隨機數Б，如果Б≤α，則jm=jmΔ，反之則維持jm不變。

為了降低克里金方程組的維度以提高計算效率，針對jmi的條件數據，只選取其他道的第i個采樣點的波阻抗1mi,2mi,…,j-1mi，而忽略其他道的第i個采樣點以外的波阻抗。引入這種近似合乎地質認識，因為jm基于層序地層網格采樣，同一采樣點處于同一層序地層。由于不同道、同一層序地層之間的波阻抗的相關性較強，而不同道、不同層序地層之間的波阻抗相關性較弱，因此可以忽略。在實際地震數據反演時，還可以結合垂向變程和水平變程，進一步縮減jmi條件數據的有效個數。

以上為某一道的算法流程。對于整個三維工區而言，須結合序貫模擬和改進的MCMC算法(圖2)：

圖2 反演流程

(1)從三維工區中隨機選擇一道j，結合水平變程，將當前道鄰域內已反演道以及當前道的確定性反演結果作為條件數據；

(2)利用序貫Gibbs同位協同擾動模擬與擴展

M-H算法對目標分布p(jm|1m,2m,…,j-1m,jO,jξ)抽樣，獲取當前道的隨機反演結果jm；

(3)將jm在地震網格重采樣；

(4)重復步驟(1)～(3)，直到完成三維工區所有道的反演。

6 模型試算

利用模型數據(圖3)驗證所提隨機反演算法的效果。由于該模型數據為一維數據，因此不用考慮橫向變差，只需考慮波阻抗的期望、方差以及垂向變差。通過對時間域波阻抗曲線進行10Hz低通濾波獲取波阻抗的期望，在此基礎上計算波阻抗的方差。圖4為垂向實驗變差與垂向擬合變差，其中變差函數模型由0.48倍的指數模型與0.52倍的高斯模型混合而成，兩個模型的垂向變程均為3ms。

圖3 模型數據

圖4 垂向實驗變差(黑點)與垂向擬合變差(紅線)

圖5為反演結果。由圖可見：①確定性反演結果與實測曲線整體吻合度較高，但是垂向分辨率較為了更逼近真實情況，圖a選自一段實際測井曲線，經過Backus濾波(500Hz)再線性插值，得到圖b(采樣間隔為1ms)；在此基礎上計算法向反射系數，并與主頻為30Hz的雷克子波褶積，并加入高斯隨機噪聲得到圖c低，僅為實測曲線的一種近似平滑(圖5a)；②隨機反演的高頻細節與實測數據并不完全吻合(圖5b、圖5c)；③基于確定性反演協同約束的隨機反演(圖5b)的“灰色條帶”寬度明顯小于常規隨機反演(圖5c)，說明前者通過引入確定性反演協同約束，有效降低了隨機反演高頻成分的不確定性。造成以上現象的原因為：由于合成記錄中含有隨機噪聲，且反演中無測井數據約束，合成記錄的帶限特性決定地震數據無法約束隨機反演的高頻成分(高頻成分僅僅來自隨機模擬的數學實現)，因此高頻成分具有不確定性。

圖5 反演結果

為了進一步量化隨機反演的不確定性，提出不確定性指示公式

(16)

式(16)的計算結果表明，確定性反演協同約束的隨機反演的高頻成分不確定性僅為244364.55，遠小于常規隨機反演的434247.63。

7 實際資料試算

實際資料取自中國陸上M區塊的過A井的地震剖面，目的層巖性主要為碎屑巖，反演時窗從層位T1到T4，共3個層段(圖6)。首先根據解釋的4個層位構建層序地層網格，由于區內沉積相對平穩，地層厚度變化不大，因此3個層段的沉積模式均為平行頂、底。

圖6 過A井的地震剖面

A井數據未參與約束，僅用于估算波阻抗的期望、方差和垂向變差。建立波阻抗期望的方法與建立確定性反演低頻模型的方法一致，即對A井的波阻抗曲線沿層位橫向插值并由0～10Hz的低通濾波獲取(圖7)。根據擬合的高斯分布(圖8)求得T1-T2、T2-T3、T3-T4層段的波阻抗標準差，分別為867225、990082、898831kg·m-2·s-1。利用A井數據統計了T1-T2、T2-T3以及T3-T4層段的波阻抗垂向實驗變差(圖9)，變差函數模型由0.8倍的指數型變差函數與0.2倍的高斯型變差函數混合而成，得到T1-T2、T2-T3、T3-T4層段的垂向變程，分別為2.5、4.5、5.0ms。通過統計波阻抗確定性反演結果的層段計算均方根屬性，得到水平實驗變差(圖10)，變差函數模型與垂向變差一致，得到T1-T2、T2-T3、T3-T4層段的水平變程，分別為1000、1400、800m。

圖7 確定性反演的低頻模型

圖8 T1-T2(左)、T2-T3(中)、T3-T4(右)層段的波阻抗統計直方圖及擬合的高斯分布(紅線)

圖9 T1-T2(上)、T2-T3(中)以及T3-T4(下)層段的波阻抗垂向實驗變差(藍圈)及垂向擬合變差(紅線)

圖10 T1-T2(上)、T2-T3(中)以及T3-T4(下)層段的波阻抗水平實驗變差(藍圈)及水平擬合變差(紅線)

圖11為波阻抗確定性反演剖面，可見確定性反演結果更忠實于地震數據本身，其垂向分辨率與地震數據基本一致。圖12和圖13分別為隨機地震反演的實現A和實現B。由圖可見：①與確定性反演結果(圖11)相比，隨機地震反演的垂向分辨率更高，在水平變差以及層序地層網格的雙重約束下，橫向連續且自然，井旁道反演結果與實測曲線吻合度較高(圖12、圖13)。但是由于隨機地震反演結果的高頻成分來自隨機模擬和變差函數，當參與約束的測井數據較少時，其高頻成分的不確定性較大。②雖然兩次實現地震數據參與約束的程度一致，但是兩者的儲層整體特征差異較大(圖12、圖13)，并且與圖11的整體特征也有較大差異。

圖11 波阻抗確定性反演剖面

圖14和圖15分別為確定性反演協同約束的隨機地震反演的實現C和實現D。由圖可見：與圖12、圖13之間的差異相比，圖14和圖15的差異相對減小；圖14、圖15的整體特征更接近圖11，但前兩者的垂向分辨率更高。這是由于序貫Gibbs同位協同擾動模擬是借助同位協同簡單克里金而非簡單克里金確定隨機反演的先驗解空間。相比常規隨機反演，確定性反演協同約束的隨機地震反演的某個采樣點的波阻抗隨機反演結果不僅與鄰近采樣點的波阻抗值相關，還與當前采樣點的波阻抗確定性反演結果相關，在一定程度上增強了波阻抗隨機反演結果與確定性反演結果之間的相關性，從而降低隨機反演結果高頻成分的不確定性。

圖12 波阻抗隨機反演實現A

圖13 波阻抗隨機反演實現B

圖14 確定性反演協同約束的隨機地震反演實現C

圖15 確定性反演協同約束的隨機地震反演實現 D

圖16為隨機反演四個不同實現的信噪比，可見信噪比值約為10，說明地震數據對反演結果的約束程度基本相同，但這四個實現展現的儲層細節存在一定差異，從而驗證了地震反演具有多解性。

圖16 隨機反演四個不同實現的信噪比

8 結論

(1)相比確定性反演，隨機地震反演可以產生高分辨率的反演結果，其中垂向變差影響隨機反演的垂向分辨率，橫向變差影響隨機反演的橫向連續性。

(2)與均勻地震網格相比，由于層序地層網格融入了構造和沉積模式等信息，因此更適合變差函數橫向約束隨機反演；借助重采樣矩陣，在層序地層網格進行抽樣模擬，在均勻地震網格進行褶積正演，整個反演過程既滿足構造和沉積模式的約束，同時又符合地球物理原理。

(3)常規隨機反演的高頻成分來自隨機模擬、變差函數以及約束井數據等。當工區內參與約束的測井數據較少時，常規隨機反演高頻成分的不確定性較大。通過引入確定性反演的協同約束，可進一步限定候選解的解空間，增強波阻抗隨機反演結果與確定性反演結果的相關性，進而降低隨機反演高頻成分不確定性。實際數據試算表明，通過對比隨機反演的四個不同實現，驗證了所提算法的有效性。

附錄A 重采樣矩陣

(A-1)

(A-2)

因此，重采樣矩陣為Ns行、N列的稀疏矩陣

(A-3)

附錄B 同位協同簡單克里金計算克里金估值和克里金方差

同位協同簡單克里金是同位協同克里金與簡單克里金的結合[16]，其克里金權重由線性方程組

(B-1)

(B-2)

和

(B-3)