蘆永明，張偉1,,3*

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

0 引言

研究表明在地殼、上地幔和地核內(nèi)部廣泛地存在地震各向異性(Shearer and Orcutt,1985;Tanimoto and Anderson,1985;Song and Richards,1996).它對(duì)地震波傳播的振幅和走時(shí)有很大的影響(Tsvankin and Thomsen,1995;Xu and Mao，2018)，如果忽略了各向異性的影響會(huì)造成基于走時(shí)的成像和反演的誤差.因此，準(zhǔn)確計(jì)算各向異性介質(zhì)中走時(shí)是必要的.在各向同性介質(zhì)中，體波只有壓縮波(P波)和剪切波(S波)，并且相速度(走時(shí)梯度方向)和群速度(射線方向)方向一致.但是在各向異性介質(zhì)中，體波有三種類型的波：qP波、qSV波和qSH波.每種波以各自的傳播速度和極化方向傳播，并且相速度和群速度方向不再一致(黃光南等，2015)，簡(jiǎn)單應(yīng)用各向同性走時(shí)計(jì)算方法有可能引起違背時(shí)間因果性的問(wèn)題，使得求解初至走時(shí)變得比較困難.

對(duì)于各向異性介質(zhì)，快速行進(jìn)法只能解決一些簡(jiǎn)單的情況(Bin Waheed and Alkhalifah，2017)，并且這類方法在各向異性介質(zhì)中，如果從波前上最小走時(shí)點(diǎn)擴(kuò)展波前面，可能會(huì)違背時(shí)間的因果條件.而快速掃描法不需要存儲(chǔ)和追蹤波前，可以克服該問(wèn)題，在求解向異性介質(zhì)中的走時(shí)方面得到了廣泛的應(yīng)用.Kao等(2004)用 Lax-Friedrichs 差分格式求解VTI介質(zhì)中的哈密爾頓雅可比方程.對(duì)該方法，人工黏性參數(shù)的選取會(huì)極大影響方法的精度和速度.為了提高該方法的精度和效率，Zhang等(2006)用高階加權(quán)的非震蕩方案改進(jìn)Lax-Friedrichs格式.Bin Waheed等(2015)、Bin Waheed和Alkhalifah(2017)提出將qP波的走時(shí)求解分解成橢圓項(xiàng)和非橢圓項(xiàng).先迭代求解橢圓項(xiàng)，然后將非橢圓項(xiàng)看成高階的誤差迭代地?cái)M合.該方法的優(yōu)勢(shì)是可以擴(kuò)展到方位各向異性介質(zhì)中，但是求解需要更多的迭代次數(shù).Le Bouteiller等(2018)用快速掃描法求解二維TI介質(zhì)中qP波的走時(shí)，他們用間斷伽遼金方法求解時(shí)間依賴的哈密爾頓雅可比方程.為了使快速掃描法適用于強(qiáng)各向異性介質(zhì)，Han等(2017)發(fā)展了基于慢度四次方程的快速掃描法計(jì)算二維TTI介質(zhì)中qP波的走時(shí).他們通過(guò)構(gòu)造局部解把qP波和qSV波的解耦慢度四次方程轉(zhuǎn)換成走時(shí)的四次方程數(shù)值求解.對(duì)走時(shí)四次方程求導(dǎo)得到三次方程，求出三次方程的根作為四次方程的求解子區(qū)間，在每個(gè)子區(qū)間內(nèi)使用二分法尋找可能的走時(shí)解.Huang等(2020)、Huang和Luo(2020)也發(fā)展了快速掃描法計(jì)算二維TI介質(zhì)中qP、qSV和qSH波的走時(shí).他們首先離散化各向異性介質(zhì)的程函方程，并轉(zhuǎn)換為走時(shí)的四次方程.為了求解這個(gè)方程，首先根據(jù)費(fèi)馬原理確定待求點(diǎn)可能的走時(shí)解范圍，然后根據(jù)程函方程求導(dǎo)以后的根將這個(gè)范圍劃分為一些只包含一個(gè)解的子區(qū)間.因此，可以用試位法在每個(gè)子區(qū)間尋找走時(shí)解.

