旅行時(shí)線性插值射線追蹤提高計(jì)算精度和效率的改進(jìn)方法

葉 佩，李慶春

1.中國(guó)國(guó)土資源航空物探遙感中心，北京 100083

2.長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院，西安 710054

0 引言

射線理論及射線方法是研究地震波傳播理論的重要途徑之一，射線理論也經(jīng)常被用于研究不均勻介質(zhì)和復(fù)雜構(gòu)造中的地震波傳播問(wèn)題。利用旅行時(shí)信息進(jìn)行層析成像最常用的正演方法就是射線追蹤，該方法還被廣泛用于偏移、地震波反演及其他正演模擬中，計(jì)算地震波的旅行時(shí)和波前［1－11］在實(shí)際應(yīng)用中得到了不斷的發(fā)展與完善。旅行時(shí)線性插值（linear travel－time interpolation，LTI）射線追蹤法是Asawaka［4］提出來(lái)的一種基于線性假設(shè)的網(wǎng)格單元擴(kuò)展方法，該方法已廣泛應(yīng)用于層析成像，具有計(jì)算速度快、精度高、易于實(shí)現(xiàn)等優(yōu)點(diǎn)。近年來(lái)，國(guó)內(nèi)很多學(xué)者對(duì)該方法做了大量的研究。如：張建中等［6］推導(dǎo)出了旅行時(shí)非線性插值及變速介質(zhì)中旅行時(shí)插值計(jì)算公式；為了改善該方法在近場(chǎng)旅行時(shí)的精度，降低累計(jì)誤差，聶建新等［12］提出了旅行時(shí)二次／線性聯(lián)合插值（簡(jiǎn)稱 QLTI）；彭直興等［13］在震源的近場(chǎng)采用非線性插值計(jì)算旅行時(shí)，在遠(yuǎn)場(chǎng)仍用線性插值的計(jì)算方法，改善了傳統(tǒng)的LTI法因線性插值引起的誤差；為了改善旅行時(shí)的計(jì)算精度和計(jì)算效率，孫章慶等［14－15］利用窄帶技術(shù)來(lái)模擬波前，并在局部采用線性插值計(jì)算走時(shí)，隨后，又將該方法擴(kuò)展到復(fù)雜地表走時(shí)計(jì)算。由于傳統(tǒng)的LTI法在進(jìn)行旅行時(shí)計(jì)算時(shí)射線方向考慮不全，計(jì)算的節(jié)點(diǎn)旅行時(shí)不一定最小，導(dǎo)致追蹤的射線路徑不一定滿足最小旅行時(shí)。張東等［16－17］提出在向前處理中采用多方向的循環(huán)計(jì)算，使所計(jì)算的節(jié)點(diǎn)旅行時(shí)達(dá)到最??；同年，又將該方法擴(kuò)展到三維。梅勝全等［18］用次震源波前面掃描算法代替旅行時(shí)插值并通過(guò)判定條件進(jìn)行快速插值的方法來(lái)提高計(jì)算速度。

Asawaka在文獻(xiàn)［4］中指出，LTI算法向前處理過(guò)程中計(jì)算節(jié)點(diǎn)旅行時(shí)，應(yīng)先按列掃描再進(jìn)行按行掃描，但其并沒(méi)有對(duì)按行掃描進(jìn)行進(jìn)一步的分析，也沒(méi)有給出計(jì)算實(shí)例說(shuō)明若不這樣做會(huì)存在什么樣的問(wèn)題。國(guó)內(nèi)學(xué)者在對(duì)該算法進(jìn)行研究時(shí)，也沒(méi)有提到在按列掃描之后還要進(jìn)行按行掃描，相反讓人誤以為L(zhǎng)TI算法向前處理只需按列掃描即可。但實(shí)際上，在向前處理過(guò)程中僅按列或按行掃描均沒(méi)有考慮到逆向傳播的射線，在處理速度變化劇烈的介質(zhì)時(shí)，會(huì)導(dǎo)致所追蹤出來(lái)的射線并不滿足最小旅行時(shí)。因此，筆者通過(guò)實(shí)例分析，說(shuō)明LTI算法在向前處理過(guò)程中應(yīng)采用全方位循環(huán)掃描的方法，并給出了其具體實(shí)現(xiàn)步驟；在保證計(jì)算精度的同時(shí)，在邊界加入次生節(jié)點(diǎn)的方法可以大大提高該方法的計(jì)算效率。

