基于Chebyshev走時(shí)逼近的三維多次反射射線計(jì)算

孫建國，苗 賀

吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，長春 130026

0 引言

射線走時(shí)和射線路徑是構(gòu)成零階幾何射線場(亦稱為幾何光學(xué)場)的兩個(gè)基本要素。其中，射線走時(shí)對應(yīng)高頻波場的相位，射線路徑對應(yīng)高頻波場的能量傳播軌跡。因此，自從射線理論被引入到地震學(xué)中以來，對射線走時(shí)和射線路徑計(jì)算方法的研究就一直沒有間斷過，尤其是在20世紀(jì)80年代以后，在反射地震偏移和反演研究的推動下[1-2]，對射線走時(shí)和射線路徑的研究一度成為研究熱點(diǎn)，有很多研究者和工程師參與到其中，提出了很多方法和技術(shù)[3-9]。形成這種局面的主要原因是基于射線理論的偏移和反演方法具有對輸入數(shù)據(jù)體觀測裝置適應(yīng)性強(qiáng)和計(jì)算速度快等優(yōu)點(diǎn)[1-2]。再有，在偏移成像中的射線追蹤與常規(guī)的兩點(diǎn)射線追蹤有所不同，只須設(shè)定射線起點(diǎn)的位置，而不須設(shè)定射線終點(diǎn)的位置。這是因?yàn)榛谏渚€理論的偏移在具體實(shí)施上是個(gè)一點(diǎn)對多點(diǎn)過程，任意一個(gè)成像點(diǎn)上的射線走時(shí)和射線方向均可利用相鄰點(diǎn)上的相關(guān)信息通過插值得到。然而，由于這種只考慮初值的射線追蹤方法不能在計(jì)算過程中隨時(shí)插入由非均勻介質(zhì)中的聚焦和反聚焦現(xiàn)象所產(chǎn)生的新射線(散射射線)，從而會在成像區(qū)域形成射線覆蓋盲區(qū)。因此，在20世紀(jì)90年代以后，無論是工業(yè)界還是學(xué)術(shù)界均傾向于采用波前構(gòu)建法來替代常規(guī)的初值射線追蹤法，并逐步建立了基于波前構(gòu)建的Kirchhoff型偏移和波束偏移[1-2]。但是，波前構(gòu)建法在實(shí)現(xiàn)時(shí)需要利用復(fù)雜的數(shù)據(jù)結(jié)構(gòu)，而且對于多路徑走時(shí)的存儲也比較復(fù)雜[10-14]。所以，在20世紀(jì)90年代以后，基于程函方程數(shù)值分析的射線走時(shí)計(jì)算方法逐步受到關(guān)注和重視[1-5]。

在歷史上，利用程函方程數(shù)值解計(jì)算地震波走時(shí)的基本思想早已有之，但是由于沒有相應(yīng)的需求而一直不為人所關(guān)注。1988年，Vidale首先將前人的基本思想付諸實(shí)踐[3-4]，提出了基于程函方程有限差分(FD)數(shù)值解的地震波走時(shí)計(jì)算方法。與上文提到的初值方法相比，F(xiàn)D法具有計(jì)算速度快和沒有計(jì)算盲區(qū)等特點(diǎn)，從而可以作為一種理想的地震波走時(shí)計(jì)算方法替代傳統(tǒng)的射線追蹤法。然而，在無人為干預(yù)時(shí)，利用有限差分等基于網(wǎng)格的程函方程數(shù)值分析方法只能沿著波的前進(jìn)方向計(jì)算走時(shí)，而不能沿著波的回折方向計(jì)算走時(shí)。換句話說，在無人為干預(yù)時(shí)，利用任何一種基于程函方程數(shù)值分析的方法均只能計(jì)算沿著波的前進(jìn)方向的透射(波)走時(shí)，而不能同時(shí)計(jì)算由該透射波所激發(fā)的反射(波)走時(shí)，更不能計(jì)算由該透射波所激發(fā)的、其傳播方向經(jīng)過多次改變的多次透射(波)和多次反射(波)走時(shí)。針對這個(gè)問題，Rawlinson等[5]在2004年提出了一種在人為干預(yù)下利用程函方程數(shù)值分析方法計(jì)算多次透射和多次反射走時(shí)的分區(qū)多級方案，在這種方案中，走時(shí)計(jì)算采用了基于窄帶技術(shù)的快速推進(jìn)法(fast marching method，F(xiàn)MM)[5]。為計(jì)算多次反射走時(shí)，Rawlinson等[5]將每一個(gè)帶有光滑彎曲界面的非均勻?qū)右暈橐粋€(gè)獨(dú)立的計(jì)算單元，而每一個(gè)層面均被視為窄帶(層面窄帶)。在層面窄帶中，快速推進(jìn)過程被反復(fù)重置為初始狀態(tài)(重新初始化)，并根據(jù)預(yù)定的傳播方向重新進(jìn)行推進(jìn)計(jì)算。為了計(jì)算射線路徑，Rawlinson等[5]采用了逆向梯度法，即首先利用有限差分網(wǎng)格節(jié)點(diǎn)上的走時(shí)值計(jì)算走時(shí)梯度，然后再根據(jù)梯度方向計(jì)算射線路徑的局部坐標(biāo)，最后將所有的局部坐標(biāo)點(diǎn)連成一體，得到射線路徑。

