蔣星達(dá)，張偉，陳曉非

1 南方海洋科學(xué)與工程廣東省實驗室(廣州)，廣州 511458 2 南方科技大學(xué)深圳市深遠(yuǎn)海油氣勘探技術(shù)重點實驗室，廣東深圳 518055 3 南方科技大學(xué)地球與空間科學(xué)系，廣東深圳 518055

0 引言

地下巖石受構(gòu)造應(yīng)力以及礦物定向排列等因素影響形成的周期性沉積、薄互層、定向裂縫等，可以表現(xiàn)出明顯的各向異性特征(Johnston and Christensen，1995；Helbig and Thomsen，2005).隨著油氣勘探對成像的精度要求提高，地下介質(zhì)的各向異性性質(zhì)越來越引起關(guān)注(Tsvankin et al.，2010；Thomsen，2014).以垂直方向為對稱軸的橫向各向同性層狀VTI介質(zhì)已經(jīng)被廣泛應(yīng)用于地震疊前時間偏移成像(Douma and De Hoop，2006；黃中玉等，2009；Jenner，2011)和非常規(guī)油氣微地震監(jiān)測(Li et al.，2013，2014；Perez et al.，2016；Michel and Tsvankin，2017；Gajek and Malinowski，2021；Kim et al.，2021)中.以上成像處理過程中均需要準(zhǔn)確計算VTI介質(zhì)中的地震波走時和射線路徑.各向異性介質(zhì)中地震波能量以群速度方向(群角)傳播，但是程函方程卻是以相速度角度(相角)為自變量，導(dǎo)致各向異性介質(zhì)中發(fā)展的各種射線追蹤算法需要搜索群角對應(yīng)的相角，影響計算效率(趙愛華等，2006；Sadri and Riahi，2010；Bai et al.，2013；Wang，2014；趙后越和張美根，2014；楊穎航等，2018；Li et al.，2020).探索適用于疊前時間偏移和微地震監(jiān)測的層狀VTI介質(zhì)高效射線追蹤技術(shù)對提高油氣勘探數(shù)據(jù)處理效率和微地震實時監(jiān)測具有重要意義.

基于初值問題的打靶法是經(jīng)典的射線追蹤算法.打靶法通過二分法修改射線出射角度，直到能夠準(zhǔn)確獲得震源到檢波點的射線路徑.打靶法原理簡單，易于實現(xiàn)，在各向同性介質(zhì)走時計算中獲得廣泛應(yīng)用(徐濤等，2004).打靶法的一個難點在于長偏移距時出射角度的微小變化可能導(dǎo)致射線路徑出現(xiàn)巨大擾動.田玥和陳曉非(2005)以射線在高速層位中的水平慢度代替射線出射角度有效的解決了這個問題，F(xiàn)ang和Chen(2019)進(jìn)一步將該方法擴展到具有速度梯度的一維速度模型，但目前該方法還未見到在強各向異性介質(zhì)中的應(yīng)用(Perez et al.，2016).最短路徑法基于圖論，利用網(wǎng)格節(jié)點表示速度模型，通過計算鄰近節(jié)點之間的走時信息選擇連接震源點到檢波點的最小走時節(jié)點連線作為射線路徑(Nakanishi and Yamaguchi，1986；Moser，1991).最短路徑法適合于復(fù)雜介質(zhì)，可以同時追蹤多個檢波點的初至波、反射波、折射波等各種震相走時，在各向同性和各向異性介質(zhì)中均獲得大量的應(yīng)用(Fischer and Lees，1993；劉洪等，1995；張美根等，2006a；唐小平和白超英，2009；Bai et al.，2009，2010；趙后越和張美根，2014；Li et al.，2020).最短路徑法的計算結(jié)果精度依賴于網(wǎng)格密度，節(jié)點連線形成的入射角和出射角不能連續(xù)變化，當(dāng)網(wǎng)格稀疏時，會產(chǎn)生不滿足Snell定律的“之”字型射線路徑.提高網(wǎng)格節(jié)點的密度可以提升最短路徑方法的計算精度，但會顯著增加計算時間(Dijkstra，1959；Moser，1991；劉洪等，1995；尤建軍等，2006；張美根等，2006a).張美根等(2006b)和黃國嬌等(2020)提出僅在速度分界面處設(shè)置網(wǎng)格節(jié)點的界面源法，減少了網(wǎng)格節(jié)點的個數(shù)，提高了計算效率.張建中等(2004)提出動態(tài)網(wǎng)格射線追蹤方法，網(wǎng)格節(jié)點的走時和射線路徑的計算不受固定節(jié)點位置的約束，在不增加節(jié)點個數(shù)的情況下顯著提升了射線追蹤算法的精度(李培明等，2013；桑運云等，2013，2014).在各向異性介質(zhì)中，最短路徑法計算節(jié)點之間的走時需要獲得兩點連線方向的群速度，計算給定群角的群速度顯著增加了計算時間(白海軍等，2011)，而前人發(fā)展出的各種群速度與群角近似公式(Byun et al.，1989；Sena，1991；趙愛華和丁志峰，2005)只適合于弱各向異性介質(zhì).

