动态网络最短路径射线追踪算法中向后追踪方法的改进*1

卢江波，方　志

(湖南大学 土木工程学院，湖南 长沙　410082)



动态网络最短路径射线追踪算法中向后追踪方法的改进*1

卢江波，方志†

(湖南大学 土木工程学院，湖南 长沙410082)

动态网络最短路径射线追踪算法中的向后追踪方法能够解决线性走时插值算法(LTI)向后追踪过程不稳定的问题，但是其计算效率较低.综合利用节点次级源的位置信息以及波的传播规律，提出了改进方法，排除了动态网络最短路径射线追踪算法向后追踪过程中存在的大量冗余计算.数值算例表明，改进的向后追踪方法具有较高的计算效率，是动态网络最短路径射线追踪算法中向后追踪方法的几倍至几十倍；若将改进后的向后追踪方法应用于动态网络最短路径射线追踪改进算法，则该算法的计算效率将提高一倍左右.

射线追踪；线性走时插值；向后追踪方法；计算效率；初至波射线追踪

射线追踪技术在地震层析成像以及混凝土超声波射线层析成像等领域具有重要作用．目前常用射线追踪方法主要有两点射线追踪算法(包括试射法以及弯曲法)[1-3]、有限差分解程函方程法[4-6]、最短路径法[7-10]以及LTI(LinearTravel-timeInterpolation)射线追踪算法[11-22]等.其中，LTI射线追踪算法因其计算精度较高、计算速度较快且适用于任意复杂的速度介质模型，在地震层析成像等领域得到了广泛应用.但是LTI原算法[11]存在两个问题：在向前计算节点最小走时时，不能正确追踪逆向传播的射线，相关节点不能得到正确的最小走时[16-22]；在向后追踪接收点射线路径时，存在不能正确追踪接收点射线路径的可能，算法的稳定性存在不足.

文献[17-18]将波前扩展方式与LTI算法基本方程相结合，提出了动态网络最短路径射线追踪算法，该算法在向前计算各节点的最小走时时，从震源点开始，采用波前扩展的方式逐点计算各个节点上的最小走时，改变了LTI原算法的计算方式，确保各节点能得到最小走时；在向后追踪接收点射线路径时，基于互换原理考虑了接收点所有可能的射线路径，确保算法能正确追踪接收点的射线路径.

动态网络最短路径射线追踪算法虽然能够解决LTI原算法存在的两个问题，但是其计算效率偏低.文献[22]基于波的传播规律提出了动态网络最短路径射线追踪改进算法，改进并提高了动态网络最短路径射线追踪算法向前计算节点最小走时这一步骤的计算效率，但是，该改进方法的向后追踪方法仍然采用动态网络最短路径射线追踪算法中的方法，计算效率依然偏低.

针对动态网络最短路径射线追踪算法中的向后追踪方法存在计算效率低的问题，本文提出了改进方法.首先，在向前计算节点最小走时这一步骤中，不仅计算各节点的最小走时，而且还记录各节点次级源的位置；然后，在向后追踪接收点的射线路径时，利用各节点次级源的位置信息以及波的传播规律对算法进行改进，降低算法的计算量，提高算法的计算效率.

1　LTI算法基本方程的推导及接收点次级源的选定

如图1所示，射线通过AB节段的D点到达C节点.A,B节点的走时分别为TA和TB，节段AB长为L，单元慢度为s.以A点为原点建立局部坐标系，C,D点的局部坐标分别为(xc,yc)、(r,0).现推导通过AB节段，C节点能得到的最小走时，以及C节点取得最小走时时D点在局部坐标系中坐标r的计算公式[11-12]．

图1　通过节段AB到达C点的计算示意图Fig.1　Schematic of computing from segment AB to C

由LTI算法的线性走时假定可得D点的走时：

(1)

根据D点的走时，以及C,D点在局部坐标系中的坐标和单元慢度，得到C点走时：

