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

盧江波+方志

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1）首先將接收點分為以下3種情況，然后對不同的情況采用不同的策略求解接收點的次級源：

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

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

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

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

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

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

同理得到圖4（d）中R3點的次級源S.

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

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

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

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

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

3 數(shù)值算例

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

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

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

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

4 結(jié) 論

針對LTI原算法中向后追蹤方法存在的無限循環(huán)以及可能不能得到正確射線路徑的問題，動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤算法中基于互換原理提出的向后追蹤方法能夠予以有效的解決.但是該算法存在較多的無效計算.本文根據(jù)模型中節(jié)點次級源的位置信息以及波的傳播規(guī)律，提出了改進(jìn)的向后追蹤方法.數(shù)值結(jié)果表明，改進(jìn)后的向后追蹤方法，其計算效率較之改進(jìn)前的方法有較大程度的提高；此外，本文提出的向后追蹤改進(jìn)方法能將動態(tài)網(wǎng)絡(luò)最短路徑射線追蹤改進(jìn)算法的計算效率提高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ì)中的快速兩點間射線追蹤方法[J]. 地震學(xué)報，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)格的有限差分走時計算方法[J]. 地震學(xué)報， 2004，26（6）：644-650.

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

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

[9] 趙愛華，徐濤.提高規(guī)則網(wǎng)格最短路徑方法反射波走時計算精度的走時校正技術(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]趙改善，郝守玲，楊爾皓，等.基于旅行時線性插值的地震射線追蹤算法[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]聶建新，楊慧珠.地震波旅行時二次/線性聯(lián)合插值法[J].清華大學(xué)學(xué)報：自然科學(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é)報， 2002，24（4）： 105-108.

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

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

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

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

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

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