一种用于三维ATI介质的弹性波场矢量分解方法及系统

一种用于三维ati介质的弹性波场矢量分解方法及系统
技术领域
1.本发明涉及地球物理数据处理和弹性参数反演领域，特别是关于一种用于三维ati(arbitrary transverse isotropy，任意横向各向同性)介质的弹性波场矢量分解方法及系统。


背景技术：

2.实际地下介质普遍存在各向异性特征，通常由定向排列的垂直裂缝或具有周期性的薄互层引起。当地震波穿过这些介质时，由于波长大于裂缝宽度或者薄层厚度，常出现地震各向异性特征。地层各向异性对称轴沿垂直方向时，可采用vti(verticaltransverse isotropy，垂直横向各向同性)模型近似。若经过强烈的构造运动，对称轴方向发生改变，不再沿着垂直方向，当倾斜对称轴在水平面的投影具有一定方位时，各向异性介质可以采用任意倾向的ti介质进行近似，此时对称轴的空间取向可以为任意方向，简写为ati。面对这种越来越复杂的各向异性油气储层，现有的地震勘探方法已经无法满足当前高精度油气勘探的需求。
3.目前，工业界常用的地震偏移方法主要包括射线类和波动方程类两类，基于弹性波动方程的逆时偏移方法由于不受倾角和偏移孔径的限制，可以较好地保存反映各向异性地层结构的地震波振幅和相位等信息，所以逐渐应用到实际勘探生产中。该方法可以在地下重建p波场和s波场，更符合真实的地震波传播过程，应用成像条件可以得到不同波模式的成像结果。但是在弹性各向异性地层中，耦合传播的p波场和s波场将产生串扰噪声，降低成像剖面的质量。特别是针对复杂的三维ati介质来说，更需要研究一种高效精确的弹性波场解耦方法，以提高弹性逆时偏移的成像精度。
4.传统弹性波场分解基于亥姆霍兹定理，通过旋度和散度算子实现p波场和s波场的分离，但是该算子认为波的偏振方向和传播方向一样，所以只能应用在各向同性介质中。面对各向异性介质，通过求解克里斯托弗方程可以计算实际波的偏振方向，进而实现p波场和s波场的分离。许多学者也基于这个思路，研究了弹性波场的矢量分解方法，通过在二维非均匀vti介质中，利用非稳态滤波实现p波场和s波场的分离。另外一种在均匀vti介质中的弹性波场分解方法，使分解前后的波场具有相同的振幅、相位和物理单位。随后低秩近似的引入大大简化了数值计算的过程，同时，有学者将该方法推广到二维非均匀vti介质中。
5.然而，当面对更复杂的各向异性介质时，例如三维ati介质，这些适用于vti模型的弹性波场分解方法由于无法考虑对称轴的倾角和方位信息，不能对弹性波场进行有效分离。另外，随着实际资料采集数量的增加，这些分解方法的计算成本依然很高，始终无法满足实际三维生产的要求。


技术实现要素：