2006年，de Kool等[6]將上述計(jì)算方法推廣到三維。2017年，孫章慶等[7]將這種多級計(jì)算方法與群推進(jìn)法組合在一起，提出了一種針對復(fù)雜山地多波型走時(shí)計(jì)算的多級次群推進(jìn)迎風(fēng)混合法。類似地，唐小平等[8]將多級計(jì)算方法與最短路徑算法組合在一起，提出了一種基于最短路徑法的、針對三維層狀介質(zhì)中多次波的走時(shí)與射線路徑計(jì)算方法。在理論上，分區(qū)多級計(jì)算的基本思想可以用于任何一種基于程函方程數(shù)值分析的地震波走時(shí)計(jì)算方法，因此存在有無限多的組合方式。

無論是在Rawlinson等的開創(chuàng)性工作中，還是在后續(xù)的跟蹤完善性工作中，走時(shí)梯度都是直接利用網(wǎng)格點(diǎn)走時(shí)值來計(jì)算的。由此而產(chǎn)生的問題是：①梯度計(jì)算的精度和穩(wěn)定性受有限差分網(wǎng)格常數(shù)的制約，并且隨計(jì)算區(qū)域和網(wǎng)格密度的變化而變化，尤其是在復(fù)雜(曲)邊界附近或是在變網(wǎng)格中更是如此；②梯度計(jì)算的精度和穩(wěn)定性受單個(gè)網(wǎng)格點(diǎn)走時(shí)精度的制約。為解決這個(gè)問題，張東等[15]采用B樣條插值多項(xiàng)式，將對離散數(shù)據(jù)的數(shù)值求導(dǎo)問題轉(zhuǎn)化為對分片光滑函數(shù)的求導(dǎo)問題。2014年，張婷婷等[16]將B樣條插值用于三維層狀介質(zhì)中的多次波走時(shí)和多次波射線路徑計(jì)算，為在三維條件下研究多次反射提供了一個(gè)新的計(jì)算途徑。然而，根據(jù)張婷婷等[17]的計(jì)算結(jié)果，基于B樣條插值多項(xiàng)式的計(jì)算方法會在速度突變區(qū)域內(nèi)產(chǎn)生一定的誤差，因而難以在速度劇烈變化區(qū)域內(nèi)得到滿意的計(jì)算結(jié)果。本文針對這個(gè)問題：①利用Chebyshev正交多項(xiàng)式逼近走時(shí)場，使分區(qū)多級計(jì)算的誤差在最小二乘的意義下達(dá)到最小，從而避開速度突變所帶來的走時(shí)突變問題；②對Chebyshev逼近函數(shù)求導(dǎo)，得到射線路徑在所考慮區(qū)域內(nèi)的局部坐標(biāo)。

1 基于程函方程有限差分解的分區(qū)多級走時(shí)計(jì)算

根據(jù)Rawlinson等[5]的工作，分區(qū)多級計(jì)算的主要步驟是：①在模型中劃分出反射區(qū)和透射區(qū)。②將反射(透射)區(qū)內(nèi)的網(wǎng)格點(diǎn)劃分為接收點(diǎn)、活動點(diǎn)和遠(yuǎn)離點(diǎn)(圖1)，其中：接收點(diǎn)位于窄帶以外的上風(fēng)區(qū)，是已經(jīng)完成走時(shí)計(jì)算的點(diǎn)；活動點(diǎn)位于窄帶之內(nèi)，是當(dāng)前正在進(jìn)行走時(shí)計(jì)算的點(diǎn)；遠(yuǎn)離點(diǎn)位于窄帶以外的下風(fēng)區(qū)，是還沒有進(jìn)行走時(shí)計(jì)算的點(diǎn)。③利用公式