1 LTI法的基本原理

LTI方法存在一個(gè)線性假設(shè)前提，如圖1所示：設(shè)A（xA，zA）、B（xB，zB）兩點(diǎn)的坐標(biāo)已知，A、B 兩點(diǎn)的旅行時(shí)tA、tB也是已知的，則線段A、B之間任意一點(diǎn)D（xD，zD）的旅行時(shí)tD可由A、B兩點(diǎn)的旅行時(shí)tA、tB線性插值得到，即

其中：tA為旅行時(shí)較小的點(diǎn)，即tB≥tA為線段AB長(zhǎng)度；r為A到D的距離，滿足

圖1 LTI法旅行時(shí)線性插值幾何關(guān)系Fig.1 Ttraveltime linear interpolation geometric relationship of LTI method

若射線從A，B之間的D 點(diǎn)穿過(guò)到達(dá)C點(diǎn)，則C點(diǎn)的旅行時(shí)tC可表示為

其中：v為介質(zhì)的波傳播速度；φ為線段AB 與AC之間的夾角。

式（2）中C點(diǎn)的旅行時(shí)是關(guān)于r的函數(shù)。根據(jù)最短旅行時(shí)原理，射線從D點(diǎn)傳到C點(diǎn)應(yīng)滿足

聯(lián)立公式（2）、（3）可以解出

將式（4）代入式（2）可得C點(diǎn)最小旅行時(shí)為

由式（4）可得D點(diǎn)的坐標(biāo)為

該方法可以分為向前處理和向后處理兩步：向前處理是從炮點(diǎn)開(kāi)始計(jì)算炮點(diǎn)到各網(wǎng)格邊界節(jié)點(diǎn)的最短旅行時(shí)，可由式（5）求得；向后處理是從接收點(diǎn)開(kāi)始，確定所有炮點(diǎn)-接收點(diǎn)的射線路徑與各單元網(wǎng)格邊界的交點(diǎn)，可由式（6）求得［4，19］。

2 LTI算法存在的問(wèn)題及其改進(jìn)

2.1 復(fù)雜介質(zhì)中的逆向傳播

LTI算法在向前處理過(guò)程中若只按列或按行掃描來(lái)計(jì)算節(jié)點(diǎn)旅行時(shí)沒(méi)有考慮到逆向傳播的射線。以按列掃描為例，如圖2a所示，S點(diǎn)表示炮點(diǎn)所在位置，R點(diǎn)是其中的一個(gè)計(jì)算節(jié)點(diǎn)，位于炮點(diǎn)右方。LTI算法在向前處理過(guò)程中對(duì)R點(diǎn)最小旅行時(shí)的計(jì)算只考慮了炮點(diǎn)S，計(jì)算點(diǎn)R所在左半平面的所有射線；而在復(fù)雜的速度結(jié)構(gòu)中，射線極有可能先繞到R點(diǎn)的右半平面再傳到R點(diǎn)，如圖2b所示。同理，當(dāng)計(jì)算點(diǎn)位于炮點(diǎn)左邊，也存在同樣的問(wèn)題。這樣會(huì)導(dǎo)致在橫向速度變化劇烈的介質(zhì)中（即出現(xiàn)逆向傳播），計(jì)算出來(lái)的節(jié)點(diǎn)旅行時(shí)不是最小，向后處理中追蹤出來(lái)的射線路徑也不滿足最短旅行時(shí)原理，降低了該方法對(duì)復(fù)雜模型的適應(yīng)能力。

圖2 逆向傳播的射線示意圖Fig.2 Schematic illustrations of counterpropagating rays

圖3為Asakawa和Kawanaka的一個(gè)經(jīng)典模型［4］，是一個(gè)三層速度模型，兩低速層中夾有一高速層。模型剖分網(wǎng)格單元大小為0.1m×0.1m。模型參數(shù)及觀測(cè)系統(tǒng)如圖3a所示，圖3b是只按列掃描LTI算法追蹤的初至射線路徑，與文獻(xiàn)［4］的結(jié)果吻合。

