一种粘滞各向异性介质地震波数值模拟方法、装置及设备

1.本文涉及地球物理技术领域，具体涉及一种粘滞各向异性介质地震波数值模拟方法、装置及设备。


背景技术：

2.在油气探测开发中，地震勘探技术发挥着极其重要的作用。经过数十年的开采，来自简单圈闭构造的油气资源几近枯竭，工业界的勘探重心也逐渐由浅层转向深层、由常规地带转向极端地带、由构造油气藏转向隐蔽的岩性油气藏。因此，发展针对复杂介质的地震波场高精度数值模拟技术已经成为当今工业界的共识。其中，地层对地震波能量的吸收衰减和弹性参数随方位的变化(即各向异性)是高精度波场模拟不可忽略的重要因素。
3.已有研究表明，地下介质的很多性质都呈现出与方向或方位的相关性，即通常说的各向异性。现有技术的研究主要集中在速度各向异性上，却忽略了吸收衰减各向异性对地震波传播的影响，现有技术中提出了一个各向异性时间分数阶粘滞弹性波方程，在该方程的刚度矩阵中，分别含有速度各向异性和衰减各向异性参数，所使用的速度各向异性参数是thomsen系数，衰减各向异性参数是由一个形式上类似的q矩阵来表征的，但是由于其中包含的时间分数阶算子在求解过程中需要存储大量的当前时刻以前的波场值，因此导致该方程需要较大的运算能力才能计算出准确的地震波变化，从而提高了数据计算的成本，进而导致上述方程实用性较差，很难在实际生产开发中运行。因此亟需一种能能够同时准确描述速度和衰减各向异性、同时又便于数值计算方程，提高地震波数值模拟的效率。


技术实现要素：

