CN106646597A - 基于弹簧网络模型的正演模拟方法及装置 - Google Patents

Abstract

本发明公开一种基于弹簧网络模型的正演模拟方法及装置，具体包括：设定正演模拟采用的时间采样间隔和空间采样间隔；根据所述空间采样间隔对待模拟的地质模型进行离散化，获得数字化地质模型；所述数字化地质模型为矩形网格形式或长方体网格形式；根据所述空间采样间隔确定所述数字化地质模型的各个方向的弹性参数；根据震源子波函数、时间采样间隔及所述数字化地质模型的各个方向的弹性参数，更新所述数字化地质模型的各个时刻的波场值。本发明克服了以往对网格形状的限制，不仅可以适用于正方形或者正方体网格，也适用于矩形或者长方体网格，可以满足实际生产对模拟网格的要求，且能根据实际情况优选时间和空间采样间隔以提高计算效率。

Description

基于弹簧网络模型的正演模拟方法及装置

技术领域

本发明涉及石油地震勘探技术领域，特别涉及一种基于弹簧网络模型的正演模拟方法及装置。

背景技术

地震正演模拟是地震资料采集、处理和地质解释的基础，是降低地震多解性的重要手段。目前，地震波数值模拟的主流方法都是通过近似求解波动方程来实现的，比如常用的有限差分法、有限元法等，其固有缺点是必须假设地质体是相对均匀而连续的单元，以确保描述地震波传播的偏微分方程具有可解性。然而，随着地质体非均匀性的增强，当介质中存在强间断面或者含有多相流体时，这一问题变得更加突出。因此，怎样将实际复杂的地下孔隙介质，用简单理想化且不失一般性的模型来描述，一直是地震波场正演模拟研究的重要课题与方向。

弹簧网络模型作为一种微观力学模型，通过与速度Verlet结合可用于模拟弹性波的传播过程。由于该方法是从最基本的牛顿运动定律和胡克定律出发来进行波场模拟，不会受到波动方程的近似假设条件限制，因而能够适应任意复杂的地下介质模型。

弹簧网络模型(Lattice Spring Model，简称LSM)是一种从微观角度研究介质弹塑性的方法，其相关原理可追溯到20世纪初Born等学者提出的晶体动力学理论。早期的Born弹簧模型，相邻质点之间只考虑中心力(Central-Force)作用，这种力也称为线弹簧作用力。这种模型模拟的介质泊松比固定为0.25，其应用受到限制。1939年，Kirkwood在Born弹簧模型的基础上引入了一种键弯曲作用力(Bond-Bending Force)，即角弹簧作用力，使得模拟材料的泊松比能够在一定范围内变化，该模型因此被称为键弯曲模型(Bond-Bending Model)。1954年，Born和Huang合著的图书《Dynamical Theory of CrystalLattices》比较全面地介绍了这样一种微观模型，为该模型用于研究晶体的弹性、内聚能和晶格动力学等打下了坚实基础。之后，Kittel、Lax和Keating等人对Born弹簧模型作了进一步地发展与完善，保证了形变后晶体的应变能量计算公式具有旋转不变性。

Born弹簧模型在早期主要用于研究介质的弹性性能与断裂行为。1984年，美国学者Grest和Webman提出了将固体离散成一系列由弹簧与节点连接成的立方体元的想法，并将其用于研究渗透集群的震动模式。Kantor和Webman(1984)、Arbabi和Sahimi(1988)先后通过将渗流系统分解成微观的键弯曲模型，研究了其弹性性能。Hassold和Srolovitz(1989)、Murat等(1992)通过Born弹簧模型模拟得到了介质的破裂统计特性。Starzewski等(1996)采用Spring Network Model研究了复合材料和晶体的伸缩性和断裂性。Ladd等(1997)进一步研究了线弹簧与角弹簧的综合模型，并最先使用了“Lattice Spring Model”一词。Buxton等(2001)，Chung等(2002)采用LSM描述了非均匀介质的弹塑性形变和脆性破裂效应。

2000年，爱尔兰的两位学者Toomey和Bean采用一种离散粒子模型，通过将介质离散成微观球体紧密排列的正六边形构造单元，借助虎克定律和速度Verlet算法，实现了泊松体中弹性波的模拟，该方法被称为离散粒子方法(Discrete Particle Scheme，简称DPS)，也有学者称其为弹性网格方法(Elastic Lattice Method，简称ELM)。其后，O’Brien和Bean对相关理论作了进一步推广，解除了泊松比的限制，并在火山震源机制、弹性波模拟等方面得到了应用。

此外，Yim和Sohn(2000)构造了一种新的质点弹簧模型(Mass Spring LatticeModel，简称MSLM)，将传统的弹性波方程中的弹性常数用弹簧的弹性常数替换，采用有限差分方法实现了复杂介质中的超声波模拟。2008年，O’Brien通过在LSM中加入粘弹性弹簧模拟了地震波在粘弹性介质中的传播。Pazdniakou和Adler(2012)对LSM模型的相关理论作了比较全面的介绍，并尝试将其用于模拟干燥孔隙介质中的弹性波。Zhao等(2011)在原始LSM的基础上增加了一种独特的剪切弹簧并应用Distinct Lattice Spring Model(DLSM)研究了介质的弹性和动态破裂性质。Zhu等(2011)借助DLSM模拟了P波在不同岩石交界处的传播过程。2012年，Zhao等对DLSM模型的演化作了简要的概述，并提出了一种采用高阶DLSM模拟介质弹性的方法。同年，Zhao和Khalili(2012)对DLSM的GPU并行加速算法进行了研究。

综上所述，目前基于弹簧网络模型或者类似方法的正演模拟手段，主要集中于正方形网格或者正三角形网格情形，而在实际生产中，为了适应特定的观测系统，常常需要将网格剖分成矩形或长方体，以节省工作量以满足特定方向的计算精度，但现有技术中尚欠缺适应矩形网格或长方体网格的问题研究。

发明内容

本发明实施例提供了一种基于弹簧网络模型的正演模拟方法，将常规的仅仅适用于正方形网格或者正三角形网格的弹簧网络模型，推广到了适用于矩形网格或长方体网格，能够更好地适应实际生产需要，该方法包括：

