馬琦琦，孫贊東

中國石油大學(北京)地球物理學院，北京 102249

0 引言

地震波數(shù)據(jù)中蘊含了豐富的有助于揭示地下巖層及其孔隙流體特征的信息。隨著勘探難度的加大以及勘探技術的進步，基于地震資料的地下巖性探測以及孔隙流體識別已成為現(xiàn)階段油氣勘探的關鍵技術[1]。在高品質(zhì)的地震數(shù)據(jù)基礎上，儲層、流體識別的精度主要依賴于選取的彈性參數(shù)對于巖性和流體的敏感度以及提取彈性參數(shù)方法的可靠性[2-4]。

國內(nèi)外學者利用彈性參數(shù)對儲集層識別及流體檢測方面做了大量的研究。Gregory和Domenico[5-6]提出了利用縱橫波速度比VpVs識別流體；Fatti等[7]利用泥巖基線進行了烴類流體的檢測；Goodway等[8-9]闡述了拉梅參數(shù)λ、μ和密度ρ的組合(例如，λμ、λρ、μρ)可以有效地識別儲層和流體；Russell等[10]結合了孔隙彈性理論提出了Russell流體因子用來識別流體；印興耀等[11]對基于Russell流體因子F和縱波阻抗Ip的兩項彈性阻抗反演進行了研究；Du等[12]指出Russell流體因子可以有效識別烴類流體。

疊前數(shù)據(jù)相對于疊后數(shù)據(jù)，包含更多的地下介質(zhì)信息[13]。疊前AVA/AVO(振幅隨入射角度/偏移距變化)反演技術是提取隱藏在疊前地震數(shù)據(jù)中的巖性、流體信息的重要途徑[14-15]。許多學者基于縱波反射地震數(shù)據(jù)進行了疊前反演的研究，然而縱波反射系數(shù)對橫波信息和密度是不夠敏感的[16]。因此，由于擁有更豐富的彈性信息，基于縱波和橫波信息的疊前聯(lián)合反演相對僅采用縱波數(shù)據(jù)的疊前縱波反演，可以獲得更加優(yōu)質(zhì)的反演結果[17-19]。Stewart[20]首先提出了疊前聯(lián)合反演方法，并利用縱、橫波數(shù)據(jù)提取了縱波速度反射系數(shù) ΔVpVp和橫波速度反射系數(shù) ΔVsVs；Buland和Omre[21]利用馬爾科夫鏈的方法實現(xiàn)了聯(lián)合反演；Viere和Landro[22]提出利用最小二乘估計算法求解PP和PS反射系數(shù)的近似線性表達式；Rabben等[23]利用Metropolis-Hastings算法實現(xiàn)了聯(lián)合反演；張廣智等[24]推導了基于脆性因子、泊松比和密度的近似方程，并利用聯(lián)合反演提取了以上彈性參數(shù)；Lu等[25]采用泰勒展開和均值漂移法提高了疊前聯(lián)合反演結果的準確性；杜炳毅等[26]推導了基于Gassmann流體項、剪切模量和密度的轉(zhuǎn)換波公式，實現(xiàn)了基于這三個參數(shù)的聯(lián)合反演；王彥飛等[16]提出了基于粒子濾波提供先驗信息的L1范數(shù)約束的算法來進行聯(lián)合反演，進一步提高了反演的精度。

