顧 雯 印興耀 巫芙蓉 李 坤 翟浩杰

(①中國石油大學(華東)地球科學與技術(shù)學院，山東青島 266555；②東方地球物理公司研究院，河北涿州 072750；③北京中恒利華石油技術(shù)研究所,北京 100102)

0 引言

從前人研究成果看，薄儲層預(yù)測從地質(zhì)綜合研究(沉積環(huán)境、成巖作用、分布規(guī)律等)和疊后定量研究逐步走向疊前多參數(shù)定量研究，地震反演技術(shù)是薄儲層定量預(yù)測的最主要方法之一[2]。高分辨率疊前地震反演是定量評價頁巖氣層彈性、物性及有利甜點區(qū)空間展布范圍的重要方法。基于統(tǒng)計學理論的隨機反演方法或概率化反演算法在獲取高分辨率模型參數(shù)最優(yōu)解的同時，可以有效地評價反演結(jié)果和先驗信息的不確定性，能更好地定量解釋儲層參數(shù)。以貝葉斯反演框架為例，業(yè)界針對模型參數(shù)先驗信息構(gòu)建、后驗概率密度分布構(gòu)建、隨機采樣算法優(yōu)選進行了大量研究。Eide等[3]從概率角度出發(fā)，假設(shè)先驗概率密度函數(shù)以及和井震有關(guān)的條件似然函數(shù)均符合高斯分布，基于貝葉斯公式推導(dǎo)出與儲層參數(shù)有關(guān)的后驗概率密度函數(shù)，通過對后驗概率密度函數(shù)多次采樣最終得到反演結(jié)果。Escobar等[4]基于貝葉斯理論和Zoeppritz近似方程求取縱橫波阻抗的后驗分布，用新的序貫算法對后驗概率進行抽樣，最終提出一種層序網(wǎng)格內(nèi)的超快速疊前隨機反演方法。為了加快隨機反演計算速度，Contreras等[5]引入馬爾科夫鏈—蒙特卡洛(MCMC)技術(shù)，通過聯(lián)合模擬縱波和橫波阻抗將其擴展為對多個部分疊加道集的反演，將隨機反演推向疊前。Smith等[6]利用MCMC方法獲取模型參數(shù)的后驗分布，利用MCMC方法逼近貝葉斯統(tǒng)計分析中出現(xiàn)的積分運算。大量文獻研究成果表明，基于協(xié)克里金、序貫高斯、模擬退火等地質(zhì)統(tǒng)計學算法的隨機反演能夠隨機模擬薄儲層[7-11]。近年來，基于地質(zhì)統(tǒng)計學反演的地震波形指示反演方興未艾[12-15]，由于地震波形的變化反映了沉積環(huán)境和巖性組合的空間變化，將這種變化定義為儲層構(gòu)型，可以利用波形變化間接推斷儲層空間的相變特征。該技術(shù)在空間結(jié)構(gòu)化數(shù)據(jù)指導(dǎo)下不斷尋優(yōu)，利用地震波形的橫向變化和空間距離代替變差函數(shù)進行樣本選樣，使反演結(jié)果在橫向上更符合地質(zhì)沉積規(guī)律。

基于疊后地震波形指示反演在薄儲層預(yù)測方面的優(yōu)勢，以及疊前彈性參數(shù)較疊后彈性參數(shù)信息更豐富、對儲層的敏感性更高的特點，本文基于巖石物理資料，應(yīng)用地震驅(qū)動+儲層構(gòu)型約束的高精度疊前隨機反演方法，探尋定量表征優(yōu)質(zhì)薄頁巖的技術(shù)，為深層頁巖氣地質(zhì)甜點預(yù)測提供技術(shù)支撐。

1 方法原理