本文綜合打靶法和最短路徑法的優(yōu)點提出了一種層狀VTI介質(zhì)高效射線追蹤方法.該方法不需要尋找特定群角對應(yīng)的群速度，而是以相角出發(fā)，從震源點出射一系列射線，基于打靶法計算射線同各層速度邊界的交點位置和走時.將射線交點作為網(wǎng)格節(jié)點，基于節(jié)點信息計算檢波點走時信息和射線路徑.當(dāng)檢波點周圍的網(wǎng)格密度不足時，通過自適應(yīng)加密增加檢波點周圍的網(wǎng)格節(jié)點個數(shù)，最終利用插值和動態(tài)網(wǎng)格方法(張建中等，2004)獲得檢波點的走時和射線路徑.與傳統(tǒng)VTI介質(zhì)最短路徑方法相比，本方法網(wǎng)格節(jié)點設(shè)置靈活，射線路徑不受節(jié)點位置限制，在內(nèi)存和CPU占用資源少的基礎(chǔ)上可以快速同時獲得多個檢波點的射線路徑和走時信息.該方法適用于任意強度層狀VTI介質(zhì)，計算效率高，在需要計算直達(dá)波、反射波的微地震監(jiān)測、垂直地震剖面(Vertical Seismic Profiling, VSP)成像、井間地震成像、疊前偏移等油氣勘探領(lǐng)域具有廣泛的應(yīng)用價值.

1 方法描述

本文提出了一種根據(jù)稀疏網(wǎng)格節(jié)點計算水平層狀VTI介質(zhì)中震源到檢波點射線路徑和走時的高效射線追蹤方法.該方法流程如圖1所示，主要包括網(wǎng)格模型建立和檢波點射線追蹤計算兩個步驟.在建立網(wǎng)格模型過程中，我們根據(jù)從震源出射的稀疏相角精確計算VTI介質(zhì)中的群速度和群角，根據(jù)群角方向計算出射射線與各層位的交點作為稀疏網(wǎng)格節(jié)點，并獲得各個節(jié)點的精確走時.若檢波點周圍網(wǎng)格節(jié)點走時差距較大，則在檢波點周圍局部空間內(nèi)加密網(wǎng)格節(jié)點，精確計算各節(jié)點的走時信息.在檢波點射線追蹤計算過程中，通過三次樣條插值算法獲得檢波點的走時信息，并利用動態(tài)網(wǎng)格法反傳快速計算檢波點到震源的射線路徑信息.

1.1 VTI介質(zhì)中的群速度和群角

地震波在各向異性介質(zhì)中傳播時會分裂為qP、qSV、qSH三種波形.VTI介質(zhì)中相速度與相角具有如下關(guān)系(Thomsen, 1986)：

