CN117607957B - 基于等效慢度的快速推进法的地震波走时求解方法及系统 - Google Patents

Abstract

本发明公开了基于等效慢度的快速推进法的地震波走时求解方法及系统，该方法包括：获取在初始离散网格内的震源点位置，计算震源点位置与初始离散网格内最近的网格点的偏移向量，根据偏移向量对初始离散网格进行移动，得到目标网格；建立在目标网格内空间点的走时、震源点到空间点的距离和等效慢度的关系方程，并根据关系方程得到基于等效慢度的程函方程；对程函方程进行离散，得到求解空间点对应的等效慢度的目标方程组；根据目标方程组和关系方程，计算目标网格内多个空间点各自对应的走时；根据多个地震波走时各自对应的等效慢度，得到在初始离散网格内多个网格点各自对应的当前地震波走时。本发明提高了地震波走时的计算精度。

Description

基于等效慢度的快速推进法的地震波走时求解方法及系统

技术领域

本发明涉及地球物理勘查技术领域，尤其涉及一种基于等效慢度的快速推进法的地震波走时求解方法、系统、终端及存储介质。

背景技术

地震波的走时是指地震波从震源激发后传播至地下介质中任意一点所需要的时间。地震波的走时与地下介质的速度分布密切相关，因此，走时计算在地球物理勘查技术领域的层析成像、深度偏移、地震定位等方面具有非常重要的应用。

快速推进法是一种快速的有限差分求解地震波走时的方法。然而，由于波前面在震源附近曲率较大，利用快速推进法求解震源附近的走时，往往存在较大的误差，并且这种误差会逐渐累积，从而影响整体的走时计算精度，另外，震源点的坐标有时并未与离散的网格点完全重合，这时采用传统方法计算的走时场，会产生系统误差，影响最终的计算结果。

因此，现有技术还有待于改进和发展。

发明内容

本发明的主要目的在于提供一种基于等效慢度的快速推进法的地震波走时求解方法及系统，旨在解决现有技术中采用快速推进法由于震源点的坐标与离散网格点未完全重合而产生系统误差，以及走时场在震源附近曲率较大，导致计算的地震波走时的精度较低的问题。

本申请实施例第一方面提供一种基于等效慢度的快速推进法的地震波走时求解方法，所述基于等效慢度的快速推进法的地震波走时求解方法包括如下步骤： 获取在初始离散网格内的震源点位置，计算所述震源点位置与所述初始离散网格内最近的网格点的偏移向量，根据所述偏移向量对所述初始离散网格进行移动，得到目标网格；其中，所述震源点位置在所述目标网格内的空间点上；建立在目标网格内空间点的走时、等效慢度以及震源点到所述空间点的距离的关系方程，并根据所述关系方程得到基于等效慢度的程函方程；对所述程函方程进行离散，得到用于求解所述空间点对应的所述等效慢度的目标方程组；根据所述目标方程组和所述关系方程，计算所述目标网格内多个所述空间点各自对应的地震波走时。

在一些实现方式中，所述关系方程中的所述等效慢度等于在所述目标网格内所述震源点到所述空间点的走时与所述震源点到所述空间点的距离的商；所述根据所述关系方程得到基于等效慢度的程函方程，具体包括：获取用于计算地震波走时的基本方程；将所述关系方程中的等效慢度和距离替代所述基本方程中的走时，得到基于等效慢度的程函方程。

在一些实现方式中，所述目标方程组包括离散方程和差分算子方程组；所述对所述程函方程进行离散，得到用于求解所述空间点对应的所述等效慢度的目标方程组，具体包括：采用迎风格式对所述程函方程进行求解，得到离散方程；将所述关系方程中的所述等效慢度和所述距离代入所述离散方程中的差分算子，得到差分算子方程组。

在一些实现方式中，所述差分算子方程组包括一阶精度方程和二阶精度方程；所述将所述关系方程中的所述等效慢度和所述距离代入所述离散方程中的差分算子，得到差分算子方程组，具体包括：确定所述空间点的走时的精度阶位，其中，所述精度阶位包括一阶精度和二阶精度；将所述关系方程中的所述等效慢度和所述距离代入所述离散方程中所述一阶精度对应的差分算子，得到一阶精度方程；将所述关系方程中的所述等效慢度和所述距离代入所述离散方程中所述二阶精度对应的差分算子，得到二阶精度方程。

在一些实现方式中，所述目标网格内多个所述空间点各自对应的走时包括初始走时和中间走时；所述根据所述目标方程组和所述关系方程，计算所述目标网格内多个所述空间点各自对应的地震波走时，具体包括：获取所述目标网格中的震源点位置；根据所述震源点位置、所述目标方程组和所述关系方程，计算所述震源点位置对应的初始走时和所述震源点位置周围的多个所述空间点各自对应的中间走时；其中，所述震源点位置周围的多个所述空间点为筛选空间点，所述筛选空间点中已剔除过作为所述震源点位置的所述空间点；根据多个所述空间点各自对应的所述中间走时，确定多个所述空间点中的目标空间点，并将所述目标空间点作为新的震源点位置计算新的中间走时，直至所有所述空间点各自对应的走时均已计算。

在一些实现方式中，所述根据所述震源点位置、所述目标方程组和所述关系方程，计算所述震源点位置对应的初始走时和所述震源点位置周围的多个所述空间点各自对应的中间走时，具体包括：根据所述目标方程组，计算所述震源点位置对应的初始等效慢度，并根据所述初始等效慢度和所述关系方程得到所述初始走时；根据所述目标方程组，计算所述震源点位置周围的多个所述空间点各自对应的中间等效慢度，并根据各所述中间等效慢度和所述关系方程得到多个所述空间点各自对应的中间走时。

