最短路径射线追踪方法及其改进

井间地震层析成像

第18卷　第1期　　　　　　　　　　　　地　球　物　理　学　进　展　　　　　　　　　　　　Vol.18　No.12003年3月(146～150)　　　　　　　　　　PROGRESS　IN　GEOPHYSICS　　　　

　　　　　　March　2003

最短路径射线追踪方法及其改进

张建中1,3,　陈世军2,3,　余大祥3

(1.厦门大学电子工程系,厦门,361005;　2.中国科学院地质与地球物理研究所,北京,100029;　3.胜利油田,东营257000)

摘　要　综述了用网络最短路径算法求解地震射线追踪问题的原理、方法技术以及存在问题和改进措施.特别介绍了作者在最短路径算法基础上,提出的动态网络最短路径地震射线追踪方法.该方法先采集从炮点到整个模型所有节点上的初至旅行时,其中,在一个单元内,对相邻每对已计算出最小旅行时的节点进行线性插值,并利用Fermat原理计算未知节点的最小旅行时;然后,利用同样的方法,从接收点开始,反向追踪炮点到接收点的射线路径.该方法能适于各种复杂的非均匀介质,极大地提高了射线追踪的精度.关键词　最短路径算法,射线追踪,动态网络

中图分类号　P315　　　　　　　　文献标识码　A　　　　　　　文章编号　100422903(2003)0120146205

Improvementofshortestpathraytracingmethod

ZHANGJian2zhong1,　CHENShi22,　Da2xiang3

(1.CollegeofComputerandInformationEngineering,Xiamen,,;.ofandGeophysics,ChineseAcademyofSciences,Beijing,;3.Field,,,)

Abstract　Thealgorithminnetworksarediscussed,therelevantdraw2backsinwhichareEspecially,wedevelopashortestpathraytracingmethodwithdynamicnetworksbasedonpathalgorithm.Thewavefronttimesaresampledatthenodesawayfromthesourcethroughouttheen2tiremodel.Thefirstarrivaltimesatnodesinacellareexpressedwithlinearinterpolationbetweentraveltimesobtainedaforehandateachpairoftwoneighboringnodes,andthendeterminedbyFermatprinciple.Afterthefirstarrivaltimesbeingsampledthroughoutthemodel,theraypathsfromthesourcetoeachreceiverareobtainedbackwardawayfromthereceivertosourceinthesamewayastheforwardsampling.Thisalgorithmisregardlessofmodelcomplexity,andhasahigherprecisionforcomputationofseismictravel2timesandraypaths.

Keywords　shortestpathalgorithm,raytracing,dynamicnetworks

0　引　　言

不同领域的很多科技问题都可归结为网络理论和网络中的最短路径问题.交通图就是一个典型的例子,其网络由城镇和连接这些城镇的道路构成.人们常常寻求从一个城镇到另一个城镇或环绕所有城镇的最短路径.Fermat原理指出,地震波沿着旅行时最小的路径传播.地震射线追踪的目的就是寻找震源到各个接收点的地震波传播的最小旅行时路径,这对于研究地震波的传播规律以及地震层析成像具有十分重要的意义.波场的数值计算需要对地下复杂介质进行离散化处理,即把介质划分成若干小单元,这样单元边界上的节点以及彼此相邻节点的连

收稿日期　2002207224;　修回日期　2002210210.

线就构成了一个网络.若把震源点和接收点也作为

网络的节点,并与其相邻的其它节点相连,我们就可以借鉴网络中的最短路径算法求出从震源点到包括接收点在内的其它所有节点的最小旅行时路径.由于网络理论没有空间维数的概念,二维和三维射线追踪的最短路径算法是一样的.

最短路径方法起源于网络理论[1],首次由Nakanishi和Yamaguchi[2]应用于地震射线追踪中.Moser[3]以及Klimes和Kvasnicha[4]对最短路径射线

追踪方法及其性能进行了详细研究.最短路径射线追踪方法是用网络节点之间的最小旅行时连线近似地震射线路径的.当网络较稀时,计算的射线路径及其旅行时有较大的误差,且误差的大小与网络节点

