2015, Vol. 30
波诱导流体流动是含流体孔隙介质中波衰减和频散的主要机制.研究表明，流体Biot流动(全局流动)和局部喷射流动是波诱导流体流动的两种重要形式，对含流体孔隙介质中弹性波的传播有着重要影响.Biot(1956，1962)考虑波诱导的流体全局流动，最先建立了流体饱和多孔隙介质中波传播的双相介质理论，奠定了孔隙介质理论研究的基础.然而，在很多情况下，Biot理论预测的波能量衰减和速度频散比实际观测值低，不能对许多实验结果和实际数据给出合理地解释.Dvorkin和Nur(1993)认为当波在孔隙介质中传播时Biot流动和喷射流动同时存在，它们作为一个耦合过程共同对波的衰减和频散产生影响，他们基于孔隙各向同性一维问题，将这两种固-流相互作用的力学机制有机地结合起来，提出了统一的Biot-Squirt(BISQ)模型，并将BISQ模型的预测结果、Biot理论的预测结果与实验观测数据进行比较，发现BISQ模型的预测结果更为准确.之后，国内外一些学者对BISQ模型进行了广泛深入地研究，将该模型由各向同性推广到各向异性(Parra，1997;Yang and Zhang， 2002)，由弹性推广到粘弹性(Cheng， et al.， 2002；张智等，2005；Nie et al., 2008)，由流体饱和推广到流体非饱和(Nie et al., 2004)，研究了基于BISQ模型的孔隙介质中的波衰减和频散(Yang and Zhang， 2000；李红星等，2009；聂建新等，2010；张显文等，2010)、波场数值模拟(杨宽德等， 2002a，2002b，2011；刘财等，2007；李红星等，2007)、储层参数反演(聂建新等，2004；高凤霞等，2013)等问题，形成了较为系统和完善的BISQ理论体系.
地震波场数值模拟技术是认识地震波在地下介质中的传播规律，建立地震响应与储层参数之间关系的重要手段.数值模拟基于BISQ模型的双相介质中的地震波场是认识含流体孔隙介质中地震波传播规律的重要方法.杨宽德等(2002a，2011)基于BISQ模型先后适用FCT有限差分法和FCT紧致有限差分法模拟了Biot流动和喷射流动作用下的弹性波传播；杨宽德等(2002b)、李红星等(2007)分别采用低阶交错网格有限差分法和高阶交错网格有限差分法对基于BISQ模型的双相横向各向同性介质中的波场进行了数值模拟；Wang等(2008)研究了基于BISQ机制的三维正交各向异性介质中弹性波的波场数值模拟问题.然而，以上波场数值模拟工作都是在时间域进行的(刘财等，2005；郭智奇等，2007；皮红梅等，2009).
在Biot流动和喷射流动两种机制耦合作用下，双相介质中的波动方程是一个复变系数偏微分方程组(Yang and Zhang， 2002)，方程的部分系数是与频率相关的复数值.实际地震波是由多个单一频率的简谐波合成的复合波，因而，当地震波在基于BISQ模型的双相介质中传播时，不同频率的谐波具有不同的响应特征，导致地震波的传播特性具有频率相关性.因而，若要通过数值求解波动方程，实现基于BISQ模型的双相介质中弹性波传播的数值模拟，需要采用频率域波场正演模拟的方法.时间域的波场正演方法仅适用于波动方程的系数不随频率变化的情况，即仅仅适用于低、高频极限情况，而对中间频段，波动方程系数是与频率相关的，因而就必须在频率域对于每个频率分别求解波动方程，以获得每一频率下的地震波场相应特征.
频率域正演模拟方法是按照频率片对空间网格进行整体求解方程组，计算误差分配到每一个网格点上，并且各个频率片之间是独立计算的，不存在累积误差，是一种特别适合于并行计算的正演方法.采用频率域正演模拟方法进行孔隙介质中的波传播数值模拟相比于时间域方法而言，有以下几点优势：
(1)频率-空间域方程不存在刚性问题，无需将传播方程分裂为刚性部分和非刚性部分；
(2)可以精确地进行全频段的波场数值模拟；
(3)易于求解更为复杂的孔隙介质理论的波动方程(例如，基于孔隙粘弹性理论的波动方程)；
(4)可以为孔隙介质中弹性波频率域全波形反演问题的研究奠定基础.
本文采用Min等(2000)的25点频率-空间域优化差分算子，在保证精度的前提下，该算法将单个横波波长范围内的网格点数减小到3.3个，正演模拟精度高，压制网格频散效果好.本文将采用25点频率-空间域有限差分法数值求解频率-空间域基于BISQ模型的双相各向同性介质的弹性波传播方程，从而研究地震频段内基于BISQ模型的双相各向同性介质中的地震波传播规律.
1 基于BISQ模型的双相介质的波传播方程
基于BISQ模型的双相各向同性介质的弹性波方程，在二维二分量情况下可简化成如下形式(杨宽德等，2002a)：
，为有效应力之孔隙弹性系数张量(Ks为固体体积模量)，P 为流体压力，ρ1= 1－φρs(φ为孔隙度，ρs为固体骨架密度)，ρ2=φρf(ρf为流体密度)，ρ12=－ρa(ρa为固-流耦合附加密度)，ρ22=φρf+ρa，η为流体黏滞系数，k为渗透率.根据Parra(1997)、Yang和Zhang(2002)的研究可知，方程组(1)中的流体压力P为
确定)， 
，ω为角频率).事实上，基于Biot模型的双相介质的波传播方程也可写成(1)式的形式，只是流体压力 P 的表达式中不包含喷射流系数，即相当于喷射流系数S(ω)=1，而这恰好是喷射流系数在高频极限(ω→∞)时的取值，换句话讲，Biot模型可以看成BISQ模型中喷射流动处于非松弛状态的情况，而在其它频率，从喷射流系数S(ω)的计算公式可以看出，它是一个与频率相关的复数值.由此可见，基于Biot模型的双相各向异性介质波传播方程是一个实常系数偏微分方程组，而基于BISQ模型的双相各向异性介质波传播方程则是一个复变系数偏微分方程组.2 频率-空间域有限差分数值解法
2.1 频率-空间域
将(2)式带入方程组(1)式，并进行傅立叶变换，从而得到BISQ模型的双相各向同性介质的频率-空间域波动方程：
和
是频率域位移波场值，G(ω)是频率域震源项.2.2 频率-空间域优化差分算子
优化25点差分格式的构造主要包括三个方面：加权平均算子、平均加速度项和优化系数.在频率-空间域波动方程(4.8)式中，包含了3种12个空间二阶微分项：
(i=1、2、3、4，A1=ux、A2=uz、A3=Ux、A4=Uz).相应的空间差分算子表示为
频率-空间域波动方程中还有多个加速度项，为了逼近加速度项要用到25个网格点，25点加速度项的差分算子如下：
图 1展示了4种25点差分算子网格点的分布情况(Min et al., 2000)，网格中心点的坐标为(i，j).图 1a展示了点(i，j)处二阶导数项∂2/∂x2的25 点差分算子模型.首先以每一行的中心点为基础建构两个二阶中心差分，网格间距分别为Δx和2Δx，权系数分别为c和d；然后再根据该行与点(i，j)距离的不同，对各行进行加权，权系数分别为b1、b2、b3.∂2/∂z2，∂2/∂x∂z求导与此类似，差分系数分布如图 1b和图 1c所示.
 | 图 1 25点有限差分格式的网格点分布情况 ；(d)加速度项. Fig. 1 Electromagnetic25 grid point distribution of finite difference scheme |