在一些实现方式中，所述根据多个所述空间点各自对应的所述中间走时，确定多个所述空间点中的目标空间点，并将所述目标空间点作为新的震源点位置计算新的中间走时，直至所有所述空间点各自对应的走时均已计算，具体包括：根据多个所述空间点各自对应的所述中间走时，确定多个所述中间走时中的最小值走时；其中，多个所述中间走时为已计算过的全部走时中剔除了作为过所述震源点位置对应的走时；将所述最小值走时对应的所述空间点作为目标空间点，并将所述目标空间点作为新的震源点位置； 根据所述新的震源点位置、所述目标方程组和所述关系方程，计算所述新的震源点位置周围的未计算过的所述空间点对应的新的中间走时，直至所有所述空间点均已计算过走时。

本申请实施例第二方面还提供一种基于等效慢度的快速推进法的地震波走时求解系统，其中，所述基于等效慢度的快速推进法的地震波走时求解系统包括：网格确定模块，用于获取在初始离散网格内的震源点位置，计算所述震源点位置与所述初始离散网格内最近的网格点的偏移向量，根据所述偏移向量对所述初始离散网格进行移动，得到目标网格；其中，所述震源点位置在所述目标网格内的空间点上；方程建立模块，用于建立在目标网格内空间点的走时、等效慢度以及震源点到所述空间点的距离的关系方程，并根据所述关系方程得到基于等效慢度的程函方程；方程离散模块，用于对所述程函方程进行离散，得到用于求解所述空间点对应的所述等效慢度的目标方程组；走时计算模块，用于根据所述目标方程组和所述关系方程，计算所述目标网格内多个所述空间点各自对应的地震波走时；当前走时确定模块，用于根据多个所述地震波走时各自对应的等效慢度，得到在所述初始离散网格内多个网格点各自对应的当前地震波走时。

本申请实施例第三方面还提供一种终端，其中，所述终端包括：存储器、处理器及存储在所述存储器上并可在所述处理器上运行的基于等效慢度的快速推进法的地震波走时求解程序，所述基于等效慢度的快速推进法的地震波走时求解程序被所述处理器执行时实现如上所述的基于等效慢度的快速推进法的地震波走时求解方法的步骤。

本申请实施例第四方面还提供一种计算机可读存储介质，其中，所述计算机可读存储介质存储有基于等效慢度的快速推进法的地震波走时求解程序，所述基于等效慢度的快速推进法的地震波走时求解程序被处理器执行时实现如上所述的基于等效慢度的快速推进法的地震波走时求解方法的步骤。

有益效果：本发明提供了基于等效慢度的快速推进法的地震波走时求解方法及系统，该方法通过将离散网格移动形成目标网格的空间的与震源点坐标重合，建立某一点的走时、点到震源点的直线距离和等效慢度的关系，并根据该关系将以走时的方程转换成基于等效慢度的方程，等效慢度具有更小的曲率半径，从而在进行方程计算求解走时阶段，使得计算得到的走时具有较高的精度，最后对应计算初始离散网格点上的走时，避免由于震源点未完全与离散网格重合所产生的系统误差，进而提高地震波走时的计算精度。

附图说明

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

图1是本申请基于等效慢度的快速推进法的地震波走时求解方法的较佳实施例的流程图；

图2是本申请基于等效慢度的快速推进法的地震波走时求解方法的较佳实施例中快速推进法所采用的窄带技术示意图；

图3是本申请基于等效慢度的快速推进法的地震波走时求解方法的较佳实施例中快速推进法地震波走时计算的具体实施步骤流程示意图；

图4是本申请基于等效慢度的快速推进法的地震波走时求解方法的较佳实施例中二维常速度梯度模型；

图5是本申请基于等效慢度的快速推进法的地震波走时求解系统的较佳实施例的结构示意图；

图6为本申请终端的较佳实施例的结构示意图。

附图标记说明：

10、基于等效慢度的快速推进法的地震波走时求解系统；100、网格确定模块；200、方程建立模块；300、方程离散模块；400、走时计算模块；500、当前走时确定模块。

具体实施方式

为使本发明的目的、技术方案及效果更加清楚、明确，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，所描述的实施例仅为本发明的可能的技术实现，并非全部实现可能。基于本发明中的实施例，本领域技术人员完全可以结合本发明的实施例，在没有进行创造性劳动的情况下得到其他实施例，而这些实施例也在本发明的保护范围之内。

首先，对本发明实施例中涉及的名词进行介绍：地震波走时，是指地震波从震源点传递到观测点所经过的时间。由于地震纵波和横波的传播速度以及所经之处的物质组成各不相同，因此地震波的走时会受到这些因素影响而有所差异。在地震学领域，研究地震波的运动学信息，特别是计算地震波的走时，是至关重要的。

现有技术中，快速推进法是一种快速的有限差分求解地震波走时的方法，在求解地震波走时时，由于波前面在震源附近曲率较大，采用一阶或二阶差分格式进行离散往往存在较大误差，然而，采用高阶差分格式（三阶或以上）很难提高精度，这主要是由于快速推进法采用的特殊求解方式（迎风格式及窄带技术）所决定的，快速推进法通常采用二阶差分格式，当所需节点不足时（如震源附近），差分格式会自动降为一阶格式进行离散求解。