圖3 Askawa經(jīng)典模型試驗(yàn)Fig.3 Classic model test of Askawa

為說(shuō)明射線逆向傳播的問(wèn)題，將圖3a的模型及觀測(cè)系統(tǒng)整體順時(shí)針旋轉(zhuǎn)90°，如圖4a所示，重新追蹤出的初至波射線路徑見(jiàn)圖4b。對(duì)比圖3b，4b中的射線路徑可以看出，旋轉(zhuǎn)后追蹤的1、2兩條射線路徑有誤。由表1也可以看出，旋轉(zhuǎn)后僅用按列掃描的LTI算法計(jì)算出來(lái)的1、2兩條射線的旅行時(shí)不是最小。由圖3b中的射線路徑可知，正確情況下圖4b中標(biāo)示1、2的2條射線應(yīng)該是由左側(cè)的高速層折射回來(lái)。究其原因是只按列掃描的LTI法在向前處理過(guò)程中，計(jì)算震源左側(cè)單元節(jié)點(diǎn)最小旅行時(shí)只考慮了來(lái)自右半平面的射線，并沒(méi)有考慮來(lái)自左半平面的射線。同理，只按行掃描也存在射線考慮不全的問(wèn)題。因此，筆者認(rèn)為只按列或按行掃描的LTI法在向前處理過(guò)程并沒(méi)有考慮來(lái)自各個(gè)方向的射線，只有采用全方位循環(huán)的方法才能確保計(jì)算得到的各節(jié)點(diǎn)的旅行時(shí)滿足最小。

圖4 LTI法向前處理按列掃描及其缺陷Fig.4 Forward processing by column scanning and its defects of LTI method

表1 Asakawa模型旋轉(zhuǎn)前后旅行時(shí)計(jì)算比較Table1 Traveltime comparison of Asakawa’s model calculated before and after rotation

在具體實(shí)施中，先按LTI算法按列掃描計(jì)算各節(jié)點(diǎn)的旅行時(shí)，再按行掃描的方法來(lái)重新計(jì)算各網(wǎng)格節(jié)點(diǎn)旅行時(shí)；將按行掃描計(jì)算的節(jié)點(diǎn)旅行時(shí)與按列掃描計(jì)算的節(jié)點(diǎn)旅行時(shí)相比較，取其中旅行時(shí)最小的值作為網(wǎng)格節(jié)點(diǎn)的旅行時(shí)。這樣就可以確保在向前處理過(guò)程中考慮來(lái)自各個(gè)方向的射線，即所計(jì)算出來(lái)的節(jié)點(diǎn)旅行時(shí)是最小的。

全方位循環(huán)法向前處理實(shí)現(xiàn)步驟如下。

1）先用LTI算法按列掃描計(jì)算出整個(gè)區(qū)域各網(wǎng)格節(jié)點(diǎn)的旅行時(shí)［4］。

2）除震源所在網(wǎng)格外，重新計(jì)算震源所在行各節(jié)點(diǎn)的旅行時(shí)（按行掃描），如圖5a所示，并與原先按列掃描計(jì)算的結(jié)果對(duì)比，取其中最小值作為節(jié)點(diǎn)的旅行時(shí)。

3）由下到上逐個(gè)計(jì)算震源所在行上面所有節(jié)點(diǎn)的旅行時(shí)：

①計(jì)算各單元水平邊界及垂直邊界上各節(jié)點(diǎn)的旅行時(shí)，這時(shí)只考慮來(lái)自下邊界的射線（圖5b），并與原先按列計(jì)算的比較，取其中旅行時(shí)最小的；

②從左往右，利用單元左邊界上的已知點(diǎn)計(jì)算右邊界和上邊界節(jié)點(diǎn)的旅行時(shí)（圖5c），并與節(jié)點(diǎn)已有值比較，取其中旅行時(shí)最小的；

③從右往左，利用單元右邊界上的已知點(diǎn)計(jì)算左邊界和上邊界節(jié)點(diǎn)的旅行時(shí)（圖5d），并與節(jié)點(diǎn)已有值比較，取其中旅行時(shí)最小的。