常規(guī)疊前反演主要通過巖石物理建模技術(shù)計算彈性參數(shù)建立參數(shù)模型，再根據(jù)AVO響應(yīng)特征建立道集與彈性參數(shù)的關(guān)系，利用Zoeppritz近似方程或者彈性阻抗方程建立目標函數(shù)迭代求解參數(shù)模型。針對深層優(yōu)質(zhì)頁巖定量預(yù)測技術(shù)難題，本文在疊后波形指示反演基礎(chǔ)上，充分利用疊前多參數(shù)相似特征和巖相組合的指示作用，利用波形相似性和空間距離多參數(shù)優(yōu)選模擬樣本，以建立樣本空間(圖1)。疊后波形指示反演根據(jù)波形相似性優(yōu)選統(tǒng)計樣本，對比預(yù)測道地震波形與所有已知井旁道地震波形，優(yōu)選最相似的若干井樣本，再比較這些井的不同頻段的濾波曲線，尋找共性結(jié)構(gòu)特征并建立初始模型[8]。在疊前，該方法利用道集波形和AVO特征，基于道集波形相似性、AVO特征和空間距離的三變量優(yōu)選方法提取結(jié)構(gòu)相似的井數(shù)據(jù)作為空間估值樣本，建立待判別道集初始模型，然后在貝葉斯框架下開展MCMC模擬，補充確定性疊前反演缺失的高頻成分。該模擬算法的核心是根據(jù)實際概率分布得到統(tǒng)計意義上正確的隨機樣點分布，實現(xiàn)全局優(yōu)化的多個等概率模擬結(jié)果[16-17]。由鉆/測井數(shù)據(jù)、巖相地質(zhì)數(shù)據(jù)及已有的確定性反演數(shù)據(jù),結(jié)合地層框架模型建立的目的層概率密度函數(shù)和縱橫向變差函數(shù)獲得地質(zhì)統(tǒng)計學信息，降低了MCMC模擬的不確定性，在反演過程中加入地震數(shù)據(jù)約束，進一步降低了地質(zhì)統(tǒng)計學反演的不確定性，可保證等概率結(jié)果的相似性。

圖1 疊前波形特征選樣分析示意圖(AVO屬性平面圖)vP、vS、Den分別代表縱波速度、橫波速度、密度

1.1 疊前波形指示樣本集的建立

首先，通過疊前道集特征矩陣建立測井曲線樣本數(shù)據(jù)集，其中疊前道集特征矩陣是由道集波形特征動態(tài)矩陣、AVO特征指示因子矩陣和距離權(quán)重因子矩陣組成，描述了待預(yù)測點與樣本點的井旁道波形相似度、AVO特征相似度和空間距離等3個井震相關(guān)性特征。道集波形特征動態(tài)矩陣由對比預(yù)測點與各樣本點的近、中、遠分角度疊加地震波形獲得，AVO特征指示因子矩陣由對比預(yù)測點與各樣本點的AVO特征獲得，距離權(quán)重因子矩陣由對比預(yù)測點與各樣本點的距離獲得。其次，在小波域進行頻率分解，提取樣本集曲線的共性結(jié)構(gòu)作為初始模型。對于不同的共性結(jié)構(gòu)相關(guān)截止頻率l應(yīng)用

(1)

在道集優(yōu)化的基礎(chǔ)上計算道集波形特征參數(shù)，利用奇異值分解進行波形分類，建立地震道集波形與不同結(jié)構(gòu)特征的井上彈性參數(shù)曲線樣本的映射關(guān)系。首先建立初始樣本集

D=UΣV*

(2)

式中：U和V分別為n×n階地震波形數(shù)據(jù)正交矩陣和m×m階井點屬性正交矩陣，V*為V的共軛轉(zhuǎn)置；Σ為n×m階非負實數(shù)對角矩陣，表征井上曲線樣本與地震波形的相關(guān)性。定義U為疊前道集矩陣，則

U=[αβγ]

(3)

式中：α為采集數(shù)據(jù)波形特征參數(shù)；β為AVO特征因子參數(shù)；γ為距離加權(quán)因子參數(shù)。