本文提出一種可以將慢度方程退化為二次方程形式的計(jì)算qSV波和qSH波的初至走時(shí)的快速掃描算法.本研究通過(guò)選取恰當(dāng)?shù)木W(wǎng)格模板，發(fā)現(xiàn)待求點(diǎn)的慢度分量中有一個(gè)是已知的，可以通過(guò)相鄰點(diǎn)的走時(shí)近似計(jì)算得到.因此，qP波和qSV波解耦的慢度四次方程可以簡(jiǎn)化為二次方程解析求解.相比于傳統(tǒng)的慢度方程或各向異性程函方程轉(zhuǎn)換為走時(shí)四次方程(Han et al.,2017;Huang et al.,2020;Huang and Luo,2020)，本文提出的方法既包含了該類方法的優(yōu)勢(shì)，又極大地提高了計(jì)算效率.對(duì)于qSH波，它的慢度方程是二次的，可以直接用解析方法求解.

本文的結(jié)構(gòu)安排如下：首先推導(dǎo)出各向異性介質(zhì)中qSV波和qSH波的慢度方程.然后，在構(gòu)建的三角形網(wǎng)格局部解中，建立起慢度和走時(shí)之間的關(guān)系，把qSV波的慢度方程簡(jiǎn)化為二次方程解析地求解.qSH波的慢度方程是二次方程，可以解析求解.接著，給出因果性測(cè)試并總結(jié)用快速掃描算法計(jì)算二維VTI介質(zhì)中qSV波和qSH波走時(shí)的流程.最后，用兩個(gè)二維數(shù)值實(shí)驗(yàn)驗(yàn)證了本文方法在VTI模型上準(zhǔn)確計(jì)算走時(shí)的有效性，并給出了結(jié)論.

1 基本理論

本文提出的求解二次方程的快速掃描法組成如下：根據(jù)相鄰點(diǎn)的走時(shí)并利用慢度方程構(gòu)建局部解更新待求點(diǎn)的走時(shí)，通過(guò)求解群速度分量測(cè)試因果性.如果局部解無(wú)解，那么沿著固定方向更新走時(shí).最后，將上述步驟結(jié)合Gauss-Seidel迭代沿著可選方向進(jìn)行掃描更新走時(shí).

1.1 局部解構(gòu)造以及走時(shí)和慢度的關(guān)系

在二維VTI介質(zhì)中，從Christoffel方程導(dǎo)出的解耦qP波和qSV波的慢度四次方程有如下形式：

(1)

這里的ε，δ表示Thomsen(1986)參數(shù)，α0和β0表示對(duì)稱軸方向的P波和S波速度，px和pz分別表示水平方向和垂直方向的慢度分量.為了表示方便，公式(1)可以簡(jiǎn)化為如下形式：

(2)

(3)

SqP和SqSV分別表示qP波和qSV波的慢度方程.對(duì)于qSH波，慢度方程表示如下：

(4)

γ是橫波各向異性和橫波分裂強(qiáng)度的參數(shù).圖1展示了qP波、qSV波和qSH波的慢度曲面.慢度方程的詳細(xì)導(dǎo)出參見(jiàn)附錄A.

圖1 (a)qP波和qSV波的慢度曲面，內(nèi)側(cè)曲線代表qP波的，外側(cè)曲線代表qSV波的;(b)qSH波的慢度曲面Fig.1 (a)The coupled slowness surface of qP and qSV waves.The inner curve is the qP-wave slowness surface,and the outer one is the qSV-wave slowness surface;(b)The slowness surface of qSH wave

目前常用的兩種局部解的模型如圖2a和2b(Qian et al.,2007)所示.局部解是基于慢度方程用相鄰點(diǎn)的走時(shí)求解待求點(diǎn)的走時(shí)，也就是圖2所示中C點(diǎn)的走時(shí).為此，需要在待求點(diǎn)C建立走時(shí)和慢度的關(guān)系.首先以圖2a三角形單元0為例.在一個(gè)三角網(wǎng)格中，如果A點(diǎn)和B點(diǎn)的走時(shí)已知，CA和CB方向的慢度分量可以表示為：

(5)

這里的TA和TB分別表示A點(diǎn)和B點(diǎn)的走時(shí)，TC表示待求點(diǎn)C的走時(shí)，h表示網(wǎng)格間距.通過(guò)幾何關(guān)系，CA和CB方向的慢度也可以由水平和垂直方向的慢度P=(px,pz)映射得到：

(6)

其中Q是映射矩陣.聯(lián)合公式(5)和(6)，可以得到如下關(guān)系式：

(7)

這里

(8)

α=∠ACB=45°，如圖2a所示.其他三角形網(wǎng)格單元的局部解構(gòu)造類似于單元0.

