CN116774292A - 一种地震波走时确定方法、系统、电子设备及存储介质 - Google Patents

Abstract

本发明公开一种地震波走时确定方法、系统、电子设备及存储介质，涉及地震数据处理领域，该方法包括：首先利用快速行进法求解程函方程获取初始走时值；然后结合初始走时值，利用间断有限元求解程函方程的弱形式，确定地震波走时，提高走时精度。这两个关键步骤组合构成“高精度求解技术方案”的核心。利用快速行进法求解程函方程，虽然与快速扫描法相比有一定的效率差距，但仍然是高效率的方法之一，且实践中精度高于快速扫描方法，因此能够为间断有限元法提供网格单元每个节点正确的方向信息。最后使用二阶走时基函数，实现网格单元内走时计算精度的提高。

Description

一种地震波走时确定方法、系统、电子设备及存储介质

技术领域

本发明涉及地震数据处理领域，特别是涉及一种地震波走时确定方法、系统、电子设备及存储介质。

背景技术

地球物理方法中，地震方法在海洋或陆地上应用广泛，主要通过振动传感器阵列系统观测人工地震的地下地层反射波场信息，然后通过地震成像和反演流程，获取地下介质的构造或者属性信息，帮助地球物理学家和地质学家勘察海洋或者陆地地下矿产资源。地震波传播过程中，走时是非常关键的信息，很多地震成像和反演方法需要进行走时计算，程函方程是描述走时的基本数学模型。因此，如何高效准确地求解程函方程成为了地球物理学家的重要研究工作。

程函方程与波动方程波场描述不同，主要是描述波前运动，其可进一步简化为射线方程组，描述波前面上少数点的行进过程（射线），然后利用射线追踪方法龙格库塔递推求解，但射线方法常常会遇到阴影区问题，可能造成某些区域射线无法到达。有限差分解是计算走时最有效的方法之一(Vidale，1990；Popovici，1991），可称之为Vidale类方法，该方法能够在规则模型网格上实施，可以覆盖整个空间。快速行进法(Sethian，1996，1999，2002)被广泛用于求解程函方程，该方法包括迎风差分、窄带技术和堆排序技术，迎风差分和窄带技术用于地震波传播过程中波前扩展的过程，堆排序运用二叉树堆从全局中选取最小走时值。快速行进方法具有无条件稳定，计算精确度高，算法效率高的特点，已经广泛用于地震波传播的走时计算。快速扫描方法（Zhao，2007)也被广泛应用于走时的计算，它最大的优势是计算效率比快速行进更高，但实践中，快速扫描法比快速行进方法的精度相对比较低。

发明内容

本发明的目的是提供一种地震波走时确定方法、系统、电子设备及存储介质，以提高计算地震波走时的精度。

为实现上述目的，本发明提供了如下方案：

一种地震波走时确定方法，包括：

获取地震波的震源位置点和观测位置点；

将所述震源位置点和所述观测位置点之间的区域离散为（I-1）×（J-1）个网格；I为水平方向网格节点数量；J为垂直方向网格节点数量；

建立地震波的波前在任一所述网格节点的程函方程；所述程函方程是根据地震波的波前在任一所述网格节点的走时和慢度模型确定的；所述慢度模型为根据所述地震波的波前在任一所述网格节点的地震速度模型确定的；

基于所述程函方程，利用逆风差分法和快速行进算法，确定所有网格节点的初始走时；

基于所有网格节点的初始走时以及每个网格节点的程函方程，利用间断伽辽金法，确定地震波走时。

可选地，基于所述程函方程，利用逆风差分法和快速行进算法，确定所有网格节点的初始走时，具体包括：

利用逆风差分法，确定所述程函方程的差分方程；

基于所述差分方程，利用快速行进算法，确定所有网格节点的初始走时。

可选地，基于所有网格节点的初始走时以及每个网格节点的程函方程，利用间断伽辽金法，确定地震波走时，具体包括：

根据当前网格节点的程函方程，利用间断有限元方法，确定每个有限元单元的弱形式方程；

根据所述弱形式方程和加权函数系数，确定代数方程；

根据当前网格节点的初始走时，对所述代数方程进行求解，得到当前有限元单元的中心点平均值、横向梯度和纵向梯度；

根据当前有限元单元的中心点平均值、横向梯度和纵向梯度，确定当前有限元单元的走时；

根据所有有限元单元的走时，确定所述地震波走时。

可选地，根据当前有限元单元的中心点平均值、横向梯度和纵向梯度，确定当前有限元单元的走时，具体包括：