(2)

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

(3)

根据费马原理，TC对r的一阶偏导数应当满足等于零的条件，并记Δ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时的情况相同.

文献[8]给出了最短路径射线追踪算法中次级源的确定方法，考虑到LTI算法中的射线可通过单元边界上任意点进行传播，这与最短路径射线追踪算法中射线只能通过节点进行传播不同，因此最短路径射线追踪算法中次级源的确定方法不完全适合LTI算法.为了更好地描述本文提出的改进方法，将LTI算法中次级源的确定方法规定如下：通过接收点所在单元各节段(不包括接收点所在节段)可得到多个接收点走时，其中走时最小值对应的点即为接收点的次级源，如果走时最小值对应多个点，取距离接收点最近的点作为接收点的次级源.如图1所示，若接收点C通过AB节段D点得到的走时小于通过单元其它节段得到的走时，那么D点即为接收点C的次级源.

2　LTI原始算法向后追踪过程存在的问题及动态网络最短路径射线追踪算法存在的不足及改进

LTI原始算法在向后追踪接收点射线路径时，首先逐点计算接收点所在单元中各节点走时与节点至接收点的走时之和，然后选出最小的走时之和以及相应的节点，最后将单元中包含该节点的节段作为接收点次级源的可能区域[11].事实上，接收点的次级源并不一定在这些节段中.此外，LTI原算法在确定接收点次级源的可能区域时，并未排除接收点所在的节段，算法可能会陷入无限循环.因为通过接收点所在节段线性插值得到的接收点走时，可能比通过单元其它节段得到的走时更小，那么接收点取最小走时对应的点可能为接收点本身，算法可能会进入无限循环.

以图2所示模型为例，模型尺寸为3m×3m，单元大小为1m×1m，单元边界划分为2个节段，模型上层、中层以及下层单元的速度分别为525m/s,515m/s和505m/s，震源S的x,y坐标分别为1.5和3．0，接收点R的x,y坐标分别为2.05和0；图中虚线为采用LTI原始算法中的向后追踪方法得到的射线路径，需要说明的是，计算这条射线路径时，已经排除了接收点的次级源位于接收点所在节段的情况，实线为根据动态网络最短路径射线追踪算法中的向后追踪方法[17-18]得到的射线路径.

虚线为LTI算法向后追踪结果，实线为动态网络 最短路径射线追踪结果图2　LTI原始算法向后追踪方法存在的问题Fig.2　Problems of LTI algorithm’s backward tracing method

从图2可以看出，在III号单元中，两种方法的计算结果一样，接收点R的次级源均为R1.在II号单元中，对于LTI原始算法的向后追踪方法，首先须确定II号单元中节点走时与节点至接收点R1走时之和最小的节点，经计算确定为d节点，然后确定cd或de节段为接收点R1次级源的可能区域；对于动态网络最短路径射线追踪算法中的向后追踪方法，接收点R1次级源所在区域为ef节段.计算结果表明：通过ef节段计算得到的R1点的走时(T=0.005 416 05s)要小于通过cd或de节段的走时(T=0.005 422 57s)，接收点R1的次级源不在cd或de节段内.造成这一问题的根源在于:LTI算法假定射线路径可以经过单元边界的任意一点，在向后追踪过程中必须考虑所有可能的射线路径，并根据费马原理选择走时最短的那条路径[11]，而LTI原始算法的具体追踪方法却没有考虑接收点所在单元所有节段的射线，而是采用了一种“简化”方法确定接收点射线路径或次级源的可能区域，由图2所示模型可以看出，这种“简化”方法存在不足.

此外，在确定接收点次级源时，若不排除接收点所在的节段，那么在计算新接收点R1的次级源时，由于通过节段cd插值得到的走时，比通过ed节段得到的走时小，计算得到的“次级源”为R1本身，程序将陷入无限循环.