圖1 本文高效射線追蹤方法流程圖Fig.1 The flowchart of our high-efficient raytracing method

(1)

其中[vP,vSV,vSH]表示qP、qSV、qSH波相速度值(圖2)，[α0,β0,ε,δ,γ]是表征VTI介質(zhì)性質(zhì)的Thomsen參數(shù).α0為qP波垂向傳播相速度；β0為qS波垂向傳播相速度；ε為qP波水平和垂直方向速度變化率大小，用來度量qP波各向異性強度；γ為qSH波水平和垂直方向速度變化率大小，用來度量qSH波各向異性強度;δ為qP波垂向各向異性變化的快慢程度(Tsvankin, 1997).

圖2 各向異性介質(zhì)波前示意圖(Thomsen，1986)其中φ為群角，θ為相角.Fig.2 Schematic diagram of the wavefront in anisotropic media (From Thomsen, 1986)φ is group angle, θ is phase angle.

相角和相速度是地震波沿波前面垂直方向變化所展現(xiàn)的表觀性質(zhì)，實際上的地震波能量是以群速度大小沿著群角方向傳播.VTI介質(zhì)中，群速度/群角和相速度/相角具有如下關(guān)系：

(2)

其中[V,φ]為群速度與群角.將相速度公式(1)代入(2)式，即可獲得qP、qSV、qSH波給定相角時的群速度和群角表達(dá)式(具體表達(dá)式見附錄A).

由公式(2)可知，強各向異性VTI介質(zhì)中群速度和群角之間沒有直接表達(dá)式，而是通過相角連接.給定相角時可以通過公式(A1)和(A2)計算出群速度和群角，但是如果計算特定群角方向的群速度，則需要先搜索確定相應(yīng)的相角，才能通過公式(A1)計算出群速度.盡管可以在走時計算前通過掃描相角建立群角和群速度查找表，提高傳統(tǒng)最短路徑方法節(jié)點之間走時計算的效率，但建立高分辨率的查找表仍然對計算效率具有明顯影響.

1.2 網(wǎng)格模型建立

與傳統(tǒng)最短路徑法設(shè)置固定網(wǎng)格節(jié)點不同，本文以從震源按不同相角出射的一系列射線與界面的交點作為網(wǎng)格節(jié)點(圖3).

圖3 直達(dá)波網(wǎng)格節(jié)點分布示意圖其中S代表震源，k代表層位序號.Fig.3 Schematic diagram of grid nodes in direct waveS is source, k is layer number.

1.2.1 確定射線出射相角分布

如果在震源處采用等相角分布出射射線，有可能在遠(yuǎn)偏移距處界面節(jié)點間距過大.為此，本文采用了一種簡易的不均勻出射相角分布：以震源與相鄰層位上等間距分布的點連線構(gòu)成的角度作為出射射線相角(圖4)，相應(yīng)的相角角度為：

圖4 射線出射相角計算示意圖Fig.4 Schematic diagram of the calculation of ray emergence phase angle

(3)

h為震源與鄰近層位的垂直距離，Rn為層位上的某一個點距離震源的水平偏移距.Rn=nΔR，ΔR為震源相鄰層位節(jié)點間隔，n為出射相角個數(shù).需要注意的是，以上公式只是為了給出一個非均勻的出射相角分布，真實的射線路徑按照相角對應(yīng)的群角出射.

1.2.2 射線在層狀VTI介質(zhì)中傳播

射線在不同層位傳播過程中相角遵循Snell定律，

(4)

其中下標(biāo)k代表層位序號(圖3)，[θk+1,θk]為同一層位界面的入射和出射相角，[vk+1(θk+1),vk(θk)]為對應(yīng)的入射和出射相速度，p為射線參數(shù).當(dāng)入射相角θk+1已知時，根據(jù)公式(1)和(4)可以確定相速度vk+1(θk+1)和射線參數(shù)p，出射相角θk可以推導(dǎo)獲得(具體推導(dǎo)過程請見附錄B).