(1)

計(jì)算所考慮活動點(diǎn)上走時(shí)梯度的模。式中：Si,j,k為網(wǎng)格慢度；|t|i,j,k為走時(shí)梯度的模；和分別表示點(diǎn)(i,j,k)處走時(shí)在對應(yīng)方向上的n階向前和向后差分算子。當(dāng)n=1時(shí)，

(2)

式中，h為網(wǎng)格間距。④將最小走時(shí)值賦予所考慮的活動點(diǎn)，并將其標(biāo)記為接收點(diǎn)，移出窄帶。⑤將距所考慮活動點(diǎn)最近的遠(yuǎn)離點(diǎn)移入窄帶，重復(fù)③—④所描述的計(jì)算過程，直到全部遠(yuǎn)離點(diǎn)均變?yōu)榻邮拯c(diǎn)。

圖2給出了上述分區(qū)多級計(jì)算步驟的進(jìn)一步說明。①首先由震源開始，計(jì)算反射(透射)面及以上區(qū)域內(nèi)所有網(wǎng)格點(diǎn)處的走時(shí)(圖2a)。②將界面(反射(透射)面)重置為窄帶(圖2b)。③為計(jì)算反射波，將界面(反射(透射)面)以上的點(diǎn)設(shè)為遠(yuǎn)離點(diǎn)，按照上文中的計(jì)算步驟重啟FMM(圖2c)。為計(jì)算透射波，將界面(反射(透射)面)以下層內(nèi)的點(diǎn)設(shè)為遠(yuǎn)離點(diǎn)，按照上文中的計(jì)算步驟重啟FMM(圖2d)。重復(fù)圖2所描述的計(jì)算過程，即可得到任意級次的多次反射和多次透射。

圖1 反射(透射)區(qū)內(nèi)的接收點(diǎn)(黑色圓點(diǎn))、活動點(diǎn)(灰色圓點(diǎn))和遠(yuǎn)離點(diǎn)(白色圓點(diǎn))Fig.1 receiver point (black), active point (gray) and far point (white) in reflection (transmission) region

據(jù)文獻(xiàn)[4]修改。圖2 分區(qū)多級走時(shí)計(jì)算過程示意圖Fig.2 Partition multi-level travel time calculation process

2 基于Chebyshev多項(xiàng)式的離散走時(shí)逼近

Chebyshev逼近是基于Chebshev多項(xiàng)式的一種逼近方法，是為解決零偏差問題而提出的一種基于正交多項(xiàng)式的最佳逼近方法[18]。在歷史上，Chebyshev逼近和Chebshev多項(xiàng)式由俄羅斯數(shù)學(xué)家P L Chebyshev在1854年提出，其最初目的是解決連桿設(shè)計(jì)問題[19]。根據(jù)其數(shù)學(xué)表現(xiàn)形式可將Chebyshev提出的多項(xiàng)式分為兩種，即第一類和第二類Chebyshev多項(xiàng)式。其中，第一類Chebyshev多項(xiàng)式起源于多倍角余弦和正弦函數(shù)展開，用Tn(x)表示。具體地講，Tn(x)=cos[ncos-1(x)]。第二類Chebyshev多項(xiàng)式用Un(x)表示，與Tn(x)的關(guān)系為Un(x)=[dTn+1(x)/dx]/(n+1)。在實(shí)際應(yīng)用中，利用下列遞推公式計(jì)算Tn(x)和Un(x)，即T0(x)=1、T1(x)=x、Tn+1(x)=2xTn(x)-Tn-1(x)、U0(x)=1、U1(x)=2x、Un+1(x)=2xUn(x)-Un-1(x)。

