盧江波 方 志

(中國長沙410082湖南大學(xué)土木工程學(xué)院)

一種線性走時插值射線追蹤改進算法*

(中國長沙410082湖南大學(xué)土木工程學(xué)院)

針對線性走時插值算法(LTI)不能正確追蹤逆向傳播射線的問題， 目前已提出多種改進算法， 如擴張收縮LTI算法、 循環(huán)計算LTI算法、 動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法等， 但這些算法的計算效率普遍偏低． 在分析各種改進LTI算法的優(yōu)劣后， 本文提出了改進動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法. 該改進算法依據(jù)波的傳播規(guī)律以及LTI算法的基本方程， 排除動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法中大量冗余節(jié)點計算, 并采用傳統(tǒng)的二叉樹堆排序算法對波前陣列節(jié)點進行管理． 數(shù)值算例表明， 本文提出的改進算法具有較高的計算效率， 其計算效率是動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法的4.5—30倍， 是原始LTI算法的2—6.5倍； 當動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法采用堆排序算法時， 改進算法的計算效率是其3.5—15倍．

射線追蹤 線性走時插值 改進算法 波前擴展 計算效率

引言

射線追蹤技術(shù)在地震層析成像以及混凝土超聲波射線層析成像等領(lǐng)域具有重要作用. 目前常用的射線追蹤方法主要有兩點射線追蹤算法(包括試射法)(Julian, Gubbins, 1977； 徐濤等, 2004； 田玥， 陳曉非， 2005)、 彎曲法(Thurber, Ellsworth, 1980; Xuetal, 2006, 2010)、 有限差分解程函方程法(Vidale, 1988; Qinetal, 1992; 李振春等， 2004)、 最短路徑法(Moser, 1991; 趙愛華等, 2000； Zhaoetal, 2004; Baietal, 2010; 趙愛華， 徐濤， 2012)和LTI (linear travel-time interpolation)射線追蹤算法(Asakawa, Kawanaka, 1993; 趙改善等， 1998; Cardarelli, Cerreto, 2002; 黃靚， 黃政宇, 2002; 聶建新， 楊慧珠, 2003; 張建中等, 2003, 2004; 黃靚， 2008； 張東等, 2009; 盧江波， 方志， 2014)等． 其中， LTI射線追蹤算法因其計算精度較高、 計算速度較快， 且適用于任意復(fù)雜的速度介質(zhì)模型， 在地震層析成像等領(lǐng)域得到廣泛應(yīng)用． 但原始LTI算法(Asakawa, Kawanaka, 1993)所采用的由震源向模型邊界外推的計算方式， 不能正確追蹤逆向傳播的射線， 影響地震層析成像的精度．

針對這一問題， 不少研究人員提出了多種改進算法. 例如： 黃靚和黃政宇(2002)及黃靚(2008)提出了擴張收縮LTI算法， 盧江波和方志(2014)提出了擴張收縮LTI改進算法， 張東等(2009)提出了循環(huán)計算LTI算法； 此外， 張建中等(2003， 2004)將波前擴展方式與LTI基本方程相結(jié)合， 提出了動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法. 這些算法中， 擴張收縮LTI算法、 擴張收縮LTI改進算法以及循環(huán)計算LTI算法都需要迭代計算， 對于復(fù)雜的速度模型以及網(wǎng)格劃分較為精細的模型， 其迭代次數(shù)較多， 計算效率偏低. 動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法能夠一次計算出所有節(jié)點的走時， 具有相對較高的計算效率， 且也能有效解決上述問題， 但該算法在進行波前擴展時， 所選擇的插值線段不夠合理， 無效重復(fù)計算較多， 計算量大; 此外， 該算法采用快速排序法及插入排序法管理波前陣列節(jié)點， 效率較低．

本文針對動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法存在的不足對其進行改進. 首先依據(jù)波的傳播規(guī)律以及LTI算法的基本方程， 選擇更為合理的插值線段， 并排除該算法中的無效重復(fù)計算； 然后采用傳統(tǒng)的二叉樹堆排序算法對波前陣列節(jié)點進行管理， 以提高該算法的計算效率．