通過(guò)①、②、③步計(jì)算，就可求得震源所在行上面所有節(jié)點(diǎn)的最小旅行時(shí)。

4）同理，從上到下可以計(jì)算出震源所在行下面各節(jié)點(diǎn)的旅行時(shí)，并與按列計(jì)算的旅行時(shí)比較取其中最小的。最后可得到整個(gè)模型區(qū)域內(nèi)所有節(jié)點(diǎn)的最小旅行時(shí)。

圖6為L(zhǎng)TI算法全方位循環(huán)對(duì)圖4a追蹤的結(jié)果。與圖3b相比，追蹤的射線其路徑和理論值是一樣的，尤其是標(biāo)示為1、2的兩條射線，其旅行時(shí)與理論值一致（表1）。由試算結(jié)果可知，LTI算法全方位循環(huán)法在向前處理過(guò)程中可以確保計(jì)算出來(lái)的節(jié)點(diǎn)旅行時(shí)最小，從而能夠避免僅按列或按行掃描未考慮逆向傳播的射線導(dǎo)致在橫向變化劇烈的模型中所追蹤出來(lái)的初至射線路徑不滿足最小旅行時(shí)原理的問(wèn)題。

2.2 LTI算法的計(jì)算效率

LTI法是基于線性假設(shè)的，模型網(wǎng)格剖分的精細(xì)程度直接影響射線旅行時(shí)的計(jì)算精度：在一定范圍內(nèi)，網(wǎng)格剖分越精細(xì)，射線旅行時(shí)精度越高；但網(wǎng)格剖分得越精細(xì)，在計(jì)算時(shí)，所耗費(fèi)的時(shí)間也就越長(zhǎng)。因此在保證計(jì)算精度同時(shí)，如何提高計(jì)算效率也是目前需解決的問(wèn)題。

為了保證LTI法的計(jì)算精度，同時(shí)提高該方法的計(jì)算效率，筆者提出在網(wǎng)格單元邊界加入次生節(jié)點(diǎn)的技巧。

圖5 向前處理按行掃描計(jì)算節(jié)點(diǎn)最小旅行時(shí)的實(shí)現(xiàn)過(guò)程Fig.5 Implementation process of minimum travel－time of grid node calculated by line scanning in the forward processing

圖6 LTI算法全方位循環(huán)法模型射線追蹤路徑結(jié)果圖Fig.6 Ray traces calculated by all－round cycling of LTI method

圖7a是改進(jìn)前LTI算法網(wǎng)格節(jié)點(diǎn)示意圖，圖7b為改進(jìn)后在網(wǎng)格邊界插入次生節(jié)點(diǎn)的示意圖?？梢钥闯?，在剖分精度相同的情況下，傳統(tǒng)LTI計(jì)算點(diǎn)數(shù)多。改進(jìn)后只在主節(jié)點(diǎn)上給定速度值，邊界上插入的次生節(jié)點(diǎn)速度可由單元邊界上相鄰兩主節(jié)點(diǎn)的速度線性插值得到；而原來(lái)的做法（圖7a）則需在每個(gè)節(jié)點(diǎn)都賦速度值。因此在將該方法用于層析成像時(shí)，插入次生節(jié)點(diǎn)的LTI法在不影響計(jì)算精度的前提下所需反演的參數(shù)更少，更有利于反演。

為了驗(yàn)證加入次生節(jié)點(diǎn)后LTI法的計(jì)算精度，建立一均勻速度模型。模型參數(shù)及觀測(cè)系統(tǒng)如圖8a所示。圖8b顯示的是網(wǎng)格剖分個(gè)數(shù)為10×10，單元邊界次生節(jié)點(diǎn)數(shù)為99×99，即節(jié)點(diǎn)間距為1m情況下的初至波射線路徑；表2列出了傳統(tǒng)的LTI法與網(wǎng)格邊界插入次生節(jié)點(diǎn)的LTI法在節(jié)點(diǎn)間距為5m和1m兩種情況下各道所接收到的初至波旅行時(shí)；圖9給出了不同節(jié)點(diǎn)間距下插入次生節(jié)點(diǎn)前后LTI法的計(jì)算耗時(shí)。

