范曉輝,楊景陽(yáng),單 博

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

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

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

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

1 方法原理

1.1 AVO正演公式

ZOEPPRITZ[2]基于平面波理論,根據(jù)斯奈爾(Snell)定理、界面兩側(cè)位移和應(yīng)力的連續(xù)性,建立了用位移振幅比表示的反射透射系數(shù),對(duì)Zoeppritz方程組進(jìn)行簡(jiǎn)化推導(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è)矩陣中的速度、角度項(xiàng)均滿(mǎn)足Snell定律。則有:

(2)

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

(3)

簡(jiǎn)化矩陣形式為:

(4)

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

(5)

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

(6)

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

(7)

對(duì)(7)式泰勒展開(kāi)可得:

(8)

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

同樣展開(kāi)角度項(xiàng),則有:

(9)

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

(10)

其中,

A1=B1=T+1

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

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

(11)

其中,

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

(12)

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

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

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

(13)

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

(14)

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

(15)

其中,

(16)

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

(17)

對(duì)(17)式右端余弦項(xiàng)泰勒展開(kāi)

(18)

(19)

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

(20)

將(8)式代入(20)式

(21)

(22)

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

(23)

其中,

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

(24)

(23)式可表示為:

(25)

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

2 正演模擬

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

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

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

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

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

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

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

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

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

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

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

3 反演分析

3.1 一維層狀模型模擬

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

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

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

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

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

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

3.2 反演流程

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

d=Wr

(26)

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

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

(27)

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

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

(28)

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

(29)

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

(30)

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

3.3 應(yīng)用實(shí)例

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

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

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

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

4 結(jié)論

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

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

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

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

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

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