1 LTI算法的基本方程

如圖1所示， 射線經(jīng)單元下邊界的AB節(jié)段到達C節(jié)點， 射線與AB的交點為D．A點及B點走時分別為tA和tB， 節(jié)段AB長為L， 單元慢度為s，C點距A點的水平及豎向距離均為已知， 分別為xC，yC． 建立以A點為原點的局部坐標系， 確定A，B，C，D點的局部坐標． 現(xiàn)推導(dǎo)C點走時以及交點D距A點長度r的計算公式(Asakawa, Kawanaka, 1993; 趙改善等， 1998)．

由線性追蹤算法的基本假設(shè)可得D點的走時為

圖1 經(jīng)過邊界AB到達C點的射線路徑圖Fig.1 The diagram of ray path from segment AB to C

(1)

根據(jù)D點的走時， 結(jié)合單元慢度以及各點的局部坐標等條件， 可得C點的走時為

(2)

將式(1)代入式(2)得

(3)

根據(jù)費馬原理， tC對r的一階偏導(dǎo)數(shù)應(yīng)滿足等于零的條件， 即(設(shè)Δt=tB-tA)

(4)

當L2s2-Δt2>0時， 解方程(4)可得

(5)

若r≥0且r≤L, 則

(6)

若r<0或r>L， 則計算r=0和r=L時的tC值， 并取兩者較小值作為最終tC值．

當L2s2-Δt2≤0時， tC的計算方法與r<0或r>L時的情況相同．

2 LTI算法及動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法存在的不足及改進

圖2 LTI算法存在的問題 Fig.2 Problems of LTI algorithm

對于這一問題， 張建中等(2003， 2004)提出的動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法能予以有效解決. 但是該算法在進行波前擴展時， 選擇的插值線段不夠合理， 無效重復(fù)計算較多． 現(xiàn)以圖3所示均質(zhì)模型為例， 對該算法的基本步驟及其計算策略進行詳細說明， 同時依據(jù)波的傳播規(guī)律以及LTI算法的基本方程對該算法進行改進， 改進算法的基本步驟及其計算策略如圖4所示. 具體步驟及改進分析如下:

1) 首先計算震源S所在單元中所有節(jié)點的走時， 如圖3a及圖4a所示， 然后將這些節(jié)點加入波前陣列中用于下一步計算．

2) 從波前陣列中找出走時最小的節(jié)點作為當前擴展點， 假定節(jié)點為A1， 然后對A1的所有鄰點(S，A2，A12)一一進行分析， 如圖3b--d及圖4b--d所示， 并根據(jù)鄰點的不同狀態(tài)采用不同的計算策略， 具體如下：

① 鄰點未計算過的策略. 當前擴展點不與該鄰點形成插值段， 而是以當前擴展點為震源， 對單元內(nèi)還未做過擴展點的節(jié)點進行計算． 對鄰點A12進行分析時， 即采用這一策略， 如圖3b所示． 事實上， 對于鄰點未計算過的情況， 只需利用當前擴展點A1的走時信息對該鄰點進行計算(圖4b). 考慮波沿單元邊界的速度與通過單元內(nèi)部的速度不同， 將Ⅱ號單元中的節(jié)點分為A1A4、A1A14邊界上的節(jié)點以及其它節(jié)點； 對于其它節(jié)點， 顯然在A12做擴展點， 并與A1形成插值段時進行計算更為合理有效； 對于A1A4、A1A14邊界上的節(jié)點， 與最短路徑法的計算方法相同， 只需計算節(jié)點A2以及A12， 而節(jié)點A2的計算則可放在對鄰點A2的分析中, 因此， 此時只需對節(jié)點A12進行計算．