在前人研究的基礎上，本文提出了基于雙約束的流體因子F、巖性因子μρ和密度ρ的聯(lián)合反演方法。首先推導了基于F、μρ和ρ的縱波和轉(zhuǎn)換波反射系數(shù)方程，并參考致密砂巖儲層建立正演參數(shù)模型對新推導的公式進行精度分析，為實現(xiàn)F、μρ的直接提取奠定了基礎，彈性參數(shù)的直接提取可以有效避免間接計算引入的累積誤差。在實踐中，由于地震資料質(zhì)量不高、子波提取不準等問題，反演本身是一個不適定的問題，雖然PP-PS聯(lián)合反演在一定程度上降低了解的不適定性，但是仍有必要在目標函數(shù)中增加適當?shù)恼齽t化約束來提高反演的穩(wěn)定性。因此，本文借助貝葉斯理論引入模型參數(shù)的先驗分布構建正則化項，以提高反演的穩(wěn)定性，并在目標函數(shù)中加入改進的低頻約束項，不僅得到了可靠的低頻信息，還提高了反演的魯棒性。由于反演的目標函數(shù)是有關模型參數(shù)的非線性公式，通常假定背景縱橫波速度比為常數(shù)，而這種方法往往會降低反演精度，本文引入了重加權迭代算法來求解方程，在不斷的迭代過程中更新背景值，以提高反演結果的精度。最后，將提出的新方法應用到了模型資料和實際資料中。對比結果表明該方法具有一定的可行性和實用性，提供了一種可靠的流體和巖性因子的直接提取方法。

1 方法原理

隨著勘探程度的不斷加深，非常規(guī)油氣勘探成為了全球油氣資源勘探的重要組成部分，其中，致密儲層的勘探是非常規(guī)油氣領域研究的熱點問題之一[27]。Zhou和Hilterman[28]通過對大量實測數(shù)據(jù)分析，認為參數(shù)Russell流體因子F對固結碎屑巖儲層的適用性最高，可以滿足對具有低孔、低滲的特點、固結程度較高的致密儲層孔隙流體的敏感指示。不同的流體具有不同的F值，因此根據(jù)流體項F可以區(qū)分干燥和飽和巖石以及飽和巖石中的流體類型；剪切模量μ代表了儲層骨架的特性，表征了巖石的巖性特征，密度ρ也能顯示儲層的巖性特征，Goodway[8]提出兩者的乘積μρ可以做為優(yōu)質(zhì)的巖性因子，有效地進行儲層的識別。因此，利用縱橫波疊前聯(lián)合反演獲取高精度的F和μρ，可以有助于解決致密儲層識別和流體檢測等問題。

1.1 FMD公式推導

疊前聯(lián)合反演的理論基礎是Zoeppritz公式，但是因為其表達形式復雜，許多學者對其進行了簡化近似。Russell等人[29]基于孔隙彈性介質(zhì)理論推導了基于Gassmann流體項f的縱波反射系數(shù)Rpp的線性近似方程，具體表達式如下：

上式中為縱波反射角和透射角的平均角度；γsat為上下層介質(zhì)飽和巖石縱橫波速度比的平均值；表示干巖石的縱橫波速度比的平方，該值的具體計算方法可以基于實驗室測量或者利用Gassmann方程的計算得到[10,30]，具體到實際應用來說，應該根據(jù)研究區(qū)條件靈活選擇估算方法，從而得到適用于研究區(qū)的數(shù)值；和Δμ分別為上下層介質(zhì)剪切模量的平均值和差；f還可以表示為β2M，其中，β為Boit系數(shù)(抽空條件下每個單位體積變化中孔隙體積的變化)，M為模量(在不改變孔隙體積的前提下，把流體壓入地層孔隙所需要的壓力)，f和Δf分別為上下層介質(zhì)Gassmann流體項的平均值和差。

Russell在Gassmann流體項基礎上又提出了Russell流體因子F，其表達式為

公式(2)中Ip，Is分別為縱波阻抗和橫波阻抗。

為了研究基于F、μρ和ρ的縱橫波聯(lián)合反演方法，首先根據(jù)孔隙介質(zhì)理論推導基于這三個參數(shù)的的縱波以及轉(zhuǎn)換波公式。針對縱波Rpp公式，我們首先從公式(1)的密度項提取公式(3)和公式(4)

則剩下的密度項為

公式(6)即為基于F、μρ和ρ的縱波近似系數(shù)表達式，在本文中稱作FMD公式。

Aki和Richards[31]在入射角度較小且反射界面兩側(cè)的彈性參數(shù)變化不大的假設基礎上推導了轉(zhuǎn)換波Rps反射系數(shù)表達式：

