求解二次方程的二維VTI介質(zhì)qSV波和qSH波走時(shí)快速掃描算法

蘆永明，張偉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).它對地震波傳播的振幅和走時(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)，簡單應(yīng)用各向同性走時(shí)計(jì)算方法有可能引起違背時(shí)間因果性的問題，使得求解初至走時(shí)變得比較困難.

對于各向異性介質(zhì)，快速行進(jìn)法只能解決一些簡單的情況(Bin Waheed and Alkhalifah，2017)，并且這類方法在各向異性介質(zhì)中，如果從波前上最小走時(shí)點(diǎn)擴(kuò)展波前面，可能會(huì)違背時(shí)間的因果條件.而快速掃描法不需要存儲(chǔ)和追蹤波前，可以克服該問題，在求解向異性介質(zhì)中的走時(shí)方面得到了廣泛的應(yīng)用.Kao等(2004)用 Lax-Friedrichs 差分格式求解VTI介質(zhì)中的哈密爾頓雅可比方程.對該方法，人工黏性參數(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)勢是可以擴(kuò)展到方位各向異性介質(zhì)中，但是求解需要更多的迭代次數(shù).Le Bouteiller等(2018)用快速掃描法求解二維TI介質(zhì)中qP波的走時(shí)，他們用間斷伽遼金方法求解時(shí)間依賴的哈密爾頓雅可比方程.為了使快速掃描法適用于強(qiáng)各向異性介質(zhì)，Han等(2017)發(fā)展了基于慢度四次方程的快速掃描法計(jì)算二維TTI介質(zhì)中qP波的走時(shí).他們通過構(gòu)造局部解把qP波和qSV波的解耦慢度四次方程轉(zhuǎn)換成走時(shí)的四次方程數(shù)值求解.對走時(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í)的快速掃描算法.本研究通過選取恰當(dāng)?shù)木W(wǎng)格模板，發(fā)現(xiàn)待求點(diǎn)的慢度分量中有一個(gè)是已知的，可以通過相鄰點(diǎn)的走時(shí)近似計(jì)算得到.因此，qP波和qSV波解耦的慢度四次方程可以簡化為二次方程解析求解.相比于傳統(tǒng)的慢度方程或各向異性程函方程轉(zhuǎn)換為走時(shí)四次方程(Han et al.,2017;Huang et al.,2020;Huang and Luo,2020)，本文提出的方法既包含了該類方法的優(yōu)勢，又極大地提高了計(jì)算效率.對于qSH波，它的慢度方程是二次的，可以直接用解析方法求解.

本文的結(jié)構(gòu)安排如下：首先推導(dǎo)出各向異性介質(zhì)中qSV波和qSH波的慢度方程.然后，在構(gòu)建的三角形網(wǎng)格局部解中，建立起慢度和走時(shí)之間的關(guān)系，把qSV波的慢度方程簡化為二次方程解析地求解.qSH波的慢度方程是二次方程，可以解析求解.接著，給出因果性測試并總結(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í)，通過求解群速度分量測試因果性.如果局部解無解，那么沿著固定方向更新走時(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表示對稱軸方向的P波和S波速度，px和pz分別表示水平方向和垂直方向的慢度分量.為了表示方便，公式(1)可以簡化為如下形式：

(2)

(3)

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

(4)

γ是橫波各向異性和橫波分裂強(qiáng)度的參數(shù).圖1展示了qP波、qSV波和qSH波的慢度曲面.慢度方程的詳細(xì)導(dǎo)出參見附錄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ā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)對應(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)

