标准线性固体模型的稳定性条件数值解的计算方法及系统

技术领域

本发明涉及地震勘探基础应用技术领域，尤其涉及一种标准线性固体模型 的有限差分解的稳定性条件数值解的计算方法及系统。

背景技术

地震数值模拟是地震勘探和地震学的重要基础，同时也是了解复杂介质中 地震波传播规律的重要工具，其作用贯穿于整个地震采集、处理和解释中。随 着地震勘探开发的深入，常规的弹性介质理论难以满足实际介质需求。地震波 在实际地层中传播时，能量和相位都发生改变，直接影响地震资料的分辨率， 实际地层表现出一定的黏弹特性，因此通过对黏弹介质中的地震波进行数值模 拟，来研究和分析地震波传播过程中的衰减特征，对实际地震资料分辨率的提 高非常有意义。但是确保黏弹介质数值模拟稳定(即时间步长的确定)是一个急 需解决的问题，这关系到模拟算法的成败，也关系到模拟效率的高低。

目前主要运用黏弹固体模型来表征实际介质的粘滞特性，应用最广的黏弹 固体模型主要包括：Kelvin-Voigt固体模型、Maxwell固体模型和标准线性固 体模型(也称为标准线性黏弹固体模型)。其中，Kelvin-Voigt固体模型不能 考虑应力作用下应变的突然变化，也不能表示应力消失后的剩余应变，Maxwell 固体模型不具备蠕变特征，两者都不足以描述大多数黏弹介质的特征，但是标 准线性固体模型可以同时考虑具备应变突然变化和剩余应变及蠕动特征，因 此，标准线性固体模型能更加全面表征实际地层的黏弹特性，相比较而言更符 合实际情况，所以进行标准线性固体模型的黏弹介质数值模拟是非常有必要 的。

现有技术中标准线性固体模型的黏弹介质数值模拟主要是全三维、全波场 的模拟，此种方法的计算量非常大，从而产生的费用(主要体现在电力、存储设 备、人力等资源的消耗)是巨大的，数值模拟的效率低，因此有必要提供一种计 算方法定量化地给出标准线性固体模型数值模拟时间步长(即标准线性固体模型 的有限差分解的稳定性条件数值解)，从而更高效、更成功地指导数值模拟的进 行。

发明内容

本发明所要解决的技术问题是现有技术中采用全三维、全波场的模拟方法 产生的费用高、数值模拟效率低的缺陷。

为了解决上述技术问题，本发明提供了一种标准线性固体模型的稳定性条 件数值解的计算方法及系统。

本发明的技术方案为：

一种标准线性固体模型的稳定性条件数值解的计算方法，包括：

构建标准线性固体模型，并使所述标准线性固体模型包括彼此串联的第一 弹性体和第二弹性体、以及与所述第一弹性体并联的阻尼器；

确定所述标准线性固体模型有限差分解的稳定性条件的状态传递矩阵；

获取多组参数集，并使每组参数集中包括所述第一弹性体的弹性系数、所 述第二弹性体的弹性系数、所述阻尼器的黏滞系数、频率和介质密度；

在设定的空间差分精度和空间网格步长下，依次利用各组所述参数集并通 过计算机计算使得所述状态传递矩阵的特征值的模小于1的时间步长；

确定计算得出的时间步长为标准线性固体模型的稳定性条件数值解。

优选的是，所述确定所述标准线性固体模型有限差分解的稳定性条件的状态 传递矩阵包括：

确定所述标准线性固体模型的本构方程 $(M_1+M_3)p+M_2\frac{\partial p}{\partial t}=M_1{M}_3ϵ+M_2{M}_3\frac{\partial ϵ}{\partial t},$其中p为总应力，ε为总应变，弹性 系数为M₁的所述第一弹性体的应变与黏滞系数为M₂的所述阻尼器的应变相等， M₃为所述第二弹性体的弹性系数；

根据总应变ε与质点位移(u,v,w)间的关系方程以及所述 本构方程，得到第一方程：

$(M_1+M_3)p+M_2\frac{\partial p}{\partial t}=M_1{M}_3(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z})+M_2{M}_3\frac{\partial }{\partial t}(\frac{\partial u}{{\partial }_x}+\frac{\partial v}{{\partial }_y}+\frac{{\partial }_w}{{\partial }_z});$

对所述第一方程的左右两边分别对时间求二次偏导数，得到第二方程：

$\left(\begin{array}{c}(M_1+M_3)\frac{{\partial }^{2}p}{{\partial t}^2}+M_2\frac{{\partial }^{3}p}{{\partial t}^3}=M_1{M}_3(\frac{{\partial }^{3}u}{{\partial x\partial t}^2}+\frac{{\partial }^{3}v}{{\partial x\partial t}^2}+\frac{{\partial }^{3}w}{{\partial z\partial t}^2})\\ +M_2{M}_3\frac{\partial }{\partial t}(\frac{{\partial }^{3}u}{{\partial x\partial t}^2}+\frac{{\partial }^{3}v}{{\partial y\partial t}^2}+\frac{{\partial }^{3}w}{{\partial z\partial t}^2})\end{array}\right);$

利用所述第二方程和声波的纳维尔方程$\left(\begin{array}{c}\frac{{\partial }^{3}u}{{\partial x\partial t}^2}=\frac{1}{\rho }\frac{{\partial }^{2}P}{{\partial x}^2}\\ \frac{{\partial }^{3}v}{{\partial y\partial t}^2}=\frac{1}{\rho }\frac{{\partial }^{2}P}{{\partial y}^2}\\ \frac{{\partial }^{3}w}{{\partial z\partial t}^2}=\frac{1}{\rho }\frac{{\partial }^{2}P}{{\partial z}^2}\end{array}\right),$得到第三方 程：