圖2 二維網(wǎng)格示意圖(a)C點(diǎn)周圍有八個(gè)三角形單元；(b)C點(diǎn)周圍有四個(gè)三角形單元(Qian et al.,2007)；(c)模板(a)對(duì)應(yīng)的二階情況.P代表入射到C點(diǎn)的慢度向量.Fig.2 The mesh and stencils of the scheme in the 2D case(a)There are the eight-triangle stencils around point C;(b)There are the four-triangle stencils around point C (Qian et al.,2007);(c)Second order case for stencil (a).P is a slowness vector pointing to C.

再以圖2b中三角形單元0為例，如果A點(diǎn)和B點(diǎn)的走時(shí)已知，CA和CB方向的慢度分量和水平垂直方向的慢度分量相等，直接可以表示為：

(9)

通過(guò)以上分析構(gòu)建出了兩種常用局部解中走時(shí)和慢度的關(guān)系(公式(7)和(9)).

1.2 求解qSV波和qSH波的慢度方程

Han等(2017)使用圖2a中的局部解建立走時(shí)和慢度的關(guān)系，然后直接把公式(7)代入公式(1)得到走時(shí)的四次方程，并數(shù)值求解.Huang等(2020)、Huang和Luo(2020)使用圖2b中的局部解，通過(guò)將各向異性程函方程轉(zhuǎn)換為走時(shí)的四次方程，數(shù)值迭代求解.本研究發(fā)現(xiàn)，在第一種局部解中有水平慢度分量和垂直慢度分量只有一個(gè)未知這樣的特征，可以避免求解四次方程.先簡(jiǎn)化局部解公式(7)：

(10)

類似地，對(duì)于圖2a中局部解的二階格式，可以構(gòu)造如圖2c中的局部解，得到如下的慢度分量和走時(shí)關(guān)系：

(11)

那么，有

(12)

這里由于選取網(wǎng)格的對(duì)角點(diǎn)都為C點(diǎn)，水平或者豎直的是B點(diǎn)(圖2a所示).所以在每個(gè)三角形單元中必然有(xC-xB)或者(zC-zB)為0，這說(shuō)明不論一階還是二階格式，都有px和pz其中一個(gè)是已知的，可以由相鄰點(diǎn)的走時(shí)計(jì)算.還是以單元0為例，其中：xC-xB=0和zB-zA=0，代入公式(10)可以得到：

(13)

這里px為已知，所以公式(2)中的B和C也為已知，pz為待求，那么SqSV就簡(jiǎn)化為一個(gè)二次方程，直接可以解析求解如下：

(14)

公式(14)直接給出可能的根，而數(shù)值求解四次方程的方法需要判斷根的子區(qū)間和數(shù)值迭代求根(迭代次數(shù)從幾次到數(shù)十次)，所以在計(jì)算效率方面本文方法更高.本文只考慮初至波，因此TC的走時(shí)范圍為：

(15)

1.3 因果性測(cè)試和固定點(diǎn)方向的更新

如圖3所示，在各向異性介質(zhì)中，相速度方向(走時(shí)梯度)和群速度方向(射線方向)不再像在各向同性介質(zhì)中保持一致.因此，不能再用時(shí)間的梯度確定更新下一步的方向，只能使用群速度的方向.對(duì)于每個(gè)空間點(diǎn)，完成局部解的求解之后，還需要添加因果條件的測(cè)試確保得到的走時(shí)解要符合因果性條件，才可以更新待求點(diǎn)的走時(shí).求解公式(14)之后，由對(duì)應(yīng)的慢度向量計(jì)算得到如下的群速度向量：

圖3 二維均勻各向異性介質(zhì)中群速度(紅色線)和相速度(藍(lán)色線)的示意圖(Tsvankin，2005)θ表示相速度角度，ψ表示群速度角度.Fig.3 Group-velocity direction (red line)and phase-velocity direction (blue line)in a 2D anisotropic homogeneous media (Tsvankin,2005)θ denotes the phase-velocity angle,and ψ denotes the group-velocity angle.

(16)

圖4 因果性測(cè)試圖中vg是一個(gè)群速度向量，落在了三角形ABC中，那么對(duì)應(yīng)的走時(shí)就滿足因果性.Fig.4 Causality testThe vg is a group-velocity vector,which falls in triangular mesh ABC,and the related traveltime satisfies the causality.

當(dāng)三角形單元0的局部解無(wú)解或者不滿足因果性，假設(shè)地震波是直接通過(guò)線段AC或者線段BC傳播到C點(diǎn)，那么走時(shí)的更新公式如下：

(17)