动态网络最短路径射线追踪算法中基于互换原理提出的向后追踪方法[17-18]，考虑了来自接收点所在单元中除接收点所在节段外所有节段的射线，确保了接收点能得到其次级源.但是，该方法计算量大，计算效率低.为了解决这一问题，本文首先在向前计算节点最小走时的步骤中，建立了一个数组，专门用于记录节点次级源的位置信息.然后利用各节点次级源的位置信息以及波的传播规律对动态网络最短路径射线追踪算法中的向后追踪方法进行改进.现以图3所示模型为例，对该向后追踪方法的基本步骤及计算策略进行说明，同时对向后追踪改进方法进行阐述，改进前后的向后追踪方法及计算策略分别如图3和图4所示.具体步骤和分析如下：

1)首先将接收点分为以下3种情况，然后对不同的情况采用不同的策略求解接收点的次级源：

①接收点位于单元内部，如图3(a)及图4(a)所示的R点.此时利用LTI算法的相关公式计算接收点所在单元各节段至接收点的走时，然后从中选出走时最小值对应的点，如图3(b)及图4(b)中的R1点，这个点就是接收点的次级源.

②接收点位于单元边界上，但是非单元边界上的节点.如图3(b),(d)中的R1和R3点，此时除接收点所在的节段外，单元中的其它节段均要计算从该节段至接收点的最小走时，然后从中选出走时最小值对应的点，如图3(c),(e)中的R2和S点.事实上，对于接收点位于单元边界但非单元节点的情况，除非该接收点为原始接收点(图3和图4中的R)，否则并不需要对接收点进行全方位的计算.可以采取如下两个步骤确定接收点次级源的可能区域．

图3　动态网络最短路径射线追踪算法向后追踪方法的基本步骤和计算策略Fig.3　Basic steps and computing strategy of the dynamic networks tracing algorithm’s backward tracing method

首先，利用接收点所在节段端点的次级源位置信息确定两个定位点．定位点的确定方法为：若节点端点的次级源不在单元内，或者次级源在单元内但与节点端点处于同一边界，则定位点为节段端点本身，其它情况定位点为节点端点的次级源.如图4(b)所示，接收点R1所在节段的端点为A1和A2，易知接收点R1的次级源位于单元I中，现以端点A1的次级源位置信息为例说明定位点的确定方法，若A1的次级源不在单元I中，比如A1的次级源为单元II中的b点，则由A1的次级源位置信息确定的定位点为A1本身；若A1的次级源与A1处于同一边界，比如次级源位于A1A4或者A1A6边界，则定位点为A1本身；对于其它情况，比如A1的次级源为单元I中的c点，则定位点为A1的次级源c.在图4(b)中，假定A1节点的次级源为A3节点，A2节点的次级源为a1点，则由A1，A2节点的次级源位置信息可以确定两个定位点，分别为A1节点本身以及点a1.

同理得到图4(d)中R3点的次级源S.

③接收点位于单元边界上，且为单元边界上的节点.如图3(c)中的R2节点，此时需要在R2节点所处的几个单元中执行情况②的计算，如图3(c).

显然，由于在向前计算节点最小走时的过程中，已经记录了各节点的次级源，因此，在改进方法中，可以直接得到接收点的次级源，如图4(c).

2)以步骤1)中获得的次级源为新的接收点，重复步骤1)，直至新的接收点为震源.

3)依次连接各接收点得到R点的初至波射线路径，如图3(e)和图4(e)所示.

对比图3与图4可以看出，改进后的向后追踪方法计算量明显减少.

图4　向后追踪改进方法的基本步骤和计算策略Fig.4　Basic steps and computing strategy of the backward tracing method which presented in this study

3　数值算例