② 鄰點已計算過但未做過擴展點(不確定該鄰點是否已得到理論最小走時)的策略. 當前擴展點與該鄰點形成走時插值段， 并應(yīng)用LTI基本方程計算插值段所在單元中還未做過擴展點的節(jié)點． 對鄰點A2進行分析時， 即采用這一策略， 具體如圖3c所示． 事實上， 對于鄰點已計算過但未做過擴展點的情況， 只需利用當前擴展點A1的走時信息對該鄰點進行計算(圖4c). 由LTI基本方程式(3)可知， 在鄰點A2計算過但未做過擴展點的情況下， 節(jié)點A1與A2形成插值段的意義不明確． 而且節(jié)點A2做擴展點時， 節(jié)點A2與A1也會形成插值段， 此時， 節(jié)點A1、A2均確定得到最小走時． 顯然， 與未做過擴展點的鄰點A2形成插值段既不合理也無必要． 但考慮節(jié)點A2的理論射線路徑可能經(jīng)過節(jié)點A1(如Ⅱ號單元的速度遠大于Ⅰ號單元)， 因此， 在不形成插值段的情況下， 還需利用節(jié)點A1的走時信息對節(jié)點A2進行計算．

③ 鄰點已做過擴展點(可以確定該鄰點已得到理論最小走時)的策略. 當前擴展點與該鄰點形成走時插值段， 并應(yīng)用LTI基本方程計算插值段所在單元中還未做過擴展點的節(jié)點． 對鄰點S進行分析時， 即采用這一策略， 具體如圖3d所示． 實際上， 在當前擴展點與已做過擴展點的鄰點形成插值段時， 并不需要對插值段所在單元中所有未做過擴展點的節(jié)點進行計算， 此時僅需考慮位于當前擴展點一側(cè)且不與插值段在同一邊界的節(jié)點(圖4d). 而且在對這些節(jié)點進行計算時， 還可以通過規(guī)定節(jié)點的計算順序， 并增加一些判斷條件， 進一步排除無效重復(fù)計算. 具體為： 先計算位于插值段A1S對邊A4A7， 且與當前擴展點A1處于同一位置的節(jié)點A4， 然后再依次計算其它節(jié)點A3、A2； 在節(jié)點計算前對節(jié)點進行判斷， 以節(jié)點A4為例， 若節(jié)點已獲得的走時小于當前擴展點鄰點的走時與節(jié)點至當前擴展點的走時之和, 即tA4

a) 計算范圍的確定. 首先， 對于與插值段處于同一邊界的節(jié)點， 顯然只需考慮與當前擴展點相鄰的節(jié)點， 當鄰點已做過擴展點時(如圖3d中的節(jié)點S)， 不需要計算， 若鄰點未做過擴展點， 則可按照上面的計算策略①或②進行計算， 此時無需考慮； 然后， 對位于鄰點一側(cè)且不與插值段處于同一邊界的節(jié)點(A5—A9)， 可知在節(jié)點S已做過擴展點的情況下， 若A5—A9節(jié)點的理論射線路徑過A1S插值段， 則交點必然為節(jié)點S． 現(xiàn)對節(jié)點S的左鄰點A11進行分析， 若節(jié)點A11的理論最小走時小于或等于節(jié)點S的理論最小走時， 則節(jié)點A6—A9的理論射線路徑不過A1S插值段， 因為這些節(jié)點通過A11計算得到的走時較通過A1S插值段更??； 若A11的理論最小走時大于節(jié)點S的理論最小走時(圖3d)， 則當A11做擴展點時， 節(jié)點S必已做過擴展點，A11將與節(jié)點S形成插值段， 通過對A11一側(cè)的節(jié)點進行計算， 節(jié)點A6—A9能得到最小走時. 無論哪種情況均不用利用A1S插值段對節(jié)點A6—A9進行計算． 對于節(jié)點A5， 可放在節(jié)點S做擴展點時進行計算， 考慮節(jié)點S做擴展點時， 節(jié)點A11和A1可能均未做過擴展點， 而節(jié)點A11和A1做擴展點時， 又只考慮位于節(jié)點A11或A1一側(cè)的節(jié)點，A5將不能得到最小走時， 因此在節(jié)點S為當前擴展點、 而其鄰點A11和A1均未做過擴展點的情況下， 除對節(jié)點A11和A1進行計算外， 還需增加節(jié)點A5的計算(具體見計算策略④)．

