基于散射矩陣分解的反射系數(shù)二階近似

龔誠(chéng)誠(chéng) 吳國(guó)忱*② 單俊臻

(①中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，山東青島 266580；②海洋國(guó)家實(shí)驗(yàn)室海洋礦產(chǎn)資源評(píng)價(jià)與探測(cè)技術(shù)功能實(shí)驗(yàn)室，山東青島 266580)

0 引言

AVO技術(shù)利用振幅信息研究巖性、檢測(cè)油氣，經(jīng)過數(shù)十年的發(fā)展，已經(jīng)形成了一套完整而嚴(yán)謹(jǐn)?shù)睦碚擉w系。Zoeppritz[1]研究了平面波入射在水平界面上的縱橫波反射、透射情況，建立了Zoeppritz方程。如何利用地震反射振幅隨入射角變化的信息檢測(cè)油氣成為地球物理學(xué)家關(guān)注的問題，但是由于Zoeppritz方程的結(jié)構(gòu)過于復(fù)雜且物理含義并不直觀，不易進(jìn)行數(shù)值計(jì)算，一直沒有得到很好的應(yīng)用。Muskat等[2]首先提出了AVO概念，即“平面波的反射和透射系數(shù)是入射角的函數(shù)”；Koefoed[3]將泊松比與反射系數(shù)相聯(lián)系，給出了精確Zoeppritz方程的近似公式，使AVO技術(shù)得到快速發(fā)展；Bortfeld[4]在界面兩側(cè)巖性參數(shù)變化小的假設(shè)下得到了區(qū)分流體、固體的近似公式；Richards等[5]在假設(shè)介質(zhì)參數(shù)變化小的情況下根據(jù)位移和應(yīng)力連續(xù)條件推導(dǎo)出較直觀的反射/透射系數(shù)近似表達(dá)式；Aki等[6]給出了反射/透射系數(shù)精確表達(dá)式，并利用射線參數(shù)與角度的關(guān)系給出了弱參數(shù)變化下的反射/透射系數(shù)近似式。Ostrander[7]利用AVO技術(shù)識(shí)別“亮點(diǎn)”型含油氣砂巖，使AVO 技術(shù)用于實(shí)際生產(chǎn)。在前人的基礎(chǔ)上，人們從不同角度推導(dǎo)出各類近似表達(dá)式，表征縱波反射系數(shù)與不同巖石物理參數(shù)間的關(guān)系。如Shuey[8]給出了具有不同角度項(xiàng)的關(guān)于泊松比的縱波反射系數(shù)近似式；Smith 等[9]利用Gardner公式，采用加權(quán)分析提出了速度隨密度變化的經(jīng)驗(yàn)公式；Hilterman[10]在 Shuey近似的基礎(chǔ)上，給出了關(guān)于泊松比的另一種近似式；鄭曉東[11]、楊紹國(guó)等[12]推出了Zoeppritz方程的近似冪級(jí)數(shù)表達(dá)式；Fatti等[13]發(fā)布了關(guān)于相對(duì)波阻抗的縱波反射系數(shù)近似式；Goodway等[14]提出了關(guān)于拉梅參數(shù)的縱波反射系數(shù)近似式；Xu等[15]引入剪切模量、體積模量和壓縮模量簡(jiǎn)化了縱波反射系數(shù)；Gray等[16]給出了關(guān)于剪切模量、體積模量和壓縮模量相對(duì)變化量的近似表達(dá)式；Russell等[17]在多孔流體飽和巖石的前提下提出了關(guān)于Russell流體因子的縱波反射系數(shù)近似表達(dá)式；宗兆云等[18]推出了關(guān)于楊氏模量和泊松比縱波反射系數(shù)近似方程。