设定正演模拟采用的时间采样间隔和空间采样间隔；

根据所述空间采样间隔对待模拟的地质模型进行离散化，获得数字化地质模型；所述数字化地质模型为矩形网格形式或长方体网格形式；

根据所述空间采样间隔确定所述数字化地质模型的各个方向的弹性参数；

根据震源子波函数、时间采样间隔及所述数字化地质模型的各个方向的弹性参数，更新所述数字化地质模型的各个时刻的波场值。

本发明实施例提供了一种基于弹簧网络模型的正演模拟装置，将常规的仅仅适用于正方形网格或者正三角形网格的弹簧网络模型，推广到了适用于矩形网格或长方体网格，能够更好地适应实际生产需要，该装置包括：

参数设定模块，用于设定正演模拟采用的时间采样间隔和空间采样间隔；

离散模块，用于根据所述空间采样间隔对待模拟的地质模型进行离散化，获得数字化地质模型；所述数字化地质模型为矩形网格形式或长方体网格形式；

弹性参数确定模块，用于根据所述空间采样间隔确定所述数字化地质模型的各个方向的弹性参数；

波场值确定模块，用于根据震源子波函数、时间采样间隔及所述数字化地质模型的各个方向的弹性参数，更新所述数字化地质模型的各个时刻的波场值。

在本发明实施例中，设定正演模拟采用的时间采样间隔和空间采样间隔；根据所述空间采样间隔对待模拟的地质模型进行离散化，获得数字化地质模型；所述数字化地质模型为矩形网格形式或长方体网格形式；根据所述空间采样间隔确定所述数字化地质模型的各个方向的弹性参数；根据震源子波函数、时间采样间隔及所述数字化地质模型的各个方向的弹性参数，更新所述数字化地质模型的各个时刻的波场值。采用本发明方法可以将常规的仅仅适用于正方形网格或者正三角形网格的弹簧网络模型，推广到了适用于更一般性的矩形网格或长方体网格，这样能够更好地与实际野外地震资料采集观测系统采集的数据进行匹配，满足了实际生产中对于不同方向的计算精度以及计算效率的需求。

附图说明

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

图1为本发明实施例提供的一种基于弹簧网络模型的正演模拟方法的处理流程图；

图2为本发明实施例提供的一种基于弹簧网络模型的正演模拟方法具体处理流程图；

图3为本发明实施例提供的二维(D2Q8)弹簧网络模型示意图；

图4为本发明实施例提供的三维(D3Q18)弹簧网络模型示意图；

图5a)为本发明实施例提供的LSM与FDM在网格比(dz/dx＝5/8)下的P波的归一化相速度频散曲线对比图；

图5b)为本发明实施例提供的LSM与FDM在网格比(dz/dx＝5/8)下的S波的归一化相速度频散曲线对比图；

图5c)为本发明实施例提供的LSM与FDM在网格比(dz/dx＝1)下的P波的归一化相速度频散曲线对比图；

图5d)为本发明实施例提供的LSM与FDM在网格比(dz/dx＝1)下的S波的归一化相速度频散曲线对比图；

图5e)为本发明实施例提供的LSM与FDM在网格比(dz/dx＝8/5)下的P波的归一化相速度频散曲线对比图；

图5f)为本发明实施例提供的LSM与FDM在网格比(dz/dx＝8/5)下的S波的归一化相速度频散曲线对比图；

图6a)为本发明实施例提供的二阶空间精度的FDM模拟得到的波场(v_x)切片图；

图6b)为本发明实施例提供的二阶空间精度的FDM模拟得到的波场(v_z)切片图；

图6c)为本发明实施例提供的四阶空间精度的FDM模拟得到的波场(v_x)切片图；

图6d)为本发明实施例提供的四阶空间精度的FDM模拟得到的波场(v_z)切片图；

图6e)为本发明实施例提供的LSM模拟得到的波场(v_x)切片图；

图6f)为本发明实施例提供的LSM模拟得到的波场(v_z)切片图；

图7a)为本发明实施例提供的二维层状介质情形下的界面之上的接收点处的LSM与二阶及四阶精度FDM的单道地震记录对比图(v_x)；

图7b)为本发明实施例提供的二维层状介质情形下的界面之下的接收点处的LSM与二阶及四阶精度FDM的单道地震记录对比图(v_x)；

图7c)为本发明实施例提供的二维层状介质情形下的界面之上的接收点处的LSM与二阶及四阶精度FDM的单道地震记录对比图(v_z)；

图7d)为本发明实施例提供的二维层状介质情形下的界面之下的接收点处的LSM与二阶及四阶精度FDM的单道地震记录对比图(v_z)；

图8a)为本发明实施例提供的Marmousi模型的FDM模拟得到的波场(v_x)快照图；

图8b)为本发明实施例提供的Marmousi模型的FDM模拟得到的波场(v_z)快照图；

图8c)为本发明实施例提供的Marmousi模型的LSM模拟得到的波场(v_x)快照图；

图8d)为本发明实施例提供的Marmousi模型的LSM模拟得到的波场(v_z)快照图；

图9为本发明实施例提供的一种基于弹簧网络模型的正演模拟装置的结构示意图；

图10为本发明实施例提供的一种基于弹簧网络模型的正演模拟装置的一种具体结构示意图。

具体实施方式

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

现有的基于弹簧网络模型或者类似原理的地震波场正演模拟方法，几乎全都是采用正方形或者正三角形网格，这样的方法在某种意义上不能满足实际野外地震资料采集的需求。基于此，本发明提出一种基于弹簧网络模型的正演模拟方法，以弥补前人提出方法的不足之处，从而更好地与实际地震资料进行匹配，并达到提高数值计算效率的目的。

图1是本发明实施例提供的一种基于弹簧网络模型的正演模拟方法的处理流程图，如图1所示，该方法包括：

步骤101、设定正演模拟采用的时间采样间隔和空间采样间隔。

步骤102、根据所述空间采样间隔对待模拟的地质模型进行离散化，获得数字化地质模型；所述数字化地质模型为矩形网格形式或长方体网格形式。

步骤103、根据所述空间采样间隔确定所述数字化地质模型的各个方向的弹性参数。

步骤104、根据震源子波函数、时间采样间隔及所述数字化地质模型的各个方向的弹性参数，更新所述数字化地质模型的各个时刻的波场值。

