馬琦琦，孫贊東

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

0 引言

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

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

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

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

1 方法原理

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

1.1 FMD公式推導

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

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

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

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

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

則剩下的密度項為

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

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

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

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

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

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

1.2 FMD公式精度分析

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

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

1.3 疊前聯合反演方法

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

其中，

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

其中，

同理，轉換波公式可以改寫為

圖1 模型1縱波反射系數(a)與轉換波反射系數(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)即可得到待求參數，但是在實際運算中公式(13)和公式(14)有較強的不適定性，因此本文借助貝葉斯反演[32]，將模型參數的先驗信息引入反演的正則化項中，可以有效地改善反演問題的不適定性問題。根據貝葉斯理論，模型參數的后驗概率密度分布等于

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

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

圖2 模型2縱波反射系數(a)與轉換波反射系數(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=分別為縱波數據和橫波數據噪聲的協方差；，分別為縱波數據和橫波數據噪聲的均方差；I為單位矩陣。

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

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

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

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

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

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

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

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

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

2 模型測試

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

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

圖3 理論井數據及初始模型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波數據反演結果與實測井數據對比(不含噪音)Fig. 7 Comparison of inversion results with real logging data base on PP-wave inversion (noise free).

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

圖9 基于PP波數據反演結果與實測井數據對比(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 反演結果相關系數統計表Table 2 Comparison of the correlation coefficients of the inversion results

3 實際資料測試

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

圖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

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

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

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

4 結論

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