孫建國(guó)，苗 賀

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

0 引言

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

在歷史上，利用程函方程數(shù)值解計(jì)算地震波走時(shí)的基本思想早已有之，但是由于沒(méi)有相應(yīng)的需求而一直不為人所關(guān)注。1988年，Vidale首先將前人的基本思想付諸實(shí)踐[3-4]，提出了基于程函方程有限差分(FD)數(shù)值解的地震波走時(shí)計(jì)算方法。與上文提到的初值方法相比，F(xiàn)D法具有計(jì)算速度快和沒(méi)有計(jì)算盲區(qū)等特點(diǎn)，從而可以作為一種理想的地震波走時(shí)計(jì)算方法替代傳統(tǒng)的射線追蹤法。然而，在無(wú)人為干預(yù)時(shí)，利用有限差分等基于網(wǎng)格的程函方程數(shù)值分析方法只能沿著波的前進(jìn)方向計(jì)算走時(shí)，而不能沿著波的回折方向計(jì)算走時(shí)。換句話說(shuō)，在無(wú)人為干預(yù)時(shí)，利用任何一種基于程函方程數(shù)值分析的方法均只能計(jì)算沿著波的前進(jìn)方向的透射(波)走時(shí)，而不能同時(shí)計(jì)算由該透射波所激發(fā)的反射(波)走時(shí)，更不能計(jì)算由該透射波所激發(fā)的、其傳播方向經(jīng)過(guò)多次改變的多次透射(波)和多次反射(波)走時(shí)。針對(duì)這個(gè)問(wèn)題，Rawlinson等[5]在2004年提出了一種在人為干預(yù)下利用程函方程數(shù)值分析方法計(jì)算多次透射和多次反射走時(shí)的分區(qū)多級(jí)方案，在這種方案中，走時(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)過(guò)程被反復(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ì)算方法與群推進(jìn)法組合在一起，提出了一種針對(duì)復(fù)雜山地多波型走時(shí)計(jì)算的多級(jí)次群推進(jìn)迎風(fēng)混合法。類(lèi)似地，唐小平等[8]將多級(jí)計(jì)算方法與最短路徑算法組合在一起，提出了一種基于最短路徑法的、針對(duì)三維層狀介質(zhì)中多次波的走時(shí)與射線路徑計(jì)算方法。在理論上，分區(qū)多級(jí)計(jì)算的基本思想可以用于任何一種基于程函方程數(shù)值分析的地震波走時(shí)計(jì)算方法，因此存在有無(wú)限多的組合方式。

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

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

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

(1)

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

(2)

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

圖2給出了上述分區(qū)多級(jí)計(jì)算步驟的進(jìn)一步說(shuō)明。①首先由震源開(kāi)始，計(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ì)算過(guò)程，即可得到任意級(jí)次的多次反射和多次透射。

圖1 反射(透射)區(qū)內(nèi)的接收點(diǎn)(黑色圓點(diǎn))、活動(dòng)點(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ū)多級(jí)走時(shí)計(jì)算過(guò)程示意圖Fig.2 Partition multi-level travel time calculation process

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

Chebyshev逼近是基于Chebshev多項(xiàng)式的一種逼近方法，是為解決零偏差問(wèn)題而提出的一種基于正交多項(xiàng)式的最佳逼近方法[18]。在歷史上，Chebyshev逼近和Chebshev多項(xiàng)式由俄羅斯數(shù)學(xué)家P L Chebyshev在1854年提出，其最初目的是解決連桿設(shè)計(jì)問(wèn)題[19]。根據(jù)其數(shù)學(xué)表現(xiàn)形式可將Chebyshev提出的多項(xiàng)式分為兩種，即第一類(lèi)和第二類(lèi)Chebyshev多項(xiàng)式。其中，第一類(lèi)Chebyshev多項(xiàng)式起源于多倍角余弦和正弦函數(shù)展開(kāi)，用Tn(x)表示。具體地講，Tn(x)=cos[ncos-1(x)]。第二類(lèi)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ù)展開(kāi)為Chebyshev級(jí)數(shù)等等。因此，為在最小平方意義下逼近來(lái)自于快速推進(jìn)法的走時(shí)場(chǎng)，在下文中選取Chebyshev多項(xiàng)式為逼近函數(shù)。如圖3所示，在諸多正交多項(xiàng)式中，Chebyshev多項(xiàng)式的方差最小，即逼近誤差亦為最小。

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

(3)

展開(kāi)后得到

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)及方差對(duì)比圖(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ù)對(duì)應(yīng)的正交項(xiàng)；M,N,L分別為x,y,z方向的網(wǎng)格節(jié)點(diǎn)數(shù)。

一旦完成了公式(4)所描述的逼近過(guò)程，即可根據(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為速度。

對(duì)比圖4a、b可見(jiàn)：Chebyshev逼近法所得到的射線路徑在源點(diǎn)附近優(yōu)于B樣條插值法；在其他區(qū)域，兩者所得射線路徑在整體上相差不大。為便于進(jìn)一步對(duì)比，取y=25.0 m處的xz剖面進(jìn)行對(duì)比，其結(jié)果示于圖5，可以看出，Chebyshev逼近法所給出的走時(shí)具有較小的相對(duì)誤差。另外，在一臺(tái)CPU為酷睿i5的臺(tái)式機(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%。因此，綜合來(lái)看，在這個(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í)相對(duì)誤差對(duì)比。圖5 速度連續(xù)變化模型射線路徑側(cè)面圖Fig.5 Continuous velocity model ray path side view

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

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

3.3 多次透射一次反射

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

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

a.4層模型中的透射等時(shí)線；b.第四層底邊界的一次反射等時(shí)線；c.第四層底邊界的一次反射射線路徑；d.4層模型中所有4個(gè)層面上的反射射線路徑；e.根據(jù)d圖計(jì)算的走時(shí)相對(duì)誤差。圖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可見(jiàn)：高速區(qū)射線路徑偏向?qū)用娣较?，而低速區(qū)射線偏向?qū)用娴姆ň€方向。這與Snell定律所要求的完全一致，說(shuō)明基于Chebyshev逼近的射線路徑計(jì)算是正確的。

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

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

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

4 結(jié)論

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

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

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

(

)：

[1] 孫建國(guó). Kirchhoff型偏移理論的研究歷史、研究現(xiàn)狀與發(fā)展趨勢(shì)展望:與光學(xué)繞射理論的類(lèi)比、若干新結(jié)果、新認(rèn)識(shí)以及若干有待于解決的問(wè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] 孫建國(guó). 高頻漸近散射理論及其在地球物理場(chǎng)數(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] 孫章慶，孫建國(guó)，王雪秋，等. 三維復(fù)雜山地多級(jí)次群推進(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] 石秀林,孫建國(guó),孫輝,等. 基于波前構(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] 石秀林，孫建國(guó)，孫輝，等．基于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] 孫建國(guó)，何洋. 基于波前構(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ù)興, 孫建國(guó), 王坤. 速度模型光滑分析[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ù)興，孫建國(guó)，孫章慶. 波前構(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í)場(chǎng)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.