作者简介　张建中,男,1963年生,陕西白水人,博士,副教授,现在厦门大学计算机与信息工程学院从事信号处理和地球物理方面的科研工作.

曾获省级科技进步一等奖,在国内外期刊发表学术论文30余篇(Email:zhangjianzhong98@jzzhang@).

井间地震层析成像

1471期　　　　　　　　　　　　张建中,等:最短路径射线追踪方法及其改进

的分布和数量密切相关.Klimes和Kvasnicha[4]及

Zhang和Toksoz[5]等通3过改变节点的分布,Moser等[6]、Fischer和Lees[7]及VanAvendonk等[8]通过引入其它方法进一步修正射线路径,来提高计算效率,减小计算误差.为了克服现有方法中波只能从一个节点到达另一个节点的限制,我们建立了动态网络,并在一个单元内,利用节点最小旅行时插值和Fer2mat原理来求任意节点的局部最小旅行时及其子震源位置,提高了射线追踪的精度.

之间的旅行时为它们之间的欧氏距离与其平均慢度之积.图1所示(二维),一个节点与其周围16个节点直接相连,这样,一个震源或子震源就可以有16个出射方向和16个与之对应的计算节点;若是三维情形,网络节点分布在立方体单元的角点上,震源则需按26个方向出射.这种设置便于绘制旅行时和速度的等值线图,更重要的是可以引入速度界面,进行反射波的射线追踪[3].在实际应用中,Moser等[6]以及Klimes和Kvasnicha[4]均采用了这种模型

.

1　最短路径射线追踪方法原理

最短路径射线追踪方法基于地震波传播的Fer2mat最小旅行时原理和网络理论中的最短路径算法来实现.把地下介质离散成若干小单元体,并在各单元边界上设置一些节点,地下速度模型就可表示成由这些节点以及它们之间的连线所形成的网络.网络中的每一个节点只能与彼此相邻的节点连接.相邻节点之间的连接权等于地震波沿其连线传播的旅行时.一条路径是由相互连接的节点序列组成的,.原理,.1.1　速度模型和网络的建立

在最短路径射线追踪算法中,介质一般都被离散化成规则的矩形单元,网络结构和与之对应的速度模型却有下列不同的类型

.

图2　模型Ⅱ节点设置及出射射线分布

点线是单元边界;空心点和点划线表示按等间距设置的节点和相应的出射射线(a1=a3=0.25d,a2=0.50d);实心点和实线表示按等出射角设置的节点和相应的出射射线(b=b3=0.29289d,b2=0.41421d)

Fig.2　NodesandraysdistributionformodelⅡ

图1　模型Ⅰ速度和网络结构

速度离散值分布在各节点上,任一节点与其周围的16个节点直接相连形成网络.

Fig.1　VelocityandnetworksstructureformodelⅠ

模型Ⅰ　网络节点设置在矩形单元的四个角点上,整个介质的速度或慢度用这些角点上的离散采样值表示,每个节点都可以与相邻的节点相连,相邻节点

模型Ⅱ　网络节点按一定规则分布在矩形单元的边

界上(不含单元角点),在每一个矩形单元内,速度或慢度为常数.只有在两个节点之间不存在单元边界时,这两个节点才可以连接,如图2所示.相连节点之间的旅行时由它们之间的欧氏距离乘以它们所在单元的慢度得到.这种模型由常速度单元组成,使两个相连节点之间的旅行时的计算简单而快速,并且在一个单元内的射线追踪是精确的.Moser[3]及Fis2cher和Lees[7]按照等间距规则把节点均匀分布在网格上,每个单元边界上设置二个节点的节点分布如图2中的空心点所示,若单元边界的长度(即网格距)为d,那么节点与单元角点及节点之间的距离则为a1=a3=0.25d,a2=0.5d.Zhang和Toksoz[5]按照等角原则(即震源与节点连线之间的夹角相等)设

井间地震层析成像

置节点,如图2中的实心点所示,若单元边界的长度为d,这时节点与单元角点及节点之间的距离分别为b1=0.29289d,b2=0.41421d,b3=0.29289d.对于折射波的计算,节点的等角分布要比等间隔分布减小50%的误差.

1.2　节点最小旅行时的计算和射线路径的确定

2　存在问题及一些改进方法

