傾斜層狀TI介質反射波旅行時快速計算

張建中 安 全* 于建明 陳 龍

(①內蒙古自治區地震局，內蒙古呼和浩特 010080； ②東方地球物理公司新興物探開發處，河北涿州 072751)

0 引言

地球內部廣泛存在各向異性介質[1-5]。TI介質是目前地震勘探理論研究和資料處理中經常使用的各向異性介質模型[6-9]。在反射波地震勘探中，非均勻各向異性介質反射旅行時一般由多方位、多炮檢距的射線追蹤計算。盡管已有的各向異性射線追蹤正演模擬足夠快，但當用于反射波反演時，因需反復計算多方位、多炮檢距旅行時，使反演過程非常費時[10-13]。對于最大排列長度與反射深度比接近1(簡稱常規排列)，介質具有不太復雜結構的CMP(Common Middle Point)反射時距關系可由測線方位上的NMO(Normal Moveout)速度描述[13-16]，因而可通過該反射時距關系計算反射波旅行時。這樣，NMO速度的計算則成為反射波旅行時快速計算的基礎。

在各向異性介質中，隨射線傳播角度而變化的速度使NMO速度方法的求解變得異常復雜。已有的計算各向異性NMO速度多數對應于簡單的TI模型，如：Thomsen[17]給出了水平界面VTI介質中同類波(非轉換波)的NMO速度解；Tsvankin[15，18]導出了水平界面HTI和TI介質對稱軸方位與界面走向一致時同類波的NMO速度解；Grechka等[19]推導了三維傾斜界面VTI介質同類波的NMO速度解；姚陳[20]導出了三維傾斜界面對稱軸任意取向的TI介質同類波NMO速度解。這些已有的NMO速度解都是針對單層介質，而由地面觀測得到的NMO速度一般是地下多層介質速度的等效結果。如果界面傾斜，無法再利用Dix公式[21]或廣義Dix方程[13]建立地表觀測的等效NMO速度和層速度的聯系。為此，Grechka等[22]提出了NMO速度面的概念，通過對沿零炮檢距射線路徑上不同層的速度面的Dix型平均化(Dix型平均算法)，得到水平地表上的等效NMO速度，該算法需要數值求解由介質彈性常數表示的Christoffel方程。

本文在現有的單層三維傾斜界面對稱軸任意取向TI介質NMO速度面解析解[23-24]的基礎上，運用Dix型平均算法[22]求解多層介質水平地表的等效NMO速度，進而運用由NMO速度描述的CMP反射時距關系實現反射波旅行時快速算法。與Grechka等[22]的數值計算NMO速度面不同，上述NMO速度面解析解完全由介質的各向異性參數和界面參數描述。

本文所述快速算法在探究NMO速度與地層各向異性參數及界面特征的關系、垂向非均勻TI介質反射波旅行時反演等方面有實際應用價值。

1 等效NMO速度的計算方法

1.1 NMO速度面解析解簡述

考慮多層傾斜界面對稱軸任意空間取向的TI介質，其常規排列同類波反射時距關系可由隨測線方位變化的NMO速度描述[13-16](假設每一炮檢對的反射波旅行時是單值)

(1)

式中：t是測線方位角為φ、炮檢距為X時同類反射波旅行時；t0為零炮檢距雙程反射旅行時；vnmo(φ)是測線方位角為φ時的NMO速度。

對于單層TI介質，利用姚陳[20]導出的NMO速度解，可由式(1)快速計算反射旅行時。對于多層傾斜TI介質，因來自地下不同層反射波旅行時t是其傳播路徑上多層介質的綜合反映，與t相對應的vnmo(φ)也是多層介質的綜合反映，因此需要先求得vnmo(φ)，再由式(1)快速計算來自不同界面的反射波旅行時。由于傾斜界面的存在，使等效NMO速度和層速度無法再通過Dix公式[21]或廣義Dix方程[13]建立聯系，因此本文運用由Grechka等[22]提出的NMO速度面概念，即把NMO速度做為三維空間方向單位矢量的函數，得到單層傾斜界面TI介質的NMO速度面，再運用Dix型平均算法計算水平地表上的等效NMO速度。