差分算子中的Δx和Δz是空间网格步长，ai、bi以及c、d、e、f是25点频率-空间域正演算法的加权系数，为了得到上述差分算子中的优化加权系数，构造其目标函数，使离散模型和连续模型的相速度基本相等.这是一个非线性优化问题，由Gauss-Newton优化方法求解(Min et al., 2000)，由此，可以寻找到优化系数使得差分方程相速度与波动方程相速度尽可能地接近，从而得到25点频率-空间域正演算法的最优化加权系数.
将式(4)至(7)式带入频率-空间域波动方程组(3)式，并将式中的其他变量离散化，则可得到离散化的基于BISQ模型双相各向同性介质的频率-空间域波动方程，以(3a)式为例，其离散化方程见附录A.
相应的(3b、3c、3d)都有类似的采用25点优化差分算子的差分方程，它们构成放置点(i，j)处的差分方程组，且对于网格中每一个点都能建立一个这样的方程组.通过求解该差分方程组，可解得频率ω下的每一网格点的波场值，且对每一个频率ω都能解得所有网格点的波场值.然后对每一个网格点所有频率下的波场值进行IFFT变换，得到该网格点在各个时间点上的波场值，而每个时间点上所有网格的波场值组成波场的时间切片.最后对时间切片抽取数据，得到合成记录.
2.3 PML吸收边界条件
边界条件是有限差分算法的重要组成部分.本文采用的PML吸收边界条件，是由Berenger(1994)最早提出，吴国忱等(2007)、刘璐等(2013)成功地将其应用到频率-空间域有限差分算法中，核心思想是在计算区域四周加上完全匹配层起到吸收边界的作用.我们设定在匹配层中的吸收因子形式为
其中，αj=2πa0f0(lj/LPML)2称为衰减系数，a0为常数，f0为子波主频，lj为波在匹配层中的传播距离，LPML为完全匹配层厚度.a0和LPML需因模型不同而变，以吸收效果最佳为准则，本文根据模型实验效果，取a0=1.79、LPML=20. 3 数值模拟算例及分析
3.1 均匀双相各向同性介质中的地震波传播
为了详细分析地震波在基于BISQ模型的双相各向同性介质中的传播规律，设计一个单层均匀双相各向同性介质模型，模型大小为2000 m×2000 m，介质物性参数分别为：干燥情况下固体骨架的拉梅常数为λ=6.26 GPa，μ=14.61 GPa，孔隙度φ=0.1，固相密度ρs=2650 kg/m3，固相体积模量Ks=38 GPa，渗透率k=0.5×10-15 m2，特征喷射流动长度R=0.002，流体密度ρf=1000 kg/m3，流体声速v0=1500 m/s，固流耦合密度ρa=420 kg/m3.其他计算参数选择为：空间采样间隔Δx=Δz=10 m，震源子波采用Ricker子波，震源位置(1000m，1000m)，主频为25 Hz，模拟中震源为固相x方向集中力源.
基于BISQ模型的双相介质的地震波场数值模拟，考虑理想相界(η=0 Pa·s)和粘滞相界(η=0.001 Pa·s)两种情况，我们分别研究这两种情况下的地震波传播特征.同时，为了突出Biot流动和喷射流动耦合作用对地震波传播的影响，我们将基于BISQ模型的波场数值模拟结果与基于Biot模型的波场数值模拟结果进行对比.在前文中我们已经阐述了Biot模型可以看成BISQ模型的高频极限，因而可以直接将基于BISQ模型的双相介质波传播方程中的喷射流系数Sj(ω)=1，而其它参数不变，采用同样的求解方式实现基于Biot模型的双相介质中弹性波传播的数值模拟.
图 2-5给出了理想相界(η=0)情况下，基于BISQ模型和基于Biot模型的双相各向同性介质地震波场数值模拟结果.图 2显示的是基于BISQ模型的双相各向同性介质中固相和流相垂直位移分量和水平位移分量25 Hz的单频波场，图 3给出了相应的300 ms时时间域波场传播快照；而图 4显示的是基于Biot模型的双相各向同性介质中固相和流相垂直位移分量和水平位移分量25 Hz的单频波场，图 5是相应的300 ms时时间域波场传播快照.从图 2和4中可以看出，单频波场中波前面呈圆形从震源点向外传播，同时，时间域的波场快照中的各种波型的波前面也都呈以震源为中心的圆形，这充分体现了介质的各向同性性质.从图 3中可以观测到快P波和S波，没有观测到慢P波，而在图 5中，除快P波和S波外，还可观测到明显的慢P波，而且慢P波表现出很强的频散，尤其在能量较强的流相分量上.这表明，在理想相界情况下，在Biot流动作用下，慢P波呈传播模式，而在Biot流动和喷射流动耦合作用下，慢P波呈扩散模式，在波场快照中无法观测到，即在两种流体流动机制作用下，慢P波表现出更强的衰减性.同时，通过对比图 3和5，可以发现，快P波在基于BISQ模型的双相各向同性介质中传播时，其传播速度明显小于在基于Biot模型的双相各向同性介质中传播时的速度，而S波的波速在基于两种模型的双相介质中传播时则无明显差异，这是由于当喷射流动处于非松弛状态时，孔隙介质的胀缩性将变差，但剪切性没有明显变化.同时这也说明喷射流动对压缩波的影响较大，而对剪切波的影响较小，这与Dvorkin和Nur(1993)的理论研究结果是一致的.
 | 图 2 η=0时BISQ模型的固相(a)和流相(b)的垂直分量(c)和水平分量(d)的25 Hz单频波场Fig. 2 The solid phase of the BISQ model(a) and fluid(b)of the vertical component(c) and the horizontal component(d)25 Hz monochromatic wave field when η=0 |
 | 图 3 η=0时BISQ模型的固相(a)和流相(b)的垂直分量(c)和水平分量(d)300 ms时的波场传播快照Fig. 3 The solid state BISQ model(a) and fluid(b)of the vertical component(c) and the horizontal component(d)300 ms wave field propagating snapshots when η=0 |
 | 图 4 η=0时Biot模型的固相(a)和流相(b)的垂直分量(c)和水平分量(d)的25 Hz单频波场Fig. 4 The solid state Biot model(a) and flow(b)the vertical component(c) and the horizontal component(d)25 Hz monochromatic wave field when η=0 |
 | 图 5 η=0时Biot模型的固相(a)和流相(b)的垂直分量(c)和水平分量(d)300 ms时的波场传播快照Fig. 5 The solid state Biot model(a) and fluid(b)of the vertical component(c) and the horizontal component(d)300 ms wave field propagating snapshots when η=0 |