公式(7)中，為轉(zhuǎn)換波反射角和透射角的平均角度；、和分別為上下層介質(zhì)縱波速度、橫波速度和密度的平均值；ΔVp、ΔVs和Δρ分別為上下層介質(zhì)縱波速度、橫波速度和密度的差。已知橫波速度與剪切模量和密度的關系式為根據(jù)微分性質(zhì)可以得到則公式(7)可以改寫為公式(8)

首先，從公式(8)的密度項提取公式(9)

則公式(8)中剩下的密度項為

則公式(8)可以改寫為FMD方程的轉(zhuǎn)換波反射系數(shù)表達式：

1.2 FMD公式精度分析

為了驗證FMD公式的合理性，本文參考某致密砂巖儲層的實測井數(shù)據(jù)，設計了兩套模型對新推導的公式(6)、(11)進行精度分析。模型參數(shù)取值如表1所示，其中，模型1的兩套砂巖層填充的孔隙流體均為水，上覆巖層的孔隙度為10%，下伏巖層的孔隙度為5%；模型2中兩套巖層的孔隙度均為10%，但是上覆巖層的孔隙流體為氣，下伏巖層的孔隙流體填充的是水。分別利用精確的Zoeppritz公式、Aki＆Richards近似公式、Russell近似式以及本文新推導的FMD公式求取縱波反射系數(shù)，以及利用精確的Zoeppritz公式、Aki＆Richards近似公式和FMD公式求取轉(zhuǎn)換波反射系數(shù)，并進行對比。圖1和圖2分別是根據(jù)模型1和模型2數(shù)據(jù)得到的不同縱波反射系數(shù)、轉(zhuǎn)換波反射系數(shù)曲線對比圖。其中，黑色實線代表精確Zoeppritz公式，紅色實線代表Aki＆Richards公式，藍色實線代表Russell公式，綠色虛線代表新推導的FMD公式。在圖中可見，利用新公式求得的兩套模型的縱波和轉(zhuǎn)換波反射系數(shù)在35°內(nèi)時與精確Zoeppritz公式得到的反射系數(shù)值誤差非常小，并且在模型一中的求得的橫波反射系數(shù)比Aki＆Richards公式的精度要高一些。因此，F(xiàn)MD公式的精度在中、小角度入射情況下可以滿足聯(lián)合反演的要求，是直接提取相關彈性參數(shù)的理論基礎。

表1 致密砂巖氣儲層模型參數(shù)Table 1 Parameters of gas reservoirs in tight sandstone

1.3 疊前聯(lián)合反演方法

求解公式(6)和公式(11)時需要針對縱波數(shù)據(jù)以及轉(zhuǎn)換波數(shù)據(jù)各聯(lián)立M(M取大于等于3的整數(shù))個方程，其中每個入射角θi(i=1,2…M)對應一個方程，假設每道有M個入射角度，N個采樣點，設則公式(6)可以寫作

其中，

為了表達簡便，公式(12)可以寫為

其中，

同理，轉(zhuǎn)換波公式可以改寫為

圖1 模型1縱波反射系數(shù)(a)與轉(zhuǎn)換波反射系數(shù)(b)對比圖Fig. 1 Comparison of PP wave re flection coef ficients(a) and PS wave re flection coef ficients(b) of model 1

公式(14)中

其中，

則通過求解公式(13)、(14)即可得到待求參數(shù)，但是在實際運算中公式(13)和公式(14)有較強的不適定性，因此本文借助貝葉斯反演[32]，將模型參數(shù)的先驗信息引入反演的正則化項中，可以有效地改善反演問題的不適定性問題。根據(jù)貝葉斯理論，模型參數(shù)的后驗概率密度分布等于

其中，P(m)為模型參數(shù)的先驗分布；為似然函數(shù)，表征道集數(shù)據(jù)的噪聲分布；當后驗概率分布函數(shù)的形狀不變的情況下，P(dpp)P(dps)為常數(shù)，可以忽略。

在面對實際資料時，疊前道集數(shù)據(jù)含有一定程度的噪聲，本文假設噪聲是互不相關的，均值為0且服從高斯分布的似然函數(shù)，縱、橫波道集數(shù)據(jù)的噪聲似然函數(shù)可表示為

