CN110824558B - 一种地震波数值模拟方法 - Google Patents

Abstract

本发明公开一种地震波数值模拟方法。该方法包括：获取地震波速度模型；确定地震波速度模型适应的剖分节点径向半径范围；确定地震波速度模型中的地震波速度对应的剖分节点径向半径；确定地震波速度模型的多个剖分节点、边界内点、权重因子和地震波场传播角；确定各剖分节点对应的最邻近的n个剖分节点；确定最优化形变参数；确定所有剖分节点对应的差分加权系数；根据各剖分节点、边界内点、权重因子、地震波场传播角和所有剖分节点对应的差分加权系数采用地震波场计算公式进行地震波场延拓，确定数值模拟的地震波场。本发明可以获得高精度的地震波数值模拟结果，有效提升地震数据处理的精度及稳定性。

Description

一种地震波数值模拟方法

技术领域

本发明涉及地震学领域，特别是涉及一种地震波数值模拟方法。

背景技术

地震波数值模拟是通过数学手段模拟地下介质对地震波响应规律的一种技术。利用地震波数值模拟可以使地球物理学家更加直观地理解地震波在地下介质中的传播规律和特征，对地震的预防和研究提供指导，也可以为观测地震数据的处理与解释提供理论依据。它不仅是地震学的重要理论基础，也是整个地震数据处理流程的技术核心，尤其在偏移成像及速度建模的方法流程中占据至关重要的位置。因此，高精度的地震波数值模拟方法对地震数据处理技术在实际数据中的广泛应用至关重要。

随着计算机技术的发展，地震波数值模拟方法得到了广泛的发展。基于对波动方程的求解方式不同，地震波数值模拟方法主要包括谱形式、弱形式和强形式。与弱形式中的有限元、边界元和谱元法以及强形式中的伪谱法和有限体积法相比，有限差分法通过使用一组系数对当前点附近的函数值进行加权来近似偏导数。由于其内存需求小、精度高、易于实现等明显优势，在地震波数值模拟中仍被广泛应用。然而，标准的有限差分法到目前为止在空间和时间采样间隔上达到了极限。欠采样会降低精度，造成数值频散，而过采样会导致存储量和计算量增加。在这种情况下，对任意形状异常或速度变化较大的低速异常进行非常精细的离散化，虽然满足稳定条件，但对于速度较高的深部区域，将导致结果不准确或强烈波动。此外，在笛卡尔网格中精确离散复杂的几何结构(包括崎岖的地形和复杂的地下断层结构)是难以实现的。近年来，变网格有限差分法有了较大的发展。变网格有限差分方法根据精细网格和粗网格之间过渡的处理方式，大致可分为基于插值的、基于坐标变换的和基于差分近似的变网格有限差分方法。然而，基于插值的变网格有限差分法容易引起虚假反射，同时存在不稳定问题；基于坐标变换的变网格有限差分法会遇到曲面边界条件的稳定性问题；基于差分近似的变网格有限差分方法存在网格剖分不灵活、虚假散射强、稳定性差等问题。

发明内容

针对现有技术存在的上述缺陷，提供一种简单灵活、稳定性强、精度高的地震波数值模拟方法。

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

一种地震波数值模拟方法，包括：

获取地震波速度模型；

根据所述地震波速度模型，确定地震波速度模型适应的剖分节点径向半径范围；

根据所述剖分节点径向半径范围，确定所述地震波速度模型中的地震波速度对应的剖分节点径向半径；

根据各所述地震波速度对应的剖分节点径向半径和所述地震波速度模型，采用快速节点剖分法，得到所述地震波速度模型的多个剖分节点、边界内点、权重因子和地震波场传播角；

根据各所述剖分节点采用快速邻域搜索算法，得到各所述剖分节点对应的最邻近的n个剖分节点；

确定最优化形变参数；

根据所述最优化形变参数、各所述地震波速度对应的剖分节点径向半径和各所述剖分节点对应的最邻近的n个剖分节点，确定所有剖分节点对应的差分加权系数；