對于qP波和qSV波,

(5)

對于qSH波，

(6)

以上公式中的介質(zhì)參數(shù)為k層的值，這里為了表達(dá)簡潔沒有顯示下標(biāo).在qP波和qSV波中選擇符合0°≤θk≤90°的正實數(shù)作為最終的出射相角，進(jìn)而可以計算該層位射線傳播的群角和群速度(附錄A).當(dāng)遇到qSV波多值射線問題時，我們選擇最小的走時(即按照最大速度傳播)作為射線到達(dá)檢波點時間.

1.2.3 網(wǎng)格節(jié)點位置和走時

確定出射相角后可以計算對應(yīng)的群角以及群速度，進(jìn)而可以計算該層射線與相鄰層位的交點，確定相鄰層位節(jié)點的位置和走時(圖5).該位置和走時可以表示為：

圖5 射線傳播路徑和網(wǎng)格節(jié)點位置計算Fig.5 The computation of grid node position and ray propagation path

(7)

1.2.4 局部網(wǎng)格節(jié)點加密

在獲得各層網(wǎng)格節(jié)點位置后，由于檢波點與網(wǎng)格節(jié)點位置基本不會正好重合，因此需要根據(jù)網(wǎng)格節(jié)點走時和射線路徑計算檢波點的走時和射線路徑信息.當(dāng)檢波點兩端網(wǎng)格節(jié)點走時差距較大時，我們選擇在兩個節(jié)點對應(yīng)的兩條射線之間添加更多的射線，加密檢波點周圍節(jié)點的密度.具體步驟表示為：

(1)獲取檢波點兩側(cè)最近節(jié)點之間的走時差Δt，判斷Δt是否大于規(guī)定時間差ΔT；

(3)重新獲取檢波點兩側(cè)最近節(jié)點走時差Δt′，若Δt′>ΔT，則重復(fù)步驟(2)；若Δt′≤ΔT，則停止.

通過以上操作，可以保證檢波點最終的走時誤差在小于ΔT的精度范圍.

1.3 檢波點走時和射線路徑計算

檢波點周圍的網(wǎng)格節(jié)點如圖6所示.根據(jù)檢波點所在層位的網(wǎng)格節(jié)點時間，檢波點D的走時tr,D可以利用三次樣條函數(shù)插值獲得：

tr,D=aD+bD(xr,D-xr,A2)+cD(xr,D-xr,A2)2

+dD(xr,D-xr,A2)3，(8)

系數(shù)[aD,bD,cD,dD]由A2、B2界面所有節(jié)點位置和走時獲得.

圖6 檢波點D的射線路徑示意圖CD為需要計算的射線路徑；A1、A2、B1、B2為檢波點兩端網(wǎng)格節(jié)點；r為層位.Fig.6 Schematic diagram of the raypath of receiver DCD is the raypath to be calculated; A1, A2, B1, B2 are grid nodes on both sides of the receiver; r is the layer number.

傳統(tǒng)最短路徑法通常以檢波點同網(wǎng)格節(jié)點連接具有最小走時的路徑作為射線路徑，但是存在“之”字形路徑問題.張建中等(2004)提出的動態(tài)網(wǎng)絡(luò)射線追蹤方法不受固定節(jié)點限制，基于Fermat原理確定出不位于網(wǎng)格節(jié)點上的射線路徑，本文應(yīng)用該方法實現(xiàn)VTI介質(zhì)的更準(zhǔn)確的射線路徑計算，其原理可表示為：