1.2 基于貝葉斯框架的隨機反演

針對反演優(yōu)化問題，本文在貝葉斯框架下構(gòu)建了待反演模型參數(shù)的后驗概率密度分布(PDF)及目標泛函。貝葉斯推斷為在當前觀測地震數(shù)據(jù)的情況下，結(jié)合觀測數(shù)據(jù)之前對待反演模型參數(shù)的先驗知識計算模型參數(shù)PDF的統(tǒng)計學習方法，是一種將先驗信息和似然函數(shù)轉(zhuǎn)化為后驗信息的方法，其數(shù)學基礎(chǔ)是貝葉斯定理[18-19]。假設(shè)未知模型參數(shù)服從某種類型的概率密度分布，即所謂的先驗知識，它的引入是對帶限地震數(shù)據(jù)的修正，貝葉斯公式為

(4)

式中：p(m|d)為待反演模型參數(shù)的PDF；p(d|m)為觀測地震數(shù)據(jù)的似然函數(shù)；p(m)為待反演模型參數(shù)的先驗分布。貝葉斯反演是從已知數(shù)據(jù)d中估計彈性參數(shù)模型m的PDF的過程。假設(shè)觀測地震數(shù)據(jù)的噪聲分布服從高斯概率密度分布

(5)

式中：G為地震數(shù)據(jù)的正演矩陣；Cd為彈性阻抗數(shù)據(jù)體的協(xié)方差矩陣。通過測井數(shù)據(jù)統(tǒng)計分析待反演模型參數(shù)的先驗信息，得到模型參數(shù)的先驗概率密度分布

(6)

式中：Cm為模型參數(shù)的協(xié)方差矩陣；μm為模型參數(shù)的先驗均值。將式(6)代入貝葉斯公式，可以求取模型參數(shù)的p(m|d)

(7)

其中

對式(7)兩邊取對數(shù)，得到等價的目標泛函

(8)

令lnO(m|d)=0，得到反演方程

(9)

利用式(9)可以求出貝葉斯推斷中的最大后驗概率解，可以采用迭代重加權(quán)最小二乘算法求解，也可以利用蒙特卡洛仿真模擬和馬爾科夫鏈模型獲取模型參數(shù)的隨機解。由于基于統(tǒng)計學和隨機采樣的反演算法可以有效提高反演分辨率，因此本文利用貝葉斯概率解評價模型參數(shù)的最優(yōu)解和不確定性。此外，本文基于疊前高精度反演結(jié)果，通過建立道集與彈性參數(shù)曲線的樣本模型尋找待預(yù)測區(qū)與已知樣本井數(shù)據(jù)的共性結(jié)構(gòu)(截距、梯度、振幅、相位等)，在相似共性結(jié)構(gòu)的道集組合約束下獲得高分辨率的反演結(jié)果(圖2)。

圖2 疊前多參數(shù)相似性高精度反演流程

MCMC采用的算法為Metropolis-Hasting算法，通過該算法不斷對由波形優(yōu)選樣本建立的地震波形指示的初始模型進行隨機擾動，從而得到貝葉斯推斷中的最大后驗概率解，取多次可行實現(xiàn)的均值作為期望值輸出，具體步驟如下。

(3)產(chǎn)生隨機數(shù)μ，其中μ～U(0,1)；

(10)

(6)令i=i+1，返回步驟(2)、迭代直至平穩(wěn)狀態(tài)(總迭代次數(shù)為N)。

可以證明以此種方式構(gòu)建的接受概率是滿足細致平衡條件的，算法實現(xiàn)流程見圖3。

圖3 算法流程

1.3 地質(zhì)甜點參數(shù)(TOC)模擬