下面我们考察在粘滞相界情况下，基于BISQ模型和基于Biot模型的双相各向同性介质中地震波的传播特征.这里取粘滞系数η的值为0.001 Pa·s.图 6-9给出了粘滞相界情况下，基于BISQ模型和基于Biot模型的双相各向同性介质地震波场数值模拟结果.图 6显示的是基于BISQ模型的双相各向同性介质中固相和流相垂直位移分量和水平位移分量25 Hz的单频波场，图 7给出了相应的300 ms时时间域波场传播快照；而图 8显示的是基于Biot模型的双相各向同性介质中固相和流相垂直位移分量和水平位移分量25 Hz的单频波场，图 9是相应的300 ms时时间域波场传播快照.同时，为了更直观地对比，在粘滞相界情况下，基于两种模型的双相各向同性介质中地震波传播的差异性，取(350 m，350 m)处质点的x分量振动记录进行分析.图 10给出了基于两种模型的双相各向同性介质中地震波场数值模拟得到的质点(350 m，350 m)处的x分量振动记录对比结果.
 | 图 6 η=0.001 Pa·s时BISQ模型的固相(a)和流相(b)的垂直分量(c)和水平分量(d)的25 Hz单频波场Fig. 6 The solid state BISQ model(a) and fluid(b)of the vertical component(c) and the horizontal component(d)25 Hz monochromatic wave field when η=0.001 |
 | 图 7 η=0.001 Pa·s时BISQ模型的固相(a)和流相(b)的垂直分量(c)和水平分量(d)300 ms时的波场传播快照Fig. 7 The solid state BISQ model(a) and fluid(b)of the vertical component(c) and the horizontal component(d)300 ms wave field propagating snapshots when η=0.001 |
 | 图 8 η=0.001 Pa·s时Biot模型的固相(a)和流相(b)的垂直分量(c)和水平分量(d)的25 Hz单频波场Fig. 8 The solid state Biot model(a) and fluid(b)of the vertical component(c) and the horizontal component(d)25 Hz monochromatic wave field when η=0.001 |
 | 图 9 η=0.001 Pa·s时Biot模型的固相(a)和流相(b)的垂直分量(c)和水平分量(d)300 ms时的波场传播快照Fig. 9 The solid state Boit model(a) and fluid(b)of the vertical component(c) and the horizontal component(d)300 ms wave field propagating snapshots when η=0.001 |
 | 图 10 η=0.001 Pa·s时基于BISQ模型和基于Biot模型的双相各向同性介质中质点(350 m，350 m)处水平分量记录对比图Fig. 10 Based on BISQ model and Biot model of particle two phase isotropic medium based on(350 m，350 m)contrast diagram of horizontal component recording when η=0.001 |