$\left(\begin{array}{c}(M_1+M_3)\frac{{\partial }^{2}p}{{\partial t}^2}+M_2\frac{{\partial }^{3}p}{{\partial t}^3}=\frac{M_1{M}_3}{\rho }(\frac{{\partial }^{2}p}{{\partial x}^2}+\frac{{\partial }^{2}p}{{\partial y}^2}+\frac{{\partial }^{2}p}{{\partial z}^2})\\ +\frac{M_2{M}_3}{\rho }\frac{\partial }{\partial t}(\frac{{\partial }^{2}p}{{\partial x}^2}+\frac{{\partial }^{2}p}{{\partial y}^2}+\frac{{\partial }^{2}p}{{\partial z}^2})\end{array}\right),$其中ρ为所述介 质密度；

将所述第三方程中的应力取空间傅里叶变换，得到第四方程： $(M_1+M_3)\frac{{\partial }^2\hat{p}}{{\partial t}^2}+M_2\frac{{\partial }^3\hat{p}}{{\partial t}^3}=-\frac{M_1{M}_3}{\rho }{k}^2\hat{p}-\frac{M_2{M}_3}{\rho }{k}^2\frac{\partial \hat{p}}{\partial t},$其中为总应力p的空间傅 里叶变换，k为波数；

对所述第四方程中的时间偏导数用差分近似，得到第五方程：

$\left(\begin{array}{c}{\hat{p}}^{n+1}=\frac{[{3M}_2+2(M_1+M_3)\mathrm{\Delta t}-M_1{M}_3{\mathrm{\Delta t}}^3{k}^2/\rho -M_2{M}_3{\mathrm{\Delta t}}^2{k}^2/\rho ]}{[M_2+(M_1+M_3)\mathrm{\Delta t}]}{\hat{p}}^n\\ -\frac{[M_2{M}_3{\mathrm{\Delta t}}^2{k}^2/\rho -{3M}_2-(M_1+M_3)\mathrm{\Delta t}]}{[M_2+(M_1+M_3)\mathrm{\Delta t}]}{\hat{p}}^{n-1}+\frac{M_2}{[M_2+(M_1+M_3)\mathrm{\Delta t}]}{\hat{p}}^{n-2}\end{array}\right),$其中 分别为第n-2,n-1,n,n+1时刻的值，并且所述波数k满足： 在所述空间差分精度为2N的情况下，x,y,z三个 方向上的空间网格步长分别为Δx,Δy,Δz，a_l为对应所述空间差分精度2N的空间 差分系数，Δt为所述时间步长；

根据所述第五方程，得到所述标准线性固体模型有限差分解的稳定性条件的 状态传递矩阵$A=\left(\begin{array}{ccc}a& b& c\\ 1& 0& 0\\ 0& 1& 0\end{array}\right),$其中：

$a=\frac{[{3M}_2+2(M_1+M_3)\mathrm{\Delta t}-M_1{M}_3{\mathrm{\Delta t}}^3{k}^2/\rho -M_2{M}_3{\mathrm{\Delta t}}^2{k}^2/\rho ]}{[M_2+(M_1+M_3)\mathrm{\Delta t}]};$

$b=-\frac{[M_2{M}_3{\mathrm{\Delta t}}^2{k}^2/\rho -{3M}_2-(M_1+M_3)\mathrm{\Delta t}]}{[M_2+(M_1+M_3)\mathrm{\Delta t}]};$

$c=\frac{M_2}{[M_2+(M_1+M_3)\mathrm{\Delta t}]}.$

优选的是，所述方法还包括：在计算得出的所述时间步长后，根据所述参 数集计算所述标准线性固体模型的品质因子，并根据所述参数集和所述品质因 子计算所述标准线性固体模型的介质速度。

优选的是，在设定的空间差分精度和空间网格步长下，依次利用各组所述 参数集，运用安装在所述计算机上的Matlab仿真软件编程计算使得所述状态传 递矩阵的特征值的模小于1的时间步长。

一种标准线性固体模型的稳定性条件数值解的计算系统，包括：

模型构建单元，用于构建标准线性固体模型，并使所述标准线性固体模型 包括彼此串联的第一弹性体和第二弹性体、以及与所述第一弹性体并联的阻尼 器；

状态传递矩阵确定单元，用于确定所述标准线性固体模型有限差分解的稳 定性条件的状态传递矩阵；

参数集获取单元，用于获取多组参数集，并使每组参数集中包括所述第一 弹性体的弹性系数、所述第二弹性体的弹性系数、所述阻尼器的黏滞系数、频 率和介质密度；

时间步长计算单元，用于在设定的空间差分精度和空间网格步长下，依次 利用各组所述参数集并通过计算机计算使得所述状态传递矩阵的特征值的模小 于1的时间步长；

稳定性条件数值解确定单元，用于确定计算得出的时间步长为标准线性固 体模型的稳定性条件数值解。

优选的是，所述状态传递矩阵确定单元包括：

本构方程确定单元，用于确定所述标准线性固体模型的本构方程 $(M_1+M_3)p+M_2\frac{\partial p}{\partial t}=M_1{M}_3ϵ+M_2{M}_3\frac{\partial ϵ}{\partial t},$其中p为总应力，ε为总应变，弹性 系数为M₁的所述第一弹性体的应变与黏滞系数为M₂的所述阻尼器的应变相等， M₃为所述第二弹性体的弹性系数；

第一方程确定单元，用于根据总应变ε与质点位移(u,v,w)间的关系方程 以及所述本构方程确定单元确定的本构方程，得到第一方 程：

