0 引言

近年来，分数阶微积分被广泛应用于反常扩散、信号处理与控制、流体力学、图像处理、软物质研究、地震分析、黏弹阻尼器等领域.相较于整数阶导数，分数阶微分算子可以更简洁准确地描述具有历史依赖性和空间全域相关性的复杂力学和物理过程.分数阶导数在数学上的一个主要特点是它的非局部性，即当前状态与过去所有状态有关.

国内外大多数学者通过组合机械元件来描述黏弹性模型，如Kelvin-Voigt模型(聂建新等，2012)、Maxwell模型等.但是，将此类模型用于黏弹地震波衰减模拟方面，会存在对高频成分的衰减过量与不足的问题，并不能真实地表达黏弹介质的物理和力学特征(孙成禹和印兴耀，2007；孙成禹等，2010).品质因子Q能够度量介质对地震波能量的吸收衰减，实验观测数据表明，在地震勘探频带范围内，固体的品质因子基本不随波的频率变化，是固体介质的固有属性，但吸收系数与频率成正比(Aki and Richards, 1980).鉴于此，许多研究者将多个标准线性固体并联得到的广义标准线性固体(Zener模型)，在给定的频带内拟合黏弹介质中Q值为常数的特征，这种广义标准线性固体被广泛用来构建近似常Q黏弹模型，并使用一系列方法开展了数值模拟(Liu et al., 1976；Emmerich and Korn, 1987；Carcione et al., 1988；Zhu et al., 2013；郭振波等，2014)，但出于对计算效率的考虑，这些广义线性体模型在实际地震正演模拟、偏移成像等方面中使用相对较少.郭智奇等(2008)用Carcione记忆变量自身迭代方法来代替应力应变之间的卷积运算，较卷积运算计算量与计算复杂度确实减少许多，精度也较高，但是对于较为复杂的方程在大计算域、长时间历程情况下计算量和存储量也是相当大的.除这些近似常Q模型与记忆变量自身迭代方法之外，Kjartansson建立了一种基于数学模式假定的完全常Q模型，可对非完全弹性介质吸收频散关系予以定量描述，基于Kjartansson理论(Kjartansson，1979)，Carcione(Carcione et al., 2002；Carcione, 2009)使用Grünwald-Letnikov近似法进行了常Q黏声波和黏弹波的数值模拟，Qiao等(2016)基于Kjartansson常Q理论采用分数阶导数对黏弹各向异性单相介质中波场进行了数值模拟.

波在含流体孔隙介质中传播时，宏观Biot流动机制和微观Squirt流动机制作为一个耦合过程共同对地震波的传播和衰减产生影响(Biot, 1941, 1956a, b).Dvorkin和Nur(1993)从孔隙各向同性一维问题出发，提出统一处理两种衰减机制的Biot-Squirt(BISQ)模型，兼顾了流体在轴向的自由流动(体现Biot的全局流的流动)及在径向的流动(体现Squirt的局部流的流动)，很好地解释了实验中观测到的高频散和强衰减现象.Parra(1997)从频率域孔隙弹性波方程出发，将一维各向同性的BISQ模型推广到横向各向同性孔隙介质情况，Yang和Zhang(2002)将其推广到高维横向各向同性和一般各向异性情况.此后，国内外众多学者又对BISQ模型的基本理论和应用范围进一步完善和发展(Parra，2000；杨宽德等，2002；聂建新等，2004；李红星等，2007；Yang et al., 2007; 聂建新等，2010；杨宽德等，2011).但是，BISQ模型中含有一个表征喷射流动的微观尺度参数—特征喷流长度，含虚数单位，其物理意义不明确、实验测得相对困难，不利于时间域正演计算.Diallo和Appel(2000)提出了一种不依赖于特征喷流长度参数的改进BISQ模型，改进的BISQ模型克服了BISQ模型中特征喷流长度的缺陷，使该模型具有较大的实用前景.虽然众多国内外学者对BISQ模型进行了诸多改进，但经典的Biot理论和改进后的BISQ理论均假设双相介质的固体骨架为理想弹性体，未考虑真实地下介质固体骨架的黏弹性(聂建新等，2012).

另一方面，改进的BISQ理论也未考虑实际地下介质中孔隙流体的流变性，Liu和Katsube(1990)指出由于孔隙流体的黏滞性，孔隙介质中固、液相之间既存在相对平动，又存在相对转动，因此在含黏滞流体孔隙介质的固液分界面会出现慢横波，它的存在也是孔隙介质中地震波传播能量耗散的原因之一, 故在进行双相介质中地震波场模拟时，孔隙黏滞流体的流变性也应该被考虑进来.

