范曉輝,楊景陽,單 博

(中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院,山東青島 266580)

地球物理反演的本質(zhì)是依據(jù)地震數(shù)據(jù)獲取我們所需地層的物性、彈性參數(shù),揭示地下地層空間展布和含油氣性等特征。OSTRANDER[1]初次提出AVO技術(shù),揭示地震反射波振幅反映儲層特征的奧秘,由于Zoeppritz方程[2]精確求解的復(fù)雜性,許多學(xué)者開展了基于Zoeppritz方程的AVO方法研究。AKI等[3]在《定量地震學(xué)》中,對Zoeppritz方程進行了簡化近似。拉梅參數(shù)(λ,μ)是儲層預(yù)測和流體識別中最常用的彈性參數(shù),分別表征巖石的壓縮性和剛性,相較縱、橫波速度和密度參數(shù)具有更好的預(yù)測及識別能力[4]。GOODWAY等[5]將Fatti近似式改寫為模量組合和密度表示的形式,通過提取縱、橫波阻抗相對變化量求取縱、橫波阻抗信息,進而依據(jù)公式計算λρ和μρ,但其本質(zhì)仍是一種間接反演方法,而間接計算產(chǎn)生的累積誤差會降低預(yù)測精度[6-7]。XU等[8]提出利用體積模量K和拉梅常數(shù)λ,μ表征的完全隱含速度信息的近似方程。在密度相對變化較小時,能夠利用兩項參數(shù)反映入射角不超過45°的反射特性。GRAY等[9]基于Aki-Richards近似方程提出了由體積模量/壓縮模量、剪切模量和密度的相對變化率表示的反射系數(shù)近似公式,實現(xiàn)利用疊前地震數(shù)據(jù)直接反演出λ和μ,可以更加準確地識別巖性和流體。CONNOLLY[10]率先提出彈性阻抗的概念,能夠獲得縱、橫波速度和密度信息,進而計算其它彈性參數(shù)。王保麗等[11]基于Gray近似方程,推導(dǎo)了由λ、μ和ρ三參數(shù)表達式表示的彈性阻抗方程,并進行標準化處理,可以直接從彈性阻抗數(shù)據(jù)中提取拉梅參數(shù)。宗兆云等[12]基于平面波入射假設(shè)條件,推導(dǎo)了基于楊氏模量Y、泊松比σ和密度D的YPD反射系數(shù)近似方程,實現(xiàn)目標參數(shù)的直接反演,能夠有效識別頁巖氣“甜點”。

URSIN等[13]認為利用直接從常規(guī)縱波地震數(shù)據(jù)反演縱、橫波速度和密度3個參數(shù)難度較大。張春濤等[14]認為利用縱波地震數(shù)據(jù)一般很難反演得到密度比。疊前AVO反演能夠利用更多地下真實數(shù)據(jù),其本質(zhì)仍是不適定問題,一般通過刪除第三項提高反演問題的穩(wěn)定性[15]。SMITH等[16]利用Gardner經(jīng)驗關(guān)系式[17]改寫Aki-Richards近似方程,消除密度項,雖然穩(wěn)定性有所提高,但經(jīng)驗關(guān)系式與實際地層的吻合情況嚴重影響反演的精確度。

在Zoeppritz方程的基礎(chǔ)上,提出一種關(guān)于剪切模量和密度的AVO直接反演方法,該方法直接提取剪切模量和密度,避免間接計算引起的誤差。本方法僅涉及兩參數(shù)反演,且未引入經(jīng)驗關(guān)系式,具有比常規(guī)三參數(shù)反演更好的穩(wěn)定性。通過模型數(shù)據(jù)試算與方法對比,驗證了此方法正、反演問題的有效性與穩(wěn)定性。最后,將其應(yīng)用于實際疊前地震數(shù)據(jù),以驗證方法在巖性識別與儲層預(yù)測中的正確性。

1 方法原理

1.1 AVO正演公式

