葉 佩，李慶春

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

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

0 引言

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

Asawaka在文獻(xiàn)［4］中指出，LTI算法向前處理過程中計(jì)算節(jié)點(diǎn)旅行時(shí)，應(yīng)先按列掃描再進(jìn)行按行掃描，但其并沒有對按行掃描進(jìn)行進(jìn)一步的分析，也沒有給出計(jì)算實(shí)例說明若不這樣做會(huì)存在什么樣的問題。國內(nèi)學(xué)者在對該算法進(jìn)行研究時(shí)，也沒有提到在按列掃描之后還要進(jìn)行按行掃描，相反讓人誤以為LTI算法向前處理只需按列掃描即可。但實(shí)際上，在向前處理過程中僅按列或按行掃描均沒有考慮到逆向傳播的射線，在處理速度變化劇烈的介質(zhì)時(shí)，會(huì)導(dǎo)致所追蹤出來的射線并不滿足最小旅行時(shí)。因此，筆者通過實(shí)例分析，說明LTI算法在向前處理過程中應(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長度；r為A到D的距離，滿足

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

若射線從A，B之間的D 點(diǎn)穿過到達(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)開始計(jì)算炮點(diǎn)到各網(wǎng)格邊界節(jié)點(diǎn)的最短旅行時(shí)，可由式（5）求得；向后處理是從接收點(diǎn)開始，確定所有炮點(diǎn)-接收點(diǎn)的射線路徑與各單元網(wǎng)格邊界的交點(diǎn)，可由式（6）求得［4，19］。

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

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

LTI算法在向前處理過程中若只按列或按行掃描來計(jì)算節(jié)點(diǎn)旅行時(shí)沒有考慮到逆向傳播的射線。以按列掃描為例，如圖2a所示，S點(diǎn)表示炮點(diǎn)所在位置，R點(diǎn)是其中的一個(gè)計(jì)算節(jié)點(diǎn)，位于炮點(diǎn)右方。LTI算法在向前處理過程中對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)左邊，也存在同樣的問題。這樣會(huì)導(dǎo)致在橫向速度變化劇烈的介質(zhì)中（即出現(xiàn)逆向傳播），計(jì)算出來的節(jié)點(diǎn)旅行時(shí)不是最小，向后處理中追蹤出來的射線路徑也不滿足最短旅行時(shí)原理，降低了該方法對復(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ù)及觀測系統(tǒng)如圖3a所示，圖3b是只按列掃描LTI算法追蹤的初至射線路徑，與文獻(xiàn)［4］的結(jié)果吻合。

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

為說明射線逆向傳播的問題，將圖3a的模型及觀測系統(tǒng)整體順時(shí)針旋轉(zhuǎn)90°，如圖4a所示，重新追蹤出的初至波射線路徑見圖4b。對比圖3b，4b中的射線路徑可以看出，旋轉(zhuǎn)后追蹤的1、2兩條射線路徑有誤。由表1也可以看出，旋轉(zhuǎn)后僅用按列掃描的LTI算法計(jì)算出來的1、2兩條射線的旅行時(shí)不是最小。由圖3b中的射線路徑可知，正確情況下圖4b中標(biāo)示1、2的2條射線應(yīng)該是由左側(cè)的高速層折射回來。究其原因是只按列掃描的LTI法在向前處理過程中，計(jì)算震源左側(cè)單元節(jié)點(diǎn)最小旅行時(shí)只考慮了來自右半平面的射線，并沒有考慮來自左半平面的射線。同理，只按行掃描也存在射線考慮不全的問題。因此，筆者認(rèn)為只按列或按行掃描的LTI法在向前處理過程并沒有考慮來自各個(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í)，再按行掃描的方法來重新計(jì)算各網(wǎng)格節(jié)點(diǎn)旅行時(shí)；將按行掃描計(jì)算的節(jié)點(diǎn)旅行時(shí)與按列掃描計(jì)算的節(jié)點(diǎn)旅行時(shí)相比較，取其中旅行時(shí)最小的值作為網(wǎng)格節(jié)點(diǎn)的旅行時(shí)。這樣就可以確保在向前處理過程中考慮來自各個(gè)方向的射線，即所計(jì)算出來的節(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é)果對比，取其中最小值作為節(jié)點(diǎn)的旅行時(shí)。

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