根据各所述剖分节点、所述边界内点、所述权重因子、所述地震波场传播角和所述所有剖分节点对应的差分加权系数采用地震波场计算公式进行地震波场延拓，确定数值模拟的地震波场。

可选的，所述根据所述地震波速度模型，确定地震波速度模型适应的剖分节点径向半径范围，具体包括：

根据所述地震波速度模型利用稳定性条件和频散关系，确定地震波速度模型适应的剖分节点径向半径范围；

所述稳定条件为：

其中，r为剖分节点径向半径，v为地震波速度，τ为时间步长，a_j为差分系数，k为差分阶数，j表示距离差分中心点的网格数；

所述的频散关系具体为：

其中，v_min为地震波速度模型中的最小速度，f_max为地震波震源频带中的最大频率。

可选的，所述根据所述剖分节点径向半径范围，确定所述地震波速度模型中的地震波速度对应的剖分节点径向半径，具体包括：

根据所述剖分节点径向半径范围，确定每个速度值适应的剖分节点径向半径范围的中值；

根据所述速度值适应的剖分节点径向半径范围的中值，确定所述速度值对应的剖分节点径向半径。

可选的，所述根据各所述地震波速度对应的剖分节点径向半径和所述地震波速度模型，采用快速节点剖分法，得到所述地震波速度模型的多个剖分节点、边界内点、权重因子和地震波场传播角，具体包括：

根据各所述地震波速度对应的剖分节点径向半径和所述地震波速度模型，采用快速节点剖分法进行节点剖分，得到计算域和边界域的剖分节点在已知地震波速度模型中对应的真实空间位置(x,z)；

根据所述真实空间位置(x,z)，从所述地震波速度模型提取相同空间位置处的地震波速度v；

根据所述地震波速度对应的剖分节点径向半径，确定所述地震波速度v对应的剖分节点径向半径r；所述真实空间位置(x,z)、所述地震波速度v和所述剖分节点径向半径r构成所述地震波速度模型的剖分节点(x,z,v,r)；

根据所述剖分节点(x,z,v,r)，确定边界内点(x_inner,z_inner,v_inner,r_inner)；

获取边界域的厚度H_b和边界域内各个剖分节点到计算域最外层剖分节点的垂直距离D_b；

根据所述边界域的厚度H_b和所述边界域内各个剖分节点到计算域最外层剖分节点的垂直距离D_b，确定边界域内各剖分节点的权重因子w；

根据所述边界域内各个剖分节点相对于计算域的空间位置，确定所述边界域内各个剖分节点的地震波场传播角θ。

可选的，所述根据所述边界域的厚度H_b和所述边界域内各个剖分节点到计算域最外层剖分节点的垂直距离D_b，确定边界域内各剖分节点的权重因子w，具体包括：

根据所述边界域的厚度H_b和所述边界域内各个剖分节点到计算域最外层剖分节点的垂直距离D_b采用公式w＝D_b/H_b，计算边界域内各剖分节点的权重因子w。

可选的，所述根据所述边界域内各个剖分节点相对于计算域的空间位置，确定所述边界域内各个剖分节点的地震波场传播角θ，具体包括：

根据边界域内各个剖分节点相对于计算域的空间位置，将所述边界域的各个剖分节点的空间位置进行分类，得到八类位置，所述八类位置分别是上边、左上边、左变、左下变、下边、右下边、右边和右上边；

对所述八类位置分别给于0°、45°、90°、135°、180°、225°、270°和315°的地震波场传播角度，所有的地震波场传播角度构成了边界域内各个剖分节点的地震波场传播角θ。

可选的，所述根据所述最优化形变参数、各所述地震波速度对应的剖分节点径向半径和各所述剖分节点对应的最邻近的n个剖分节点，确定所有剖分节点对应的差分加权系数，具体包括：