为了缓解快速推进法在震源附近模拟精度不足的问题，通常会在震源附近加密网格，也就是通过自动加密技术（震源附近采用小网格，远离震源处采用大网格）提高走时的求解精度，但该方法对精度的提升有限，并且这种方式增加了实现的困难，降低了计算效率。另外，震源点的坐标有时并未与离散的网格点完全重合，这时采用传统方法计算的走时场，会产生系统误差，影响最终的计算结果。

下面参考附图描述本申请实施例的基于等效慢度的快速推进法的地震波走时求解方法及系统。针对上述提到的相关技术中采用快速推进法由于震源点的坐标与离散网格点未完全重合而产生系统误差，以及走时场在震源附近曲率较大，导致计算的地震波走时的精度较低的问题，本申请提供了一种基于等效慢度的快速推进法的地震波走时求解方法，在该方法中，通过将离散网格移动形成目标网格的空间的与震源点坐标重合，建立某一点的走时、点到震源点的直线距离和等效慢度的关系，并根据该关系将以走时的方程转换成基于等效慢度的方程，等效慢度具有更小的曲率半径，从而在进行方程计算求解走时阶段，使得计算得到的走时具有较高的精度，最后通过插值计算初始离散网格点上的走时，避免由于震源点未完全与离散网格重合所产生的系统误差，进而提高地震波走时的计算精度。由此，解决了相关技术中采用快速推进法在震源附近计算的地震波走时由于震源点的坐标与离散网格点未完全重合而产生系统误差，以及走时场在震源附近曲率较大，导致计算的地震波走时的精度较低的技术问题。

本发明的发明构思为：由于在实际非均匀介质中，射线的传播路径是弯曲的，而走时是对慢度的曲线积分，本发明采用“以直代曲”的思想，利用点到源点的直线路径来代替真实的曲线路径，等效慢度代替真实慢度，因此走时的计算可以由慢度沿真实曲线路径的积分转换为点到源点的直线距离与等效慢度的乘积，即，等效慢度为走时与点到震源点直线距离的商。该处理方式是将走时场的复杂度或者曲率信息部分转移到具有解析解的距离项上，使得等效慢度具有更小的曲率半径，并且等效慢度场的分布更加光滑。因此，在使用一阶或二阶差分格式对等效慢度进行离散求解，在获得等效慢度以后，利用距离乘以等效慢度可转换得到该点的走时，当震源点并未完全与离散网格重合时，计算会产生系统误差，需要对网格进行移动产生新的离散网格，从而使得震源点完全与新的离散网格的某一离散点完全重合，在计算出所有新网格的走时后，通过插值计算原始网格点上的走时，避免由于震源点未完全与离散网格重合所产生的系统误差，从而使计算的地震波走时能够获得更高的精度。

本发明提高了地震波走时求解的精度（较于现有技术相同网格），并且提高了地震波走时求解的效率（较于现有技术相同精度）。

下面以具体地实施例对本发明的技术方案进行详细说明。下面这几个具体的实施例可以相互结合，对于相同或相似的概念或过程可能在某些实施例不再赘述。

本发明较佳实施例所述的基于等效慢度的快速推进法的地震波走时求解方法，如图1所示，所述基于等效慢度的快速推进法的地震波走时求解方法包括以下步骤：

在步骤S101中，获取在初始离散网格内的震源点位置，计算所述震源点位置与所述初始离散网格内最近的网格点的偏移向量，根据所述偏移向量对所述初始离散网格进行移动，得到目标网格；其中，所述震源点位置在所述目标网格内的空间点上。

具体的，震源点位置为震源点的空间坐标，对原始网格按偏移向量移动获得新的离散网格，使震源与新的网格的某一网格点完全重合，即为目标网络（用于计算的新的计算网络），空间点为目标网络内的网格点，该空间点与初始离散网格的网格点不重合。

需要说明的是，在进行走时计算时，经常遇到震源点并未与离散网格点完全重合，此时，利用传统方法或等效慢度方法都会产生系统误差，本申请实施例在原始离散网络获得震源点位置坐标以后，求出震源点最近的网格点的偏移向量，然后把原始网格按偏移向量移动，产生一个新的网格，通过震源点与新网格的某一网格点完全重合，从而可以消除走时计算的系统误差，以提高走时的计算精度。

在步骤S102中，建立在目标网格内空间点的走时、等效慢度以及震源点到所述空间点的距离的关系方程，并根据所述关系方程得到基于等效慢度的程函方程。

在一种实现方式中，所述关系方程中的所述等效慢度为在所述目标网格内所述空间点的走时与所述震源点到所述空间点的距离的商；获取用于计算地震波走时的基本方程；将所述关系方程中的等效慢度和距离替代所述基本方程中的走时，得到基于等效慢度的程函方程。

具体的，等效慢度为空间一点的走时与该点到震源点直线距离的商，即关系方程为：

；（1）

其中，是震源点到空间任意一点（即空间点）的直线距离，/>是点处的等效慢度，/>为走时场，/>为空间点的坐标，/>分别为三维坐标中对应三个方向的位置；

求解地震波走时的基本方程为：

；（2）

其中，为慢度，/>为梯度符号；

把方程（1）带入到方程（2）中消除走时，从而可以得到由等效慢度所表征的程函方程为：

；（3）

需要说明的是，根据关系方程将基本方程转换为基于等效慢度的程函方程，在程函方程（3）中，和/>具有解析，因此，这两项的求取是准确的，主要误差来自于/>的求解，通常情况下，等效慢度是非常光滑的，尤其是在震源附近，例如，在均匀介质中，等效慢度是常数，本质上，等效慢度方法将走时场的复杂度或曲率部分转移到具有解析解的距离项上，在震源附近，/>具有较大的曲率，但可以通过解析方法求解，因此不会损失精度，而等效慢度本身在震源附近具有相对光滑的分布，因此，当采用一阶或二阶差分格式时可获得更高的精度。