$(M_1+M_3)p+M_2\frac{\partial p}{\partial t}=M_1{M}_3(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z})+M_2{M}_3\frac{\partial }{\partial t}(\frac{\partial u}{{\partial }_x}+\frac{\partial v}{{\partial }_y}+\frac{{\partial }_w}{{\partial }_z});$

第二方程确定单元，用于对所述第一方程确定单元得到的第一方程的左右 两边分别对时间求二次偏导数，得到第二方程：

$\left(\begin{array}{c}(M_1+M_3)\frac{{\partial }^{2}p}{{\partial t}^2}+M_2\frac{{\partial }^{3}p}{{\partial t}^3}=M_1{M}_3(\frac{{\partial }^{3}u}{{\partial x\partial t}^2}+\frac{{\partial }^{3}v}{{\partial x\partial t}^2}+\frac{{\partial }^{3}w}{{\partial z\partial t}^2})\\ +M_2{M}_3\frac{\partial }{\partial t}(\frac{{\partial }^{3}u}{{\partial x\partial t}^2}+\frac{{\partial }^{3}v}{{\partial y\partial t}^2}+\frac{{\partial }^{3}w}{{\partial z\partial t}^2})\end{array}\right);$

第三方程确定单元，用于利用所述第二方程确定单元得到的第二方程和声波 的纳维尔方程$\left(\begin{array}{c}\frac{{\partial }^{3}u}{{\partial x\partial t}^2}=\frac{1}{\rho }\frac{{\partial }^{2}P}{{\partial x}^2}\\ \frac{{\partial }^{3}v}{{\partial y\partial t}^2}=\frac{1}{\rho }\frac{{\partial }^{2}P}{{\partial y}^2}\\ \frac{{\partial }^{3}w}{{\partial z\partial t}^2}=\frac{1}{\rho }\frac{{\partial }^{2}P}{{\partial z}^2}\end{array}\right),$得到第三方程：

$\left(\begin{array}{c}(M_1+M_3)\frac{{\partial }^{2}p}{{\partial t}^2}+M_2\frac{{\partial }^{3}p}{{\partial t}^3}=\frac{M_1{M}_3}{\rho }(\frac{{\partial }^{2}p}{{\partial x}^2}+\frac{{\partial }^{2}p}{{\partial y}^2}+\frac{{\partial }^{2}p}{{\partial z}^2})\\ +\frac{M_2{M}_3}{\rho }\frac{\partial }{\partial t}(\frac{{\partial }^{2}p}{{\partial x}^2}+\frac{{\partial }^{2}p}{{\partial y}^2}+\frac{{\partial }^{2}p}{{\partial z}^2})\end{array}\right),$其中ρ为所述介 质密度；

第四方程确定单元，用于将所述第三方程确定单元得到的第三方程中的应 力取空间傅里叶变换，得到第四方程： $(M_1+M_3)\frac{{\partial }^2\hat{p}}{{\partial t}^2}+M_2\frac{{\partial }^3\hat{p}}{{\partial t}^3}=-\frac{M_1{M}_3}{\rho }{k}^2\hat{p}-\frac{M_2{M}_3}{\rho }{k}^2\frac{\partial \hat{p}}{\partial t},$其中为总应力p的空间傅 里叶变换，k为波数；

第五方程确定单元，用于对所述第四方程确定单元得到的第四方程中的时 间偏导数用差分近似，得到第五方程：

$\left(\begin{array}{c}{\hat{p}}^{n+1}=\frac{[{3M}_2+2(M_1+M_3)\mathrm{\Delta t}-M_1{M}_3{\mathrm{\Delta t}}^3{k}^2/\rho -M_2{M}_3{\mathrm{\Delta t}}^2{k}^2/\rho ]}{[M_2+(M_1+M_3)\mathrm{\Delta t}]}{\hat{p}}^n\\ -\frac{[M_2{M}_3{\mathrm{\Delta t}}^2{k}^2/\rho -{3M}_2-(M_1+M_3)\mathrm{\Delta t}]}{[M_2+(M_1+M_3)\mathrm{\Delta t}]}{\hat{p}}^{n-1}+\frac{M_2}{[M_2+(M_1+M_3)\mathrm{\Delta t}]}{\hat{p}}^{n-2}\end{array}\right),$其中 分别为第n-2n-1nn+1时刻的值，并且所述波数k满足： 在所述空间差分精度为2N的情况下，x,y,z三个 方向上的空间网格步长分别为Δx,Δy,Δz，a_l为对应所述空间差分精度2N的空间 差分系数，Δt为所述时间步长；

状态传递矩阵子确定单元，用于根据所述第五方程确定单元得到的第五方 程，得到所述标准线性固体模型有限差分解的稳定性条件的状态传递矩阵 $A=\left(\begin{array}{ccc}a& b& c\\ 1& 0& 0\\ 0& 1& 0\end{array}\right),$其中：

$a=\frac{[{3M}_2+2(M_1+M_3)\mathrm{\Delta t}-M_1{M}_3{\mathrm{\Delta t}}^3{k}^2/\rho -M_2{M}_3{\mathrm{\Delta t}}^2{k}^2/\rho ]}{[M_2+(M_1+M_3)\mathrm{\Delta t}]};$

$b=-\frac{[M_2{M}_3{\mathrm{\Delta t}}^2{k}^2/\rho -{3M}_2-(M_1+M_3)\mathrm{\Delta t}]}{[M_2+(M_1+M_3)\mathrm{\Delta t}]};$

$c=\frac{M_2}{[M_2+(M_1+M_3)\mathrm{\Delta t}]}.$