常規(guī)AVO反演利用精確Zoeppritz方程的一階近似式，更加適用于弱介質(zhì)變化的反射界面、小角度或小炮檢距反射問題，其假設(shè)條件導(dǎo)致計(jì)算誤差[19]，不利于準(zhǔn)確提取密度參數(shù)[20]，不能很好地預(yù)測(cè)復(fù)雜儲(chǔ)層或中深部?jī)?chǔ)層，且無法充分利用近臨界角數(shù)據(jù)。因此，針對(duì)復(fù)雜儲(chǔ)層或中深部?jī)?chǔ)層，需要改進(jìn)反射系數(shù)近似方程，使其適應(yīng)近臨界角入射，充分利用大角度地震數(shù)據(jù)提高反演精度。Xu等[21]基于Aki-Richards近似方法給出了具有歸一化速度、密度二階項(xiàng)的近似公式；Wang[22]對(duì)縱波反射系數(shù)精確方程的射線參數(shù)表達(dá)形式進(jìn)行泰勒展開，得到了關(guān)于垂直慢度、射線參數(shù)的二階近似公式；Ramos[23]對(duì)縱波反射系數(shù)精確方程進(jìn)行三階泰勒展開，給出了適用于強(qiáng)反射界面的橫波反射系數(shù)近似公式；Ursin等[24]利用上、下行波特征向量矩陣及矩陣對(duì)稱性對(duì)反射/透射系數(shù)進(jìn)行泰勒二階展開，得到二階近似公式；Charles[25]提出最優(yōu)Zoeppritz近似，通過研究縱波反射系數(shù)精確值與縱橫波速度和密度相對(duì)變化量各階泰勒展開的關(guān)系得到擬線性近似，提高了近臨界角處的反射系數(shù)反演精度。

基于反射系數(shù)精確表達(dá)式利用Aki-Richards近似方法進(jìn)行高階近似過于繁瑣[22-23]，利用特征向量化簡(jiǎn)物理意義不夠明確[24]。為此，本文直接從P-SV平面波入射/散射矩陣出發(fā)，給出了一種求取散射矩陣高階近似的方法，即利用擾動(dòng)思想將散射矩陣分解為背景矩陣與一階、二階擾動(dòng)矩陣，求取縱波反射系數(shù)背景項(xiàng)與一階、二階擾動(dòng)項(xiàng)，推得縱波反射系數(shù)的二階近似公式。模型對(duì)比分析表明，所推公式在中高角度乃至近臨界角入射情況下具有較高的精度，對(duì)密度參數(shù)的敏感性更高，為充分利用大炮檢距地震數(shù)據(jù)準(zhǔn)確地反演物性參數(shù)提供了基礎(chǔ)。

1 基于散射矩陣分解的高階近似

1.1 散射矩陣分解

Aki等[6]研究了P-SV平面波在固體—固體分界面上的反射與透射問題，得到了入射與散射P-SV平面波完整系統(tǒng)(圖1)，并基于此系統(tǒng)推導(dǎo)了完整的散射矩陣。

(1)

M=Mb+ΔM1+ΔM2

(2)

N=Nb+ΔN1+ΔN2

(3)

R=Rb+ΔR1+ΔR2

(4)

式中：Mb、Nb、Rb分別為M、N、R的背景矩陣，由各向同性介質(zhì)參數(shù)構(gòu)成；ΔM1、ΔN1、ΔR1分別為M、N、R的一階擾動(dòng)矩陣，由介質(zhì)參數(shù)的一階變化構(gòu)成；ΔM2、ΔN2、ΔR2分別為M、N、R的二階擾動(dòng)矩陣，由介質(zhì)參數(shù)的二階變化構(gòu)成。

圖1 入射與散射的P-SV平面波完整系統(tǒng)

1.2 縱波反射系數(shù)二階近似定量表達(dá)(附錄A)

對(duì)式(1)變形得到如下矩陣表達(dá)式

MR=N

(5)

首先，僅考慮系數(shù)矩陣M、N的背景項(xiàng)Mb、Nb時(shí)，可以推得散射矩陣的背景項(xiàng)Rb

(6)

當(dāng)某種波型入射時(shí)，由于背景介質(zhì)均為各向同性，故背景介質(zhì)中不存在該波型的反射以及其他波型的反射與透射，在Rb中相同波型的透射系數(shù)為1，其他波型的反射、透射系數(shù)為0。

假設(shè)M、N與R僅由背景矩陣與一階擾動(dòng)矩陣構(gòu)成，則式(5)可表示為

(Mb+ΔM1)(Rb+ΔR1)=(Nb+ΔN1)

(7)

對(duì)式(7)進(jìn)行整理，忽略一階以上高階項(xiàng)，得

ΔR1=(Mb)-1(ΔN1-ΔM1Rb)

(8)

假設(shè)M、N與R由背景矩陣與一階擾動(dòng)矩陣、二階擾動(dòng)矩陣共同構(gòu)成，則式(5)可表示為

(Mb+ΔM1+ΔM2)(Rb+ΔR1+ΔR2)

=Nb+ΔN1+ΔN2

(9)

對(duì)式(9)進(jìn)行整理，忽略二階以上高階項(xiàng)，得

ΔR2=(Mb)-1(ΔN2-ΔM1ΔR1-ΔM2Rb)

(10)