最短路径射线追踪方法是用离散的网络节点之

间的最小旅行时连线近似地震射线路径的.只有当网络节点很密时,才能较好地逼近实际波传播路径,这势必需要大量的计算时间和计算机内存,特别是三维问题的计算量令人难以接受.但是当节点较稀时,射线常常呈“之”字形路径,计算出的旅行时将比实际旅行时系统偏大,且在波传播方向上节点越少,误差越大.另一个问题是,当网络稀疏时,特别对于速度变化平缓的区域,在两个点之间常常会有几条等时最短路径.其中可能有一条能较好地近似真实射线路径,而其它的却不能.Dijkstra算法[1]不能区别这些路径,选取那一条具有任意性,有时不能得到接近真实射线的路径.所以说,这种射线追踪问题具有一定的不确定性,在计算过程中,需要施加一些约束,,.2.1,Snell定律对获得的射线路径进行修正[7].对射线上三个相邻节点,固定两端点,用Snell定律改变中间点的位置,使该小段射线的旅行时达到最小.沿射线对每三个相邻节点的射线段进行优化,不断迭代,直到整个射线的旅行时几乎不变时(小于0.01%)为至.这些方法能使射线更圆滑,改善了射线的旅行时,但由于不能使射线的移动超出一个网格单元,从而不能在较大程度上改善射线空间位置的误差.

对于速度变化很小的模型,上述第二个问题是很常见的,为了解决该问题,人们必须寻求一种方法,从两点之间存在的几条相同旅行时路径中,确定最接近真实射线的路径.传统的做法是在单元边界上增加节点数,但实际上随着节点的增加,该问题会更严重.Fischer和Lees在Dijkstra算法中引入直射线追踪的约束[7].当遇到两点之间出现几条相同旅行时路径时,把它们与直线路径的旅行时相比较,如果直线路径的旅行时小,就用直线路径取代原路径,再用Snell定律优化方法微调直线路径.2.2　混合射线追踪方法

Moser等[6]把Dijkstra算法得到的射线路径作为初值,再用经典的两点射线追踪方法—弯曲法对每条射线进行改善.这种混合射线追踪方法被Papaza2chos和Nolet[12]及VanAvendonk等[13]成功地用于三维层析反演中.VanAvendonk等[8]对这种混合方法

给出震源点的旅行时,从震源点开始,用Dijk2stra算法[1]逐点算出从震源到达各个网络节点上的

最小旅行时.在整个过程中,把所有网络节点的集合

N分成四个子集.P:已经获得最小旅行时的节点,

即已作过子波源的节点集合;Q:已计算出从P中至少一个节点波传来的旅行时,但还没有作子震源的节点集合;R:在N中除去P、Q后剩余节点的集合;F(i):与震源或子震源直接相连的节点集合,每一步仅计算这些节点的旅行时.具体算法为:①初始化.P:=Φ,t(s)=0,s为震源,Φ表示空集;Q:=

N,t(i)=∞,i∈N.②选择.在Q中选择旅行时最