對(duì)于qSV波的AC和BC固定方向的群速度求解，本文使用迭代方法.在二維VTI介質(zhì)中，相速度和群速度之間的關(guān)系可以表示為(Berryman,1979)：

(18)

k表示波數(shù)，其中的兩項(xiàng)可以分別展開(kāi)為

(19)

在公式(19)中用ψ表示群速度的角度，θ表示相速度的角度，V(θ)表示對(duì)應(yīng)相角下的相速度，可以由Christoffel方程求得.借助于(18)式，可以得到如下公式：

(20)

因?yàn)橹磺?5°,135°，225°,275°這些固定方向的群速度，直接將這些群角代入公式(20)數(shù)值求解.例如求解45°方向群速度，tan45°=1，選取相速度的區(qū)間為[0°,89°]，如圖5所示，可以用數(shù)值方法尋找相角使得FG小于給定的閾值，對(duì)應(yīng)的相角再代回公式(19)求解對(duì)應(yīng)的群速度.對(duì)于存在多個(gè)解的情況(圖5b)，只選擇最小的群速度即可.

圖5 公式(20)在群速度等于45°的曲線(a)Thomsen參數(shù)為：α0=3 km·s-1,β0=1.5 km·s-1,ε=0.3,和δ=0.1.這種情況下，可以看到計(jì)算得到的相速度角度近似地為40°，可以直接利用公式(19)計(jì)算群速度角度.(b)Thomsen參數(shù)為：α0=3 km·s-1,β0=1.5 km·s-1,ε=0.3,和δ=-0.1.這種情況下，可以看到計(jì)算得到的相速度角度近似地為25°和60°，我們選擇利用公式(19)計(jì)算得到的最小的群速度.Fig.5 The curves of Eq.(20)in the group direction of 45°(a)The Thomsen parameters are:α0=3 km·s-1,β0=1.5 km·s-1,ε=0.3,and δ=0.1.In this case,we can see that the calculated phase angle is approximately 40°.Therefore,we can compute the group velocity using the calculated phase angle by Eq.(19).(b)The Thomsen parameters are:α0=3 km·s-1,β0=1.5 km·s-1,ε=0.3,and δ=-0.1.In this case,we can see that the calculated phase angles are approximately 25° and 60°.We choose the one corresponding to the minimal group velocity by Eq.(19).

對(duì)于qSH波，本文使用如下的公式求固定方向的群速度(Tsvankin,2005):

(21)

2 算法流程

本文提出的快速掃描法計(jì)算qSV波和qSH波走時(shí)的流程如下：

1)對(duì)震源進(jìn)行初始化，然后將所有的空間網(wǎng)格點(diǎn)賦予一個(gè)比計(jì)算走時(shí)大很多的大值.

2)Gauss-Seidel 迭代

a)根據(jù)每一次的掃描方向應(yīng)用局部解求解可能的走時(shí).掃描方向如下：

i=1∶nx,j=1∶nz;

i=nx∶1,j=1∶nz;

i=1∶nx,j=nz∶1;

i=nx∶1,j=nz∶1.

b)對(duì)于qSV波，驗(yàn)證求得的走時(shí)是否與 qSV 波相關(guān)，并測(cè)試因果性.qSH波過(guò)程類似.

c)如果局部解無(wú)解或者求出的解不滿足因果性，那么用固定方向更新走時(shí).

對(duì)于每個(gè)空間點(diǎn)，在掃描完所有方向以后，選擇其中最小的走時(shí)作為計(jì)算結(jié)果.

3)計(jì)算所有網(wǎng)格點(diǎn)的本次迭代與上一次迭代走時(shí)的誤差，從中選取最大的誤差，判斷誤差是否到達(dá)閾值.如果沒(méi)有，那么繼續(xù) Gauss-Seidel 迭代，直到誤差達(dá)到設(shè)置的閾值為止.

3 數(shù)值算例

在這部分中，我們給出兩個(gè)二維的算例.首先，用本文提出的方法計(jì)算走時(shí)，然后用波動(dòng)方程計(jì)算qSV波前.選取一定的時(shí)刻對(duì)比二者的匹配性，以此驗(yàn)證本文方法計(jì)算的走時(shí)的正確性.因?yàn)槟壳皼](méi)有二維qSH波的模擬程序，不做波前和走時(shí)的對(duì)比，只用均勻模型的解析解和數(shù)值解做對(duì)比.