优选的是，所述系统还包括：

品质因子计算单元，用于在所述时间步长计算单元计算得出所述时间步长 后，根据所述参数集获取单元获取的参数集计算所述标准线性固体模型的品质 因子；

介质速度计算单元，用于在所述时间步长计算单元计算得出所述时间步长 后，根据所述参数集获取单元获取的参数集和所述品质因子计算所述标准线性 固体模型的介质速度。

优选的是，所述时间步长计算单元，具体用于在设定的空间差分精度和空 间网格步长下，依次利用各组所述参数集，运用安装在所述计算机上的Matlab 仿真软件编程计算使得所述状态传递矩阵的特征值的模小于1的时间步长。

与现有技术相比，上述方案中的一个或多个实施例可以具有如下优点或有 益效果：

应用本发明实施例提供的标准线性固体模型的稳定性条件数值解的计算方 法及系统，从标准线性固体模型出发，进而得到标准线性固体模型有限差分解 的稳定性条件的状态传递矩阵，通过选取相应的时间步长使状态传递矩阵的特 征值的模小于1来满足黏弹数值模拟的稳定，从而能够定量化确定更符合实际 黏弹介质的标准线性固体模型数值模拟的时间步长，指导数值模拟顺利完成， 同时也使模拟的效率得以保障，最终运用模拟结果来解决实际地震勘探中的一 些地质问题。相比于现有技术中采用全三维、全波场的模拟方法进行标准线性 固体模型的黏弹介质数值模拟，本发明在确定时间步长后可确保基于标准线性 固体黏弹介质数值模拟的顺利完成，保证效率，节约成本，因此在实际工作 中，此项发明效果是非常显著的。

本发明的其它特征和优点将在随后的说明书中阐述，并且部分地从说明书 中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通 过在说明书、权利要求书以及附图中所特别指出的结构来实现和获得。

附图说明

附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发 明的实施例共同用于解释本发明，并不构成对本发明的限制。在附图中：

图1示出了本发明实施例标准线性固体模型的稳定性条件数值解的计算方 法的流程图；

图2示出了标准线性固体模型的结构示意图；

图3示出了本发明实施例中确定所述标准线性固体模型有限差分解的稳定 性条件的状态传递矩阵的方法的流程图；

图4示出了应用本发明实施例所述的计算方法，针对不同组合的参数集确 定的品质因子，时间步长与介质速度之间的关系曲线；

图5示出了应用本发明实施例所述的计算方法，针对不同组合的参数集确 定的介质速度，时间步长与品质因子之间的关系曲线；

图6示出了本发明实施例标准线性固体模型的稳定性条件数值解的计算系 统的结构示意图；

图7示出了本发明实施例中状态传递矩阵确定单元的结构示意图。

具体实施方式

以下将结合附图及实施例来详细说明本发明的实施方式，借此对本发明如 何应用技术手段来解决技术问题，并达成技术效果的实现过程能充分理解并据 以实施。需要说明的是，只要不构成冲突，本发明中的各个实施例以及各实施 例中的各个特征可以相互结合，所形成的技术方案均在本发明的保护范围之 内。

为解决现有技术中采用全三维、全波场的模拟方法进行标准线性固体模型 的黏弹介质数值模拟，产生的费用高，模拟效率低的缺陷，本发明实施例提供了 一种标准线性固体模型的稳定性条件数值解的计算方法，该方法能够定量化确 定更符合实际黏弹介质的标准线性固体模型数值模拟的时间步长，指导数值模 拟顺利完成，同时也使模拟的效率得以保障，节约了成本。

如图1所示，是本发明实施例标准线性固体模型的稳定性条件数值解的计 算方法的流程图，所述标准线性固体模型的稳定性条件数值解的计算方法，包 括以下步骤：

步骤101：构建标准线性固体模型，并使所述标准线性固体模型包括彼此串 联的第一弹性体和第二弹性体、以及与所述第一弹性体并联的阻尼器。

具体地，标准线性固体模型的结构示意图可参照图2，该标准线性固体模型 包括第一弹性体1、第二弹性体3和阻尼器2，其中第一弹性体1与第二弹性体 3串联连接，并且第一弹性体1和阻尼器2并联连接，并且第一弹性体1的弹性 系数表示为M1,第二弹性体3的弹性系数表示为M₃，阻尼器2的黏滞系数表示 为M₂，第一弹性体1和阻尼器2的应变相等，均表示为ε1，第二弹性体3的应 变表示为ε2。

步骤102：确定所述标准线性固体模型有限差分解的稳定性条件的状态传递 矩阵。

具体地，标准线性固体模型有限差分解的稳定性条件的状态传递矩阵的确 定方法，将在本发明的下一实施例中进行详细地说明。

步骤103：获取多组参数集，并使每组参数集中包括所述第一弹性体1的弹 性系数、所述第二弹性体3的弹性系数、所述阻尼器2的黏滞系数、频率和介质 密度。

具体地，在本步骤中，准备多组参数集，以备在步骤104中进行时间步长 的计算。每组参数集中包括5个参数，分别为第一弹性体1的弹性系数M1、所 述第二弹性体3的弹性系数M₃、所述阻尼器2的黏滞系数M₂、频率ω(这里 ω＝2πf)和介质密度ρ。

步骤104：在设定的空间差分精度和空间网格步长下，依次利用各组所述参 数集并通过计算机计算使得所述状态传递矩阵的特征值的模小于1的时间步 长。