鉴于双相介质理论的相对不完善，以往研究中所建立的模型的适用范围将受到一定限制.本文旨在将分数阶微分算子的记忆特性引入双相介质理论来描述固体骨架的黏弹性、并综合考虑孔隙黏滞流体的流变性，固体骨架各向异性.针对双相孔隙储层岩石，提出一种适用范围更为广泛的分数阶黏弹双相孔隙介质模型，并对其进行波场数值模拟与特征分析.为此，本文基于含分数阶时间导数的常Q黏弹本构关系，将改进的BISQ模型(Diallo and Appel, 2000)、分数阶常Q黏弹双相介质理论与孔隙黏滞流体流变性(Liu and Katsube, 1990)结合，提出了一种新的基于分数阶时间导数黏弹骨架本构关系的含黏滞流体双相VTI介质模型，并推导出其波传播方程.该方程各物理参量均具有明确物理意义，然后采用精度较高、实现快速、占用内存较小的高阶交错网格有限差分法求解该方程.另外，因分数阶微分算子具有全局记忆特性，对长时间历程及大计算域分数阶常Q黏弹波动方程数值模拟时计算量和存储量极大.为此，对于分数阶时间导数采用短时记忆有限差分法离散，结合上述的高阶交错网格有限差分法和短时记忆有限差分法，我们实现了该类介质的地震波场数值模拟，并且具体分析了不同黏滞相界、不同品质因子组及双层介质模型情况下波的传播特征，与此同时研究了波在介质界面处的反射、透射及波形转换现象.最后，对地震波在基于分数阶时间导数常Q黏弹本构关系的含黏滞流体双相VTI介质中的传播特征进行了总结.

1 常Q分数阶波传播方程 1.1 建立模型的假设前提

(1) 孔隙只含有一种各向同性、黏滞性和可压缩的流体，且为饱和态.

(2) 孔隙中流体可以在波传播方向和垂直传播方向发生流动，当波传播时，流体能可以从软孔隙中喷出，泵入刚性孔隙中，流体渗流遵守广义达西定律.孔隙流体与岩石基质之间无化学作用，且不考虑热弹性效应.

(3) 介质由于岩性、孔隙性、构造应力等因素而具有各向异性.

1.2 常Q分数阶黏骨架各向异性本构关系

地下介质固体骨架的线性黏弹性，在应力-应变的本构关系上，可表示为松弛函数与应变对时间导数的卷积形式(Christensen，1982)，公式为

(1)

其中，i、j、k、l=1, 2, 3，*表示卷积运算，C_ijkl为四阶松弛张量，σ_ij和ε_kl分别为应力张量和应变张量，对于VTI各向异性介质有：

(2)

其中，c₁₂⁰(t)=c₁₁⁰(t)-2c₆₆⁰(t)，根据Qiao等(2016)c₁₁⁰(t)与c₃₃⁰(t)是两个独立的松弛函数，分别对应P波模量的水平分量和垂直分量，即：

其中，λ_//和λ_⊥、μ_//和μ_⊥分别表示平行于地层和垂直于地层的VTI介质拉梅常数，c₁₃⁰(t)=λ_⊥(t)，剪切模量c₄₄⁰(t)=μ_*(t)，c₆₆⁰(t)=μ_//(t)，其中μ_*(t)是连接两个方向的剪切模量.

将(2)式代入(1)式，在二维x-z平面，得：

(3)

其中：

H(t)为Heaviside阶跃函数，ω₀=1/t₀为参考频率(Carcione，2002).

根据分数阶导数理论(Caputo et al., 1971)，式中卷积运算转化为分数阶时间导数，(3)式变为

(4)

这里D^2γ=∂^2γ/∂t^2γ是2γ阶时间分数阶微分算子，且有：

其分数阶导数阶数定义为

这里Q_P和Q_S分别是P波和S波的品质因子，根据Zhu和Tsvankin(2006)有：

另外，参考模量M_P//、M₁₃、M_P⊥、M_μ可分别由参考频率为ω₀时的平行于地层的P波速度v_P//、垂直于地层的P波速度v_P⊥和S波速度v_S表示为

其中，ρ_s为固体骨架密度，C₁₁⁰、C₁₂⁰、C₁₃⁰、C₃₃⁰、C₄₄⁰、C₆₆⁰为干燥情况下固体骨架的弹性刚度矩阵张量C⁰中的元素.

