動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法中向后追蹤方法的改進(jìn)

盧江波+方志

摘 要：動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法中的向后追蹤方法能夠解決線性走時(shí)插值算法（LTI）向后追蹤過程不穩(wěn)定的問題，但是其計(jì)算效率較低.綜合利用節(jié)點(diǎn)次級(jí)源的位置信息以及波的傳播規(guī)律，提出了改進(jìn)方法，排除了動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法向后追蹤過程中存在的大量冗余計(jì)算.數(shù)值算例表明，改進(jìn)的向后追蹤方法具有較高的計(jì)算效率，是動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法中向后追蹤方法的幾倍至幾十倍；若將改進(jìn)后的向后追蹤方法應(yīng)用于動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤改進(jìn)算法，則該算法的計(jì)算效率將提高一倍左右.

關(guān)鍵詞：射線追蹤；線性走時(shí)插值；向后追蹤方法；計(jì)算效率；初至波射線追蹤

中圖分類號(hào)：P631 文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：1674-2974（2016）05-0106-07

Abstract：The backward tracing method of the shortest path ray tracing algorithm with dynamic networks can solve the unstability problem in the backward tracing procedure of the LTI （Linear Travel-time Interpolation） algorithm， but the computational efficiency of the method is low. This study presented an improved method on backward tracing. According to the location information of the secondary sources for the nodes and the law of wave propagation， a large number of redundancy calculation are excluded in the backward tracing of the dynamic networks tracing algorithm. The numerical examples show that the improved method exhibits the higher computational efficiency. The calculation efficiency of the improved method is several times that of the backward tracing method of the dynamic networks tracing algorithm. When the improved method is applied to the improved algorithm of the shortest path ray tracing with dynamic networks， the computational efficiency of the algorithm can be increased by about 100 %.

Key words：ray tracing； linear traveltime interpolation； improved algorithm； backward tracing； computational efficiency；first arrival ray tracing

