宗兆云，印興耀，張 峰，吳國忱

1 中國石油大學(xué)（華東）地球科學(xué)與技術(shù)學(xué)院，青島 266555

2 中國石油大學(xué)（北京）地球物理與信息工程學(xué)院，北京 102249

1 引 言

地震波縱橫波速度、巖石縱橫波模量、拉梅參數(shù)、體積模量及泊松比等彈性參數(shù)在油氣儲(chǔ)層預(yù)測和流體識別中發(fā)揮了重要作用.疊前地震反演是從疊前地震資料中獲取此類彈性參數(shù)的主要方法.由于正演算法不同，疊前地震反演可以分為基于射線追蹤和Zoeppritz方程及其近似的AVO反演和基于全波波動(dòng)方程的全波形反演.地震全波形反演不僅能夠利用地震波傳播的走時(shí)和振幅特性，還能夠自動(dòng)考慮地震波傳播過程中的反射、透射、散射、繞射、多次波、波型轉(zhuǎn)換等問題，而不局限于只利用反射波信息，充分利用地震波場的全波形特征，實(shí)現(xiàn)介質(zhì)彈性參數(shù)直接反演，但由于其計(jì)算效率和分辨能力有限，目前主要用于地震速度建模.相比全波形反演，AVO反演能夠綜合測井等多尺度資料的信息，其在儲(chǔ)層預(yù)測和流體識別方面的高效性及高分辨率特征，一直備受人們青睞[1-3］.

Zoeppritz方程是AVO反演的基礎(chǔ)，但由于該方程具有較強(qiáng)的非線性，人們更多地采用其各種近似方程進(jìn)行 AVO反演[4-5］.Aki和Richard[6］基于微擾理論和界面兩側(cè)介質(zhì)參數(shù)弱變化近似條件下推導(dǎo)得到基于縱橫波速度和密度的縱波反射系數(shù)近似方程.Shuey[7］建立了基于泊松比、橫波速度和密度的反射系數(shù)近似方程，該方程在含氣儲(chǔ)層AVO分析中發(fā)揮了重要作用，但由于該方程中泊松比項(xiàng)具有非線性特征，很少用于AVO反演.為提高參數(shù)反演穩(wěn)定性，F(xiàn)atti等[8］推導(dǎo)了基于縱橫波阻抗和密度的反射系數(shù)近似方程，該方程能夠更穩(wěn)定地反演得到縱橫波阻抗參數(shù).為直接獲取流體敏感參數(shù)，Goodway等[9］和 Russell等[10］分別推導(dǎo)了基于拉梅參數(shù)和流體因子的反射系數(shù)近似方程，并實(shí)現(xiàn)了從疊前地震資料中直接提取拉梅參數(shù)和流體因子.

基于反射系數(shù)近似方程，F(xiàn)atti等[8］用權(quán)疊加的方法得到了縱波阻抗的相對變化量和橫波阻抗的相對變化量，進(jìn)而計(jì)算得到絕對縱橫波阻抗.Simmons和Backus[11］應(yīng)用線性化的反演方法求取縱波阻抗、橫波阻抗和密度反射系數(shù)，反演過程中采用Gardner[12］的縱波速度與密度的經(jīng)驗(yàn)關(guān)系式和Castagna[13］的關(guān)于縱波與橫波速度的經(jīng)驗(yàn)方程.Buland和Omre[14］研究了基于貝葉斯理論的線性AVO反演方法，直接反演得到縱波速度、橫波速度和密度參數(shù)，但該方法沒有建立參數(shù)之間的關(guān)系.Yin等人[15］在考慮了參數(shù)去相關(guān)情況下提出了點(diǎn)約束稀疏脈沖疊前地震反演，直接反演得到縱橫波速度和密度參數(shù).通過對反射系數(shù)方程轉(zhuǎn)換到對數(shù)域，并進(jìn)行簡單推導(dǎo)，Connolly[16］首先提出彈性阻抗（EI）的概念，并有學(xué)者認(rèn)為EI在抗噪能力方面比疊前 AVO 反 演 有 優(yōu) 勢[17］.Yin 等[18］提 出 利 用Connolly彈性阻抗方程從三個(gè)角度反演結(jié)果中通過計(jì)算提取縱、橫波速度和密度參數(shù)的方法，進(jìn)而可以計(jì)算得到縱、橫波阻抗，拉梅常數(shù)，泊松比等豐富多樣的巖性參數(shù)，從而進(jìn)行地下儲(chǔ)層的展布情況及含油氣性預(yù)測.Wang等[19］針對Zoeppritz方程的Gray近似[20］作了彈性阻抗方程的推導(dǎo)，并實(shí)現(xiàn)了從EI反演中直接提取拉梅常數(shù)，降低了間接計(jì)算帶來的累積誤差.

隨著常規(guī)油氣資源獲取難度的增大，全球油氣勘探已進(jìn)入了常規(guī)與非常規(guī)油氣并重發(fā)展的時(shí)代.頁巖氣作為全球重要接替資源，其勘探開發(fā)技術(shù)成為目前人們的研究熱點(diǎn)和難點(diǎn).頁巖氣是指賦存于頁巖中的非常規(guī)氣.儲(chǔ)集空間以裂縫為主并可以吸附氣和水溶氣形式賦存，在裂縫發(fā)育帶可獲較高產(chǎn)量.頁巖儲(chǔ)層的基質(zhì)滲透率較低，需要裂縫才能形成工業(yè)產(chǎn)能.除本身的天然裂縫外，在開發(fā)過程中還應(yīng)考慮頁巖儲(chǔ)層在加砂壓裂改造時(shí)，是否易于被改造.脆性越大，越易被改造.地震巖石物理研究表明，楊氏模量和泊松比能夠較好地表征巖石的脆性，評價(jià)頁巖儲(chǔ)層的造縫能力，頁巖氣“甜點(diǎn)”具有高楊氏模量和低泊松比特征[21］.本文推導(dǎo)得到一種基于楊氏模量和泊松比的新的反射系數(shù)近似方程，并建立了一種穩(wěn)定獲取楊氏模量和泊松比的疊前地震直接反演方法.

2 YPD反射系數(shù)近似方程推導(dǎo)

1980年，Aki和Richards提出了基于縱橫波速度和密度的Zoeppritz近似公式為

其中，α、β分別為縱、橫波速度，ρ為密度，θ為入射角，Δα／α、Δβ／β、Δρ／ρ分別是縱、橫波速度反射系數(shù)和密度反射系數(shù).

由方程（1）出發(fā)，Zong等[23］建立反射系數(shù)與縱橫波模量和密度的關(guān)系式如下：

其中，M為縱波模量，與介質(zhì)抗壓縮性和硬度直接相關(guān)，體現(xiàn)儲(chǔ)層骨架和流體信息.μ為橫波模量，與介質(zhì)抗剪切性和剛度直接相關(guān)，體現(xiàn)儲(chǔ)層骨架信息[24］.ΔM／M、Δμ／μ分別為縱、橫波模量反射系數(shù).

在各向同性介質(zhì)中，縱、橫波模量與楊氏模量和泊松比關(guān)系為

令

隨著我國信息化技術(shù)的不斷發(fā)展，許多銀行在日常的運(yùn)營過程中也都已經(jīng)運(yùn)用計(jì)算機(jī)技術(shù)來形成了自身的數(shù)據(jù)庫管理系統(tǒng)，其對于銀行管理水平的提升也有著非常重要的意義。但是目前許多銀行在信息化技術(shù)的應(yīng)用以及建設(shè)過程中還存在有比較多的安全問題，并直接威脅到了該銀行的持續(xù)穩(wěn)定發(fā)展。針對這一問題，也就要求各銀行能夠在已有基礎(chǔ)上來對信息化管理系統(tǒng)以及數(shù)據(jù)庫系統(tǒng)進(jìn)行不斷的更新與優(yōu)化，來促使自身的信息管理能力得到進(jìn)一步的提升。只有在結(jié)合客戶實(shí)際需求的基礎(chǔ)上進(jìn)行信息的準(zhǔn)確把握跟及時(shí)決策，才能夠給用戶進(jìn)行全方位金融服務(wù)的有效提供。

則由方程（3）和方程（5）得

方程（7）兩邊同除以縱波模量可得

同理，由方程（4）和方程（6）得

令

由方程（5）、（6）和（10）可得

又

其中，σ1和σ2分別為上下層介質(zhì)泊松比.同理，

由于

方程（14）和（15）分別代入方程（12）和（13）可得

令

則

將方程（8）、（9）、（11）、（16）、（17）、（18）和（19）代入方程（2）可得

方程（20）建立了縱波反射系數(shù)與楊氏模量反射系數(shù)、泊松比反射系數(shù)及密度反射系數(shù)的線性關(guān)系.我們稱之為YPD近似方程.

3 YPD近似方程精度分析

利用Goodway等[9］根據(jù)實(shí)測資料給出的含氣砂巖模型，我們對利用YPD近似方程、精確Zoeppritz方程和Aki-Richards近似方程計(jì)算反射系數(shù)的精度進(jìn)行了分析，模型參數(shù)見表1.

表1 Goodway（1997）[9］三層含氣砂巖與頁巖模型Table 1 Model for gas-bearing sandstone and shale designed by Goodway（1997）[9］

入射介質(zhì)波阻抗小于透射介質(zhì)波阻抗的界面稱為正波阻抗界面，反之稱為負(fù)波阻抗界面.上述模型中，上覆頁巖、下伏含油氣砂巖的反射界面就是負(fù)波阻抗界面，而上覆砂巖、下伏頁巖的反射界面為正波阻抗界面.以上述模型為基礎(chǔ)，分別用精確的Zoeppritz方程、Aki-Richards近似方程、YPD近似方程計(jì)算了不同界面處的反射系數(shù)及近似方程與精確方程的殘差.圖1和圖2分別為正負(fù)波阻抗界面對比分析，圖1a、圖2a為分別用精確的Zoeppritz方程、Aki-Richards近似方程、YPD近似方程計(jì)算得到的反射系數(shù)隨入射角的變化，圖1b、圖2b為Aki-Richards近似方程、YPD近似方程計(jì)算得到反射系數(shù)與精確方程計(jì)算得到的反射系數(shù)的差值隨入射角的變化.由圖可知，基于YPD近似方程得到的反射系數(shù)在入射角為40°左右仍與精確Zoeppritz方程有較好的近似，能夠適用于較大角度入射情況下反射系數(shù)的準(zhǔn)確求解.

4 YPD-AVA疊前地震反演

以方程（20）為基礎(chǔ)，在貝葉斯反演框架下，假設(shè)待反演楊氏模量、泊松比及密度反射系數(shù)服從柯西分布，該分布假設(shè)可以最大限度提高反演分辨率[25-26］，假設(shè)似然函數(shù)服從高斯分布，同時(shí)在反演目標(biāo)函數(shù)中加入初始模型約束，并通過初始模型建立各道去相關(guān)矩陣，消除待反演參數(shù)間的互相關(guān)性，建立了一種YPD-AVA疊前地震反演方法，具體反演流程如圖3所示.

圖3 YPD-AVA疊前地震反演流程圖Fig.3 Flow chart for YPD-AVA pre-stack seismic inversion

貝葉斯框架下，待反演參數(shù)θ的后驗(yàn)概率分布為

其中：θ為待反演參數(shù)，x為觀測樣本，p（θ｜x）為后驗(yàn)概率密度，p（θ）為先驗(yàn)概率密度，p（x｜θ）為似然函數(shù).

貝葉斯估計(jì)由帶噪聲的疊前部分角度疊加道集數(shù)據(jù)D估計(jì)待反演參數(shù)（楊氏模量、泊松比與密度的反射系數(shù)）矩陣R′.假設(shè)似然函數(shù)服從高斯分布為

相比高斯分布，待反演參數(shù)的柯西分布在提高反演分辨能力等方面更有優(yōu)勢，因此，我們假設(shè)待反演參數(shù)服從柯西分布：

最大化后驗(yàn)概率分布（24）得到初始目標(biāo)函數(shù)為：

在目標(biāo)函數(shù)（25）的基礎(chǔ)上加入初始模型約束，目標(biāo)函數(shù)變?yōu)椋?/p>

其中：

λY，λPR和λρ分別為楊氏模量、泊松比和密度約束系數(shù)，P為

方程（26）具有弱非線性，可采用反復(fù)加權(quán)最小二乘（IRLS）[27］算法進(jìn)行優(yōu)化求解.

5 模型驗(yàn)證

為驗(yàn)證YPD-AVA疊前地震反演的可行性和穩(wěn)定性，本文采用實(shí)際單井模型進(jìn)行驗(yàn)證.模型測試中采用主頻為40Hz的雷克子波.圖4b中藍(lán)色曲線為某工區(qū)通過實(shí)際測井資料計(jì)算并通過深時(shí)轉(zhuǎn)換得到的時(shí)間域楊氏模量、泊松比和密度曲線，采用精確Zoeppritz方程進(jìn)行正演得到角度域疊前角度道集，然后利用本文提出的YPD-AVA疊前地震反演方法實(shí)現(xiàn)楊氏模量、泊松比和密度參數(shù)反演.圖4a和圖4b分別為無噪聲情況下角度域合成地震記錄及楊氏模量、泊松比和密度參數(shù)反演結(jié)果，其中紅色曲線為反演結(jié)果，反演結(jié)果表明，無噪聲情況下，該方法能夠獲取與真實(shí)值基本吻合的楊氏模量、泊松比和密度，驗(yàn)證了方法的可行性.為進(jìn)一步驗(yàn)證方法的穩(wěn)定性，對合成道集加不同程度的服從高斯分布的隨機(jī)噪聲，圖5—圖6依次是信噪比為5∶1和2∶1情況下的角度域合成地震記錄與反演結(jié)果，可以看出，在加噪情況下，楊氏模量和泊松比反演結(jié)果仍與真實(shí)值有較高的吻合度，但密度反演結(jié)果在局部有些不合理，這是因?yàn)槊芏软?xiàng)對反射系數(shù)值的影響較小，在使用AVO近似公式進(jìn)行反演時(shí)，很難在含噪聲情況下獲取合理的密度參數(shù)[28-29］.

6 實(shí)際資料試處理

實(shí)際資料來自中國陸上某勘探工區(qū)，進(jìn)行疊前地震反演之前，需要對地震數(shù)據(jù)進(jìn)行保幅處理，包括精細(xì)的波前擴(kuò)散補(bǔ)償、震源組合與檢波器組合效應(yīng)的校正、反Q濾波、地表一致性處理、疊前去噪處理、去除多次波等，并假設(shè)處理后的層間多次波、各向異性的影響可以忽略不計(jì)[30］.筆者利用本文提出的YPD-AVA疊前地震反演方法對不同角度疊加剖面進(jìn)行了實(shí)際資料試處理.本文反演方法輸入的是部分角度疊加剖面，首先根據(jù)已知的速度信息，將疊前CMP或CRP偏移距道集轉(zhuǎn)換到角度域[14］，然后，根據(jù)實(shí)際入射角范圍，進(jìn)行分角度疊加，提高資料信噪比.圖7a、7b和7c分別是小、中、大三個(gè)角度部分疊加剖面，圖中白色橢圓區(qū)域?yàn)闅獠匕l(fā)育區(qū).圖8為反演得到的楊氏模量剖面，在含氣區(qū)域呈高值異常顯示，圖9為反演得到的泊松比剖面，在含氣區(qū)域呈低值異常顯示，該反演結(jié)果與實(shí)際鉆遇結(jié)果吻合較好，驗(yàn)證了方法的有效性.

7 結(jié) 論

楊氏模量和泊松比能夠較好地表征巖石的脆性，評價(jià)頁巖儲(chǔ)層的造縫能力，頁巖氣“甜點(diǎn)”具有高楊氏模量和低泊松比特征，利用疊前地震反演獲取楊氏模量和泊松比參數(shù)成為利用疊前地震資料進(jìn)行頁巖氣“甜點(diǎn)”識別的重要手段.本文在平面縱波入射條件下，由Aki-Richards近似出發(fā)，推導(dǎo)得到基于楊氏模量、泊松比和密度的Zoeppritz近似公式（YPD近似方程），奠定了疊前反演獲取楊氏模量和泊松比的理論基礎(chǔ).在YPD近似方程基礎(chǔ)上建立了一種穩(wěn)定反演楊氏模量和泊松比的方法，模型測試表明，該方法在信噪比比較低的情況下仍能得到合理的反演結(jié)果，實(shí)際資料試處理驗(yàn)證了該方法在實(shí)際生產(chǎn)中的可行性，提供了一種從地震數(shù)據(jù)中直接提取楊氏模量和泊松比的疊前地震反演方法.與常規(guī)AVO方法類似，該方法建立在水平層狀假設(shè)基礎(chǔ)之上，且在中等入射角情況下具有較高精度.地震巖石物理研究表明，在實(shí)際頁巖氣地震識別中，除進(jìn)行地層巖石的脆性評價(jià)外，還需通過地層應(yīng)力參數(shù)等評價(jià)頁巖氣地層可壓裂性，因此下一步有必要建立地震反射系數(shù)與地層應(yīng)力等參數(shù)間的直接關(guān)系，并發(fā)展合理的疊前反演方法，以便更好地進(jìn)行頁巖氣藏地震識別與評價(jià).

（References）

[1]彭蘇萍，高云峰，楊瑞召等.AVO探測煤層瓦斯富集的理論探討和初步實(shí)踐——以淮南煤田為例.地球物理學(xué)報(bào)，2005，48（6）：1475-1486.Peng S P，Gao Y F，Yang R Z，et al.Theory and application of AVO for detection of coalbed methane——A case from the Huainan coalfield.Chinese J.Geophys.（in Chinese），2005，48（6）：1475-1486.

[2]牛濱華，文鵬飛，溫寧等.基于BSR的AVO正演估算水合物含量方法的研究.地球物理學(xué)報(bào)，2006，49（1）：143-152.Niu B H，Wen P F，Wen N，et al.Estimation of hydrate concentration based on AVO modeling of BSR.Chinese J.Geophys.（in Chinese），2006，49（1）：143-152.

[3]彭真明，李亞林，魏文閣等.粒子濾波非線性AVO反演方法.地球物理學(xué)報(bào)，2008，51（4）：1218-1225.Peng Z M，Li Y L，Wei W G，et al.Nonlinear AVO inversion using particle filter.Chinese J.Geophys.（in Chinese），2008，51（4）：1218-1225.

[4]李景葉，陳小宏，郝振江等.多波時(shí)移地震AVO反演研究.地球物理學(xué)報(bào)，2005，48（4）：902-908.Li J Y，Chen X H，Hao Z J，et al.A study on multiple timelapse seismic AVO inversion.Chinese J.Geophys.（in Chinese），2005，48（4）：902-908.

[5]陳建江，印興耀.基于貝葉斯理論的AVO三參數(shù)波形反演.地球物理學(xué)報(bào)，2007，50（4）：1251-1260.Chen J J，Yin X Y.Three-parameter AVO waveform inversion based on Bayesian theorem.Chinese J.Geophys.（in Chinese），2007，50（4）：1251-1260.

[6]Aki K，Richards P G.Quantitative Seismology（2ed）.San Francisco：W.H.Freeman，1980：158-164.

[7]Shuey R T.A simplification of the Zoeppritz equations.Geophysics，1985，50（4）：609-614.

[8]Fatti J L，Smith G C，Vail P J，et al.Detection of gas in sandstone reservoirs using AVO analysis：A 3-D seismic case history using the Geostack technique.Geophysics，1994，59（9）：1362-1376.

[9]Bill G，Chen T W，Downton J.Improved AVO fluid detection and lithology discrimination using Lame′s petrophysical parameters：“λρ”，“μρ”，&“λ／μfluid stack”，from P and S inversions.SEG Technical Program Expanded Abstracts，1997，16（1）：183-186.

[10]Russell B H，Gray D，Hampson D P.Linearized AVO and poroelasticity.Geophysics，2011，76（3）：C19-C29.

[11]Simmons J L，Backus M M.Waveform-based AVO inversion and AVO prediction-error.Geophysics，1996，61（6）：1575-1588.

[12]Gardner G H F，Gardner L W G，Gregory A R.Formation velocity and density-the diagnostic basics for stratigraphic traps.Geophysics，1974，39（6）：770-780.

[13]Castagna J，Batzle M，Eastwood R L.Relationships between compressional-wave and shear-wave velocities in clastic silicate rocks.Geophysics，1985，50（4）：571-581.

[14]Arild B， More H.Bayesian linearized AVO inversion.Geophysics，2003，68（1）：185-198.

[15]Yin X Y，Yang P J，Zhang G Z.A novel prestack AVO inversion and its application.SEG Technical Program Expanded Abstracts，2008，27（1）：2041-2045.

[16]Connolly P.Elastic impedance.The Leading Edge，1999，18（4）：438-452.

[17]Guillaume C.AVO inversion and elastic impedance.SEG Technical Program Expanded Abstracts，2000，19（1）：142-145.

[18]Yin X Y，Yuan S H，Zhang F C.Rock elastic parameters calculated from elastic impedance.CPS／SEG Technical Program Expanded Abstracts，2004：538-542.

[19]Wang B L，Yin X Y，Zhang F C.Laméparameters inversion based on elastic impedance and its application.Applied Geophysics，2006，3（3）：174-178.

[20]Gray D， Verttas D G C.Elastic inversion for Lame parameters.SEG Technical Program Expanded Abstracts，2002，21（1）：213-216.

[21]Sena A，Gabino C，Chesser K，et al.Seismic reservoir characterization in resource shale plays： “Sweet spot”discrimination and optimization of horizontal well placement.SEG Technical Program Expanded Abstracts，2011，30（1）：1744-1748.

[22]Goodway B，Perez M，Varsek J，et al.Seismic petrophysics and isotropic-anisotropic AVO methods for unconventional gas exploration.The Leading Edge，2010，29（12）：1500-1508.

[23]Zong Z Y，Yin X Y，Wu G C.Robust AVO inversion for elastic modulus and its application in fluid factor calculation.SEG Technical Program Expanded Abstracts，2011，30（1）：406-411.

[24]Mavko G， Mukerji T，Dvorkin J.The Rock Physics Handbook：Tools for Seismic Analysis in Porous Media.Cambridge：Cambridge University Press，2009：21-25.

[25]Wubshet A，Sacchi M D.High-resolution three-term AVO inversion by means of a Trivariate Cauchy probability distribution.Geophysics，2011，76（3）：R43-R55.

[26]楊培杰，印興耀.非線性二次規(guī)劃貝葉斯疊前反演.地球物理學(xué)報(bào)，2008，51（6）：1876-1882.Yang P J，Yin X Y.Non-linear quadratic programming Bayesian prestack inversion.Chinese J.Geophys.（in Chinese），2008，51（6）：1876-1882.

[27]Daubechies I，Devore R，F(xiàn)ornasier M，et al.Iteratively reweighted least squares minimization for sparse recovery.Communications on Pure and Applied Mathematics，2010，63（1）：1-38.

[28]Lines L R.Density contrast is difficult to determine from AVO.Canadian Journal of Exploration Geophysics，1998.

[29]Downton J E.Seismic parameter estimation from AVO inversion[Ph.D.thesis］.Canada：University of Calgary，2005：87.

[30]宗兆云，印興耀，吳國忱.基于疊前地震縱橫波模量直接反演的流體檢測方法.地球物理學(xué)報(bào)，2012，55（1）：284-292.Zong Z Y，Yin X Y，Wu G C.Fluid identification method based on compressional and shear modulus direct inversion.Chinese J.Geophys.（in Chinese），2012，55（1）：284-292.