根据所述最优化形变参数和所述径向基函数，得到优化的径向基函数；

根据所述优化的径向基函数、所述地震波速度对应的剖分节点径向半径和所述剖分节点对应的最邻近的n个剖分节点，确定各所述剖分节点对应的最邻近的n个剖分节点的差分加权系数；

根据各所述剖分节点对应的最邻近的n个剖分节点的差分加权系数，得到所有剖分节点对应的差分加权系数；

所述所有剖分节点对应的差分加权系数具体为：

其中，Φ_opt为优化的径向基函数，x(x,z)表示地震波速度模型剖分节点的空间坐标矢量，x表示空间坐标矢量的水平分量，z表示空间坐标矢量的垂直分量，下标“0”表示参与当前地震波速度模型剖分节点差分运算的周围第0个最邻近剖分节点，下标“n”表示参与当前地震波速度模型剖分节点差分运算的周围第n个最邻近剖分节点；b_i为第i个剖分节点周围n个最邻近剖分节点对应的差分加权系数，i＝0,1...,n,....n+5，L[]表示偏导数算子。

可选的，所述根据各所述剖分节点、所述边界内点、所述权重因子、所述地震波场传播角和所述所有剖分节点对应的差分加权系数采用地震波场计算公式进行地震波场延拓，确定数值模拟的地震波场，具体包括：

采用公式

确定计算域和边界域所有剖分节点对应的第一类地震波场；

其中，U为每个剖分节点对应的第一类地震波场，t为地震波场传播时间，v为地震波速度，τ为时间步长，b_i为第i个剖分节点周围n个最邻近剖分节点对应的差分加权系数；

确定计算域和边界域所有剖分节点对应的第二类地震波场；

根据所述权重因子对所有剖分节点的所述第一类地震波场和所述第二类地震波场进行加权运算，得到数值模拟的地震波场W；

所述数值模拟的地震波场W具体为：

其中，U为每个剖分节点对应的第一类地震波场，V为每个剖分节点对应的第二类地震波场。

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

1)本发明方法为一种地震波数值模拟方法，与常规地震波数值模拟方法相比，本发明方法简单灵活、稳定性强、精度高。

2)本发明基于速度模型进行灵活网格剖分，自主选取各个网格点位置的差分步长及差分系数，无需人为干扰，可以获得高精度的地震波数值模拟结果，利用该方法可以有效提升地震数据处理的精度及稳定性，为地球浅部油气资源勘探和深部结构探测提供重要的理论和技术支撑。

3)本发明可以广泛用于地震学及勘探地震学领域中，对涉及地震波数值模拟的地球物理研究具有重要理论和应用价值。

附图说明

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

图1为本发明地震波数值模拟方法流程图；

图2是双层介质模型；

图3是图2所示双层介质模型的局部网格剖分图；

图4是图2所示双层介质模型的地震波场快照剖面；

图5是图2所示双层介质模型的地震记录剖面；

图6是本发明提供的Marmousi-2模型；

图7是利用本发明方法所得的图6所示Marmousi-2模型的局部网格剖分图；

图8是图6所示Marmousi-2模型的地震波场快照剖面。

具体实施方式

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

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

图1为本发明地震波数值模拟方法流程图。一种地震波数值模拟方法，包括：

步骤101：获取地震波速度模型；

步骤102：根据所述地震波速度模型，确定地震波速度模型适应的剖分节点径向半径范围，具体包括：

根据所述地震波速度模型利用稳定性条件和频散关系，确定地震波速度模型适应的剖分节点径向半径范围；

所述稳定条件为：

其中，r为剖分节点径向半径，v为地震波速度，τ为时间步长，a_j为差分系数，k为差分阶数，j表示距离差分中心点的网格数；

所述的频散关系具体为：

其中，v_min为地震波速度模型中的最小速度，f_max为地震波震源频带中的最大频率。

步骤103：根据所述剖分节点径向半径范围，确定所述地震波速度模型中的地震波速度对应的剖分节点径向半径，具体包括：