为了对比改进前后向后追踪方法的计算效率，建立了尺寸为2 500m×600m的二维模型，如图5所示，其中x值在500～2 000m之间且y值在200～400m之间的区域为低速区，速度为500m/s，其余速度均为4 000m/s，震源S位于模型上表面的正中间(1 250，0)，单元尺寸为5m×5m，单元边界划分段数为4,10,20等3种情况，模型下表面的每个单元布置一个接收点，共500个，部分接收点的射线路径如图6所示，分别记录两种向后追踪方法追踪所有接收点的射线路径所耗费的总时间.此外，为了说明向后追踪方法计算效率的提高对整个算法的影响，将改进后的向后追踪方法应用于动态网络最短路径射线追踪改进算法[22]，然后比较应用改进方法前后动态网络最短路径射线追踪改进算法的总耗时(不包括算法的前处理过程).计算机CPU主频为3.4GHz，计算结果如表1所示.

x/m图5　速度模型Fig.5　 Velocity model表1　计算效率对比Tab.1　Computational efficiency comparisons

注：算法A为动态网线最短路径射线追踪改进算法；算法B为动态网线最短路径射线追踪改进算法与本文提出的向后追踪改进方法相结合的一种算法.

由表1可知，向后追踪改进方法的计算效率较高，是改进前向后追踪方法的几倍至几十倍.并且能将动态网络最短路径射线追踪改进算法的计算效率提高1倍左右.

x/m(a)动态网络最短路径射线追踪算法的向后追踪结果

x/m(b)改进后的向后追踪结果图6　部分接收点的射线追踪结果Fig.6　Ray path tracing results of a part of receive points

为验证向后追踪改进方法对复杂模型的有效性，采用向后追踪改进方法对Marmousi速度模型进行射线追踪，其中节点最小走时的计算采用的是动态网络最短路径射线追踪改进算法.模型尺寸为9 192m×2 904m，震源S位于模型左上角(0，0)，单元尺寸为24m×24m，单元边界划分为10段，共设置5个接收点R1～R5，分别位于模型上表面的3 600,4 800m,6 000m,7 200m和8 400m，射线追踪结果如图7(a)所示，作为对照，本文给出了单元边界划分为30段时最短路径法的射线追踪结果，如图7(b)所示.两种算法下，各接收点根据射线追踪结果计算的走时如表2所示.计算结果表明，向后追踪改进方法对于复杂的速度模型同样有效.

x/m(a)向后追踪改进方法的追踪结果

x/m(b)最短路径法的追踪结果图7　Marmousi速度模型的射线追踪结果Fig.7　Ray path tracing results of theMarmousi velocity model

表2　各接收点的最小走时Tab.2　The minimum travel time of the receive points

注：改进算法为动态网线最短路径射线追踪改进算法与本文提出的向后追踪改进方法相结合的一种算法．

4　结　论

针对LTI原算法中向后追踪方法存在的无限循环以及可能不能得到正确射线路径的问题，动态网络最短路径射线追踪算法中基于互换原理提出的向后追踪方法能够予以有效的解决.但是该算法存在较多的无效计算.本文根据模型中节点次级源的位置信息以及波的传播规律，提出了改进的向后追踪方法．数值结果表明，改进后的向后追踪方法，其计算效率较之改进前的方法有较大程度的提高；此外，本文提出的向后追踪改进方法能将动态网络最短路径射线追踪改进算法的计算效率提高1倍左右.

[1]JULIANBR，GUBBINSD.Three-dimensionalseismicraytracing[J].JGeophys,1977,43(1/2):95-113.

[2]田玥，陈晓非.水平层状介质中的快速两点间射线追踪方法[J]. 地震学报,2005,27(2):147-154.

TIANYue，CHENXiao-fei.Arapidandaccuratetwo-pointraytracingmethodinhorizontallylayeredvelocitymodel[J].ActaSeismologicaSinica,2005,27(2):147-154. (InChinese)

[3]XUTao，ZHANGZhong-jie，GAOErgen，et al.Segmentallyiterativeraytracingincomplex2Dand3Dheteroge-neousblockmodels[J].BulletinoftheSeismologicalSocietyofAmerica, 2010,100(2): 841-850.