射線追蹤技術(shù)在地震層析成像以及混凝土超聲波射線層析成像等領(lǐng)域具有重要作用.目前常用射線追蹤方法主要有兩點(diǎn)射線追蹤算法（包括試射法以及彎曲法）[[1-3]、有限差分解程函方程法[[4-6]、最短路徑法[[7-10]以及LTI（Linear Travel-time Interpolation）射線追蹤算法[[11-22]等.其中，LTI射線追蹤算法因其計(jì)算精度較高、計(jì)算速度較快且適用于任意復(fù)雜的速度介質(zhì)模型，在地震層析成像等領(lǐng)域得到了廣泛應(yīng)用.但是LTI原算法[[11]存在兩個(gè)問題：在向前計(jì)算節(jié)點(diǎn)最小走時(shí)時(shí)，不能正確追蹤逆向傳播的射線，相關(guān)節(jié)點(diǎn)不能得到正確的最小走時(shí)[[16-22]；在向后追蹤接收點(diǎn)射線路徑時(shí)，存在不能正確追蹤接收點(diǎn)射線路徑的可能，算法的穩(wěn)定性存在不足.

文獻(xiàn)[17-18]將波前擴(kuò)展方式與LTI算法基本方程相結(jié)合，提出了動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法，該算法在向前計(jì)算各節(jié)點(diǎn)的最小走時(shí)時(shí)，從震源點(diǎn)開始，采用波前擴(kuò)展的方式逐點(diǎn)計(jì)算各個(gè)節(jié)點(diǎn)上的最小走時(shí)，改變了LTI原算法的計(jì)算方式，確保各節(jié)點(diǎn)能得到最小走時(shí)；在向后追蹤接收點(diǎn)射線路徑時(shí)，基于互換原理考慮了接收點(diǎn)所有可能的射線路徑，確保算法能正確追蹤接收點(diǎn)的射線路徑.

動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法雖然能夠解決LTI原算法存在的兩個(gè)問題，但是其計(jì)算效率偏低.文獻(xiàn)[22]基于波的傳播規(guī)律提出了動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤改進(jìn)算法，改進(jìn)并提高了動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法向前計(jì)算節(jié)點(diǎn)最小走時(shí)這一步驟的計(jì)算效率，但是，該改進(jìn)方法的向后追蹤方法仍然采用動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法中的方法，計(jì)算效率依然偏低.

針對(duì)動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法中的向后追蹤方法存在計(jì)算效率低的問題，本文提出了改進(jìn)方法.首先，在向前計(jì)算節(jié)點(diǎn)最小走時(shí)這一步驟中，不僅計(jì)算各節(jié)點(diǎn)的最小走時(shí)，而且還記錄各節(jié)點(diǎn)次級(jí)源的位置；然后，在向后追蹤接收點(diǎn)的射線路徑時(shí)，利用各節(jié)點(diǎn)次級(jí)源的位置信息以及波的傳播規(guī)律對(duì)算法進(jìn)行改進(jìn)，降低算法的計(jì)算量，提高算法的計(jì)算效率.

1 LTI算法基本方程的推導(dǎo)及接收點(diǎn)

文獻(xiàn)[8]給出了最短路徑射線追蹤算法中次級(jí)源的確定方法，考慮到LTI算法中的射線可通過單元邊界上任意點(diǎn)進(jìn)行傳播，這與最短路徑射線追蹤算法中射線只能通過節(jié)點(diǎn)進(jìn)行傳播不同，因此最短路徑射線追蹤算法中次級(jí)源的確定方法不完全適合LTI算法.為了更好地描述本文提出的改進(jìn)方法，將LTI算法中次級(jí)源的確定方法規(guī)定如下：通過接收點(diǎn)所在單元各節(jié)段（不包括接收點(diǎn)所在節(jié)段）可得到多個(gè)接收點(diǎn)走時(shí)，其中走時(shí)最小值對(duì)應(yīng)的點(diǎn)即為接收點(diǎn)的次級(jí)源，如果走時(shí)最小值對(duì)應(yīng)多個(gè)點(diǎn)，取距離接收點(diǎn)最近的點(diǎn)作為接收點(diǎn)的次級(jí)源.如圖1所示，若接收點(diǎn)C通過AB節(jié)段D點(diǎn)得到的走時(shí)小于通過單元其它節(jié)段得到的走時(shí)，那么D點(diǎn)即為接收點(diǎn)C的次級(jí)源.

2 LTI原始算法向后追蹤過程存在的問題及動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法存在的不足及改進(jìn)

LTI原始算法在向后追蹤接收點(diǎn)射線路徑時(shí)，首先逐點(diǎn)計(jì)算接收點(diǎn)所在單元中各節(jié)點(diǎn)走時(shí)與節(jié)點(diǎn)至接收點(diǎn)的走時(shí)之和，然后選出最小的走時(shí)之和以及相應(yīng)的節(jié)點(diǎn)，最后將單元中包含該節(jié)點(diǎn)的節(jié)段作為接收點(diǎn)次級(jí)源的可能區(qū)域[[11].事實(shí)上，接收點(diǎn)的次級(jí)源并不一定在這些節(jié)段中.此外，LTI原算法在確定接收點(diǎn)次級(jí)源的可能區(qū)域時(shí)，并未排除接收點(diǎn)所在的節(jié)段，算法可能會(huì)陷入無限循環(huán).因?yàn)橥ㄟ^接收點(diǎn)所在節(jié)段線性插值得到的接收點(diǎn)走時(shí)，可能比通過單元其它節(jié)段得到的走時(shí)更小，那么接收點(diǎn)取最小走時(shí)對(duì)應(yīng)的點(diǎn)可能為接收點(diǎn)本身，算法可能會(huì)進(jìn)入無限循環(huán).

以圖2所示模型為例，模型尺寸為3 m×3 m，單元大小為1 m×1 m，單元邊界劃分為2個(gè)節(jié)段，模型上層、中層以及下層單元的速度分別為525 m/s，515 m/s和505 m/s，震源S的x，y坐標(biāo)分別為1.5和3.0，接收點(diǎn)R的x，y坐標(biāo)分別為2.05和0；圖中虛線為采用LTI原始算法中的向后追蹤方法得到的射線路徑，需要說明的是，計(jì)算這條射線路徑時(shí)，已經(jīng)排除了接收點(diǎn)的次級(jí)源位于接收點(diǎn)所在節(jié)段的情況，實(shí)線為根據(jù)動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法中的向后追蹤方法[[17-18]得到的射線路徑.

從圖2可以看出，在III號(hào)單元中，兩種方法的計(jì)算結(jié)果一樣，接收點(diǎn)R的次級(jí)源均為R1.在II號(hào)單元中，對(duì)于LTI原始算法的向后追蹤方法，首先須確定II號(hào)單元中節(jié)點(diǎn)走時(shí)與節(jié)點(diǎn)至接收點(diǎn)R1走時(shí)之和最小的節(jié)點(diǎn)，經(jīng)計(jì)算確定為d節(jié)點(diǎn)，然后確定cd或de節(jié)段為接收點(diǎn)R1次級(jí)源的可能區(qū)域；對(duì)于動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法中的向后追蹤方法，接收點(diǎn)R1次級(jí)源所在區(qū)域?yàn)閑f節(jié)段.計(jì)算結(jié)果表明：通過ef節(jié)段計(jì)算得到的R1點(diǎn)的走時(shí)（T=0.005 416 05 s）要小于通過cd或de節(jié)段的走時(shí)（T=0.005 422 57 s），接收點(diǎn)R1的次級(jí)源不在cd或de節(jié)段內(nèi).造成這一問題的根源在于：LTI算法假定射線路徑可以經(jīng)過單元邊界的任意一點(diǎn)，在向后追蹤過程中必須考慮所有可能的射線路徑，并根據(jù)費(fèi)馬原理選擇走時(shí)最短的那條路徑[[11]，而LTI原始算法的具體追蹤方法卻沒有考慮接收點(diǎn)所在單元所有節(jié)段的射線，而是采用了一種“簡(jiǎn)化”方法確定接收點(diǎn)射線路徑或次級(jí)源的可能區(qū)域，由圖2所示模型可以看出，這種“簡(jiǎn)化”方法存在不足.

此外，在確定接收點(diǎn)次級(jí)源時(shí)，若不排除接收點(diǎn)所在的節(jié)段，那么在計(jì)算新接收點(diǎn)R1的次級(jí)源時(shí)，由于通過節(jié)段cd插值得到的走時(shí)，比通過ed節(jié)段得到的走時(shí)小，計(jì)算得到的“次級(jí)源”為R1本身，程序?qū)⑾萑霟o限循環(huán).

動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法中基于互換原理提出的向后追蹤方法[[17-18]，考慮了來自接收點(diǎn)所在單元中除接收點(diǎn)所在節(jié)段外所有節(jié)段的射線，確保了接收點(diǎn)能得到其次級(jí)源.但是，該方法計(jì)算量大，計(jì)算效率低.為了解決這一問題，本文首先在向前計(jì)算節(jié)點(diǎn)最小走時(shí)的步驟中，建立了一個(gè)數(shù)組，專門用于記錄節(jié)點(diǎn)次級(jí)源的位置信息.然后利用各節(jié)點(diǎn)次級(jí)源的位置信息以及波的傳播規(guī)律對(duì)動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法中的向后追蹤方法進(jìn)行改進(jìn).現(xiàn)以圖3所示模型為例，對(duì)該向后追蹤方法的基本步驟及計(jì)算策略進(jìn)行說明，同時(shí)對(duì)向后追蹤改進(jìn)方法進(jìn)行闡述，改進(jìn)前后的向后追蹤方法及計(jì)算策略分別如圖3和圖4所示.具體步驟和分析如下：

1）首先將接收點(diǎn)分為以下3種情況，然后對(duì)不同的情況采用不同的策略求解接收點(diǎn)的次級(jí)源：

①接收點(diǎn)位于單元內(nèi)部，如圖3（a）及圖4（a）所示的R點(diǎn).此時(shí)利用LTI算法的相關(guān)公式計(jì)算接收點(diǎn)所在單元各節(jié)段至接收點(diǎn)的走時(shí)，然后從中選出走時(shí)最小值對(duì)應(yīng)的點(diǎn)，如圖3（b）及圖4（b）中的R1點(diǎn)，這個(gè)點(diǎn)就是接收點(diǎn)的次級(jí)源.

②接收點(diǎn)位于單元邊界上，但是非單元邊界上的節(jié)點(diǎn).如圖3（b），（d）中的R1和R3點(diǎn)，此時(shí)除接收點(diǎn)所在的節(jié)段外，單元中的其它節(jié)段均要計(jì)算從該節(jié)段至接收點(diǎn)的最小走時(shí)，然后從中選出走時(shí)最小值對(duì)應(yīng)的點(diǎn)，如圖3（c），（e）中的R2和S點(diǎn).

事實(shí)上，對(duì)于接收點(diǎn)位于單元邊界但非單元節(jié)點(diǎn)的情況，除非該接收點(diǎn)為原始接收點(diǎn)（圖3和圖4中的R），否則并不需要對(duì)接收點(diǎn)進(jìn)行全方位的計(jì)算.可以采取如下兩個(gè)步驟確定接收點(diǎn)次級(jí)源的可能區(qū)域.

首先，利用接收點(diǎn)所在節(jié)段端點(diǎn)的次級(jí)源位置信息確定兩個(gè)定位點(diǎn).定位點(diǎn)的確定方法為：若節(jié)點(diǎn)端點(diǎn)的次級(jí)源不在單元內(nèi)，或者次級(jí)源在單元內(nèi)但與節(jié)點(diǎn)端點(diǎn)處于同一邊界，則定位點(diǎn)為節(jié)段端點(diǎn)本身，其它情況定位點(diǎn)為節(jié)點(diǎn)端點(diǎn)的次級(jí)源.如圖4（b）所示，接收點(diǎn)R1所在節(jié)段的端點(diǎn)為A1和A2，易知接收點(diǎn)R1的次級(jí)源位于單元I中，現(xiàn)以端點(diǎn)A1的次級(jí)源位置信息為例說明定位點(diǎn)的確定方法，若A1的次級(jí)源不在單元I中，比如A1的次級(jí)源為單元II中的b點(diǎn)，則由A1的次級(jí)源位置信息確定的定位點(diǎn)為A1本身；若A1的次級(jí)源與A1處于同一邊界，比如次級(jí)源位于A1A4或者A1A6邊界，則定位點(diǎn)為A1本身；對(duì)于其它情況，比如A1的次級(jí)源為單元I中的c點(diǎn)，則定位點(diǎn)為A1的次級(jí)源c.在圖4（b）中，假定A1節(jié)點(diǎn)的次級(jí)源為A3節(jié)點(diǎn)，A2節(jié)點(diǎn)的次級(jí)源為a1點(diǎn)，則由A1，A2節(jié)點(diǎn)的次級(jí)源位置信息可以確定兩個(gè)定位點(diǎn)，分別為A1節(jié)點(diǎn)本身以及點(diǎn)a1.

然后沿著單元邊界連接兩個(gè)定位點(diǎn)，其中不包含接收點(diǎn)的那條路徑即為接收點(diǎn)次級(jí)源的可能區(qū)域.如圖4（b）所示，沿單元邊界連接定位點(diǎn)A1和a1可以得到兩條路徑：A1A3A4a1和A1A6A7a1，其中路徑A1A3A4a1不包含接收點(diǎn)R1，因此確定該路徑為接收點(diǎn)次級(jí)源的可能區(qū)域.

得到次級(jí)源的可能區(qū)域后，再確定這個(gè)區(qū)域內(nèi)的節(jié)段，圖4（b）中可能區(qū)域A1A3A4a1內(nèi)的節(jié)段為A1A3，A3A4，A4A5，計(jì)算從這些節(jié)段至接收點(diǎn)的最小走時(shí)，其中走時(shí)最小的點(diǎn)即為接收點(diǎn)的次級(jí)源，如圖4（c）中的R2.

同理得到圖4（d）中R3點(diǎn)的次級(jí)源S.

③接收點(diǎn)位于單元邊界上，且為單元邊界上的節(jié)點(diǎn).如圖3（c）中的R2節(jié)點(diǎn)，此時(shí)需要在R2節(jié)點(diǎn)所處的幾個(gè)單元中執(zhí)行情況②的計(jì)算，如圖3（c）.

顯然，由于在向前計(jì)算節(jié)點(diǎn)最小走時(shí)的過程中，已經(jīng)記錄了各節(jié)點(diǎn)的次級(jí)源，因此，在改進(jìn)方法中，可以直接得到接收點(diǎn)的次級(jí)源，如圖4（c）.

2）以步驟1）中獲得的次級(jí)源為新的接收點(diǎn)，重復(fù)步驟1），直至新的接收點(diǎn)為震源.

3）依次連接各接收點(diǎn)得到R點(diǎn)的初至波射線路徑，如圖3（e）和圖4（e）所示.

對(duì)比圖3與圖4可以看出，改進(jìn)后的向后追蹤方法計(jì)算量明顯減少.

3 數(shù)值算例

為了對(duì)比改進(jìn)前后向后追蹤方法的計(jì)算效率，建立了尺寸為2 500 m×600 m的二維模型，如圖5所示，其中x值在500～2 000 m之間且y值在200～400 m之間的區(qū)域?yàn)榈退賲^(qū)，速度為500 m/s，其余速度均為4 000 m/s，震源S位于模型上表面的正中間（1 250，0），單元尺寸為5 m×5 m，單元邊界劃分段數(shù)為4，10，20等3種情況，模型下表面的每個(gè)單元布置一個(gè)接收點(diǎn)，共500個(gè)，部分接收點(diǎn)的射線路徑如圖6所示，分別記錄兩種向后追蹤方法追蹤所有接收點(diǎn)的射線路徑所耗費(fèi)的總時(shí)間.此外，為了說明向后追蹤方法計(jì)算效率的提高對(duì)整個(gè)算法的影響，將改進(jìn)后的向后追蹤方法應(yīng)用于動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤改進(jìn)算法[[22]，然后比較應(yīng)用改進(jìn)方法前后動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤改進(jìn)算法的總耗時(shí)（不包括算法的前處理過程）.計(jì)算機(jī)CPU主頻為3.4 GHz，計(jì)算結(jié)果如表1所示.

由表1可知，向后追蹤改進(jìn)方法的計(jì)算效率較高，是改進(jìn)前向后追蹤方法的幾倍至幾十倍.并且能將動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤改進(jìn)算法的計(jì)算效率提高1倍左右.

為驗(yàn)證向后追蹤改進(jìn)方法對(duì)復(fù)雜模型的有效性，采用向后追蹤改進(jìn)方法對(duì)Marmousi速度模型進(jìn)行射線追蹤，其中節(jié)點(diǎn)最小走時(shí)的計(jì)算采用的是動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤改進(jìn)算法.模型尺寸為9 192 m×2 904 m，震源S位于模型左上角（0，0），單元尺寸為24 m×24 m，單元邊界劃分為10段，共設(shè)置5個(gè)接收點(diǎn)R1～R5，分別位于模型上表面的3 600，4 800 m，6 000 m，7 200 m和8 400 m，射線追蹤結(jié)果如圖7（a）所示，作為對(duì)照，本文給出了單元邊界劃分為30段時(shí)最短路徑法的射線追蹤結(jié)果，如圖7（b）所示.兩種算法下，各接收點(diǎn)根據(jù)射線追

蹤結(jié)果計(jì)算的走時(shí)如表2所示.計(jì)算結(jié)果表明，向后追蹤改進(jìn)方法對(duì)于復(fù)雜的速度模型同樣有效.

4 結(jié) 論

針對(duì)LTI原算法中向后追蹤方法存在的無限循環(huán)以及可能不能得到正確射線路徑的問題，動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法中基于互換原理提出的向后追蹤方法能夠予以有效的解決.但是該算法存在較多的無效計(jì)算.本文根據(jù)模型中節(jié)點(diǎn)次級(jí)源的位置信息以及波的傳播規(guī)律，提出了改進(jìn)的向后追蹤方法.數(shù)值結(jié)果表明，改進(jìn)后的向后追蹤方法，其計(jì)算效率較之改進(jìn)前的方法有較大程度的提高；此外，本文提出的向后追蹤改進(jìn)方法能將動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤改進(jìn)算法的計(jì)算效率提高1倍左右.

參考文獻(xiàn)

[1] JULIAN B R，GUBBINS D.Three-dimensional seismic ray tracing[J].J Geophys，1977，43（1/2）：95-113.

[2] 田玥，陳曉非.水平層狀介質(zhì)中的快速兩點(diǎn)間射線追蹤方法[J]. 地震學(xué)報(bào)，2005，27（2）：147-154.

[3] XU Tao，ZHANG Zhong-jie，GAO Ergen，et al. Segmentally iterative ray tracing in complex 2D and 3D heteroge-neous block models[J]. Bulletin of the Seismological Society of America， 2010，100（2）： 841-850.

[4] VIDALE J.Finite-difference calculation of travel times[J].Bulletin of the Seismological Society of America， 1988，78（6）： 2062-2076.

[5] QIN Fu-hao， LUO Yi， OLSENl K B， et al. Finite-difference solution of the eikonal equation along expanding wavefronts[J].Geophysics， 1992，57（3）：478-487.

[6] 李振春，劉玉蓮，張建磊，等.基于矩形網(wǎng)格的有限差分走時(shí)計(jì)算方法[J]. 地震學(xué)報(bào)， 2004，26（6）：644-650.

[7] MOSER T J. Shortest path calculation of seismic rays[J]. Geophysics， 1991，56（1）： 59-67.

[8] 劉洪，孟凡林，李幼銘.計(jì)算最小走時(shí)和射線路徑的界面網(wǎng)全局方法[J].地球物理學(xué)報(bào)， 1995，38（6）： 823-832.

[9] 趙愛華，徐濤.提高規(guī)則網(wǎng)格最短路徑方法反射波走時(shí)計(jì)算精度的走時(shí)校正技術(shù)[J]. 地球物理學(xué)進(jìn)展， 2012，27（5）：1854-1862.

[10]BAI Chao-yin，HUANG Guo-jiao，ZHAO Rui. 2-D/3-D irregular shortest-path ray tracing for multiple arrivals and its applications[J]. Geophysical Journal International， 2010，183（3）： 1596-1612.

[11]ASAKAWA E，KAWANAKA T. Seismic ray tracing using linear traveltime interpolation[J]. Geophysical Prospecting， 1993，41（1）： 99-111.

[12]趙改善，郝守玲，楊爾皓，等.基于旅行時(shí)線性插值的地震射線追蹤算法[J].石油物探， 1998，37（2）：14-24.

[13]CARDARELLI E，CERRETO A. Ray tracing in elliptical anisotropic media using the linear traveltime interpolation （LTI） method applied to traveltime seismic tomography[J]. Geophysical Prospecting， 2002，50（1）： 55-72.

[14]聶建新，楊慧珠.地震波旅行時(shí)二次/線性聯(lián)合插值法[J].清華大學(xué)學(xué)報(bào)：自然科學(xué)版 ， 2003，43（11）：1495-1498.

NIE Jian-xin， YANG Hui-zhu.Quadratic/linear travel timeinterpolationof seismic ray tracing[J].J Tsinghua Univ ：Sci&Tech， 2003，43（11）：1495-1498.（In Chinese）

[15]ZHANG Jian-zhong， HUANG Yue-qin， SONG Lin-ping，et al.Fast and accurate 3-D ray tracing using bilinear traveltime interpolation and the wave front group marching[J]. Geophysical Journal International， 2011，184（3）： 1327-1340.

[16]黃靚，黃政宇.線性插值射線追蹤的改進(jìn)方法[J]. 湘潭大學(xué)自然科學(xué)學(xué)報(bào)， 2002，24（4）： 105-108.

[17]張建中，陳世軍，余大祥.最短路徑射線追蹤方法及其改進(jìn)[J].地球物理學(xué)進(jìn)展， 2003，18（1）：146-150.

[18]張建中，陳世軍，徐初偉.動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤[J]. 地球物理學(xué)報(bào)， 2004，47（5）： 899-904.

[19]黃靚.混凝土超聲波層析成像的理論方法和試驗(yàn)研究[D]. 長(zhǎng)沙： 湖南大學(xué)土木工程學(xué)院， 2008：33-36.

[20]張東，謝寶蓮，楊艷，等.一種改進(jìn)的線性走時(shí)插值射線追蹤算法[J]. 地球物理學(xué)報(bào)， 2009，52（1）： 200-205.

[21]盧江波，方志.線性走時(shí)插值射線追蹤算法的改進(jìn)[J].湖南大學(xué)學(xué)報(bào)：自然科學(xué)版， 2014，41（1）：39-44.

[22]盧江波，方志.一種線性走時(shí)插值射線追蹤改進(jìn)算法[J].地震學(xué)報(bào)， 2014，36（6）：1089-1100.