在步骤S103中，对所述程函方程进行离散，得到用于求解所述空间点对应的所述等效慢度的目标方程组。

在一种实现方式中，所述目标方程组包括离散方程和差分算子方程组；采用迎风格式对所述程函方程进行求解，得到离散方程；将所述关系方程中的所述等效慢度和所述距离代入所述离散方程中的差分算子，得到差分算子方程组。

具体的，采用迎风格式对方程（3）进行求解，其离散形式表示的离散方程为：

；（4）

其中，为求取最大值函数，/>为网格点处/>的慢度值，E为等效慢度，、/>、/>、/>、/>和/>为空间导数算子，下标/>，/>，/>分别代表/>，/>，方向的索引，上标/>，/>，/>分别代表沿着三个方向的导数，上标中的正负号分别代表前向差分和后向差分。

在一种实现方式中，所述差分算子方程组包括一阶精度方程和二阶精度方程；确定所述空间点的走时的精度阶位，其中，所述精度阶位包括一阶精度和二阶精度；将所述关系方程中的所述等效慢度和所述距离代入所述离散方程中所述一阶精度对应的差分算子，得到一阶精度方程；将所述关系方程中的所述等效慢度和所述距离代入所述离散方程中所述二阶精度对应的差分算子，得到二阶精度方程。

具体的，对等效慢度的一阶精度方程包括一阶精度的前向差分算子和一阶精度的后向差分算子；一阶精度的前向差分算子可以表示为：

；（5）

一阶精度的后向差分算子可以表示为：

；（6）

对等效慢度的二阶精度方程包括二阶精度的前向差分算子和二阶精度的后向差分算子；二阶精度的前向差分算子可以表示为：

；（7）

二阶精度的后向差分算子可以表示为：

；（8）

可以理解的是，为/>方向的采样间隔，/>是点/>处的等效慢度，是点/>处的等效慢度，/>是点/>处的等效慢度，/>是点处的等效慢度，/>是点/>处的等效慢度。参见图2，下标/>，/>，/>分别表示震源点/>，/>，/>右侧的空间点（网格点），下标/>，/>，/>分别表示震源点/>，/>，/>左侧的空间点（网格点），下标/>，/>，/>分别表示震源点/>，/>，/>右侧两次的空间点（网格点），下标/>，/>，/>分别表示震源点/>，/>，/>左侧两次的空间点（网格点）。在未示出的每个精度的差分算子公式中，下标/>，/>，/>分别表示震源点/>，/>，/>上侧的空间点（网格点），下标/>，/>，/>分别表示震源点/>，/>，/>下侧的空间点（网格点）；下标/>，/>，/>分别表示震源点/>，/>，/>前侧的空间点（网格点），下标/>，/>，/>分别表示震源点/>，/>，/>后侧的空间点（网格点）；同理/>方向的采样间隔为/>，/>方向的采样间隔为/>。

求解方程（4）时，并非完全稳定的，也即并非都能求出该点的等效慢度，然后再转换成走时。而对方程（3）的离散求解是无条件稳定的，从而可以转换到传统方法，利用传统方法求得该点的走时，然后再转换成等效慢度，而继续求解别的点的等效慢度时，仍通过方程（4）求解，仅在方程（4）无解时转换采用方程（3）离散求解。因此，当求解方程（4）无解时，转而对方程（3）离散求解进行替代，其离散形式的方程为：

，（9）

具体的，对走时的一阶精度方程包括一阶精度的前向差分算子和一阶后向差分算子可以表示为：一阶精度的前向差分算子可以表示为：

；（10）

一阶精度的后向差分算子可以表示为

；（11）

值得说明的是，本发明实施例在震源附近空间点不足以及等效慢度方法在某些点无解需要转到传统方法求解的情况下选用一阶精度进行计算，在其他情况下均采用二阶精度计算，从而提高地震波走时的计算精度。

需要说明的是，方程（4）通过迎风格式利用相邻网格点的已知走时信息计算新的走时，需要保证节点（空间点）的更新顺序是由小到大的（走时）。因此，快速推进算法采用了窄带技术，如图2所示。

本发明通过离散等效慢度代替传统的离散走时，计算出等效慢度以后，再转换成走时，通过这种方式避免了震源奇异性的问题，能够大幅提高震源附近走时计算的精度。

在步骤S104中，根据所述目标方程组和所述关系方程，计算所述目标网格内多个所述空间点各自对应的地震波走时。

在一种实现方式中，所述目标网格内多个所述空间点各自对应的走时包括初始走时和中间走时；获取所述目标网格中的震源点位置；根据所述震源点位置、所述目标方程组和所述关系方程，计算所述震源点位置对应的初始走时和所述震源点位置周围的多个所述空间点各自对应的中间走时；其中，所述震源点位置周围的多个所述空间点（即窄带中的空间点）为筛选空间点，所述筛选空间点中已剔除过作为所述震源点位置的所述空间点；根据多个所述空间点各自对应的所述中间走时，确定多个所述空间点中的目标空间点，并将所述目标空间点作为新的震源点位置计算新的中间走时，直至所有所述空间点各自对应的走时均已计算。

可以理解的是，初始走时为震源点位置处对应空间点的走时（一个），中间走时为目标网格中处除了震源点位置处对应空间点以外的空间点的走时（多个）。