利用公式，确定当前有限元单元的走时；其中，/>为中心点平均值；u_ij为横向梯度；v_ij为纵向梯度；d_x为横坐标方向采样间隔；d_y为纵坐标方向采样间隔；x为当前有限单元中任意一点的横向坐标；y为当前有限单元中任意一点的纵向坐标；x_i为当前有限单元的中线点位置的横向坐标；y_j为当前有限单元的中心点位置的纵向坐标。

一种地震波走时确定系统，包括：

数据获取模块，用于获取地震波的震源位置点和观测位置点；

离散模块，用于将所述震源位置点和所述观测位置点之间的区域离散为（I-1）×（J-1）个网格；I为水平方向网格节点数量；J为垂直方向网格节点数量；

程函方程建立模块，用于建立地震波的波前在任一所述网格节点的程函方程；所述程函方程是根据地震波的波前在任一所述网格节点的走时和慢度模型确定的；所述慢度模型为根据所述地震波的波前在任一所述网格节点的地震速度模型确定的；

初始走时确定模块，用于基于所述程函方程，利用逆风差分法和快速行进算法，确定所有网格节点的初始走时；

地震波走时确定模块，用于基于所有网格节点的初始走时以及每个网格节点的程函方程，利用间断伽辽金法，确定地震波走时。

一种电子设备，包括：存储器及处理器，所述存储器用于存储计算机程序，所述处理器运行所述计算机程序以使所述电子设备执行上述的地震波走时确定方法。

一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现上述的地震波走时确定方法。

根据本发明提供的具体实施例，本发明公开了以下技术效果：

本发明的地震波走时确定方法、系统、电子设备及存储介质，首先利用快速行进法求解程函方程获取初始走时值；然后结合初始走时值，利用间断有限元求解程函方程的弱形式，确定地震波走时，提高走时精度。这两个关键步骤组合构成“高精度求解技术方案”的核心。利用快速行进法求解程函方程，虽然与快速扫描法相比有一定的效率差距，但仍然是高效率的方法之一，且实践中精度高于快速扫描方法，因此能够为间断有限元法提供网格单元每个节点正确的方向信息。步最后使用二阶走时基函数，实现网格单元内走时计算精度的提高。

附图说明

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

图1为本发明提供的地震波走时确定方法流程图；

图2为本发明地震波走时确定方法在实际应用中的流程图；

图3为不同传统方法确定地震波走时的结果对比图；

图4为本发明方法与传统方法确定地震波走时的结果对比图；

图5为本发明方法与传统方法确定地震波走时的结果对比的局部放大图。

具体实施方式

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

本发明的目的是提供一种地震波走时确定方法、系统、电子设备及存储介质，以提高计算地震波走时的精度。

本发明在规则网格上，在程函方程波前快速行进法基础上，利用间断有限元的间断特性和网格单元中二阶走时表示形式，解决网格单元内计算精度不足问题，提出了一种基于程函方程波前行进的地震波走时高精度求解技术，从而强化快速行进法的高精度优势。

为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。

如图1和图2所示，本发明提供的地震波走时确定方法，包括：

步骤101：获取地震波的震源位置点和观测位置点。

在实际应用中，获取地震速度模型，用v（x，y）表示，x、y是地震波水平和垂直坐标。

获取所有震源点位置s_x和s_y，下面以（s_x，s_y）表示某一震源点坐标位置。

获取网格间隔d_x和d_y，那么基于d_x和d_y可得到离散网格（i，j)，i=0,1,...,I-1，j=0,1,...,J-1，即有I×J个网格节点，有(I-1)×(J-1)矩形网格。x=i×d_x，y=j×d_y。

步骤102：将所述震源位置点和所述观测位置点之间的区域离散为（I-1）×（J-1）个网格；I为水平方向网格节点数量；J为垂直方向网格节点数量。

步骤103：建立地震波的波前在任一所述网格节点的程函方程；所述程函方程是根据地震波的波前在任一所述网格节点的走时和慢度模型确定的；所述慢度模型为根据所述地震波的波前在任一所述网格节点的地震速度模型确定的。

在实际应用中，根据输入地震速度模型v（x，y），得到慢度模型s（x,y）=1.0/v（x,y）。速度模型中，地震波行进的波前在（x，y）位置的走时T和（x，y）位置的慢度模型s的关系表示为，该方程为程函方程，如公式（1），实质表示地震波波前面上走时梯度的模，表示如下：

（1）。

步骤104：基于所述程函方程，利用逆风差分法和快速行进算法，确定所有网格节点的初始走时。

作为一种可选地实施方式，步骤104，具体包括：

S1：利用逆风差分法，确定所述程函方程的差分方程。