由于無法直接反演頁巖氣甜點參數(shù)(TOC、含氣量、孔隙度等物性參數(shù))，目前主要基于單參數(shù)或多參數(shù)回歸分析建立預(yù)測模型預(yù)測甜點參數(shù)。這種方法往往預(yù)測精度較低，且多解性強，對于指導(dǎo)頁巖氣評層選區(qū)意義不大。波形的變化特征能有效指示巖相組合的彈性參數(shù)變化特征，利用其波形特征優(yōu)選模擬樣本建立樣本空間(也稱波形變差函數(shù))，在貝葉斯框架下開展隨機模擬，補充確定性疊前反演縱橫波速度比時缺失的高頻成分[20-21]。

地質(zhì)統(tǒng)計學由區(qū)域變量抽取的樣本值估算方差，通過變程反映變量的影響范圍。經(jīng)典變差函數(shù)中滯后距h代表區(qū)域空間兩點之間的距離，變差函數(shù)值γ(h)代表區(qū)域空間兩點參數(shù)值的差(圖4)。經(jīng)典變差函數(shù)描述的是參數(shù)隨實際區(qū)域空間距離而變化的規(guī)律。假設(shè)有樣本A、B，則滯后距h=|AB|。

圖4 變差函數(shù)相關(guān)性分析示意圖

對于一維區(qū)域空間，即A=(x1)，B=(x2)，則

|AB|=|x1-x2|

(11)

對于二維區(qū)域空間，即A=(x1，y1)，B=(x2，y2)，則

(12)

對于三維區(qū)域空間，即A=(x1，y1，z1)，B=(x2，y2，z2)，則

(13)

通過計算波形空間兩點之間的距離得到變差函數(shù)滯后距h。假設(shè)波形空間是四維空間，波形空間中的點自然是四維點。設(shè)有兩點的地震波形，質(zhì)心、均值、方差、變方差分別為C、A、V、Vv，則

(14)

式中下標1、2表示地震波序號。甜點參數(shù)模擬思想基于子波不變假設(shè)，認為地震彈性參數(shù)差異與井特征參數(shù)結(jié)果變化密切相關(guān)，即地震彈性參數(shù)隨井位而變化。根據(jù)質(zhì)心、均值、方差、變方差等四種屬性描述地震特征參數(shù)差異的特征向量，利用測井參數(shù)的統(tǒng)計變量與井的特征參數(shù)垂向變差函數(shù)表征井垂向結(jié)構(gòu)變化相對地震差異的變化量。最后統(tǒng)計預(yù)測道地震特征參數(shù)的特征向量，利用波形變差函數(shù)模擬預(yù)測道(井)的特征參數(shù)實現(xiàn)特征參數(shù)模擬。該模擬方法有效結(jié)合地震、地質(zhì)和測井信息，利用地震信息指導(dǎo)井參數(shù)高頻模擬，利用波形橫向變化與井特征參數(shù)相對變化關(guān)系，建立特征向量變差函數(shù)模型，實現(xiàn)井震協(xié)同高頻模擬。

2 模型驗證

為了驗證所提方法的有效性和對薄層的縱、橫向識別能力，首先對比薄砂層與厚砂層的疊前道集特征，選取A1、A4井開展正演分析。由正演模擬的近、中、遠道集(圖5)可見：道集波形特征發(fā)生明顯變化；兩口井在零炮檢距處均表現(xiàn)為強振幅、正極性AVO特征，振幅呈減小趨勢，但隨著入射角增大，兩口井AVO特征逐漸發(fā)生變化，且遠道集更能體現(xiàn)波形的差異性(圖6)。

圖5 A1井(上)和A4井(下)正演道集A1井同時存在薄砂層與厚砂層，A4井只存在厚砂層。結(jié)論中藍色代表泥巖，黃色代表砂巖

圖6 A1和A4井AVO響應(yīng)特性

進一步設(shè)計了砂泥巖薄互層地質(zhì)模型(圖7)。為了便于開展波形相似性反演，在該模型上建立5口虛擬井(A1～A5)提取不同的儲層特征。利用主頻為40～100Hz的零相位雷克子波作為震源進行正演模擬，選取低頻(40Hz)、中頻(70Hz)和高頻(100Hz)地震正演剖面開展不同反演方法實驗。