整理二階擾動(dòng)矩陣及其內(nèi)部參數(shù)，經(jīng)過計(jì)算可以得到ΔN2=ΔM2Rb，則式(10)可以進(jìn)一步簡(jiǎn)化為

ΔR2=(Mb)-1(-ΔM1ΔR1)

(11)

(12)

分別計(jì)算散射矩陣的背景項(xiàng)、一階擾動(dòng)項(xiàng)、二階擾動(dòng)項(xiàng)并相加，整理可得入射、散射P-SV平面波完整系統(tǒng)下的縱波反射系數(shù)RPP關(guān)于縱、橫波速度與密度的二階近似定量表達(dá)式

(13)

式中： Δρ、ΔVP、ΔVS分別為界面兩側(cè)介質(zhì)密度、縱波速度、橫波速度的差值；ρ、VP、VS分別為界面兩側(cè)介質(zhì)密度、縱波速度、橫波速度的平均值；θP、θS分別為界面兩側(cè)下、上行縱(θP1、θP2)、橫波入射角(θS1、θS2)的平均值；p為射線參數(shù)。

1.3 橫波因子cosθS

觀察式(13)發(fā)現(xiàn)，基于散射矩陣得到的縱波反射系數(shù)的二階項(xiàng)受橫波入射角、透射角的平均值θS影響，但在實(shí)際縱波反演方法中不存在θS。因此，文中將cosθS稱為橫波因子，并由反射界面的彈性參數(shù)近似表示。

由Snell定律可知

(14)

假設(shè)界面兩側(cè)介質(zhì)參數(shù)變化較小，入射角、反射角與透射角均為實(shí)數(shù)，則界面兩側(cè)介質(zhì)參數(shù)可由p的一階擾動(dòng)項(xiàng)及其平均值的形式表示，即

VS1=VS-ΔVS/2

VS2=VS+ΔVS/2

VP1=VP-ΔVP/2

VP2=VP+ΔVP/2

θS1=θS-ΔθS/2

θS2=θS+ΔθS/2

θP1=θP-ΔθP/2

θP2=θP+ΔθP/2

將各參數(shù)代入式(14)，得

(15)

進(jìn)一步化簡(jiǎn)并忽略高階擾動(dòng)項(xiàng)，得

(16)

橫波因子可由實(shí)測(cè)縱、橫波速度與入射角定量表示為

(17)

將式(17)代入式(13)，可得到縱波反射系數(shù)的另一種表達(dá)形式

(18)

式(18)為化簡(jiǎn)后的關(guān)于縱、橫波速度以及密度的縱波反射系數(shù)二階近似式。與式(13)相比，用橫縱波速度比的形式表示射線參數(shù)p與橫波因子cosθS，使近似公式的參數(shù)形式得到統(tǒng)一，便于選取實(shí)際反演參數(shù)，更加適用于AVO反演。

2 縱波反射系數(shù)二階近似方程精度分析

2.1 二階項(xiàng)對(duì)精度的影響

若忽略式(18)中的二階擾動(dòng)項(xiàng)，便得到關(guān)于縱、橫波速度以及密度的縱波反射系數(shù)一階近似方程

(19)

(20)

本文借鑒Wang[22]的四種模型，分別利用不同近似表達(dá)式給出了不同巖石界面下的縱波反射系數(shù)曲線。這四種模型分別表示頁(yè)巖/砂巖、頁(yè)巖/石灰?guī)r(或白云巖)、硬石膏/砂巖和硬石膏/石灰?guī)r(或白云巖)界面，模型參數(shù)如表1所示。

表1 模型參數(shù)

圖2為RPP曲線(未考慮橫波因子)。由圖可見：①當(dāng)為頁(yè)巖/砂巖(圖2a)、頁(yè)巖/石灰?guī)r(或白云巖)(圖2b)界面時(shí)，上層介質(zhì)的波阻抗小于下層介質(zhì)(正阻抗界面)，Aki-Richards近似方程僅在小角度入射時(shí)得到的RPP具有較高精度，在中高角度處會(huì)產(chǎn)生明顯誤差，而二階近似方程在中高角度仍能取得很好的近似效果，且本文得到的二階近似方程(式(18))在近臨界角處仍能取得很高的近似效果；②由于頁(yè)巖/砂巖(圖2a)界面兩側(cè)密度與縱波速度變化小、橫波速度變化大，頁(yè)巖/石灰?guī)r(或白云巖)(圖2b)界面兩側(cè)密度與速度變化趨勢(shì)相反，因此后者的二階近似方程精度更高，表明二階RPP反映了橫波的影響，即變化小的橫波速度提高了二階RPP的精度；③當(dāng)為硬石膏/砂巖(圖2c)、硬石膏/石灰?guī)r(或白云巖)(圖2d)界面時(shí)，上層介質(zhì)波阻抗大于下層介質(zhì)(負(fù)阻抗界面)，在小角度入射時(shí)三種近似方程得到的RPP均具很高的精度，但在中高角度處僅有二階近似方程可以滿足精度要求，且本文推導(dǎo)的二階近似方程(式(18))所產(chǎn)生的誤差最??；④由于硬石膏中的密度、縱波速度、橫波速度遠(yuǎn)大于砂巖與石灰?guī)r，介質(zhì)參數(shù)變化大，進(jìn)行射線參數(shù)與橫波因子處理時(shí)，在大角度情況下會(huì)產(chǎn)生較大誤差，二階近似方程的精度隨著入射角的增大而逐漸變低。