具体的，根据所述目标方程组，计算所述震源点位置对应的初始等效慢度，并根据所述初始等效慢度和所述关系方程得到所述初始走时；根据所述目标方程组，计算所述震源点位置周围的多个所述空间点各自对应的中间等效慢度，并根据各所述中间等效慢度和所述关系方程得到多个所述空间点各自对应的中间走时。

进一步，把相应的差分算子带入到方程（5）-（8）的一个带入到方程（4）中，即可求解在空间点处的等效慢度。当方程（4）求解无解时，需要把相应的差分算子带入到方程（10）-（11）的一个到方程（9）中，即可求解该空间点处的走时，然后计算该点的等效慢度。可以理解的是，空间点可以是三维下的点，也可以是二维下的点。

需要说明的是，如图2所示，快速推进算法采用了窄带技术以保证节点（空间点）的走时更新顺序是由小到大的；窄带技术是把所有需要计算的网格节点分为三类：激活点（Alive points）、近点（Close points）和远点（Far points），激活点是已经完成计算并固定的节点，近点是已经计算但未固定的节点，远点是尚未计算的节点，激活点构成了上风（Upwind）区域，近点构成了窄带（Narrow band），远点则构成了下风（Downwind）区域，也就是说，当前的各空间点为目标网格内处于上风区域的激活点、处于窄带的近点及处于下风区域的远点中的一种。

具体的，根据多个所述空间点各自对应的所述中间走时，确定多个所述中间走时中的最小值走时；其中，多个所述中间走时为已计算过的全部走时中剔除作为过所述震源点位置对应的走时；将所述最小值走时对应的所述空间点作为目标空间点，并将所述目标空间点作为新的震源点位置；根据所述新的震源点位置（即子震源点）、所述目标方程组和所述关系方程，计算所述新的震源点位置周围的未计算过的所述空间点对应的新的中间走时，重复上述步骤（重复子震源点周围空间点的等效慢度和走时计算），直至所有所述空间点均已计算过走时。

进一步，快速推进算法从已知值的上风区域开始，系统地向下风方向构建走时。参见图3，该过程可以概括为：首先计算并确定震源点的走时，然后计算其周围（前后、上下、左右）网格节点的走时，并将它们加入到一个窄带（即已计算走时的空间点的集合）中；接着从窄带中选择走时最小的节点（空间点），将其固定为一个新的子震源，并计算其周围节点的走时；如果节点已在窄带中，则更新其走时并保留较小的走时；如果节点不在窄带中，则将其加入窄带。然后从更新后的窄带中再选择走时最小的节点，固定其走时，并计算其周围节点的走时。重复这个过程，直到所有网格节点的走时都被确定。

在步骤S105中，根据多个所述地震波走时各自对应的等效慢度，得到在所述初始离散网格内多个网格点各自对应的当前地震波走时。

可以理解的是，目标网络（即新的计算网格）的所有的等效慢度以及地震波走时求取完成以后，通过双线性插值（二维）或三线性插值（三维）得到初始离散网络内各个网格点（即原始网格点）对应的当前等效慢度，并根据当前等效慢度和距离计算得到各个网格点各自对应的当前地震波走时。

可以理解的是，双线性插值是一种在两个方向分别进行一次线性插值的方法，具体来说，双线性插值利用原图像中的四个点来计算新图像中的某一个点。三线性插值是一种在三个方向分别进行一次线性插值的方法，具体来说，三线性插值利用原图像中的八个点来计算新图像中的某一个点。

本发明可以采用更大的网格进行等效慢度的计算，然后再转换成走时，可大幅提高走时计算的效率。

下面结合一具体应用场景对本发明的具体实现进行描述，在该具体应用场景中，基于等效慢度的快速推进算法从具有已知值的上风区域开始，系统地向下风方向构建走时：

步骤S1，获取震源点的空间坐标，计算震源点与最近的网格点的偏移向量，然后对原始网格按偏移向量进行移动，从而产生新的计算网格。

步骤S2，初始化新网格中所有网格点为远点，走时初始化为一个大值（例如10e⁶）；

步骤S3，初始化震源点的等效慢度为该点的慢度值，距离为零，通过将距离乘以等效慢度，可以得该点的走时值，震源点被加入到窄带中；

步骤S4，利用最小堆排序获取窄带内走时最小的点，取出并标记为激活点；

步骤S5，将该激活点作为新的子震源，根据方程（5）-（8）及（4）计算其周围（前后、上下、左右）未激活节点的等效慢度，然后，乘以距离以获得周围未激活点的走时；

步骤S6，如果节点已在窄带中，则更新其走时，并保留较小的走时，同时，也更新该点的等效慢度。如果节点不在窄带中，则将其加入到窄带中；

步骤S7，重复执行步骤S4-S6，直到所有新网格节点的等效慢度和走时都被确定。

步骤S8，利用新网格点的等效慢度对原始网格点进行插值，获得原始网格点的等效慢度，然后转换成走时。

为了验证基于等效慢度的快速推进法求解地震波走时方法的有效性。本发明还构建了一个二维常速度梯度模型，该二维常速度梯度模型的走时场和等效慢度场具有解析解，从而可以用于验证数值模拟方法的精度和效率。