圖7 砂泥巖薄互層地質(zhì)模型

圖8為40Hz、70Hz和100Hz地震正演、稀疏脈沖反演和多參數(shù)波形相似性反演剖面。由圖可見：

①40Hz地震正演剖面由于分辨率較低，完全無法識別薄互層砂體，由于砂體組合不同，地震波形差異也較大(圖8a上)；70Hz地震正演剖面只能識別疊置砂體(圖8b上)；100Hz地震正演剖面可以識別每一個薄砂體(圖8c上)。②稀疏脈沖反演剖面與地震正演剖面的分辨率特征相似，如40Hz稀疏脈沖反演剖面完全無法識別薄互層砂體(圖8a中)，70Hz稀疏脈沖反演剖面可以分辨疊置的兩套砂體(圖8b中)，100Hz稀疏脈沖反演剖面可以分辨第3套薄砂體(圖8c中)。③多參數(shù)波形相似性反演剖面的分辨率明顯高于稀疏脈沖反演(圖8a下、圖8b下、圖8c下)，并且在70Hz時就能完全識別薄砂體(圖8b下)。

圖8 40Hz(a)、70 Hz(b)和100Hz(c)地震正演 (上)、稀疏脈沖反演(中)和多參數(shù)波形相似性反演(下)剖面

泥巖背景中由上至下發(fā)育4組薄互層砂體，其中前兩套砂體為疊置砂體，第三套砂體最薄(厚度為4m)，砂巖、泥巖速度分別為3500、2800 m/s，密度分別為2.65、2.26g/cm3

3 應(yīng)用實例

3.1 研究區(qū)地質(zhì)背景

3.2 疊前道集優(yōu)化處理

多參數(shù)波形相似性反演技術(shù)的關(guān)鍵是道集優(yōu)化處理、敏感參數(shù)體生成、樣本選取、截止頻率確定等流程?；谛律傻寞B前敏感參數(shù)體，再利用其波形相似性和空間距離多參數(shù)優(yōu)選模擬樣本，建立波形與彈性參數(shù)曲線的樣本模型，在貝葉斯框架下開展隨機模擬，補充確定性疊前反演缺失的高頻成分，尋找待預(yù)測區(qū)與已知樣本的共性結(jié)構(gòu)(振幅、頻率、相位等)，在相似共性結(jié)構(gòu)的道集組合約束下反演，從而獲得高分辨率的薄層反演成果。

研究區(qū)疊前時間偏移CRP道集上存在明顯的隨機噪聲和遠炮檢距剩余時差，前者影響疊前道集的信噪比，后者導(dǎo)致道集不平、AVO規(guī)律不明確。由于疊前高精度反演對道集質(zhì)量要求較高，為此制定了道集優(yōu)化流程。首先，在利用射線追蹤方法將炮檢距轉(zhuǎn)換為入射角時需要質(zhì)控和優(yōu)化地震速度，使其滿足射線追蹤的需求；其次，在道集優(yōu)化中需考慮隨機噪聲壓制與多次波問題；最后，在去除噪聲的道集上，針對遠道進行剩余時差校正。圖9為拉平前、后道集。由圖可見：①由于道集中存在較多的隨機噪聲且遠炮檢距道集不平，尤其在遠炮檢距處，原始道集與合成記錄的相關(guān)性較差，O3w的AVO規(guī)律也較差(圖9左)。②利用Radon變換和移動積分道集拉平技術(shù)優(yōu)化處理后的道集資料明顯提高了信噪比，同時大幅提高了疊前道集的質(zhì)量(圖9右)，為后續(xù)高精度反演提供了可靠的道集資料；利用優(yōu)化后的道集進一步開展井震標定，優(yōu)化后的道集與合成記錄的相關(guān)性得到明顯改善，同時O3w底部AVO曲線與井上正演道集的AVO曲線更吻合。