圖2 模型2縱波反射系數(shù)(a)與轉(zhuǎn)換波反射系數(shù)(b)對比圖Fig. 2 Comparison of PP wave re flection coef ficients(a) and PS wave re flection coef ficients(b) of model 2

其中，Cnp=，Cns=分別為縱波數(shù)據(jù)和橫波數(shù)據(jù)噪聲的協(xié)方差；，分別為縱波數(shù)據(jù)和橫波數(shù)據(jù)噪聲的均方差；I為單位矩陣。

模型參數(shù)的先驗分布可以分為單變量和多變量的，考慮到本文同時提取的三個參數(shù)的相關性，選取多變量分布可以降低由于三參數(shù)之間的相關性造成的反演病態(tài)問題。高斯分布僅可以產(chǎn)生一致性加權系數(shù)，這樣會影響反演結果的稀疏性，而柯西分布可在求解過程中可以得到非一致加權系數(shù)，產(chǎn)生稀疏效果，更具有地質(zhì)意義。因此，本文采用多變量柯西分布[33]作為模型參數(shù)的先驗分布，具體公式為

其中，ψ為相關矩陣，可以通過最大期望估計法得到[33]。Di為3N×3N的矩陣，其表達式為

把公式(16)、(17)、(18)帶入到公式(15)中，可以得到模型參數(shù)的后驗概率密度為

對公式(20)進行代數(shù)變換，便可以得到目標函數(shù)F(m)，其中目標函數(shù)為

由于低頻地震數(shù)據(jù)采集成本巨大，大多數(shù)的地區(qū)仍以常規(guī)采集為主，而常規(guī)采集中地震數(shù)據(jù)的低頻信息受噪聲影響嚴重，導致低頻信息無法從地震數(shù)據(jù)中獲得[34]，因此反演過程中需要合并可靠的低頻信息。低頻信息不僅提供了反演參數(shù)的整體趨勢，還增加了頻帶寬度，從而提高了反演結果的分辨率。因此，為了獲得低頻信息，還需要在目標函數(shù)中增加低頻約束。常用的低頻約束T如公式(22)所示

其中，PF，Pμρ，Pρ為積分矩陣，PFm，Pμρm，Pρm分別表示待求參數(shù)F、μρ、ρ的自然對數(shù)；LF，Lμρ，Lρ分別為真實的F、μρ、ρ自然對數(shù)的低頻成分，可以由初始模型進行濾波得到；λF，λμρ，λρ分別為F、μρ、ρ的正則項的權系數(shù)。以F為例，

通常，在實際操作中，低頻約束權系數(shù)的選擇通?；诠こ處煹慕?jīng)驗，通過不斷修改加權系數(shù)來嘗試獲得最佳的反演結果。通過公式(22)可見，常規(guī)的低頻約束是使反演得到的彈性參數(shù)的自然對數(shù)逼近模型的低頻成分，這種約束方法對加權系數(shù)非常敏感，加權系數(shù)的選擇不當則會影響反演結果的精度。因此本文在目標函數(shù)中引入了改進的低頻約束T′，表達式如下所示

其中，HF，Hμρ，Hρ分別為彈性參數(shù)的低通濾波矩陣?？梢钥吹礁倪M的低頻約束項使反演得到的彈性參數(shù)自然對數(shù)的低頻成分逼近模型的低頻成分[35]，這種方法對加權系數(shù)相對不敏感，不僅降低了人為因素的干擾，還在反演過程中降低了對中、高頻信息的影響，從而提高了反演精度。加入低頻約束項后，新得到的目標函數(shù)為：

最后，可采用迭代重加權最小二乘來求解目標函數(shù)，在迭代過程中可以不斷更新公式(24)中的得到F、μρ、ρ的反射系數(shù)，然后利用道積分即可得到F、μρ、ρ。

2 模型測試