作為一種正交多項(xiàng)式，Chebyshev多項(xiàng)式與其他正交多項(xiàng)式相比，例如與Legendre、Laguerre和Hermite等正交多項(xiàng)式相比具有許多獨(dú)特的性質(zhì)[9]。事實(shí)上，Chebyshev多項(xiàng)式特別適合于進(jìn)行最佳平方逼近。此外，在一定條件下可將連續(xù)函數(shù)展開為Chebyshev級數(shù)等等。因此，為在最小平方意義下逼近來自于快速推進(jìn)法的走時(shí)場，在下文中選取Chebyshev多項(xiàng)式為逼近函數(shù)。如圖3所示，在諸多正交多項(xiàng)式中，Chebyshev多項(xiàng)式的方差最小，即逼近誤差亦為最小。

現(xiàn)在利用第一類Chebyshev多項(xiàng)式Tn(x)來構(gòu)造3維離散走時(shí)τ(x,y,z)的近似解析表達(dá)式t(x,y,z)。為此目的，引入系數(shù)序列a′m，m=0,1,…,8，使得

(3)

展開后得到

t(x,y,z)=a0+a1z+a2(2z2-1)+a3y+

a4yz+a5y(2z2-1)+a6(2y2-1)+

a7z(2y2-1)+a8(2y2-1)(2z2-1)+

a9x+a10xz+a11x(2z2-1)+a12xy+

a13xyz+a14xy(2z2-1)+a15x(2y2-1)+

a16(2y2-1)xz+a17x(2y2-1)(2z2-1)+

a18(2x2-1)+a19(2x2-1)z+

a20(2x2-1)(2z2-1)+

a21(2x2-1)(2y2-1)+

a22(2x2-1)y+a23(2x2-1)yz+

圖3 三種數(shù)據(jù)4種多項(xiàng)式逼近效果(a、b、c)及方差對比圖(d)Fig.3 Four polynomial approximations for three kinds of data (a, b, c) and variance comparison chart (d)

a25(2x2-1)(2y2-1)z+

a26(2x2-1)(2y2-1)(2z2-1)。

(4)

式中，ap(p∈0,1,2,…,26)是利用最小二乘法得到的插值系數(shù)，即

(5)

式中：τ(xi,yj,zk)為網(wǎng)格節(jié)點(diǎn)上的走時(shí)值；φp(xi,yj,zk)(p∈0,1,2,…,26)為逼近函數(shù)中與系數(shù)對應(yīng)的正交項(xiàng)；M,N,L分別為x,y,z方向的網(wǎng)格節(jié)點(diǎn)數(shù)。

一旦完成了公式(4)所描述的逼近過程，即可根據(jù)走時(shí)函數(shù)t(x,y,z)的梯度函數(shù)t(x,y,z)求取當(dāng)前追蹤點(diǎn)的射線方向。根據(jù)定義，

t(x,y,z)=

(6)

3 數(shù)值計(jì)算實(shí)例

3.1 常梯度速度模型

首先考慮圖4中所示的常梯度速度模型。模型尺度為x×y×z=500 m×500 m×250 m，網(wǎng)格間距為10.0 m。震源點(diǎn)坐標(biāo)是(1.0 m，1.0 m，1.0 m)；接收點(diǎn)均勻地分布在模型表面上，其間距為20.0 m；速度分布函數(shù)為v=500+60z。式中：z為深度；v為速度。

對比圖4a、b可見：Chebyshev逼近法所得到的射線路徑在源點(diǎn)附近優(yōu)于B樣條插值法；在其他區(qū)域，兩者所得射線路徑在整體上相差不大。為便于進(jìn)一步對比，取y=25.0 m處的xz剖面進(jìn)行對比，其結(jié)果示于圖5，可以看出，Chebyshev逼近法所給出的走時(shí)具有較小的相對誤差。另外，在一臺CPU為酷睿i5的臺式機(jī)(CPU 頻率為2 450 MHz；內(nèi)存為6 G RAM)上計(jì)算時(shí)，B樣條插值法所用的CPU時(shí)間為5.191 s，而Chebyshev逼近法所用的 CPU時(shí)間為4.205 s，其計(jì)算效率相差約20%。因此，綜合來看，在這個(gè)模型中，Chebyshev逼近法要優(yōu)于B樣條插值法。

a.Chebyshev插值法；b.B樣條插值法。圖4 速度連續(xù)變化模型射線路徑圖Fig.4 Ray path diagram of continuous velocity model

