線性走時插值射線追蹤算法的改進＊

盧江波，方 志

（湖南大學 土木工程學院，湖南 長沙 410082）

射線追蹤技術(shù)在地震層析成像以及混凝土超聲波射線層析成像等領(lǐng)域具有重要作用.目前常用射線追蹤方法主要有有限差分解程函方程法［1－2］、最短路徑法［3－4］以及LTI（Linear Travel－time Inter－polation）射線追蹤算法［5－6］等.實驗表明［5］，LTI算法在走時計算以及射線路徑追蹤上比其它方法（如有限差分解程函方程法）更為快速、精確.

LTI射線追蹤算法由Asakawa等人提出［5］，該算法以走時線性變化為前提，分兩步進行.第一步，向前計算最小走時：先將模型劃分為若干規(guī)則的單元，并將各單元邊界劃分為若干節(jié)段，然后，根據(jù)走時線性變化的假定以及走時最小原理（Fermat原理），求出從發(fā)射點到接收點的最小走時；第二步，根據(jù)得到的最小走時，反向追蹤射線路徑.但是該算法在向前計算走時時，沒有考慮射線的逆向傳播，不能追蹤回波［9］，進而影響射線追蹤的精度.針對這一問題，不少學者提出了改進算法，張東等人提出了循環(huán)計算LTI改進算法［7］，王浩全提出了交叉掃描LTI改進算法［8］，黃靚等人提出了擴張－收縮掃描改進算法［9－10］.在考慮射線的逆向傳播時，文獻［7－8］采用從發(fā)射點所在列向外逐列逆向掃描的方式，文獻［9－10］則采用從邊界列（行）向發(fā)射點所在列（行）進行逆向掃描的方式.在考慮射線的逆向傳播上，文獻［9－10］提出的改進算法更為合理，但是該算法存在計算效率低、收斂速度慢的問題.因為在LTI算法中，接收點要得到其最小走時，需其理論射線路徑與模型單元邊界所有相交節(jié)點都已得到其最小走時，這些相交節(jié)點的最小走時可能來自列掃描，也可能來自行掃描.根據(jù)本文作者對文獻［5，9－10］算法步驟的理解以及對文獻［9－10］中數(shù)值算例的研究，發(fā)現(xiàn)文獻［9－10］不僅增加了從邊界進行逆向傳播的射線，即文獻［9－10］中的收縮掃描過程，而且將文獻［5］中既考慮鄰列也考慮來自鄰行入射射線的逐列掃描方式改為逐列（行）掃描只考慮鄰列（行）入射射線的方式，這點可由文獻［9］算例模型1的計算結(jié)果得出：若文獻［9］的逐列（行）掃描過程既考慮鄰列又考慮鄰行的入射射線，則模型1經(jīng)一次擴張掃描即得接收點的最小走時，此時相對誤差應為0.247‰，而不應為按行列分開掃描方式經(jīng)一次擴張掃描得到的計算結(jié)果3.3‰①.文獻［9］采用行列分開掃描的方式使得文獻［9－10］在能正確追蹤回波的同時，接收點為獲取其最小走時也進行了較多的無效掃描，降低了算法的計算效率.由以上分析可知，若采用行列交叉掃描的方式對擴張收縮掃描算法進行改進，將有效提高擴張收縮掃描算法的計算效率.

同時，由于繞射波以及回波的存在，文獻［5］算法步驟3對每列都進行水平邊界節(jié)點最小走時搜索，其意義并不明確.此外，由于擴張收縮掃描算法相比文獻［5］增加了收縮掃描過程，在逐列掃描過程中能夠考慮上行或下行首波的最小走時，所以文獻［5］中的逐行掃描可以省略.最后，由于存在收縮掃描過程，在計算豎向邊界各節(jié)點最小走時時，若按文獻［5］算法步驟5的計算方法，則存在著重復無效掃描.因此，也有必要對文獻［5］的擴張掃描過程進行簡化和改進.