在离散网格中，快速行进方法通过逆风差分的方法求解程函方程，微分形式用下列差分替代，公式（1）的差分方程表示如下：

（2）。

其中，为在节点（i，j）沿x方向的n阶前向差分，/>为在节点（i，j）沿x方向的n阶前后差分，/>为在节点（i，j）沿y方向的n阶前向差分，/>分别是在节点（i，j）沿y方向的n阶后向差分，n是差分阶数，这里取二阶差分，边界处降为一阶差分。一阶差分和二阶差分表示形式如下：

一阶差分表示为：

（3）。

T_i,j为网格节点（i,j)处的旅行时；T_i-1,j为网格节点（i-1,j)处的旅行时；T_i+1,j为网格节点（i+1,j)处的旅行时；T_i，j-1为网格节点（i,j-1)处的旅行时；T_i，j+1为网格节点（i,j+1)处的旅行时。

二阶差分表示为：

（4）。

公式（2）表示声波的波场是从最小值传播到大值。其中，T_i-2,j为为网格节点（i-2,j)处的旅行时；T_i+2,j为网格节点（i+2,j)处的旅行时；T_i，j-2为网格节点（i,j-2)处的旅行时；T_i，j+2为网格节点（i,j+2)处的旅行时。

S2：基于所述差分方程，利用快速行进算法，确定所有网格节点的初始走时。

在实际应用中，快速行进方法计算的走时就是由最小值（震源点处）向外扩展构建方程式|▽T|_i，j的解。求解步骤如下：

（1）初始化所有网格节点的走时为T_i,j=+∞，都作为未处理节点。

（2）根据震源点位置s_x和s_y和网格间隔d_x和d_y，确定炮点在X坐标方向的坐标所在网格点索引isx=INT（s_x/d_x）和炮点在Y坐标方向的坐标所在网格点索引jsy=INT（s_y/d_y），其中，INT表示四舍五入取整。然后，根据isx和isy和震源点速递v(s_x,s_y)，求取临近节点的走时：

（5）。

（6）。

（7）。

（8）。

其中，T_isx,jsy为网格节点（isx,jsy)处的旅行时；T_isx,jsy+1为网格节点（isx,jsy+1)处的旅行时；T_isx+1,jsy+1为网格节点（isx+1,jsy+1)处的旅行时；T_isx+1,jsy为网格节点（isx+1,jsy)处的旅行时；震源点临近节点划为已处理节点区。

（3）计算处理节点区波前的临近未处理节点（作为待处理节点）的走时，根据上述公式（2）形成一元二次方程组，得到临近节点的走时，划入已处理节点区。

（4）重复（3）步骤直到完成所有网格节点的走时计算，最终得到网格节点的初始走时，i∈[0,I-1],j∈[0,J-1]。/>就是差分方程的|▽T|_i，j。

在已知震源位置和地震速度模型的基础上，该步骤能够快速得到所有离散网格点初始走时，在每个网格单元上为步骤105提供每个网格节点计算的初始走时信息，主要计算每个节点波前的方向信息。

步骤105：基于所有网格节点的初始走时以及每个网格节点的程函方程，利用间断伽辽金法，确定地震波走时。

作为一种可选地实施方式，步骤105，具体包括：

根据当前网格节点的程函方程，利用间断有限元方法，确定每个有限元单元的弱形式方程。

根据所述弱形式方程和加权函数系数，确定代数方程。

根据当前网格节点的初始走时，对所述代数方程进行求解，得到当前有限元单元的中心点平均值、横向梯度和纵向梯度。

根据当前有限元单元的中心点平均值、横向梯度和纵向梯度，确定当前有限元单元的走时。具体包括：

利用公式，确定当前有限元单元的走时；其中，为中心点平均值；u_ij为横向梯度；v_ij为纵向梯度；d_x为横坐标方向采样间隔；d_y为纵坐标方向采样间隔；x为当前有限单元中任意一点的横向坐标；y为当前有限单元中任意一点的纵向坐标；x_i为当前有限单元的中线点位置的横向坐标；y_j为当前有限单元的中心点位置的纵向坐标。

根据所有有限元单元的走时，确定所述地震波走时。

在实际应用中，方程的间断有限元弱形式表示如下：

（9）。

其中，为有限元单元左边界位置；/>为有限元单元右边界位置；/>为有限元单元上边界位置；/>为有限元单元下边界位置，/>为有限元单元左边界位置的加权函数系数；/>有限元单元右边界位置的加权函数系数。左边的后四项主要是边界的通量信息求取。

在x方向上，。

在y方向上，。