具体地，所述设定的空间差分精度表示为2N，在所述空间差分精度为2N 的情况下，x,y,z三个方向上的空间网格步长分别表示为Δx,Δy,Δz。在设定所 述空间差分精度2N和空间网格步长下，依次将各组参数集输入至计算机内， 通过计算机计算使得所述状态传递矩阵的特征值的模小于1的时间步长。这 里，需要使用计算机计算时间步长的原因是，通过步骤102得到的标准线性固 体模型有限差分解的稳定性条件的状态传递矩阵特别复杂，无法直接求取该状 态传递矩阵的特征值表达式，从而也就无法计算使得该状态传递矩阵的模小于1 的时间步长。

步骤105：确定计算得出的时间步长为标准线性固体模型的稳定性条件数值 解。

在本实施例中，应用本发明实施例提供的标准线性固体模型的稳定性条件 数值解的计算方法，从标准线性固体模型出发，进而得到标准线性固体模型有 限差分解的稳定性条件的状态传递矩阵，通过选取相应的时间步长使状态传递 矩阵的特征值的模小于1来满足黏弹数值模拟的稳定，从而能够定量化确定更 符合实际黏弹介质的标准线性固体模型数值模拟的时间步长，指导数值模拟顺 利完成，同时也使模拟的效率得以保障，最终运用模拟结果来解决实际地震勘 探中的一些地质问题。相比于现有技术中采用全三维、全波场的模拟方法进行 标准线性固体模型的黏弹介质数值模拟，本方法在确定时间步长后可确保基于 标准线性固体黏弹介质数值模拟的顺利完成，保证效率，节约成本，因此在实 际工作中，此项发明效果是非常显著的。

进一步地，如图3所示，是本发明实施例中确定所述标准线性固体模型有限 差分解的稳定性条件的状态传递矩阵的方法的流程图，该方法包括以下步骤：

步骤201：首先得出控制标准线性固体模型的方程为：

$p=M_1{ϵ}_1+M_2\frac{{\partial ϵ}_1}{\partial t}=M_3{ϵ}_2---\left(1\right)$

ε＝ε₁+ε₂(2)

又因为ε₂＝p/M₃，将该公式代入(2)式，再联合(1)式，消去ε₁和ε₂后， 得到应力与应变的关系(即标准线性固体模型的本构方程)为：

$(M_1+M_3)p+M_2\frac{\partial p}{\partial t}=M_1{M}_3ϵ+M_2{M}_3\frac{\partial ϵ}{\partial t}---\left(3\right)$

其中p为总应力，ε为总应变，ε₁为弹性系数为M₁的第一弹性体1的应变、 以及黏滞系数为M₂的阻尼器2的应变，ε₂为弹性系数为M₃的第二弹性体3的 应变。

步骤202：总应变ε与质点位移(u,v,w)间的关系方程为：

$ϵ=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}---\left(4\right)$

将(4)式代入上述本构方程(即式(3))，得到第一方程(即式(5))：

$(M_1+M_3)p+M_2\frac{\partial p}{\partial t}=M_1{M}_3(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z})+M_2{M}_3\frac{\partial }{\partial t}(\frac{\partial u}{{\partial }_x}+\frac{\partial v}{{\partial }_y}+\frac{{\partial }_w}{{\partial }_z})---\left(5\right)$

步骤203：对所述第一方程(即式(5))的左右两边分别对时间t求二次偏 导数，得到第二方程(即式(6))：

$\left(\begin{array}{c}(M_1+M_3)\frac{{\partial }^{2}p}{{\partial t}^2}+M_2\frac{{\partial }^{3}p}{{\partial t}^3}=M_1{M}_3(\frac{{\partial }^{3}u}{{\partial x\partial t}^2}+\frac{{\partial }^{3}v}{{\partial x\partial t}^2}+\frac{{\partial }^{3}w}{{\partial z\partial t}^2})\\ +M_2{M}_3\frac{\partial }{\partial t}(\frac{{\partial }^{3}u}{{\partial x\partial t}^2}+\frac{{\partial }^{3}v}{{\partial y\partial t}^2}+\frac{{\partial }^{3}w}{{\partial z\partial t}^2})\end{array}\right)---\left(6\right)$

步骤204：利用所述第二方程(即式(6))和声波的纳维尔方程：

$\left(\begin{array}{c}\frac{{\partial }^{3}u}{{\partial x\partial t}^2}=\frac{1}{\rho }\frac{{\partial }^{2}P}{{\partial x}^2}\\ \frac{{\partial }^{3}v}{{\partial y\partial t}^2}=\frac{1}{\rho }\frac{{\partial }^{2}P}{{\partial y}^2}\\ \frac{{\partial }^{3}w}{{\partial z\partial t}^2}=\frac{1}{\rho }\frac{{\partial }^{2}P}{{\partial z}^2}\end{array}\right)---\left(7\right)$

得到第三方程(即式(8))：

$\left(\begin{array}{c}(M_1+M_3)\frac{{\partial }^{2}p}{{\partial t}^2}+M_2\frac{{\partial }^{3}p}{{\partial t}^3}=\frac{M_1{M}_3}{\rho }(\frac{{\partial }^{2}p}{{\partial x}^2}+\frac{{\partial }^{2}p}{{\partial y}^2}+\frac{{\partial }^{2}p}{{\partial z}^2})\\ +\frac{M_2{M}_3}{\rho }\frac{\partial }{\partial t}(\frac{{\partial }^{2}p}{{\partial x}^2}+\frac{{\partial }^{2}p}{{\partial y}^2}+\frac{{\partial }^{2}p}{{\partial z}^2})\end{array}\right)---\left(8\right)$

其中，ρ为介质密度。

步骤205：将所述第三方程(即式(8))中的应力取空间傅里叶变换，得 到第四方程(即式(9))：