a .Chebyshev插值法；b.B樣條插值法；c.B樣條插值(黑色)與Chebyshev多項(xiàng)式逼近(灰色)走時(shí)相對誤差對比。圖5 速度連續(xù)變化模型射線路徑側(cè)面圖Fig.5 Continuous velocity model ray path side view

3.2 Marmousi模型數(shù)值計(jì)算

為檢驗(yàn)在上文中所提出的方法在復(fù)雜模型中的計(jì)算精度，首先利用三維Marmousi模型進(jìn)行對比研究。圖6給出了三維Marmousi模型中一個(gè)剖面上的對比結(jié)果，不難看出：運(yùn)動學(xué)射線追蹤雖然能準(zhǔn)確地計(jì)算出射線路徑，但是存在射線路徑盲區(qū)，從而不能求出模型內(nèi)的全部走時(shí)；與此相反，在Chebyshev逼近法中，射線路徑是通過走時(shí)來計(jì)算的，而且是程函方程的有限差分解，存在于整個(gè)速度模型之中，因此，在Chebyshev逼近法中沒有射線覆蓋盲區(qū)問題。值得指出的是，雖然利用基于程函方程有限差分法的Chebyshev逼近法能算出整個(gè)模型中的射線路徑，但是考慮到過多的射線條數(shù)會影響對比效果，所以在圖6中只在射線追蹤法的覆蓋區(qū)內(nèi)進(jìn)行了對比。

3.3 多次透射一次反射

黑色射線為運(yùn)動學(xué)射線追蹤結(jié)果；灰色射線為Chebyshev逼近法結(jié)果；白色等值線為用FMM得到的走時(shí)等值線(等時(shí)線)。圖6 Marmousi模型射線路徑對比Fig.6 Marmousi model ray path comparison

圖7a、b分別給出4層水平層狀模型中的透射走時(shí)和在第4層底邊界的反射走時(shí)。圖7c給出的是在第4層底邊界利用Chebyshev逼近法得到的一次反射射線路徑，其中透射部分的射線均垂直穿過在圖7a中所給出的等時(shí)線，反射部分的射線均垂直穿過圖7b中所給出的等時(shí)線。圖7d是利用Chebyshev逼近法求得的所有4個(gè)層面上的入射與反射射線路徑。圖7e給出的是Chebyshev逼近法與樣條插值法之間關(guān)于射線路徑計(jì)算的相對誤差對比，可以看出，Chebyshev逼近法的精度要高于B樣條插值法的精度。

a.4層模型中的透射等時(shí)線；b.第四層底邊界的一次反射等時(shí)線；c.第四層底邊界的一次反射射線路徑；d.4層模型中所有4個(gè)層面上的反射射線路徑；e.根據(jù)d圖計(jì)算的走時(shí)相對誤差。圖7 4層水平層狀模型多次反射透射路徑圖Fig.7 Multiple reflection transmission ray path diagram of 4-level horizontal layered model

a.Marmousi模型分層與剖面路徑圖(據(jù)文獻(xiàn)[8])；b.Marmousi三維路徑圖。圖8 Marmousi模型多次反射路徑圖 Fig.8 Marmousi model multiple reflection ray path diagram

3.4 層間多次反射

為考察Chebyshev逼近法計(jì)算多次反射的效果，首先選取Marmousi模型中的一個(gè)小部分，構(gòu)成局部層狀介質(zhì)(圖8)。由圖8可見：高速區(qū)射線路徑偏向?qū)用娣较?，而低速區(qū)射線偏向?qū)用娴姆ň€方向。這與Snell定律所要求的完全一致，說明基于Chebyshev逼近的射線路徑計(jì)算是正確的。

鑒于在圖8所示模型中，Chebyshev逼近和B樣條插值之間的誤差對比分析與圖7e中的結(jié)果類似，故在此不再進(jìn)行專門的討論。

其次，考察高速度反差條件下Chebyshev逼近法的計(jì)算效果。為此選用二維Sigsbee模型，其計(jì)算結(jié)果展示在圖9中。為了保證能在鹽丘底部形成有效的多次反射，將震源置于遠(yuǎn)離高速體的位置，并以模型的底界面為一次反射面。首先，利用多級次FMM計(jì)算出整個(gè)模型中的多次反射走時(shí)，再從模型頂界面開始進(jìn)行逆向梯度計(jì)算，進(jìn)而得到如圖9所示的射線路徑。為方便起見，只計(jì)算了模型底界面和鹽丘底界面之間的2次反射射線路徑。盡管如此，也可以清晰地看出在高速體存在時(shí)的層間多次反射始終滿足Snell定律，且由高速體向低速區(qū)域呈“掃把”形出射。這進(jìn)一步地說明了Chebyshev逼近法計(jì)算一次和多次反射射線路徑的有效性。