ZOEPPRITZ[2]基于平面波理論,根據(jù)斯奈爾(Snell)定理、界面兩側(cè)位移和應(yīng)力的連續(xù)性,建立了用位移振幅比表示的反射透射系數(shù),對Zoeppritz方程組進行簡化推導(dǎo):

(1)

式中:RPP、TPP、RPS、TPS分別為縱波反射系數(shù)、縱波透射系數(shù)、轉(zhuǎn)換橫波反射系數(shù)和轉(zhuǎn)換橫波透射系數(shù);vP1、vP2、vS1、vS2分別為界面兩側(cè)的縱、橫波速度;α1、β1、α2、β2分別為界面兩側(cè)的縱、橫波的反射和透射角;ρ1、ρ2分別為界面上、下地層的密度。(1)式左右兩側(cè)矩陣中的速度、角度項均滿足Snell定律。則有:

(2)

式中:P為射線參數(shù)。記x=sinα1,a=vP2/vP1,b=vS2/vS1,c=ρ2/ρ1,γ1=vS1/vP1,代入(1)式整理可得:

(3)

簡化矩陣形式為:

(4)

主要對縱波反射系數(shù)進行重新推導(dǎo),根據(jù)克萊姆法可知:

(5)

det表示求取矩陣的行列式,其中:

(6)

式中:A11是將A中的第一列用B替代所得到的矩陣。借鑒AVO公式常用各類參數(shù)相對變化率表示反射系數(shù)的方式,引入橫波、密度的相對變化率,記RS=(vS2-vS1)/(vS2+vS1),Rρ=(ρ2-ρ1)/(ρ2+ρ1),可對(3)式中參數(shù)進行轉(zhuǎn)換:

(7)

對(7)式泰勒展開可得:

(8)

式中:Rμ=(μ2-μ1)/(μ2+μ1),表示界面上、下剪切模量的相對變化率;T=γ1/γ2,表示界面上、下地層橫、縱波速度比的比值,可通過測井曲線獲得。

同樣展開角度項,則有:

(9)

將(8)式、(9)式代入(5)式,為保持精度,角度項保留至4次方項,整理可得:

(10)

其中,

A1=B1=T+1

(10)式為由剪切模量、密度的相對變化率表示的反射系數(shù)方程,記為μ-ρ方程,對其進行化簡得:

RPP(θ)=ARμ+BRρ+C

(11)

其中,

假定疊前地震數(shù)據(jù)有m個入射角,每一道地震數(shù)據(jù)有n個時間采樣點,結(jié)合褶積模型,反射系數(shù)與子波褶積得到合成地震記錄,即正演方程為:

(12)

式中:d(θ)為合成地震數(shù)據(jù);W(θ)為子波矩陣;RPP(θ)為PP波反射系數(shù);上標均代表對應(yīng)時間采樣點;下標均代表對應(yīng)入射角。

1.2 與Aki-Richards近似方程的比較

AKI等[3]在《定量地震學(xué)》中給出了Zoeppritz方程的經(jīng)典近似式為:

(13)

由于vP=vS/γ,結(jié)合縱波速度與剪切模量的關(guān)系,得:

(14)

將(14)式代入(13)式,Aki-Richards近似方程可改寫為:

(15)

其中,

(16)

對(16)式兩端取正弦得:

(17)

對(17)式右端余弦項泰勒展開

(18)

(19)

將(18)式和(19)式代入(17)式且忽略高階項,則有:

(20)

將(8)式代入(20)式

(21)

(22)

將(21)式和(22)式代入(15)式,保留相對變化率的一階近似,有:

(23)

其中,

考慮Aki-Richards近似方程成立的假設(shè)條件為相鄰兩層介質(zhì)的彈性參數(shù)變化較小,則

(24)

(23)式可表示為:

(25)

比較(10)式和(25)式可以發(fā)現(xiàn),假設(shè)相鄰兩層介質(zhì)的彈性參數(shù)變化較小,Aki-Richards近似方程在外推過程中,γ的平均值用界面上層的γ替代,從而得到Aki-Richards近似方程等價于μ-ρ方程,但μ-ρ方程在推導(dǎo)過程中沒有進行平均值的近似,避免了這部分近似的誤差。

