基于旅行时线性插值的地震射线追踪算法
用于地震射线追踪的旅行时线性插值法
用于地震射线追踪的旅行时线性插值法Asaka.,E;蒋录全【期刊名称】《石油物探译丛》【年(卷),期】1993(000)004【摘要】本文提出了一种称为旅行时线性插值(LTI)的新的射线追踪法。
对于二维速度结构来说,用该方法计算旅行时积追踪射线路径比其它常规方法更为快速精确。
LTI 公式仅用于二维单元模型边界上旅行时和射线路径的计算。
为此认为射线路径在一个匀速单元中总为直线。
这种近似法在层析成像分析中较为适用。
LTI 算法分两步进行。
第一步计算所有单元边界上的旅行时,第二步对全部震源-接收点对射线路径进行追踪。
假设单元边界上任一点的旅行时为其相邻两离散点算出的旅行时的线性插值,这时可以费马原理为准则从若干可能的情况中选定正确的旅行时和射线路径。
LTI 法与激发法及程函方程有限差分法(FDM)的数值计算结果对比表明,用LTI 法计算旅行时和射线路径其有速度快和精度高的优点。
LTI 法可看作是常规的程函方程有限差分法的高级形式,因为用LTI 可以直接导出 FDM 公式。
在推导过程中,从理论上可以看出 LTI 法要比 FDM 法精确。
并且使用LTI 法可以避免Vidale 法中由于速度突变产生的数值不稳定现象。
【总页数】9页(P1-9)【作者】Asaka.,E;蒋录全【作者单位】不详;不详【正文语种】中文【中图分类】P631.4【相关文献】1.起伏地表地震波旅行时混合网格线性插值射线追踪计算方法 [J], 王琦;朱盼;叶佩;李勤;李庆春2.基于旅行时线性插值的地震射线追踪算法 [J], 赵改善;郝守玲3.地震波旅行时非线性/线性联合插值法 [J], 彭直兴;张芬;周熙襄;钟本善4.旅行时线性内插地震射线追踪 [J], Asaka.,E;严又生5.地震射线追踪的线性走时扰动插值法 [J], 李同宇;张建中因版权原因，仅展示原文概要，查看原文内容请购买。
基于快速推进迎风双线性插值法的三维地震波走时计算
基于快速推进迎风双线性插值法的三维地震波走时计算孙章庆;孙建国;岳玉波;江兆南【摘要】三维地震波走时计算技术是三维地震反演、层析成像、偏移成像等诸多地震数据处理技术中非常重要的正演计算工具.为了获得精度高且兼顾效率的三维走时计算方法:首先,在常规双线性插值公式推导过程中,充分利用平面波双线性假设的结论,获得了二元极小值超越方程的解析解,进而推导出了准确的局部走时计算公式,同时构造性地证明了该计算公式满足地震波的传播规律和Eikonal方程;其次,引入迎风差分的基本思想,提出迎风双线性插值的局部走时计算策略,该计算策略能简化算法、提高效率且保证无条件稳定性;然后,将上述计算公式和迎风双线性插值策略与常规快速推进法中的窄带技术结合,获得了一种新的基于快速推进迎风双线性插值法的三维地震波走时计算方法;最后,通过精度和效率分析检验了新算法的精度、效率和正确性,并通过计算实例验证了算法在面对复杂介质时的稳定性和有效性.【期刊名称】《地球物理学报》【年(卷),期】2015(058)006【总页数】13页(P2011-2023)【关键词】三维;平面波双线性假设;迎风双线性插值;窄带技术;地震波走时计算【作者】孙章庆;孙建国;岳玉波;江兆南【作者单位】吉林大学地球探测科学与技术学院,长春 130026;国土资源部应用地球物理综合解释理论开放实验室—波动理论与成像技术实验室,长春130026;吉林大学地球探测科学与技术学院,长春 130026;国土资源部应用地球物理综合解释理论开放实验室—波动理论与成像技术实验室,长春130026;东方地球物理公司物探技术研究中心,河北涿州 072751;重庆市地质矿产勘查开发局208水文地质工程地质队,重庆408300【正文语种】中文【中图分类】P6311 引言地震波走时是描述地震波运动学特征的一个很重要的物理参数，其刻画了地震波由震源点激发后到达地下介质空间中各个位置坐标点所需耗消的走时，而该走时和地下介质的速度参数的分布情况密切相关，所以地震波走时计算方法常常被应用于走时反演、层析成像、偏移成像、速度分析等一些重要的地震数据处理技术中.地震波走时计算的精度、效率、稳定性等性质直接影响着这些地震数据处理的效果和效率（Sun，1998，1999，2000，2004；井西利等，2007；瞿辰等，2007），同时考虑到三维地震勘探正在全面大范围展开的大背景，所以对三维地震波走时计算进行研究具有很重要的理论意义和实际价值.在走时计算方面，目前采用的主要方法有两点射线追踪法（Julian and Gubbins，1977）、最短路径射线追踪法（Moser，1991；张美根等，2006）、波前构建法（Vinje et al.，1993）、走时插值法（Asakawa and Kawanaka，1993；Wang and Ma，1999）、有限差分法（Vidale，1988；Sethian and Popovici，1999；Sun et al.，2011）、辛几何算法（秦孟兆和陈景波，2000；高亮等，2000）等.其中，两点射线追踪法在早期天然地震数据处理中起着非常重要的作用，但在现代三维地震数据处理中其计算效率相对较低.最短路径射线追踪法具有不错的走时计算精度，但由于其需要设置很多网格线上的计算节点，所以在求解三维问题时其需要耗费相对更大的计算和存储成本.波前构建法能同时计算走时和射线路径，并且其计算精度也相对不错，不过在三维算法时，其中的网格定位、插值问题等都相对非常复杂.实际上，走时插值法和有限差分法是目前被广泛采用的三维走时算法，其中有限差分法具有计算效率高、局部走时算法简洁等优势，但其计算精度有限，而具有更高计算精度的高阶有限差分需要花费相对更多的计算成本.走时插值法具有不错的计算精度，但是其算法的整体实现策略的效率相对有限差分法低很多.因此，本文将采用走时插值法作为局部走时算法，而在整体实现策略上采用有限差分算法中的快速算法.目前，三维走时插值法主要采用B样条插值法（张东等，2013）、双线性插值法（李培明等，2013；刘锋等，2012；梅胜全等，2010；张东等，2009）等算法.在双线插值法中，有两方面关键算法.一方面是局部走时算法，其包括计算公式和局部计算策略等技术环节；另一方面是算法的整体实现策略，也即算法的整体实现步骤.在局部走时算法方面，目前普遍认为双线性插值法求解极小值的方程为一个超越方程，无法求出解析解（其可能解分布在一条双曲线上），所以一些学者分别采用网格界面剖分法（张东等，2009）、快速插值法（李培明等，2013；梅胜全等，2010）、最速下降法（刘锋等，2012）来获得一个约束条件下接近于真实解的一个近似解.而在走时计算的整体实现步骤方面，目前普遍采用“向前处理”的按行列扫描计算（张东等，2009）、最短路径射线追踪算法的实现步骤（李培明等，2013；梅胜全等，2010）等.针对以上两个方面的关键技术，笔者充分利用平面波双线性假设的结论，推导出了双线性法求解极小值超越方程的解析计算公式，同时构造性地证明了该解析计算公式的正确性和相对应的地震波传播的物理规律.此外，笔者还借鉴快速推进法中的迎风差分格式的构建思想，提出了迎风双线性插值的局部走时计算的实现策略，并以快速推进法中的窄带技术（这是一种满足地震波传播基本规律的、计算效率高且灵活的波前扩展算法）作为算法的整体实现策略.最后的精度分析和计算实例表明了本文算法的正确性和在面对复杂介质时的稳定性和有效性.2 算法的描述为了获得一种精度高且兼顾效率的三维走时计算方法，本文将提出一种快速推进迎风双线性插值法的三维地震波走时算法.如图1所示，这是一种网格算法，和以往基于网格的双线插值法、有限差分法、最短路径追踪法等类似，该算法主要包括三个方面的核心内容：局部计算公式、局部实现策略以及整体实现策略.图1 算法的描述Fig.1 The description of the algorithm如图1所示，局部计算公式的建立可描述为如下问题：在计算的某一个时刻一个平面上的点A、B、C、D的走时值已知，且该平面内的走时值为双线性分布，则如何构建计算与该平面邻近的网格节点F的走时值的计算公式；局部实现策略则可以描述为如下问题：在计算过程中，与走时值待求的F点相邻近的、类似于平面ABCD的、走时值已知且分布满足双线性分布的平面可能会有很多种情况，而在不同情况时应该分别采用怎样对应的不同的处理方法；最后，整体实现策略可以描述为如下问题：从源点S处给定的初始条件出发，应该采用怎样的实现步骤或过程来逐次地计算整个计算空间内所有网格节点的走时值.与如上描述的三个问题相对应，这实际上涉及到算法的局部计算公式、算法的局部实现策略、算法的整体实现策略，所以接下来本文将分别详细阐述这三方面问题.最后，还将通过算法的精度分析和计算实例来验证算法的精度、正确性和有效性.3 算法的局部计算公式：双线性插值解析公式为了计算三维地震波走时，首先需要推导建立局部走时计算公式.如上所述，局部走时计算公式的推导可描述为：如图2所示，设点A、B、C、D的走时值已知，分别为：tA、tB、tC、tD；正方体网格的网格间距为h；假设平面ABCD上的走时值为双线性分布；求F点的走时值tF的计算公式.图2 算法的局部计算公式：双线性插值公式Fig.2 The local computationformulas of algorithm：the bilinear interpolation formulas为了推导计算tF的公式，首先以点A为坐标原点，向量AB方向为x轴正方向，向量AD方向为y轴正方向，向量AF方向为z轴正方向，建立局部笛卡尔坐标系，则有点A、B、C、D、F的坐标分别为A（0，0，0），B（h，0，0），C（h，h，0），D（0，h，0），F（0，0，h）；其次，假定到达F点的射线经过了走时值分布为已知的平面ABCD内的任意一点E，设E点的坐标为 E（x，y，0），因为E 点被限定位于平面ABCD内，则有0≤x≤h，0≤y≤h.过E点作平行于y轴的直线，该直线与线段AB、CD分别相交于E1 和E2 点，则点E1、E2 的坐标分别为E1（x，0，0）、E2（x，h，0）.根据平面ABCD上的走时值为双线性分布的假设，线段AB、CD及E1E2上的走时值均为线性分布，则有对上述（1）式进行整理可得其中k1 ＝（tB－tA）/h，k2 ＝（tD －tA）/h，k3 ＝（tA＋tC－tB－tD）/h2.设在局部正方体网格单元内地震波传播速度为均匀，并基于上述到达F点的射线经过平面ABCD内的E（x，y，0）点的假定，可得其中s为当前正方体网格单元内的地震波传播速度的倒数（即慢度）.将式（2）代入（3）中，则有根据Fermat原理，只要确定能使得（4）式取最小值的E（x，y，0）点的位置，即可获得计算tF的公式，所以对（4）式tF分别关于变量x，y求偏导数，并令该偏导数为零，即分析方程组（5）可知：只要求出该方程组的解，即可获得E（x，y，0）点的位置坐标，然后再将其带回到（4）式即可获得计算tF的公式.但是，一些研究已表明方程组（5）实际上是一个超越方程组，无法求出相应的解析解（其可能解分布在一条双曲线上），所以一些学者分别采用网格界面剖分法（张东等，2009）、快速插值法（李培明等，2013；梅胜全等，2010）、最速下降法（刘锋等，2012）来获得一个约束条件下接近于真实解的一个近似解.与以往的双线插值法的近似处理不同，在此笔者将充分利用双线性假设的结论，来获取方程组（5）的精确解析解，并构造性地证明计算tF的公式的正确性以及其背后满足的地震波传播波的物理事实.根据走时插值法的提出者Asakawa的线性假设（Asakawa and Kawanaka，1993），实际上即是平面波的假设，其阐述的是如果当前传播的为平面波，则该波传播经过的二维空间中的走时分布在二维空间中的任意一条网格线上均为线性，而在三维空间中的一个平面上的走时分布则均为双线性.同样基于平面波的假设，并假定在局部的相对较小的一个网格单元内地震波的传播速度为匀速.如图3所示，则无论从那个方向（假定箭头方向为平面波的任意一个传播方向）途经平面ABCD均有tD－tA＝tC－tB，因为如图3所示，设该平面波与平面ABCD相交于如图3所示的一系列平行的虚线，其中经过点A、B、C、D的虚线分别位于平面波tA、tB、tC、tD 时刻的波前上，因为这些波前是相互平行的，并且网格四边形ABCD为一正方形，所以必有tA时刻的波前传到tD时刻花费的走时等于tB时刻的波前传到tC时刻花费的走时，也即是tD－tA＝tC－tB. 图3 平面波假设分析Fig.3 The analysis of plane wave assumption基于上述双线性假设，即平面波假设的进一步结论tD－tA＝tC－tB，可以得出k3＝0.将该结论带回到（5）式中，则（5）可简化为再分析（6）式，若在k1、k2均不为零的情况下，则用下式除以上式再整理有：y ＝k2x/k1，即方程组（6）的解实际上是位于一条直线上（而并非一条双曲线上），并且将该直线代回方程组（6）中的任意一个方程，均可以获得解析的x、y的解.不过要获得各种情况下的方程组（6）的解，还得分情况讨论k1、k2取不同值时，方程组（6）的解的情况：①k1＝0且k2＝0时此时，将k1、k2代回方程组（6）可得：x＝y＝0，再将该结论代回（4）式中可得计算tF的公式为②k1＝0且k2＜0时在此种情况下，将k1＝0代入方程组（6）可得x＝0，再将该结论代回方程组（6）中可将其简化为融入0≤y≤h的限定条件，求解（8）式可得其中Δt＝tD－tA，且需满足：（此是基于0≤y≤h的限定条件）.综合x＝0的结论，可知 E点的坐标为0），再将其代入（4）式中即可获得此种情况下计算tF的公式：③k2＝0且k1＜0时这种情况与情况②相似，不同之处仅在于交换了变量x与y的角色，所以同理情况②的推导过程可得：其中：Δt＝tB－tA，且需满足：（此是基于0≤x≤h的限定条件）.同时E点的坐标为计算tF 的公式为：④k1＜0且k2＜0这种情况下，直接求解方程组（6），即将上述结论代入方程组（6）的上式，并融入0≤x≤h、0≤y≤h的限定条件，则可得方程（6）的解：其中，Δt1＝tB－tA、Δt2＝tD－tA，且需满足：0≤（此是基于0≤x≤h、0≤y≤h 的限定条件）.将（14）式代回（5）式中，即可获得计算tF的公式：同时，此时E点的坐标为：⑤k1＞0或k2＞0这种情况下，因为上述将E点限定于平面ABCD内的限定条件有0≤x≤h、0≤y≤h，所以很容易分析出，此时方程组（6）无解.这种无解的情况表面上表明：上述计算公式可能出现不稳定的情况.不过，这种所谓的不稳定的情况，在将上述算法融入一定的局部实现策略和算法的整体实现策略后，是不可能出现的，关于该问题将在后续内容中阐述.经过以上的推导，实际上已经获得了当k1、k2取不同值的情况时，点E（x，y，0）的位置坐标和tF的计算公式.如图4所示，下面将构造性地证明这些公式的正确性，以及它们对应的不同情况下的地震波传播的物理事实.图4 双线性插值公式的分析Fig.4 The analysis of bilinear interpolation formulas上述情况①：当k1＝0且k2＝0时，即有tB＝tA，tD＝tA，也即tA＝tB＝tD，该式对应的物理事实为点A、B、D位于地震波传播过程中的同一个等时面上.此外，根据不在同一条直线上的三点确定一个平面的数学事实和上述基于的平面波假设，则有与A、B、D点位于同一平面上的C点也位于该等时面上，所以可以进一步得出tA＝tB＝tC＝tD，也即平面ABCD为平面波的一个等时面.同时，根据局部笛卡尔坐标系可知，直线AF与平面ABCD垂直，再根据射线路径和等时面相互垂直的物理规律，可以得出在这种情况下平面波沿着z轴正方向传播，而经过平面ABCD到达F点的射线与平面ABCD相交于A点，所以此时自然有：x＝y＝0，tF＝tA＋sh，也就是说上述情况①得的方程组（6）的解和tF的计算公式符合在这种情况下平面波的传播规律，即是正确的.上述情况②：当k1＝0且k2＜0时，即有tB＝tA，tD＜tA，此时再根据双线性假设，因为tB＝tA，所以平面ABCD上的走时分布沿x方向没有变化，而仅沿着y方向有变化，所以此时的插值问题实际上已简化为二维插值问题，而E点则被固定到距F更近的线段AD上，即有x＝0.此外，基于二维空间中Asakawa提出的走时线性插值法同样可以获得其中Δt＝tD－tA.再分析tF的计算公式，实际上其可以变形为该方程实际上是二维Eikonal方程（∂t/∂z）2＋（∂t/∂y）2＝s2在网格上以A点为中心的差分离散方程，也就是说此处得到的tF的计算公式，在局部上满足Eikonal方程，所以其是正确的.同时上述0的限定条件还能保证tF的计算公式不会出现负数开平方的不稳定情况.上述情况③：当k2＝0且k1＜0时，此时的情况与情况②是类似的，也可以简化为一个二维插值问题，不同之处仅在于坐标方向x与y交换了角色，所以在此不再重复阐述.上述情况④：当k1＜0且k2＜0时，即有tB＜tA，tD＜tA，在这种情况下插值问题才真真变为一个三维双线性插值问题，在此条件下方程组（6）有解.求解该方程组即可获得E点的坐标为x＝－hΔt1/同时将它们代回（5）式得到tF的计算公式：tF＝tA＋其中Δt1＝tB－tA、Δt2＝tD－tA. 同样其可以变形为该方程实际上是三维 Eikonal方程（∂t/∂x）2 ＋（∂t/∂y）2 ＋（∂t/∂z）2＝s2在网格上以A点为中心的差分离散方程，也就是说此处得出的tF的计算公式在局部上是满足Eikonal方程的，所以其是正确的.此外，再仔细分析E点的坐标可得：y＝kx，其中k＝Δt2/Δt1，也就是说，实际上E点是位于坐标平面xAy内的一条直线上，该直线的斜率为k.k的取值实际上和已知条件tA、tB、tD的取值情况有关：k＝Δt2/Δt1 ＝（tD－tA）/（tB－tA）.如图4所示，tD ＜tB时有k＞1，此时直线途径ACD半区域，也即E点位于走时分布较小（因为双线假设和tD＜tB）的区域.这实际上与Fermat原理和地震波传播的基本规律是相符合的，即地震波总是从走时值相对更小的区域传向下一个区域.而当tB＜tD时有k＜1，此时有如上同样的结论.最后当tB＝tD时，则k＝1，此时E点位于正方形ABCD的对角线AC 上.综上所述，在情况④下得出的E点的坐标位置是满足地震波的传播规律和Fermat原理的，tF的计算公式在局部上满足Eikonal方程.同时0≤2s2 h2/3的限定条件，还能保证tF的计算公式不会出现负数开平方的不稳定情况.上述情况⑤：当k1＞0或k2＞0时，即有tB＞tA，tD＞tA，此时上述的结论为无解.实际上，再回顾方程组（6），此时并非方程组（6）无解，只是解的情况为x＜0，y＜0，即E点不在平面ABCD内.实际上，为了公式推导的方便和局部网格计算的限定，将E点限定于走时值已知的平面ABCD内，而在该限定条件下x＜0，y＜0的解不能满足该限定条件，所以致使方程组（6）无解.事实上，x＜0，y＜0的解的情况说明：在这种情况下E点位于其他的邻近的网格单元内.同时tA＜tB、tA＜tD，再加上双线性假设，暗含了其他相邻的网格单元上的分布的走时值会比网格单元ABCD上的更小.而根据Fermat原理，此时计算tF的公式则应该通过其他网格单元上推导获得，也即在计算时应该采用其他网格单元进行计算.该问题实际上引伸出了局部走时计算策略的问题：如图5所示，在进行三维局部走时计算时，包围被算点F周围的类似于ABCD的网格单元实际上有24个，从几何的角度讲，到达F点的射线必与用于封闭F点的类似于ABCD的24个网格单元中的一个相交，而该交点正是我们上述需要求解的E点.当然要是在24个网格单元中逐次搜寻计算E点，必定是一种繁琐且计算效率低下的算法.为了解决该问题，此处借鉴迎风差分的思想，提出了一种非常简洁的迎风双线性插值的局部实现策略，下面将详细阐述该策略.4 算法的局部实现策略：迎风双线性插值法如图5所示，在三维走时计算的某一时刻，当前需要计算的是点F的走时值tF.如图2所示，上述推导了网格单元ABCD走时分布为已知的情况下点F的走时值tF 的计算公式.但实际上，在被算点F周围的24个类似于ABCD的网格单元中，走时分布为已知的网格单元可能有很多个，那么究竟应该采用那个网格单元作为已知条件是合理的和稳定的？在此，引入迎风差分的基本思想（Sethian and Popovici，1999），该思想是指地震波总是从地震波走时值分布更小的区域向走时值未知的区域传播，而在局部数值计算实现时就应该迎着走时值分布更小的区域去做差分代替微分的近似.迎风差分格式的基本思想实际上是Huygens原理和Fermat原理的综合利用，主要体现在其隐含的将被算点周围的走时值最小的网格节点当作子震源点，也正因为此迎风差分格式具有无条件稳定性.同样借鉴该思想，在此提出迎风双线性插值法，该方法的核心借用迎风差分格式选取最小邻近走时点的方法，选取被算点周围的最小走时分布网格单元，然后再采用类似于图2所示的推导的公式进行局部走时计算.总体上讲，实际就是要构建相应的迎风双线性插值公式.如图5a所示的局部正方体网格，以点F为中心的大正方体的六个面分别由以A1、A2、A3、A4、A5、A6为中心的4个类似于图2中的ABCD的网格单元组成，这样就总共有24个类似于ABCD的网格单元.在当前计算点F的走时值tF的时刻，这24个网格单元上的26个网格节点的走时值可能为已知，可能为未知.由后续将阐述的算法的整体实现策略的初始化步骤（详见后续阐述内容的“算法的整体实现步骤：快速推进法”部分）可知，若网格节点的走时值为未知（即没有被计算出来），其值为固定的一个非常大的数（这个数为人为规定的，其值远大于所有网格节点最终可能算出的走时值），若网格节点的走时值被计算完成，其值即为最终的走时计算结果.在一般情况下，当前被算点F周围的26个网格节点的走时值为一部分被计算完成，一部分未开始计算（如图6a所示）.而在被计算完成的部分里面，先完成的走时值一定是小于后完成的走时值（该结论由后续内容的窄带技术的阐述可知），也就是说与点F邻近的26个网格节点均有一个走时值，而用于计算tF的网格单元的走时分布应为最小.具体挑选最小走时分布区域的方法为：①确定与点F直接相邻的6个网格节点的走时值最小的网格节点Amin，其中有 tAmin ＝min（tA1，tA2，tA3，tA4，tA5，tA6）；② 在以Amin为中心的平面上分别确定其他两个方向上的最小走时网格节点Bmin，Dmin，若设Amin＝A5则分别有tBmin ＝min（tB51，tB52），tDmin ＝ min（tD51，tD52）；③Amin，Bmin，Dmin所在的网格单元即为最终的类似于公式推导过程中的网格单元ABCD，与该网格单元相对应，tF的计算公式相应地变为图5 算法的局部实现策略：迎风双线性插值Fig.5 The local implementation strategy of algorithm：the upwind bilinear interpolation其中Δt1＝tBmin－tAmin，Δt2＝tDmin－tAmin.公式（16）即为最终的迎风双线性插值公式，与常规的双线性插值公式相比，该公式实际上是一种优化的公式.该公式隐含着一个优化的局部实现策略，而该优化的局部实现策略将Huygens原理和Fermat原理融入到算法的局部计算过程中，使得算法是在满足地震波传播的基本规律的基础上实现的.与迎风差分格式类似，此处的迎风双线性插值法也具有无条件稳定性.5 算法的整体实现策略：窄带技术上述内容分别建立了算法的局部走时计算公式和局部实现策略，实际上要想完成整个计算区域的地震波走时计算，还需要拟定一个算法的整体实现步骤，即算法的整体实现策略.算法的整体实现策略往往决定着算法的稳定性和计算效率.基于对计算效率、无条件稳定性、对迎风双线性插值的适应性等因素的综合考虑，在此引入快速推进法中的窄带技术（Sethian and Popovici，1999）作为算法的整体实现策略.如图6b所示，窄带技术的核心思想是采用一个包含当前所有波前点的窄带来近似模拟地震波前，通过窄带的扩展演化来模拟地震波前的传播过程，其充分考虑了地震波前传播过程中的熵守恒理论，同时此处因为迎风差分思想的融入，所以窄带技术具有无条件稳定性.窄带技术在实现过程中将计算空间中的网格节点分为三种类型：完成点，其走时计算已经完成的点，在图6b中用黑色填充的点表示；窄带点，当前窄带内的点，也即近似波前上的网格节点，在图6b中用灰色填充的点表示；远离点，窄带扩展还未到达的点，也即还未进行走时计算的网格节点，在图6b中用白色填充的点表示.图6 算法的整体实现策略：窄带技术Fig.6 The global implementation strategy of algorithm：the narrow band technique窄带技术的实现过程可以分为初始化和扩展计算两个部分.在初始化时，震源点为唯一的完成点，其走时值为零；与震源点直接相邻的正方体网格单元内的网格节点为初始窄带点，它们构成初始窄带，即初始波前，它们的走时值通过局部的解析公式求取；其它所有剩下的计算空间中的网格节点为远离点，其初始值设置为一个相对很大的数（远大于最终可能计算出来的最大的走时值即可，这主要是为了实现上。
基于双线性插值的三维地震波旅行时计算
广
Ｉ ， ”
图 １ 平 面 ＸＺ内矩 形 单 元 旅 行 时 分 布 图 Ｏ
收 稿 日期 ： ０ ８— ３—２ ２０ ０ ８
作者简介 ： 彭直兴 （９ ３ ， 汉族 ） 重庆万州人 ， 师 ， 士研究生 ， １７ 一） 男（ ， 讲 博 主要从事矿产普查与勘探 、 地球探测 与信息技术研究。
解 ２ Ｄ和 ３ Ｄ程 函方 程 ¨ ’Ｊ他 应用 了地 震波 沿平 一 ． ， ＂ 面波前传 播 的理论 。虽然 这两 种方 法 在射 线 间 的插
化的， 么， 那 单元 内任一 点的旅行 时可用 单元 角点 （ 网
格节点 ）上 的旅行 时通 过双线 性 插值来 表 示 。 于 三 对
Ｖｏ ． ０ Ｎｏ ５ １３ ． ０ｃｔ ２ ８ ． ｏｏ
文 章 编 号 ：１０ ２ ３ （０ ８ ０ ０ ８ ０ ０ ０— ６４ ２ ０ ）５— ０ ５— ３
基 于双 线性 插值 的三 维 地 震 波 旅 行 时 计 算
彭直兴 ， 忠民 沈
（ 成都理工大学能源学院 ， 四川 成都 ６ ０ ５ ） １０ ９
ｍ ｔｏ ） 往往 会 收 敛 到 一 个 局部 最 小 旅 行 时路 径 ， ｅ ｄ， ｈ 难 以找 到全局 最 小解 ， 不 能 处 理 射 线 盲 区 以及 临 也 界折 射波 和绕射 波对应 的最 小旅 行 时路 径 Ｊ 。近
１ 旅 行 时双 线 性 插 值 基 本 原 理
假设在任 意 四边形单元 内， 旅行 时分 布是 线性 变
推 导 了相应 的计 算公 式。通过数值模拟计 算、 与解析解的对 比， 不仅 证明 了该 方法原理 的正确性 ， 而且还证 明 了该 方
法 原 理 具 有较 高 的计 算精 度 。 关 键 词 ： 线 追 踪 ； 震 波旅 行 时 ； 线 性 插 值 ； 维 结构 ； 值 模 拟 射 地 双 三 数
射线追踪算法
射线追踪算法
射线追踪的算法有多种,比如
最短路径法,
有限差分方程法,
旅行时线性插值法(LTI,Li near Traveltime Interpolation)等.
反演求解的方法有
反投影法(BPT),
代数重建法(ART),
同时迭代重建法(SIRT),
奇异值分解法(SVD),
最小二乘QR分解法(LSQR)等多种算法。
射线追踪: 地震波沿着一条传播时间最短的路径进行传播。
该方法把模型离散成均匀的正方形单元，旅行时和射线路径的确定只与单元边界上的点有关。
该方法把模型离散成均匀的正方形单元，旅行时和射线路径的确定只与单元边界上的点有关。
追踪出来的射线也不同。
随着循环次数的增多，出现了明显的回波现象。
射线追踪的理论基础是,在高频近似条件下,地震波场的主能量沿射线轨迹传播.传统的射线追踪方法,通常意义上包括初值问题的试射法(Shootingmethod)和边值问题的弯曲法(Bendingmethod).试射法根据由源发出的一束射线到达接收点的情况对射线出射角及其密度进行调整,最后由最靠近接收点的两条射线走时内插求出接收点处走时.弯曲法则是从源与接收点之间的一条假想初始路径开始,根据最小走时准则对路径进行扰动,从而求出接收点处的走时及射线路径。
起伏地表地震波旅行时混合网格线性插值射线追踪计算方法
起伏地表地震波旅行时混合网格线性插值射线追踪计算方法王琦;朱盼;叶佩;李勤;李庆春【摘要】为了提高在起伏地表和复杂构造条件下地震波旅行时的计算精度和效率,对完全矩形网格剖分的旅行时线性插值(LTI)方法进行了改进,提出用矩形网和不规则四边形网相结合的方法离散速度模型,对适用于矩形网格的局部旅行时计算公式进一步推导,得到适用于混合网格的局部旅行时计算公式,并证明了该公式可稳定求解.结合分区多步计算技术将该方法扩展为基于混合网格的分区多步LTI方法.数值计算结果表明,该方法不但可以灵活地处理剧烈起伏的地表和构造复杂的速度界面,而且旅行时和射线路径的计算结果能保持较高的精度.【期刊名称】《石油地球物理勘探》【年(卷),期】2018(053)001【总页数】12页(P35-46)【关键词】旅行时线性插值法;混合网格;起伏地表;多次波;转换波【作者】王琦;朱盼;叶佩;李勤;李庆春【作者单位】长安大学地质工程与测绘学院,陕西西安710054;长安大学地质工程与测绘学院,陕西西安710054;中国国土资源航空物探遥感中心,北京100083;西安科技大学地质与环境学院,陕西西安710054;长安大学地质工程与测绘学院,陕西西安710054【正文语种】中文【中图分类】P6311 引言旅行时线性插值(Linear traveltime interpolation 简称LTI)是Asakawa等[1]提出的基于线性假设网格单元扩展射线追踪算法。
由于该方法计算速度快、原理简单，得到了很大发展。
赵改善等[2]结合界面二次源法将该方法推广于追踪反射波旅行时；聂建新等[3]将旅行时二次插值与线性插值方法联合，降低了累积误差；张赛民等[4]用抛物线插值取代线性插值，改善了因线性引起的误差；张东等[5]为了求得最小走时，在向前处理过程中采用了多方向循环的计算方法。
近年来多次波越来越引起人们关注，对多次波的正确模拟和认识是压制多次波干扰的前提。
2 射线追踪法合成地震记录
初始入射角：方位角
2-3 试射法射线追踪
• 起伏界面试射法——问题2
射线在交点处透射线（或反射线）的方向计算问题
计算透射方向
2-3 试射法射线追踪
• 起伏界面试射法——问题2
射线在交点处透射线（或反射线）的方向计算问题
2-3 试射法射线追踪
• 起伏界面试射法——实例
2-3 试射法射线追踪
• 试射法追踪收敛方法
2-5 二维层状模型描述
• 问题分析
替代方法：封闭结构模型（研究生课程）
2-6 程序设计及演示
• 程序说明
2-6 程序设计及演示
• 程序说明
2-6 程序设计及演示
• 程序说明
2-6 程序设计及演示
• 程序说明
2-6 程序设计及演示
• 程序流程图
程序—第1页
2-6 程序设计及演示
程序—第2页
• 起伏界面试射法——问题1
射线与起伏界面的交点计算问题
思路：将问题简化为卡尔 丹公式进行计算
2-3 试射法射线追踪
• 起伏界面试射法——问题1
射线与起伏界面的交点计算问题
将目标公式转换为卡尔丹——问题1
射线与起伏界面的交点计算问题
思路
2-3 试射法射线追踪
• 起伏界面试射法——模型——界面描述
实际处理
2-3 试射法射线追踪
• 起伏界面试射法——问题分解
差异
2-3 试射法射线追踪
• 起伏界面试射法——问题1
射线与起伏界面的交点计算问题
第1种情况
2-3 试射法射线追踪
• 起伏界面试射法——问题1
射线与起伏界面的交点计算问题
第2种情况
第3种情况
利用线性旅行时插值射线追踪计算近地表模型初至波走时
利用线性旅行时插值射线追踪计算近地表模型初至波走时涂齐催;刘怀山【摘要】利用线性旅行时插值射线追踪对近地表模型进行正演计算,可以快捷、准确地获得初至波走时和射线路径.由于该算法计算的初至波不局限于折射波,因此很好地解决了浅层折射勘探中的低速"隐蔽层"问题;而且,由于该算法是基于网格划分和线性插值,因此它不仅可以追踪任意复杂介质的初至波,而且可以使得追踪的初至波射线路径逼于真实,避免了同类算法直接连接网络节点形成射线路径的缺陷(路径过于弯折,计算走时偏大).将LTI算法同其他几种算法的追踪结果进行的对比和分析表明,LTI算法在计算初至波走时和射线路径方面较其他算法更为精准、稳定,是一种有效的射线追踪方法.【期刊名称】《物探与化探》【年(卷),期】2006(030)002【总页数】6页(P148-153)【关键词】线性旅行插值;射线追踪;近地表;初至波;折射波;屏蔽层;网格剖分【作者】涂齐催;刘怀山【作者单位】中国海洋大学,海洋地球科学学院,山东,青岛,266003;中国海洋大学,海洋地球科学学院,山东,青岛,266003【正文语种】中文【中图分类】P631.4随着东部油田和陆上油田勘探程度的不断提高，西部沙漠、戈壁、黄土高原、山地和海上滩浅海[1]等复杂地区的油气勘探正成为主要勘探目标，在这些地区进行地震勘探静校正成为影响勘探效果的关键技术之一。
要想尽可能解决静校正问题，必须进行复杂近地表结构成像研究，这是一项很有意义的工作；而作为反演成像基础的正演模拟在地球物理学中一直占有很重要的地位。
因此，在进行反演成像之前，有必要进行正演模拟试验和分析。
对于近地表结构，长期以来，基本上是利用层状介质模型的折射波初至来进行研究的，并发展了一些层状介质模型的折射波正演方法。
徐涛等研究了复杂介质折射波射线追踪方法[2]，但在出现“屏蔽层”(中间低速层)时，此方法失去作用。
刘康和等[3]对此作了分析和讨论，建议实际工作中应结合初至折射和浅层反射走时来进行综合反演和解释，但浅层勘探的反射信号用来反演和成像并不可靠，原因是在信噪比低的情况下很难对其精确拾取。
基于改进的双线性旅行时插值的三维射线追踪
物探 化探 计 算技 术
２１年３ ００ 月
文章编 号 ：ｌ０ — ｌ４ （ ０ ０ Ｏ —０ ５ — ０ ｌ ７ ９ ２ １ ）２ １２ ６
基 于 改进 的双 线 性 旅行 时插值 的 三维射 线 追踪
梅胜全， 飞， 邓 钟本善， 周熙襄
（ 成都理 工 大学 摘 信 息工程 学院 ， 四川 成都 ６０５ ） １０９
线人 射 点 ， 图 ２ ａ ， 个最 小旅 行 时插 值 就是 检 如 （ ）这
波点 的初 至旅 行 时 。
各单元的插值 ， 直至激发 点， 得到一条完整 的初至 射线 路径 。
１１ 向前 处理 ： ． 计算 模型 节点 时 间场
（ ）根据当前交点的位置， ２ 利用旅行时插值算
法 ， 续找 出下 一个 网格 单 元 的交 点 如 图 ２ ｂ 所 继 （）
决 了三 维射 线 追踪算 法的计 算低 效率 ， 线精度 不 高的 问题 。 射
关 键词 ：线性旅行 插值 ； 射线 追踪 ； 网格剖分 ； 至波旅行 时 初
中图分类 号 ：Ｐ６ １４ ４ ３ ． ３
文 献标 识码 ：Ａ 巨大 障碍 。那 么 ， 如何 有效地 解决 三维地震 走 时计
பைடு நூலகம்
２期
梅 胜全 等 ： 于改进 的 双线性旅 行 时插 值 的三维射 线 追踪 基
１３ ５
（ 激发 点单元 示意 ａ ）
（） ｂ 沿六个面波前扩展
（） ｃ 计算 出次震源点
图 １ 向 前 处理 过 程
Ｆ ｇ １ Ｆ ｒ ｒ ｒ ｃ ｓ ｉｇ ｐ ｏ ｅ ｕ ｅ ｉ ． ｏ ｗａｄ ｐ ｏ ｅ ｓｎ ｒ ｃ ｄ ｒ
０ 前 言
基于Curvelet变换的缺失地震数据插值方法
２０１１年４月第４６卷 第２期 ＊北京市昌平区府学路１８号中国石油大学（北京）ＣＮＰＣ物探重点实验室９０４室，１０２２４９。
Ｅｍａｉｌ：ｇｕｏｃｈａｎｇ．ｌｉｕ＠ｙａｈｏｏ．ｃｏｍ本文于２０１０年７月１５日收到，最终修改稿于２０１１年１月１９日收到。
本项研究得到国家８６３课题（２００６ＡＡ０９Ａ１０２－０９）和国家重大专项课题（２００８ＺＸ０５０２３－００５－００４）联合资助。
·处理技术·文章编号：１０００－７２１０（２０１１）０２－０２３７－１０基于Ｃｕｒｖｅｌｅｔ变换的缺失地震数据插值方法刘国昌＊ 陈小宏 郭志峰 刘华锋 高建军（中国石油大学（北京）油气资源与探测国家重点实验室，北京１０２２４９）刘国昌，陈小宏，郭志峰，刘华锋，高建军．基于Ｃｕｒｖｅｌｅｔ变换的缺失地震数据插值方法．石油地球物理勘探，２０１１，４６（２）：２３７～２４６摘要 Ｃｕｒｖｅｌｅｔ变换具有局部性、多尺度性和多方向性，在处理非稳态地震信号中具有较大的优势，可以利用Ｃｕｒｖｅｌｅｔ变换的压缩特性重建缺失的地震数据。
本文首先分析了基于稀疏变换的缺失地震数据插值的基本原理，在反问题正则化理论框架下，针对Ｌ２范数约束和Ｌ１范数约束条件，分析了两种约束的差异，着重阐述了基于Ｃｕｒｖｅｌｅｔ变换的Ｌ１范数约束的插值方法，其优点在于对非线性同相轴不需要分窗口处理，并将凸集投影算法（ＰＯＣＳ）引入到Ｃｕｒｖｅｌｅｔ变换的插值方法中，通过采用按指数规律衰减的阈值参数加快了迭代收敛的速度。
理论和实际算例验证了Ｃｕｒｖｅｌｅｔ变换插值方法的有效性。
关键词 地震数据 重建 Ｃｕｒｖｅｌｅｔ变换 凸集投影 Ｌ２范数 Ｌ１范数中图分类号：Ｐ６３１ 文献标识码：Ａ１ 引言由于地震采集的原因以及经济条件的限制，地震数据往往是不规则或者稀疏的，因此地震数据规则化和道间插值在地震资料处理和成像中有着重要的作用［１］，如基于多道褶积的自由表面多次波衰减（ＳＲＭＥ）、波动方程偏移等地震资料处理和成像技术尤其需要规则、完全的地震数据体［２～５］。
地震波走时和射线的有限差分计算
地震波走时和射线的有限差分计算
赵成斌;许云
【期刊名称】《华北地震科学》
【年(卷),期】1995(013)004
【摘要】以往都是采用射线追踪的方法计算地震波的走时和射线，但是当速度模型复杂时这种方法存在一些问题。
本文提出另一种计算地震波走时和射线的方法。
该方法从程函方程出发，利用互换原理和Ｆｅｒｍａｔ原理计算出各种波的到时和射线。
解决了射线追踪方法存在的问题。
为地震波走时和射线的计算以及地震波走时反演开辟了一条新途径。
【总页数】6页(P41-46)
【作者】赵成斌;许云
【作者单位】不详;不详
【正文语种】中文
【中图分类】P315.3
【相关文献】
1.非均匀介质中地震波走时与射线路径快速计算技术 [J], 赵爱华;张中杰;王光杰;王辉
2.地震波走时的有限差分法计算 [J], 朱金明;王丽燕
3.二维复杂介质中地震波走时和射线的计算方法及其应用 [J], 潘纪顺;张先康;赵成斌;张景超
4.不同速度模型条件下地震波走时的射线追踪计算 [J], 嘉世旭
5.有限差分法地震波走时计算 [J], 王华忠;谢海兵
因版权原因，仅展示原文概要，查看原文内容请购买。
地震射线追踪的线性走时扰动插值法
地震射线追踪的线性走时扰动插值法李同宇;张建中【摘要】线性走时插值(LTI)方法假设离散模型单元边界上地震波走时呈线性变化,则单元边界上任意点的走时可通过相邻离散网格节点上走时的线性插值表示.而实际上,走时沿单元边界并非线性变化,当离散单元较大时,线性假设会导致较大计算误差.针对此问题,本文采用线性走时扰动插值方法(LTPI),将单元边界点的实际走时分解为等效匀速介质中的参考走时和走时扰动(后者远小于前者);在离散单元边界上,假设走时扰动线性变化,同时参考走时保持非线性变化,避免了LTI方法的弊端.文中将复杂介质离散成不规则单元,推导了适用于二维不规则单元的线性走时扰动插值公式,形成一种基于LTPI方法的透射波射线追踪方法.不同模型的测试结果表明,相比LTI方法,LTPI方法对复杂介质具有更强的适应性,计算的波前走时和射线路径具有更高计算精度和更强稳定性,在满足一定精度要求的情况下具有更高计算效率.【期刊名称】《石油地球物理勘探》【年(卷),期】2018(053)006【总页数】10页(P1165-1174)【关键词】波前走时;射线追踪;走时扰动;走时插值;不规则单元【作者】李同宇;张建中【作者单位】海底科学与探测技术教育部重点实验室,山东青岛266100;中国海洋大学海洋地球科学学院,山东青岛266100;海底科学与探测技术教育部重点实验室,山东青岛266100;青岛海洋科学与技术试点国家实验室海洋矿产资源评价与探测技术功能实验室,山东青岛266061;中国海洋大学海洋地球科学学院,山东青岛266100【正文语种】中文【中图分类】P6311 引言地震射线追踪及走时计算是地震学的基本问题之一，被应用于层析成像、地震定位、偏移以及地震数据采集设计等领域[1，2]。
传统的射线追踪方法主要针对两点射线追踪问题，可分为试射法[3]和弯曲法[4]。
这两种方法都存在一定局限性：前者需不断试验射线入射角，计算效率较低；后者易于陷入局部最优解。
对旅行时进行抛物型插值的地震射线追踪方法
对旅行时进行抛物型插值的地震射线追踪方法
对旅行时进行抛物型插值的地震射线追踪方法
提出一种对旅行时进行抛物线插值的地震射线追踪方法(简称PTI 方法),它比基于旅行时线性插值方法(简称LTI方法)计算结果更准确.PTI 和LTI方法都是基于2D网格单元模型,用于计算地震波的旅行时和射线路径.首先介绍了相关方法的一些基本概念.旅行时和射线路径都是在网格边界上进行计算的,因此,射线路径在同一恒速网格内是直线.其计算过程有两步.第一步,计算旅行时,第二步追踪射线路径.然后给出了LTI算法的基本公式.因为在炮点网格内可能存在折射波,文章也相应导出了其公式.最后详细推导了PTI方法的公式.通过模型试算对比说明,用PTI方法较LTI算法更精确、更有效,PTI方法是一种很有发展前途的地震射线追踪算法.
作者：张赛民周竹生陈灵君周惠群ZHANG Sai-min ZHOU Zhu-sheng CHEN Ling-jun ZHOU Hui-qun 作者单位：张赛民,周竹生,周惠群,ZHANG Sai-min,ZHOU Zhu-sheng,ZHOU Hui-qun(中南大学信息物理工程学院,410083长沙)
陈灵君,CHEN Ling-jun(中石化中南分公司研究院,410007长沙) 刊名：地球物理学进展ISTIC PKU英文刊名：PROGRESS IN GEOPHYSICS 年，卷(期)：2007 22(1) 分类号：P315.3 关键词：PTI 法 LTI法射线追踪旅行时射线路径网格单元。
线性走时插值射线追踪算法的改进
线性走时插值射线追踪算法的改进卢江波;方志【摘要】在LTI(Linear Travel-time Interpolation)射线追踪算法基础上提出的扩张收缩扫描算法能正确追踪直达波、绕射波和回波的射线路径，但其存在计算效率低、收敛速度慢的问题。
采用交叉扫描方式对扩张收缩扫描算法进行改进，并由此提出了基于交叉扫描方式的扩张收缩扫描改进算法。
理论分析及数值模拟结果表明：改进算法在保留了原扩张收缩扫描算法所有优点的同时，具有更高的计算效率；当模型网格尺寸划分较细时，改进算法在计算效率上的优势更为显著。
%The extension & compaction scanning algorithm based on LTI (Linear Travel-time Interpo-lation)ray tracing algorithm can properly trace the ray path of the direct wave,diffraction wave and echo 's,but it has the following problems:low computational efficiency and slow convergencespeed.Aiming at these problems and using the cross-scan mode to improve the algorithm,an improved extension &compac-tion scanning algorithm based on the cross-scan mode was presented.Theoretical analysis and the results of numerical calculation show that the improved algorithm not only maintains all the advantages of the ex-tension & compaction scanning algorithm,but also has a higher computational efficiency.When the model is mesh more refined,the advantage of the improved algorithm is more significant in computational efficien-cy.【期刊名称】《湖南大学学报（自然科学版）》【年(卷),期】2014(000)001【总页数】6页(P39-44)【关键词】射线追踪;LTI算法;交叉扫描;扩张;收缩扫描法;改进算法【作者】卢江波;方志【作者单位】湖南大学土木工程学院，湖南长沙 410082;湖南大学土木工程学院，湖南长沙 410082【正文语种】中文【中图分类】TU317射线追踪技术在地震层析成像以及混凝土超声波射线层析成像等领域具有重要作用.目前常用射线追踪方法主要有有限差分解程函方程法［1－2］、最短路径法［3－4］以及LTI（Linear Travel－time Inter－polation）射线追踪算法［5－6］等.实验表明［5］，LTI算法在走时计算以及射线路径追踪上比其它方法（如有限差分解程函方程法）更为快速、精确.LTI射线追踪算法由Asakawa等人提出［5］，该算法以走时线性变化为前提，分两步进行.第一步，向前计算最小走时：先将模型划分为若干规则的单元，并将各单元边界划分为若干节段，然后，根据走时线性变化的假定以及走时最小原理（Fermat原理），求出从发射点到接收点的最小走时；第二步，根据得到的最小走时，反向追踪射线路径.但是该算法在向前计算走时时，没有考虑射线的逆向传播，不能追踪回波［9］，进而影响射线追踪的精度.针对这一问题，不少学者提出了改进算法，张东等人提出了循环计算LTI改进算法［7］，王浩全提出了交叉扫描LTI改进算法［8］，黄靓等人提出了扩张－收缩扫描改进算法［9－10］.在考虑射线的逆向传播时，文献［7－8］采用从发射点所在列向外逐列逆向扫描的方式，文献［9－10］则采用从边界列（行）向发射点所在列（行）进行逆向扫描的方式.在考虑射线的逆向传播上，文献［9－10］提出的改进算法更为合理，但是该算法存在计算效率低、收敛速度慢的问题.因为在LTI算法中，接收点要得到其最小走时，需其理论射线路径与模型单元边界所有相交节点都已得到其最小走时，这些相交节点的最小走时可能来自列扫描，也可能来自行扫描.根据本文作者对文献［5，9－10］算法步骤的理解以及对文献［9－10］中数值算例的研究，发现文献［9－10］不仅增加了从边界进行逆向传播的射线，即文献［9－10］中的收缩扫描过程，而且将文献［5］中既考虑邻列也考虑来自邻行入射射线的逐列扫描方式改为逐列（行）扫描只考虑邻列（行）入射射线的方式，这点可由文献［9］算例模型1的计算结果得出：若文献［9］的逐列（行）扫描过程既考虑邻列又考虑邻行的入射射线，则模型1经一次扩张扫描即得接收点的最小走时，此时相对误差应为0.247‰，而不应为按行列分开扫描方式经一次扩张扫描得到的计算结果3.3‰①.文献［9］采用行列分开扫描的方式使得文献［9－10］在能正确追踪回波的同时，接收点为获取其最小走时也进行了较多的无效扫描，降低了算法的计算效率.由以上分析可知，若采用行列交叉扫描的方式对扩张收缩扫描算法进行改进，将有效提高扩张收缩扫描算法的计算效率.同时，由于绕射波以及回波的存在，文献［5］算法步骤3对每列都进行水平边界节点最小走时搜索，其意义并不明确.此外，由于扩张收缩扫描算法相比文献［5］增加了收缩扫描过程，在逐列扫描过程中能够考虑上行或下行首波的最小走时，所以文献［5］中的逐行扫描可以省略.最后，由于存在收缩扫描过程，在计算竖向边界各节点最小走时时，若按文献［5］算法步骤5的计算方法，则存在着重复无效扫描.因此，也有必要对文献［5］的扩张扫描过程进行简化和改进.基于此，本文在前人研究的基础上提出了基于交叉扫描方式的扩张－收缩扫描新算法，该算法以扩张－收缩扫描算法为基础，结合文献［5］在逐列扫描过程中进行交叉扫描的思想，改进扩张－收缩扫描算法中逐列扫描的具体计算方法，提高了算法的计算效率，减少了迭代次数；同时，保留了原扩张收缩扫描算法的所有优点.1 LTI射线追踪算法的基本方程如图1所示，射线通过单元下边界的AB节段到达C节点，射线与AB的交点为D.A点及B点走时分别为TA和TB，节段AB长为L，单元慢度为s（速度的倒数），C点距A点的水平及竖向距离均为已知，分别为x，y.建立以A点为原点的局部坐标系，确定A，B，C，D点的局部坐标.现推导C点走时以及交点D 距A 点长度r的计算公式［5－6］.由线性追踪算法的基本假设可得D点的走时：根据D点的走时，结合单元慢度以及各点的局部坐标等条件，可得C点走时：将式（1）代入式（2）得根据费马原理，TC对r的一阶偏导数应当满足等于零的条件，即（设ΔT＝TB－TA）当L2s2－ΔT2 ＞0时，解方程（4）可得① 注：事实上，模型1接收点在按行列分开扫描的方式下，经扩张－收缩－再扩张后能够得到其最小走时，文献［9］给出的是一次扩张扫描结果，这应该是文献［9］作者的疏忽.若r≥0且r≤L，则若r＜0或r＞L，则计算r＝0和r＝L时的TC值，并取两者较小值作为最终TC值. 当L2s2－ΔT2≤0时，TC的计算方法与r＜0或r＞L时的情况一样.图1 经过节段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个节段；发射点S位于2号单元下边界的中点，接收点为a；线段为a节点的理论射线路径，其与模型各单元边界的交点分别为d，c和b.由理论射线路径可以看出，节点d，b只能通过扩张扫描的行扫描得到其最小走时，而节点a，c则需通过扩张扫描的列扫描得到其最小走时，所以接收点a要得到其最小走时，则模型的第3列必须先完成一次行列扫描，再对第4列进行一次列扫描.而扩张收缩扫描算法［9－10］在其扩张过程采用行列分开扫描的方式，所以a 节点在第一次扩张扫描中不能得到其最小走时.扩张－收缩扫描算法的第一次扩张扫描过程如图2所示（假定发射点所在列各节点的走时计算已经完成），其中图2（a）为扩张扫描的列扫描结果，图2（b）为完成所有列扫描后再进行行扫描的结果.从图2（b）可以看出a节点没有得到其最小走时.对于a节点，此次扩张扫描为无效扫描.同时，从图2（b）还可以看出，除a节点外还有其它节点也进行了无效扫描.易知，当模型网格划分得更多时，将有更多节点经历无效扫描.为了减少无效扫描，提高算法的计算效率.本文采用行列交叉扫描的方式对算法［9－10］进行改进，即在扩张（收缩）扫描过程，行列扫描不再分开进行，而是在其逐列（行）扫描时就同时进行行和列的扫描.仍以图2所示模型为例，采用改进后的算法对其进行一次扩张扫描，同样假定发射点所在列各节点的走时计算已经完成.图2 扩张收缩扫描算法的扩张扫描过程及结果Fig.2 Process and results of expansion scan of expansion－contraction scanning algorithm图3（a）为改进算法在扩张扫描时对第3列进行交叉扫描后的结果，可看出，b 和c节点在这次行列交叉扫描过程中都得到了最小走时，进而在第4列进行交叉扫描后，a节点能得到其最小走时（图3（b））.图3 改进算法的扩张扫描过程及结果Fig.3 Scanning process and results of improved algorithm对比图2（b）和图3（b），可以看出，改进算法进行的无效扫描更少，具有更高的计算效率，同时，由于改进算法仅改变了扫描顺序，两种算法的单次扩张扫描计算量不变.同时，对比文献［5］的一次扩张扫描过程，可以看出，本文提出的交叉扫描算法仅保留了逐列扫描过程中的行扫描，去掉了逐列扫描后附加的逐行扫描过程，并对文献［5］的逐列扫描过程进行了简化.需要说明的是，本文的交叉扫描与扩张收缩扫描在概念上是不同的，前者针对的是一次扩张收缩扫描中行列的计算次序，而后者针对的是整个单元网格的计算次序，一个为微观操作，另一个为宏观操作.下文中，除“扩张扫描”和“收缩扫描”中的“扫描”为宏观操作外，其余“扫描”均指微观操作.3 改进算法的基本步骤改进算法仍分两步进行，第一步向前计算走时，第二步，向后追踪射线路径，其中第二步与文献［9－10］相同，不赘述.改进算法的基本步骤如下：1）如图4所示，模型为3×3的网格，发射点S位于第2行第2列.以发射点所在的行和列为轴，将网格划分为四个象限，其中发射点所在行和列为各象限共有部分.图4（a）和图4（b）中的阴影部分分别为模型的第一和第二象限.图4 模型网格及象限划分示意图Fig.4 Model grid and quadrant division schematic diagram2）计算发射点S所在单元边界上各节点的走时，并记录次级源.假定此模型在各单元的边界上均只划分两个节段（图5），则各单元均有8个节点.根据单元节点以及发射点在整体坐标系中的坐标，可求得发射点S所在单元边界上各节点的走时，并将S点记为各节点的次级源.3）计算发射点S所在列所有节点的走时并记录次级源.首先计算图6中与发射点单元上边界相连的EJLG单元各节点的走时.以I点的走时计算为例（仅考虑通过下边界GE到达I节点的射线），易知，满足I 点最小走时要求的射线可能来自GE中的任一节段.此时，根据第1节给出的计算公式和计算方法分别计算出射线通过GF节段和FE节段时I点的最小走时，取两个最小走时的较小值作为I点的最小走时，并记录相应的次级源.以同样步骤求出EJLG单元其它节点的走时，完成该单元的计算.这种通过单元下边界节点走时计算其它节点走时的过程，称为向上扫描.其它通过上边界、左边界、右边界的情况分别称为向下扫描、向右扫描以及向左扫描.图5 改进算法的基本步骤2Fig.5 The basic steps 2of improved algorithm图6 改进算法的基本步骤3Fig.6 The basic steps 3of improved algorithm然后，以同样的方式向上逐个计算其它单元的节点走时，并记录相应的次级源.最后，从发射点单元开始，向下扫描，最终完成发射点S所在列所有节点走时计算以及次级源记录.4）扩张扫描.对步骤1划分的四个象限采用不同的扫描策略，具体如下.象限Ⅰ：从靠近发射点所在列的右侧列开始，向右逐列进行向右向上交叉扫描，直至模型的最后一列.以靠近发射点所在列的右侧列为例说明交叉扫描算法：如图7所示，先对此列单元进行向右扫描，然后再向上扫描，在扫描过程中，若扫描得到的节点走时比原节点走时小，则更新此节点的走时，并记录对应的次级源，否则，原节点走时以及相应次级源保持不变.象限Ⅱ，Ⅲ，Ⅳ的扫描也均从发射点所在列开始，向左（象限Ⅱ，Ⅲ）或向右（象限Ⅳ）逐列进行向左向上（象限Ⅱ）、向左向下（象限Ⅲ）以及向右向下（象限Ⅳ）交叉扫描，直至模型的边界.图7 改进算法的基本步骤4Fig.7 The basic steps 4of improved algorithm5）收缩扫描.四个象限的收缩扫描均从模型的竖向边界开始，向左（象限Ⅰ，Ⅳ）或向右（象限Ⅱ，Ⅲ）逐列进行向左向下（象限Ⅰ）、向左向上（象限Ⅳ）、向右向下（象限Ⅱ）或向右向上（象限Ⅲ）交叉扫描，直至发射点所在列.6）重复步骤4）和5），直到模型中所有节点的走时均不再改变为止.至此得到所有节点的最小走时以及相应的次级源.对比文献［5］的逐列扫描过程，本文省去了对每列都进行水平边界节点最小走时的搜索，原因如前言所述，在绕射波以及回波存在的情况下，搜索水平边界节点最小走时的意义并不明确.同时，为了确定逐列扫描过程中行扫描的顺序，本文参考了均匀介质模型各节点的理论射线路径，如发射点右上方单元各节点的最小走时均来自各自单元的左边界或下边界，据此将模型划分为四个象限，这也与仅改变文献［9］行列扫描顺序的结果一致.最后，本文按象限划分扫描区域进行逐列扫描的方式较文献［5］的逐列扫描过程更简洁，程序编制也更为容易.4 数值算例为了验证本文提出的改进算法在计算效率上的优越性，给出三个数值算例.考虑到改进算法仅改变了扫描顺序，不改变单次扩张收缩扫描计算量，故选取计算收敛所需迭代次数为对比参数.算例1 如图8，模型尺寸为3m×5m，各单元边长及速度分别为1m和3 500m／s，单元边界划分为2个节段，发射点S位于模型的下边界，且距模型左边界1.5m.采用扩张－收缩扫描算法［9－10］及本文改进算法进行计算.结果显示：扩张－收缩扫描算法需进行4次迭代才收敛，而改进算法只需要2次迭代即收敛.图8 算例1计算模型Fig.8 Calculation model of example 1算例2 将算例1模型细划为30×50的网格，单元边长为0.1m.如图9所示，发射点位置不变.单元边界划分为2个节段.计算结果显示：扩张收缩扫描算法需进行31次迭代，而改进算法仍只需2次迭代.图9 算例2计算模型Fig.9 Calculation model of example 2现将1 500号单元在两种算法下的扫描过程及结果列于表1，其中，由于算法收敛的条件，原算法及改进算法的最后一次迭代结果均与前一次的迭代结果相同.表1 1 500号单元在两种算法下的扫描过程及结果Tab.1 Scanning process and results of No.1 500under the two algorithms?算例1及算例2的计算结果表明：对于均质模型，改进算法相比于原算法具有更快的收敛速度，且当模型网格尺寸划分较细时，改进算法在计算效率上的优势更显著.算例3 如图10，模型尺寸为20m×20m，单元边长为1m；单元速度有两种，其中白色单元的速度为3 500m／s，黑色单元的速度为100m／s，单元边界划分为2个节段，发射点S和接收点R分别位于模型的左下角和下边界距发射点6m位置.计算结果显示：改进算法只需7次迭代即收敛，而原算法则需17次迭代；同时，接收点R在两种算法下的计算最小走时均为8.439 6ms，与理论最小走时相差0.19%.图10 算例3计算模型Fig.10 Calculation model of example 3算例3的计算结果表明：改进算法保留了扩张收缩扫描算法能正确处理射线逆向传播的优点，并且具有更快的收敛速度.5 结语理论分析以及数值算例表明，本文提出的基于交叉扫描方式的扩张－收缩扫描算法，不仅具有原扩张－收缩扫描算法的所有优点，而且在不增加单次扩张－扫描计算量的前提下，通过改变算法的扫描顺序，提高了算法的计算效率，加快了算法的收敛速度.特别当模型网格划分数较多时，改进算法在计算效率方面的优势更为显著.参考文献［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］刘洪，孟凡林，李幼铭.计算最小走时和射线路径的界面网全局方法［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－hao，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.（InChinese）。
直达波旅行时线性插值算法
直达波旅行时线性插值算法黄翼坚;朱光明;敦敏【摘要】文中首先简要介绍了初至波旅行时LTI算法,并分析了该法的计算结果包含非直达波旅行时的原因,由此提出通过对射线方向进行控制实现直达波旅行时计算的方法.直达波旅行时计算步骤与初至波旅行时LTI算法的计算步骤基本一样,只有两个不同之处:其一,旅行时赋初值不同,对震源点以外的节点,初至波旅行时LTI算法赋予大干该节点初至波旅行时的值,而直达波旅行时LTI算法则赋予大干该节点直达波旅行时的值;其二,在每次更新节点旅行时的过程中,一旦发现有波从某个低速地层传向某个高速地层,则在以后的计算中,判断射线是否是从该高速层指向该低速层,如果是则舍去该射线的旅行时.平层速度模型和具有起伏速度界面的多层模型算例说明,本文提出的直达波旅行时LTI迭代算法简单而有效,反射波射线追踪算例说明基于该算法可以进行大角度反射波射线追踪.【期刊名称】《石油地球物理勘探》【年(卷),期】2010(045)002【总页数】5页(P225-229)【关键词】初至波旅行时;直达波旅行时;线性插值;算法;费马原理;数值模拟【作者】黄翼坚;朱光明;敦敏【作者单位】长安大学地质工程与测绘学院,陕西西安,710054;长安大学地质工程与测绘学院,陕西西安,710054;东方地球物理公司物探技术研究中心,河北涿州,072751【正文语种】中文复杂速度模型中地震波旅行时的计算是基尔霍夫积分法叠前深度偏移及层析成像的一个关键步骤。
自20世纪80年代末以来,国内外许多学者研究了初至波旅行时的计算方法。
例如,Vidale[1]提出了初至波旅行时程函方程有限差分解法;Moser[2]提出了最短路径射线追踪及初至波旅行时计算方法;A sawaka和 Kawanaka[3]提出了初至波旅行时线性插值(L TI)算法;刘洪等 [4]提出了界面网全局旅行时计算技术;王华忠等[5]提出了动态规划法旅行时计算方法;赵连锋等[6]提出了有序波前重建法射线追踪技术;秦义龙等[7]利用单频双程波动方程计算初至波旅行时。
基于射线追踪技术计算地震定位中震源轨迹的改进方法
基于射线追踪技术计算地震定位中震源轨迹的改进方法赵爱华;丁志峰;白志明【期刊名称】《地球物理学报》【年(卷),期】2015(058)009【摘要】使用震源轨迹确定震源位置不仅稳健而且直观,但当介质复杂时震源轨迹难以给出解析解.基于最小走时树射线追踪技术计算震源轨迹的方法(以轨迹所在的残差场中残差最小的点(初始点)至残差较小的点(震源轨迹代表点)的射线路径表示震源轨迹)适用于复杂速度模型,但尚不能正确计算由多段组成的震源轨迹,同时兼顾计算轨迹的完整性和精细性较为困难,计算参数设置烦琐不适于大批量数据的自动处理.针对该方法存在的问题,本文对其进行了改进:(1)采用一种“削皮”算法选取震源轨迹所经过的模型单元的节点作为轨迹代表点;(2)将残差较小的区域作为震源轨迹计算区域(该区域依轨迹分布自适应地划分为若干个连通区域),从未计算的轨迹代表点中选取残差最小者作为射线路径初始点,利用最小走时树算法依次计算所有连通区域内的震源轨迹;(3)通过去掉较短的不再分叉的射线路径使震源轨迹更为精细.虚拟和真实事件的算例表明,改进方法有效克服了原方法的不足,可便捷地计算复杂速度模型中事件的震源轨迹,计算的轨迹精细且较完整.【总页数】14页(P3272-3285)【作者】赵爱华;丁志峰;白志明【作者单位】中国地震局地球物理研究所,北京 100081;中国地震局地球物理研究所,北京 100081;中国科学院地质与地球物理研究所,北京 100029【正文语种】中文【中图分类】P315【相关文献】1.射线追踪技术在城市环境场强预测计算中的应用 [J], 李超峰;焦培南;聂文强2.复杂介质地震定位中震源轨迹的计算 [J], 赵爱华;丁志峰;孙为国;王椿镛3.复杂介质地震定位中震源轨迹的计算 [J], 赵爱华;丁志峰;孙为国;王椿镛4.三维复杂速度模型中地震事件震源轨迹的计算 [J], 赵爱华5.一种基于插值技术高精度计算稀疏网格地震定位中震源轨迹的方法 [J], 赵爱华因版权原因，仅展示原文概要，查看原文内容请购买。