6.针对上述问题，本发明的目的是提供一种用于三维ati介质的弹性波场矢量分解方法及系统，能够考虑对称轴的倾角和方位信息，且满足实际三维生产要求。
7.为实现上述目的，本发明采取以下技术方案：第一方面，提供一种用于三维ati介质的弹性波场矢量分解方法，包括：
8.对观测坐标系进行旋转，使得观测坐标系的垂向坐标轴平行于各向异性介质的倾斜对称轴方向，得到新观测坐标系；
9.基于设置的三维ati模型的弹性参数和thomsen参数，计算三维vti介质条件下的刚度矩阵，建立三维vti介质弹性波动方程，进行波场延拓；
10.在新观测坐标系下，采用三维vti介质弹性波动方程描述三维ati介质，确定空间域的三维ati分解算子；
11.在每一次时间迭代过程中，对弹性波场进行正反傅里叶变换，并基于设置的三维ati模型的弹性参数和thomsen参数，得到空间域的辅助波场；
12.根据空间域的三维ati分解算子和辅助波场，进行弹性波场分解，得到p波场和s波场。
13.进一步地，所述三维ati模型的弹性参数包括纵波速度v
p
、横波速度vs和地层密度ρ，所述thomsen参数包括纵波各向异性ε、纵横波各向异性过度参数δ和横波各向异性γ。
14.进一步地，所述刚度矩阵和thomson参数之间的关系为：
15.x
33
＝ρv
p2
16.c
44
＝ρv
s2
17.c
11
＝(1+2ε)ρv
p2
18.c
66
＝(1+2γ)ρv
s2
[0019][0020]
其中，c
11
、c
13
、c
33
、c
44
、c
66
为刚度矩阵中的独立参数；ρ为密度；v
p
、vs为沿着介质各向异性对称轴方向的p波速度和s波速度。
[0021]
进一步地，所述空间域的三维ati分解算子为：
[0022][0023]
其中，分别为对x轴、y轴和z轴的空间一阶偏导；r为系数；θ为倾斜对称轴在新观测坐标系中x-z平面内沿z轴顺时针方向的夹角；为倾斜对称轴在新观测坐标系中x-y平面内的投影沿x轴顺时针方向的夹角。
[0024]
进一步地，所述在每一次时间迭代过程中，对弹性波场进行正反傅里叶变换，并基于设置的三维ati模型的弹性参数和thomsen参数，得到空间域的辅助波场，包括：
[0025]
根据得到的三维ati分解算子，构建满足三维ati介质的泊松方程；
[0026]
在每一次时间迭代过程中，对弹性波场进行正反傅里叶变换，并基于设置的三维ati模型的弹性参数和thomsen参数以及构建的满足三维ati介质的泊松方程，计算空间域的辅助波场。
[0027]
进一步地，所述空间域的辅助波场w为：
[0028][0029]
其中，r为系数，u1、u2、u3、u4、u5、u6和u7分别为弹性波场u进行傅里叶变换后得到的波数域形式u构建的参数u1、u2、u3、u4、u5、u6和u7分别进行反傅里叶变换后得到的对应的空间域形式。
[0030]
进一步地，所述p波场u
p
和s波场us为：
[0031][0032][0033]
第二方面，提供一种用于三维ati介质的弹性波场矢量分解系统，包括：
[0034]
坐标系旋转模块，用于对观测坐标系进行旋转，使得观测坐标系的垂向坐标轴平行于各向异性介质的倾斜对称轴方向，得到新观测坐标系；
[0035]
刚度矩阵计算模块，用于基于设置的三维ati模型的弹性参数和thomsen参数，计算三维vti介质条件下的刚度矩阵，建立三维vti介质弹性波动方程，进行波场延拓；
[0036]
三维ati介质描述模块，用于在新观测坐标系下，采用三维vti介质弹性波动方程描述三维ati介质，确定空间域的三维ati分解算子；
[0037]
辅助波场确定模块，用于在每一次时间迭代过程中，对弹性波场进行正反傅里叶变换，并基于设置的三维ati模型的弹性参数和thomsen参数，得到空间域的辅助波场；
[0038]
弹性波场分解模块，用于根据空间域的三维ati分解算子和辅助波场，进行弹性波场分解，得到p波场和s波场。
[0039]
第三方面，提供一种处理设备，包括计算机程序指令，其中，所述计算机程序指令被处理设备执行时用于实现上述用于三维ati介质的弹性波场矢量分解方法对应的步骤。
[0040]
第四方面，提供一种计算机可读存储介质，所述计算机可读存储介质上存储有计算机程序指令，其中，所述计算机程序指令被处理器执行时用于实现上述用于三维ati介质的弹性波场矢量分解方法对应的步骤。
[0041]
本发明由于采取以上技术方案，其具有以下优点：
[0042]
1、在复杂三维各向异性介质中，本发明同时考虑各向异性以及倾角和方位角的影响，基于特征分析推导得到适用于三维ati介质的矢量波场分解，通过数值实例的结果，说明本发明能够有效得到分解干净的p波场和s波场。
[0043]
2、本发明构建泊松方程后，能够避免直接通过求解泊松方程计算三维ati介质的辅助波场，而采用一种新的方式，大大提高了计算效率以及在三维ati介质中弹性波场矢量分解的应用。
[0044]
3、本发明通过给定不同的倾角和方位角参数，带入推导的三维ati分解算子的公式，适用于大多数的ti介质，也可以退化到三维各向同性的情况，根据介质的弹性参数对本构方程的影响，针对不同弹性介质特征进行三维弹性波分解计算。
[0045]
4、本发明可以应用在其它的地球物理数据处理工作中，例如分离经过倾斜薄互层间多次反射的p波和s波，以及分离经过裂隙发育页岩气储层的混合弹性波等等，应用领域
及其广泛。
[0046]
综上所述，本发明可以广泛应用于地球物理数据处理和弹性参数反演领域中。
附图说明
[0047]
通过阅读下文优选实施方式的详细描述，各种其他的优点和益处对于本领域普通技术人员将变得清楚明了。附图仅用于示出优选实施方式的目的，而并不认为是对本发明的限制。在整个附图中，用相同的附图标记表示相同的部件。在附图中：
[0048]
图1是本发明一实施例提供的方法流程示意图。
具体实施方式
[0049]
下面将参照附图更详细地描述本发明的示例性实施方式。虽然附图中显示了本发明的示例性实施方式，然而应当理解，可以以各种形式实现本发明而不应被这里阐述的实施方式所限制。相反，提供这些实施方式是为了能够更透彻地理解本发明，并且能够将本发明的范围完整地传达给本领域的技术人员。
[0050]
应理解的是，文中使用的术语仅出于描述特定示例实施方式的目的，而无意于进行限制。除非上下文另外明确地指出，否则如文中使用的单数形式“一”、“一个”以及“所述”也可以表示包括复数形式。术语“包括”、“包含”、“含有”以及“具有”是包含性的，并且因此指明所陈述的特征、步骤、操作、元件和/或部件的存在，但并不排除存在或者添加一个或多个其它特征、步骤、操作、元件、部件、和/或它们的组合。文中描述的方法步骤、过程、以及操作不解释为必须要求它们以所描述或说明的特定顺序执行，除非明确指出执行顺序。还应当理解，可以使用另外或者替代的步骤。
[0051]
尽管可以在文中使用术语第一、第二、第三等来描述多个元件、部件、区域、层和/或部段，但是，这些元件、部件、区域、层和/或部段不应被这些术语所限制。这些术语可以仅用来将一个元件、部件、区域、层或部段与另一区域、层或部段区分开。除非上下文明确地指出，否则诸如“第一”、“第二”之类的术语以及其它数字术语在文中使用时并不暗示顺序或者次序。因此，以下讨论的第一元件、部件、区域、层或部段在不脱离示例实施方式的教导的情况下可以被称作第二元件、部件、区域、层或部段。
[0052]
需要说明的是：本发明实施例中提到的p波为一种纵波，粒子振动方向和波前进方平行，在所有地震波中，前进速度最快，也最早抵达。p波能在固体、液体或气体中传递。s波为一种横波，前进速度仅次于p波，粒子振动方向垂直于波的前进方向。
[0053]
本发明实施例提供的用于三维ati介质的弹性波场矢量分解方法及系统，考虑了复杂各向异性参数以及对称轴倾角和方位角，能够高效率地获得弹性波场分解得到的p波场和s波场。
[0054]
实施例1
[0055]
如图1所示，本实施例提供一种用于三维ati介质的弹性波场矢量分解方法，包括以下步骤：
[0056]
1)根据观测坐标系下三维ati介质的倾角和方位角信息，对观测坐标系进行旋转，使得观测坐标系的垂向坐标轴平行于各向异性介质的倾斜对称轴方向，得到新观测坐标系。
[0057]
具体地，此时，原观测坐标系下的三维ati介质可看为新观测坐标系下的三维vti介质。
[0058]
2)基于设置的三维ati模型的弹性参数和thomsen参数，计算三维vti介质条件下的刚度矩阵c
ijkl
，建立三维vti介质弹性波动方程，进行波场延拓，具体为：
[0059]
2.1)基于prem地球物理模型(初步地球参考模型)，设置三维ati模型(任意倾向横向各向同性模型)的弹性参数和thomsen参数(汤姆森参数)。
[0060]
具体地，三维ati模型的弹性参数包括纵波速度v
p
、横波速度vs和地层密度ρ，thomsen参数包括纵波各向异性ε、纵横波各向异性过度参数δ和横波各向异性γ。
[0061]
具体地，本发明可以参照地球物理勘探中熟知的prem地球物理模型来设置参数，也不仅限于该模型，例如可以参照智利和阿根廷中部区域的远震p波层析成像模型或iasp91地球物理模型。
[0062]
2.2)基于设置的三维ati模型的弹性参数和thomsen参数，计算三维vti介质条件下的刚度矩阵c
ijkl
。
[0063]
具体地，根据弹性张量的voigt对称性，四阶张量c
ijkl
中的下角标ij(kl)＝11、22、33、12、23、31可分别对应为i(j)＝1、2、3、4、5、6，于是原本81个独立常数可降为36个，记为刚度矩阵c
ijkl
。
[0064]
刚度矩阵和thomson参数之间的关系为：
[0065]c33
＝ρv
p2
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0066]c44
＝ρv
s2
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(2)
[0067]c11
＝(1+2ε)ρv
p2
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0068]c66
＝(1+2γ)ρv
s2
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0069][0070]
其中，c
11
、c
13
、c
33
、c
44
、c
66
为刚度矩阵中的独立参数；ρ为密度；v
p
、vs为沿着介质各向异性对称轴方向的p波速度和s波速度。
[0071]
2.3)根据计算的刚度矩阵c
ijkl
，建立三维vti介质弹性波动方程，进行波场延拓，模拟波场传播的过程。
[0072]
具体地，将计算的刚度矩阵c
ijkl
带入三维vti介质弹性波动方程中建立针对该模型的弹性波动方程，输入震源的三分量位移数据u
x
、uy、uz，从零时刻开始，模拟波场传播的过程。
[0073]
建立的弹性波动方程为：
[0074][0075][0076][0077]
其中，u
x
、uy、uz为三分量位移；c
11
、c
13
、c
33
、c
44
、c
66
为刚度矩阵中的独立参数；分别为二阶和混合空间偏导。
[0078]
3)在新观测坐标系下，采用三维vti介质弹性波动方程描述三维ati介质，确定空间域的三维ati分解算子。
[0079]
具体地，在新观测坐标系下，采用三维vti介质弹性波动方程描述三维ati介质，针对三维ati介质，定义θ为倾斜对称轴在新观测坐标系中x-z平面内沿z轴顺时针方向的夹角，为倾斜对称轴在新观测坐标系中x-y平面内的投影沿x轴顺时针方向的夹角，则空间域的三维ati分解算子为：
[0080][0081]
其中，分别为对x轴、y轴和z轴的空间一阶偏导；r为系数，且：
[0082][0083]
4)在每一次时间迭代过程中，对弹性波场进行正反傅里叶变换，并基于设置的三维ati模型的弹性参数和thomsen参数，得到空间域的辅助波场，具体为：
[0084]
4.1)基于亥姆霍兹定理，根据得到的三维ati分解算子构建满足三维ati介质的泊松方程：
[0085][0086]
其中，u为弹性波场，w为辅助波场。
[0087]
具体地，根据亥姆霍兹定理，任何弹性波场均可分解为一个无源场和一个无旋场，即若已知三维ati分解算子则需要通过计算一个辅助波场w，进行弹性波场分解得到分解的p波场和s波场。
[0088]
4.2)在每一次时间迭代过程中，对弹性波场u进行傅里叶变换，得到弹性波场u的波数域形式u：
[0089]
u＝-(r1k
x2
+r2k
y2
+r3k
z2
+r4k
xky
+r5k
xkz
+r6k
ykz
)w
ꢀꢀꢀꢀꢀꢀꢀꢀ
(12)
[0090]
其中，w为辅助波场w的频率域表示形式，r1、r2、r3、r4、r5、r6均为系数，且：
[0091][0092][0093]
r3＝sin2θ+r2cos2θ
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(15)
[0094]
r4＝2 sinθcosθ(1-cos2θ-r2sin2θ)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(16)
[0095][0096][0097]
4.3)根据弹性波场u的波数域形式u，构建参数u1、u2、u3、u4、u5、u6和u7：
[0098]
[0099][0100][0101][0102][0103][0104][0105]
其中，k
x
、ky和kz分别为x、y和z方向的波数。
[0106]
4.4)对参数u1、u2、u3、u4、u5、u6和u7进行反傅里叶变换，得到对应的空间域形式u1、u2、u3、u4、u5、u6和u7。
[0107]
4.5)基于设置的三维ati模型的弹性参数和thomsen参数，根据构建的三维ati介质的泊松方程以及得到的参数u1、u2、u3、u4、u5、u6和u7及其对应的空间域形式u1、u2、u3、u4、u5、u6和u7，计算空间域的辅助波场w。
[0108]
具体地，为避免直接求解该泊松方程，提高计算效率，采用快速算法确定辅助波场w，在空间域的辅助波场w的计算公式为：
[0109][0110]
其中，u1、u2、u3、u4、u5、u6、u7分别为参数u1、u2、u3、u4、u5、u6、u7的时间空间域表示形式；r为弹性系数。且：
[0111][0112]
5)采用矢量分解法，根据得到的空间域的三维ati分解算子和空间域的辅助波场w，进行弹性波场分解，得到p波场u
p
和s波场us：
[0113][0114][0115]
其中，
·
表示点乘积；
×
表示叉乘积。
[0116]
上述步骤3)的具体推导过程为：
[0117]
将三维vti介质弹性波动方程变换到频率波数域写为矩阵形式：
[0118]
[0119][0120]
其中，为三分量位移u
x
、uy、uz傅里叶变换后的形式；ω为角频率；k
x
、ky、kz为三个对称轴方向的波数。
[0121]
在频散关系(v为相速度)的基础上，将等式(30)两边同除波数的模后，简写为：
[0122][0123]
其中，a为公式(30)中左侧的3
×
3矩阵，为列向量。
[0124]
对公式(31)进行特征分析，矩阵a的特征值为ρω2，可以表示不同波的相速度矩阵a的特征向量为表示波的偏振方向。
[0125]
因此在空间域的偏振方向公式为：
[0126][0127][0128][0129]
其中，a1、a2、a3分别为p波偏振方向、sh波偏振方向和sv波偏振方向。
[0130]
三维vti分解算子公式为：
[0131][0132]
针对三维ati介质，将观测坐标系进行旋转，使得旋转后的新观测坐标系的z轴沿倾斜对称轴方向，在三维空间下的坐标旋转公式为：
[0133][0134]
构建新观测坐标系下的空间偏导数因此得到三维任意倾向ti
介质的分解算子计算公式为上述公式(9)。
[0135]
上述步骤4)的具体推导过程为：
[0136]
根据泊松方程求解空间域的辅助波场w，对上述公式(11)展开并进行傅里叶变换至频率波数域得到上述公式(12)。
[0137]
令系数因此上述公式(12)简化为：
[0138]-f(ε,δ)u＝w
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(37)
[0139]
基于弱各向异性介质的假设，将系数f(ε,δ)在ε＝0和δ＝0处进行taylor(泰勒)展开为：
[0140][0141]
其中，r、e均为弹性系数，且弹性系数r为上述公式(27)，弹性系数e为：
[0142][0143]
将系数f(ε,δ)带入上述公式(37)整理得到：
[0144][0145]
代表了频率域中与不同波数进行计算后的弹性波场，最后变换至空间域，得到辅助波场w的计算公式(11)。
[0146]
实施例2
[0147]
本实施例提供一种用于三维ati介质的弹性波场矢量分解系统，包括：
[0148]
坐标系旋转模块，用于根据观测坐标系下三维ati介质的倾角和方位角信息，对观测坐标系进行旋转，使得观测坐标系的垂向坐标轴平行于各向异性介质的倾斜对称轴方向，得到新观测坐标系。
[0149]
刚度矩阵计算模块，用于基于设置的三维ati模型的弹性参数和thomsen参数，计算三维vti介质条件下的刚度矩阵，建立三维vti介质弹性波动方程，进行波场延拓。
[0150]
三维ati介质描述模块，用于在新观测坐标系下，采用三维vti介质弹性波动方程描述三维ati介质，确定空间域的三维ati分解算子。
[0151]
辅助波场确定模块，用于在每一次时间迭代过程中，对弹性波场进行正反傅里叶变换，并基于设置的三维ati模型的弹性参数和thomsen参数，得到空间域的辅助波场。
[0152]
弹性波场分解模块，用于根据空间域的三维ati分解算子和辅助波场，进行弹性波场分解，得到p波场和s波场。
[0153]
实施例3
[0154]
本实施例提供一种与本实施例1所提供的用于三维ati介质的弹性波场矢量分解方法对应的处理设备，处理设备可以是用于客户端的处理设备，例如手机、笔记本电脑、平
板电脑、台式机电脑等，以执行实施例1的方法。
[0155]
所述处理设备包括处理器、存储器、通信接口和总线，处理器、存储器和通信接口通过总线连接，以完成相互间的通信。存储器中存储有可在处理设备上运行的计算机程序，处理设备运行计算机程序时执行本实施例1所提供的用于三维ati介质的弹性波场矢量分解方法。
[0156]
在一些实现中，存储器可以是高速随机存取存储器(ram：random access memory)，也可能还包括非不稳定的存储器(non-volatile memory)，例如至少一个磁盘存储器。
[0157]
在另一些实现中，处理器可以为中央处理器(cpu)、数字信号处理器(dsp)等各种类型通用处理器，在此不做限定。
[0158]
实施例4
[0159]
本实施例提供一种与本实施例1所提供的用于三维ati介质的弹性波场矢量分解方法对应的计算机程序产品，计算机程序产品可以包括计算机可读存储介质，其上载有用于执行本实施例1所述的用于三维ati介质的弹性波场矢量分解方法的计算机可读程序指令。
[0160]
计算机可读存储介质可以是保持和存储由指令执行设备使用的指令的有形设备。计算机可读存储介质例如可以是但不限于电存储设备、磁存储设备、光存储设备、电磁存储设备、半导体存储设备或者上述的任意组合。
[0161]
上述各实施例仅用于说明本发明，其中各部件的结构、连接方式和制作工艺等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。