1.3 二维孔隙黏滞流体压力表达式

地震波在含黏滞流体孔隙介质中传播时，孔隙或裂缝刚性程度的不同会使孔隙流体产生压力梯度，这样地震波就诱导孔隙流体流动使压力重新达到平衡.在双相孔隙介质理论中，孔隙黏滞流体压力表达式的确定是建立波传播方程的前提.当不考虑孔隙黏滞流体流变性时，根据刘财等(2013)提出的方法，将基于改进BISQ机制的HTI三维波动孔隙流体压力表达式变为二维VTI介质中二维波动孔隙流体压力，公式为

(5)

其中第一项表示与x方向波动有关的Biot流动和Squirt流动对总孔流体压力的贡献，第二项表示与z方向波动有关的两种流体流动对总孔隙流体压力的贡献.其中，W=ϕ(U-u)表示相对固、流两相位移矢量，U=(U₁，U₃)^T和u=(u₁，u₃)^T分别代表流相和固相位移矢量，而u_i和U_i(i=1，3)则分别代表固相和流相质点在笛卡尔坐标系两个正交方向x和z上的位移分量；ϕ为孔隙度；β为双相介质压缩率系数；γ_i(t)=α_ii(t)/3+(i=1，2，3)，其中，α_ii为有效应力之孔隙黏弹性系数张量α(t)=diag(α₁₁(t)，α₂₂(t)，α₃₃(t))的张量元素.对于双相VTI介质来说，有效应力的孔隙黏弹性系数与双相黏弹介质松弛函数之间的关系为

将W和γ_i(t)(i=1，3)的表达式代入(5)式，并整理得：

(6)

式(6)即为以固相和流相位移分量为基本未知量的基于改进BISQ机制的二维双相VTI介质中的孔隙黏滞流体压力表达式.

流体微元应力张量S与孔隙流体压力p′成正比，满足关系：

(7)

如果储层含黏滞流体，例如：稠油，在波场模拟时，应考虑黏滞流体流变性.波在含黏滞流体孔隙介质中传播会产生流体振动速度梯度，进而产生剪切应变，导致液相中既存在剪切应力又存在附加法相应力，由牛顿内摩擦定律可知，黏性摩擦力产生的应力偏量T_ij为

(8)

其中，η为流体黏滞系数；V_i(V_j)为流相质点在正交方向x、z上的速度分量；下标i、j=(1, 3)分别表示x、z方向分量，文中后续公式中i、j意义同此(卢明辉等，2009).

由式(7)和式(8)可得含黏滞性流体孔隙介质液相的流变学本构关系为

(9)

令p为考虑黏滞流体流变性时孔隙黏滞流体压力，则有：

(10)

推得考虑流体流变性时基于改进BISQ机制的孔隙黏滞流体压力为

(11)

1.4 波传播方程

与Biot模型、BISQ模型和改进BISQ模型不同的是，本文采用基于Kjartansson常Q理论的分数阶时间导数黏弹本构关系来描述双相介质固体骨架的黏弹性，并引入流变学本构关系表征孔隙流体的黏滞力学行为.因此，从Yang和Zhang(2002)提出的含BISQ机制及固流耦合各向异性的孔隙弹性波动理论出发，给出基于分数阶时间导数常Q黏弹本构关系的含黏滞流体双相VTI介质的黏弹骨架本构关系、几何方程、运动方程以及二维波传播方程.

黏骨架本构方程为

(12)

其中，σ为总应力张量，e(t)为固体骨架应变张量，孔隙黏滞流体压力p由(11)式给出.

几何方程为

(13)

运动方程为

(14a)

(14b)

其中σ_ij表示总应力张量分量；ρ₁=(1-ϕ)ρ_s，ρ₂=ϕρ_f，ρ_f为孔隙流体密度；ρ_22i=ρ₂-ρ_12i，ρ_12i=-ρ_ai，ρ_ai为x_i方向的固-流耦合附加质量密度；流体阻抗系数r_ij=(k_ij)^-1, k_ij为渗透率张量k的元素，在VTI介质中k=diag(k₁₁, k₁₁, k₃₃).

综合式(11)、(12)、(13)、(14)即可推导出分数阶时间导数常Q黏弹本构关系的含黏滞流体双相VTI介质波传播方程，基本未知量为固相和流相位移u_i和U_i.在二维x-z平面内波传播方程是由以下5个方程构成的含卷积运算的偏微分方程组，公式为