小的节点i,i∈Q.③更替.计算从i点传到j点的旅行时,若该值比原值小,,保持原值不变.即,Ft=(j),t(i)+dijF(,t(j=t(i)+dij.

将节点i从QP.④迭代判断.如果P=N或Q=Φ停止迭代,否则转向②.确定震源点到某一点的波传播射线路径,是从该点开始向震源逆向搜寻处理的过程.在从震源点开始扩展波前的过程中,利用指针数组记录了各个最小旅行时节点前一级子震源的位置.确定射线路径时,从接收点开始,依次找出前一级子震源,直到震源点为止,把震源点和找出的各子震源节点以及接收点顺次相连,就得到相应的射线路径.

1.3　最小旅行时节点的确定

在上述算法的迭代过程中,每一步都要从Q中选择出旅行时为最小的点作为下一个子震源,这是最耗时的一步.在Dijkstra算法[1]中,所用比较方法的计算时间为O(n2)量级(n为节点总数).如果把所有旅行时从大到小完全排队,仍然不能改善这种量级的时间消耗.Johnson[9]、Gallo和Pallotlino[10]采用堆排序法,将所需时间减少至O(nlogn)量级.常旭和刘伊克[11]通过对各种排序算法进行的对比研究认为,对于无序数列,平均性能最好的算法是快速排序算法,而对于正序数列,插入排序算法所耗时间可以减到O(n)量级.因此,综合使用快速排序算法和插入算法,能提高寻找最小旅行时节点的速度.

井间地震层析成像

进行了详细研究.混合方法能在相对低的计算成本

上获得较精确的射线路径及其旅行时,但下列条件必须成立:首先,最短路径方法得到的解必须尽量接近实际射线路径,为射线弯曲方法提供一个正确的可行解;其次,弯曲射线的迭代过程必须收敛到稳定的解.前者要求较高精确的最短路径解,后者要求射线弯曲要有足够的迭代时间.

在Dijkstra算法[1]中,随着波前面的不断扩展,要不断由已知点的旅行时计算到达其相邻节点的最小旅行时.在一个单元内,我们利用旅行时插值和Fermat原理来求任意节点的局部最小旅行时及其子震源位置.

如图4所示,在慢度为s的矩形单元内,如果A、B为已知旅行时的点,其座标分别为(xa,za)、

(xb,zb),旅行时分

3　动态网络最短路径方法

上述最短路径算法的网络是始终固定不变的,这使得地震波只能从一个节点到达另一个节点,任意两点之间的射线路径也只能由其间的有关节点直接连接而成,各节上射线的出射角不能随入射角连续变化.为了克服这些问题,我们提出了动态网络最短路径射线追踪方法

.

别为t(a)和t(b),现在来求从A、B之

间某一点通过到达D图4　

由二个节点旅行时求另一

点的最小旅行时和射节点旅行时及射线的几何关系线路径.设射线从A、Fig.4　GofgettingminimalB之间的,DfromnodesA记(cd,)t(d),且有xc=xa.

假设在A、B之间旅行时线性变化,那么可以得到由A、B点的旅行时求C点旅行时的线性插值公

图3　节点设置及可能接收的穿过单元边界的最小旅行时射线

(a)接收点不在单元角点上;(b)接收点在单元角点上.

虚线为单元边界,实心点表示节点,实线表示可能从节点及两节点之间任意一点通过到达接收点的最短旅行时路径.

式,进而得到D点的旅行时表达式,其未知量是C点的坐标.根据Fermat原理,对该式求极小,就可得到C点的位置和经过C点到达D点的最小旅行时分别为

zc=zd-(xd-xa)[t(b)-t(a)]

Fig.3　Nodessettingandpossibleminimaltraveltimerayscomefromcelledges

×{s2(zb-za)2-[t(b)-t(a)]2}

t(d)=t(a)+

-1/2

,(1)

[t(b)-t(a)]+

zb-za

(2)

3.1　节点分布和动态网络

介质被离散成矩形单元,每个单元的速度为常数,节点等间距分布在单元边界(包括单元角点)上,如图3所示.这里设置的节点指的是需要计算最小旅行时的点,不一定是网络节点,而网络节点则是在波前传播过程中即时确定的.在一个单元内,我们不只是把已知最小旅行时的节点作为子震源,求出未知节点上的旅行时,而是利用Fermat原理,在已知最小旅行时的节点所在边界上,求出以最短旅行时路径传至某个节点的子震源位置(即网络节点)及其最小旅行时.这样求出的射线路径不是单元边界上离散节点的连线,而是穿过单元边界上正好满足最小旅行时条件的那一点与到达点的连线.3.2　节点旅行时的求取

2221/2

{s(zb-za)-[t(b)-t(a)]}.

zb-za

由于射线从A、B之间通过,必须满足约束条件

()2()222

22≤[t(b)-t(a)]/s(zd-zb)+(xd-xa)(z-z)2(z-z)2

≤.(zd-za)2+(xd-xa)2

(3)

该式是在最小旅行时条件下,射线穿过A、B之间的正确位置到达D点的约束条件.只有在此条件满足的情况下,才能用(1)和(2)式分别计算出射线通过的准确位置C的坐标和到达D点的最小旅行时.3.3　射线路径的确定

在本方法中,由于射线路径可以从单元边界上的任何点(不限于边界上的已知旅行时节点)通过,

井间地震层析成像

150　　　　　　　　　　　　　　地　球　物　理　学　进　展　　　　　　　　　　　　　　　18卷

所以不能用记录前一级子震源位置的办法来确定射线路径,而是利用计算出的单元边界上节点的最小旅行时,根据互换原理,从接收点开始,反向确定满足Fermat原理的最小旅行时点的位置,直至进行到震源点所在单元为止.具体步骤为①在接收点所在单元,利用该单元边界上的所有两两相邻的已知旅行时节点,由(1)和(2)式计算出单元边界上至接收

点的各局部最小旅行时的点,从中选出旅行时为最小的点,此点便是所要求的初至射线路径与该单元边界的交点,如图5(a)所示.②以该交点为新的接收点,重复步骤①,直至进行到震源所在单元为止,如图5(b)所示.③把震源点和找出的所有交点及接收点依次相连,如图5(c),就得到了初至射线路径

.

图5　确定波传播射线路径的过程

虚线表示单元边界,星点表示震源点,圆点表示已知旅行时节点,三角形点表示接收点及实际射线与单元边界的交点,.

Fig.5　Processof]Mcalculationofseismicrays[J].Geophysics,

1991,56(1):59～67.

[4]　KlimesL,KvasnichaM.32Dnetworkraytracing[J],Geophys.J.

Int.,1994,116:726～738.

[5]　ZhangJ,ToksozMN.Nonlinearrefractiontraveltimetomography[J].

Geophysics,1998,63:1726～1737.

[6]　MoserTJ,VanEckT,NoletG.hypocenterdeterminationinstrongly

heterogeneousearthmodelsusingtheshortestpathmethod[J].J.Geo2phys.Res.,1992,97:6563～6572.

[7]　FischerR,LeesJM.Shortestpathraytracingwithsparsegraphs[J].

Geophysics,1993,58:987～996.

[8]　VanAvendonkHJA,HardingAJ,OrcuttJA,HolbrookWS.Hy2

bridshortestpathandraybendingmethodfortraveltimeandraypathcalculations[J].Geophysics,1982,66(2):648～653.

[9]　JohnsonDB.Efficientalgorithmsforshortestpathsinsparsenetworks

[J].JournaloftheACM,1977,24:1～13.

[10]　GalloG,PallottinoS.Shortestpathmethods:Aunifyingapproach

[J].MathematicalProgrammingStudy,1986,26:38～64.

[11]　常旭,刘伊克.地震正反演与成像[M].北京:华文出版社,

2001.

[12]　PapazachosC,NoletG.PandSdeep

velocitystructureoftheHel2

lenicareaobtainedbyrobustnon2linearinversionoftraveltimes[J].J.Geophys.Res.,1997,102:8349～8367.

[13]　VanAvendonkHJA,HardingAJ,OrcuttJA,McClainJS.A22

DtomographicstudyoftheClippertontransformfaults[J].J.Geo2phys.Res.,1998,103:17885～17899.

4　结　　语

法,线的模拟计算,尤其适于二维、三维地震层析成像.随着网络节点的增加,所用的计算时间也快速增加,特别对三维问题需占用很多的计算机时间和存贮空间;如减少单元边界上的节点数,会产生不精确的射线路径和旅行时.因此,在满足一定精度要求的情况下,提高效率是最短路径算法首要的发展方向.混合射线追踪方法能得到精确解,其效率不但取决于最短路径法,还与射线弯曲法及震源检波点对的数目有关.我们的动态网络最短路径射线追踪方法是对最短路径算法的发展,与原方法相比,在单元网格相同的情况下,具有较高的精度;在相同精度要求的情况下,能大大提高效率和减小计算机存贮空间.参　考　文　献(References):

[1]　DijkstraEW.Anoteontwoproblemsinconnectionwithgraphs[J].

Numer.Math.,1959,1:269～271.

[2]　NakanishiI,YamaguchiK.Anumericalexperimentonnonlinearim2

agereconstructionfromfirst2arrivaltimesfortwo2dimensionalislandarcstructure[J].J.Phys.Earth,1986,34:195～201.

本文来源：https://www.bwwdw.com/article/jxfq.html