在测试时，设置二维常速度梯度模型大小为10 km x 6 km，震源位于水平位置5km、深度200 m处。在震源位置，速度为2000 m/s。设置速度梯度为0.4。图4展示了该二维常速度梯度模型、走时以及等效慢度的解析解，具体的，图4中的a）为在不同深度下，地点对应速度场的解析解；图4中的b）为在不同深度下，地点对应走时场的理论值；图4中的c）为在不同深度下，地点对应等效慢度的理论值。可以看出，走时场在震源附近具有较大的曲率，而等效慢度的曲率较小，等效慢度的变化也相对平缓，因此，利用差分算子进行离散时，等效慢度法具有较高的精度。

本发明一实施例选用200米x 200米的网格间距进行走时场的计算，传统方法计算出的走时与理论值的最大绝对误差为31.4 ms，而采用等效慢度法计算出的走时场与理论值的最大绝对误差为1.0 ms，计算的精度提升了一个数量级。本发明另一实施例选用100米x 100米的网格间距进行走时场的计算，传统方法计算出的走时与理论值的最大绝对误差为17.1 ms，而采用等效慢度法计算出的走时场与理论值的最大绝对误差为0.31 ms。本发明另一实施例选用50米x 50米的网格间距进行走时场的计算，传统方法计算出的走时与理论值的最大绝对误差为8.4 ms，而采用等效慢度法计算出的走时场与理论值的最大绝对误差为0.096 ms。本发明另一实施例选用25米x 25米的网格间距进行走时场的计算，传统方法计算出的走时与理论值的最大绝对误差为4.2 ms，而采用等效慢度法计算出的走时场与理论值的最大绝对误差为0.032 ms。

通过上述结果，可以观测出，传统方法在网格间距较大时，存在较大的误差，该误差主要由于震源奇异性引起的。当网格间距减半时，传统方法的误差大约降低为原来的误差的一半。对于等效慢度法，即使采用较大的网格，仍能保持较高的精度。并且当网格间距减半时，等效慢度法的误差大约降低为原来误差的1/3，因此，等效慢度法具有更快的收敛率。

因此，本发明可以在相同的网格间距下，采用等效慢度法获得更高的精度。或者在相同精度下，通过等效慢度法，可以采用更大的网格进行计算，从而获得更高的效率。

其次参照附图描述根据本申请实施例提出的基于等效慢度的快速推进法的地震波走时求解系统。

图5是本申请实施例的基于等效慢度的快速推进法的地震波走时求解系统的方框示意图。

如图5所示，该基于等效慢度的快速推进法的地震波走时求解系统10包括：网格确定模块100、方程建立模块200、方程离散模块300、走时计算模块400和当前走时确定模块500。

具体地，网格确定模块100，用于获取在初始离散网格内的震源点位置，计算所述震源点位置与所述初始离散网格内最近的网格点的偏移向量，根据所述偏移向量对所述初始离散网格进行移动，得到目标网格；其中，所述震源点位置在所述目标网格内的空间点上；

方程建立模块200，用于建立在目标网格内空间点的走时、震源点到所述空间点的距离和等效慢度的关系方程，并根据所述关系方程得到基于等效慢度的程函方程。

方程离散模块300，用于对所述程函方程进行离散，得到求解所述空间点对应的所述等效慢度的目标方程组

走时计算模块400，用于根据所述目标方程组和所述关系方程，计算所述目标网格内多个所述空间点各自对应的走时。

当前走时确定模块500，用于根据多个所述地震波走时各自对应的等效慢度，得到在所述初始离散网格内多个网格点各自对应的当前地震波走时.

需要说明的是，前述对基于等效慢度的快速推进法的地震波走时求解方法实施例的解释说明也适用于该实施例的基于等效慢度的快速推进法的地震波走时求解系统，此处不再赘述。

根据本申请实施例提出的基于等效慢度的快速推进法的地震波走时求解系统，能够提高地震波走时的计算精度，且保证地震波走时的计算效率。

图6为本申请实施例提供的终端的结构示意图。该终端可以包括：

存储器501、处理器502及存储在存储器501上并可在处理器502上运行的计算机程序。

处理器502执行程序时实现上述实施例中提供的基于等效慢度的快速推进法的地震波走时求解方法。

进一步地，终端还包括：

通信接口503，用于存储器501和处理器502之间的通信。

存储器501，用于存放可在处理器502上运行的计算机程序。

存储器501可能包含高速RAM存储器，也可能还包括非易失性存储器（non-volatile memory），例如至少一个磁盘存储器。

如果存储器501、处理器502和通信接口503独立实现，则通信接口503、存储器501和处理器502可以通过总线相互连接并完成相互间的通信。总线可以是工业标准体系结构（Industry Standard Architecture，简称为ISA）总线、外部设备互连（PeripheralComponent，简称为PCI）总线或扩展工业标准体系结构（Extended IndustryStandardArchitecture，简称为EIS）总线等。总线可以分为地址总线、数据总线、控制总线等。为便于表示，图6中仅用一条粗线表示，但并不表示仅有一根总线或一种类型的总线。

可选地，在具体实现上，如果存储器501、处理器502及通信接口503，集成在一块芯片上实现，则存储器501、处理器502及通信接口503可以通过内部接口完成相互间的通信。

处理器502可能是一个中央处理器（Central Processing Unit，简称为CPU），或者是特定集成电路（Application Specific Integrated Circuit，简称为ASIC），或者是被配置成实施本申请实施例的一个或多个集成电路。

本实施例还提供一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现如上的基于等效慢度的快速推进法的地震波走时求解方法。

本申请一个实施例提供一种计算机程序产品，包括计算机程序，所述计算机程序被处理器执行时实现如本申请图1所对应的实施例中任意实施例提供的基于等效慢度的快速推进法的地震波走时求解方法。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本申请的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不必须针对的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任一个或N个实施例或示例中以合适的方式结合。此外，在不相互矛盾的情况下，本领域的技术人员可以将本说明书中描述的不同实施例或示例以及不同实施例或示例的特征进行结合和组合。