(15a)

(15b)

(15c)

(15d)

(15e)

根据方程（4），对方程（15）进行进一步推导，写成速度应力式为

(16a)

(16b)

(16c)

(16d)

(16e)

(16f)

(16g)

(16h)

其中，v_x、v_z分别代表x和z固相速度分量.式中其他变量为

2 有限差分数值解法 2.1 高阶交错网格有限差分格式

有限差分法(FDM)作为一种常用的地震波数值模拟方法，最早被用于地震波场数值模拟，相较于伪谱法(PSM)，实现简单便捷，占用内存小，可以较好地适应剧烈变化的复杂地下介质.交错网格有限差分方法由于将空间变量的导数在响应的网格点之间的半程上计算，可以在不增加计算量和内存的前提下提高精度，有效地减少传统有限差分法的网格频散.

2.2 短时记忆长度

分数阶导数波动方程的数值模拟需要极大的计算量和存储量，尤其对长时间或大计算域的模拟.为了减少计算量、提高计算效率，Podlubny(1999)提出了“短时记忆算法”，通过截断算子长度来达到减少运算量和存储量的目的，仅考虑当前时刻t以前且对当前时刻影响较大的有限时间区间[t-L, t]，L称为记忆长度.

对于时间分数阶导数，有多种定义，对于本文算例由于都是同样实数集上的算子，理论上它们是等价的，但是从数值计算角度考虑，本文采用Grünwald-Letnikov定义，公式为

(17)

式中，Δt为时间步长，定义：

这样τ(2γ, m)有递推关系：

当m=0时，τ(2γ, 0)=1.使用短时记忆中心差分近似(17)式，当Δt→0时，在t=mΔt时离散化有：

(18)

其中，C_n为相应的高阶交错网格差分算子系数.

2.3 稳定性条件和吸收边界条件

显式有限差分数值模拟，需要考虑计算过程的稳定性，由于分数阶偏微分方程数值计算的稳定性较为复杂，至今缺乏系统的研究.本文通过数值试算得出本文所建立分数阶偏微分方程的稳定性条件，其近似表达式为

(19)

