本发明属于石油勘探技术领域，具体涉及一种黏弹各向异性双相介质区域变网格求解算子。
背景技术：
地震波波场正演模拟越来越受到地球物理学者们的重视，它的地位也显得越来越重要。其中，地球介质具有黏弹性、各向异性和双相性，目前针对地震波波场正演模拟大都基于以上其中一种或两种性质考虑。但是，这种基于一种或两种性质模拟得到的地震波数据与基于三种性质模拟出的多分量的地震波数据在流体识别与描述、裂缝分布估计、岩性估计、各向异性分析、亮点反射设别等方面都有所差距，基于三种性质模拟出的多分量的地震波数据更加有效。
虽然基于三种性质模拟出的多分量的地震波数据更加有效，但是全地下构造均采用黏弹各向异性双相介质弹性波正演模拟方法进行地震波波场模拟时，模拟计算量巨大，现有的计算机硬件难以承受。因此，本发明提出了一种黏弹各向异性双相介质区域变网格求解算子，既同时考虑了地球介质的黏弹性、各向异性和双相性，又只针对目标区域进行弹性纵横波场的模拟，从而降低模拟的计算量。
技术实现要素：
本发明的目的是为克服上述现有技术的不足，提供一种黏弹各向异性双相介质区域变网格求解算子。
为实现上述目的，本发明采用如下技术方案：
一种黏弹各向异性双相介质区域变网格求解算子，包括以下步骤：
步骤1：输入速度场及各种参数场，建立观测系统；
步骤2：选择目标区域进行局部网格加密；
步骤3：在全局区域采用黏声各向异性介质进行波场延拓；
步骤4：在目标区域采用黏弹各向异性双相介质进行波场延拓；
步骤5：在目标区域的边界处进行黏声各向异性和黏弹各向异性双相波场传递；
步骤6：应用人工边界条件对边界反射进行吸收；
步骤7：输出炮记录和波场快照。
优选的，所述步骤2中，在全局区域进行粗网格剖分，在目标区域进行细网格剖分。
优选的，所述步骤2中，全局区域的粗网格尺寸为目标区域细网格尺寸的奇数倍。
优选的，所述步骤3中，在全局区域采用黏声各向异性波动方程(1)进行p波传播；
方程(1)中，p表示各向异性声波波场；q表示各向异性声波辅助波场；vpz表示对称轴方向的相速度；ε表示用以表征纵波各向异性强度的thomson参数，δ表示用以表征动校正速度的thomson参数；t表示时间；fx表示水平分量的震源；fz表示垂直分量的震源；x表示空间水平分量的空间坐标；z表示空间垂直分量的空间坐标；τ表示松弛时间，τ由公式(2)求得：
τ＝τε/τσ-1(2)
公式(2)中，τσ为应力松弛时间，τε为应变松弛时间；其中，τσ和τε与纵波品质因子qp的关系如公式(3)、公式(4)所示；
公式(3)和公式(4)中，ω代表角频率。
优选的，所述步骤3中全局区域采用的黏声各向异性波动方程(1)由公式(5)递推求得；
公式(5)中，n表示差分计算次数；f和f-1分别为傅里叶变换与傅里叶反变换；δt是时间步长；kx为水平分量的波数；kz为垂直分量的波数；
公式(6)～公式(9)中，δx表示水平分量的空间网格步长；δz表示垂直分量的空间网格步长；c2,0和c2,m为二阶偏导数的有限差分系数；m为积分变量；m为差分阶数；i为水平分量的空间坐标计数，j为垂直分量的空间坐标计数。
优选的，所述步骤4中，在目标区域采用黏弹各向异性双相介质一阶速度应力方程(10)进行多分量波场延拓；
方程(10)中，ux为固相水平分量的质点速度；uz为固相垂直分量的质点速度；ux为流相水平分量的质点速度；uz为流相垂直分量的质点速度；τxx为水平分量的正应力；τzz为垂直分量的正应力；τxz表示切应力；s表示流相正应力；q为固体与流体之间相耦合弹性系数；r为流相弹性系数；b表示耗散系数；t表示时间；x表示空间水平分量的空间坐标；z表示空间垂直分量的空间坐标；l1表示体积模量的标准线性固体个数；l2表示剪切模量的标准线性固体个数；
方程(10)中，d1、d2和d3由公式(11)求得：
公式(11)中，ρ11表示单位体积中固体相对流体运动时固体部分总的等效质量，ρ22表示单位体积中流体相对固体运动时流体部分总的等效质量，ρ12表示单位体积中流体和固体之间的质量耦合系数；
方程(10)中，
公式(12)中，
方程(10)中，e11、e55由公式(13)求得：
公式(13)中，ε表示用以表征纵波各向异性强度的thomson参数；ρ表示密度；vp表示纵波速度；vs表示横波速度；
方程(10)中，
公式(14)中，λ为与纵横波速度和密度有关的拉梅常数；μ为与横波速度和密度有关的拉梅常数；
优选的，所述步骤4中，目标区域采用的黏弹各向异性双相介质一阶速度应力方程(10)由公式(16)递推求得；
公式(16)中，δt1为目标区域的时间采样间隔，δt1与时间步长δt之间的关系如公式(17)所示；
δt＝(2k 1)δt1(17)
公式(17)中，k为任一正整数；
公式(16)中，dx表示水平分量的空间一阶偏导数，dz表示垂直分量的空间一阶偏导数，dx和dz由公式(18)求得：
公式(18)中，u为波场值，可表示ux、uz、ux、uz、τxx、τzz、τxz、
优选的，所述步骤5中，在目标区域的左右边界采用公式(20)所示的边界条件进行声压p与应力τxx的波场传递；
-p＝τxx(20)
在目标区域的上下边界采用公式(21)进行声压p与应力τzz的波场传递；
-p＝τzz(21)。
优选的，所述步骤5中，目标区域边界处全局区域延拓的黏声各向异性介质声波波场与目标区域延拓的黏弹各向异性双相介质弹性波波场互相传递的步骤如下：
步骤51：将目标区域上下边界上处于粗网格点上的声压p传递给对应的细网格点上的τxx，粗网格点上的uz传递给对应细网格点上的uz；将目标区域左右边界上处于粗网格点上的声压p传递给对应的细网格点上的τzz，粗网格点上的ux传递给对应细网格点上的ux；
步骤52：利用插值公式计算目标区域边界上其他细网格点上的波场值；
步骤53：采用降阶方法求解过渡区域的波场值；
步骤54：在目标区域其他细网格点黏弹各向异性双相介质波动方程(16)进行波场延拓；
步骤55：将目标区域的过渡区域上位于细网格点上的波场值传递给对应的粗网格；
步骤56：将目标区域上下边界上位于粗网格位置上的细网格点波场值τxx传递给粗网格点波场值p，将目标区域左右边界上位于粗网格位置上的细网格点波场值τzz传递给粗网格点波场值p，位于粗网格位置上的细网格点波场值ux和uz传递给粗网格点波场值ux和uz。
优选的，所述步骤6中，在全局区域采用pml边界条件对波场的衰减进行吸收，对目标区域的边界采用cfs-pml边界条件进行波场延拓；
使用拉伸函数将速度-应力方程变换到拉伸坐标域；在频率域中，2d拉伸函数的具体表述形式：
式中，σx代表cfs-pml中的水平分量衰减系数；σz代表cfs-pml中的垂直分量衰减系数；εx为水平分量的拉伸函数；εz为垂直分量的拉伸函数；ω为角频率；i表示复数虚单位；
笛卡尔坐标系(x，z)与拉伸坐标系
将方程(24)代入方程(10)，得到拉伸坐标系下一阶速度-应力方程的表达式：
其中，
本发明的有益效果是：
本发明黏弹各向异性双相介质区域变网格求解算子，同时考虑了地下介质的黏弹性、各向异性和双相性，在全局区域采用黏声各向异性介质的声波算子进行波场延拓，在目标区域采用黏弹各向异性双相介质的弹性波算子进行延拓，在目标区域的边界上进行稳定的声波和弹性波场传递，并在目标区域采用细网格剖分，提高了地震波场的模拟准确度并提高了计算效率。
附图说明
构成本申请的一部分的说明书附图用来提供对本申请的进一步理解，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。
图1是本发明黏弹各向异性双相介质区域变网格求解算子的结构示意图；
图2(a)为实施例中黏弹各向异性双相介质模型的双相介质参数；
图2(b)为实施例中黏弹各向异性双相介质模型的速度密度参数；
图2(c)为实施例中黏弹各向异性双相介质模型的品质因子参数；
图2(d)为实施例中黏弹各向异性双相介质模型的各向异性参数；
图3(a)为采用本发明模拟得到的地震波正演模拟炮记录的固相水平分量图；
图3(b)为采用本发明模拟得到的地震波正演模拟炮记录的固相垂直分量图；
图3(c)为采用本发明模拟得到的地震波正演模拟炮记录的流相水平分量图；
图3(d)为采用本发明模拟得到的地震波正演模拟炮记录的流相垂直分量图；
图4(a)为采用本发明模拟得到的波场快照的固相水平分量图；
图4(b)为采用本发明模拟得到的波场快照的固相垂直分量图；
图4(c)为采用本发明模拟得到的波场快照的流相水平分量图；
图4(d)为采用本发明模拟得到的波场快照的流相垂直分量图；
图5(a)为只考虑了各向异性、双相性的各向异性双相介质地震波正演模拟炮记录的固相水平分量图；
图5(b)为只考虑了各向异性、双相性的各向异性双相介质地震波正演模拟炮记录的固相垂直分量图；
图5(c)为只考虑了各向异性、双相性的各向异性双相介质地震波正演模拟炮记录的流相水平分量图；
图5(d)为只考虑了各向异性、双相性的各向异性双相介质地震波正演模拟炮记录的流相垂直分量图；
图6(a)为只考虑了黏弹性、双相性的黏弹各向同性双相介质地震波正演模拟炮记录的固相水平分量图；
图6(b)为只考虑了黏弹性、双相性的黏弹各向同性双相介质地震波正演模拟炮记录的固相垂直分量图；
图6(c)为只考虑了黏弹性、双相性的黏弹各向同性双相介质地震波正演模拟炮记录的流相水平分量图；
图6(d)为只考虑了黏弹性、双相性的黏弹各向同性双相介质地震波正演模拟炮记录的流相垂直分量图。
具体实施方式
应该指出，以下详细说明都是例示性的，旨在对本申请提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本申请所属技术领域的普通技术人员通常理解的相同含义。
需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。
下面结合附图和实施例对本发明进一步说明。
如图1所示，一种黏弹各向异性双相介质区域变网格求解算子，包括以下步骤：
步骤1：输入速度场及各种参数场，建立观测系统；
具体地，输入纵波速度场、横波速度场、密度场、耗散系数场、流相弹性系数场、固体与流体之间相耦合弹性系数场、纵波品质因子场qp、横波品质因子场qs、各向异性参数场；
步骤2：在全局区域进行粗网格剖分，在目标区域进行细网格剖分，从而实现整个计算区域的非均匀网格划分；
具体地，全局区域的粗网格尺寸为目标区域细网格尺寸的奇数倍，如3、5、7等倍数；
步骤3：在全局区域采用黏声各向异性介质进行波场延拓；
具体地，所述步骤3中，在全局区域采用黏声各向异性波动方程(1)进行p波传播；
方程(1)中，p表示各向异性声波波场；q表示各向异性声波辅助波场；vpz表示对称轴方向的相速度；ε表示用以表征纵波各向异性强度的thomson参数，δ表示用以表征动校正速度的thomson参数；t表示时间；fx表示水平分量的震源；fz表示垂直分量的震源；x表示空间水平分量的空间坐标；z表示空间垂直分量的空间坐标；τ表示松弛时间，τ由公式(2)求得：
τ＝τε/τσ-1(2)
公式(2)中，τσ为应力松弛时间，τε为应变松弛时间；其中，τσ和τε与纵波品质因子qp的关系如公式(3)、公式(4)所示；
公式(3)和公式(4)中，ω代表角频率。
具体地，所述步骤3中全局区域采用的黏声各向异性波动方程由公式(5)递推求得；
公式(5)中，n表示差分计算次数；f和f-1分别为傅里叶变换与傅里叶反变换；δt是时间步长；kx为水平分量的波数；kz为垂直分量的波数；
公式(6)～公式(9)中，δx表示水平分量的空间网格步长；δz表示垂直分量的空间网格步长；c2,0和c2,m为二阶偏导数的有限差分系数；m为积分变量；m为差分阶数；i为水平分量的空间坐标计数，j为垂直分量的空间坐标计数。
步骤4：在目标区域采用黏弹各向异性双相介质进行波场延拓；
具体地，所述步骤4中，在目标区域采用黏弹各向异性双相介质一阶速度应力方程(10)进行多分量波场延拓；
方程(10)中，ux为固相水平分量的质点速度；uz为固相垂直分量的质点速度；ux为流相水平分量的质点速度；uz为流相垂直分量的质点速度；τxx为水平分量的正应力；τzz为垂直分量的正应力；τxz表示切应力；s表示流相正应力；q为固体与流体之间相耦合弹性系数；r为流相弹性系数；b表示耗散系数；t表示时间；x表示空间水平分量的空间坐标；z表示空间垂直分量的空间坐标；l1表示体积模量的标准线性固体个数；l2表示剪切模量的标准线性固体个数；
方程(10)中，d1、d2和d3由公式(11)求得：
公式(11)中，ρ11表示单位体积中固体相对流体运动时固体部分总的等效质量，ρ22表示单位体积中流体相对固体运动时流体部分总的等效质量，ρ12表示单位体积中流体和固体之间的质量耦合系数；
方程(10)中，
公式(12)中，
方程(10)中，e11、e55由公式(13)求得：
公式(13)中，ε表示用以表征纵波各向异性强度的thomson参数；ρ表示密度；vp表示纵波速度；vs表示横波速度；
方程(10)中，
公式(14)中，λ为与纵横波速度和密度有关的拉梅常数；μ为与横波速度和密度有关的拉梅常数；
具体地，所述步骤4中，目标区域采用的黏弹各向异性双相介质一阶速度应力方程(10)由公式(16)递推求得；
公式(16)中，δt1为目标区域的时间采样间隔，δt1与时间步长δt之间的关系如公式(17)所示；
δt＝(2k 1)δt1(17)
公式(17)中，k为任一正整数；
公式(16)中，dx表示水平分量的空间一阶偏导数，dz表示垂直分量的空间一阶偏导数，dx和dz由公式(18)求得：
公式(18)中，u为波场值，可表示ux、uz、ux、uz、τxx、τzz、τxz、
步骤5：在目标区域的边界处进行黏声各向异性和黏弹各向异性双相波场传递；
具体地，所述步骤5中，在目标区域的左右边界采用公式(20)所示的边界条件进行声压p与应力τxx的波场传递；
-p＝τxx(20)
在目标区域的上下边界采用公式(21)进行声压p与应力τzz的波场传递；
-p＝τzz(21)
具体地，所述步骤5中，目标区域边界处全局区域延拓的黏声各向异性介质声波波场与目标区域延拓的黏弹各向异性双相介质弹性波波场互相传递的步骤如下：
步骤51：将目标区域上下边界上处于粗网格点上的声压p传递给对应的细网格点上的τxx，粗网格点上的uz传递给对应细网格点上的uz；将目标区域左右边界上处于粗网格点上的声压p传递给对应的细网格点上的τzz，粗网格点上的ux传递给对应细网格点上的ux；
步骤52：利用插值公式计算目标区域边界上其他细网格点上的波场值；
步骤53：采用降阶方法求解过渡区域的波场值；如采用10阶差分精度，在边界区域逐渐采用8阶、6阶、4阶、2阶差分精度；
步骤54：在目标区域其他细网格点黏弹各向异性双相介质波动方程(16)进行波场延拓；
步骤55：将目标区域的过渡区域上位于细网格点上的波场值传递给对应的粗网格；
步骤56：将目标区域上下边界上位于粗网格位置上的细网格点波场值τxx传递给粗网格点波场值p，将目标区域左右边界上位于粗网格位置上的细网格点波场值τzz传递给粗网格点波场值p，位于粗网格位置上的细网格点波场值ux和uz传递给粗网格点波场值ux和uz。
在进行波场传递时，直接一对一传递会引起不稳定，因此本发明采用的是九点加权法，利用粗网格点周围若干个细网格点波场值计算该点相应的粗网格点波场值；变时间步长的基本原理是在全局区域采用大尺度时间步长，在目标区域采用小尺度时间步长，边界区域的小尺度时间步长利用双线性插值计算得到；
步骤6：应用人工边界条件对边界反射进行吸收；
具体地，所述步骤6中，在全局区域采用pml边界条件对波场的衰减进行吸收，对目标区域的边界采用cfs-pml边界条件进行波场延拓；
使用拉伸函数将速度-应力方程变换到拉伸坐标域；在频率域中，2d拉伸函数的具体表述形式：
式中，σx代表cfs-pml中的水平分量衰减系数；σz代表cfs-pml中的垂直分量衰减系数；εx为水平分量的拉伸函数；εz为垂直分量的拉伸函数；ω为角频率；i表示复数虚单位；
笛卡尔坐标系(x，z)与拉伸坐标系
将方程(24)代入方程(10)，得到拉伸坐标系下一阶速度-应力方程的表达式：
其中，
步骤7：输出炮记录和波场快照。
本发明公开了一种黏弹各向异性双相介质区域变网格求解算子，该方法在全局区域采用黏声各向异性介质的声波算子进行波场延拓，在目标区域采用黏弹各向异性双相介质的弹性波算子进行延拓，在目标区域的边界上进行稳定的声波和弹性波场传递，并在目标区域采用细网格剖分，提高了地震波场的模拟准确性并提高了计算效率。
实施例：
本发明黏弹各向异性双相介质区域变网格求解算子，应用于黏弹各向异性双相介质模型，取得了理想的计算效果。
首先输入速度场及各种参数场，建立观测系统，其中实施例中黏弹各向异性双相介质模型如图2(a)～图2(d)所示，图2(a)为黏弹各向异性双相介质模型的双相介质参数，图2(b)为黏弹各向异性双相介质模型的速度密度参数，图2(c)为黏弹各向异性双相介质模型的品质因子参数，图2(d)为黏弹各向异性双相介质模型的各向异性参数。
然后采用本发明黏弹各向异性双相介质区域变网格求解算子进行地震波正演模拟，模型大小为3km×3km。在实施例中，全局网格间距为10m×10m，全局时间采样间隔为1ms，其中，全局区域的粗网格尺寸为目标区域细网格尺寸的3倍；记录时间为1.5s，炮点位于地表位置激发，水平位置为1500m，检波点均匀分布于地表处，共300个检波点。
采用本发明模拟得到的地震波正演模拟炮记录如图3(a)～图3(d)所示，其中图3(a)为固相水平分量，图3(b)为固相垂直分量，图3(c)为流相水平分量，图3(d)为流相垂直分量；
采用本发明模拟得到的波场快照图如图4(a)～图4(d)所示，从图中可以看出，地震波在目标区域以弹性纵横波形式传播，在其他区域以声学纯纵波形式传播，波场快照证明了本发明方法的正确性。
为了作为对比，将采用本发明模拟得到的地震波正演模拟炮记录与其他地震波正演模拟炮记录进行对比，具体如下：
图5(a)～图5(d)为只考虑了各向异性、双相性的各向异性双相介质地震波正演模拟炮记录，其中图5(a)为固相水平分量，图5(b)为固相垂直分量，图5(c)为流相水平分量，图5(d)为流相垂直分量；
图6(a)～图6(d)为只考虑了黏弹性、双相性的黏弹各向同性双相介质地震波正演模拟炮记录，其中图6(a)为固相水平分量，图6(b)为固相垂直分量，图6(c)为流相水平分量，图6(d)为流相垂直分量；
通过将3(a)～图3(d)与图5(a)～图5(d)对比可以看出，采用本发明模拟得到的地震波正演模拟炮记录相比于各向异性双相介质地震波正演模拟炮记录，地震波的形态基本一致，但存在明显的地震波场衰减，说明本发明除了能准确模拟地震波的各向异性和双相性之外，还能模拟出地下介质中的黏弹衰减性；通过将3(a)～图3(d)与图6(a)～图6(d)进行对比可以看出，采用本发明模拟得到的地震波正演模拟炮记录相比于黏弹各向同性双相介质地震波正演模拟炮记录，振幅能量一致，不同的是，本发明能够将各向异性引起的波形准确地模拟出来。因此可以说明本发明方法同时考虑了地下介质的黏弹、各向异性和双相介质的性质，得到了准确的正演模拟结果。
本发明提出的一种黏弹各向异性双相介质区域变网格求解算子，既同时考虑了地下介质的黏弹性、各向异性和双相性，又只针对目标区域进行弹性纵横波场的模拟；为后续的复杂地质构造成像、反演工作提供地震波场计算基础。
上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。
1.一种黏弹各向异性双相介质区域变网格求解算子，其特征是，包括以下步骤：
步骤1：输入速度场及各种参数场，建立观测系统；
步骤2：选择目标区域进行局部网格加密；
步骤3：在全局区域采用黏声各向异性介质进行波场延拓；
步骤4：在目标区域采用黏弹各向异性双相介质进行波场延拓；
步骤5：在目标区域的边界处进行黏声各向异性和黏弹各向异性双相波场传递；
步骤6：应用人工边界条件对边界反射进行吸收；
步骤7：输出炮记录和波场快照。
2.如权利要求1所述的黏弹各向异性双相介质区域变网格求解算子，其特征是，所述步骤2中，在全局区域进行粗网格剖分，在目标区域进行细网格剖分。
3.如权利要求2所述的黏弹各向异性双相介质区域变网格求解算子，其特征是，所述步骤2中，全局区域的粗网格尺寸为目标区域细网格尺寸的奇数倍。
4.如权利要求3所述的黏弹各向异性双相介质区域变网格求解算子，其特征是，所述步骤3中，在全局区域采用黏声各向异性波动方程(1)进行p波传播；
方程(1)中，p表示各向异性声波波场；q表示各向异性声波辅助波场；vpz表示对称轴方向的相速度；ε表示用以表征纵波各向异性强度的thomson参数，δ表示用以表征动校正速度的thomson参数；t表示时间；fx表示水平分量的震源；fz表示垂直分量的震源；x表示空间水平分量的空间坐标；z表示空间垂直分量的空间坐标；τ表示松弛时间，τ由公式(2)求得：
τ＝τε/τσ-1(2)
公式(2)中，τσ为应力松弛时间，τε为应变松弛时间；其中，τσ和τε与纵波品质因子qp的关系如公式(3)、公式(4)所示；
公式(3)和公式(4)中，ω代表角频率。
5.如权利要求4所述的黏弹各向异性双相介质区域变网格求解算子，其特征是，所述步骤3中全局区域采用的黏声各向异性波动方程(1)由公式(5)递推求得；
公式(5)中，n表示差分计算次数；f和f-1分别为傅里叶变换与傅里叶反变换；δt是时间步长；kx为水平分量的波数；kz为垂直分量的波数；
公式(6)～公式(9)中，δx表示水平分量的空间网格步长；δz表示垂直分量的空间网格步长；c2,0和c2,m为二阶偏导数的有限差分系数；m为积分变量；m为差分阶数；i为水平分量的空间坐标计数，j为垂直分量的空间坐标计数。
6.如权利要求5所述的黏弹各向异性双相介质区域变网格求解算子，其特征是，所述步骤4中，在目标区域采用黏弹各向异性双相介质一阶速度应力方程(10)进行多分量波场延拓；
方程(10)中，ux为固相水平分量的质点速度；uz为固相垂直分量的质点速度；ux为流相水平分量的质点速度；uz为流相垂直分量的质点速度；τxx为水平分量的正应力；τzz为垂直分量的正应力；τxz表示切应力；s表示流相正应力；q为固体与流体之间相耦合弹性系数；r为流相弹性系数；b表示耗散系数；t表示时间；x表示空间水平分量的空间坐标；z表示空间垂直分量的空间坐标；l1表示体积模量的标准线性固体个数；l2表示剪切模量的标准线性固体个数；
方程(10)中，d1、d2和d3由公式(11)求得：
公式(11)中，ρ11表示单位体积中固体相对流体运动时固体部分总的等效质量，ρ22表示单位体积中流体相对固体运动时流体部分总的等效质量，ρ12表示单位体积中流体和固体之间的质量耦合系数；
方程(10)中，
公式(12)中，
方程(10)中，e11、e55由公式(13)求得：
公式(13)中，ε表示用以表征纵波各向异性强度的thomson参数；ρ表示密度；vp表示纵波速度；vs表示横波速度；
方程(10)中，
公式(14)中，λ为与纵横波速度和密度有关的拉梅常数；μ为与横波速度和密度有关的拉梅常数；
7.如权利要求6所述的黏弹各向异性双相介质区域变网格求解算子，其特征是，所述步骤4中，目标区域采用的黏弹各向异性双相介质一阶速度应力方程(10)由公式(16)递推求得；
公式(16)中，δt1为目标区域的时间采样间隔，δt1与时间步长δt之间的关系如公式(17)所示；
δt＝(2k 1)δt1(17)
公式(17)中，k为任一正整数；
公式(16)中，dx表示水平分量的空间一阶偏导数，dz表示垂直分量的空间一阶偏导数，dx和dz由公式(18)求得：
公式(18)中，u为波场值，可表示ux、uz、ux、uz、τxx、τzz、τxz、
8.如权利要求7所述的黏弹各向异性双相介质区域变网格求解算子，其特征是，所述步骤5中，在目标区域的左右边界采用公式(20)所示的边界条件进行声压p与应力τxx的波场传递；
-p＝τxx(20)
在目标区域的上下边界采用公式(21)进行声压p与应力τzz的波场传递；
-p＝τzz(21)。
9.如权利要求8所述的黏弹各向异性双相介质区域变网格求解算子，其特征是，所述步骤5中，目标区域边界处全局区域延拓的黏声各向异性介质声波波场与目标区域延拓的黏弹各向异性双相介质弹性波波场互相传递的步骤如下：
步骤51：将目标区域上下边界上处于粗网格点上的声压p传递给对应的细网格点上的τxx，粗网格点上的uz传递给对应细网格点上的uz；将目标区域左右边界上处于粗网格点上的声压p传递给对应的细网格点上的τzz，粗网格点上的ux传递给对应细网格点上的ux；
步骤52：利用插值公式计算目标区域边界上其他细网格点上的波场值；
步骤53：采用降阶方法求解过渡区域的波场值；
步骤54：在目标区域其他细网格点黏弹各向异性双相介质波动方程(16)进行波场延拓；
步骤55：将目标区域的过渡区域上位于细网格点上的波场值传递给对应的粗网格；
步骤56：将目标区域上下边界上位于粗网格位置上的细网格点波场值τxx传递给粗网格点波场值p，将目标区域左右边界上位于粗网格位置上的细网格点波场值τzz传递给粗网格点波场值p，位于粗网格位置上的细网格点波场值ux和uz传递给粗网格点波场值ux和uz。
10.如权利要求9所述的黏弹各向异性双相介质区域变网格求解算子，其特征是，所述步骤6中，在全局区域采用pml边界条件对波场的衰减进行吸收，对目标区域的边界采用cfs-pml边界条件进行波场延拓；
使用拉伸函数将速度-应力方程变换到拉伸坐标域；在频率域中，2d拉伸函数的具体表述形式：
式中，σx代表cfs-pml中的水平分量衰减系数；σz代表cfs-pml中的垂直分量衰减系数；εx为水平分量的拉伸函数；εz为垂直分量的拉伸函数；ω为角频率；i表示复数虚单位；
笛卡尔坐标系(x，z)与拉伸坐标系
将方程(24)代入方程(10)，得到拉伸坐标系下一阶速度-应力方程的表达式：
其中，