其中，为有限元单元右边界位置的旅行时差分（类似间断点沿y方向导数）；/>为有限元单元右边界位置的右边旅行时；/>为有限元单元右边界位置的左边旅行时；/>为有限元单元下边界位置的旅行时差分（类似间断点沿x方向导数）；/>为有限元单元下边界位置的下边的旅行时；/>为有限元单元下边界位置的上边旅行时。

（10）。

（11）。

其中，为右边相邻有限单元左边界的旅行时；/>为当前有限单元右边界的旅行时；/>为下边相邻有限单元上边界的旅行时；为当前有限单元下边界的旅行时。

这是间断有限元推导的基本步骤，利用加权余量方法形成弱形式方程，表示矩形单元（element）里的积分。这两需要进一步设定的参数包括：

（1）设定网格单元上波前走时近似解：

（12）。

公式（12）是二阶基函数。，u_ij，v_ij是有限元单元（element）上的中心点平均值、横向梯度及纵向梯度，是待求变量。

（2）设定加权函数系数w。间断有限元中通常使用的权函数与基函数(形函数或者插值函数）一致，所以该系数是向量系数中的一个分量。

（3）通量系数α_l,ij,α_r,ij,α_t,ij,α_d,ij，是依赖于本单元（element）及其邻近单元的数值解的局部常数。根据地震波走时的传播特点，单元的求解需要利用相邻单元的走时信息，主要是依赖于单元之间信息关系；l，r，t和b分别表示单元的左（left）、右（right）、上（top）和底（bottom）。在弱形式中，通过α_l,ij,α_r,ij,α_t,ij,α_d,ij来确定左右上下四个方向上的相邻单元走时是否参与计算，从而保证数值解符合初至走时沿最短路径扩展等时波前面的因果性。

弱形式中的通量系数α_l,ij,α_r,ij,α_t,ij,α_d,ij，实质上相当于单元节点上的波前走时的方向，可以表示为：

（13）。

其中，这里的，u_ij，v_ij将由/>计算得到。

设定上述所需要的参数后。将w选择，其中分量代入间断有限元弱形式的公式（9）后，可以得到三个代数方程。代数方程通过进一步的简化可以得到便于计算的形式：

（14）。

其中，为横向通量常系数均值；/>为纵向通量常系数均值；/>为单元通量常系数差；/>为单元通量常系数差值及横向通量常系数差值之和；/>为单元通量常系数差值及纵向通量常系数差值之和；R_1,ij为权值1弱形式积分已知值；R_2,ij为权值/>弱形式积分已知值；R_3,ij为权值/>弱形式积分已知值。

（15）。

（16）。

方程（14）就是求解单元走时未知数，u_ij，v_ij的方程组，然后通过公式（12）得到所有网格点的间断有限元解。

基于快速扫描法的间断有限元算法具体步骤为：

（1）在步骤104基础上得到，i∈[0,I-1],j∈[0,J-1]后，在四个方向上进行扫描，求解各个有限元单元的走时。分别按照四个方向扫描一遍：

。

求解每一个有限元单元的走时解，有限元单元求解如下步骤（2）。

（2）由于方程（4）包含未知量u_ij，v_ij。通过走时得到：

（17）。

作为/>，u_ij，v_ij，代入公式（13），进一步求解方程（14）的系数，最后求解方程组得到/>，u_ij，v_ij的解。然后，基于公式（12）得到网格点上的解。

该部分间断有限元主要利用间断有限元的间断特性（每个网格单元有自己独立的基函数）和二阶走时表示能力，实现网格单元内走时计算精度的提高。

为验证本发明所提的程函方程快速行进间断有限元求解方法的有效性，设计一个均匀介质模型，速度为1000m/s，水平和垂直采样间隔都20m，震源置于（1000m，1000m）位置。用传统的快速扫描法（即图3中点线）和快速行进法（图3中虚线）与解析解（图3中实线）比较，快速行进法明显由于快速扫描法，这里表示步骤104能够提高较高精度的初始波前走时，能够为本发明高精度求解技术方案提供基础。

在上述模型来测试本发明（快速行进法+间断有限元），结果如图4中点线，与传统快速行进法（图4中虚线）和解析解（图4中实线）比较，它们都接近真实解。放大后，如图5所示，本发明的结果（点线）明显好于传统的快速行进法（虚线）。该例子的发明效果非常明显，通过与解析解（正确解）比较，本发明地震波场波前走时求解明显好于传统方法。

实施例二

为了执行上述实施例一对应的方法，以实现相应的功能和技术效果，下面提供一种地震波走时确定系统，包括：

数据获取模块，用于获取地震波的震源位置点和观测位置点。

离散模块，用于将所述震源位置点和所述观测位置点之间的区域离散为（I-1）×（J-1）个网格；I为水平方向网格节点数量；J为垂直方向网格节点数量。

程函方程建立模块，用于建立地震波的波前在任一所述网格节点的程函方程；所述程函方程是根据地震波的波前在任一所述网格节点的走时和慢度模型确定的；所述慢度模型为根据所述地震波的波前在任一所述网格节点的地震速度模型确定的。

初始走时确定模块，用于基于所述程函方程，利用逆风差分法和快速行进算法，确定所有网格节点的初始走时。

地震波走时确定模块，用于基于所有网格节点的初始走时以及每个网格节点的程函方程，利用间断伽辽金法，确定地震波走时。

实施例三

本发明提供了一种电子设备，包括：存储器及处理器，所述存储器用于存储计算机程序，所述处理器运行所述计算机程序以使所述电子设备执行实施例一的地震波走时确定方法。

实施例四

本发明提供了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现实施例一的地震波走时确定方法。

本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。

本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (7)

1.一种地震波走时确定方法，其特征在于，包括：

获取地震波的震源位置点和观测位置点；

将所述震源位置点和所述观测位置点之间的区域离散为（I-1）×（J-1）个网格；I为水平方向网格节点数量；J为垂直方向网格节点数量；

建立地震波的波前在任一所述网格节点的程函方程；所述程函方程是根据地震波的波前在任一所述网格节点的走时和慢度模型确定的；所述慢度模型为根据所述地震波的波前在任一所述网格节点的地震速度模型确定的；

基于所述程函方程，利用逆风差分法和快速行进算法，确定所有网格节点的初始走时；

基于所有网格节点的初始走时以及每个网格节点的程函方程，利用间断伽辽金法，确定地震波走时。

2.根据权利要求1所述的地震波走时确定方法，其特征在于，基于所述程函方程，利用逆风差分法和快速行进算法，确定所有网格节点的初始走时，具体包括：

利用逆风差分法，确定所述程函方程的差分方程；

基于所述差分方程，利用快速行进算法，确定所有网格节点的初始走时。

3.根据权利要求1所述的地震波走时确定方法，其特征在于，基于所有网格节点的初始走时以及每个网格节点的程函方程，利用间断伽辽金法，确定地震波走时，具体包括：

根据当前网格节点的程函方程，利用间断有限元方法，确定每个有限元单元的弱形式方程；

根据所述弱形式方程和加权函数系数，确定代数方程；

根据当前网格节点的初始走时，对所述代数方程进行求解，得到当前有限元单元的中心点平均值、横向梯度和纵向梯度；

根据当前有限元单元的中心点平均值、横向梯度和纵向梯度，确定当前有限元单元的走时；

根据所有有限元单元的走时，确定所述地震波走时。

4.根据权利要求3所述的地震波走时确定方法，其特征在于，根据当前有限元单元的中心点平均值、横向梯度和纵向梯度，确定当前有限元单元的走时，具体包括：

利用公式，确定当前有限元单元的走时；其中，/>为中心点平均值；u_ij为横向梯度；v_ij为纵向梯度；d_x为横坐标方向采样间隔；d_y为纵坐标方向采样间隔；x为当前有限单元中任意一点的横向坐标；y为当前有限单元中任意一点的纵向坐标；x_i为当前有限单元的中线点位置的横向坐标；y_j为当前有限单元的中心点位置的纵向坐标。

5.一种地震波走时确定系统，其特征在于，包括：

数据获取模块，用于获取地震波的震源位置点和观测位置点；

离散模块，用于将所述震源位置点和所述观测位置点之间的区域离散为（I-1）×（J-1）个网格；I为水平方向网格节点数量；J为垂直方向网格节点数量；

程函方程建立模块，用于建立地震波的波前在任一所述网格节点的程函方程；所述程函方程是根据地震波的波前在任一所述网格节点的走时和慢度模型确定的；所述慢度模型为根据所述地震波的波前在任一所述网格节点的地震速度模型确定的；

初始走时确定模块，用于基于所述程函方程，利用逆风差分法和快速行进算法，确定所有网格节点的初始走时；

地震波走时确定模块，用于基于所有网格节点的初始走时以及每个网格节点的程函方程，利用间断伽辽金法，确定地震波走时。

6.一种电子设备，其特征在于，包括：存储器及处理器，所述存储器用于存储计算机程序，所述处理器运行所述计算机程序以使所述电子设备执行权利要求1-4任一项所述的地震波走时确定方法。

7.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现权利要求1-4任一项所述的地震波走时确定方法。