具体实施时，在得到一个地质模型对其进行正演模拟时，首先需要设定正演模拟所要采用的时间采样间隔和空间采样间隔(步骤101)，然后需要根据预设的空间采样间隔对待模拟的地质模型(模拟地质区域)进行离散化，获得数字化地质模型(步骤102)。进行离散化的方法是：对于二维(2D)数值模拟问题，将待模拟的地质模型剖分为矩形网格(D2Q8模型，见图3)，或者，对于三维(3D)数值模拟问题，将待模拟的地质模型剖分为长方体网格(D3Q18模型，见图4)，二维及三维模型的网格边长可以相等，也可以不相等。然后将待模拟的地质模型离散成由离散网格节点与线弹簧构成的弹性元，相邻离散网格节点之间通过线弹簧连接；相邻离散网格节点之间的距离为预设的空间采样间隔。矩形网格(D2Q8模型)形式的数字化地质模型由弹性面元构成；长方体网格(D3Q18模型)形式的数字化地质模型由弹性体元构成。

当然，本发明方法对正方形网格或者正方体网格同样适用。

具体实施时，在获得数字化地质模型后，首先需要根据预设的空间采样间隔求取数字化地质模型的各个方向的弹性参数。假设，对数字化地质模型(离散模型)施加外力或者扰动后，相应位置的线弹簧会发生形变，弹性面元(2D)或者弹性体元(3D)的中心节点i受到的线弹簧的总作用力F_i可以按照如下公式计算：

其中，φ_i为第i个弹性元的能量密度，k_j为中心节点i与相邻节点j之间的线弹簧的弹性常数，u_i为节点的位移，x_ij为节点i指向节点j的向量，为归一化的方向向量，∑代表对中心节点的所有邻点求和，符号“·”代表向量内积运算；

对于二维D2Q8模型，上式中，E为弹性面元的面积，N＝8；

对于三维D3Q18模型，上式中，E为弹性体元的体积，N＝18。

线弹簧的弹性参数是一个与弹簧长度有关的常数，在离散模型中具体表现为与邻点位置或者方向有关的常数，在计算弹性体元或者弹性面元时，用到了邻近的8个节点(D2Q8模型)或者18个节点(D3Q18模型)上的波场值，且不同邻点与中心节点之间的线弹簧的弹性常数不相等。假设相等长度的线弹簧具有相同的弹性参数时，其计算公式如下：

①对于二维D2Q8模型，各个方向弹性参数的计算公式如下：

且有，

其中，Δx和Δz分别为沿着x轴和z轴方向的网格长度，ρ为介质的质量密度，V_P为纵波传播速度，k₁₀、k₂₀和k₃₀分别为以中心节点为顶点，沿着x轴正负方向、z轴正负方向以及对角正负方向的线弹簧弹性参数，如图3所示。

②对于三维D3Q18模型，各个方向弹性常数的计算公式如下：

且有，

其中，Δx、Δy和Δz分别为沿着x轴、y轴和z轴方向的网格长度，ρ为介质的质量密度，V_P为纵波传播速度，k₁₀₀、k₂₀₀和k₃₀₀分别为以中心节点为顶点，沿着x轴、y轴和z轴方向的线弹簧弹性参数，k₄₀₀、k₅₀₀和k₆₀₀分别代表沿着xoy、xoz和yoz坐标平面对角方向的线弹簧弹性参数，如图4所示。

具体实施时，计算完模拟区域内各个节点的总的作用力后，根据震源子波函数、时间采样间隔及数字化地质模型的各个方向的弹性参数，利用速度Verlet算法更新数字化地质模型的各个时刻的波场值(步骤104)。该方法不是基于有限差分方法(FDM)或者有限元法(FEM)中常用的弹性波动方程，而是依据最基本的牛顿运动定律和胡克定律来模拟弹性波在固体介质中的传播，具体实施过程中，引用流体力学中经常用到的速度Verlet算法，用于更新模拟区域内各个时刻的波场值，波场值包括位移x_i、速度v_i和加速度a_i：

其中，Δt为模拟采用的时间采样间隔，F_i为中心节点i所受到的线弹簧的总作用力，m_i为中心节点i的质量，χ为粘滞项系数，用于衰减节点的震动能量，这里取χ＝0。

上述速度Verlet算法是一个时间二阶计算精度的数值算法，而弹簧网络模型在数值模拟时的空间计算精度取决于中心节点i所受作用力的计算精度。

具体实施时，如图2所示，该方法还包括：

步骤105：根据所述时间采样间隔和空间采样间隔确定所述数字化地质模型的相速度频散曲线。

具体的，在二维数值模拟情况下，针对一般性的矩形网格(Δx≠Δz)，数字化地质模型的P波和S波的相速度频散曲线按照如下公式确定：

中间变量计算公式如下：

上式中

且有

其中，q_p为LSM模拟的波场对应的P波的相速度频散曲线，q_s为LSM模拟的波场对应的S波的相速度频散曲线，r＝Δz/Δx为z轴方向与x轴方向上的网格边长之比，Δt为时间采样间隔，V_P和V_S分别为纵波和横波速度，θ为平面波与x轴正方向的夹角，λ为波数。

需要指出的是，上述针对一般性矩形网格的相速度频散曲线，对于正方形网格(Δx＝Δz)情形同样适用。

倘若仅仅考虑线弹簧作用力时，纵横波速度之比(V_P/V_S)为也就是模拟介质的泊松比在三维情况下为0.25。

具体实施时，如图2所示，该方法还包括：

步骤106：对所述时间采样间隔和空间采样间隔进行验证，确定时间采样间隔和空间采样间隔的选择是否合理。具体的，根据弹簧网络模型的稳定性条件和相速度频散曲线，验证时间采样间隔和空间采样间隔是否合理。当所述时间采样间隔不满足弹簧网络模型的稳定性条件，且所述空间采样间隔不满足相速度频散精度要求时，重新设定时间采样间隔和空间采样间隔，对重新设定的时间采样间隔和空间采样间隔进行再次验证，直到重新设定的时间采样间隔满足弹簧网络模型的稳定性条件，且重新设定的空间采样间隔满足相速度频散精度要求为止。