在本算例中，選取均勻的VTI介質(zhì)，P波的速度為3000 m·s-1，S波的速度為2000 m·s-1.模型的網(wǎng)格橫向大小為301，縱向大小為301，網(wǎng)格間距為10 m，震源位于(1.5，1.5)km 處.使用兩組各向異性參數(shù)測(cè)試.第一組：各向異性參數(shù)為ε=0.3，δ=0.1.計(jì)算結(jié)果如圖6a所示，選取t=0.35 s時(shí)刻，抽取對(duì)應(yīng)的走時(shí)等時(shí)線和波動(dòng)方程模擬的波場(chǎng)進(jìn)行匹配，圖中的結(jié)果表明兩者重合較好.為了定量評(píng)估誤差，圖6b展示了解析解和數(shù)值解的對(duì)比，圖中的藍(lán)色曲線表示解析解得到的走時(shí)結(jié)果，紅色曲線表示本文方法數(shù)值求解得到的走時(shí)結(jié)果，可以觀察到兩者吻合.第二組：介質(zhì)各向異性參數(shù)為ε=0.3，δ=-0.1.計(jì)算的結(jié)果如圖6c，選取t=0.35 s時(shí)刻，抽取對(duì)應(yīng)的走時(shí)等時(shí)線和波動(dòng)方程模擬的波場(chǎng)進(jìn)行匹配.對(duì)于該組參數(shù)，從波動(dòng)方程模擬的結(jié)果看到，在特定的方向有qSV波交叉重疊的現(xiàn)象.對(duì)應(yīng)的群速度的求解存在多個(gè)值，本文在求解的過(guò)程中只選取最小的值，可以從圖6中觀察到這些特定方向的走時(shí)和波前匹配有誤差，其他地方的走時(shí)結(jié)果和波前重合.兩組參數(shù)的計(jì)算結(jié)果驗(yàn)證了本文方法在均勻介質(zhì)準(zhǔn)確求解qSV波走時(shí)的有效性.對(duì)于qSH波，各向異性參數(shù)設(shè)置為γ=0.2；圖7展示了VTI介質(zhì)中的走時(shí)計(jì)算結(jié)果.圖7中的藍(lán)色曲線表示解析解得到的走時(shí)結(jié)果，紅色曲線表示本文方法數(shù)值求解得到的走時(shí)結(jié)果，可以觀察到兩者吻合，驗(yàn)證了該方法求解qSH波走時(shí)的有效性.

圖6 (a)ε=0.3,δ=0.1下的走時(shí)和波前的匹配；(b)數(shù)值解和解析解的對(duì)比；(c)ε=0.3,δ=-0.1下的走時(shí)和波前的匹配Fig.6 (a)The traveltime contours overlaid on the corresponding snapshots with anisotropic parameters:ε=0.3,δ=0.1;(b)Comparison of the numerical and analytical solutions;(c)The traveltime contours overlaid on the corresponding snapshots with anisotropic parameters ε=0.3,δ=-0.1

圖7 qSH 波的走時(shí)計(jì)算結(jié)果Fig.7 Traveltime results of qSH wave

再將本方法應(yīng)用于抽取的部分BP模型測(cè)試算法的有效性.圖8展示了模型的P波速度、參數(shù)ε和δ.設(shè)置S波的速度為P波速度除以1.5.該模型大小為橫向1200個(gè)網(wǎng)格點(diǎn)，縱向901個(gè)網(wǎng)格點(diǎn)，網(wǎng)格間距為10 m，震源位于(6,4.5)km.圖9a和9b展示了在VTI情況下t=0.48 s和t=0.98 s時(shí)刻計(jì)算得到的qSV波走時(shí)和波前的匹配，從圖中可以觀察到這兩個(gè)時(shí)刻計(jì)算得到的走時(shí)(紅色等時(shí)線)和對(duì)應(yīng)時(shí)刻的波場(chǎng)模擬得到的波前吻合，驗(yàn)證了本文方法在VTI介質(zhì)計(jì)算走時(shí)的準(zhǔn)確性.表1是計(jì)算效率的對(duì)比，可以看出，達(dá)到相同的誤差，原始四次方程迭代求解需要的快速掃描計(jì)算迭代次數(shù)比本文方法多一次，計(jì)算時(shí)間約為本文方法的5倍.圖9c展示了本文方法和數(shù)值求解四次方程方法(Han et al.,2017)的計(jì)算結(jié)果對(duì)比，二者一致，因?yàn)閮烧叨际窃伎刂扑拇温确匠痰慕?本算例驗(yàn)證了本文方法在VTI介質(zhì)計(jì)算走時(shí)的高效性和準(zhǔn)確性.圖10展示了在VTI情況下計(jì)算得到的qSH波走時(shí)，這里設(shè)置γ參數(shù)等于ε.圖10中包含了原始網(wǎng)格(紅色虛線)計(jì)算結(jié)果和加密1倍網(wǎng)格(藍(lán)色實(shí)線)計(jì)算結(jié)果的對(duì)比，觀察到二者吻合.本算例說(shuō)明即使對(duì)于BP這樣復(fù)雜的各向異性模型，本文提出的方法也得到了較好的效果.