此外，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括至少一个该特征。在本申请的描述中，“N个”的含义是至少两个，例如两个，三个等，除非另有明确具体的限定。

流程图中或在此以其他方式描述的任何过程或方法描述可以被理解为，表示包括一个或N个用于实现定制逻辑功能或过程的步骤的可执行指令的代码的模块、片段或部分，并且本申请的优选实施方式的范围包括另外的实现，其中可以不按所示出或讨论的顺序，包括根据所涉及的功能按基本同时的方式或按相反的顺序，来执行功能，这应被本申请的实施例所属技术领域的技术人员所理解。

在流程图中表示或在此以其他方式描述的逻辑和/或步骤，例如，可以被认为是用于实现逻辑功能的可执行指令的定序列表，可以具体实现在任何计算机可读存储介质中，以供指令执行系统、装置或设备（如基于计算机的系统、包括处理器的系统或其他可以从指令执行系统、装置或设备取指令并执行指令的系统）使用，或结合这些指令执行系统、装置或设备而使用。就本说明书而言，“计算机可读存储介质”可以是任何可以包含、存储、通信、传播或传输程序以供指令执行系统、装置或设备或结合这些指令执行系统、装置或设备而使用的装置。计算机可读存储介质的更具体的示例（非穷尽性列表）包括以下：具有一个或N个布线的电连接部（电子装置），便携式计算机盘盒（磁装置），随机存取存储器（RAM），只读存储器（ROM），可擦除可编辑只读存储器（EPROM或闪速存储器），光纤装置，以及便携式光盘只读存储器（CDROM）。另外，计算机可读存储介质甚至可以是可在其上打印所述程序的纸或其他合适的介质，因为可以通过对纸或其他介质进行光学扫描，接着进行编辑、解译或必要时以其他合适方式进行处理来以电子方式获得所述程序，然后将其存储在计算机存储器中。

应当理解，本申请的各部分可以用硬件、软件、固件或它们的组合来实现。在上述实施方式中，N个步骤或方法可以用存储在存储器中且由合适的指令执行系统执行的软件或固件来实现。如，如果用硬件来实现和在另一实施方式中一样，可用本领域公知的下列技术中的任一项或他们的组合来实现：具有用于对数据信号实现逻辑功能的逻辑门电路的离散逻辑电路，具有合适的组合逻辑门电路的专用集成电路，可编程门阵列（PGA），现场可编程门阵列（FPGA）等。

本技术领域的普通技术人员可以理解实现上述实施例方法携带的全部或部分步骤是可以通过程序来指令相关的硬件完成，所述的程序可以存储于一种计算机可读存储介质中，该程序在执行时，包括方法实施例的步骤之一或其组合。

此外，在本申请各个实施例中的各功能单元可以集成在一个处理模块中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个模块中。上述集成的模块既可以采用硬件的形式实现，也可以采用软件功能模块的形式实现。所述集成的模块如果以软件功能模块的形式实现并作为独立的产品销售或使用时，也可以存储在一个计算机可读取存储介质中。

上述提到的存储介质可以是只读存储器，磁盘或光盘等。尽管上面已经示出和描述了本申请的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本申请的限制，本领域的普通技术人员在本申请的范围内可以对上述实施例进行变化、修改、替换和变型。

最后应说明的是：以上各实施例仅用以说明本申请的技术方案，而非对其限制；尽管参照前述各实施例对本申请进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本申请各实施例技术方案的范围。

Claims (7)

1.一种基于等效慢度的快速推进法的地震波走时求解方法，其特征在于，所述基于等效慢度的快速推进法的地震波走时求解方法包括：

获取在初始离散网格内的震源点位置，计算所述震源点位置与所述初始离散网格内最近的网格点的偏移向量，根据所述偏移向量对所述初始离散网格进行移动，得到目标网格；其中，所述震源点位置在所述目标网格内的空间点上；

建立在所述目标网格内空间点的走时、等效慢度以及震源点到所述空间点的距离的关系方程，并根据所述关系方程得到基于等效慢度的程函方程；

对所述程函方程进行离散，得到用于求解所述空间点对应的所述等效慢度的目标方程组；

根据所述目标方程组和所述关系方程，计算所述目标网格内多个所述空间点各自对应的地震波走时；

根据多个所述地震波走时各自对应的等效慢度，得到在所述初始离散网格内多个网格点各自对应的当前地震波走时；

所述关系方程中的所述等效慢度为在所述目标网格内所述震源点到所述空间点的走时与所述震源点到所述空间点的距离的商；

所述根据所述关系方程得到基于等效慢度的程函方程，具体包括：

获取用于计算地震波走时的基本方程；

将所述关系方程中的等效慢度和距离替代所述基本方程中的走时，得到基于等效慢度的程函方程；

所述目标方程组包括离散方程和差分算子方程组；

所述对所述程函方程进行离散，得到用于求解所述空间点对应的所述等效慢度的目标方程组，具体包括：

采用迎风格式对所述程函方程进行求解，得到离散方程；

将所述关系方程中的所述等效慢度和所述距离代入所述离散方程中的差分算子，得到差分算子方程组；

所述差分算子方程组包括一阶精度方程和二阶精度方程；

所述将所述关系方程中的所述等效慢度和所述距离代入所述离散方程中的差分算子，得到差分算子方程组，具体包括：