[4]VIDALEJ.Finite-differencecalculationoftraveltimes[J].BulletinoftheSeismologicalSocietyofAmerica, 1988,78(6): 2062-2076.

[5]QINFu-hao,LUOYi,OLSENlKB, et al.Finite-differencesolutionoftheeikonalequationalongexpandingwavefronts[J].Geophysics, 1992,57(3):478-487.

[6]李振春，刘玉莲，张建磊,等.基于矩形网格的有限差分走时计算方法[J]. 地震学报, 2004,26(6):644-650.

LIZhen-chun,LIUYu-lian,ZHANGJian-lei,et al.Finite-differencecalculationoftraveltimesbasedonrectangulargrid[J].ActaSeismologicaSinica, 2004,26(6):644-650. (InChinese)

[7]MOSERTJ.Shortestpathcalculationofseismicrays[J].Geophysics, 1991,56(1): 59-67.

[8]刘洪，孟凡林，李幼铭.计算最小走时和射线路径的界面网全局方法[J].地球物理学报, 1995,38(6): 823-832．

LIUHong，MENGFang-lin，LIYou-ming.Theinterfacegridmethodforseekingglobalminimuntraveltimeandcorrespondentraypath[J]ActaGeophysicaSinicaChinJGeophys-ch,1995,38(6):823-832.(InChinese)

[9]赵爱华，徐涛.提高规则网格最短路径方法反射波走时计算精度的走时校正技术[J]. 地球物理学进展, 2012,27(5):1854-1862.

ZHAOAi-hua，XUTao.Atraveltimecorrectiontechniqueforimprovingtheaccuracyofreflectionwavetraveltimeswiththeshortestpathmethodbasedonaregulargrid[J].ProgressinGeophys, 2012,27(5):1854-1862.(InChinese)

[10]BAIChao-yin，HUANGGuo-jiao，ZHAORui. 2-D/3-Dirregularshortest-pathraytracingformultiplearrivalsanditsapplications[J].GeophysicalJournalInternational, 2010,183(3): 1596-1612.

[11]ASAKAWAE，KAWANAKAT.Seismicraytracingusinglineartraveltimeinterpolation[J].GeophysicalProspecting, 1993,41(1): 99-111.

[12]赵改善，郝守玲，杨尔皓,等.基于旅行时线性插值的地震射线追踪算法[J].石油物探, 1998,37(2):14-24.

ZHAOGai-shan,HAOShou-ling,YANGEr-hao,et al.Seismicraytracingalgorithmbasedonthelineartraveltimeinterpolation[J].GeophysicalProspectingforPetroleum，1998,37(2):14-24.(InChinese)

[13]CARDARELLIE，CERRETOA.Raytracinginellipticalanisotropicmediausingthelineartraveltimeinterpolation(LTI)methodappliedtotraveltimeseismictomography[J].GeophysicalProspecting, 2002,50(1): 55-72.

[14]聂建新，杨慧珠.地震波旅行时二次/线性联合插值法[J].清华大学学报:自然科学版 , 2003,43(11):1495-1498.

NIEJian-xin,YANGHui-zhu.Quadratic/lineartraveltimeinterpolationofseismicraytracing[J].JTsinghuaUniv:Sci&Tech, 2003,43(11):1495-1498.(InChinese)

[15]ZHANGJian-zhong,HUANGYue-qin,SONGLin-ping,et al.Fastandaccurate3-Draytracingusingbilineartraveltimeinterpolationandthewavefrontgroupmarching[J].GeophysicalJournalInternational, 2011,184(3): 1327-1340.

[16]黄靓，黄政宇.线性插值射线追踪的改进方法[J]. 湘潭大学自然科学学报, 2002,24(4): 105-108.

HUANGLiang，HUANGZheng-yu.Animprovedmethodoflinearinterpolationraytracing[J].NaturalScienceJournalofXiangtanUniversity, 2002,24(4): 105-108. (InChinese)