$(M_1+M_3)\frac{{\partial }^2\hat{p}}{{\partial t}^2}+M_2\frac{{\partial }^3\hat{p}}{{\partial t}^3}=-\frac{M_1{M}_3}{\rho }{k}^2\hat{p}-\frac{M_2{M}_3}{\rho }{k}^2\frac{\partial \hat{p}}{\partial t}---\left(9\right)$

其中，为总应力p的空间傅里叶变换，k为波数。

步骤206：对所述第四方程(即式(9))中的时间偏导数用差分近似，得 到第五方程(即式(10))：

$\left(\begin{array}{c}{\hat{p}}^{n+1}=\frac{[{3M}_2+2(M_1+M_3)\mathrm{\Delta t}-M_1{M}_3{\mathrm{\Delta t}}^3{k}^2/\rho -M_2{M}_3{\mathrm{\Delta t}}^2{k}^2/\rho ]}{[M_2+(M_1+M_3)\mathrm{\Delta t}]}{\hat{p}}^n\\ -\frac{[M_2{M}_3{\mathrm{\Delta t}}^2{k}^2/\rho -{3M}_2-(M_1+M_3)\mathrm{\Delta t}]}{[M_2+(M_1+M_3)\mathrm{\Delta t}]}{\hat{p}}^{n-1}+\frac{M_2}{[M_2+(M_1+M_3)\mathrm{\Delta t}]}{\hat{p}}^{n-2}\end{array}\right)---\left(10\right)$

其中分别为第n-2,n-1,n,n+1时刻的值，并且所述波数 k满足：在所述空间差分精度为2N的情况下，

x,y,z三个方向上的空间网格步长分别为Δx,Δy,Δz，a_l为对应所述空间差分精度 2N的空间差分系数，Δt为所述时间步长。

步骤207：根据所述第五方程(即式(10))，得到所述标准线性固体模型 有限差分解的稳定性条件的状态传递矩阵$A=\left(\begin{array}{ccc}a& b& c\\ 1& 0& 0\\ 0& 1& 0\end{array}\right),$其中：

$a=\frac{[{3M}_2+2(M_1+M_3)\mathrm{\Delta t}-M_1{M}_3{\mathrm{\Delta t}}^3{k}^2/\rho -M_2{M}_3{\mathrm{\Delta t}}^2{k}^2/\rho ]}{[M_2+(M_1+M_3)\mathrm{\Delta t}]};$

$b=-\frac{[M_2{M}_3{\mathrm{\Delta t}}^2{k}^2/\rho -{3M}_2-(M_1+M_3)\mathrm{\Delta t}]}{[M_2+(M_1+M_3)\mathrm{\Delta t}]};$

$c=\frac{M_2}{[M_2+(M_1+M_3)\mathrm{\Delta t}]}.$

在本实施例中，通过上述步骤可以确定标准线性固体模型的有限差分解的 稳定性条件的状态传递矩阵，另外根据本实施例中涉及的所有的公式可以看 出，该状态传递矩阵实际上是关于以下9个参数的表达式：M₁(第一弹性体1 的弹性系数)、M₂(阻尼器2的黏滞系数)、M₃(第二弹性体3的弹性系数)、 ω(频率)、ρ(介质密度)、Δx,Δy,Δz(空间网格步长)、2N(空间差分精 度)、al(空间差分系数)和Δt(时间步长)，也就是说，可以在设定上述前8 个参数的基础上，借助计算机计算得出满足下列要求的时间步长Δt：得到的时 间步长需要使得状态传递矩阵的特征值的模小于1。

进一步地，所述方法还包括以下步骤：在计算得出的所述时间步长后，根 据所述参数集计算所述标准线性固体模型的品质因子，并根据所述参数集和所 述品质因子计算所述标准线性固体模型的介质速度。

下面结合具体实例、以及图4和图5说明本实施例，设定空间网格步长 Δx＝Δy＝Δz＝5m，经查二阶导数2N各阶精度的权系数值(据孙成禹，2007)，空 间差分精度的空间差分系数依次分别为：

a₀＝-2.7222222,a_-1＝a₁＝1.5000000,a_-2＝a₂＝-0.15000000,a_-3＝a₃＝0.011111111。

标准线性固体模型的品质因子的计算公式：

$Q=\frac{{M_1}^2}{M_2{M}_3\omega }+\frac{M_1}{M_2\omega }+\frac{M_2\omega }{M_3}$

标准线性固体模型的介质速度的计算公式：

$v^{\prime 2}=\frac{2(Q^2+1)[(M_1+M_3)M_1{M}_3+{M_2}^2{\omega }^2{M}_3]}{[{(Q^2+1)}^{1/2}+Q]Q[{(M_1+M_3)}^2+{\left(M_2\omega \right)}^2]\rho },$当Q→∞时，变为弹性介 质情形：上述介质速度的计算公式变为：

$v^{\prime 2}=v^2=\frac{(M_1+M_3)M_1{M}_3+{M_2}^2{\omega }^2{M}_3}{[{(M_1+M_3)}^2+{\left(M_2\omega \right)}^2]\rho }.$

通过设置不同参数集(M₁,M₂,M₃,ω,ρ)，使得品质因子Q的取值分别为1、 10、100、1000的情况下，时间步长Δt与介质速度v的定量关系(如图4所示)； 或者通过设置不同参数集(M₁,M₂,M₃,ω,ρ)，使得介质速度v的取值分别为 500m/s、1000m/s、1500m/s、2000m/s、3000m/s、4000m/s和5000m/s的 情况下，时间步长Δt与品质因子Q的定量关系(如图5所示)。由图4和图5 可知，在设置特定参数集的情况下，可以找到相应的时间步长。另外，从图4 和图5中还可以看出，随着介质速度的增大，稳定性条件变差，要求的时间步 长变小。