根據(jù)圖6，已知節(jié)點A1、B1的二維坐標(biāo)和走時分別為[xr+1,A1,zr+1,A1,tr+1,A1]，[xr+1,B1,zr+1,B1,tr+1,B1].在已知檢波點(或射線路徑與層位的交點)D的坐標(biāo)和走時[xr,D,zr,D,tr,D]和其射線傳播路徑節(jié)點C的縱坐標(biāo)zr+1,C條件下，C點的走時tr+1,C可以通過A1、B1界面節(jié)點走時三次樣條插值獲得：

tr+1,C=aC+bC(xr+1,C-xr+1,A1)+cC(xr+1,C

-xr+1,A1)2+dC(xr+1,C-xr+1,A1)3，(9)

系數(shù)[aC,bC,cC,dC]由A1、B1界面所有節(jié)點位置和走時獲得.

tr+1,C和tr,D關(guān)系可以表示為：

(10)

sCD為射線CD方向的群慢度，可根據(jù)A1、B1界面所有已知節(jié)點的群慢度插值獲得：

sCD=aCD+bCD(xr+1,C-xr+1,A1)+cCD(xr+1,C-xr+1,A1)2

+dCD(xr+1,C-xr+1,A1)3，(11)

系數(shù)[aCD,bCD,cCD,dCD]由A1、B1界面所有節(jié)點位置和群慢度獲得.

(1)確定檢波點所在的層位以及兩側(cè)網(wǎng)格節(jié)點，根據(jù)公式(9)獲得檢波點的走時信息，計算檢波點與相鄰層位的射線路徑交點；

(2)以該層位射線路徑交點為新的“檢波點”，重復(fù)以上步驟(1)獲得下一個層位的射線交點.根據(jù)新的射線交點繼續(xù)計算，直到到達(dá)震源所在層位；

(3)將以上各層獲得的射線交點與震源和檢波點依次連接，獲得該檢波點在VTI模型中的射線路徑信息.

反射波網(wǎng)格節(jié)點的確定與直達(dá)波相似(圖7)，除卻在反射層位依據(jù)反射定律射線沿反射角傳播外，各層節(jié)點位置以及走時信息均可以通過公式(7)獲得.進(jìn)而根據(jù)插值方法計算獲得各個檢波點的走時信息和射線路徑.

圖7 反射波網(wǎng)格節(jié)點示意圖Fig.7 Schematic diagram of grid nodes in reflected wave

2 數(shù)值實驗

為了驗證本文算法的有效性以及實用性，我們設(shè)計了一個含高速夾層的強各向異性多層VTI介質(zhì)模型(圖8).位于地面和井中的震源以及檢波器用來模擬疊前時間偏移和微地震井中監(jiān)測的數(shù)據(jù)采集方式.兩個震源的位置分別為(500 m, 1200 m)和(500 m, 0 m).井中檢波器共20個，水平位置為1000 m，垂直方向沿290 m至670 m按20 m等間隔排列.地面共有40個檢波器，按照50 m等間隔均勻水平排列.

圖8 VTI介質(zhì)模型(a) 震源點(500 m, 1200 m)與井中檢波器； (b) 震源點(500 m, 0 m)與地面檢波器； (c) 各層Thomsen參數(shù)大小.Fig.8 VTI medium model(a) Source (500 m, 1200 m) and the receivers in the downhole well; (b) Source (500 m, 0 m) and the receivers on the surface; (c) The Thomsen parameters in each layer.

2.1 單個檢波器走時和射線路徑計算

