高精度地震波走时射线追踪方法与流程

本发明属于地震监测技术领域，涉及一种高精度地震波走时射线追踪方法。

背景技术：

传统的射线追踪方法包括试射法和弯线法，试射法是从震源出发，以一定的角度射出一条射线，射线根据snell定律进行传播，然后根据射线在地面上的落点和实际接收点的情况不断调整初射角度，最后由最靠近接收点的两条射线走时差值求出接收点的走时和射线路径。弯线法是先假设震源和接收点之间存在一条初始的射线路径，然后根据费马原理对射线路径进行改动，从而求出接收点的最小走时和射线路径。这两种方法时常面临陷于局部解的困扰，有时候不能求得全局最小解。

LTI算法是日本科学家Asakawa于1993年提出的一种射线追踪方法，该方法分为向前-向后计算2个基本步骤，并以线性插值为基本假设。向前计算的目的是为了得到各个网格节点上的最小走时，向后计算是在向前计算的基础上，由检波点根据走时最小的原理逐步追踪到炮点。

射线追踪的基本理论是，在高频近似的条件下，地震波的主要能量沿着射线路径传播。而射线追踪的目的就是在给定速度模型(m行n列的网格，每个网格内的速度为常数，不同网格间的速度可能不一样)、震源点位置和接收点位置时，求出一条射线的传播路径，使得从震源点到接收点的传播时间(走时)最小。

旅行时线性插值(Linear Traveltime Interpolation，LTI)法射线追踪具有较高的计算精度和效率，常用于二维速度模型的射线追踪。然而，该方法在考虑射线的来向时只考虑了一个半平面，即炮点相对于网格节点所在的半平面，没有考虑出现回折波的可能；对于地质情况复杂时，不能适应复杂速度模型的射线追踪。

下面通过建模来说明，给出一种速度模型，如图1所示，模型为1000m＊1000m的正方形，网格大小为10m＊10m，其中0～400m和600m～1000m深度上为1000m/s的低速体，中间400m～600m为1500m/s的高速体，该模型最先由Asakawa建立，用来验证射线追踪的准确性，本文称之为Asakawa模型。现在以(0，200)的位置为震源，(1000，0)、(1000,100)、……、(1000，900)的位置为接收点，得到的射线路径如图1中的曲线所示。需要注意的是射线3和射线4两条射线是经过高速体的折射后到达的。现在把模型顺时针旋转90°，震源和接收点相对于模型的位置不变，仍然用传统LTI算法进行射线追踪，结果如图2所示。对比图1和图2可以发现，射线3和射线4两条射线的路径不一样，图1中的两条射线是经过高速体的折射到达接收点的，而图2中的两条射线是直接由震源点到达接收点的。表1给出了图1和图2中射线走时的对比，可以看出，Asakawa模型旋转后射线3和射线4两条射线走时变大，说明传统LTI算法在面对复杂模型时计算精度不够。

表1：图1和图2中射线走时对比

技术实现要素：

为实现上述目的，本发明提供一种高精度地震波走时射线追踪方法，解决了现有技术中由于在考虑射线来向的问题上存在缺陷，算法的适应能力不强的问题。

本发明所采用的技术方案是，高精度地震波走时射线追踪方法，按照以下步骤进行：

步骤1，计算得到整个区域网格节点上的走时；

步骤2，除震源所在网格外，重新计算震源所在行网格上节点的走时，若比原来的小就替代；

步骤3，对于震源所在行之上的各行，先以网格下边界为插值线段，重新计算网格上边界上的2个节点的走时，若果比原来的小就替代；

步骤4，再以网格左边界为插值线段，重新计算网格右边界上2个节点的走时，若比原来的小就替代；

步骤5，再以网格右边界为插值线段，重新计算网格左边界上2个节点的走时，若比原来的小就替代；

步骤6，重复步骤3-5，直至到达空间内的第一行。

进一步的，所述步骤1中，计算得到整个区域网格节点上的走时具体按照以下步骤进行：

步骤a，计算震源所在网格上4个节点的走时；

步骤b，计算震源所在列网格上各节点的最小走时，包括以网格下边界为插值线段计算上方2个节点的最小走时，以及以网格上边界为插值线段计算下方2个节点的最小走时，以此类推，直到计算到第一行或最后一行；

步骤c，计算炮点所在列右侧一列网格节点走时，以网格左边界为插值线段计算网格右方2个节点的最小走时，对于重复计算的节点，若计算所得值比原来的小就替代；

步骤d，计算炮点所在列右侧一列网格节点走时，由下往上，以网格下边边界为插值线段，计算网格上方2个节点的最小走时，如果比原来的值小就替代；

步骤e，计算炮点所在列右侧一列网格节点走时，由上往下，以网格上边界为插值线段计算网格下边界2个节点的最小走时，如果比原来的值小就替代；

步骤f，重复步骤c-e，一直到最右边的一列网格；

对于炮点左侧网格节点最小走时的计算，采用和右侧相对应的方式得到。

本发明的有益效果是：在传统LTI算法的基础上进行了改进，考虑了网格节点射线从各个方向传来的可能，实验结果表明：改进后的LTI方法具有较强的适应性，网格节点走时计算精度得到提高，相应的射线路径也更加准确。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是Asakawa模型由传统LTI算法追踪出的射线路径图。