式中Δx和Δz为横向和纵向空间网格步长，Δt为时间间隔，δ=0.872为稳定性常数，v=max{v_P//, v_P⊥, v_S}.

为了消除人工边界的虚反射，本文采用简便且实用的Cerjan衰减边界(Cerjan et al., 1985)，即沿人工边界向外扩充N个网格点，使进入扩充区域内的波场通过乘以下面的衰减因子，公式为

(20)

逐渐衰减为零，其中i为吸收层内节点序号，a为衰减系数，通过多次试验确定其最佳值.

3 数值模拟与分析

为了考察本文所建模型的有效性，我们选取不同的模型进行波场数值模拟研究.

3.1 模型1：不同相界单层模型

为了详细分析含黏滞流体常Q双相VTI介质中各种弹性波的传播特性以及流体黏滞性系数对波传播特性的影响，设计理想相界单层模型及流体黏滞性不同的单层模型.模型大小为3000 m×3000 m，介质物性参数分别为：干燥情况下各向同性背景孔隙介质的v_P//=5600 m·s^-1、v_P⊥=5167 m·s^-1、v_S=3126 m·s^-1、ϕ=0.2、ρ_s=2500 kg·m^-3、ρ_f=1040 kg·m^-3，固相体积模量K_s=80 GPa，x和z方向的固流耦合密度ρ_a1=420 kg·m^-3、ρ_a3=450 kg·m^-3, x和z方向渗透率k₁₁=6×10^-10 m²·s^-1、k₃₃=1×10^-9 m²·s^-1，近似理想相界流体黏滞系数为：η=1×10^-10 Pa·s，四种不同的流体黏滞系数分别为η₁=1×10^-10 Pa·s、η₂=0.001 Pa·s、η₃=0.004 Pa·s、η₄=0.01 Pa·s，Q_P=60，Q_S=30，ω₀=80 Hz，短时记忆长度L=0.15 s.其他计算参数选择为：空间采样间隔Δx=Δz=10 m，时间步长为Δt=5×10^-4s，t=0.2 s.震源子波采用Ricker子波，震源位置(1500 m，1500 m)，主频为40 Hz，加载方式为x方向集中力源.

图 1给出了近似理想相界情况下含黏滞流体双相VTI介质中固相和流相速度在t=0.15 s时刻的波场快照.分析近似理想相界时波场数值模拟结果：(1)从图 1可以看出，波场快照中均包含3种不同的波，它们由外向内依次为：1快准P波、2准SV波和3慢准P波，固相和流相波场快照中慢准P波均清晰可见，且能量较强，这与无耗散方程所描述的慢P波传播特征相同；(2)对比固相波场与流相波场可见，这三种波引起的固相和流相质点振动相位相同，但快准P波及准SV波引起的固相质点振动略强于流相质点，慢准P波引起的流相质点振动明显强于固相质点，这些现象与Biot模型中质点振动特性相同(刘洋和李承楚，2000)；(3)固相和流相x、z分量中的快准P波、准SV波和慢准P波波前面均呈椭圆形，准SV波的波前面出现波面尖角现象；故在含黏滞流体常Q黏弹双相VTI介质中仍存在横波分裂，体现了黏弹双相各向异性介质中固体骨架各向异性所引起的波场特征.在图中并未观测到慢横波，是由于在地震频带内，慢横波的速度很低，振幅很弱，对实际地震勘探数据处理的影响很小，但它们的存在以及对其传播特征的研究，将进一步丰富和完善双相介质黏弹性波传播理论.

图 3给出了不同流体黏滞系数情况下固相速度水平分量在t=0.15 s时波场快照.从图 3波场快照中可以看出：随着流体黏滞性系数增大，快准P波振幅水平方向两端衰减明显，而垂直方向两端增强明显；准SV波尖角现象越来越不明显，波前面呈椭圆型；慢准P波振幅衰减.为了更明显显示不同流体黏滞性系数对含黏滞流体常Q黏弹双相VTI介质中地震波传播特性的影响，图 3给出了网格点坐标为(150，200)处不同流体黏滞性系数固相质点垂直速度分量的时间记录.从图 3中可以看出：随着流体黏滞性系数变大，快准P波的相位延迟，准SV波的相位提前，振幅逐渐均变小，但快准P波振幅减小程度比准SV波振幅减小程度明显；快准P波的波长变窄，而准SV波波长几乎无变化.这是由于孔隙流体黏滞性变强，在压力梯度下在孔隙及裂缝间流动变缓，起到阻碍纵波传播作用，而随着孔隙流体黏滞性变强，使得质点间剪切方向联系作用变强，进而横波相位提前.

3.2 模型2：不同品质因子组黏滞相界单层模型

为了考察黏滞相界模型中不同品质因子情况下波的传播特征，取流体黏滞系数η=0.001 Pa·s, 模型大小、震源子波类型，震源位置、震源主频、空间采样间隔、时间步长、模拟时间和加载方式均同模型1.波场数值模拟选择四组不同的品质因子参数，第一组：Q_P=60，Q_S=30；第二组：Q_P=30，Q_S=15；第三组：Q_P=10，Q_S=5；第四组：Q_P=6，Q_S=3，其他介质参数与模型1中的参数一致.

图 4分别给出了四组不同同品质因子组情况下黏滞相界模型中固相速度水平分量在t=0.15 s时刻的波场快照.从图 4中可以看出，随着品质因子的减小，波形无明显变化，振幅衰减逐渐减小.同样为了更明显显示不同品质因子组对含黏滞流体常Q黏弹双相VTI介质中地震波传播特性的影响，图 5给出了网格点坐标为(150, 200)处不同品质因子组固相质点垂直速度分量的时间记录.从图中可以看出随着品质因子的减小，快准P波的左侧波峰相位提前，右侧波峰相延迟，导致波长变宽，波谷相位延迟；准SV波的左侧波谷相位提前，右侧波谷相位延迟，导致波长变宽，波峰相位提前；二者振幅均明显减小.这是由于孔隙介质固体骨架品质因子减小，表现为黏弹性增强.随着固体骨架的黏弹性增强，受力之后弛豫现象增强，因此，波形变宽，振幅减小.

3.3 模型3：双层介质模型

为了分析含黏滞流体常Q黏弹双相VTI介质中波在介质分界面处的反射、透射及波形转换现象，设计了一个大小为3000 m×3000 m双层介质模型，界面深度为1700 m，ρ_f=1040 kg·m^-3，η=0.001 Pa·s.模型上层介质参数为：干燥情况下各向同性背景孔隙介质的v_P//=3700 m·s^-1、v_P⊥=3400 m·s^-1、v_S=1600 m·s^-1，Q_P=30，Q_S=15，ρ_S=2650 kg·m^-3；模型下层介质参数为：干燥情况下各向同性背景孔隙介质的v_P//=5400 m·s^-1、v_P⊥=5060 m·s^-1、v_S=3000 m·s^-1，Q_P=5，Q_S=3，其他模型参数及波场数值模拟计算的震源参数和加载方式，空间采样间隔、时间步长和模拟时间均同模型1.选取固、液相x分量t=0.3 s时的波场快照及地震记录来分析3种波在介质分界面处的反射、透射及转换特征，见图 6.

分析图 6可知：快准P波、准SV波和慢准P波在介质分界面处均发生了反射、透射及波形转换现象；在波场快照中可同时观测到多种波，这些波大致可以分为5类，参见(刘财等，2013).除此之外，由于固相上下层的品质因子差异，在固相波场中可以明显观测到每种波的上下层振幅差异，但是在液相波场中可见固相上下层品质因子差异对慢P波、快准P波产生的透射转换准SV波及慢准P波产生的透射转换快准P波几乎无影响，观测不到振幅差异.这是由于在震源附近固体骨架黏弹性对流相影响较小而对固相影响较大，当波继续传播，固相的衰减也将影响流相波场衰减，因此流相下层离震源较远处波场也呈衰减状态.双层含黏滞流体常Q黏弹双相VTI介质模型中的波场信息比较丰富，每种波的入射都将在分界面处产生多种波，并且波的反射、透射及转换都是在上下层不同品质因子组的情况下发生的.深入研究含黏滞流体常Q黏弹双相介质模型中地震波在分界面处的反射、透射及转换现象，分析各种波的传播特征，能够加深对实际地下黏弹性双相介质中波传播机理的认识，进而有助于利用多波地震资料对含黏滞流体黏弹双相介质油气藏进行更为深入的研究，具有重要的理论和实际意义.

4 结论

将固体骨架分数阶时间导数常Q黏弹各向异性理论、孔隙流体黏滞性理论、改进的BISQ模型及双相介质理论相结合，建立了一种更具有普适性的含黏滞流体分数阶时间导数常Q黏弹骨架孔隙介质模型，即基于分数阶时间导数常Q黏弹本构关系的含黏滞流体双相VTI介质模型，推导了介质的二维速度-应力波传播方程，使用短时记忆方法的中心差分计算分数阶时间导数，通过对模拟结果进行分析，得出以下结论：

(1) 分数阶导数考虑黏弹介质的全局记忆特性，使用基于短时记忆方法的中心差分近似分数阶时间导数，在保证一定精度的前提下提高了计算效率，解决了分数阶导数计算耗时长和储存量大的问题，波场模拟结果表明以上方法可以准确模拟常Q黏弹各向异性介质中地震波传播过程，模拟精度较高.

(2) 随着流体黏滞性系数增大，孔隙流体黏滞性变强，在压力梯度下在孔隙及裂缝间流动变缓，起到阻碍纵波传播作用，而随着孔隙流体黏滞性变强，使得质点间联系作用变强，进而横波相位提前.

(3) 随着品质因子的减小，孔隙介质固体骨架品质因子表现为黏弹性增强.随着固体骨架的黏弹性增强，受力之后弛豫现象增强，因此，波形变宽，振幅减小.

(4) 在双层介质中，上下层介质孔隙流体黏滞性系数相同，固体骨架品质因子组不同的情况下，快准P波、准SV波和慢准P波在双层介质的分界面处均发生反射、透射和波型转换现象，波场信息丰富，在固相波场中可以明显观测到每种波的上下层振幅差异，但是在液相波场中可见固相上下层品质因子差异对慢P波、快准P波产生的透射转换准SV波及慢准P波产生的透射转换快准P波几乎无影响，观测不到振幅差异.这是由于在震源附近固体骨架黏弹性对流相影响较小而对固相影响较大，当波继续传播，固相的衰减也将影响流相波场衰减.

本文提出了基于分数阶时间导数常Q黏弹本构关系的含黏滞流体双相VTI介质模型，并进行了波场模拟，波场数值模拟结果证实了其合理性和有效性，为与其他衰减模型以及实际黏弹各向异性介质地震波的传播机理的对比提供了参考，为从含黏滞流体常Q黏弹双相各向异性理论角度研究储层的地震响应奠定了理论基础，同时也为进一步研究黏弹双相介质中黏弹地震波反演提供正演理论基础.