圖9 Sigsbee模型剖面Fig.9 Sigsbee model section

4 結(jié)論

針對利用B樣條插值進(jìn)行走時(shí)梯度計(jì)算所存在的問題(在速度突變時(shí)會出現(xiàn)較大的計(jì)算誤差)，本文提出利用Chebyshev多項(xiàng)式替代B樣條對走時(shí)場進(jìn)行逼近，并在此基礎(chǔ)上結(jié)合多級快速推進(jìn)法(FMM)建立了一種計(jì)算三維多次反射射線的方法。根據(jù)簡單的水平層狀介質(zhì)和復(fù)雜的Marmousi和Sigsbee模型中的數(shù)值計(jì)算結(jié)果可以得到下列結(jié)論：

1)利用Chebyshev多項(xiàng)式逼近走時(shí)可以有效地解決離散數(shù)據(jù)求導(dǎo)過程中出現(xiàn)的不穩(wěn)定問題，以及B樣條插值法在速度突變時(shí)所引起的大誤差問題，并能較好適應(yīng)分區(qū)多級次算法得到的多次反射、透射走時(shí)場。

2)Chebyshev逼近法在簡單和復(fù)雜模型中均能保證所計(jì)算出的射線路徑具有較高的精度，能夠滿足偏移成像的要求。

(

)：

[1] 孫建國. Kirchhoff型偏移理論的研究歷史、研究現(xiàn)狀與發(fā)展趨勢展望:與光學(xué)繞射理論的類比、若干新結(jié)果、新認(rèn)識以及若干有待于解決的問題[J]. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版),2012,42(5):1521-1552.

Sun Jianguo. The History, the State of the Art and the Future Trend of The Research of Kirchhoff-Type Migration Theory: A Comparison with Optical Diffraction Theory, Some New Understanding and Some Problems to be Solved[J]. Journal of Jilin University (Earth Science Edition), 2012, 42(5): 1521-1552.

[2] 孫建國. 高頻漸近散射理論及其在地球物理場數(shù)值模擬與反演成像中的應(yīng)用:研究歷史與研究現(xiàn)狀概述以及若干新進(jìn)展[J]. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版)，2016，46(4)：1231-1259.

Sun Jianguo. High-Frequency Asymptotic Scattering Theories and Their Applications in Numerical Modeling and Imaging of Geophysical Fields: An Overview of the Research History and the State-of-the-Art, and Some New Developments[J]. Journal of Jilin University (Earth Science Edition), 2016, 46(4): 1231-1259.

[3]Vidale J E. Finite-Difference Calculations of Travel-Times[J]. Bulletin of the Seismological Society of America, 1988, 78: 2062-2076.

[4]Vidale J E. Finite-Difference Calculations of Travel-times in Three Dimensions[J]. Geophysics, 1990, 55: 521-526.

[5] Rawlinson N, Sambridge M. Multiple Reflection and Transmission Phases in Complex Layered Media Using a Multistage Fast Marching Method[J]. Geophysics, 2004, 69(5): 1338-1350.

[6] de Kool M, Rawlinson N, Sambridge M. A Practical Grid-Based Method for Tracking Multiple Refraction and Reflection Phases in Three-Dimensional Heterogeneous Media[J]. Geophysical Journal International, 2006, 167(1): 253-270.

[7] 孫章慶，孫建國，王雪秋，等. 三維復(fù)雜山地多級次群推進(jìn)迎風(fēng)混合法多波型走時(shí)計(jì)算[J]. 地球物理學(xué)報(bào)，2017，60(5)：1861-1873.

Sun Zhangqing, Sun Jianguo, Wang Xueqiu, et al. Computation of Multiples Seismic Traveltime in Mountainous Areas with Complex 3D Conditions Using Multistage Group Marching up Wind Hybrid Method[J]. Chinese Journal of Geophysics, 2017, 60(5): 1861-1873.

[8] 唐小平,白超英. 最短路徑算法下三維層狀介質(zhì)中多次波追蹤[J]. 地球物理學(xué)報(bào)，2009，52(10)：2635-2643.

