标准线性固体模型的稳定性条件数值解的计算方法及系统的制作方法
【技术领域】
[0001] 本发明涉及地震勘探基础应用技术领域，尤其涉及一种标准线性固体模型的有限 差分解的稳定性条件数值解的计算方法及系统。
【背景技术】
[0002] 地震数值模拟是地震勘探和地震学的重要基础，同时也是了解复杂介质中地震波 传播规律的重要工具，其作用贯穿于整个地震采集、处理和解释中。随着地震勘探开发的深 入，常规的弹性介质理论难W满足实际介质需求。地震波在实际地层中传播时，能量和相位 都发生改变，直接影响地震资料的分辨率，实际地层表现出一定的黏弹特性，因此通过对黏 弹介质中的地震波进行数值模拟，来研究和分析地震波传播过程中的衰减特征，对实际地 震资料分辨率的提高非常有意义。但是确保黏弹介质数值模拟稳定（即时间步长的确定） 是一个急需解决的问题，送关系到模拟算法的成败，也关系到模拟效率的高低。
[0003] 目前主要运用黏弹固体模型来表征实际介质的粘滞特性，应用最广的黏弹固体模 型主要包括；Kelvin-Voigt固体模型、Maxwell固体模型和标准线性固体模型（也称为标准 线性黏弹固体模型）。其中，Kelvin-Voigt固体模型不能考虑应力作用下应变的突然变化， 也不能表示应力消失后的剩余应变，Maxwell固体模型不具备蠕变特征，两者都不足W描述 大多数黏弹介质的特征，但是标准线性固体模型可W同时考虑具备应变突然变化和剩余应 变及蠕动特征，因此，标准线性固体模型能更加全面表征实际地层的黏弹特性，相比较而言 更符合实际情况，所W进行标准线性固体模型的黏弹介质数值模拟是非常有必要的。
[0004] 现有技术中标准线性固体模型的黏弹介质数值模拟主要是全Η维、全波场的模 拟，此种方法的计算量非常大，从而产生的费用（主要体现在电力、存储设备、人力等资源 的消耗）是巨大的，数值模拟的效率低，因此有必要提供一种计算方法定量化地给出标准 线性固体模型数值模拟时间步长（即标准线性固体模型的有限差分解的稳定性条件数值 解），从而更高效、更成功地指导数值模拟的进行。
【发明内容】
[0005] 本发明所要解决的技术问题是现有技术中采用全Η维、全波场的模拟方法产生的 费用高、数值模拟效率低的缺陷。
[0006] 为了解决上述技术问题，本发明提供了一种标准线性固体模型的稳定性条件数值 解的计算方法及系统。
[0007] 本发明的技术方案为：
[0008] -种标准线性固体模型的稳定性条件数值解的计算方法，包括：
[0009] 构建标准线性固体模型，并使所述标准线性固体模型包括彼此串联的第一弹性体 和第二弹性体、W及与所述第一弹性体并联的阻尼器；
[0010] 确定所述标准线性固体模型有限差分解的稳定性条件的状态传递矩阵；
[0011] 获取多组参数集，并使每组参数集中包括所述第一弹性体的弹性系数、所述第二 弹性体的弹性系数、所述阻尼器的黏滞系数、频率和介质密度；
[0012] 在设定的空间差分精度和空间网格步长下，依次利用各组所述参数集并通过计算 机计算使得所述状态传递矩阵的特征值的模小于1的时间步长；
[0013] 确定计算得出的时间步长为标准线性固体模型的稳定性条件数值解。
[0014] 优选的是，所述确定所述标准线性固体模型有限差分解的稳定性条件的状态传递 矩阵包括：
[0015] 确定所述标准线性固体模型的本构方程
其中P为总应力，ε为总应变，弹性系数为Ml的 所述第一弹性体的应变与黏滞系数为M2的所述阻尼器的应变相等，Ms为所述第二弹性体的 弹性系数；
[0016] 根据总应变ε与质点位移（u，v，w)间的关系方程
W及所述本 构方程，得到第一方程：
[0017]
[0018] 对所述第一方程的左右两边分别对时间求二次偏导数，得到第二方程：
[0019]
[0020] 利用所述第二方程和声波的纳维尔方程
得到第Η方程：
?
[0021] 其中Ρ为所述介质密 度；
[0022] 将所述第S方程中的应力取空间傅里叶变换，得到第四方程：
其中#为总应力Ρ的空间傅里叶变换， k为波数；
[0023] 对所述第四方程中的时间偏导数用差分近似，得到第五方程：
[0024]
其 中护声"、；ri分别为第n-2, n-1，η, n+1时刻的梦值，并且所述波数k满足：
在所述空间差分精度为2N的情况下，X，y，Z Η个方向上的空 间网格步长分别为Δχ，Ay, Δζ，曰1为对应所述空间差分精度2Ν的空间差分系数，At为 所述时间步长；
[0025] 根据所述第五方程，得到所述标准线性固体模型有限差分解的稳定性条件的状态 传递矩阵
其中：
[0029] 优选的是，所述方法还包括；在计算得出的所述时间步长后，根据所述参数集计算 所述标准线性固体模型的品质因子，并根据所述参数集和所述品质因子计算所述标准线性 固体模型的介质速度。
[0030] 优选的是，在设定的空间差分精度和空间网格步长下，依次利用各组所述参数集， 运用安装在所述计算机上的Matlab仿真软件编程计算使得所述状态传递矩阵的特征值的 模小于1的时间步长。
[0031] 一种标准线性固体模型的稳定性条件数值解的计算系统，包括：
[0032] 模型构建单元，用于构建标准线性固体模型，并使所述标准线性固体模型包括彼 此串联的第一弹性体和第二弹性体、W及与所述第一弹性体并联的阻尼器；
[0033] 状态传递矩阵确定单元，用于确定所述标准线性固体模型有限差分解的稳定性条 件的状态传递矩阵；
[0034] 参数集获取单元，用于获取多组参数集，并使每组参数集中包括所述第一弹性体 的弹性系数、所述第二弹性体的弹性系数、所述阻尼器的黏滞系数、频率和介质密度；
[0035] 时间步长计算单元，用于在设定的空间差分精度和空间网格步长下，依次利用各 组所述参数集并通过计算机计算使得所述状态传递矩阵的特征值的模小于1的时间步长；
[0036] 稳定性条件数值解确定单元，用于确定计算得出的时间步长为标准线性固体模型 的稳定性条件数值解。
[0037] 优选的是，所述状态传递矩阵确定单元包括：
[003引本构方程确定单元，用于确定所述标准线性固体模型的本构方程
其中P为总应力，ε为总应变，弹性系数为Ml的 所述第一弹性体的应变与黏滞系数为M2的所述阻尼器的应变相等，Ms为所述第二弹性体的 弹性系数；
[0039] 第一方程确定单元，用于根据总应变ε与质点位移（u，v，w)间的关系方程
W及所述本构方程确定单元确定的本构方程，得到第一方程：
[0040]
[0041] 第二方程确定单元，用于对所述第一方程确定单元得到的第一方程的左右两边分 别对时间求二次偏导数，得到第二方程：
[0042]
[0043] 第Η方程确定单元，用于利用所述第二方程确定单元得到的第二方程和声波的纳 维尔方程
得到第Η方程：
>
[0044] 其中Ρ为所述介质密 度；
[0045] 第四方程确定单元，用于将所述第Η方程确定单元得到的第Η方程中的应力取空 间傅里叶变换，得到第四方程：
其中# 为总应力Ρ的空间傅里叶变换，k为波数；
[0046] 第五方程确定单元，用于对所述第四方程确定单元得到的第四方程中的时间偏导 数用差分近似，得到第五方程：
[0047]
其中r-气护1.护扩'' 分别为第η-化-lnn+1时刻的#值，并且所述波数k满足：
在所述空间差分精度为2N的情况下，X，y，Z Η个方向上的空 间网格步长分别为Δχ，Ay, Δζ，曰1为对应所述空间差分精度2Ν的空间差分系数，At为 所述时间步长；
[0048] 状态传递矩阵子确定单元，用于根据所述第五方程确定单元得到的第五方程，得 到所述标准线性固体模型有限差分解的稳定性条件的状态传递矩阵
痒中：
[0052] 优选的是，所述系统还包括：
[0053] 品质因子计算单元，用于在所述时