b) 節(jié)點計算前的判斷條件. 以圖4d中的節(jié)點A4為例進行說明. 顯然， 鄰點S的走時為插值段A1S上的最小走時， 節(jié)點A4到插值段A1S的最短距離為A1A4. 根據(jù)LTI基本方程式(1)和(2)可知， 若tA4

c) 節(jié)點計算后的判斷條件. 假定A4滿足費馬原理的射線過當前擴展點A1. 顯然， 節(jié)點A3和A2滿足費馬原理的射線也過當前擴展點A1. 如果A2—A4的理論射線路徑過A1S插值段， 那么這些節(jié)點的理論射線路徑必然沿A1A4邊界傳播， 與最短路徑法的計算方法相同， 只需對A1的鄰點A2進行計算即可. 而節(jié)點A2的計算則可放在對鄰點A2的分析中. 因此， 當A4滿足費馬原理的射線過當前擴展點A1時，A4后面的節(jié)點A3和A2均不用計算．

以上分析針對的是當前擴展點為單元端節(jié)點的情況， 對于當前擴展點為中間節(jié)點的情況亦有同樣的結(jié)論. 以節(jié)點A11為當前擴展點、S為其已做過擴展點的鄰點為例， 節(jié)點的計算順序從A6到A9， 若某節(jié)點滿足費馬原理的射線過當前擴展點A11， 假定為A6， 則根據(jù)LTI基本方程可以證明，A7—A9滿足費馬原理的射線也過當前擴展點A11． 如果A7—A9的理論射線路徑過A11S插值段， 那么這些節(jié)點的理論射線路徑必然過節(jié)點A11． 現(xiàn)對節(jié)點A11的左鄰點A10進行分析. 如果節(jié)點A10的理論最小走時小于或等于節(jié)點A11的理論最小走時， 則節(jié)點A7—A9的理論射線路徑必然不過A11S插值段， 因為這些節(jié)點通過A10計算得到的走時較通過A11計算的走時更?。?如果節(jié)點A10的理論最小走時大于節(jié)點A11的理論最小走時， 那么當節(jié)點A10做擴展點時， 節(jié)點A11必然已做過擴展點，A10與A11將形成插值段， 通過對A10一側(cè)的節(jié)點進行計算，A7—A9可得到最小走時； 無論哪種情況，A6后面的節(jié)點A7—A9均不用計算．

④ 如果當前擴展點為單元的中間節(jié)點， 且其兩個鄰點均未做過擴展點， 則需利用當前擴展點的走時信息， 對位于當前擴展點對邊， 且與當前擴展點處于同一位置的節(jié)點進行計算． 以圖4d為例， 如果I號單元右邊界上的中間節(jié)點A2為當前擴展點， 且其兩個鄰點A1及A3均未做過擴展點， 則需利用A2的走時信息對A9及A15進行計算． 需要特別說明的是， 這一條計算規(guī)則是改進算法特有的. 動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法沒有這一規(guī)則．

對于以上4種計算策略， 所有節(jié)點在計算后， 都需進行如下處理： 若節(jié)點之前未被計算過， 則更新其走時并將其加入波前陣列中； 若節(jié)點之前已計算過， 而新計算出的走時小于原來的走時， 則更新其走時． 在分析完當前擴展點的所有鄰點后， 將當前擴展點從波前陣列中刪除．

3) 重復(fù)步驟2)， 直到模型中所有節(jié)點均做過擴展點， 此時所有節(jié)點均得到最小走時， 圖3e、 圖4e為最終狀態(tài)．