圖1給出的觀測系統示意圖中，設S為激發點，R為接收點，其中心點CMP即為坐標系原點O,測線SR在三維空間中的方向單位矢量為L，可將NMO速度vnmo(L)表示為[22]

圖1 觀測系統示意圖

(2)

式中：L=(sinθcosφ，sinθsinφ，cosθ)是單位矢量，其中θ和φ為L的傾角和方位角；U是一個3×3對稱矩陣

(3)

顯然，式(2)表示了vnmo(L)在三維空間中所有可能方向的矢徑，所有矢徑端點構成了NMO速度面；當θ=90°，式(2)表示的則是水平面上僅隨φ變化的vnmo(φ)。對于單層傾斜界面、對稱軸任意空間取向、任意各向異性強度的TI介質，文獻[23-24]將U的各元素顯式地表示為零炮檢距射線的相速度v及其對相角γ的導數dv/dγ、群角γg以及介質對稱軸方向c和反射界面法線方向r的解析函數。其中相速度v是對稱軸方向的qP(qS)波速度vP0(vS0)、各向異性參數(如ε、δ等)和相角γ的函數。c和r可表示為

(4)

式中：θc、φc分別為對稱軸的傾角和方位角；θr、φr分別為傾斜界面法線的傾角和方位角。

顯然，零炮檢距射線相角γ與c、r關系為

cosγ=c·r=sinθcsinθrcos(φc-φr)+cosθccosθr

(5)

S′、R′為S、R在水平地表的投影，P為界面上相反射點，G為界面上群反射點；γ為OP與c的夾角，φ為OS′與x1軸夾角，θ為OS與x3軸夾角；當炮檢距X=0時，零炮檢距射線OP與r方向一致，OG表示零炮檢距群反射路徑，與c的夾角為γg

而γ與γg滿足[17]

(6)

(7)

式中vg為群傳播速度，有

(8)

在得到零炮檢距射線相角γ后，將相速度v及式(5)～式(8)代入U的解析解，可求得群反射點G及零炮檢距雙程反射旅行時t0。

此外，U的元素滿足[24]

(9)

式中

(10)

根據文獻[22]，U滿足

(11)

式中：τ0是從零炮檢距射線反射點到CMP的單程旅行時；p=(p1，p2，p3)為在零炮檢距反射點激發、CMP位置的慢度矢量。對單一均勻層，Uk,m可由給定零炮檢距射線慢度矢量求解Christoffel方程獲得[22]。與此相比較，由給定的單層傾斜界面對稱軸任意空間取向TI介質參數，應用U解析解直接計算NMO速度面更為簡潔、快速。

特別指出的是，Grechka在理論上證實了均勻任意各向異性介質中同類反射波的NMO速度面為一橢圓柱面，其軸線平行于零炮檢距射線方向[22，24]。從幾何意義上理解，NMO速度面與三維空間內任一平面(與NMO速度面軸線平行的平面除外)的交線都為橢圓。由單層TI介質NMO速度面解析解[23-24]計算多層介質中來自反射界面同類波的水平地表上等效NMO速度，涉及到沿零炮檢距射線路徑上計算每一單層的層間Uint(l)(l為層序號，目標層為第1層向地面算起)，以及計算從第1層到第l層的等效U(l)。

1.2 單層Uint的計算

若l=1，零炮檢距射線方向矢量直接取反射界面法線空間取向(正入射原理)，從而計算單層的Uint(1)以及零炮檢距群反射點到CMP的單程旅行時τ0,1，此時Uint(1)亦是等效U(1)。若l>1，除目標層下界面(反射界面)外，零炮檢距射線與其路徑上的其他界面一般并不垂直，此時需要該層中零炮檢距射線方向參數代替式(5)中的界面法線參數θr和φr，然后再計算層l的Uint(l)。每層中的射線方向參數可由零炮檢距射線追蹤確定，同時獲得射線在每層中的單程旅行時τ0,l。

1.3 U與空間任一平面的交線W的計算[22]

NMO速度面與空間任一平面D的交線為該平面內的橢圓。平面D的法線單位矢量z由傾角θD和方位角φD表示為

z=(sinθDcosθD,sinθDsinφD,cosθD)

(12)

由該矢量，在平面D內可以找到兩個相互垂直單位矢量b(1)和b(2)