根据所述剖分节点径向半径范围，确定每个速度值适应的剖分节点径向半径范围的中值；

根据所述速度值适应的剖分节点径向半径范围的中值，确定所述速度值对应的剖分节点径向半径。

基于步骤102所确定的地震波速度模型适应的剖分节点径向半径范围，对已知的地震波速度模型中的每一个不同的速度值，选取与该速度适应的剖分节点径向半径范围的中值，作为该速度对应的剖分节点径向半径范围，获得每一个不同的地震波速度对应的剖分节点径向半径r。

步骤104：根据各所述地震波速度对应的剖分节点径向半径和所述地震波速度模型，采用快速节点剖分法，得到所述地震波速度模型的多个剖分节点、边界内点、权重因子和地震波场传播角，具体包括：

根据各所述地震波速度对应的剖分节点径向半径和所述地震波速度模型，采用快速节点剖分法进行节点剖分，得到计算域和边界域的剖分节点在已知地震波速度模型中对应的真实空间位置(x,z)；

根据所述真实空间位置(x,z)，从所述地震波速度模型提取相同空间位置处的地震波速度v；

根据所述地震波速度对应的剖分节点径向半径，确定所述地震波速度v对应的剖分节点径向半径r；所述真实空间位置(x,z)、所述地震波速度v和所述剖分节点径向半径r构成所述地震波速度模型的剖分节点(x,z,v,r)；

根据所述剖分节点(x,z,v,r)，确定边界内点(x_inner,z_inner,v_inner,r_inner)；

获取边界域的厚度H_b和边界域内各个剖分节点到计算域最外层剖分节点的垂直距离D_b；

根据所述边界域的厚度H_b和所述边界域内各个剖分节点到计算域最外层剖分节点的垂直距离D_b，确定边界域内各剖分节点的权重因子w；

根据所述边界域内各个剖分节点相对于计算域的空间位置，确定所述边界域内各个剖分节点的地震波场传播角θ。

所述根据所述边界域的厚度H_b和所述边界域内各个剖分节点到计算域最外层剖分节点的垂直距离D_b，确定边界域内各剖分节点的权重因子w，具体包括：

根据所述边界域的厚度H_b和所述边界域内各个剖分节点到计算域最外层剖分节点的垂直距离D_b采用公式w＝D_b/H_b，计算边界域内各剖分节点的权重因子_w。

所述根据所述边界域内各个剖分节点相对于计算域的空间位置，确定所述边界域内各个剖分节点的地震波场传播角θ，具体包括：

根据边界域内各个剖分节点相对于计算域的空间位置，将所述边界域的各个剖分节点的空间位置进行分类，得到八类位置，所述八类位置分别是上边、左上边、左变、左下变、下边、右下边、右边和右上边；

对所述八类位置分别给于0°、45°、90°、135°、180°、225°、270°和315°的地震波场传播角度，所有的地震波场传播角度构成了边界域内各个剖分节点的地震波场传播角θ。

基于已知的地震波速度模型和确定的地震波速度对应的剖分节点径向半径，利用快速节点剖分法逐次对计算域最外层的地震波速度、计算域内的地震波速度和边界域的地震波速度进行节点剖分，获得每一个计算域和边界域的剖分节点在已知地震波速度模型中对应的真实空间位置(x,z)；依据此空间位置，从已知的地震波速度模型中提取相同空间位置处的地震波速度v，基于步骤103确定的剖分节点径向半径，获得此地震波速度v对应的剖分节点径向半径r；所述真实空间位置(x,z)、地震波速度v和剖分节点径向半径r，构成了已知地震波速度模型的剖分节点(x,z,v,r)，该剖分节点包含两部分：计算域的剖分节点(x_c,z_c,v_c,r_c)和边界域的剖分节点(x_b,z_b,v_b,r_b)；然后，由所述的剖分节点(x,z,v,r)中提取边界域的各个剖分节点(x_b,z_b,v_b,r_b)对应的向内剖分节点(x_inner,z_inner,v_inner,r_inner)，此剖分节点称为边界内点；之后，依据边界域的厚度H_b和边界域内各个剖分节点到计算域最外层剖分节点的垂直距离D_b计算边界域内各剖分节点的权重因子w；根据边界域内各个剖分节点相对于计算域的空间位置，计算边界域内各个剖分节点的地震波场传播角θ；最后，存储上述参数，具体包括已知地震波速度模型的剖分节点(x,z,v,r)，边界内点(x_inner,z_inner,v_inner,r_inner)，权重因子w，地震波场传播角θ；