从图 7和9可以看出，在粘滞相界情况下，基于Biot模型和BISQ模型的双相各向同性介质中的慢P波均由于强衰减作用而无法在波场快照中观测到，波场快照中只能观测到快P波和S波，而且这两种波在基于两种模型的双相介质中传播时的速度差别与理想相界情况下相同，这些特征在图 10中更为直观，而且从图中我们可以看出在基于BISQ模型的双相各向同性介质中传播的快P波的振幅小于基于Biot模型的双相各向同性介质中的快P波，这说明Biot流动和喷射流动两种机制耦合作用造成了更大的快P波衰减.
3.2 双相各向同性介质中的地震波反射和透射特征
为了研究BISQ模型所描述的双相各向同性介质中地震波在介质分界面处的反射和透射特征，设计一个如下双层介质模型：模型网格数为200×200，网格大小为10 m×10 m，界面位于1000 m处，震源点位于(1000 m，900 m).上、下层均为双相各向同性介质，上层介质固体骨架的弹性参数为λ=5.18 GPa，μ=6.12 GPa，下层介质固体骨架弹性参数为λ=22.72 GPa，μ=16.12 GPa.两层介质中流体粘滞系数η=0.震源类型以及其他介质参数与以上均匀双相各向同性介质波场数值模拟时的情况相同.
图 11、图 12和图 13分别显示的是双层介质模型的25 Hz单频波场、300 ms时的波场快照以及地震记录.从图 11中可以看到入射波，透射波和反射波，入射波和反射波的相互干涉使得反射界面附近波场产生扰动；通过反傅氏变换将频率域波场转换到时间域，得到波场快照和合成地震记录.在图 12和13中，可以观测到快P波、反射快P波、S波、反射S波、快P波与S波之间的反射转换波、透射快P波、透射S波、快P波与S波之间的透射转换波.可以看出，由于Biot流动和喷射流动两种机制的耦合作用造成慢P波的强衰减，在波场快照和合成地震记录中无法观测到直达慢P波，以及与慢P波相关的反射、透射、转换波类，波场特征类似于单相介质的情况.
 | 图 11 双层双相各向同性介质模型固相垂直分量(a)和水平分量(b)的25 Hz单频波场Fig. 11 Double two phase isotropic medium model of solid vertical component(a) and the horizontal component(b)25 Hz monochromatic wave field |
 | 图 12 双层双相各向同性介质模型固相垂直分量(a)和水平分量(b)300 ms时的波场快照Fig. 12 Double two phase isotropic medium model of solid vertical component(a) and the horizontal component(b)wave field snapshot when 300 msspan> |
 | 图 13 双层双相各向同性介质模型固相垂直分量(a)和水平分量(b)的地震记录 1.直达快P波；2.直达S波；3.反射快P波；4.快P波产生的反射转换S波或S波产生的反射转换快P波；5.反射S波. Fig. 13 Double two phase isotropic medium model of solid vertical component(a) and the horizontal component(b)of the seismic record |