通過以上分析構(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中的局部解，通過將各向異性程函方程轉(zhuǎn)換為走時(shí)的四次方程，數(shù)值迭代求解.本研究發(fā)現(xiàn)，在第一種局部解中有水平慢度分量和垂直慢度分量只有一個(gè)未知這樣的特征，可以避免求解四次方程.先簡化局部解公式(7)：

(10)

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

(11)

那么，有

(12)

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

(13)

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

(14)

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

(15)

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

如圖3所示，在各向異性介質(zhì)中，相速度方向(走時(shí)梯度)和群速度方向(射線方向)不再像在各向同性介質(zhì)中保持一致.因此，不能再用時(shí)間的梯度確定更新下一步的方向，只能使用群速度的方向.對于每個(gè)空間點(diǎn)，完成局部解的求解之后，還需要添加因果條件的測試確保得到的走時(shí)解要符合因果性條件，才可以更新待求點(diǎn)的走時(shí).求解公式(14)之后，由對應(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 因果性測試圖中vg是一個(gè)群速度向量，落在了三角形ABC中，那么對應(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的局部解無解或者不滿足因果性，假設(shè)地震波是直接通過線段AC或者線段BC傳播到C點(diǎn)，那么走時(shí)的更新公式如下：

(17)

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

(18)

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

(19)

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

(20)

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

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

(21)

2 算法流程

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

1)對震源進(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)對于qSV波，驗(yàn)證求得的走時(shí)是否與 qSV 波相關(guān)，并測試因果性.qSH波過程類似.

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

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

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

3 數(shù)值算例

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

在本算例中，選取均勻的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ù)測試.第一組：各向異性參數(shù)為ε=0.3，δ=0.1.計(jì)算結(jié)果如圖6a所示，選取t=0.35 s時(shí)刻，抽取對應(yīng)的走時(shí)等時(shí)線和波動(dòng)方程模擬的波場進(jìn)行匹配，圖中的結(jié)果表明兩者重合較好.為了定量評估誤差，圖6b展示了解析解和數(shù)值解的對比，圖中的藍(lán)色曲線表示解析解得到的走時(shí)結(jié)果，紅色曲線表示本文方法數(shù)值求解得到的走時(shí)結(jié)果，可以觀察到兩者吻合.第二組：介質(zhì)各向異性參數(shù)為ε=0.3，δ=-0.1.計(jì)算的結(jié)果如圖6c，選取t=0.35 s時(shí)刻，抽取對應(yīng)的走時(shí)等時(shí)線和波動(dòng)方程模擬的波場進(jìn)行匹配.對于該組參數(shù)，從波動(dòng)方程模擬的結(jié)果看到，在特定的方向有qSV波交叉重疊的現(xiàn)象.對應(yīng)的群速度的求解存在多個(gè)值，本文在求解的過程中只選取最小的值，可以從圖6中觀察到這些特定方向的走時(shí)和波前匹配有誤差，其他地方的走時(shí)結(jié)果和波前重合.兩組參數(shù)的計(jì)算結(jié)果驗(yàn)證了本文方法在均勻介質(zhì)準(zhǔn)確求解qSV波走時(shí)的有效性.對于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ù)值解和解析解的對比；(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模型測試算法的有效性.圖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í)線)和對應(yīng)時(shí)刻的波場模擬得到的波前吻合，驗(yàn)證了本文方法在VTI介質(zhì)計(jì)算走時(shí)的準(zhǔn)確性.表1是計(jì)算效率的對比，可以看出，達(dá)到相同的誤差，原始四次方程迭代求解需要的快速掃描計(jì)算迭代次數(shù)比本文方法多一次，計(jì)算時(shí)間約為本文方法的5倍.圖9c展示了本文方法和數(shù)值求解四次方程方法(Han et al.,2017)的計(jì)算結(jié)果對比，二者一致，因?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é)果的對比，觀察到二者吻合.本算例說明即使對于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)對應(yīng)著不同時(shí)刻的波前和走時(shí)匹配；(c)本文方法和迭代求解四次方程方法(Han et al.,2017)結(jié)果對比圖.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é)果對比圖Fig.10 The traveltime comparisons of qSH wave for the encrypted and original grids

表1 BP模型的VTI情況下兩種方法的對比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ì)算效率低的問題.本文發(fā)展了一種求解二次慢度方程的快速掃描法(FSM)計(jì)算二維VTI介質(zhì)中qSV波和qSH波的初至走時(shí).首先，推導(dǎo)了qSV波和qSH波的慢度方程，并在構(gòu)建的局部解中建立起走時(shí)和慢度的關(guān)系.接著，構(gòu)建了一種八個(gè)三角形的空間局部解模板，并推導(dǎo)得到在模板中的每個(gè)三角形單元局部解中水平或者垂直慢度分量中只有一個(gè)是未知的.因此，對于qSV波可以將解耦的慢度四次方程簡化為二次方程解析求解，既包含了求解四次慢度方程方法的優(yōu)勢，又極大地提高了計(jì)算效率.求出慢度以后，再利用走時(shí)和慢度的關(guān)系得到走時(shí).之后，對局部解求出的走時(shí)進(jìn)行因果性測試，如果不滿足因果性或者無解，那么進(jìn)行固定方向的走時(shí)更新.對于qSH波，由于慢度方程是二次方程，解析求解走時(shí)的思路以及因果性測試和qSV波類似.最后，總結(jié)了本文提出的方法的工作流程，并用數(shù)值試驗(yàn)驗(yàn)證了該方法在均勻和復(fù)雜各向異性介質(zhì)中的有效性.

致謝感謝《地球物理學(xué)報(bào)》編輯部以及兩位匿名審稿人細(xì)致的評審，對提升本文的質(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)對應(yīng)的是qSH波，第二項(xiàng)對應(yīng)的是解耦的qP波和qSV波，這兩項(xiàng)是相互獨(dú)立的，做如下的定義：

(A2)

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

(A3)

對于2D中的qSH波有：

(A4)