其中，所述对边界域的地震波速度进行节点剖分具体为：基于边界域的已知的地震波速度模型和步骤103确定的地震波速度对应的剖分节点径向半径，利用快速节点剖分法，按照由内向外逐层进行的顺序，由计算域最外层节点为起始点，获得每一层边界所对应的剖分节点空间位置，地震波速度以及剖分节点径向半径；所有边界域的地震波速度模型均进行上述操作，获得边界域的剖分节点空间位置(x_b,z_b)，地震波速度v_b以及剖分节点径向半径r_b；

其中，所述的由所述的剖分节点(x,z,v,r)中提取边界域的各个剖分节点(x_b,z_b,v_b,r_b)对应的向内剖分节点(x_inner,z_inner,v_inner,r_inner)具体为：针对边界域的每一个剖分节点(x_b,z_b,v_b,r_b)，以该剖分节点空间位置为基点，沿着垂直于边界向计算域的方向，寻找距离此基点最近的一个剖分节点，所找到的距离此基点最近的一个剖分节点即为该基点对应剖分节点的向内剖分节点；对于每一个边界域的各个剖分节点(x_b,z_b,v_b,r_b)，均能找一个与之对应的向内剖分节点，所有的向内剖分节点构成了边界域的各个剖分节点(x_b,z_b,v_b,r_b)对应的向内剖分节点(x_inner,z_inner,v_inner,r_inner)。

步骤105：根据各所述剖分节点采用快速邻域搜索算法，得到各所述剖分节点对应的最邻近的n个剖分节点；

优选的，利用快速邻域搜素算法，针对每一个剖分节点，寻找其周围最近的n个剖分节点；存储每一个剖分节点对应的最邻近的n个剖分节点。

步骤106：确定最优化形变参数；

采用均匀地震波速度模型，基于上述节点剖分方法获得均匀地震波速度模型对应的剖分节点及参数；对形变参数ε进行数值采样，利用径向基函数有限差分地震波数值模拟方法进行不同形变参数下不同差分节点数的地震波数值模拟；构建与形变参数相关的目标泛函，利用最优化算法对此目标泛函进行最优化反演，使其取得极小值，此时目标泛函对应的形变参数，即为最优化形变参数ε_opt；

其中，所述的对形变参数ε进行数值采样，利用径向基函数有限差分地震波数值模拟方法进行不同形变参数下不同差分节点数的地震波数值模拟具体为：对形变参数进行等间隔采样，其中最小形变参数为0.001，最大形变参数为1.0，采样间隔为0.0001；对于每一个采样得到的形变参数，分别给定不同的差分节点数进行径向基函数有限差分地震波数值模拟方法，获得与之对应的模拟地震波场

其中，所述的目标泛函具体为：

其中n＝1,2,...,M

其中，

为模拟的均匀介质地震波场，

为解析的均匀介质地震波场，M为试验采用的最多的最邻近剖分节点个数，obj为目标泛函，

为二范数运算。

步骤107：根据所述最优化形变参数、各所述地震波速度对应的剖分节点径向半径和各所述剖分节点对应的最邻近的n个剖分节点，确定所有剖分节点对应的差分加权系数，具体包括：

根据所述最优化形变参数和所述径向基函数，得到优化的径向基函数；

