張建中 安 全* 于建明 陳 龍

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

0 引言

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

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

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

本文所述快速算法在探究NMO速度與地層各向異性參數(shù)及界面特征的關(guān)系、垂向非均勻TI介質(zhì)反射波旅行時(shí)反演等方面有實(shí)際應(yīng)用價(jià)值。

1 等效NMO速度的計(jì)算方法

1.1 NMO速度面解析解簡(jiǎn)述

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

(1)

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

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

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

圖1 觀測(cè)系統(tǒng)示意圖

(2)

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

(3)

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

(4)

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

顯然，零炮檢距射線相角γ與c、r關(guān)系為

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

(5)

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

而γ與γg滿足[17]

(6)

(7)

式中vg為群傳播速度，有

(8)

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

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

(9)

式中

(10)

根據(jù)文獻(xiàn)[22]，U滿足

(11)

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

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

1.2 單層Uint的計(jì)算

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

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

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

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

(12)

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

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

(13)

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

(14)

平面D內(nèi)的任一單位矢量b可表示為

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

(15)

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

(16)

式中

(17)

(18)

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

1.4 任一平面上等效W的計(jì)算

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

(19)

1.5 多層等效U的重構(gòu)

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

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

2 計(jì)算實(shí)例

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

表1 三層傾斜界面的TI介質(zhì)模型參數(shù)

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

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

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

(20)

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

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

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

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

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

3 結(jié)論

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

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

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