圖2 RPP曲線(未考慮橫波因子)

2.2 橫波因子對(duì)精度的影響

將橫波因子表達(dá)式(式(17))代入式(13)與式(20)，可得到RPP曲線(考慮橫波因子)(圖3)。對(duì)比未考慮橫波因子(圖2)與考慮橫波因子(圖3)的RPP曲線可知，后者的精度較高，確保了近似方程在入射角增大情況下的適用性，如本文提出的二階近似方程(式(18))曲線與精確方程曲線在近臨界角處幾乎完全擬合。當(dāng)界面兩側(cè)物性差異較大時(shí)(圖3c、圖3d)，考慮橫波因子明顯提高了近似方程的精度。

圖3 RPP曲線(考慮橫波因子)

3 密度敏感性分析

基于AVO反演可以從地震數(shù)據(jù)中獲取相應(yīng)介質(zhì)參數(shù)，小炮檢距(小入射角)地震振幅數(shù)據(jù)主要與地震波速度相關(guān)，密度參數(shù)對(duì)中高炮檢距(中高入射角)地震振幅數(shù)據(jù)更敏感[26-29]，但是由于縱波反射系數(shù)一階近似方程在中高入射角時(shí)無法獲得很高的近似效果，因此基于此類方程的AVO反演方法的密度反演精度相對(duì)有限[30]。Fatti等[13]提出了相對(duì)波阻抗的概念，并在Aki-Richards的反射系數(shù)近似式的基礎(chǔ)上推導(dǎo)了關(guān)于縱波阻抗、橫波阻抗、密度的Fatti近似式。由于縱、橫波阻抗及縱、橫波速度都與密度直接相關(guān)，其中密度對(duì)阻抗的影響最大[31]。因此為了更清晰地體現(xiàn)本文所提近似方程(式(18))在中高角度地震數(shù)據(jù)密度反演中的優(yōu)勢(shì)，將式(18)

改寫為關(guān)于縱波阻抗、橫波阻抗、密度的近似方程

(21)

圖4 Fatti近似公式RPP曲線

4 結(jié)束語

本文基于P-SV平面波入射與散射完整系統(tǒng)下的散射矩陣，給出了一種利用矩陣分解推導(dǎo)反射/透射系數(shù)高階近似的方法，將散射矩陣分解為背景矩陣與不同階數(shù)擾動(dòng)矩陣，得到了基于入射角的縱波反射系數(shù)二階近似公式。與常規(guī)AVO分析中的Aki-Richards線性近似[6]、Wang[22]的二階近似表達(dá)式相比，文中所提反射系數(shù)二階近似公式在中高角度入射時(shí)與精確方程具有良好的吻合性，且在近臨界角入射情況下也具有較高的精度，當(dāng)界面兩側(cè)橫波速度與密度變化差異適中時(shí)，近似效果可以達(dá)到最佳。因此利用本文的近似公式可以充分利用大炮檢距地震振幅數(shù)據(jù)，解決了一階近似公式未能充分利用臨界角數(shù)據(jù)的問題，提高了中高入射角的反射系數(shù)反演精度，并且為近臨界角反演提供了一種方法。利用AVO反演可以提取各種彈性、物性參數(shù)，其中密度參數(shù)對(duì)小入射角地震響應(yīng)的影響弱，在大炮檢距地震振幅數(shù)據(jù)中富含更精細(xì)的信息，常規(guī)地震方法很難獲得穩(wěn)定的密度參數(shù)估計(jì)。本文的近似公式增強(qiáng)了對(duì)密度的敏感性，提高了富含密度信息的大炮檢距地震振幅數(shù)據(jù)的利用程度，為準(zhǔn)確而穩(wěn)定地密度反演提供了基礎(chǔ)。