确定所述空间点的走时的精度阶位，其中，所述精度阶位包括一阶精度和二阶精度；

将所述关系方程中的所述等效慢度和所述距离代入所述离散方程中所述一阶精度对应的差分算子，得到一阶精度方程；

将所述关系方程中的所述等效慢度和所述距离代入所述离散方程中所述二阶精度对应的差分算子，得到二阶精度方程。

2.根据权利要求1所述的基于等效慢度的快速推进法的地震波走时求解方法，其特征在于，所述目标网格内多个所述空间点各自对应的走时包括初始走时和中间走时；

所述根据所述目标方程组和所述关系方程，计算所述目标网格内多个所述空间点各自对应的地震波走时，具体包括：

获取所述目标网格中的震源点位置；

根据所述震源点位置、所述目标方程组和所述关系方程，计算所述震源点位置对应的初始走时和所述震源点位置周围的多个所述空间点各自对应的中间走时；其中，所述震源点位置周围的多个所述空间点为筛选空间点，所述筛选空间点中已剔除过作为所述震源点位置的所述空间点；

根据多个所述空间点各自对应的所述中间走时，确定多个所述空间点中的目标空间点，并将所述目标空间点作为新的震源点位置计算新的中间走时，直至所有所述空间点各自对应的走时均已计算。

3.根据权利要求2所述的基于等效慢度的快速推进法的地震波走时求解方法，其特征在于，所述根据所述震源点位置、所述目标方程组和所述关系方程，计算所述震源点位置对应的初始走时和所述震源点位置周围的多个所述空间点各自对应的中间走时，具体包括：

根据所述目标方程组，计算所述震源点位置对应的初始等效慢度，并根据所述初始等效慢度和所述关系方程得到所述初始走时；

根据所述目标方程组，计算所述震源点位置周围的多个所述空间点各自对应的中间等效慢度，并根据各所述中间等效慢度和所述关系方程得到多个所述空间点各自对应的中间走时。

4.根据权利要求3所述的基于等效慢度的快速推进法的地震波走时求解方法，其特征在于，所述根据多个所述空间点各自对应的所述中间走时，确定多个所述空间点中的目标空间点，并将所述目标空间点作为新的震源点位置计算新的中间走时，直至所有所述空间点各自对应的走时均已计算，具体包括：

根据多个所述空间点各自对应的所述中间走时，确定多个所述中间走时中的最小值走时；其中，多个所述中间走时为已计算过的全部走时中剔除了作为过所述震源点位置对应的走时；

将所述最小值走时对应的所述空间点作为目标空间点，并将所述目标空间点作为新的震源点位置；

根据所述新的震源点位置、所述目标方程组和所述关系方程，计算所述新的震源点位置周围的未计算过的所述空间点对应的新的中间走时，直至所有所述空间点均已计算过走时。

5.一种基于等效慢度的快速推进法的地震波走时求解系统，其特征在于，所述基于等效慢度的快速推进法的地震波走时求解系统包括：

网格确定模块，用于获取在初始离散网格内的震源点位置，计算所述震源点位置与所述初始离散网格内最近的网格点的偏移向量，根据所述偏移向量对所述初始离散网格进行移动，得到目标网格；其中，所述震源点位置在所述目标网格内的空间点上；

方程建立模块，用于建立在目标网格内空间点的走时、等效慢度以及震源点到所述空间点的距离的关系方程，并根据所述关系方程得到基于等效慢度的程函方程；

所述关系方程中的所述等效慢度为在所述目标网格内所述震源点到所述空间点的走时与所述震源点到所述空间点的距离的商；

所述方程建立模块，还用于：获取用于计算地震波走时的基本方程；将所述关系方程中的等效慢度和距离替代所述基本方程中的走时，得到基于等效慢度的程函方程；

方程离散模块，用于对所述程函方程进行离散，得到用于求解所述空间点对应的所述等效慢度的目标方程组；

所述目标方程组包括离散方程和差分算子方程组；

所述方程离散模块，还用于：采用迎风格式对所述程函方程进行求解，得到离散方程；将所述关系方程中的所述等效慢度和所述距离代入所述离散方程中的差分算子，得到差分算子方程组；

所述差分算子方程组包括一阶精度方程和二阶精度方程；

所述方程离散模块，还用于：确定所述空间点的走时的精度阶位，其中，所述精度阶位包括一阶精度和二阶精度；将所述关系方程中的所述等效慢度和所述距离代入所述离散方程中所述一阶精度对应的差分算子，得到一阶精度方程；将所述关系方程中的所述等效慢度和所述距离代入所述离散方程中所述二阶精度对应的差分算子，得到二阶精度方程；

走时计算模块，用于根据所述目标方程组和所述关系方程，计算所述目标网格内多个所述空间点各自对应的地震波走时；

当前走时确定模块，用于根据多个所述地震波走时各自对应的等效慢度，得到在所述初始离散网格内多个网格点各自对应的当前地震波走时。

6.一种终端，其特征在于，所述终端包括：存储器、处理器及存储在所述存储器上并可在所述处理器上运行的基于等效慢度的快速推进法的地震波走时求解程序，所述基于等效慢度的快速推进法的地震波走时求解程序被所述处理器执行时实现如权利要求1-4任一项所述的基于等效慢度的快速推进法的地震波走时求解方法的步骤。

7.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有基于等效慢度的快速推进法的地震波走时求解程序，所述基于等效慢度的快速推进法的地震波走时求解程序被处理器执行时实现如权利要求1-4任一项所述的基于等效慢度的快速推进法的地震波走时求解方法的步骤。