對比圖3與圖4可以看出， 動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法存在較多的無效重復(fù)計算， 經(jīng)改進后， 該算法的節(jié)點計算量明顯減少． 此外， 動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法采用快速排序法和插入排序法管理波前陣列節(jié)點， 并認為這種排序方法要優(yōu)于二叉樹堆排序算法， 且兩者的運行時間分別為O(n)和O(nlgn)， 式中n表示當前波前陣列的節(jié)點數(shù)(王輝， 常旭， 2000； 張建中等*張建中等(2003， 2004)在引用王輝和常旭(2000)一文時， 認為n為總節(jié)點數(shù). 實際上是誤解了原文的意思.， 2003， 2004). 但是， 這一結(jié)論的前提是波前陣列節(jié)點始終保持從小到大的排列順序． 事實上， 采用堆排序算法對波前陣列進行管理并不要求波前陣列節(jié)點保持這種順序， 而僅要求堆結(jié)構(gòu)在插入、 更新以及刪除節(jié)點后仍保持最小堆的性質(zhì)即可， 整個過程的運行時間為O(nlgn)(Cormenetal, 2009)． 因此二叉樹堆排序算法實際上要優(yōu)于插入排序算法． 當然， 除二叉樹堆排序算法外， 還存在性能更為優(yōu)越的算法， 如Fibonacci堆、 quake堆以及Brodal隊列等(Brodal, 2013)， 但這些算法較為復(fù)雜， 實用性還存在不足． 綜合考慮， 本文仍采用傳統(tǒng)的二叉樹堆排序算法對波前陣列節(jié)點進行管理．

3 數(shù)值算例

為說明原始LTI算法存在的問題, 比較各算法的計算效率, 以及驗證動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法和本文改進算法的有效性， 建立了如圖5所示的存在高速區(qū)和低速區(qū)的連續(xù)介質(zhì)模型， 其整體尺寸為2500 m×600 m. 其中深度在400—600 m之間的區(qū)域為高速區(qū)， 速度為6000 m/s； 深度在200—400 m之間且水平方向在500—2200 m的區(qū)域為低速區(qū)， 速度為500 m/s； 其余區(qū)域的速度隨深度增加而線性增加， 深度z處的速度v(z)=v0(1+βz)， 式中v0取2000 m/s，β取0.001． 震源S位于模型的左上角， 接收點R均勻布設(shè)于模型上邊界， 間距為100 m； 單元大小為5 m×5 m， 共有6×104個單元， 單元邊界劃分段數(shù)為1．

圖5 模型速度Fig.5 Velocity model

圖6 各接收點的理論射線路徑Fig.6 Theoretical ray paths of all receiving points

圖7 原始LTI算法(a)、 動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法(b)和本文改進算法(c)的射線追蹤結(jié)果Fig.7 Ray path tracing results by the three algorithms(a) Original LTI algorithm; (b) The shortest path ray tracing algorithm with dynamic networks; (c) Improved algorithm presented in this study

圖8 原始LTI算法(a)、 動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法(b)和本文改進算法(c)的走時相對誤差a圖中x<500 m的區(qū)域沒有出現(xiàn)如b， c圖所示的等高線， 是因為該區(qū)域內(nèi)的相對誤差所對應(yīng)的顏色為淺色， 未能有效顯示所致. 圖9也存在類似的情況Fig.8 Travel time relative errors by three algorithms (a) Original LTI algorithm; (b) The shortest path ray tracing algorithm with dynamic networks; (c) Improved algorithm presented in this study. The reason why the contour in the x<500 m do not presentin Fig.8a like Figs.8b or 8c is that the color which corresponds to the value of relative error in thex<500 m was too light to present in Fig.8a. The similar situation is presented in Fig.9

圖9 原始LTI算法(a)、 動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法(b)和改進算法(c)的走時絕對誤差Fig.9 Travel time absolute error by three algorithms(a) Original LTI algorithm; (b) The shortest path ray tracing algorithm with dynamic networks; (c) Improved algorithm presented in this study