基于此，本文在前人研究的基礎(chǔ)上提出了基于交叉掃描方式的擴張－收縮掃描新算法，該算法以擴張－收縮掃描算法為基礎(chǔ)，結(jié)合文獻［5］在逐列掃描過程中進行交叉掃描的思想，改進擴張－收縮掃描算法中逐列掃描的具體計算方法，提高了算法的計算效率，減少了迭代次數(shù)；同時，保留了原擴張收縮掃描算法的所有優(yōu)點.

1 LTI射線追蹤算法的基本方程

如圖1所示，射線通過單元下邊界的AB節(jié)段到達C節(jié)點，射線與AB的交點為D.A點及B點走時分別為TA和TB，節(jié)段AB長為L，單元慢度為s（速度的倒數(shù)），C點距A點的水平及豎向距離均為已知，分別為x，y.建立以A點為原點的局部坐標系，確定A，B，C，D點的局部坐標.現(xiàn)推導C點走時以及交點D 距A點長度r的計算公式［5－6］.

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

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

將式（1）代入式（2）得

根據(jù)費馬原理，TC對r的一階偏導數(shù)應當滿足等于零的條件，即（設(shè)ΔT＝TB－TA）

當L2s2－ΔT2＞0時，解方程（4）可得

① 注：事實上，模型1接收點在按行列分開掃描的方式下，經(jīng)擴張－收縮－再擴張后能夠得到其最小走時，文獻［9］給出的是一次擴張掃描結(jié)果，這應該是文獻［9］作者的疏忽.

若r≥0且r≤L，則

若r＜0或r＞L，則計算r＝0和r＝L時的TC值，并取兩者較小值作為最終TC值.

當L2s2－ΔT2≤0時，TC的計算方法與r＜0或r＞L時的情況一樣.

圖1 經(jīng)過節(jié)段AB到達C點的射線路徑圖Fig.1 The graph of ray path from segment ABto C

在求得r值后，可以通過A點在整體坐標系中的值推出D點在整體坐標系中的值.在這里，定義D點為C點的次級源.

2 擴張收縮掃描算法存在的問題及改進

為說明擴張收縮掃描算法［9－10］存在的問題，以圖2模型為例，模型尺寸為3m×5m，單元邊長及速度分別為1m和3 500m／s；單元邊界劃分為2個節(jié)段；發(fā)射點S位于2號單元下邊界的中點，接收點為a；線段為a節(jié)點的理論射線路徑，其與模型各單元邊界的交點分別為d，c和b.

擴張－收縮掃描算法的第一次擴張掃描過程如圖2所示（假定發(fā)射點所在列各節(jié)點的走時計算已經(jīng)完成），其中圖2（a）為擴張掃描的列掃描結(jié)果，圖2（b）為完成所有列掃描后再進行行掃描的結(jié)果.從圖2（b）可以看出a節(jié)點沒有得到其最小走時.

對于a節(jié)點，此次擴張掃描為無效掃描.同時，從圖2（b）還可以看出，除a節(jié)點外還有其它節(jié)點也進行了無效掃描.易知，當模型網(wǎng)格劃分得更多時，將有更多節(jié)點經(jīng)歷無效掃描.

為了減少無效掃描，提高算法的計算效率.本文采用行列交叉掃描的方式對算法［9－10］進行改進，即在擴張（收縮）掃描過程，行列掃描不再分開進行，而是在其逐列（行）掃描時就同時進行行和列的掃描.仍以圖2所示模型為例，采用改進后的算法對其進行一次擴張掃描，同樣假定發(fā)射點所在列各節(jié)點的走時計算已經(jīng)完成.

圖2 擴張收縮掃描算法的擴張掃描過程及結(jié)果Fig.2 Process and results of expansion scan of expansion－contraction scanning algorithm

圖3（a）為改進算法在擴張掃描時對第3列進行交叉掃描后的結(jié)果，可看出，b和c節(jié)點在這次行列交叉掃描過程中都得到了最小走時，進而在第4列進行交叉掃描后，a節(jié)點能得到其最小走時（圖3（b））.