基于BISQ模型双相介质的波动方程综合考虑了两种双相介质中流体流动的形式：全局Biot流动和局部喷射流动.由于该方程是一个复变系数偏微分方程，本文首次采用了频率-空间域有限差分方法(25点)首先模拟研究了基于BISQ模型均匀双相各向同性介质中波的传播特性，并与基于Biot模型相同介质中的波做了对比，然后模拟研究了双层双相各向同性介质中波的传播规律.通过对模拟结果进行分析，得出以下结论：
(1)理想相界情况下，基于Biot模型的慢P波呈传播模式，可以在波场快照中观测到，而基于BISQ模型的慢P波呈扩散模式，在波场快照中无法观测到，即在两种流体流动机制作用下，慢P波表现出更强的衰减性.
(2)理想相界情况下，基于BISQ模型双相介质中传播的快P波的速度和振幅都小于基于Biot模型双相介质中的快P波，而S波的波速和振幅在基于两种模型的双相介质中传播时则无明显差异.这是由于当喷射流动处于非松弛状态时，孔隙介质的胀缩性将变差，但剪切性没有明显变化，同时也说明喷射流动对P波的影响较大，而对S波的影响较小.
(3)粘滞相界情况下，基于BISQ模型和Biot模型的双相各向同性介质中的慢P波均由于强衰减作用而无法在波场快照中观测到，波场快照中只能观测到快P波和S波，而且这两种波在基于两种模型的双相介质中传播时的速度和振幅差别与理想相界情况下相同.
(4)BISQ模型双层双相各向同性介质的波，由于Biot流动和喷射流动两种机制的耦合作用造成慢P波的强衰减，在波场快照和合成地震记录中无法观测到直达慢P波，以及与慢P波相关的反射、透射、转换波类，波场特征类似于单相介质的情况.
本文建立了基于BISQ机制的双相各向同性介质的频率-空间域波动方程，波场数值模拟结果证实了其合理性和有效性，可以为孔隙介质中弹性波频率域全波形反演问题的研究奠定基础.
附录A基于BISQ模型双相各向同性介质的频率-空间域波动方程，以(3a)式为例，其离散化方程为
| [1] | Berenger J P. 1994. A perfectly matched layer for the absorption of electromagnetic waves[J]. Journal of Computational Physics, 114(2): 185-200. |
| [2] | Biot M A. 1956. Theory of propagation of elastic waves in a fluid-saturated porous solid. PartⅠ: low-frequency range[J]. J. Acoust. Soc. Am., 28(2): 168-178. |
| [3] | Biot M A. 1962. Mechanics of deformation and acoustic propagation in porous media[J]. J. Appl. Phys., 33(4): 1482-1498. |
| [4] | Cheng Y F, Yang D H, Yang H Z. 2002. Biot/squirt model in viscoelastic porous media[J]. Chinese Physics Letters, 19(3): 445-448. |
| [5] | Dvorkin J, Nur A. 1993. Dynamic poroelasticity: A unified model with the squirt and the Biotmechanisms[J]. Geophysics, 58(4): 524-533. |
| [6] | Min D J, Shin C, Kwon B D, et al. 2000. Improved frequency-domain elastic wave modelingusingweighted-averaging difference operators[J]. Geophysics, 65(3): 884-895. |
| [7] | Nie J X, Yang D H, Yang H Z. 2004. Wave dispersion and attenuation in partially saturated sandstones[J]. Chinese Physics Letters, 21(3): 572-575. |
| [8] | Nie J X, Yang D H, Yang H Z. 2008. A generalized viscoelastic Biot/squirt model for clay-bearing sandstones in a wide range of permeabilities[J]. Applied Geophysics, 5(4): 249-260. |
| [9] | Parra J O. 1997. The transversely isotropic poroelastic wave equation including the Biot and the squirt mechanisms: Theory and application[J]. Geophysics, 62(1): 309-318. |
| [10] | Wang Z J, He Q D, Wang D L. 2008. The numerical simulation for a 3D two-phase anisotropic medium based on BISQ model[J]. Applied Geophysics, 5(1): 24-43. |
| [11] | Yang D H, Zhang Z J. 2000. Effects of the Biot and the squirt-flow coupling interaction on anisotropic elastic waves[J]. Chinese Science Bulletin, 45(23): 2130-2138. |
| [12] | Yang D H, Zhang Z J. 2002. Poroelastic wave equation including the Biot/squirt mechanism and the solid/fluid coupling anisotropy[J]. Wave Motion, 35(3): 223-245. |
| [13] | 高凤霞, 刘财, 冯晅,等. 2013. 几种优化方法在频率域全波形反演中的应用效果及对比分析研究[J]. 地球物理学进展, 28(4): 2060-2068, doi: 10.6038/pg20130450. |
| [14] | 郭智奇, 刘财, 杨宝俊,等. 2007. 粘弹各向异性介质中地震波场模拟与特征[J]. 地球物理学进展, 22(3): 804-810, doi: 10.3969/j.issn.1004-2903.2007.03.022. |
| [15] | 李红星, 刘财, 陶春辉. 2007. 基于横向各向同性BISQ模型的弹性波高阶交错网格有限差分数值模拟[J]. 石油地球物理勘探, 42(6): 686-693. |
| [16] | 李红星, 陶春辉, 周建平,等. 2009. 双相各向异性介质中波场速度与衰减特征分析[J]. 石油地球物理勘探, 44(4): 457-465. |
| [17] | 刘财, 李红星, 陶春辉. 2007. 基于BISQ模型的各向同性孔隙介质弹性波三维交错网格高阶有限差分数值模拟[J]. 世界地质, 26(4): 501-508. |
| [18] | 刘财, 张智, 邵志刚,等. 2005. 线性粘弹体中地震波场伪谱法模拟技术[J]. 地球物理学进展, 20(3): 640-644. |
| [19] | 刘璐, 刘洪, 刘红伟. 2013. 优化15点频率-空间域有限差分正演模拟[J]. 地球物理学报, 56(2): 644-652, doi: 10.6038/cjg20130228. |
| [20] | 聂建新, 杨顶辉, 杨慧珠. 2004. 基于非饱和多孔隙介质BISQ模型的储层参数反演[J]. 地球物理学报, 47(6): 1101-1105. |
| [21] | 聂建新, 杨顶辉, 巴晶. 2010. 含泥质低孔渗各向异性黏弹性介质中的波频散和衰减研究[J]. 地球物理学报, 53(2): 385-392, doi: 10.3969/j.issn.0001-5733.2010.02.016. |
| [22] | 皮红梅, 蒋先艺, 刘财,等. 2009.波动方程数值模拟的三种方法及对比[J]. 地球物理学进展, 24(2): 391-397, doi: 10.3969/j.issn.1004-2903.2009.02.004. |
| [23] | 吴国忱, 罗彩明, 梁楷. 2007. TTI介质弹性波频率-空间域有限差分数值模拟[J]. 吉林大学学报(地球科学版), 37(5): 1023-1033. |
| [24] | 杨宽德, 杨顶辉, 王书强. 2002a. 基于Biot-Squirt方程的波场模拟[J]. 地球物理学报, 45(6): 853-861. |
| [25] | 杨宽德, 杨顶辉, 王书强. 2002b. 基于横向各向同性BISQ方程的弹性波传播数值模拟[J]. 地震学报, 24(6): 599-606. |
| [26] | 杨宽德, 宋国杰, 李静爽. 2011. Biot流动和喷射流动耦合作用下波传播的FCT紧致差分模拟[J]. 地球物理学报, 54(5): 1348-1357, doi: 10.3969/j.issn.0001-5733.2011.05.024. |
| [27] | 张显文, 王德利, 王者江,等. 2010. 基于BISQ机制三维双相正交裂隙各向异性介质衰减及频散方位特性研究[J]. 地球物理学报, 53(10): 2452-2459, doi: 10.3969/j.issn.0001-5733.2010.10.019. |
| [28] | 张智, 刘财, 邵志刚,等. 2005. 伪谱法在常Q粘弹介质地震波场模拟中的应用效果[J]. 地球物理学进展, 20(4): 945-949, doi: 10.3969/j.issn.1004-2903.2005.04.010. |