图2是Asakawa模型旋转以后用传统LTI算法追踪出的射线路径图。

图3是传统的旅行时线性插值(LTI)原理示意图。

图4是向前计算确定各节点最小走时的实现过程图。

图5是向后计算确定射线路径的过程图。

图6是本发明实施例算法的实现过程图。

图7是Asakawa模型旋转以后用本发明方法追踪出的射线路径。

具体实施方式

下面将结合本发明实施例中，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

先详细介绍传统的LTI算法的基本原理：

如图3所示，求某点P最小走时的问题可以作如下描述：假设A和B两点的坐标和走时已知，求经过线段AB上的某点C到达点P的最小走时。设A点为坐标原点，AB长度为L，AC长度为r，P点坐标为(x，y)，根据LTI的基本假设，C点的走时T_C可以由A和B两点的走时T_A和T_B线性插值得到，即：

则P点的走时为C点的走时加上CP线段所用的时间：

式中：ΔT＝T_B-T_A，s为空间内的慢度。对T_P求关于r的导数，并令其等于0

可得：

当时，

当时，

对一阶导数继续求导，可以得到T_P关于r的二阶导数，可以看出二阶导数在这两个解处大于0，可见这两个解都是函数T_P的极小值，而T_P又是关于r的初等函数，在区间0≤r≤L内连续，对比这2个解以及AB两端点的值，就可以得出T_P在0≤r≤L区间内的最小值。

下面介绍LTI算法实现步骤

LTI分为向前计算和向后计算两个基本步骤：向前计算是从炮点开始，由公式(4)求出空间内所有网格节点的最小走时；向后计算是从接收点开始，由公式(3)确定射线路径上与各个网格边界的交点。

向前计算的步骤为：

步骤1，计算震源所在网格上4个节点的走时，如图4(a)所示；

步骤2，计算震源所在列网格上各节点的最小走时，包括以网格下边界为插值线段计算上方2个节点的最小走时，以及以网格上边界为插值线段计算下方2个节点的最小走时，以此类推，直到计算到第一行或最后一行，如图4(b)所示；

步骤3，计算炮点所在列右侧一列网格节点走时，以网格左边界为插值线段计算网格右方2个节点的最小走时，对于重复计算的节点，若计算所得值比原来的小就替代，如图4(c)所示；

步骤4，计算炮点所在列右侧一列网格节点走时，由下往上，以网格下边边界为插值线段，计算网格上方2个节点的最小走时，如果比原来的值小就替代，如图4(d)所示；

步骤5，计算炮点所在列右侧一列网格节点走时，由上往下，以网格上边界为插值线段计算网格下边界2个节点的最小走时，如果比原来的值小就替代，如图4(e)所示；

步骤6，重复步骤3-5，一直到最右边的一列网格，如图4(f)。

对于炮点左侧网格节点最小走时的计算，可采用和右侧相似的方式得到。

向后计算的过程为：

步骤1，根据公式t＝t_向前+sd计算经接收点所在网格上4个节点到达接收点的走时，找出其中最小的如图5(a)所示，其中t_向前为向前计算得到的网格上节点的走时，s为网格的慢度，d为接收点到节点距离；

步骤2，根据步骤1得到的网格节点，用公式(3)和公式(4)计算到达接收点的最小走时以及对应的插值点；

步骤3，将步骤2中求得的最小值对应的插值点作为新的接收点，重复步骤1和步骤2，得到下一个最小走时对应的插值点，如图5(b)所示；

步骤4，重复步骤3，直到求得的插值点到达炮点所在网格的边界，如图5(c)所示；

步骤4，从接收点开始，依次连接最小走时所对应的插值点，最后到达震源点，这样就得到了由震源点到接收点最小走时的路径。

因此，结合背景技术所说，传统LTI算法由于在考虑射线来向的问题上存在缺陷，导致算法的适应能力不强，现在针对这一点进行改进，提出了高精度LTI算法，实现过程如图6所示，具体按照以下步骤进行：

步骤1，先按照传统的方法计算得到整个区域网格节点上的走时；

步骤2，除震源所在网格外，重新计算震源所在行网格上节点的走时，若比原来的小就替代，如图6(a)所示；

步骤3，对于震源所在行之上的各行，先以网格下边界为插值线段，重新计算网格上边界上的2个节点的走时，若果比原来的小就替代，如图6(b)所示；

步骤4，再以网格左边界为插值线段，重新计算网格右边界上2个节点的走时，若比原来的小就替代，如图6(c)所示；

步骤5，再以网格右边界为插值线段，重新计算网格左边界上2个节点的走时，若比原来的小就替代，如图6(d)所示；

步骤6，重复步骤3-5，直至到达空间内的第一行。

对于炮点所在网格下方各行节点走时的重新计算，可采用相似的方式。

当采用本发明方法之后，重新按照图2的模型和震源-接收点关系进行射线追踪，得到的结果如图7所示，由图7可以看出，射线3和射线4两条射线也是经由高速层折射到达的接收点，这与图1中的结果相同，同时对比表2中的走时信息，改进后的走时比原来的更小，也与图1中对应接收点的走时相同，说明改进后的LTI算法比传统LTI算法得到的走时精度更高，追踪出来的射线路径也更准确。

表2：图2和图7中射线走时对比

以上所述仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内所作的任何修改、等同替换、改进等，均包含在本发明的保护范围内。