b(1)=(cosθDcosφD,cosθDsinφD,-sinθD)

(13)

b(2)=(-sinθD,cosθD,0)

(14)

平面D內的任一單位矢量b可表示為

b=b(1)cosφb+b(2)sinφb

(15)

式中φb為b相對b(1)的方位角。將b視作L，則φb=φ，將式(13)～式(15)代入式(2)，可得

(16)

式中

(17)

(18)

式(16)表示，平面D內隨方位角φ變化的NMO速度為橢圓，說明在已知等效U或Uint后，可以計算其與任一確定平面D的交線，即該平面上的等效W或層間Wint(元素Wi,j=Wj,i)。

1.4 任一平面上等效W的計算

在求得每層的Uint后，因Uint(1)等于U(1)，用式(17)可求界面l(l層上界面)與U(l)的交線W(l)，及界面l與Uint(l+1)的交線Wint(l+1)，然后由零炮檢距射線位移連續性及在界面上滿足Snell定律的Dix型平均算法[22]，獲得關于界面l+1的等效W(l+1)

(19)

1.5 多層等效U的重構

因NMO速度面軸線平行于零炮檢距射線[22-23]，因此，在得到W(l+1)后，可以通過式(9)求得相應的Uk,3(k=1,2,3)，從而實現對等效U(l+1)的重構。但在這個過程中，式(5)及式(10)中要用零炮檢距射線方向參數代替界面法線參數θr、φr。

重復1.3～1.5節的計算，可求得全部層位的等效U。最后，運用式(16)求得水平地面上的等效NMO速度橢圓vnmo(φ)，此時，式(12)中z=(0,0,1)。

2 計算實例

對于三層傾斜界面、對稱軸任意取向的TI介質模型(表1)，本文方法計算的的等效NMO速度面、等效NMO速度橢圓、CMP道集反射旅行時誤差如圖2所示，計算中P波相速度v使用了簡明的Thomsen弱各向異性近似[17]。由水平地表上第3界面P波反射等效NMO速度面(圖2a)可以看出，等效NMO速度面為一橢圓柱面(第1、第2界面P波反射等效NMO速度面亦為橢圓柱面)。圖2b為水平地表與三個界面P波反射等效NMO速度面的交線，即水平地表上獲取的三個界面反射的等效NMO速度橢圓，與常規排列下最優擬合射線追蹤計算旅行時得到的NMO速度(星號)在所有方位(間隔30°)都相符；第1、2及3界面的等效NMO速度橢圓長軸(粗直線)方位取向不同，分別為129°、118°和78°。圖2c計算了圖2b中3個橢圓上的等效NMO速度與最優擬合射線追蹤計算旅行時得到的NMO速度的相對誤差隨方位角的變化，可以看出：對應三個界面的等效NMO速度的最大相對誤差分別位于80°、25°、165°；隨著深度的增加，相對誤差也整體增大，但仍在10-3數量級內。在每一界面等效NMO速度相對誤差最大的方位上，計算了由本文方法及射線追蹤法得到的旅行時的相對誤差(圖2d)，可以看出：每層的反射旅行時相對誤差都很小(10-4～10-3數量級)，但隨炮檢距與界面深度比(X/dm)的增大，相對誤差以似指數形式增大；比較不同界面旅行時相對誤差發現，隨著界面深度增大，同一X/dm值對應的相對誤差也有所增大。

表1 三層傾斜界面的TI介質模型參數

圖2 等效NMO速度面、等效NMO速度橢圓及CMP道集旅行時精度分析(a)第3界面的等效NMO速度面；(b)水平面上三個界面的本文方法計算NMO速度橢圓與常規排列最優擬合射線追蹤計算的NMO速度(星號)對比；(c)NMO速度相對誤差隨方位角的變化曲線；(d)旅行時相對誤差隨方位角、炮檢距的變化曲線

可用本文方法計算共炮點(CSP)道集反射旅行時。設炮點處界面的深度為ds，它與測線方位角φ、炮檢距X及中心點深度dm關系為

dm(φ,X)=ds-0.5Xcosθrcos(φ-φr)

(20)