[17]张建中，陈世军，余大祥.最短路径射线追踪方法及其改进[J]．地球物理学进展, 2003,18(1):146-150．

ZHANGJian-zhong,CHENShi-jun,YUDa-xiang.Improvementofshortestpathraytracingmethod[J].ProgressinGeophys, 2003,18(1):146-150.(InChinese)

[18]张建中，陈世军，徐初伟.动态网络最短路径射线追踪[J]. 地球物理学报, 2004,47(5): 899-904．

ZHANGJian-zhong，CHENShi-jun，XUChu-wei.Amethodofshortestpathraytracingwithdynamicnetworks[J].ChineseJ.Geophys， 2004, 47(5): 899-904.(InChinese)

[19]黄靓.混凝土超声波层析成像的理论方法和试验研究[D]. 长沙: 湖南大学土木工程学院, 2008：33-36.

HUANGLiang.Methodologyandexperimentresearchonconcreteultrasoniccomputerizedtomography[D].Changsha:CollegeofCivilEngineering,HunanUniversity, 2008:33-36(InChinese)

[20]张东，谢宝莲，杨艳,等.一种改进的线性走时插值射线追踪算法[J]. 地球物理学报, 2009,52(1): 200-205.

ZHANGDong,XIEBao-lian,YANGYan,et al.Araytracingmethodbasedonimprovedlineartraveltimeinterpolation[J].ChineseJ.Geophys, 2009,52(1): 200-205.(InChinese)

[21]卢江波，方志.线性走时插值射线追踪算法的改进[J].湖南大学学报：自然科学版, 2014,41(1):39-44．

LUJiang-bo,FANGZhi.Improvementofraytracingalgorithmbasedonlineartraveltimeinterpolation[J].JournalofHunanUniversity:NaturalSciences, 2014,41(1):39-44. (InChinese)

[22]卢江波，方志.一种线性走时插值射线追踪改进算法[J].地震学报, 2014,36(6):1089-1100.

LUJiang-bo,FANGZhi.Animprovedray-tracingalgorithmbasedonlineartravel-timeinterpolation[J].ActaSeismologicaSinica,2014,36(6):1089-1100.(InChinese)

An Improved Method on Backward Tracing of the Shortest Path Raytracing Algorithm with Dynamic Networks

LU Jiang-bo，FANG Zhi†

(College of Civil Engineering，Hunan Univ，Changsha, Hunan410082，China)

ThebackwardtracingmethodoftheshortestpathraytracingalgorithmwithdynamicnetworkscansolvetheunstabilityprobleminthebackwardtracingprocedureoftheLTI(LinearTravel-timeInterpolation)algorithm,butthecomputationalefficiencyofthemethodislow.Thisstudypresentedanimprovedmethodonbackwardtracing.Accordingtothelocationinformationofthesecondarysourcesforthenodesandthelawofwavepropagation,alargenumberofredundancycalculationareexcludedinthebackwardtracingofthedynamicnetworkstracingalgorithm.Thenumericalexamplesshowthattheimprovedmethodexhibitsthehighercomputationalefficiency.Thecalculationefficiencyoftheimprovedmethodisseveraltimesthatofthebackwardtracingmethodofthedynamicnetworkstracingalgorithm.Whentheimprovedmethodisappliedtotheimprovedalgorithmoftheshortestpathraytracingwithdynamicnetworks,thecomputationalefficiencyofthealgorithmcanbeincreasedbyabout100 %.

raytracing;lineartraveltimeinterpolation;improvedalgorithm;backwardtracing;computationalefficiency；firstarrivalraytracing

1674-2974(2016)05-0106-07

2015-06-17

国家自然科学基金资助项目(51278182), National Natural Science Foundation of China(51278182)

卢江波(1987-)，男，湖南郴州人，湖南大学博士研究生

†通讯联系人，E-mail: fangzhi@hnu.edu.cn

P631

A