對比圖7a與圖6可以看出， 原始LTI算法對于大多數(shù)接收點均能正確追蹤其射線路徑， 但當接收點的理論射線路徑存在逆向傳播的射線時(如圖6中的c，d，e，f點)， 原始LTI算法不能對其進行正確的射線追蹤． 圖7b, c中各接收點的射線路徑與圖6中的理論射線路徑一致， 表明動態(tài)網(wǎng)絡(luò)射線追蹤算法和本文改進算法均能考慮逆向傳播的射線， 能夠正確處理復(fù)雜介質(zhì)模型中接收點的射線追蹤問題．

從圖8和圖9可以看出， 原始LTI算法在存在逆向傳播射線的區(qū)域, 其走時相對誤差和絕對誤差較大， 最大分別達到43.36%和297.8 ms； 而動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法和本文改進算法則不存在這一問題， 均能夠正確計算各節(jié)點的最小走時， 其最大相對誤差和絕對誤差分別為3.96%和1.5 ms．

為對比3種算法的計算效率， 以圖5所示模型為計算對象， 單元邊界劃分段數(shù)為1， 2， 4， 10， 20和30等6種情況， 分別記錄各算法的計算時間(不包括算法的前處理過程). 其中, 動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法采用插入排序法及二叉樹堆排序法對波前節(jié)點進行管理， 本文改進算法只采用二叉樹堆排序法． 計算機CPU主頻為4.3 GHz， 對比結(jié)果如表1所示.

表1 3種算法計算效率的對比Table 1 Comparison of computational efficiency for three algorithms 單位： s

從表1可以看出， 本文改進算法具有較高的計算效率， 是動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法的4.5—30倍， 是原始LTI算法的2—6.5倍； 當動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法采用堆排序算法時， 本文改進算法的計算效率是其3.5—15倍．

圖10 本文改進算法對Marmousi速度模型的初至波計算結(jié)果 (a) 初至波等時線; (b) 相對誤差Fig.10 The first arrival calculations result based on Marmousi velocity model by using the improved algorithm presented in this study (a) Isochrons of first arrivals; (b) The relative error

4 討論與結(jié)論

動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法能夠有效解決原始LTI算法不能正確追蹤逆向傳播射線的問題， 且只需一次擴展計算就能得到所有節(jié)點的最小走時， 具有較高的計算效率. 但該算法存在較多的無效重復(fù)計算， 且其波前陣列節(jié)點的管理方法效率較低. 為此本文根據(jù)波的走時信息以及LTI基本方程提出了改進動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法， 排除了動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法中大量的節(jié)點計算， 同時采用傳統(tǒng)二叉樹堆排序法管理波前陣列節(jié)點， 較大幅度提高了該算法的計算效率. 數(shù)值算例表明， 本文提出的改進算法具有較高的計算效率， 是動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法的4.5—30倍， 是原始LTI算法的2—6.5倍； 當動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法采用堆排序算法時， 本文改進算法的計算效率是其3.5—15倍． 此外， 若將本文的改進算法用于三維模型， 則效果將更為顯著.

黃靚， 黃政宇. 2002. 線性插值射線追蹤的改進方法[J]. 湘潭大學(xué)自然科學(xué)學(xué)報， 24(4): 105--108.

Huang L, Huang Z Y. 2002. An improved method of linear interpolation ray tracing[J].NaturalScienceJournalofXiangtanUniversity, 24(4): 105--108 (in Chinese).

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

Huang L. 2008.MethodologyandExperimentResearchonConcreteUltrasonicComputerizedTomography[D]. Changsha: College of Civil Engineering, Hunan University: 33--36 (in Chinese).

李振春， 劉玉蓮， 張建磊， 馬在田， 王華忠. 2004. 基于矩形網(wǎng)格的有限差分走時計算方法[J]. 地震學(xué)報， 26(6): 644--650.

Li Z C, Liu Y L, Zhang J L, Ma Z T, Wang H Z. 2004. Finite-difference calculation of traveltimes based on rectangular grid[J].ActaSeismologicaSinica, 26(6): 644--650 (in Chinese).