①計(jì)算各單元水平邊界及垂直邊界上各節(jié)點(diǎn)的旅行時(shí)，這時(shí)只考慮來自下邊界的射線（圖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í)最小的。

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

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

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

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í)間也就越長。因此在保證計(jì)算精度同時(shí)，如何提高計(jì)算效率也是目前需解決的問題。

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

圖5 向前處理按行掃描計(jì)算節(jié)點(diǎn)最小旅行時(shí)的實(shí)現(xiàn)過程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)的示意圖?？梢钥闯觯谄史志认嗤那闆r下，傳統(tǒng)LTI計(jì)算點(diǎn)數(shù)多。改進(jìn)后只在主節(jié)點(diǎn)上給定速度值，邊界上插入的次生節(jié)點(diǎn)速度可由單元邊界上相鄰兩主節(jié)點(diǎn)的速度線性插值得到；而原來的做法（圖7a）則需在每個(gè)節(jié)點(diǎn)都賦速度值。因此在將該方法用于層析成像時(shí)，插入次生節(jié)點(diǎn)的LTI法在不影響計(jì)算精度的前提下所需反演的參數(shù)更少，更有利于反演。

為了驗(yàn)證加入次生節(jié)點(diǎn)后LTI法的計(jì)算精度，建立一均勻速度模型。模型參數(shù)及觀測系統(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法絕對誤差最大值為0.391ms，在單元邊界插入次生節(jié)點(diǎn)后其絕對誤差最大值降為0.006 ms；在節(jié)點(diǎn)間距為1m的情況下，傳統(tǒng)的LTI法絕對誤差最大值為0.084ms，單元邊界插入次生節(jié)點(diǎn)后絕對誤差最大值為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ù)雜模型測試

為了綜合檢驗(yàn)改進(jìn)后LTI法對復(fù)雜模型的適應(yīng)能力，筆者引用了經(jīng)典的Marmousi模型對改進(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)。觀測系統(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算法追蹤出來的初至射線路徑，次生節(jié)點(diǎn)數(shù)為4×4。由射線路徑可以看出，在模型底部，射線沿著高速層走，符合最小旅行時(shí)射線路徑的規(guī)律。圖11是傳統(tǒng)的LTI法計(jì)算的各道初至旅行時(shí)減去改進(jìn)后的LTI法計(jì)算的各道初至旅行時(shí)得到的旅行時(shí)殘差曲線?？梢钥吹剑?jì)算出來的旅行時(shí)殘差均處于殘差為0的參考線上方，也就是說改進(jìn)后的LTI算法算出來的初至旅行時(shí)小于傳統(tǒng)的LTI法所計(jì)算出來的初至旅行時(shí)?？梢则?yàn)證改進(jìn)后的LTI算法較傳統(tǒng)的LTI法計(jì)算精度更高，且改進(jìn)后的LTI法對復(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é)果對比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算法僅按列或按行掃描沒有考慮到逆向傳播的射線，因此建議在向前處理時(shí)采用全方位循環(huán)的方法計(jì)算網(wǎng)格節(jié)點(diǎn)走時(shí)，這樣可以確保追蹤出來的射線路徑滿足最短旅行時(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］孫章慶，孫建國，韓復(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］孫章慶，孫建國，韓復(fù)興.針對復(fù)雜地形的三種地震波走時(shí)算法及對比［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ù)興，孫建國，楊昊.基于二維三次卷積插值算法的波前構(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］孫章慶，孫建國，韓復(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］孫章慶，孫建國，韓復(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］／／中國地球物理學(xué)會(huì)年會(huì).長沙：中國科學(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.