进一步地，在上述步骤104中，运用安装在所述计算机上的Matlab仿真软 件编程计算所述时间步长，值得说明的是，利用Matlab求解使得一个复杂矩阵 满足特定需求时该矩阵中某个参量的取值，是本领域技术人员常规采用的技术 手段，故在本文中不再进行展开说明。

相应地，本发明实施例还提供一种标准线性固体模型的稳定性条件数值解的 计算系统，图6示出了本发明实施例标准线性固体模型的稳定性条件数值解的计 算系统的结构示意图，如图6所示，该系统包括：

模型构建单元301，用于构建标准线性固体模型，并使所述标准线性固体模 型包括彼此串联的第一弹性体1和第二弹性体3、以及与所述第一弹性体1并联 的阻尼器2；

状态传递矩阵确定单元302，用于确定所述模型构建单元301构建的标准线 性固体模型的有限差分解的稳定性条件的状态传递矩阵；

参数集获取单元303，用于获取多组参数集，并使每组参数集中包括所述第 一弹性体1的弹性系数、所述第二弹性体3的弹性系数、所述阻尼器2的黏滞系 数、频率和介质密度；

时间步长计算单元304，用于在设定的空间差分精度和空间网格步长下，依 次利用所述参数集获取单元303获取的各组参数集并通过计算机计算使得所述 状态传递矩阵确定单元302确定的状态传递矩阵的特征值的模小于1的时间步 长；

稳定性条件数值解确定单元305，用于确定所述时间步长计算单元304计算 得出的时间步长为标准线性固体模型的稳定性条件数值解。

进一步地，如图7所示，是本发明实施例中状态传递矩阵确定单元302的结 构示意图，所述状态传递矩阵确定单元302包括：

本构方程确定单元401，用于确定所述标准线性固体模型的本构方程 $(M_1+M_3)p+M_2\frac{\partial p}{\partial t}=M_1{M}_3ϵ+M_2{M}_3\frac{\partial ϵ}{\partial t},$其中p为总应力，ε为总应变，弹性 系数为M₁的所述第一弹性体1的应变与黏滞系数为M₂的所述阻尼器2的应变相 等，M₃为所述第二弹性体3的弹性系数；

第一方程确定单元402，用于根据总应变ε与质点位移(u,v,w)间的关系方程 以及所述本构方程确定单元401确定的本构方程，得到第一 方程：

$(M_1+M_3)p+M_2\frac{\partial p}{\partial t}=M_1{M}_3(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z})+M_2{M}_3\frac{\partial }{\partial t}(\frac{\partial u}{{\partial }_x}+\frac{\partial v}{{\partial }_y}+\frac{{\partial }_w}{{\partial }_z});$

第二方程确定单元403，用于对所述第一方程确定单元402得到的第一方程 的左右两边分别对时间求二次偏导数，得到第二方程：

$\left(\begin{array}{c}(M_1+M_3)\frac{{\partial }^{2}p}{{\partial t}^2}+M_2\frac{{\partial }^{3}p}{{\partial t}^3}=M_1{M}_3(\frac{{\partial }^{3}u}{{\partial x\partial t}^2}+\frac{{\partial }^{3}v}{{\partial x\partial t}^2}+\frac{{\partial }^{3}w}{{\partial z\partial t}^2})\\ +M_2{M}_3\frac{\partial }{\partial t}(\frac{{\partial }^{3}u}{{\partial x\partial t}^2}+\frac{{\partial }^{3}v}{{\partial y\partial t}^2}+\frac{{\partial }^{3}w}{{\partial z\partial t}^2})\end{array}\right);$

第三方程确定单元404，用于利用所述第二方程确定单元403得到的第二方 程和声波的纳维尔方程$\left(\begin{array}{c}\frac{{\partial }^{3}u}{{\partial x\partial t}^2}=\frac{1}{\rho }\frac{{\partial }^{2}P}{{\partial x}^2}\\ \frac{{\partial }^{3}v}{{\partial y\partial t}^2}=\frac{1}{\rho }\frac{{\partial }^{2}P}{{\partial y}^2}\\ \frac{{\partial }^{3}w}{{\partial z\partial t}^2}=\frac{1}{\rho }\frac{{\partial }^{2}P}{{\partial z}^2}\end{array}\right),$得到第三方程：

$\left(\begin{array}{c}(M_1+M_3)\frac{{\partial }^{2}p}{{\partial t}^2}+M_2\frac{{\partial }^{3}p}{{\partial t}^3}=\frac{M_1{M}_3}{\rho }(\frac{{\partial }^{2}p}{{\partial x}^2}+\frac{{\partial }^{2}p}{{\partial y}^2}+\frac{{\partial }^{2}p}{{\partial z}^2})\\ +\frac{M_2{M}_3}{\rho }\frac{\partial }{\partial t}(\frac{{\partial }^{2}p}{{\partial x}^2}+\frac{{\partial }^{2}p}{{\partial y}^2}+\frac{{\partial }^{2}p}{{\partial z}^2})\end{array}\right),$其中ρ为所述介 质密度；