根据所述优化的径向基函数、所述地震波速度对应的剖分节点径向半径和所述剖分节点对应的最邻近的n个剖分节点，确定各所述剖分节点对应的最邻近的n个剖分节点的差分加权系数；

根据各所述剖分节点对应的最邻近的n个剖分节点的差分加权系数，得到所有剖分节点对应的差分加权系数；

所述所有剖分节点对应的差分加权系数具体为：

其中，Φ_opt为优化的径向基函数，x(x,z)表示地震波速度模型剖分节点的空间坐标矢量，x表示空间坐标矢量的水平分量，z表示空间坐标矢量的垂直分量，下标“0”表示参与当前地震波速度模型剖分节点差分运算的周围第0个最邻近剖分节点，下标“n”表示参与当前地震波速度模型剖分节点差分运算的周围第n个最邻近剖分节点；b_i为第i个剖分节点周围n个最邻近剖分节点对应的差分加权系数，i＝0,1...,n,....n+5，L[]表示偏导数算子。

根据精度需求，选取适当的差分节点数，基于步骤107，将当前剖分节点数对应的最优化形变参数ε_opt带入径向基函数Φ表达式，获得优化的径向基函数Φ_opt，利用此优化的径向基函数Φ_opt、步骤103所得的剖分节点半径及步骤105所得的每一个剖分节点对应的最邻近的n个剖分节点，计算当前节点对应的准确的差分加权系数，进而获得所有剖分节点对应的差分加权系数；

步骤108：根据各所述剖分节点、所述边界内点、所述权重因子、所述地震波场传播角和所述所有剖分节点对应的差分加权系数采用地震波场计算公式进行地震波场延拓，确定数值模拟的地震波场，具体包括：

采用公式

确定计算域和边界域所有剖分节点对应的第一类地震波场；

其中，U为每个剖分节点对应的第一类地震波场，t为地震波场传播时间，v为地震波速度，τ为时间步长，b_i为第i个剖分节点周围n个最邻近剖分节点对应的差分加权系数；

确定计算域和边界域所有剖分节点对应的第二类地震波场，计算域内所有剖分节点对应的第二类地震波场为零；具体的确定方法可参照《An absorbing boundarycondition for acoustic-wave propagation using a mesh-free method》。

根据所述权重因子对所有剖分节点的所述第一类地震波场和所述第二类地震波场进行加权运算，得到数值模拟的地震波场W；

所述数值模拟的地震波场W具体为：

其中，U为每个剖分节点对应的第一类地震波场，V为每个剖分节点对应的第二类地震波场。

本发明与现有技术相比，具有下列优点：

1)本发明与常规地震波数值模拟方法相比，方法简单灵活、稳定性强、精度高。

2)本发明基于速度模型进行灵活网格剖分，自主选取各个网格点位置的差分步长及差分系数，无需人为干扰，可以获得高精度的地震波数值模拟结果，利用该方法可以有效提升地震数据处理的精度及稳定性，为地球浅部油气资源勘探和深部结构探测提供重要的理论和技术支撑。

3)本发明可以广泛用于地震学及勘探地震学领域中，对涉及地震波数值模拟的地球物理研究具有重要理论和应用价值。

为进一步说明本发明的可行性和有效性，下面举两个实例：

实例1：

图2是双层介质速度模型，图3为双层模型局部的节点分布：其中，图3(a)是由快速节点生成方法生成，共计25462个独立节点，图3(b)是由本发明方法生成共计20701个独立节点。与这两类无网格节点相比，采用常规网格的有限差分法剖分同等大小模型需要40000个节点。将图3(a)和3(b)相比可以看出，界面附近的速度场获得了有效而均匀的剖分，证明了本发明提出的地震波数值模拟方法能够提供更合理的节点分布和更平滑的速度/晶粒半径过渡区，可以避免节点丢失或重叠。图4是图2所示双层介质模型的地震波场快照剖面：其中，图4(a)是利用常规方法所得的500ms的地震波场快照剖面，图4(b)是利用本发明方法所得的500ms的地震波场快照剖面；从图4(a)和4(b)中可以看出，本发明方法所得的地震波场具有更高的精度。图5是图2所示双层介质模型的地震记录剖面：其中，图5(a)是利用常规方法所得的地震记录剖面，图5(b)是利用本发明方法所得的地震记录剖面；从图5(a)和5(b)中也可以得出相同的结论，证明了本发明方法的可行性及有效性。