Tang Xiaoping, Bai Chaoying. Multiple Ray Tracing Within 3-D Layered Media with the Shortest Path Method[J]. Chinese Journal of Geophysics, 2009, 52(10): 2635-2643.

[9] Sun Jianguo, Sun Zhangqing, Han Fuxing. A Finite Difference Scheme for Solving the Eikonal Equation Including Surface Topography[J]. Geophysics, 2011, 76(4)：T53-T63.

[10] 石秀林,孫建國,孫輝,等. 基于波前構(gòu)建法的時(shí)間域深度偏移:Delta波包途徑[J]. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版)，2016，46(6)：1847-1854.

Shi Xiulin, Sun Jianguo, Sun Hui, et al. Depth Migration in Time Domain Using Wavefront Construction: Delta Packet Approach[J]. Journal of Jilin University (Earth Science Edition), 2016, 46(6): 1847-1854.

[11] 石秀林，孫建國，孫輝，等．基于delta波包疊加的時(shí)間域深度偏移[J]．地球物理學(xué)報(bào)，2016，59(7)：2641-2649.

Shi Xiulin, Sun Jianguo, Sun Hui, et al. The Time Domain Depth Migration by the Summation of Delta Packets[J]. Chinese Journal of Geophysics, 2016, 59(7), 2641-2649.

[12] 孫建國，何洋. 基于波前構(gòu)建的射線追蹤：一種Java實(shí)現(xiàn)[J].吉林大學(xué)學(xué)報(bào)(地球科學(xué)版)，2007, 37(4)：814-820.

Sun Jianguo, He Yang. Ray-Tracing Based on Wavefront Construction: A Java Implementation[J]. Journal of Jilin University (Earth Science Edition), 2007, 37(4): 814-820.

[13] 韓復(fù)興, 孫建國, 王坤. 速度模型光滑分析[J]. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版)，2011, 41(5): 1610-1616.

Han Fuxing, Sun Jianguo, Wang Kun. Analysis of Velocity Model Smoothing[J]. Journal of Jilin University (Earth Science Edition), 2011, 41(5): 1610-1616.

[14] 韓復(fù)興，孫建國，孫章慶. 波前構(gòu)建法研究現(xiàn)狀[J]. 地球物理學(xué)進(jìn)展，2011，26(3): 1045-1051.

Han Fuxing, Sun Jianguo, Sun Zhangqing. Research Status of the Wavefront Construction Method[J]. Progress in Geophysics, 2011, 26(3): 1045-1051.

[15] 張東，張婷婷，喬友鋒，等.三維旅行時(shí)場B樣條插值射線追蹤方法[J]. 石油地球物理勘探，2013，48(4)：559-566.

Zhang Dong, Zhang Tingting, Qiao Youfeng, et al. A Ray Tracing Method Based on B-Spline Traveltime Interpolation[J]. Oil and Gas Prospecting, 2013, 48(4)：559-566.

[16] 張婷婷，張東，邱達(dá). 三維層狀介質(zhì)中基于走時(shí)梯度的多次波射線追蹤[J]. 石油地球物理勘探，2014，49(6)：1097-1105.

Zhang Tingting, Zhang Dong, Qiu Da. Multiple Ray Tracing in 3D Layered Media Using the Traveltime Gradient Method[J]. Oil and Gas Prospecting, 2014，49(6)：1097-1105.

[17] 張婷婷，邱達(dá)，張東. 一種改進(jìn)的三維旅行時(shí)梯度射線追蹤方法[J]. 石油地球物理勘探，2016，51(5)：916-923.

Zhang Tingting, Qiu Da, Zhang Dong. A Modified 3D Traveltime Gradient Ray Tracing Method[J]. Oil and Gas Prospecting, 2016，51(5)：916-923.

[18] Press W H, Teukolsky S A, Vettering W T, et al. Numerical Recipies in FORTRAN[M]. Cambridge: Cambridge University Press, 1992.

[19] 蔣迅，王淑紅. 切比雪夫和切比雪夫多項(xiàng)式的故事[J]. 科學(xué), 2016, 68(4): 54-58.

Jiang Xun, Wang Shuhong. The Story of Chebyshev and Chebyshev Polynomials[J]. Science, 2016, 68(4): 54-58.