圖8 BP模型參數(shù)(a)P波速度模型;(b)ε模型;(c)δ模型.Fig.8 BP model parameters(a)P-wave velocity model;(b)ε model;(c)δ model.

圖9 BP 模型的VTI情況(a)和(b)對(duì)應(yīng)著不同時(shí)刻的波前和走時(shí)匹配；(c)本文方法和迭代求解四次方程方法(Han et al.,2017)結(jié)果對(duì)比圖.Fig.9 VTI case for the BP model(a)and (b)present the traveltime contours overlaid on the corresponding snapshots at different times;(c)Comparison of the results using the proposed method and the iterative method for solving quartic equation (Han et al.,2017).

圖10 qSH波的加密網(wǎng)格和原網(wǎng)格計(jì)算結(jié)果對(duì)比圖Fig.10 The traveltime comparisons of qSH wave for the encrypted and original grids

表1 BP模型的VTI情況下兩種方法的對(duì)比Table 1 The comparisons between the two methods in VTI for the BP model

4 結(jié)論

在各向異性介質(zhì)中計(jì)算初至波的走時(shí)十分重要.快速掃描法不需要存儲(chǔ)和追蹤波前面信息，在各向異性介質(zhì)初至走時(shí)計(jì)算中有著重要的應(yīng)用.基于求解慢度四次方程轉(zhuǎn)換為走時(shí)四次方程的快速掃描法適用于強(qiáng)各向異性介質(zhì)，但是存在計(jì)算效率低的問(wèn)題.本文發(fā)展了一種求解二次慢度方程的快速掃描法(FSM)計(jì)算二維VTI介質(zhì)中qSV波和qSH波的初至走時(shí).首先，推導(dǎo)了qSV波和qSH波的慢度方程，并在構(gòu)建的局部解中建立起走時(shí)和慢度的關(guān)系.接著，構(gòu)建了一種八個(gè)三角形的空間局部解模板，并推導(dǎo)得到在模板中的每個(gè)三角形單元局部解中水平或者垂直慢度分量中只有一個(gè)是未知的.因此，對(duì)于qSV波可以將解耦的慢度四次方程簡(jiǎn)化為二次方程解析求解，既包含了求解四次慢度方程方法的優(yōu)勢(shì)，又極大地提高了計(jì)算效率.求出慢度以后，再利用走時(shí)和慢度的關(guān)系得到走時(shí).之后，對(duì)局部解求出的走時(shí)進(jìn)行因果性測(cè)試，如果不滿足因果性或者無(wú)解，那么進(jìn)行固定方向的走時(shí)更新.對(duì)于qSH波，由于慢度方程是二次方程，解析求解走時(shí)的思路以及因果性測(cè)試和qSV波類似.最后，總結(jié)了本文提出的方法的工作流程，并用數(shù)值試驗(yàn)驗(yàn)證了該方法在均勻和復(fù)雜各向異性介質(zhì)中的有效性.

致謝感謝《地球物理學(xué)報(bào)》編輯部以及兩位匿名審稿人細(xì)致的評(píng)審，對(duì)提升本文的質(zhì)量有很大的幫助.感謝深圳市深遠(yuǎn)海油氣勘探技術(shù)重點(diǎn)實(shí)驗(yàn)室(項(xiàng)目編號(hào)：ZDSYS20190902093007855)和深圳市科技計(jì)劃(項(xiàng)目編號(hào)：KQTD20170810111725321)資助.

附錄A

(A1)

這里的n=(nx,ny,nz)代表波的傳播方向.在本研究中只考慮二維平面[x,z]，因此：ny=0，大括號(hào)里面的第一項(xiàng)對(duì)應(yīng)的是qSH波，第二項(xiàng)對(duì)應(yīng)的是解耦的qP波和qSV波，這兩項(xiàng)是相互獨(dú)立的，做如下的定義：

(A2)

代入公式(A1)中，可以得到如下的qP波和qSV波的解耦四次慢度方程：

(A3)

對(duì)于2D中的qSH波有：

(A4)