圖9 拉平前(左)、后(右)道集T3x為須家河組，T1f為飛仙關(guān)組

3.3 優(yōu)質(zhì)頁巖地質(zhì)甜點預(yù)測關(guān)鍵參數(shù)

3.3.1 巖石物理分析

由巖石物理敏感參數(shù)分析直方圖(圖10)看出，密度能夠較好地篩選優(yōu)質(zhì)頁巖，因此將密度作為敏感參數(shù)。進一步將頁巖分為Ⅰ、Ⅱ、Ⅲ類：密度小于2.57g/cm3的為Ⅰ類(優(yōu)質(zhì)頁巖)；密度為2.57～2.65g/cm3的為Ⅱ類(優(yōu)質(zhì)頁巖)；密度大于2.65g/cm3的為Ⅲ類。對比地質(zhì)統(tǒng)計學反演與多參數(shù)波形相似性反演的密度剖面(圖11)可見：地質(zhì)統(tǒng)計學反演結(jié)果的橫向變化不穩(wěn)定，不符合頁巖分布地質(zhì)規(guī)律，無法有效識別Ⅰ類頁巖；多參數(shù)波形相似性反演結(jié)果的縱向分辨率高，反演信息更豐富，且能有效區(qū)分Ⅰ、Ⅱ類頁巖。

圖10 龍馬溪組優(yōu)質(zhì)頁巖密度分析直方圖

圖11 地質(zhì)統(tǒng)計學反演 (上)與多參數(shù)波形相似性反演(下)的密度剖面P2l為龍?zhí)督M，P1l為梁山組，S1l2為龍二段

3.3.2 正則化因子求取

正則化可理解為對帶限觀測數(shù)據(jù)的補充，通過引入正則化乘法項，將不適定的地震反演問題轉(zhuǎn)化為適定的反問題求解，以獲取待反演模型參數(shù)的最優(yōu)解或近似解[22]。為了均衡預(yù)測樣本井點與非樣本井點彈性參數(shù)的相關(guān)性，需求取樣本井點的彈性阻抗與預(yù)測地震中頻彈性阻抗的直接匹配系數(shù)，因而設(shè)置正則化因子參數(shù)。通過正則化因子加強約束，進行地震趨勢(振幅、頻率、相位)控制，防止過度擬合。引入彈性阻抗模型作為先驗信息，可以修正反演結(jié)果，能夠滿足樣本要求而且穩(wěn)定。目標函數(shù)為

(15)

3.3.3 TOC預(yù)測

本文通過特征參數(shù)模擬直接預(yù)測頁巖氣的地質(zhì)甜點參數(shù)，提高了甜點參數(shù)的預(yù)測精度。圖12為H202井測井綜合解釋圖。由圖可見，優(yōu)質(zhì)頁巖的測井響應(yīng)表現(xiàn)為低密度、高GR、高TOC、高含氣量、高孔隙度、中高脆性等特征。圖13為過L101-H203-H202H3井TOC模擬圖。由圖可見，TOC模擬結(jié)果與測井曲線吻合度較高，橫向變化趨勢合理，其中向斜平緩區(qū)TOC相對較高。以高品質(zhì)OVT域處理的三維地震資料為基礎(chǔ)，基于巖石物理分析，將彈性參數(shù)轉(zhuǎn)化為儲層指標，利用疊前波形指示反演預(yù)測優(yōu)質(zhì)頁巖厚度和基于波形變差函數(shù)模擬方法直接預(yù)測TOC甜點參數(shù)(表1)，能準確識別水平井的油藏目標，為水平井布署提供技術(shù)支撐。

圖12 H202井測井綜合解釋圖

圖13 過L101-H203-H202H3井TOC模擬圖

表優(yōu)質(zhì)頁巖預(yù)測結(jié)果對井分析表

4 結(jié)論