先根據式(20)計算出激發點與接收點的中心點距各層底界面的深度dm(φ,X)，再計算該中心點不同界面反射波的水平地表等效NMO速度和零炮檢距反射時間，然后按照式(1)計算出對應層的CSP道集反射旅行時。

圖3為本文方法計算的三層傾斜界面TI模型CSP集反射旅行時與射線追蹤計算旅行時之間相對誤差隨X/ds的變化曲線，可以看出：X/ds<0.5時，三個界面的反射旅行時相對誤差在所有方位上都非常小；隨X/ds增大，旅行時相對誤差也在增大；第1界面的相對誤差總體上小于第2層，所有方位上相對誤差在X/ds<1.2時亦處于10-4數量級；第1界面較大相對誤差出現在60°和240°方位，第2、第3界面較大相對誤差都出現在0°和180°方位；第3界面相對誤差明顯大于第1、第2界面，在X/ds<1.2時處于10-3數量級。由圖3可知，各層全方位反射旅行時相對誤差處于10-4～10-3數量級，但隨X/ds增大，以似指數形式增大；隨著反射深度增大，相同X/ds對應的各方位旅行時相對誤差也增大。

表1所示的三層傾斜界面TI模型具有一定的代表性。由圖2、圖3的計算結果可知，由單層傾斜TI介質NMO速度面解析解，運用Grechka給出的Dix型平均算法計算等效NMO速度，利用常規排列CMP時距關系實現了CMP及CSP集反射波旅行時的快速計算，且計算得到的旅行時具有較高的精度。但隨著炮檢距與界面深度比(X/ds或X/dm)的增大，各層界面反射波旅行時相對誤差以指數形式增大，且隨著反射深度增大，相同X/ds或X/dm對應的旅行時相對誤差也增大。這些現象可由與各向異性及垂向非均勻相關的非雙曲時距特性予以解釋，同時對本文提出的快速計算反射波旅行時方法給出約束條件。

對于炮檢距與反射深度比接近1(常規排列)，介質具有不太復雜結構的反射波CMP時距關系，可由雙曲方程式(1)較好地描述[13-16]，但當這一比值大于1后，受各向異性及垂向非均勻的影響，反射波時距曲線逐漸表現出明顯的非雙曲特性[13]，此時通過雙曲方程(式(1))計算的旅行時則是不可靠的，因而旅行時相對誤差必然增大。與此相類似的原因，在相同炮檢距與界面深度比的情況下，隨著反射深度增大，反射波時距也逐漸表現出了一定程度的非雙曲性，盡管從圖2d及圖3結果認識到，這種非雙曲性在炮檢距與界面深度比不大于1時還是很小的。有鑒于此，作為依據CMP時距的快速模擬計算反射波旅行時方法，在最大排列長度與反射深度比接近1，即在常規排列下，本文方法計算結果合理、可靠[14，24]。

圖3 CSP道集反射旅行時相對誤差隨X/ds的變化曲線(a)第1界面；(b)第2界面；(c)第3界面

3 結論

本文以單層傾斜TI介質NMO速度面解析解計算等效NMO速度為基礎，利用常規排列下CMP時距關系實現了傾斜層狀TI介質模型CMP及CSP集反射波旅行時的快速計算。數值實驗表明，等效NMO速度的計算結果與常規排列下最優擬合射線追蹤計算旅行時得到的NMO速度相比，其相對誤差為10-3數量級，而由等效NMO速度運用CMP時距計算的理論旅行時與射線追蹤計算旅行時相比，其相對誤差為10-4～10-3數量級。本文的快速計算反射波旅行時方法具有較高的計算精度，適用于介質結構不太復雜、最大排列長度與反射深度比接近1的常規排列。

本文的反射波旅行時的快速計算方法，因等效NMO速度橢圓(全方位上的等效NMO速度)的計算僅需通過正入射原理追蹤一條零炮檢距射線并由單層傾斜TI介質NMO速度面解析解運用Dix平均算法求得，因此，由反射時距關系快速計算反射波旅行時要比通常多方位角、多炮檢距的射線追蹤方法效率高得多，一般情況下可能快幾個數量級，適應用于TI介質旅行時反演[10-13]。

本文TI介質反射波旅行時快速計算方法可推廣到具有光滑彎曲界面的層狀各向異性模型。