首先以震源(500 m, 1200 m)到地面單個檢波器(1500 m, 0 m)的走時和射線路徑為例展示本文方法的計算過程和結(jié)果(圖9).對震源上下相鄰的1100 m和1500 m層位，以50 m為間距計算出射射線，根據(jù)震源到每個射線點的連線計算出射相角，在震源兩側(cè)共計算得到74條射線(上下層位各37條).計算這些射線沿著各層傳播的相角，根據(jù)相角對應(yīng)的群角和群速度獲得各層速度界面的網(wǎng)格節(jié)點位置和走時.由圖9a所示，射線密度在遠(yuǎn)偏移距并未出現(xiàn)明顯擴散，動態(tài)網(wǎng)格節(jié)點分布較密集.檢波器兩側(cè)節(jié)點位置為(1471.6 m, 0 m)(1526.9 m, 0 m)，間隔超過50 m.節(jié)點走時為565.57 ms和572.90 ms，走時差為7.33 ms.檢波器并未與網(wǎng)格節(jié)點重合，走時無法直接獲得(圖9b).為了保持檢波器兩側(cè)最近節(jié)點之間的走時差小于1.00 ms，計算加密射線條數(shù)為8條(7.33/1.00≈8)，在檢波器兩側(cè)最近節(jié)點之間從震源出發(fā)添加8條射線,獲得添加射線的走時和路徑信息(圖9c和圖9d).此時，檢波器兩側(cè)節(jié)點位置分別為(1497.9 m, 0 m)(1503.8 m, 0 m)，兩側(cè)節(jié)點走時為569.03 ms和569.81 ms.因為兩側(cè)節(jié)點走時差為0.78 ms，小于預(yù)設(shè)精度1.00 ms，停止插入新的射線.通過三次樣條插值方法獲得檢波器的走時為569.30 ms，與真實值569.20 ms相對誤差小于0.1%.通過動態(tài)網(wǎng)格法逐層回傳計算獲得檢波器到震源的射線路徑，該射線路徑與檢波器真實射線路徑基本重合.射線路徑和走時的準(zhǔn)確計算驗證了本文方法的可行性.

圖9 檢波器(1500 m, 0 m)射線追蹤計算(a) 研究介質(zhì)中的網(wǎng)格節(jié)點； (b) 地面網(wǎng)格節(jié)點走時； (c) 檢波器走時計算； (d) 檢波器射線路徑計算. Node:動態(tài)網(wǎng)格節(jié)點. New:重新插入的射線. True:真實射線路徑和走時. Receiver:計算的檢波器射線路徑和走時.Fig.9 The raytracing calculation of receiver (1500 m, 0 m)(a) The grid nodes in the study medium; (b) The traveltime of the grid nodes on the surface; (c) The calculation of the traveltime of the receiver; (d) The calculation of the raypath of the receiver. Node:The nodes in dynamic networks. New: The new inserted rays. True: The true raypath and traveltime. Receiver: The calculated raypath and traveltime of the receiver.

圖11 不同方法反射波射線路徑和走時對比(a) qP波真實射線路徑與最短路徑法結(jié)果對比； (b) qP波真實射線路徑與本文高效射線追蹤法結(jié)果對比； (c) 走時對比. True:真實結(jié)果. Shortest:傳統(tǒng)最短路徑法結(jié)果. High-eff:本文高效射線追蹤方法結(jié)果.Fig.11 The comparison of raypath and traveltimes with different methods in reflected wave(a) The comparison of raypath between the true results and the shortest-path method in qP wave; (b) The comparison of raypath between the true results and our high-efficient method in qP wave; (c) The comparison between different traveltimes. True: True result. Shortest: The result of the traditional shortest-path method. High-eff: The result of our high-efficient method.

圖12 最短路徑法和本文射線追蹤方法對比(a) 不同網(wǎng)格間距兩種方法的計算時間和計算誤差對比； (b) 不同偏移距兩種方法的計算時間和計算誤差對比. Time:Shortest:最短路徑法消耗時間. Time:High-eff: 本文高效射線追蹤方法消耗時間. Misfit:Shortest: 最短路徑法計算誤差. Misfit: High-eff: 本文高效射線追蹤方法計算誤差.Fig.12 The comparison between the shortest-path method and our high-effective raytracing method(a) The comparison of calculation time and misfit in different grid spaces between the two methods; (b) The comparison of calculation time and misfit in different offsets between the two methods. Time:Shortest: The calculation time of the shortest-path method. Time:High-eff: The calculation time of our high-effective raytracing method. Misfit:Shortest: The calculation misfit of the shortest-path method. Misfit:High-eff: The calculation misfit of our high-effective raytracing method.