比如，验证时间采样间隔和空间采样间隔选择过大或过小，则重新设定时将时间采样间隔和空间采样间隔减小或增大，然后再次验证，直到时间采样间隔和空间采样间隔满足数值计算稳定性条件和数值频散精度要求。

在实际数值模拟过程中，弹簧网络模型(LSM)的稳定性条件为：

Δt＜Δd_min/V_max (14)

其中，Δt为时间采样间隔，Δd_min为最小的空间采样间隔，V_max为最大传播速度。

数值模拟中，该方法(LSM)的稳定性条件(Δt＜Δd_min/V_max)比FDM(对于2D情形：对于3D情形：更宽松，与FEM的稳定性条件(Δt＜Δd_min/V_max)在同一个水平。

最后，输出并保存计算获得的各个时刻的波场值。

下面通过比较来说明本方法的优点。

图5a)至5f)显示的是LSM与时间二阶、空间二阶精度的FDM在不同网格比下的P波或S波的归一化相速度频散曲线对比图。其中，不同网格比分别为：dz/dx＝5/8、dz/dx＝1、dz/dx＝8/5。从图5c)和5d)中可看出，当采用正方形网格(Δx＝Δz)时，两种方法的频散曲线很接近；从图5a)和5b)、5e)和5f)中可看出，而当采用矩形网格(Δx≠Δz即dz/dx＝5/8或dz/dx＝8/5)时，两种方法的数值频散特性各有优劣，总体差别不大。

接下来通过一个二维双层介质模型来验证本发明方法的有效性。模型上层与下层纵波速度分别为4000m/s与5000m/s，而横波速度取时间采样间隔Δt＝0.5ms；x轴方向与z轴方向的空间采样间隔取值为：Δx＝5m,Δz＝8m；模型大小为：5600m×5000m；震源采用主频为15Hz雷克子波。针对相同的模型参数，采用LSM、二阶(FDM-2nd)与四阶(FDM-4th)空间精度的FDM模拟得到的波场切片对比图见图6a)至6f)，相应的单道地震记录见图7a)至7d)。

根据图6a)至6f)和图7a)至7d)，不难看出，对于二维层状介质模型，LSM模拟得到的地震波场切片和单道地震记录与二阶及四阶精度的FDM吻合很好，两者的误差很小，验证了该方法的正确性。

为了更进一步说明本方法对实际复杂介质模型的适用性，我们采用Marmousi模型进行实验测试。模拟采用的时间采样间隔为0.5ms；x轴方向与z轴方向的空间采样间隔取值为：Δx＝5m,Δz＝7m；15Hz主频的雷子子波震源加载在模型中间，其坐标为(1915m，700m)。为了作对比，采用相同参数的FDM也进行了相应的波场模拟。其中，为了减弱模拟区域的边界反射波，FDM采用的是CPML吸收边界，而LSM采用的是传统的指数衰减吸收边界(ABC)，模拟得到的波场快照见图8a)至8d)。通过对比，可发现两种方法模拟得到的复杂介质模型的波场快照吻合，这说明了LSM能够适应实际复杂地质模型的弹性波正演工作。

基于同一发明构思，本发明实施例中还提供了一种基于弹簧网络模型的正演模拟装置，如下面的实施例所述。由于基于弹簧网络模型的正演模拟装置解决问题的原理与基于弹簧网络模型的正演模拟方法相似，因此基于弹簧网络模型的正演模拟装置的实施可以参见基于弹簧网络模型的正演模拟方法的实施，重复之处不再赘述。以下所使用的，术语“单元”或者“模块”可以实现预定功能的软件和/或硬件的组合。尽管以下实施例所描述的装置较佳地以软件来实现，但是硬件，或者软件和硬件的组合的实现也是可能并被构想的。

图9是本发明实施例的基于弹簧网络模型的正演模拟装置的一种结构框图，如图9所示，包括：

参数设定模块901，用于设定正演模拟采用的时间采样间隔和空间采样间隔；

离散模块902，用于根据所述空间采样间隔对待模拟的地质模型进行离散化，获得数字化地质模型；所述数字化地质模型为矩形网格形式或长方体网格形式；

弹性参数确定模块903，用于根据所述空间采样间隔确定所述数字化地质模型的各个方向的弹性参数；

波场值确定模块904，用于根据震源子波函数、时间采样间隔及所述数字化地质模型的各个方向的弹性参数，更新所述数字化地质模型的各个时刻的波场值。

下面对该结构进行说明。

具体实施时，离散模块902具体用于：

将待模拟的地质模型离散成由离散网格节点与线弹簧构成的弹性元，相邻离散网格节点之间通过线弹簧连接；相邻离散网格节点之间的距离为所述空间采样间隔；

矩形网格形式的数字化地质模型由弹性面元构成；

长方体网格形式的数字化地质模型由弹性体元构成。

具体实施时，所述弹性参数确定模块903具体用于：

对所述数字化地质模型施加外力或者扰动，相应位置的线弹簧发生形变，弹性元的中心节点i受到的线弹簧的总作用力F_i按照如下公式确定：

其中，φ_i为第i个弹性元的能量密度，k_j为中心节点i与相邻节点j之间的线弹簧的弹性常数，u_i为中心节点i的位移，u_j为相邻节点j的位移，x_ij为中心节点i指向相邻节点j的向量，为中心节点i指向相邻节点j的归一化的方向向量，∑代表对中心节点i的所有相邻节点j求和，符号“·”代表向量内积运算，N为相邻节点j的个数；中心节点的个数和弹性元的个数相同；

对于矩形网格形式的数字化地质模型，E为弹性面元的面积，N＝8；

对于长方体网格形式的数字化地质模型，E为弹性体元的体积，N＝18；

对于矩形网格形式的数字化地质模型，各个方向的弹性参数按如下公式确定：

其中，Δx和Δz分别为沿着x轴和z轴方向的网格长度，ρ为介质的质量密度，V_P为纵波传播速度，k₁₀、k₂₀和k₃₀分别为以中心节点为顶点，沿着x轴方向、z轴方向以及对角方向的线弹簧弹性常数；

对于长方体网格形式的数字化地质模型，各个方向的弹性参数按如下公式确定：