圖7 網(wǎng)格邊界插入次生節(jié)點(diǎn)前（a）后（b）布置示意圖Fig.7 Schematic illustrations of before（a）and after（b）secondary node inserted into mesh boundary

由表2可以看出：在節(jié)點(diǎn)間距為5m的情況下，傳統(tǒng)的LTI法絕對(duì)誤差最大值為0.391ms，在單元邊界插入次生節(jié)點(diǎn)后其絕對(duì)誤差最大值降為0.006 ms；在節(jié)點(diǎn)間距為1m的情況下，傳統(tǒng)的LTI法絕對(duì)誤差最大值為0.084ms，單元邊界插入次生節(jié)點(diǎn)后絕對(duì)誤差最大值為0ms。從圖9中可以看到在節(jié)點(diǎn)間距相同的情況下，傳統(tǒng)的LTI法比插入次生節(jié)點(diǎn)的LTI法耗時(shí)少，而且隨著節(jié)點(diǎn)間距逐漸減小，耗時(shí)的下降更為顯著。從試算結(jié)果可以看到：在節(jié)點(diǎn)間距相同的情況下，網(wǎng)格邊界插入次生節(jié)點(diǎn)的LTI法較傳統(tǒng)的LTI法計(jì)算精度至少可以提高一個(gè)數(shù)量級(jí)。同時(shí)，計(jì)算速度也更快，隨著節(jié)點(diǎn)間距剖分的越精細(xì)，計(jì)算耗時(shí)下降也越明顯，計(jì)算速度較傳統(tǒng)的方法可提高n倍。

圖8 網(wǎng)格邊界插入次生節(jié)點(diǎn)模型試驗(yàn)Fig.8 Model test of secondary node inserted into grid boundary

3 復(fù)雜模型測(cè)試

為了綜合檢驗(yàn)改進(jìn)后LTI法對(duì)復(fù)雜模型的適應(yīng)能力，筆者引用了經(jīng)典的Marmousi模型對(duì)改進(jìn)前后的LTI法進(jìn)行試算及誤差分析。

圖9 不同節(jié)點(diǎn)間距下插入次生節(jié)點(diǎn)前后LTI法CPU耗時(shí)Fig.9 CPU time of LTI method before and after secondary nodes inserted into grid boundary with different node spacing

圖10為Marmousi模型的試算結(jié)果。模型參數(shù)為：網(wǎng)格數(shù)384×122，網(wǎng)格大小10m×10m，速度位于網(wǎng)格節(jié)點(diǎn)。觀測(cè)系統(tǒng)設(shè)置參數(shù)：激發(fā)點(diǎn)位于模型底部（2 800，1 220），在模型頂部均勻放置77個(gè)接收點(diǎn)，接收點(diǎn)分別位于：（0，0），（50，0），（100，0），…，（3 800，0）。圖10是改進(jìn)后的LTI算法追蹤出來(lái)的初至射線路徑，次生節(jié)點(diǎn)數(shù)為4×4。由射線路徑可以看出，在模型底部，射線沿著高速層走，符合最小旅行時(shí)射線路徑的規(guī)律。圖11是傳統(tǒng)的LTI法計(jì)算的各道初至旅行時(shí)減去改進(jìn)后的LTI法計(jì)算的各道初至旅行時(shí)得到的旅行時(shí)殘差曲線。可以看到，所計(jì)算出來(lái)的旅行時(shí)殘差均處于殘差為0的參考線上方，也就是說(shuō)改進(jìn)后的LTI算法算出來(lái)的初至旅行時(shí)小于傳統(tǒng)的LTI法所計(jì)算出來(lái)的初至旅行時(shí)?？梢则?yàn)證改進(jìn)后的LTI算法較傳統(tǒng)的LTI法計(jì)算精度更高，且改進(jìn)后的LTI法對(duì)復(fù)雜模型具有更強(qiáng)的適應(yīng)能力及穩(wěn)定性。

圖10 Marmousi模型改進(jìn)LTI法初至射線路徑示意圖（★.激發(fā)點(diǎn)）Fig.10 Schematic illustrations of Marmousi model firstbreak races calculated by improved LTI method