2.2 精度對比

我們通過同VTI介質(zhì)最短路徑射線追蹤算法結(jié)果對比驗證本文算法精度.傳統(tǒng)最短路徑法沿層位界面設(shè)置固定節(jié)點，節(jié)點間隔為1.0 m，依據(jù)Dijkstra算法依次計算單個節(jié)點與相鄰層位所有節(jié)點的走時，選擇其中的最小走時節(jié)點作為該節(jié)點的地震波射線傳播路徑.由于每次計算節(jié)點之間走時需要獲得群角方向?qū)?yīng)的群速度，本文按照1.0°等間隔計算了群角在[0°，90°]區(qū)間與群速度對應(yīng)表，兩節(jié)點之間走時計算需要通過查表獲得相應(yīng)的群速度.

2.2.1 井中監(jiān)測

本測試模擬了微地震井中監(jiān)測的數(shù)據(jù)采集方式.震源(500 m, 1200 m)與垂直檢波器的射線路徑和走時如圖10所示(射線路徑以qP波為例，qSV和qSH波類似).相比于最短路徑法，本文射線追蹤算法獲得的射線路徑與準(zhǔn)確射線路徑重疊程度更高，在各個層位的交點更加接近.qP、qSV、qSH三種波的走時對比顯示本文射線追蹤算法走時曲線與真實曲線更加接近，計算精度更高.

2.2.2 地面監(jiān)測

本測試模擬了疊前時間偏移的數(shù)據(jù)采集方式.地面檢波器接收地面震源(500 m, 0 m)在介質(zhì)層位深度1000.0 m的反射波射線路徑和走時如圖11所示(以qP波為例).同直達(dá)波計算結(jié)果相同，本文方法在反射層位的反射點更加接近準(zhǔn)確反射點的位置，射線路徑與準(zhǔn)確路徑更加相似.本文方法相比于最短路徑法qP、qSV、qSH波的走時精度更高，走時誤差均小于0.1%，顯示本文射線追蹤方法在反射波計算方面同樣具有優(yōu)越性.

2.3 效率對比

我們對比了最短路徑法和本文射線追蹤方法的計算效率.兩種方法消耗時間和對應(yīng)的計算精度如圖12所示.圖12a顯示了震源(500 m, 1200 m)到檢波點(1500 m, 0 m)在不同網(wǎng)格間距下最短路徑法和本文射線追蹤方法計算直達(dá)波消耗時間和計算誤差對比.從圖中可以明顯看出，隨著節(jié)點間距的增大，最短路徑法效率逐漸提高，但是計算誤差也逐漸增大.而本文射線追蹤方法隨著初射間隔的增大(圖4中Rn)，消耗時間變化不明顯，計算誤差也一直保持小于1%.在節(jié)點間距為1 m時，最短路徑法與本文方法計算誤差相當(dāng)，但是計算消耗時間是本文方法的幾十倍.節(jié)點間距的增大提高了最短路徑法的計算效率，但是相比本文射線追蹤方法計算誤差明顯增加，在節(jié)點間距為10 m時誤差已經(jīng)超過了10%.圖12a顯示了震源(500 m, 0 m)以1000.0 m為反射層位在不同偏移距下最短路徑法節(jié)點間距為1 m時消耗時間和計算精度對比.在較小偏移距時(0～500 m)，最短路徑法網(wǎng)格范圍較小，計算消耗時間較少，同時誤差較低.隨著偏移距的增大，最短路徑法計算網(wǎng)格范圍較大，計算量增加，消耗時間增加，由于網(wǎng)格離散的特點計算走時誤差也逐漸積累.而本文射線追蹤方法隨著偏移距增大計算耗時未出現(xiàn)明顯變化，計算誤差也一直保持在0.1%左右的較高精度，證明本文方法受震源到檢波點的距離影響較小.