為了驗證本文提出的方法的可行性和抗噪性，本文選取了一組實測致密砂巖井數(shù)據(jù)進行測試。首先，在反演前對實測數(shù)據(jù)進行了Backus平均[36]處理，將數(shù)據(jù)從測井尺度轉(zhuǎn)化為地震尺度，然后又進行了時深轉(zhuǎn)換，使數(shù)據(jù)從深度域轉(zhuǎn)換到了時間域。圖3中藍色實線表示處理后的實測井數(shù)據(jù)，黑色實線表示初始模型。然后，基于精確的Zoeppritz方程正演得到了實測井數(shù)據(jù)在不同采樣時間和不同角度(5°，15°，25°和 35°)下的縱波和轉(zhuǎn)換波反射系數(shù)，利用30 Hz的雷克子波與反射系數(shù)進行褶積，得到縱波以及轉(zhuǎn)換波的角道集合成記錄，如圖4所示。為了驗證本文提出的直接提取F、μρ、ρ方法的抗噪性，對合成記錄加入了信噪比為2的高斯隨機噪音，加入噪音后的角度道集合成記錄如圖5所示。

圖6和圖7分別為利用聯(lián)合反演和僅利用縱波數(shù)據(jù)反演得到的結果與實測數(shù)據(jù)的對比，圖中藍色實線為理論值，紅色實線為反演得到的值，黑色實線為初始模型，通過對比可見，雖然僅利用縱波數(shù)據(jù)也能得到合理的反演結果，但是三個參數(shù)的反演結果精度均低于聯(lián)合反演結果的精度，因為縱波數(shù)據(jù)對橫波和密度信息不敏感，而F、μρ、ρ參數(shù)均包含這兩種信息，因此聯(lián)合反演相對縱波反演可以有效提高F、μρ、ρ參數(shù)反演的精度。圖8和圖9為信噪比為2的情況下聯(lián)合反演和縱波反演的結果，從圖中可見在含有噪聲的情況下，聯(lián)合反演相對縱波反演同樣可以得到精度更高的F、μρ、ρ的反演結果，具有較好的抗噪性。但是在實際應用中，大部分的情況下只有縱波數(shù)據(jù)，轉(zhuǎn)換波數(shù)據(jù)比較少見，因此本文利用縱波反演方法進行了基于FMD公式直接求取彈性參數(shù)以及間接計算方法的對比。其中，間接計算是利用Aki & Richards近似公式反演得到縱波速度、橫波速度、密度，然后利用它們計算得到F、μρ。圖10為僅利用信噪比為2的縱波數(shù)據(jù)直接計算以及間接算得到的反演結果對比。為了更加清晰的對比，我們求取了圖10中兩種方法得到反演結果與理論值的相關系數(shù)(表2)，其中相關系數(shù)指越大說明反演結果越趨近理論值，反演結果精度越高。從表中可以看到直接反演結果精度要高于間接計算得到的結果，因為直接反演可以有效避免間接計算引入的累積誤差。

圖3 理論井數(shù)據(jù)及初始模型Fig. 3 Theoretical well-logging data and the initial model for inversion

圖4 合成角度道集；(a) PP波；(b) PS波Fig. 4 Synthetic angle gathers. (a) PP gathers; (b) PS gathers

圖7 基于PP波數(shù)據(jù)反演結果與實測井數(shù)據(jù)對比(不含噪音)Fig. 7 Comparison of inversion results with real logging data base on PP-wave inversion (noise free).

圖8 基于PP-PS波數(shù)據(jù)反演結果與實測井數(shù)據(jù)對比(S/N=2)Fig. 8 Comparison of inversion results with real logging data base on PP-PS joint inversion (S/N=2).

圖9 基于PP波數(shù)據(jù)反演結果與實測井數(shù)據(jù)對比(S/N=2)Fig. 9 Comparison of inversion results with real logging data base on PP-wave inversion (S/N=2).

圖10 直接反演結果與間接反演結果對比 (S/N=2)；(a) 直接反演；(b)間接反演Fig. 10 Comparison of direct inversion results and indirect inversion results (S/N=2). (a) direct inversion；(b)indirect inversion

表2 反演結果相關系數(shù)統(tǒng)計表Table 2 Comparison of the correlation coefficients of the inversion results

3 實際資料測試