表2 網(wǎng)格邊界插入次生節(jié)點(diǎn)LTI法與傳統(tǒng)LTI法初至波旅行時(shí)計(jì)算結(jié)果對(duì)比Table2 First－break traveltime comparison of traditional LTI and second node inserted into grid boundary LTI ms

圖11 傳統(tǒng)LTI法與改進(jìn)后LTI法追蹤的旅行時(shí)殘差示意圖Fig.11 Schematic illustrations of traveltime residuals between traditional and improved LTI

4 結(jié)論

1）筆者用實(shí)例詳細(xì)分析了LTI算法僅按列或按行掃描沒(méi)有考慮到逆向傳播的射線，因此建議在向前處理時(shí)采用全方位循環(huán)的方法計(jì)算網(wǎng)格節(jié)點(diǎn)走時(shí)，這樣可以確保追蹤出來(lái)的射線路徑滿足最短旅行時(shí)。

2）在保證計(jì)算精度的前提下，為能有效地提高計(jì)算效率，筆者采用在網(wǎng)格邊界加入次生節(jié)點(diǎn)的措施。試算結(jié)果表明，較傳統(tǒng)的LTI法，本文所采用的LTI算法其適應(yīng)能力更強(qiáng)、計(jì)算精度及計(jì)算效率均大幅度提高。

3）本文采用2D模型，由于其高精度、高效率的特點(diǎn)，稍加擴(kuò)展后即可較好地適用三維射線追蹤及三維層析成像。

（References）：

［1］Vidale J E.Finite Difference Calculation of Traveltimes［J］.Bulletin of the Seismological Society of America，1988，78（6）：2062－2076.

［2］Vidale J E.Finite Difference Calculation of Traveltime in Three Dimensions［J］.Geophysics，1990，55（5）：521－526.

［3］Moser T J.Shortest Path Calculation of Seismic Rays［J］.Geophysics，1991，56：59－67.

［4］Asawaka E，Kawanaka T.Seismic Ray Tracing Using Linear Traveltime Interpolation ［J］. Geophysical Prospecting，1993，41：99－111.

［5］Sethian J A.A Fast Marching Level Set Method for Monotonically Advancing Fronts［J］.Proe Nat Acad Sci，1996，93（4）：1591－1595.

［6］張建中，陳世軍.復(fù)雜介質(zhì)地震初至波數(shù)值模型［J］.計(jì)算物理，2003，20（5）：429－433.Zhang Jianzhong，Chen Shijun.Numerical Modeling of Seismic First－Break in Complex Media［J］.Chinese Journal of Computational Physics，2003，20（5）：429－433.

［7］唐小平，白超英.最短路徑算法下三維層狀介質(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 J Geophys，2009，52（10）：2635－2643.

［8］孫章慶，孫建國(guó)，韓復(fù)興.復(fù)雜地表?xiàng)l件下快速推進(jìn)法地震波走時(shí)計(jì)算［J］.計(jì)算物理，2010，27（2）：281－286.Sun Zhangqing，Sun Jianguo，Han Fuxing.Traveltime Computation Using Fast Marching Method Under Complex Topographical Conditions ［J］. Chinese Journal of Computational Physics，2010，27（2）：281－286.

［9］孫章慶，孫建國(guó)，韓復(fù)興.針對(duì)復(fù)雜地形的三種地震波走時(shí)算法及對(duì)比［J］.地球物理學(xué)報(bào)，2012，55（2）：560－568.Sun Zhangqing，Sun Jianguo，Han Fuxing.The Comparison of Three Schemes for Computing Seismic Wave Traveltimes in Complex Topographical Conditions［J］.Chinese J Geophys，2012，55（2）：560－568.

［10］韓復(fù)興，孫建國(guó)，楊昊.基于二維三次卷積插值算法的波前構(gòu)建射線追蹤［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，2008，38（2）：336－340.Han Fuxing，Sun Jianguo，Yang Hao.Ray－Tracing of Wavefront Construction by Bicubic Convolution Interpolation［J］.Journal of Jilin University：Earth Science Edition，2008，38（2）：336－340.