实例2：

图6是本发明提供的Marmousi-2模型；图7是利用本发明方法所得的图6所示Marmousi-2模型的局部网格剖分图；从图7可看出本发明方法能够灵活剖分，适应于复杂地质构造以及速度剧烈变化的地质条件；图8是图6所示Marmousi-2模型的地震波场快照剖面：其中，图8(a)是利用本发明方法所得的500ms的地震波场快照剖面，图8(b)是利用本发明方法所得的1000ms的地震波场快照剖面，图8(c)是利用本发明方法所得的1500ms的地震波场快照剖面，图8(d)是利用本发明方法所得的2000ms的地震波场快照剖面；从图8中可以看出，本发明方法的地震波场能够在复杂介质中稳定传播，且具有高精度，高性噪比，低存储量的特点。进一步证明了本方面方法的有效性。

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

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

Claims (6)

1.一种地震波数值模拟方法，其特征在于，包括：

获取地震波速度模型；

根据所述地震波速度模型，确定地震波速度模型适应的剖分节点径向半径范围；

根据所述剖分节点径向半径范围，确定所述地震波速度模型中的地震波速度对应的剖分节点径向半径；

根据各所述地震波速度对应的剖分节点径向半径和所述地震波速度模型，采用快速节点剖分法，得到所述地震波速度模型的多个剖分节点、边界内点、权重因子和地震波场传播角；

根据各所述剖分节点采用快速邻域搜索算法，得到各所述剖分节点对应的最邻近的n个剖分节点；

确定径向基函数有限差分地震波数值模拟的最优化形变参数；

根据所述最优化形变参数、各所述地震波速度对应的剖分节点径向半径和各所述剖分节点对应的最邻近的n个剖分节点，确定所有剖分节点对应的差分加权系数；

根据各所述剖分节点、所述边界内点、所述权重因子、所述地震波场传播角和所述所有剖分节点对应的差分加权系数采用地震波场计算公式进行地震波场延拓，确定数值模拟的地震波场；

所述根据各所述地震波速度对应的剖分节点径向半径和所述地震波速度模型，采用快速节点剖分法，得到所述地震波速度模型的多个剖分节点、边界内点、权重因子和地震波场传播角，具体包括：

根据各所述地震波速度对应的剖分节点径向半径和所述地震波速度模型，采用快速节点剖分法进行节点剖分，得到计算域和边界域的剖分节点在已知地震波速度模型中对应的真实空间位置(x,z)；

根据所述真实空间位置(x,z)，从所述地震波速度模型提取相同空间位置处的地震波速度v；

根据所述地震波速度对应的剖分节点径向半径，确定所述地震波速度v对应的剖分节点径向半径r；所述真实空间位置(x,z)、所述地震波速度v和所述剖分节点径向半径r构成所述地震波速度模型的剖分节点(x,z,v,r)；

根据所述剖分节点(x,z,v,r)，确定边界内点(x_inner,z_inner,v_inner,r_inner)；

获取边界域的厚度H_b和边界域内各个剖分节点到计算域最外层剖分节点的垂直距离D_b；

根据所述边界域的厚度H_b和所述边界域内各个剖分节点到计算域最外层剖分节点的垂直距离D_b，确定边界域内各剖分节点的权重因子w；

根据所述边界域内各个剖分节点相对于计算域的空间位置，确定所述边界域内各个剖分节点的地震波场传播角θ。