2 正演模擬

為測試μ-ρ方程的可行性,基于CASTAGNA等[18]提出的AVO分類方法建立了4類AVO模型(表1至表4),且都利用Zoeppritz方程和μ-ρ方程(圖例以“μ-ρ”表示)、Aki-Richards近似方程和Gray近似方程計算4種模型界面的縱波反射系數(shù)并計算各類方法與精確值之間的誤差(圖1至圖4)。圖1a、圖2a、圖3a和圖4a分別是第Ⅰ、Ⅱ、Ⅲ、Ⅳ類AVO模型的正演結(jié)果;圖1b、圖2b、圖3b和圖4b分別是三類近似方法求得的第Ⅰ、Ⅱ、Ⅲ、Ⅳ類AVO模型的反射系數(shù)與Zoeppritz方程求得的準確值之間的誤差(已校正法向入射時初始誤差)。對比3種方法的正演結(jié)果可以看出,在0入射時,μ-ρ方程與Gray近似方程的正演結(jié)果均與準確值存在不同程度的誤差,說明在法向入射時兩種方法均存在直接引入彈性參數(shù)的轉(zhuǎn)化誤差。但就整體趨勢而言,μ-ρ方程正演結(jié)果在4類模型中均與實際模型趨勢吻合最好,誤差圖更加直觀地表明,誤差基本來自法向入射時計算誤差,在校正初始誤差后,本文提出的方法明顯比其它兩種方法更為精確。此外,Aki-Richards近似方程和Gray近似方程在4類模型中的正演結(jié)果具有相同的趨勢,且隨著入射角度的增加兩者逐步偏離精確Zoeppritz方程計算結(jié)果,其根本原因是Gray近似是基于Aki-Richards方程參數(shù)轉(zhuǎn)換所得,本質(zhì)上基于一種近似的不同表示,兩種方法涉及的角度項僅保留至二階,在大角度入射時易產(chǎn)生較大誤差,新方法的角度項保留至四階,能夠較好地避免此類誤差。

圖1 第Ⅰ類AVO模型不同方程計算的反射系數(shù)(a)和反射系數(shù)誤差(b)對比

圖2 第Ⅱ類AVO模型不同方程計算的反射系數(shù)(a)和反射系數(shù)誤差(b)對比

圖3 第Ⅲ類AVO模型不同方程計算的反射系數(shù)(a)和反射系數(shù)誤差(b)對比

圖4 第Ⅳ類AVO模型不同方程計算的反射系數(shù)(a)和反射系數(shù)誤差(b)對比

為驗證μ-ρ方程的穩(wěn)定性,對建立的4類AVO模型(表1至表4)進行條件數(shù)分析。圖5為4類AVO模型采用3種方法正演時,條件數(shù)隨最大入射角度的變化情況,其中最小入射角度保持為0,最大入射角度范圍為5°～40°。由圖5可以看出,當最大入射角度比較小時,3種方法的條件數(shù)都較大,表明過小角度入射時反演問題非常不穩(wěn)定,因此需要獲取較大的偏移量信息以提高反演的穩(wěn)定性。同時可以看出,μ-ρ方程所得條件數(shù)比以往3項參數(shù)的近似方程條件數(shù)成指數(shù)級降低,表明μ-ρ方程的穩(wěn)定性遠高于Aki-Richards近似方程和Gray近似方程。

表1 第Ⅰ類AVO模型參數(shù)

圖5 4類AVO模型條件數(shù)a 第Ⅰ類AVO模型; b 第Ⅱ類AVO模型; c 第Ⅲ類AVO模型; d 第Ⅳ類AVO模型

表2 第Ⅱ類AVO模型參數(shù)

表3 第Ⅲ類AVO模型參數(shù)

表4 第Ⅳ類AVO模型參數(shù)