圖3 改進算法的擴張掃描過程及結(jié)果Fig.3 Scanning process and results of improved algorithm

對比圖2（b）和圖3（b），可以看出，改進算法進行的無效掃描更少，具有更高的計算效率，同時，由于改進算法僅改變了掃描順序，兩種算法的單次擴張掃描計算量不變.同時，對比文獻［5］的一次擴張掃描過程，可以看出，本文提出的交叉掃描算法僅保留了逐列掃描過程中的行掃描，去掉了逐列掃描后附加的逐行掃描過程，并對文獻［5］的逐列掃描過程進行了簡化.

需要說明的是，本文的交叉掃描與擴張收縮掃描在概念上是不同的，前者針對的是一次擴張收縮掃描中行列的計算次序，而后者針對的是整個單元網(wǎng)格的計算次序，一個為微觀操作，另一個為宏觀操作.下文中，除“擴張掃描”和“收縮掃描”中的“掃描”為宏觀操作外，其余“掃描”均指微觀操作.

3 改進算法的基本步驟

改進算法仍分兩步進行，第一步向前計算走時，第二步，向后追蹤射線路徑，其中第二步與文獻［9－10］相同，不贅述.

改進算法的基本步驟如下：

1）如圖4所示，模型為3×3的網(wǎng)格，發(fā)射點S位于第2行第2列.以發(fā)射點所在的行和列為軸，將網(wǎng)格劃分為四個象限，其中發(fā)射點所在行和列為各象限共有部分.圖4（a）和圖4（b）中的陰影部分分別為模型的第一和第二象限.

圖4 模型網(wǎng)格及象限劃分示意圖Fig.4 Model grid and quadrant division schematic diagram

2）計算發(fā)射點S所在單元邊界上各節(jié)點的走時，并記錄次級源.假定此模型在各單元的邊界上均只劃分兩個節(jié)段（圖5），則各單元均有8個節(jié)點.根據(jù)單元節(jié)點以及發(fā)射點在整體坐標系中的坐標，可求得發(fā)射點S所在單元邊界上各節(jié)點的走時，并將S點記為各節(jié)點的次級源.

3）計算發(fā)射點S所在列所有節(jié)點的走時并記錄次級源.

首先計算圖6中與發(fā)射點單元上邊界相連的EJLG單元各節(jié)點的走時.

以I點的走時計算為例（僅考慮通過下邊界GE到達I節(jié)點的射線），易知，滿足I點最小走時要求的射線可能來自GE中的任一節(jié)段.此時，根據(jù)第1節(jié)給出的計算公式和計算方法分別計算出射線通過GF節(jié)段和FE節(jié)段時I點的最小走時，取兩個最小走時的較小值作為I點的最小走時，并記錄相應的次級源.以同樣步驟求出EJLG單元其它節(jié)點的走時，完成該單元的計算.這種通過單元下邊界節(jié)點走時計算其它節(jié)點走時的過程，稱為向上掃描.其它通過上邊界、左邊界、右邊界的情況分別稱為向下掃描、向右掃描以及向左掃描.

圖5 改進算法的基本步驟2Fig.5 The basic steps 2of improved algorithm

圖6 改進算法的基本步驟3Fig.6 The basic steps 3of improved algorithm

然后，以同樣的方式向上逐個計算其它單元的節(jié)點走時，并記錄相應的次級源.

最后，從發(fā)射點單元開始，向下掃描，最終完成發(fā)射點S所在列所有節(jié)點走時計算以及次級源記錄.

4）擴張掃描.對步驟1劃分的四個象限采用不同的掃描策略，具體如下.