為了進一步驗證FMD方程提取F和μρ參數(shù)的適用性，本文利用基于FMD方程的PP波反演和間接提取方法應用到了實際資料中。該研究區(qū)儲層呈低孔、低滲特征，為典型的致密儲層。首先，對保幅的疊前道集進行了預處理，并結合研究區(qū)實際情況，對疊前角度道集進行了分角度疊加，得到了 5°-15°，15°-25°，25°-35°的分角度疊加數(shù)據(jù)體(圖11)，圖11(a)的剖面中黑色橢圓標記處有含氣儲層發(fā)育，圖中黑色虛線所在位置為驗證井所在位置。然后，利用層位數(shù)據(jù)和井數(shù)據(jù)建立初始模型。最后，基于縱波數(shù)據(jù)，分別利用本文提出的基于FMD方程的直接反演方法和常規(guī)的間接反演方法對研究區(qū)進行F、μρ參數(shù)的提取。

圖12 反演結果對比剖面；(a) 新方法提取F；(b) 新方法提取μρ；(c) 間接計算F；(d) 間接計算μρFig. 12 Comparison of the inversion results.(a) inverted F by new method；(b) inverted μρ by new method；(c) inverted F by indirection method；(d) inverted μρ by indirection method

參數(shù)F能夠反映流體的情況，參數(shù)μρ對儲層巖性有一定的指示作用，因此在實際應用時，應將二者進行綜合對比解釋。通過巖石物理分析可知，研究區(qū)內(nèi)含氣儲層具有低F和高μρ的特征。圖12為反演結果對比圖，其中圖12(a)和圖12(b)為基于FMD方程得到的反演結果，圖12(c)和圖12(d)為間接算得的結果，通過對圖12(a)和圖12(c)的對比可見，在黑色橢圓處發(fā)育的含氣儲層中，間接計算得到的參數(shù)F剖面圖上含氣儲層響應不明顯，而利用FMD方程直接求得的參數(shù)F結果在含氣儲層處有較為強烈的響應；在圖中的黑色方框內(nèi)部，間接計算得到的F出現(xiàn)了流體響應的假象，而直接計算的結果并未存在這種假象；因此，若按照圖12(c)的結果進行油氣藏描述，會導致烴類流體預測的失??；通過對圖12(b)和圖12(d)的對比可見，利用FMD方程直接計算得到的參數(shù)μρ相對間接計算，橫向上更加連續(xù)，對于儲層的細節(jié)刻畫更加清晰(圖中紅色橢圓內(nèi)部)，有助于提高儲層預測精度。

為了更加明確地看到兩種反演結果的區(qū)別，本文在驗證井的井點處(圖11(a)黑色虛線)針對兩套反演結果分別提取了偽井曲線，并使其與Backus平均處理后的實測井曲線做對比，如圖13所示。通過偽井曲線的對比可見，利用FMD方程直接計算的方法提取的F、μρ參數(shù)與實際井曲線更加吻合。彈性參數(shù)的精度影響了儲層和流體預測的精度，因此利用本文提出的基于FMD方程直接提取F、μρ參數(shù)的方法有助于提高油氣藏描述的精度。

圖13 偽井反演結果對比；(a) F；(b) μρ；Fig. 13 Comparison of the inversion results of the pseudo well.(a) F；(b) μρ；

4 結論

本文新推導的基于F、μρ參數(shù)的縱波、轉(zhuǎn)換波反射系數(shù)方程精度滿足反演需求，為直接提取F、μρ參數(shù)提供了理論依據(jù)；將貝葉斯理論和雙項約束結合，直接從疊前數(shù)據(jù)中提取F、μρ參數(shù)是切實可行的，并且穩(wěn)定了反演過程；通過數(shù)字模型和實際資料的測試對比驗證了新方法的有效性，并且本文提出的反演方法相對直接求取方法有更高的抗噪性和精度，因此本方法具有較為廣闊的應用前景；在提取參數(shù)F時涉及干巖石縱、橫波速度比λdry的求取，而λdry參數(shù)的精度也影響了反演結果的精度，接下來還要進一步研究如何求的高精度的λdry值。