3 反演分析

3.1 一維層狀模型模擬

為驗證μ-ρ方程反演的可行性,構(gòu)建一維五層模型進行反演模擬。反演算法基于最小二乘法,模型數(shù)據(jù)如表5所示。利用Zoeppritz方程獲得模型精確數(shù)據(jù),分別利用μ-ρ方程、Aki-Richards近似方程和Gray近似方程對精確數(shù)據(jù)進行μ、ρ反演并求取各方法計算結(jié)果與精確值之間的誤差,結(jié)果如圖6和圖7所示。對比3種方法的反演結(jié)果,兩種三參數(shù)近似方程在881～1040m層段與真實數(shù)據(jù)擬合相對較好。當深度較大時,反演誤差增大,認為是反演迭代中誤差累積所致。對比3種方法的結(jié)果認為,μ-ρ方程直接反演得到的μ、ρ與真實模型的擬合誤差最小,原因在于此方法直接提取彈性參數(shù),減少了累積誤差且兩參數(shù)反演具有更高的穩(wěn)定性。Aki-Richards近似方程和Gray近似方程在層狀模型反演時仍具有相同的趨勢,再次驗證兩種近似式本質(zhì)上是基于一種近似的不同表示。

圖6 一維層狀模型μ的反演結(jié)果(a)及誤差(b)

圖7 一維層狀模型ρ的反演結(jié)果(a)及誤差(b)

表5 一維層狀模型參數(shù)

假設(shè)井數(shù)據(jù)豐富,可以通過井位速度信息求取泊松比,結(jié)合μ-ρ方程直接反演得到μ,可以通過拉梅參數(shù)與泊松比的關(guān)系獲得λ。在本次模型模擬中代入每層真實泊松比,按上述方法計算λ,并利用Gray近似方程對精確數(shù)據(jù)進行直接反演得到λ,結(jié)果如圖8 所示。由圖8可以看出,此方法求解出的λ與真實數(shù)據(jù)吻合程度更高,表明在井數(shù)據(jù)豐富條件下,μ-ρ方程能準確獲得λ。

圖8 一維層狀模型λ的計算結(jié)果(a)及誤差(b)

3.2 反演流程

在實際地震反演中,通常采用離散方程進行計算,正演公式(12)可直接簡寫為矩陣形式:

d=Wr

(26)

式中:d為疊前地震數(shù)據(jù),W為子波矩陣,r為反射系數(shù)矩陣,即(11)式。

基于最小二乘原則約束反演誤差,考慮L0范數(shù)約束下的反演問題直接求解是NP-hard[19],采用L1范數(shù)構(gòu)建稀疏脈沖反演的目標函數(shù)為:

(27)

式中:‖·‖2表示測量反演誤差的L2范數(shù),‖·‖1表示反射系數(shù)的L1范數(shù),用于測量反射系數(shù)序列的稀疏性;l1是L1范數(shù)的正則化因子。

公式(27)中,L1范數(shù)測量了反射系數(shù)序列的垂向稀疏性,并沒有考慮多道反射系數(shù)序列的橫向變化和連續(xù)性,因此引入核模約束[20],以約束r為低秩,保證r相鄰道具有相關(guān)性,其定義為矩陣r所有奇異值的和,即

(28)

式中:σ(i)為矩陣r的第i個奇異值,則(27)式可以改寫為:

(29)

式中:l2為核模約束的正則化參數(shù)。在實際數(shù)據(jù)處理中,僅通過L1范數(shù)的稀疏約束和核模約束并不能完全準確地恢復(fù)低頻分量。因此,將模型參數(shù)的先驗信息rpriori作為正則化約束項加入到目標函數(shù)中,目標函數(shù)(29)可寫為

(30)

式中:l3為先驗?zāi)Ｐ偷恼齽t化約束系數(shù)。(30)式即為最終的疊前反演目標函數(shù)?？梢酝ㄟ^迭代重加權(quán)最小二乘算法[21]進行求解。

