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

龔誠誠 吳國忱*② 單俊臻

(①中國石油大學(華東)地球科學與技術學院，山東青島 266580；②海洋國家實驗室海洋礦產資源評價與探測技術功能實驗室，山東青島 266580)

0 引言

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

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

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

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

1.1 散射矩陣分解

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

(1)

M=Mb+ΔM1+ΔM2

(2)

N=Nb+ΔN1+ΔN2

(3)

R=Rb+ΔR1+ΔR2

(4)

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

圖1 入射與散射的P-SV平面波完整系統

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

對式(1)變形得到如下矩陣表達式

MR=N

(5)

首先，僅考慮系數矩陣M、N的背景項Mb、Nb時，可以推得散射矩陣的背景項Rb

(6)

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

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

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

(7)

對式(7)進行整理，忽略一階以上高階項，得

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

(8)

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

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

=Nb+ΔN1+ΔN2

(9)

對式(9)進行整理，忽略二階以上高階項，得

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

(10)

整理二階擾動矩陣及其內部參數，經過計算可以得到ΔN2=ΔM2Rb，則式(10)可以進一步簡化為

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

(11)

(12)

分別計算散射矩陣的背景項、一階擾動項、二階擾動項并相加，整理可得入射、散射P-SV平面波完整系統下的縱波反射系數RPP關于縱、橫波速度與密度的二階近似定量表達式

(13)

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

1.3 橫波因子cosθS

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

由Snell定律可知

(14)

假設界面兩側介質參數變化較小，入射角、反射角與透射角均為實數，則界面兩側介質參數可由p的一階擾動項及其平均值的形式表示，即

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

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

(15)

進一步化簡并忽略高階擾動項，得

(16)

橫波因子可由實測縱、橫波速度與入射角定量表示為

(17)

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

(18)

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

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

2.1 二階項對精度的影響

若忽略式(18)中的二階擾動項，便得到關于縱、橫波速度以及密度的縱波反射系數一階近似方程

(19)

(20)

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

表1 模型參數

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

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

2.2 橫波因子對精度的影響

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

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

3 密度敏感性分析

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

改寫為關于縱波阻抗、橫波阻抗、密度的近似方程

(21)

圖4 Fatti近似公式RPP曲線

4 結束語

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