其中，Δx、Δy和Δz分别为沿着x轴、y轴和z轴方向的网格长度，ρ为介质的质量密度，V_P为纵波传播速度，k₁₀₀、k₂₀₀和k₃₀₀分别为以中心节点为顶点，沿着x轴、y轴和z轴方向的线弹簧弹性常数，k₄₀₀、k₅₀₀和k₆₀₀分别代表沿着xoy、xoz和yoz坐标平面对角方向的线弹簧弹性常数。

具体实施时，所述波场值确定模块904具体用于：

利用速度Verlet算法按如下公式更新所述数字化地质模型的各个时刻的波场值：

其中，x_i为位移；v_i为速度；a_i为加速度；Δt为时间采样间隔，F_i为中心节点i所受到的线弹簧的总作用力，m_i为中心节点i的质量，χ为粘滞项系数，χ＝0。

具体实施时，如图10所示，该装置还包括：相速度数值频散分析模块905，用于根据所述时间采样间隔和空间采样间隔确定所述数字化地质模型的相速度频散曲线；

所述相速度数值频散分析模块905具体用于：

按如下公式确定离散格式的弹簧网络模型的相速度频散曲线按照如下公式确定：

中间变量计算公式如下：

上式中

且有

其中，q_p为LSM模拟的波场对应的P波的相速度频散曲线，q_s为LSM模拟的波场对应的S波的相速度频散曲线，r＝Δz/Δx为z轴方向与x轴方向上的网格边长之比，Δt为时间采样间隔，V_P和V_S分别为纵波和横波速度，θ为平面波与x轴正方向的夹角，λ为波数。

具体实施时，如图10所示，该装置还包括：验证模块906，用于对所述时间采样间隔和空间采样间隔进行验证，当所述时间采样间隔不满足弹簧网络模型的稳定性条件，且所述空间采样间隔不满足相速度频散精度要求时，重新设定时间采样间隔和空间采样间隔，对重新设定的时间采样间隔和空间采样间隔进行再次验证，直到重新设定的时间采样间隔满足弹簧网络模型的稳定性条件，且重新设定的空间采样间隔满足相速度频散精度要求为止。

具体实施时，所述弹簧网络模型的稳定性条件为：

Δt＜Δd_min/V_max；

其中，Δt为时间采样间隔，Δd_min为空间采样间隔的最小值，V_max为传播速度的最大值。

综上所述，本发明实施例具有以下优点：

1.本方法在模拟过程中不依赖于传统的弹性波方程，不受其相关假设条件限制，能够适应实际复杂介质模型。

2.依据稳定性条件和相速度频散曲线优选的时间与空间采样间隔等模拟参数，保证了模拟的波场的正确性。

3.鉴于本发明提出的正演方法可以适用于矩形网格，具有一定的网格选择灵活性，因而能够满足实际生产中对于不同方向的计算精度以及计算效率的要求。

本领域内的技术人员应明白，本发明的实施例可提供为方法、系统、或计算机程序产品。因此，本发明可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本发明可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。

本发明是参照根据本发明实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。

这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的步骤。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明实施例可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (14)

1.一种基于弹簧网络模型的正演模拟方法，其特征在于，包括：

设定正演模拟采用的时间采样间隔和空间采样间隔；

根据所述空间采样间隔对待模拟的地质模型进行离散化，获得数字化地质模型；所述数字化地质模型为矩形网格形式或长方体网格形式；

根据所述空间采样间隔确定所述数字化地质模型的各个方向的弹性参数；

根据震源子波函数、时间采样间隔及所述数字化地质模型的各个方向的弹性参数，更新所述数字化地质模型的各个时刻的波场值。

2.如权利要求1所述的基于弹簧网络模型的正演模拟方法，其特征在于，根据所述空间采样间隔对待模拟的地质模型进行离散化，包括：

将待模拟的地质模型离散成由离散网格节点与线弹簧构成的弹性元，相邻离散网格节点之间通过线弹簧连接；相邻离散网格节点之间的距离为所述空间采样间隔；

矩形网格形式的数字化地质模型由弹性面元构成；

长方体网格形式的数字化地质模型由弹性体元构成。

3.如权利要求2所述的基于弹簧网络模型的正演模拟方法，其特征在于，根据所述空间采样间隔确定所述数字化地质模型的各个方向的弹性参数，包括：

对所述数字化地质模型施加外力或者扰动，相应位置的线弹簧发生形变，弹性元的中心节点i受到的线弹簧的总作用力F_i按照如下公式确定：

$F_i=\frac{\∂{\mathrm{\φ}}_i}{\∂u_i}=\frac{1}{E}\underset{j=1}{\overset{N}{\Σ}}{k}_j\[(u_j-u_i)\·{\hat{x}}_{ij}\]{\hat{x}}_{ij};$

其中，φ_i为第i个弹性元的能量密度，k_j为中心节点i与相邻节点j之间的线弹簧的弹性常数，u_i为中心节点i的位移，u_j为相邻节点j的位移，x_ij为中心节点i指向相邻节点j的向量，为中心节点i指向相邻节点j的归一化的方向向量，∑代表对中心节点i的所有相邻节点j求和，符号“·”代表向量内积运算，N为相邻节点j的个数；中心节点的个数和弹性元的个数相同；

对于矩形网格形式的数字化地质模型，E为弹性面元的面积，N＝8；

对于长方体网格形式的数字化地质模型，E为弹性体元的体积，N＝18；

对于矩形网格形式的数字化地质模型，各个方向的弹性参数按如下公式确定：