4.针对现有技术的上述问题，本文的目的在于，提供一种粘滞各向异性介质地震波数值模拟方法、装置及设备，能够提高地震波数值模拟的准确性。
5.为了解决上述技术问题，本文的具体技术方案如下：
6.一方面，本文提供一种粘滞各向异性介质地震波数值模拟方法，所述方法包括：
7.根据已知弹性介质的第一刚度矩阵，确定粘弹性介质的第一品质因子矩阵；
8.将所述第一刚度矩阵和所述第一品质因子矩阵进行简化处理，得到第二刚度矩阵和第二品质因子矩阵；
9.根据所述二刚度矩阵和所述第二品质因子矩阵，结合预设频散关系、运动平衡方程和几何方程，计算获得粘滞各向异性介质声波方程；
10.确定地下介质参数和地震波参数，并将所述地下介质参数和地震波参数带入到所述粘滞各向异性介质声波方程中，以计算获得地震波波场模拟数值。
11.作为可选地，所述第一刚度矩阵为：
[0012][0013]
并且，
[0014]
其中，v
p
表示纵波沿对称轴方向的速度，v
s
是横波速度，ρ为地下介质密度，ε和δ为速度各向异性参数，c
11
、c
13
、c
33
、c
55
和c
66
为刚度系数；
[0015]
相应地，所述第一品质因子矩阵为：
[0016][0017]
并且，
[0018]
其中，q
11
和q
33
分别对应纵波在水平方向和垂向上的品质因子，q
55
对应横波的品质因子，q
12
为品质因子剪切模量，q
12
＝q
13
,q
55
＝q
66
，ε
q
和δ
q
为衰减各向异性参数。
[0019]
进一步地，所述将所述第一刚度矩阵和所述第一品质因子矩阵进行简化处理，得到第二刚度矩阵和第二品质因子矩阵，包括：
[0020]
将横波速度和横波的品质因子设定为预设值；
[0021]
根据所述预设值，对所述第一刚度矩阵和所述第一品质因子矩阵进行简化处理，得到第二刚度矩阵和第二品质因子矩阵。
[0022]
进一步地，将所述横波速度和所述横波的品质因子均设定为零，则
[0023]
所述第二刚度矩阵为：
[0024][0025]
并且，
[0026]
其中，v
p
表示纵波沿对称轴方向的速度，v
s
是横波速度，ε和δ为速度各向异性参数；
[0027]
所述第二品质因子矩阵为：
[0028][0029]
并且，
[0030]
其中，q
11
和q
33
分别对应纵波在水平方向和垂向上的品质因子，ε
q
和δ
q
为衰减各向异性参数。
[0031]
进一步地，所述根据所述二刚度矩阵和所述第二品质因子矩阵，结合预设频散关系和运动平衡方程，计算获得粘滞各向异性介质声波方程，包括：
[0032]
确定常q模型近似频散关系和地震波运动平衡方程；
[0033]
将所述近似频散关系进行时间域变换，得到所述近似频散关系关于时间的一阶偏导数模型；
[0034]
将所述地震波运动平衡方程和几何方程带入到所述一阶偏导数模型中，并结合所述二刚度矩阵和所述第二品质因子矩阵，得到粘滞各向异性介质声波方程的一阶速度
‑
应力模型。
[0035]
进一步地，所述一阶速度
‑
应力模型为：
[0036][0037]
和为复模量，
[0038][0039][0040][0041]
其中，ω0为参考角频率，t表示时间，为拉普拉斯算子，σ
xx
和σ
zz
为t时刻x方向和z方向的地下介质质点的应力，v
x
和v
z
分别为地下介质质点沿x方向和z方向的振动速度分量，ρ为地下介质密度，γ
11
、γ
33
和γ
13
为空变分数阶阶数，γ
11
＝(1/π)arctan(1/q
11
)，γ
33
＝(1/π)arctan(1/q
33
)，γ
13
＝(1/π)arctan(1/q
13
)。
[0042]
进一步地，所述获得地震波波场模拟数值之后还包括：
[0043]
通过地下介质参数和地震波参数，结合克里斯托弗尔方程计算得到所述地震波的纵波复速度；
[0044]
根据所述纵波复速度，结合预设解析表达式，计算得到纵波品质因子参考值；
[0045]
根据地震波波场模拟数值，通过谱比法计算得到纵波品质因子估算值；
[0046]
根据所述纵波品质因子参考值和所述纵波品质因子估算值，计算获得所述粘滞各向异性介质声波方程的模拟精度。
[0047]
另一方面，本文还提供一种粘滞各向异性介质地震波数值模拟装置，所述装置包括：
[0048]
初始矩阵确定模块，用于根据已知弹性介质的第一刚度矩阵，确定粘弹性介质的第一品质因子矩阵；
[0049]
简化模块，用于将所述第一刚度矩阵和所述第一品质因子矩阵进行简化处理，得到第二刚度矩阵和第二品质因子矩阵；
[0050]
声波方程确定模块，用于根据所述二刚度矩阵和所述第二品质因子矩阵，结合预设频散关系，计算获得粘滞各向异性介质声波方程；
[0051]
模拟模块，用于确定地下介质参数和地震波参数，并将所述地下介质参数和地震
波参数带入到所述粘滞各向异性介质声波方程中，以计算获得地震波波场模拟数值。
[0052]
另一方面，本文还提供一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现如上述所述的方法。
[0053]
最后，本文还提供一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现如上述所述的方法。
[0054]
采用上述技术方案，本文所述的一种粘滞各向异性介质地震波数值模拟方法、装置及设备，根据已知弹性介质的第一刚度矩阵，确定粘弹性介质的第一品质因子矩阵；将所述第一刚度矩阵和所述第一品质因子矩阵进行简化处理，得到第二刚度矩阵和第二品质因子矩阵；根据所述二刚度矩阵和所述第二品质因子矩阵，结合预设频散关系、运动平衡方程和几何方程，计算获得粘滞各向异性介质声波方程；确定地下介质参数和地震波参数，并将所述地下介质参数和地震波参数带入到所述粘滞各向异性介质声波方程中，以计算获得地震波波场模拟数值，本文能够实现各向异性介质中地震波的速度和衰减随方位的变化特征，从而提高了地震波数值模拟的准确性。
[0055]
为让本文的上述和其他目的、特征和优点能更明显易懂，下文特举较佳实施例，并配合所附图式，作详细说明如下。
附图说明
[0056]
为了更清楚地说明本文实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本文的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
[0057]
图1示出了本文实施例提供的一种粘滞各向异性介质地震波数值模拟方法的步骤示意图；
[0058]
图2示出了本文实施例中声波方程表达式确定步骤示意图；
[0059]
图3示出了本文实施例中对声波方程精度验证步骤示意图；
[0060]
图4示出了本文实施例中声波方程的各向异性衰减声波波场；
[0061]
图5示出了本文实施例中设置震源附近为各向同性介质后所得的各向异性衰减声波波场；
[0062]
图6示出了本实施例中谱比法估算q值的观测系统；
[0063]
图7示出了本实施例中声波方程模拟精度验证对比示意图；
[0064]
图8示出了本实施例中2007bp tti模型示意图；
[0065]
图9示出了本实施例中不同模拟方式的波场快照示意图；
[0066]
图10示出了本文实施例中各向异性声波介质中的共炮点地震记录及不同方程的模拟结果示意图；
[0067]
图11示出了本实施例提供的一种粘滞各向异性介质地震波数值模拟装置的结构示意图；
[0068]
图12示出了本实施例提供的一种计算机设备的结构示意图。
[0069]
附图符号说明：
[0070]
100、初始矩阵确定模块；
[0071]
200、简化模块；
[0072]
300、声波方程确定模块；
[0073]
400、模拟模块；
[0074]
1202、计算机设备；
[0075]
1204、处理器；
[0076]
1206、存储器；
[0077]
1208、驱动机构；
[0078]
1210、输入/输出模块；
[0079]
1212、输入设备；
[0080]
1214、输出设备；
[0081]
1216、呈现设备；
[0082]
1218、图形用户接口；
[0083]
1220、网络接口；
[0084]
1222、通信链路；
[0085]
1224、通信总线。
具体实施方式
[0086]
下面将结合本文实施例中的附图，对本文实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本文一部分实施例，而不是全部的实施例。基于本文中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本文保护的范围。
[0087]
需要说明的是，本文的说明书和权利要求书及上述附图中的术语“第一”、“第二”等是用于区别类似的对象，而不必用于描述特定的顺序或先后次序。应该理解这样使用的数据在适当情况下可以互换，以便这里描述的本文的实施例能够以除了在这里图示或描述的那些以外的顺序实施。此外，术语“包括”和“具有”以及他们的任何变形，意图在于覆盖不排他的包含，例如，包含了一系列步骤或单元的过程、方法、装置、产品或设备不必限于清楚地列出的那些步骤或单元，而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。
[0088]
在地震勘探过程中，涉及到地下介质速度和衰减随方便变化的性质研究，即低下介质的各向异性，现有技术的研究主要集中在速度各向异性上，却忽略了吸收衰减各向异性对地震波传播的影响，现有技术中提出了一个各向异性时间分数阶粘滞弹性波方程，在该方程的刚度矩阵中，分别含有速度各向异性和衰减各向异性参数，所使用的速度各向异性参数是thomsen系数，衰减各向异性参数是由一个形式上类似的q矩阵来表征的，但是由于其中包含的时间分数阶算子在求解过程中需要存储大量的当前时刻以前的波场值，因此导致该方程需要较大的运算能力才能计算出准确的地震波变化，从而提高了数据计算的成本，进而导致上述方程实用性较差，很难在实际生产开发中运行。
[0089]
为了解决上述问题，本文实施例提供了一种粘滞各向异性介质地震波数值模拟方法，能够提高地震波数值模拟的准确性，同时降低了数值模拟的计算难度。图1是本文实施
例提供的一种粘滞各向异性介质地震波数值模拟方法的步骤示意图，本说明书提供了如实施例或流程图所述的方法操作步骤，但基于常规或者无创造性的劳动可以包括更多或者更少的操作步骤。实施例中列举的步骤顺序仅仅为众多步骤执行顺序中的一种方式，不代表唯一的执行顺序。在实际中的系统或装置产品执行时，可以按照实施例或者附图所示的方法顺序执行或者并行执行。具体的如图1所示，所述方法可以包括：
[0090]
s101：根据已知弹性介质的第一刚度矩阵，确定粘弹性介质的第一品质因子矩阵；
[0091]
s102：将所述第一刚度矩阵和所述第一品质因子矩阵进行简化处理，得到第二刚度矩阵和第二品质因子矩阵；
[0092]
s103：根据所述二刚度矩阵和所述第二品质因子矩阵，结合预设频散关系、运动平衡方程和几何方程，计算获得粘滞各向异性介质声波方程；
[0093]
s104：确定地下介质参数和地震波参数，并将所述地下介质参数和地震波参数带入到所述粘滞各向异性介质声波方程中，以计算获得地震波波场模拟数值。
[0094]
可以理解为，本说明书通过现有横向各向同性(vti)介质的刚度矩阵推导出vti粘弹性介质的品质因子(q)矩阵，然后为了降低q矩阵的复杂度，在对推导出的q矩阵进行简化处理，最后基于预设频散关系、运动平衡方程和几何方程可以得到粘滞各向异性介质声波方程，从而结合地下介质参数和地震波参数，模拟得到地震波波场，本文可以同时准确刻画地震波速度和衰减随方位变化的特性，同时在较小的计算量的情况下准确计算地震波在非均匀介质中的传播。
[0095]
在本说明书实施例中，由于声波方程的建立需要应力
‑
应变之间的对应关系，而应力
‑
应变关系是由与介质弹性性质密切相关的弹性刚度矩阵决定的。对于vti弹性介质，其voigt形式的弹性刚度矩阵的表达式为：
[0096][0097]
并且，
[0098]
其中，c为第一刚度矩阵，v
p
表示纵波(p波)沿对称轴方向的速度，v
s
是横波(s波)速度，ρ为地下介质密度，ε和δ为速度各向异性参数，c
11
、c
13
、c
33
、c
55
和c
66
为刚度系数，可以看出c
11
和c
33
与纵波速度相关，c
55
和c
66
与横波速度相关，c
33
可以理解为和纵波速度以及横波速度相关剪切模量。
[0099]
其中，地震波的纵波和横波可以为测量值或给定值，在模拟地震波波场过程中，可以设置一定的纵波和横波作为模拟的参数，由于地下介质是各向异性介质，因此ρ也是随着
地震波传播的位置和方位变化，ε和δ可以为thomsen系数，thomsen系数是各向异性研究中最常见的参数，一般来说，它需要岩石物理测量得到，在本说明书中可以根据地下介质的岩石学性质直接得到。
[0100]
对于弹性介质，公式(1)所示的第一刚度矩阵反映了速度随方位的变化。因此对于粘弹介质，也可以推导出一个类似的矩阵以便描述地层衰减随方向的变化，因此可以得到公式(2)中的第一品质因子(q)矩阵：
[0101][0102]
并且，
[0103]
其中，q为第一品质因子矩阵，q
11
和q
33
分别对应纵波在水平方向和垂向上的品质因子，q
55
对应横波的品质因子，q
12
为品质因子剪切模量，q
12
＝q
13
,q
55
＝q
66
，ε
q
和δ
q
为衰减各向异性参数。
[0104]
由于不同的岩性的品质因子可能不同，因此在各向异性介质中，品质因子也会随着地震波的位置或方位变化，可以通过确定介质中岩性的分布确定不同位置的品质因子的值，相应的，衰减各项异性参数ε
q
和δ
q
也可以根据实际地下介质的岩石岩性确定的，在本说明书中可以根据地下介质的岩石学性质直接得到。
[0105]
通过上述步骤得到的第一品质因子矩阵是包含纵波特征和横波特征的矩阵，即是包含地震波多个特征要素的矩阵模型，复杂的矩阵模型在后续运用和计算中会带来相当大的工作量，因此会有更高的计算性能要求，导致地震波波场模拟的性价比较低，不易大规模开展相应的研究，因此还需要对上述矩阵进行简化处理，以降低矩阵的复杂度。
[0106]
由于在各向异性介质中，针对矢量波场的研究要考虑质点的运动方向和振动方向，因此可以通过矢量波场到标量波场的近似处理，降低矩阵的复杂度。近似方法可以包括椭圆各向异性近似、弱各向异性近似以及声学近似。
[0107]
因此，在本说明书实施例中，所述将所述第一刚度矩阵和所述第一品质因子矩阵进行简化处理，得到第二刚度矩阵和第二品质因子矩阵，包括：
[0108]
将横波速度和横波的品质因子设定为预设值；
[0109]
根据所述预设值，对所述第一刚度矩阵和所述第一品质因子矩阵进行简化处理，得到第二刚度矩阵和第二品质因子矩阵。
[0110]
可以理解为，通过声学近似处理，可以将上述得到的第一刚度矩阵和第一品质因子矩阵进行简化(或退化)处理，得到维度较低的矩阵，从而降低了后续数据处理的难度。作
为优选地，可以设置与横波相关的变量值为零，即所述预设值为零(横波速度和横波的品质因子均为零。
[0111]
经过简化处理得到的第二刚度矩阵可以为：
[0112][0113]
并且，
[0114]
其中，vp表示纵波沿对称轴方向的速度，v
s
是横波速度，ε和δ为速度各向异性参数；
[0115]
相应地，所述第二品质因子矩阵可以为：
[0116][0117]
并且，
[0118]
其中，q
11
和q
33
分别对应纵波在水平方向和垂向上的品质因子，ε
q
和δ
q
为衰减各向异性参数。
[0119]
在得到第二刚度矩阵和所述第二品质因子矩阵的基础上，作为可选地，如图2所示，所述根据所述二刚度矩阵和所述第二品质因子矩阵，结合预设频散关系、运动平衡方程和几何方程，计算获得粘滞各向异性介质声波方程，包括：
[0120]
s201：确定常q模型近似频散关系和地震波运动平衡方程；
[0121]
s202：将所述近似频散关系进行时间域变换，得到所述近似频散关系关于时间的一阶偏导数模型；
[0122]
s203：将所述地震波运动平衡方程和几何方程带入到所述一阶偏导数模型中，并结合所述二刚度矩阵和所述第二品质因子矩阵，得到粘滞各向异性介质声波方程的一阶速度
‑
应力模型。
[0123]
其中，常q模型近似频散关系表示频率和速度的变化关系，其公式为：
[0124][0125]
其中，k为波数，ω为角频率，ω0为参考角频率，c为速度，γ＝(1/π)arctan(1/q)。
[0126]
将上述近似频散关系变换回时间域，得到关于时间的一阶偏导数模型：
[0127][0128]
其中，σ
ij
为t时刻地下介质质点在ij方向上的应力，为频散关系中的复模量，e
kk
为应变分量。
[0129]
地震波运动平衡方程和几何方程可以为：
[0130][0131]
其中，σ
ij
为t时刻地下介质质点在ij方向上的应力，v
i
和v
j
为地下介质质点沿i方向和j方向的振动速度分量，e
ij
为地震波在地下介质中ij方向的位移，x
i
和x
j
分别为地震波波场的空间坐标。
[0132]
通过将上述公式(7)带入到公式(6)中，结合简化后得到的第二刚度矩阵和第二品质因子矩阵，即可得到新频散关系的基于分数阶拉普拉斯算子(dfl)的粘滞各向异性介质声波方程，该方程的一阶速度
‑
应力模型为：
[0133][0134]
和为复模量，
[0135]
[0136]
其中，ω0为参考角频率，t表示时间，为拉普拉斯算子，σ
xx
和σ
zz
为t时刻x方向和z方向的地下介质质点的应力，v
x
和v
z
分别为地下介质质点沿x方向和z方向的振动速度分量，ρ为地下介质密度，γ
11
、γ
33
和γ
13
为空变分数阶阶数，γ
11
＝(1/π)arctan(1/q
11
)，γ
33
＝(1/π)arctan(1/q
33
)，γ
13
＝(1/π)arctan(1/q
13
)。
[0137]
在得到上述粘滞各向异性介质声波方程的基础上，通过确定低下介质参数(比如地下介质密度、品质因子等)，和地震波参数(纵波速度和横波速度等)，并将上述参数带入到粘滞各向异性介质声波方程从而实现了地震波波场模拟。本文是基于分数阶拉普拉斯算子的声波方程，可以同时准确刻画地震波速度和衰减随方位变化的特性，而且不包含随空间变化的分数阶拉普拉斯算子，因此可天然地准确计算地震波在非均匀介质中的传播，本文通过对声波方程的简化处理，可以实现各向异性介质声波方程的降维处理(比如三维到二维)，从而可以广泛应用于二维(2d)和三维(3d)复杂衰减各向异性介质中声波波场模拟、吸收衰减补偿逆时偏移成像及全波形反演。
[0138]
在一些其他实施例中，为了验证上述提供的粘滞各向异性介质声波方程的地震波模拟的准确性，如图3所示，所述获得地震波波场模拟数值之后还包括：
[0139]
s301：通过地下介质参数和地震波参数，结合克里斯托弗尔方程计算得到所述地震波的纵波复速度；
[0140]
s302：根据所述纵波复速度，结合预设解析表达式，计算得到纵波品质因子参考值；
[0141]
s303：根据地震波波场模拟数值，通过谱比法计算得到纵波品质因子估算值；
[0142]
s304：根据所述纵波品质因子参考值和所述纵波品质因子估算值，计算获得所述粘滞各向异性介质声波方程的模拟精度。
[0143]
可以理解为，本文通过对纵波品质因子理论数值的计算(即纵波品质因子参考值)，结合通过上述声波方程模拟得到的地震波波场进行谱比法计算得到的估算值进行比较，来确定上述声波方程的可靠性。
[0144]
其中纵波品质因子参考值通过克里斯托弗尔(christoffel)方程来计算纵波复速度，所述复速度的计算公式可以为：
[0145][0146]
中间变量a的表达式为：
[0147][0148]
其中，v
*p
为复速度，θ为传播方向与水平轴夹角的角度。
[0149]
进一步实施例中，所述预设解析表达式可以通过如下公式表示：
[0150]
[0151]
其中，q
p
为纵波品质因子参考值。
[0152]
为更好地验证所提出的衰减各向异性粘滞声波方程的精度，利用谱比法求取q值与解析解(即纵波品质因子参考值)进行对比。需要特别说明的是，所提出的粘滞各向异性介质声波方程是由各向异性粘滞弹性波方程通过声波近似得到的，虽然直接设置s波速度为零极其简便，但是波场中依然会残留有很强的横波伪影(如图4所示)，而设置震源附近为各向同性介质(将速度和衰减各向异性参数设置为零)，则可以在很大程度上去除这种伪影的干扰(如图5所示)。谱比法验证精度的具体做法是：设计一组在较大各向异性强度范围内变化的模型，利用所提出方程进行波场模拟，再通过谱比法获得其对应的q值。通过比较谱比法估计的q值与理论q值的匹配程度，即可判断出不同方案刻画衰减随方向变化的准确程度。
[0153]
在一具体实施例中，如图6所示，为所使用的观测系统，其中中心五角形代表炮点所在位置，围绕着炮点的两个同心圆上的三角形代表检波点，每个同心圆上检波点的数量为360个。利用这些检波点处的地震记录，可以通过谱比法估算出各个角度上的q值。分别求取同一角度上外侧和内侧检波点上地震数据的频谱，二者频谱比值的对数函数为：
[0154][0155]
其中，r1和r2分别为t1和t2时刻的地震波振幅谱，δt为内外两个检波器记录到地震信号的时间差。上式是以lng为截距，πδt/q为斜率k的一次函数，而所估算的q值的表达式为：
[0156]
q＝
‑
πδt/k
ꢀꢀ
(13)，
[0157]
如下表1为本文设计多个测试模型，在该组模型中，衰减各向异性强度在很大的范围内变化，其中模型1为衰减各向同性介质。其它参数为：模型大小为400
×
400个网格点，空间采样间隔10m，时间步长1ms，v
p
＝3.0km/s，thomsen速度各向异性参数为ε＝0.16，δ＝0.08，介质密度2000kg/m3。在图7中，圆形内部的黑色实线为求解公式(11)获得的理论值，圆圈所示为利用谱比法估算的q值(为便于观察，抽样为每5
°
显示一个数据)。需要说明的是，对于通过声波近似得到的方程，尽管在震源周围设置了各向同性介质，但是仍然会产生一定程度的伪影。为了排除这一干扰，逐道剔除了其中的伪横波，以此保证伪谱法估计q值的准确性。可以发现，在所有模型中，数值方法模拟的波场所估计的q值都能够很好地匹配理论值，这说明所提出的粘滞各向异性介质声波方程在刻画地层衰减随方向变化的精度是足够的。
[0158]
表1各向异性介质参数变化的测试模型
[0159]
测试模型q
11
q
13
q
33
ε
q
δ
q
品质因子(q)模型1202020000％模型23537500.430.7530％模型3504025
‑
0.5
‑
0.2250％模型4504035
‑
0.3
‑
0.6430％模型5504025
‑
0.5
‑
1.1650％模型6504015
‑
0.7
‑
1.6870％
[0160]
在另一具体实施例中，采用2007bp tti模型(不考虑倾斜角度)来验证更加复杂地质条件下所提出方法的模拟精度。图8分别展示了p波速度、p波品质因子、thomsen速度各向异性参数δ和ε，衰减各向异性参数为δ
q
＝2δ，ε
q
＝3ε。该模型大小为500
×
270个网格点，空间采样间隔10m，时间采样间隔1ms，模拟传播时间为2.5s，所采用的激发震源是主频为25hz的雷克子波。图9的a部分为各向异性弹性介质中的波场快照，图9的b部分由本文所提出的粘滞各向异性介质声波方程模拟得到的波场快照，图9的c部分由vti介质时间分数阶粘滞声波方程模拟得到的波场快照，该方程在向各向异性拓展时不需要近似，所以可以很好地描述地层衰减随方向的变化，但是由于其中包含的时间分数阶算子在数值求解过程中存在计算和存储的问题，因此不适用于大规模的数值模拟；与该方程模拟结果对比，即可验证所提出的vti介质dfl粘滞声波方程的模拟精度。图9的d部分为图9的b部分和图9的c部分相减的结果。图10为与图9对应的共炮点道集地震记录，分析对比后可以发现，与非衰减介质相比，地层衰减效应会明显使振幅能量减弱；此外，图9的d部分和图10的d部分中微弱的差值表明，所提出的粘滞各向异性介质声波方程对于复杂模型依然具有很高的数值模拟精度。
[0161]
基于同一发明构思，本文还提供一种粘滞各向异性介质地震波数值模拟装置，如图11所示，所述装置包括：
[0162]
初始矩阵确定模块100，用于根据已知弹性介质的第一刚度矩阵，确定粘弹性介质的第一品质因子矩阵；
[0163]
简化模块200，用于将所述第一刚度矩阵和所述第一品质因子矩阵进行简化处理，得到第二刚度矩阵和第二品质因子矩阵；
[0164]
声波方程确定模块300，用于根据所述二刚度矩阵和所述第二品质因子矩阵，结合预设频散关系、运动平衡方程和几何方程，计算获得粘滞各向异性介质声波方程；
[0165]
模拟模块400，用于确定地下介质参数和地震波参数，并将所述地下介质参数和地震波参数带入到所述粘滞各向异性介质声波方程中，以计算获得地震波波场模拟数值。
[0166]
通过上述装置取得的有益效果和上述方法取得的有意效果一致，本说明书实施例不做赘述。
[0167]
如图12所示，为本文实施例提供的一种计算机设备，所述计算机设备1202可以包括一个或多个处理器1204，诸如一个或多个中央处理单元(cpu)，每个处理单元可以实现一个或多个硬件线程。计算机设备1202还可以包括任何存储器1206，其用于存储诸如代码、设置、数据等之类的任何种类的信息。非限制性的，比如，存储器1206可以包括以下任一项或多种组合：任何类型的ram，任何类型的rom，闪存设备，硬盘，光盘等。更一般地，任何存储器都可以使用任何技术来存储信息。进一步地，任何存储器可以提供信息的易失性或非易失性保留。进一步地，任何存储器可以表示计算机设备1202的固定或可移除部件。在一种情况下，当处理器1204执行被存储在任何存储器或存储器的组合中的相关联的指令时，计算机设备1202可以执行相关联指令的任一操作。计算机设备1202还包括用于与任何存储器交互的一个或多个驱动机构1208，诸如硬盘驱动机构、光盘驱动机构等。
[0168]
计算机设备1202还可以包括输入/输出模块1210(i/o)，其用于接收各种输入(经由输入设备1212)和用于提供各种输出(经由输出设备1214))。一个具体输出机构可以包括呈现设备1216和相关联的图形用户接口(gui)1218。在其他实施例中，还可以不包括输入/输出模块1210(i/o)、输入设备1212以及输出设备1214，仅作为网络中的一台计算机设备。
计算机设备1202还可以包括一个或多个网络接口1220，其用于经由一个或多个通信链路1222与其他设备交换数据。一个或多个通信总线1224将上文所描述的部件耦合在一起。
[0169]
通信链路1222可以以任何方式实现，例如，通过局域网、广域网(例如，因特网)、点对点连接等、或其任何组合。通信链路1222可以包括由任何协议或协议组合支配的硬连线链路、无线链路、路由器、网关功能、名称服务器等的任何组合。
[0170]
对应于图1
‑
图3中的方法，本文实施例还提供了一种计算机可读存储介质，该计算机可读存储介质上存储有计算机程序，该计算机程序被处理器运行时执行上述方法的步骤。
[0171]
本文实施例还提供一种计算机可读指令，其中当处理器执行所述指令时，其中的程序使得处理器执行如图1至图3所示的方法。
[0172]
应理解，在本文的各种实施例中，上述各过程的序号的大小并不意味着执行顺序的先后，各过程的执行顺序应以其功能和内在逻辑确定，而不应对本文实施例的实施过程构成任何限定。
[0173]
还应理解，在本文实施例中，术语“和/或”仅仅是一种描述关联对象的关联关系，表示可以存在三种关系。例如，a和/或b，可以表示：单独存在a，同时存在a和b，单独存在b这三种情况。另外，本文中字符“/”，一般表示前后关联对象是一种“或”的关系。
[0174]
本领域普通技术人员可以意识到，结合本文中所公开的实施例描述的各示例的单元及算法步骤，能够以电子硬件、计算机软件或者二者的结合来实现，为了清楚地说明硬件和软件的可互换性，在上述说明中已经按照功能一般性地描述了各示例的组成及步骤。这些功能究竟以硬件还是软件方式来执行，取决于技术方案的特定应用和设计约束条件。专业技术人员可以对每个特定的应用来使用不同方法来实现所描述的功能，但是这种实现不应认为超出本文的范围。
[0175]
所属领域的技术人员可以清楚地了解到，为了描述的方便和简洁，上述描述的系统、装置和单元的具体工作过程，可以参考前述方法实施例中的对应过程，在此不再赘述。
[0176]
在本文所提供的几个实施例中，应该理解到，所揭露的系统、装置和方法，可以通过其它的方式实现。例如，以上所描述的装置实施例仅仅是示意性的，例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，例如多个单元或组件可以结合或者可以集成到另一个系统，或一些特征可以忽略，或不执行。另外，所显示或讨论的相互之间的耦合或直接耦合或通信连接可以是通过一些接口、装置或单元的间接耦合或通信连接，也可以是电的，机械的或其它的形式连接。
[0177]
所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部单元来实现本文实施例方案的目的。
[0178]
另外，在本文各个实施例中的各功能单元可以集成在一个处理单元中，也可以是各个单元单独物理存在，也可以是两个或两个以上单元集成在一个单元中。上述集成的单元既可以采用硬件的形式实现，也可以采用软件功能单元的形式实现。
[0179]
所述集成的单元如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本文的技术方案本质上或
者说对现有技术做出贡献的部分，或者该技术方案的全部或部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本文各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：u盘、移动硬盘、只读存储器(rom，read
‑
only memory)、随机存取存储器(ram，random access memory)、磁碟或者光盘等各种可以存储程序代码的介质。
[0180]
本文中应用了具体实施例对本文的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本文的方法及其核心思想；同时，对于本领域的一般技术人员，依据本文的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本文的限制。