象限Ⅰ：從靠近發(fā)射點所在列的右側(cè)列開始，向右逐列進行向右向上交叉掃描，直至模型的最后一列.以靠近發(fā)射點所在列的右側(cè)列為例說明交叉掃描算法：如圖7所示，先對此列單元進行向右掃描，然后再向上掃描，在掃描過程中，若掃描得到的節(jié)點走時比原節(jié)點走時小，則更新此節(jié)點的走時，并記錄對應的次級源，否則，原節(jié)點走時以及相應次級源保持不變.

象限Ⅱ，Ⅲ，Ⅳ的掃描也均從發(fā)射點所在列開始，向左（象限Ⅱ，Ⅲ）或向右（象限Ⅳ）逐列進行向左向上（象限Ⅱ）、向左向下（象限Ⅲ）以及向右向下（象限Ⅳ）交叉掃描，直至模型的邊界.

圖7 改進算法的基本步驟4Fig.7 The basic steps 4of improved algorithm

5）收縮掃描.四個象限的收縮掃描均從模型的豎向邊界開始，向左（象限Ⅰ，Ⅳ）或向右（象限Ⅱ，Ⅲ）逐列進行向左向下（象限Ⅰ）、向左向上（象限Ⅳ）、向右向下（象限Ⅱ）或向右向上（象限Ⅲ）交叉掃描，直至發(fā)射點所在列.

6）重復步驟4）和5），直到模型中所有節(jié)點的走時均不再改變?yōu)橹?至此得到所有節(jié)點的最小走時以及相應的次級源.

對比文獻［5］的逐列掃描過程，本文省去了對每列都進行水平邊界節(jié)點最小走時的搜索，原因如前言所述，在繞射波以及回波存在的情況下，搜索水平邊界節(jié)點最小走時的意義并不明確.同時，為了確定逐列掃描過程中行掃描的順序，本文參考了均勻介質(zhì)模型各節(jié)點的理論射線路徑，如發(fā)射點右上方單元各節(jié)點的最小走時均來自各自單元的左邊界或下邊界，據(jù)此將模型劃分為四個象限，這也與僅改變文獻［9］行列掃描順序的結(jié)果一致.最后，本文按象限劃分掃描區(qū)域進行逐列掃描的方式較文獻［5］的逐列掃描過程更簡潔，程序編制也更為容易.

4 數(shù)值算例

為了驗證本文提出的改進算法在計算效率上的優(yōu)越性，給出三個數(shù)值算例.考慮到改進算法僅改變了掃描順序，不改變單次擴張收縮掃描計算量，故選取計算收斂所需迭代次數(shù)為對比參數(shù).

算例1 如圖8，模型尺寸為3m×5m，各單元邊長及速度分別為1m和3 500m／s，單元邊界劃分為2個節(jié)段，發(fā)射點S位于模型的下邊界，且距模型左邊界1.5m.

采用擴張－收縮掃描算法［9－10］及本文改進算法進行計算.結(jié)果顯示：擴張－收縮掃描算法需進行4次迭代才收斂，而改進算法只需要2次迭代即收斂.

圖8 算例1計算模型Fig.8 Calculation model of example 1

算例2 將算例1模型細劃為30×50的網(wǎng)格，單元邊長為0.1m.如圖9所示，發(fā)射點位置不變.單元邊界劃分為2個節(jié)段.計算結(jié)果顯示：擴張收縮掃描算法需進行31次迭代，而改進算法仍只需2次迭代.

圖9 算例2計算模型Fig.9 Calculation model of example 2

現(xiàn)將1 500號單元在兩種算法下的掃描過程及結(jié)果列于表1，其中，由于算法收斂的條件，原算法及改進算法的最后一次迭代結(jié)果均與前一次的迭代結(jié)果相同.

表1 1 500號單元在兩種算法下的掃描過程及結(jié)果Tab.1 Scanning process and results of No.1 500under the two algorithms

算例1及算例2的計算結(jié)果表明：對于均質(zhì)模型，改進算法相比于原算法具有更快的收斂速度，且當模型網(wǎng)格尺寸劃分較細時，改進算法在計算效率上的優(yōu)勢更顯著.