第四方程确定单元405，用于将所述第三方程确定单元404得到的第三方程 中的应力取空间傅里叶变换，得到第四方程： $\left(\begin{array}{c}(M_1+M_3)\frac{{\partial }^{2}p}{{\partial t}^2}+M_2\frac{{\partial }^{3}p}{{\partial t}^3}=\frac{M_1{M}_3}{\rho }(\frac{{\partial }^{2}p}{{\partial x}^2}+\frac{{\partial }^{2}p}{{\partial y}^2}+\frac{{\partial }^{2}p}{{\partial z}^2})\\ +\frac{M_2{M}_3}{\rho }\frac{\partial }{\partial t}(\frac{{\partial }^{2}p}{{\partial x}^2}+\frac{{\partial }^{2}p}{{\partial y}^2}+\frac{{\partial }^{2}p}{{\partial z}^2})\end{array}\right),$其中为总应力p的空间傅 里叶变换，k为波数；

第五方程确定单元406，用于对所述第四方程确定单元405得到的第四方程 中的时间偏导数用差分近似，得到第五方程：

$\left(\begin{array}{c}{\hat{p}}^{n+1}=\frac{[{3M}_2+2(M_1+M_3)\mathrm{\Delta t}-M_1{M}_3{\mathrm{\Delta t}}^3{k}^2/\rho -M_2{M}_3{\mathrm{\Delta t}}^2{k}^2/\rho ]}{[M_2+(M_1+M_3)\mathrm{\Delta t}]}{\hat{p}}^n\\ -\frac{[M_2{M}_3{\mathrm{\Delta t}}^2{k}^2/\rho -{3M}_2-(M_1+M_3)\mathrm{\Delta t}]}{[M_2+(M_1+M_3)\mathrm{\Delta t}]}{\hat{p}}^{n-1}+\frac{M_2}{[M_2+(M_1+M_3)\mathrm{\Delta t}]}{\hat{p}}^{n-2}\end{array}\right),$其中 分别为第n-2,n-1,n,n+1时刻的值，并且所述波数k满足： 在所述空间差分精度为2N的情况下，x,y,z三个 方向上的空间网格步长分别为Δx,Δy,Δz，a_l为对应所述空间差分精度2N的空间 差分系数，Δt为所述时间步长；

状态传递矩阵子确定单元407，用于根据所述第五方程确定单元406得到的 第五方程，得到所述标准线性固体模型有限差分解的稳定性条件的状态传递矩阵 $A=\left(\begin{array}{ccc}a& b& c\\ 1& 0& 0\\ 0& 1& 0\end{array}\right),$其中：

$a=\frac{[{3M}_2+2(M_1+M_3)\mathrm{\Delta t}-M_1{M}_3{\mathrm{\Delta t}}^3{k}^2/\rho -M_2{M}_3{\mathrm{\Delta t}}^2{k}^2/\rho ]}{[M_2+(M_1+M_3)\mathrm{\Delta t}]};$

$b=-\frac{[M_2{M}_3{\mathrm{\Delta t}}^2{k}^2/\rho -{3M}_2-(M_1+M_3)\mathrm{\Delta t}]}{[M_2+(M_1+M_3)\mathrm{\Delta t}]};$

$c=\frac{M_2}{[M_2+(M_1+M_3)\mathrm{\Delta t}]}.$

进一步地，所述系统还包括品质因子计算单元和介质速度计算单元，其中品 质因子计算单元用于在所述时间步长计算单元304计算得出所述时间步长后，根 据所述参数集获取单元303获取的参数集计算所述标准线性固体模型的品质因 子；所述介质速度计算单元，用于在所述时间步长计算单元304计算得出所述时 间步长后，根据所述参数集获取单元303获取的参数集和所述品质因子计算所述 标准线性固体模型的介质速度。

特别地，在本发明的一优选的实施例中，所述时间步长计算单元304，具体 用于在设定的空间差分精度和空间网格步长下，依次利用各组所述参数集，运用 安装在所述计算机上的Matlab仿真软件编程计算使得所述状态传递矩阵的特征 值的模小于1的时间步长。

上述各单元的具体处理过程可参照前面本发明实施例的方法中的描述，在此 不再赘述。

在本实施例中，应用本发明实施例提供的标准线性固体模型的稳定性条件 数值解的计算系统，从标准线性固体模型出发，进而得到标准线性固体模型有 限差分解的稳定性条件的状态传递矩阵，通过选取相应的时间步长使状态传递 矩阵的特征值的模小于1来满足黏弹数值模拟的稳定，从而能够定量化确定更 符合实际黏弹介质的标准线性固体模型数值模拟的时间步长，指导数值模拟顺 利完成，同时也使模拟的效率得以保障，最终运用模拟结果来解决实际地震勘 探中的一些地质问题。相比于现有技术中采用全三维、全波场的模拟方法进行 标准线性固体模型的黏弹介质数值模拟，本系统在确定时间步长后可确保基于 标准线性固体黏弹介质数值模拟的顺利完成，保证效率，节约成本，因此在实 际工作中，此项发明效果是非常显著的。

本领域的技术人员应该明白，上述的本发明的各模块或各步骤可以用通用 的计算装置来实现，它们可以集中在单个的计算装置上，或者分布在多个计算 装置所组成的网络上，可选地，它们可以用计算装置可执行的程序代码来实 现，从而，可以将它们存储在存储装置中由计算装置来执行，或者将它们分别 制作成各个集成电路模块，或者将它们中的多个模块或步骤制作成单个集成电 路模块来实现。这样，本发明不限制于任何特定的硬件和软件结合。

虽然本发明所公开的实施方式如上，但所述的内容只是为了便于理解本发 明而采用的实施方式，并非用以限定本发明。任何本发明所属技术领域内的技 术人员，在不脱离本发明所公开的精神和范围的前提下，可以在实施的形式上 及细节上作任何的修改与变化，但本发明的专利保护范围，仍须以所附的权利 要求书所界定的范围为准。