盧江波, 方志. 2014. 線性走時插值射線追蹤算法的改進[J]. 湖南大學(xué)學(xué)報： 自然科學(xué)版, 41(1): 39--44.

Lu J B, Fang Z. 2014. Improvement of ray tracing algorithm based on linear traveltime interpolation[J].JournalofHunanUniversity:NaturalSciences, 41(1): 39--44 (in Chinese).

陸基孟， 王永剛. 2009. 地震勘探原理[M]. 第三版. 東營： 中國石油大學(xué)出版社: 47--53.

Lu J M, Wang Y G. 2009.ThePrincipleofSeismicExploration[M]. 3rd ed. Dongying: China University of Petroleum Press: 47--53 (in Chinese).

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

Nie J X, Yang H Z. 2003. Quadratic/linear travel time interpolation of seismic ray tracing[J].JournalofTsinghuaUniversity:Science&Technology, 43(11): 1495--1498 (in Chinese).

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

Tian Y, Chen X F. 2005. A rapid and accurate two-point ray tracing method in horizontally layered velocity model[J].ActaSeismologicaSinica, 27(2): 147--154 (in Chinese).

王輝， 常旭. 2000. 基于圖形結(jié)構(gòu)的三維射線追蹤方法[J]. 地球物理學(xué)報, 43(4)： 534--541.

Wang H, Chang X. 2000. 3D ray tracing method based on graphic structure[J].ChineseJournalofGeophysics, 43(4): 534--541 (in Chinese).

徐濤, 徐果明, 高爾根, 朱良保, 蔣先藝. 2004. 三維復(fù)雜介質(zhì)的塊狀建模和試射射線追蹤[J]. 地球物理學(xué)報， 47(6): 1118--1126.

Xu T, Xu G M, Gao E G, Zhu L B, Jiang X Y. 2004. Block modeling and shooting ray tracing in complex 3-D media[J].ChineseJournalofGeophysics, 47(6): 1118--1126 (in Chinese).

張東， 謝寶蓮， 楊艷， 傅相如， 秦前清. 2009. 一種改進的線性走時插值射線追蹤算法[J]. 地球物理學(xué)報， 52(1): 200--205.

Zhang D, Xie B L, Yang Y, Fu X R, Qin Q Q. 2009. A ray tracing method based on improved linear traveltime interpolation[J].ChineseJournalofGeophysics, 52(1): 200--205 (in Chinese).

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

Zhang J Z, Chen S J, Yu D X. 2003. Improvement of shortest path raytracing method[J].ProgressinGeophysics, 18(1): 146--150 (in Chinese).

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

Zhang J Z, Chen S J, Xu C W. 2004. A method of shortest path raytracing with dynamic networks[J].ChineseJournalofGeophysics, 47(5): 899--904 (in Chinese).

趙愛華， 張中杰， 王光杰， 王輝. 2000. 非均勻介質(zhì)中地震波走時與射線路徑快速計算技術(shù)[J]. 地震學(xué)報， 22(2): 151--157.

Zhao A H, Zhang Z J, Wang G J, Wang H. 2000. A new scheme for fast calculation of seismic traveltimes and ray paths in heterogeneous media[J].ActaSeismologicaSinica, 22(2): 151--157 (in Chinese).

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

Zhao A H, Xu T. 2012. A traveltime correction technique for improving the accuracy of reflection wave traveltimes with the shortest path method based on a regular grid[J].ProgressinGeophysics, 27(5): 1854--1862 (in Chinese).

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

Zhao G S, Hao S L, Yang E H, Chen W. 1998. Seismic ray tracing algorithm based on the linear traveltime interpolation[J].GeophysicalProspectingforPetroleum, 37(2): 14--24 (in Chinese).

Asakawa E, Kawanaka T. 1993. Seismic ray tracing using linear traveltime interpolation[J].GeophysicalProspecting, 41(1): 99--111.

Bai C Y, Huang G J, Zhao R. 2010. 2-D/3-D irregular shortest-path ray tracing for multiple arrivals and its applications[J].GeophysJInt, 183(3): 1596--1612.