2.根据权利要求1所述的地震波数值模拟方法，其特征在于，所述根据所述剖分节点径向半径范围，确定所述地震波速度模型中的地震波速度对应的剖分节点径向半径，具体包括：

根据所述剖分节点径向半径范围，确定每个速度值适应的剖分节点径向半径范围的中值；

根据所述速度值适应的剖分节点径向半径范围的中值，确定所述速度值对应的剖分节点径向半径。

3.根据权利要求1所述的地震波数值模拟方法，其特征在于，所述根据所述边界域的厚度H_b和所述边界域内各个剖分节点到计算域最外层剖分节点的垂直距离D_b，确定边界域内各剖分节点的权重因子w，具体包括：

根据所述边界域的厚度H_b和所述边界域内各个剖分节点到计算域最外层剖分节点的垂直距离D_b采用公式w＝D_b/H_b，计算边界域内各剖分节点的权重因子w。

4.根据权利要求1所述的地震波数值模拟方法，其特征在于，所述根据所述边界域内各个剖分节点相对于计算域的空间位置，确定所述边界域内各个剖分节点的地震波场传播角θ，具体包括：

根据边界域内各个剖分节点相对于计算域的空间位置，将所述边界域的各个剖分节点的空间位置进行分类，得到八类位置，所述八类位置分别是上边、左上边、左边、左下边、下边、右下边、右边和右上边；

对所述八类位置分别给于0°、45°、90°、135°、180°、225°、270°和315°的地震波场传播角度，所有的地震波场传播角度构成了边界域内各个剖分节点的地震波场传播角θ。

5.根据权利要求1所述的地震波数值模拟方法，其特征在于，所述根据所述最优化形变参数、各所述地震波速度对应的剖分节点径向半径和各所述剖分节点对应的最邻近的n个剖分节点，确定所有剖分节点对应的差分加权系数，具体包括：

根据所述最优化形变参数和所述径向基函数，得到优化的径向基函数；

根据所述优化的径向基函数、所述地震波速度对应的剖分节点径向半径和所述剖分节点对应的最邻近的n个剖分节点，确定各所述剖分节点对应的最邻近的n个剖分节点的差分加权系数；

根据各所述剖分节点对应的最邻近的n个剖分节点的差分加权系数，得到所有剖分节点对应的差分加权系数；

所述所有剖分节点对应的差分加权系数具体为：

其中，Φ_opt为优化的径向基函数，x(x,z)表示地震波速度模型剖分节点的空间坐标矢量，x表示空间坐标矢量的水平分量，z表示空间坐标矢量的垂直分量，下标“0”表示参与当前地震波速度模型剖分节点差分运算的周围第0个最邻近剖分节点，下标“n”表示参与当前地震波速度模型剖分节点差分运算的周围第n个最邻近剖分节点；b_i为第i个剖分节点周围n个最邻近剖分节点对应的差分加权系数，i＝0,1...,n,....n+5，L[·]表示偏导数算子。

6.根据权利要求1所述的地震波数值模拟方法，其特征在于，所述根据各所述剖分节点、所述边界内点、所述权重因子、所述地震波场传播角和所述所有剖分节点对应的差分加权系数采用地震波场计算公式进行地震波场延拓，确定数值模拟的地震波场，具体包括：

采用公式

确定计算域和边界域所有剖分节点对应的第一类地震波场；

其中，U为每个剖分节点对应的第一类地震波场，t为地震波场传播时间，v为地震波速度，τ为时间步长，b_i为第i个剖分节点周围n个最邻近剖分节点对应的差分加权系数；

确定计算域和边界域所有剖分节点对应的第二类地震波场；

根据所述权重因子对所有剖分节点的所述第一类地震波场和所述第二类地震波场进行加权运算，得到数值模拟的地震波场W；

所述数值模拟的地震波场W具体为：

其中，U为每个剖分节点对应的第一类地震波场，V为每个剖分节点对应的第二类地震波场。

Images