［11］趙昭，趙志新，徐紀(jì)人，等.高精度人工地震波層析成像技術(shù)勘探地下構(gòu)造［J］.吉林大學(xué)學(xué)報(bào)：地球科學(xué)版，2011，41（1）：362－367.Zhao Zhao，Zhao Zhixin，Xu Jiren，et al.Surveying Structures with High－Resolution Explosion Seismic Tomography［J］.Journal of Jilin University：Earth Science Edition，2011，41（1）：362－367.

［12］聶建新，楊慧珠.地震波旅行時(shí)二次／線性聯(lián)合插值法［J］.清華大學(xué)學(xué)報(bào)：自然科學(xué)版，2003，43（11）：1495－1498.Nie Jianxin，Yang Huizhu.Quadratic／Linear Travel－Time Interpolation of Seismic Ray Tracing［J］.J Tsinghua Univ：Sci ＆ Tech，2003，43（11）：1495－1498.

［13］彭直興，張芬，周熙襄，等.地震波旅行時(shí)非線性／線性聯(lián)合插值法［J］.成都理工大學(xué)學(xué)報(bào)：自然科學(xué)版，2005，32（3）：322－324.Peng Zhixing，Zhang Fen，Zhou Xixiang，et al.Nonlinear／Linear Interpolation of Seismic Traveltime［J］.Journal of Chengdu University of Techonlogy：Science ＆ Technology Edition，2005，32（3）：322－324.

［14］孫章慶，孫建國(guó)，韓復(fù)興.復(fù)雜地表?xiàng)l件下基于線性插值和窄帶技術(shù)的地震波走時(shí)計(jì)算［J］.地球物理學(xué)報(bào)，2009，52（11）：2846－2853.Sun Zhangqing， Sun Jianguo， Han Fuxing.Traveltimes Computation Using Linear Interpolation and Narrow Band Technique Under Complex Topographical Conditions［J］.Chinese J Geophys，2009，52（11）：2846－2853.

［15］孫章慶，孫建國(guó)，韓復(fù)興.基于線性插值和窄帶技術(shù)的走時(shí)計(jì)算方法［J］.石油地球物理勘探，2009，44（4）：436－441.Sun Zhangqing， Sun Jianguo， Han Fuxing.Traveltime Computation Based on Linear Interpolation and Narrow Band Technique［J］.Oil Geophysical Prospecting，2009，44（4）：436－441.

［16］張東，謝寶蓮，楊艷，等.一種改進(jìn)的線性旅行時(shí)插值射線追蹤算法［J］.地球物理學(xué)報(bào)，2009，52（1）：200－205.Zhang Dong，Xie Baolian，Yang Yan，et al.A Ray Tracing Method Based on Improved Linear Travel－Time Interpolation［J］.Chinese J Geophys，2009，52（1）：200－205.

［17］張東，傅相如，楊艷，等.基于LTI和網(wǎng)格界面剖分的三維地震射線追蹤算法［J］.地球物理學(xué)報(bào)，2009，52（9）：2370－2376.Zhang Dong，F(xiàn)u Xiangru，Yang Yan，et al.3－D Seismic Ray Tracing Algorithm Based on LTI and Partition of Grid Interface［J］.Chinese J Geophys，2009，52（9）：2370－2376.

［18］梅勝全，鄧飛，鐘本善，等.基于改進(jìn)的雙線性旅行時(shí)插值的三維射線追蹤［J］.物探化探計(jì)算技術(shù)，2010，32（2）：152－156.Mei Shengquan，Deng Fei，Zhong Benshan，et al.3－D Ray Tracing Based on Improved Bilinear Travel－Time Interpolation ［J］. Geophysical ＆ Geochemical Computing Technique，2010，32（2）：152－156.

［19］葉佩，李慶春.旅行時(shí)線性插值射線追蹤及提高計(jì)算精度和效率的改進(jìn)措施［C］／／中國(guó)地球物理學(xué)會(huì)年會(huì).長(zhǎng)沙：中國(guó)科學(xué)技術(shù)大學(xué)出版社，2011：757.Ye Pei，Li Qingchun.Improvements of Linear Traveltime Interpolation Ray Tracing for the Accuracy and Effiency ［C］／／Annual Meeting of Chinese Geophysical Society. Changsha： China Science and Technology University Press，2011：757.