3.3 應(yīng)用實例

為驗證μ-ρ方程在實際資料反演中的可行性,選取某工區(qū)的實際地震數(shù)據(jù)進行測試。該工區(qū)的疊前地震數(shù)據(jù)包括3個不同角度范圍的部分角度疊加數(shù)據(jù)體,疊加剖面如圖9所示,其中,圖9a為0°～10°疊加剖面,圖9b為8°～17°疊加剖面,圖9c為15°～25°疊加剖面。首先,對工區(qū)內(nèi)的測井資料進行預(yù)處理。然后,利用求解基于μ-ρ方程的反演目標函數(shù)實現(xiàn)剪切模量和密度參數(shù)的反演,反演結(jié)果如圖10和圖11所示。圖10為剪切模量μ的反演剖面與井位反演曲線。圖10a中紅色曲線為井的位數(shù)據(jù)經(jīng)濾波后計算的μ曲線;左側(cè)為井的巖性圖。由圖10可以看出,反演剖面的高μ區(qū)與井的砂巖段吻合較好,說明通過反演μ可以較好地識別巖性。圖10b為井位位置的反演結(jié)果。其中紅色曲線為μ-ρ方程反演結(jié)果,證實基于μ-ρ方程的反演結(jié)果與測井曲線吻合較好。圖11為密度ρ的反演剖面與井位反演曲線,井位反演結(jié)果與測井曲線基本一致。圖10和圖11表明,基于μ-ρ方程得到的反演剖面與測井解釋結(jié)果基本吻合,能得到較為精確的μ和ρ值,將其與其它參數(shù)聯(lián)合分析、構(gòu)建流體因子,可更精確地識別巖性和流體。

圖9 部分角度疊加剖面a 0°～10°疊加剖面; b 8°～17°疊加剖面; c 15°～25°疊加剖面

圖10 剪切模量μ的三維反演剖面(a)及井位反演對比(b)

圖11 密度ρ的三維反演剖面(a)及井位反演對比(b)

4 結(jié)論

拉梅參數(shù)是儲層預(yù)測和流體識別中最常用的彈性參數(shù),目前常用的提取方法是通過反演縱、橫波速度間接計算獲得,這會導(dǎo)致誤差累積。同時,利用疊前地震數(shù)據(jù)實現(xiàn)三參數(shù)的精確反演仍是一個難題。我們提出了一種基于剪切模量和密度的兩參數(shù)縱波反射系數(shù)方程,并將該方程應(yīng)用于實際資料處理,能夠較好地直接反演出剪切模量和密度。正演模擬和反演分析結(jié)果表明:

1) 從Zoeppritz方程出發(fā),直接推導(dǎo)由剪切模量和密度相對變化率表征的反射系數(shù)方程,未引入經(jīng)驗關(guān)系式,且角度項保留至四階,在較大角度入射時同樣具有較好的正演效果。避免了反演出現(xiàn)累積誤差。

2)μ-ρ方程涉及兩參量,條件數(shù)遠低于Aki-Richards方程與Gray近似方程,具有更高的穩(wěn)定性。

3) 在井資料豐富的情況下,基于μ-ρ方程精確反演得到μ,進而結(jié)合井位信息以及拉梅參數(shù)與泊松比的關(guān)系式獲得λ,實現(xiàn)三參數(shù)的獲取。

4) 引入核模約束項以及先驗信息約束項構(gòu)建反演目標函數(shù),能夠增強反演的橫向連續(xù)性并還原原始數(shù)據(jù)的頻率信息,實際資料處理證實本方法能較好地擬合井數(shù)據(jù),反演得到較為精確的μ,ρ值,并用于巖性識別。

5) 正演模擬與反演分析結(jié)果均表明,μ-ρ方程具有較高的穩(wěn)定性和準確性,基于μ-ρ方程反演得到的彈性參數(shù)可與其它參數(shù)聯(lián)合分析、構(gòu)建流體因子,可更精確地識別巖性和流體。