Brodal G S. 2013. A survey on priority queues[G]∥Space-EfficientDataStructures,Streams,andAlgorithms. Heidelberg： Springer: 150--163.

Cardarelli E, Cerreto A. 2002. Ray tracing in elliptical anisotropic media using the linear traveltime interpolation (LTI) method applied to traveltime seismic tomography[J].GeophysicalProspecting, 50(1): 55--72.

Cormen T H, Leiserson C E, Rivest R L, Stein C. 2009.IntroductiontoAlgorithms[M]. 3rd ed. London: MIT Press: 151--169.

Julian B R, Gubbins D. 1977. Three-dimensional seismic ray tracing[J].JGeophys, 43(1/2): 95--113.

Moser T J. 1991. Shortest path calculation of seismic rays[J].Geophysics, 56(1): 59--67.

Qin F, Luo Y, Olsen K B, Cai W Y, Schuster G T. 1992. Finite-difference solution of the eikonal equation along expanding wavefronts[J].Geophysics, 57(3): 478--487.

Thurber C H, Ellsworth W L. 1980. Rapid solution of ray tracing problems in heterogeneous media[J].BullSeismolSocAm, 70(4): 1137--1148.

Vidale J. 1988. Finite-difference calculation of travel times[J].BullSeismolSocAm, 78(6): 2062--2076.

Xu T, Xu G M, Gao E G, Li Y C, Jiang X Y, Luo K Y. 2006. Block modeling and segmentally iterative ray tracing in complex 3D media[J].Geophysics, 71(3): T41--T51.

Xu T, Zhang Z, Gao E, Xu G, Sun L. 2010. Segmentally iterative ray tracing in complex 2D and 3D heterogeneous block models[J].BullSeismolSocAm, 100(2): 841--850.

Zhao A H, Zhang Z J, Teng J W. 2004. Minimum travel time tree algorithm for seismic ray tracing: Improvement in efficiency[J].JGeophysEng, 1(4): 245--251.

An improved ray-tracing algorithm based on linear travel-time interpolation

(CollegeofCivilEngineering,HunanUniversity,Changsha410082,China)

In order to solver for the problem that the original LTI algorithm could not trace the reverse propagation ray, several linear travel-time interpolation (LTI for short) improved algorithms, such as extension-compaction LTI algorithm, loop computation LTI algorithm, the shortest path ray tracing algorithm with dynamic networks, have been presented, but the computational efficiency of these algorithms are low. After analyzing these improved algorithms, this paper presented a new improved shortest path ray tracing algorithm with dynamic networks. According to the law of wave propagation and the basic equation of LTI, a large number of redundancy node calculation are excluded, and the traditional binary heap sort algorithm was used to manage node of wavefront array. The numerical examples show that, the improved algorithm presented in this paper has the highest computational efficiency among all of improved algorithms; its calculation efficiency is about 4.5—30 times of the shortest path ray tracing algorithm with dynamic networks, and about 2—6.5 times of the original LTI algorithm, and about 3.5—15 times of the shortest path ray tracing algorithm with dynamic networks when the traditional binary heap sort algorithm is also used.

ray tracing; linear traveltime interpolation; improved algorithm; wavefront expansion; computational efficiency

國家自然科學(xué)基金(51278182, 51408213)資助.

2013-12-27收到初稿， 2014-06-30決定采用修改稿.

e-mail: fangzhi@hnu.edu.cn

10.3969/j.issn.0253-3782.2014.06.010.

10.3969/j.issn.0253-3782.2014.06.010

P315.3+1

A

盧江波, 方志. 2014. 一種線性走時插值射線追蹤改進算法. 地震學(xué)報, 36(6): 1089--1100.

Lu J B, Fang Z. 2014. An improved ray-tracing algorithm based on linear travel-time interpolation.ActaSeismologicaSinica, 36(6): 1089--1100. doi:10.3969/j.issn.0253-3782.2014.06.010.