算例3 如圖10，模型尺寸為20m×20m，單元邊長為1m；單元速度有兩種，其中白色單元的速度為3 500m／s，黑色單元的速度為100m／s，單元邊界劃分為2個節(jié)段，發(fā)射點S和接收點R分別位于模型的左下角和下邊界距發(fā)射點6m位置.

計算結(jié)果顯示：改進算法只需7次迭代即收斂，而原算法則需17次迭代；同時，接收點R在兩種算法下的計算最小走時均為8.439 6ms，與理論最小走時相差0.19%.

圖10 算例3計算模型Fig.10 Calculation model of example 3

算例3的計算結(jié)果表明：改進算法保留了擴張收縮掃描算法能正確處理射線逆向傳播的優(yōu)點，并且具有更快的收斂速度.

5 結(jié) 語

理論分析以及數(shù)值算例表明，本文提出的基于交叉掃描方式的擴張－收縮掃描算法，不僅具有原擴張－收縮掃描算法的所有優(yōu)點，而且在不增加單次擴張－掃描計算量的前提下，通過改變算法的掃描順序，提高了算法的計算效率，加快了算法的收斂速度.特別當模型網(wǎng)格劃分數(shù)較多時，改進算法在計算效率方面的優(yōu)勢更為顯著.

［1］VIDALE J.Finite－difference calculationoftravel times［J］.Bulletin of the Seismological Society of America，1988，78（6）：2062－2076.

［2］QIN F，LUO Y，OLSEN K B，et al.Finite－differencesolution of the eikonal equation alongexpandingwavefronts［J］.Geophysics，1992，57（3）：478－487.

［3］MOSERT J.Shortest path calculation of seismic rays［J］.Geophysics，1991，56（1）：59－67.

［4］劉洪，孟凡林，李幼銘.計算最小走時和射線路徑的界面網(wǎng)全局方法［J］.地球物理學報，1995，38（6）：823－832.LIU Hong，MENG Fan－lin，LI You－ming.The interface grid method for seeking global minimum travel－time and the correspondent ray path ［J］，Acta Geophysica Sinica，1995，38（6）：823－832.（In Chinese）

［5］ASAKAWA E，KAWANAKA T.Seismic ray tracing using linear traveltimeinterpolation［J］.Geophysical Prospecting，1993，41（1）：99－111.

［6］趙改善，郝守玲，楊爾皓，等.基于旅行時線性插值的地震射線追蹤算法［J］.石油物探，1998，37（2）：14－24.ZHAO Gai－shan，HAO Shou－lin，YANG Er－h(huán)ao，et al.Seismic ray tracing algorithm based on the linear traveltimeinterpolation［J］.Geophysical Prospecting for Petroleum，1998，37（2）：14－24.（In Chinese）

［7］張東，謝寶蓮，楊艷，等.一種改進的線性走時插值射線追蹤算法［J］.地球物理學報，2009，52（1）：200－205.ZHANG Dong，XIE Bao－lian，YANG Yan，et al.A ray tracing method based on improved linear traveltimeinterpolation［J］.Chinese Journal of Geophysics，2009，52（1）：200－205.（In Chinese）

［8］WANG Hao－quan.An improved method of linear travel－time interpolation ray tracing algorithm［J］.Acta Physica Polonica－Series A，2010，118（4）：521.

［9］黃靚，黃政宇.線性插值射線追蹤的改進方法［J］.湘潭大學學報：自然科學版，2002，24（4）：105－108.HUANG Liang，HUANG Zheng－yu.An improved method of linear interpolation ray tracing［J］.Journal of Xiangtan University：Natural Science Edition，2002，24（4）：105－108.（In Chinese）

［10］黃靚.混凝土超聲波層析成像的理論方法和試驗研究［D］.長沙：湖南大學土木工程學院，2008：33－35.HUANG Liang.Methodology and experiment research on concrete ultrasonic computerized tomography［D］.Changsha：College of Civil Engineering，Hunan University，2008：33－35.（In Chinese）