$\left\{\begin{array}{c}{k}_{10}=\left(\frac{\mathrm{\Δ}z}{\mathrm{\Δ}x}-\frac{\mathrm{\Δ}x}{3\mathrm{\Δ}z}\right){\mathrm{\ρV}}_P^2\\ k_{20}=\left(\frac{\mathrm{\Δ}x}{\mathrm{\Δ}z}-\frac{\mathrm{\Δ}z}{3\mathrm{\Δ}x}\right){\mathrm{\ρV}}_P^2\\ k_{30}=\frac{{\mathrm{\Δx}}^2+{\mathrm{\Δz}}^2}{6\mathrm{\Δ}x\mathrm{\Δ}z}{\mathrm{\ρV}}_P^2\end{array}\right.;$

其中，Δx和Δz分别为沿着x轴和z轴方向的网格长度，ρ为介质的质量密度，V_P为纵波传播速度，k₁₀、k₂₀和k₃₀分别为以中心节点为顶点，沿着x轴方向、z轴方向以及对角方向的线弹簧弹性常数；

对于长方体网格形式的数字化地质模型，各个方向的弹性参数按如下公式确定：

$\left\{\begin{array}{c}{k}_{100}=\left(\frac{\mathrm{\Δ}y\mathrm{\Δ}z}{\mathrm{\Δ}x}-\frac{\mathrm{\Δ}x\mathrm{\Δ}z}{3\mathrm{\Δ}y}-\frac{\mathrm{\Δ}x\mathrm{\Δ}y}{3\mathrm{\Δ}z}\right){\mathrm{\ρV}}_P^2\\ k_{200}=\left(\frac{\mathrm{\Δ}x\mathrm{\Δ}z}{\mathrm{\Δ}y}-\frac{\mathrm{\Δ}y\mathrm{\Δ}z}{3\mathrm{\Δ}x}-\frac{\mathrm{\Δ}x\mathrm{\Δ}y}{3\mathrm{\Δ}z}\right){\mathrm{\ρV}}_P^2\\ k_{300}=\left(\frac{\mathrm{\Δ}x\mathrm{\Δ}y}{\mathrm{\Δ}z}-\frac{\mathrm{\Δ}y\mathrm{\Δ}z}{3\mathrm{\Δ}x}-\frac{\mathrm{\Δ}x\mathrm{\Δ}z}{3\mathrm{\Δ}y}\right){\mathrm{\ρV}}_P^2\\ k_{400}=\frac{\left({\mathrm{\Δx}}^2+{\mathrm{\Δy}}^2\right)\mathrm{\Δ}z}{6\mathrm{\Δ}x\mathrm{\Δ}y}{\mathrm{\ρV}}_P^2\\ k_{500}=\frac{\left({\mathrm{\Δx}}^2+{\mathrm{\Δz}}^2\right)\mathrm{\Δ}y}{6\mathrm{\Δ}x\mathrm{\Δ}z}{\mathrm{\ρV}}_P^2\\ k_{600}=\frac{\left({\mathrm{\Δy}}^2+{\mathrm{\Δz}}^2\right)\mathrm{\Δ}x}{6\mathrm{\Δ}y\mathrm{\Δ}z}{\mathrm{\ρV}}_P^2\end{array}\right.;$

其中，Δx、Δy和Δz分别为沿着x轴、y轴和z轴方向的网格长度，ρ为介质的质量密度，V_P为纵波传播速度，k₁₀₀、k₂₀₀和k₃₀₀分别为以中心节点为顶点，沿着x轴、y轴和z轴方向的线弹簧弹性常数，k₄₀₀、k₅₀₀和k₆₀₀分别代表沿着xoy、xoz和yoz坐标平面对角方向的线弹簧弹性常数。

4.如权利要求3所述的基于弹簧网络模型的正演模拟方法，其特征在于，根据震源子波函数、时间采样间隔及所述数字化地质模型的各个方向的弹性参数，更新所述数字化地质模型的各个时刻的波场值，包括：

利用速度Verlet算法按如下公式更新所述数字化地质模型的各个时刻的波场值：

$x_i(t+\mathrm{\Δ}t)=x_i\left(t\right)+v_i\left(t\right)\mathrm{\Δ}t+\frac{a_i\left(t\right)}{2}{\mathrm{\Δt}}^2;$

$v_i(t+\frac{\mathrm{\Δ}t}{2})=v_i\left(t\right)+\frac{a_i\left(t\right)}{2}\mathrm{\Δ}t;$

$a_i(t+\mathrm{\Δ}t)=\frac{F_i(t+\mathrm{\Δ}t)}{m_i}-{\mathrm{\χv}}_i(t+\frac{\mathrm{\Δ}t}{2});$

$v_i(t+\mathrm{\Δ}t)=v_i(t+\frac{\mathrm{\Δ}t}{2})+\frac{a_i(t+\mathrm{\Δ}t)}{2}\mathrm{\Δ}t;$

其中，x_i为位移；v_i为速度；a_i为加速度；Δt为时间采样间隔，F_i为中心节点i所受到的线弹簧的总作用力，m_i为中心节点i的质量，χ为粘滞项系数，χ＝0。

5.如权利要求1所述的基于弹簧网络模型的正演模拟方法，其特征在于，还包括：根据所述时间采样间隔和空间采样间隔确定所述数字化地质模型的相速度频散曲线；

按照如下公式确定所述数字化地质模型的相速度频散曲线：

$\left\{\begin{array}{c}{q}_P=\frac{\sqrt{2}}{\mathrm{\π}\mathrm{\γ}H}\arcsin \left(\frac{\mathrm{\γ}}{2}\sqrt{\frac{\mathrm{\α}+\mathrm{\β}}{r}}\right)\\ q_S=\frac{\sqrt{2}}{\mathrm{\π}\mathrm{\γ}H}\frac{V_P}{V_S}\arcsin \left(\frac{\mathrm{\γ}}{2}\sqrt{\frac{\mathrm{\α}-\mathrm{\β}}{r}}\right)\end{array}\right.;$

中间变量计算公式如下：

$\left\{\begin{array}{c}\mathrm{\α}=\[A(\frac{1}{r}-\frac{r}{3})+B(r-\frac{1}{3r})\]+\frac{1}{3}(BC+AD)(r+\frac{1}{r})\\ \mathrm{\β}=\sqrt{{\{\[A(\frac{1}{r}-\frac{r}{3})-B(r-\frac{1}{3r})\]+\frac{1}{3}(BC+AD)(r-\frac{1}{r})\}}^2+\frac{16}{9}ABCD}\end{array}\right.;$

上式中

$\left\{\begin{array}{c}A={\sin }^2\left(\mathrm{\π}rHsin\mathrm{\θ}\right)\\ B={\sin }^2\left(\mathrm{\π}H\cos \mathrm{\θ}\right)\\ C={\cos }^2\left(\mathrm{\π}rHsin\mathrm{\θ}\right)\\ D={\cos }^2\left(\mathrm{\π}H\cos \mathrm{\θ}\right)\end{array}\right.;$

且有

$\left\{\begin{array}{c}\mathrm{\γ}=\sqrt{2}{V}_P\mathrm{\Δ}t/\mathrm{\Δ}x\\ H=\mathrm{\Δ}x/\mathrm{\λ}\end{array}\right.;$

其中，q_p为LSM模拟的波场对应的P波的相速度频散曲线，q_s为LSM模拟的波场对应的S波的相速度频散曲线，r＝Δz/Δx为z轴方向与x轴方向上的网格边长之比，Δt为时间采样间隔，V_P和V_S分别为纵波和横波速度，θ为平面波与x轴正方向的夹角，λ为波长。

6.如权利要求5所述的基于弹簧网络模型的正演模拟方法，其特征在于，还包括：对所述时间采样间隔和空间采样间隔进行验证，当所述时间采样间隔不满足弹簧网络模型的稳定性条件，且所述空间采样间隔不满足相速度频散精度要求时，重新设定时间采样间隔和空间采样间隔，对重新设定的时间采样间隔和空间采样间隔进行再次验证，直到重新设定的时间采样间隔满足弹簧网络模型的稳定性条件，且重新设定的空间采样间隔满足相速度频散精度要求为止。

7.如权利要求6所述的基于弹簧网络模型的正演模拟方法，其特征在于，所述弹簧网络模型的稳定性条件为：

Δt＜Δd_min/V_max；

其中，Δt为时间采样间隔，Δd_min为空间采样间隔的最小值，V_max为传播速度的最大值。

8.一种基于弹簧网络模型的正演模拟装置，其特征在于，包括：

参数设定模块，用于设定正演模拟采用的时间采样间隔和空间采样间隔；

离散模块，用于根据所述空间采样间隔对待模拟的地质模型进行离散化，获得数字化地质模型；所述数字化地质模型为矩形网格形式或长方体网格形式；

弹性参数确定模块，用于根据所述空间采样间隔确定所述数字化地质模型的各个方向的弹性参数；

波场值确定模块，用于根据震源子波函数、时间采样间隔及所述数字化地质模型的各个方向的弹性参数，更新所述数字化地质模型的各个时刻的波场值。

9.如权利要求8所述的基于弹簧网络模型的正演模拟装置，其特征在于，所述离散模块具体用于：

将待模拟的地质模型离散成由离散网格节点与线弹簧构成的弹性元，相邻离散网格节点之间通过线弹簧连接；相邻离散网格节点之间的距离为所述空间采样间隔；

矩形网格形式的数字化地质模型由弹性面元构成；

长方体网格形式的数字化地质模型由弹性体元构成。

10.如权利要求9所述的基于弹簧网络模型的正演模拟装置，其特征在于，所述弹性参数确定模块具体用于：

对所述数字化地质模型施加外力或者扰动，相应位置的线弹簧发生形变，弹性元的中心节点i受到的线弹簧的总作用力F_i按照如下公式确定：

$F_i=\frac{\∂{\mathrm{\φ}}_i}{\∂u_i}=\frac{1}{E}\underset{j=1}{\overset{N}{\Σ}}{k}_j\[(u_j-u_i)\·{\hat{x}}_{ij}\]{\hat{x}}_{ij};$

其中，φ_i为第i个弹性元的能量密度，k_j为中心节点i与相邻节点j之间的线弹簧的弹性常数，u_i为中心节点i的位移，u_j为相邻节点j的位移，x_ij为中心节点i指向相邻节点j的向量，为中心节点i指向相邻节点j的归一化的方向向量，∑代表对中心节点i的所有相邻节点j求和，符号“·”代表向量内积运算，N为相邻节点j的个数；中心节点的个数和弹性元的个数相同；

对于矩形网格形式的数字化地质模型，E为弹性面元的面积，N＝8；

对于长方体网格形式的数字化地质模型，E为弹性体元的体积，N＝18；

对于矩形网格形式的数字化地质模型，各个方向的弹性参数按如下公式确定：

$\left\{\begin{array}{c}{k}_{10}=\left(\frac{\mathrm{\Δ}z}{\mathrm{\Δ}x}-\frac{\mathrm{\Δ}x}{3\mathrm{\Δ}z}\right){\mathrm{\ρV}}_P^2\\ k_{20}=\left(\frac{\mathrm{\Δ}x}{\mathrm{\Δ}z}-\frac{\mathrm{\Δ}z}{3\mathrm{\Δ}x}\right){\mathrm{\ρV}}_P^2\\ k_{30}=\frac{{\mathrm{\Δx}}^2+{\mathrm{\Δz}}^2}{6\mathrm{\Δ}x\mathrm{\Δ}z}{\mathrm{\ρV}}_P^2\end{array}\right.;$

其中，Δx和Δz分别为沿着x轴和z轴方向的网格长度，ρ为介质的质量密度，V_P为纵波传播速度，k₁₀、k₂₀和k₃₀分别为以中心节点为顶点，沿着x轴方向、z轴方向以及对角方向的线弹簧弹性常数；

对于长方体网格形式的数字化地质模型，各个方向的弹性参数按如下公式确定：

$\left\{\begin{array}{c}{k}_{100}=\left(\frac{\mathrm{\Δ}y\mathrm{\Δ}z}{\mathrm{\Δ}x}-\frac{\mathrm{\Δ}x\mathrm{\Δ}z}{3\mathrm{\Δ}y}-\frac{\mathrm{\Δ}x\mathrm{\Δ}y}{3\mathrm{\Δ}z}\right){\mathrm{\ρV}}_P^2\\ k_{200}=\left(\frac{\mathrm{\Δ}x\mathrm{\Δ}z}{\mathrm{\Δ}y}-\frac{\mathrm{\Δ}y\mathrm{\Δ}z}{3\mathrm{\Δ}x}-\frac{\mathrm{\Δ}x\mathrm{\Δ}y}{3\mathrm{\Δ}z}\right){\mathrm{\ρV}}_P^2\\ k_{300}=\left(\frac{\mathrm{\Δ}x\mathrm{\Δ}y}{\mathrm{\Δ}z}-\frac{\mathrm{\Δ}y\mathrm{\Δ}z}{3\mathrm{\Δ}x}-\frac{\mathrm{\Δ}x\mathrm{\Δ}z}{3\mathrm{\Δ}y}\right){\mathrm{\ρV}}_P^2\\ k_{400}=\frac{\left({\mathrm{\Δx}}^2+{\mathrm{\Δy}}^2\right)\mathrm{\Δ}z}{6\mathrm{\Δ}x\mathrm{\Δ}y}{\mathrm{\ρV}}_P^2\\ k_{500}=\frac{\left({\mathrm{\Δx}}^2+{\mathrm{\Δz}}^2\right)\mathrm{\Δ}y}{6\mathrm{\Δ}x\mathrm{\Δ}z}{\mathrm{\ρV}}_P^2\\ k_{600}=\frac{\left({\mathrm{\Δy}}^2+{\mathrm{\Δz}}^2\right)\mathrm{\Δ}x}{6\mathrm{\Δ}y\mathrm{\Δ}z}{\mathrm{\ρV}}_P^2\end{array}\right.;$

其中，Δx、Δy和Δz分别为沿着x轴、y轴和z轴方向的网格长度，ρ为介质的质量密度，V_P为纵波传播速度，k₁₀₀、k₂₀₀和k₃₀₀分别为以中心节点为顶点，沿着x轴、y轴和z轴方向的线弹簧弹性常数，k₄₀₀、k₅₀₀和k₆₀₀分别代表沿着xoy、xoz和yoz坐标平面对角方向的线弹簧弹性常数。

11.如权利要求10所述的基于弹簧网络模型的正演模拟装置，其特征在于，所述波场值确定模块具体用于：

利用速度Verlet算法按如下公式更新所述数字化地质模型的各个时刻的波场值：

$x_i(t+\mathrm{\Δ}t)=x_i\left(t\right)+v_i\left(t\right)\mathrm{\Δ}t+\frac{a_i\left(t\right)}{2}{\mathrm{\Δt}}^2;$

$v_i(t+\frac{\mathrm{\Δ}t}{2})=v_i\left(t\right)+\frac{a_i\left(t\right)}{2}\mathrm{\Δ}t;$

$a_i(t+\mathrm{\Δ}t)=\frac{F_i(t+\mathrm{\Δ}t)}{m_i}-{\mathrm{\χv}}_i(t+\frac{\mathrm{\Δ}t}{2});$

$v_i(t+\mathrm{\Δ}t)=v_i(t+\frac{\mathrm{\Δ}t}{2})+\frac{a_i(t+\mathrm{\Δ}t)}{2}\mathrm{\Δ}t;$

其中，x_i为位移；v_i为速度；a_i为加速度；Δt为时间采样间隔，F_i为中心节点i所受到的线弹簧的总作用力，m_i为中心节点i的质量，χ为粘滞项系数，χ＝0。

12.如权利要求8所述的基于弹簧网络模型的正演模拟装置，其特征在于，还包括：相速度数值频散分析模块，用于根据所述时间采样间隔和空间采样间隔确定所述数字化地质模型的相速度频散曲线；

所述相速度数值频散分析模块具体用于：

按如下公式确定离散格式的弹簧网络模型的相速度频散曲线按照如下公式确定：

$\left\{\begin{array}{c}{q}_P=\frac{\sqrt{2}}{\mathrm{\π}\mathrm{\γ}H}\arcsin \left(\frac{\mathrm{\γ}}{2}\sqrt{\frac{\mathrm{\α}+\mathrm{\β}}{r}}\right)\\ q_S=\frac{\sqrt{2}}{\mathrm{\π}\mathrm{\γ}H}\frac{V_P}{V_S}\arcsin \left(\frac{\mathrm{\γ}}{2}\sqrt{\frac{\mathrm{\α}-\mathrm{\β}}{r}}\right)\end{array}\right.;$

中间变量计算公式如下：

$\left\{\begin{array}{c}\mathrm{\α}=\[A(\frac{1}{r}-\frac{r}{3})+B(r-\frac{1}{3r})\]+\frac{1}{3}(BC+AD)(r+\frac{1}{r})\\ \mathrm{\β}=\sqrt{{\{\[A(\frac{1}{r}-\frac{r}{3})-B(r-\frac{1}{3r})\]+\frac{1}{3}(BC+AD)(r-\frac{1}{r})\}}^2+\frac{16}{9}ABCD}\end{array}\right.;$

上式中

$\left\{\begin{array}{c}A={\sin }^2\left(\mathrm{\π}rHsin\mathrm{\θ}\right)\\ B={\sin }^2\left(\mathrm{\π}H\cos \mathrm{\θ}\right)\\ C={\cos }^2\left(\mathrm{\π}rHsin\mathrm{\θ}\right)\\ D={\cos }^2\left(\mathrm{\π}H\cos \mathrm{\θ}\right)\end{array}\right.;$

且有

$\left\{\begin{array}{c}\mathrm{\γ}=\sqrt{2}{V}_P\mathrm{\Δ}t/\mathrm{\Δ}x\\ H=\mathrm{\Δ}x/\mathrm{\λ}\end{array}\right.;$

其中，q_p为LSM模拟的波场对应的P波的相速度频散曲线，q_s为LSM模拟的波场对应的S波的相速度频散曲线，r＝Δz/Δx为z轴方向与x轴方向上的网格边长之比，Δt为时间采样间隔，V_P和V_S分别为纵波和横波速度，θ为平面波与x轴正方向的夹角，λ为波长。

13.如权利要求12所述的基于弹簧网络模型的正演模拟装置，其特征在于，还包括：验证模块，用于对所述时间采样间隔和空间采样间隔进行验证，当所述时间采样间隔不满足弹簧网络模型的稳定性条件，且所述空间采样间隔不满足相速度频散精度要求时，重新设定时间采样间隔和空间采样间隔，对重新设定的时间采样间隔和空间采样间隔进行再次验证，直到重新设定的时间采样间隔满足弹簧网络模型的稳定性条件，且重新设定的空间采样间隔满足相速度频散精度要求为止。

14.如权利要求13所述的基于弹簧网络模型的正演模拟装置，其特征在于，所述弹簧网络模型的稳定性条件为：

Δt＜Δd_min/V_max；

其中，Δt为时间采样间隔，Δd_min为空间采样间隔的最小值，V_max为传播速度的最大值。