由上述兩個例子可以看出，最短路徑法按照1 m的間隔設(shè)置網(wǎng)格節(jié)點，提高了走時計算精度，但是在研究區(qū)域范圍較大時顯著增加了計算量，計算網(wǎng)格節(jié)點之間的走時信息需要查找群角和群速度表，降低了計算效率.本文射線追蹤方法得益于稀疏網(wǎng)格和局部節(jié)點加密的特點，計算精度受初射間距影響較少，同時不用計算群角和群速度對應(yīng)關(guān)系，在較大網(wǎng)格間距和較遠(yuǎn)的偏移距時仍可以通過插值獲得高精度的走時信息.由此可看出本文射線追蹤方法相比最短路徑法計算時間耗費少，效率高，非常適于層狀VTI介質(zhì)快速射線追蹤正演計算.

3 結(jié)論

本文提出一種適用于強各向異性的水平層狀VTI介質(zhì)高效射線追蹤方法.本方法利用從震源按照不同相角出射的一系列射線與速度界面的交點作為稀疏網(wǎng)格節(jié)點，精確計算網(wǎng)格節(jié)點的位置和走時信息.在檢波點周圍自適應(yīng)增加網(wǎng)格節(jié)點密度.利用三次樣條插值法獲得檢波點的高精度走時信息，同時基于動態(tài)網(wǎng)格思想計算獲得射線路徑信息.本方法吸收了打靶法快速高效的特點，同時避免了長偏移距射線路徑巨大擾動的難題，在借鑒最短路徑法利用網(wǎng)格節(jié)點進(jìn)行多個檢波點射線追蹤的基礎(chǔ)上，避免了高密度網(wǎng)格節(jié)點引起的計算效率降低問題.本方法不受網(wǎng)格節(jié)點位置限制，可以獲得不位于網(wǎng)格節(jié)點的射線路徑，克服了節(jié)點稀疏情況下的射線路徑“之”字形缺陷.基于嚴(yán)格計算群角與群速度關(guān)系，本方法適用于任意強度VTI介質(zhì)qP、qSV、qSH波射線追蹤.本方法不用多次遍歷群角和群速度關(guān)系，在提高計算效率的基礎(chǔ)上可以高精度計算直達(dá)波、反射波等射線路徑信息.同傳統(tǒng)最短路徑方法的對比驗證了本文方法走時和射線路徑計算的高精度和高效性.本方法占用硬件資源少，耗時短，非常適合于需要多次快速正演計算的地震疊前時間偏移、井中微地震實時監(jiān)測等情況.

附錄A 群速度和群角計算公式

將qP、qSV和qSH波的相速度公式(1)代入群速度和群角公式(2)，可以得到任意強度VTI介質(zhì)群速度/群角和相角的關(guān)系：

(A1)

(A2)

(A3)

[VP,VSV,VSH]和[φP,φSV,φSH]為相應(yīng)的qP、qSV、qSH波群速度和群角.

附錄B 根據(jù)Snell定律計算各向異性介質(zhì)中的出射相角

當(dāng)入射相角θk+1已知，qP波射線參數(shù)p可以準(zhǔn)確獲得(公式(4)).根據(jù)Snell定律，射線在不同層位傳播射線參數(shù)相同，

(B1)

根據(jù)以上公式，可得出以下方程：

(B2)

(B2)式可以看作是sin2θk的一元二次方程，根的形式為：

(B3)

因此，θk的取值為

(B4)

對于qSV波，上述過程同樣成立，可獲得

(B5)

最終可以獲得與(B2)相同的方程，因此qSV波的出射角度θk與(B4)式相同.

對于qSH波，根據(jù)Snell定律可得

(B6)

可以獲得sin2θk的表達(dá)式為

因此θk的取值為

(B8)

通過以上公式，可以精確計算qP、qSV、qSH波在不同層位中的相角出射角度.