地球物理学报  2018, Vol. 61 Issue (3): 926-937   PDF    
求解双相和黏弹性介质波传播方程的间断有限元方法及其波场模拟
张金波1, 杨顶辉1 , 贺茜君2, 马啸3     
1. 清华大学数学科学系, 北京 100084;
2. 海南大学信息学院数学系, 海口 570228;
3. 西北工业大学应用数学系, 西安 710072
摘要:间断有限元(Discontinuous Galerkin:DG)方法具有低数值频散、网格剖分灵活、能模拟地震波在复杂介质中传播等优点.因此,本文将一种新的DG方法推广到双相和黏弹性等复杂介质的地震波场模拟,发展了求解Biot弹性波方程和D'Alembert介质波动方程的DG方法.首先通过引入辅助变量将Biot双相介质弹性波方程和D'Alembert介质波动方程转化为关于时间-空间的一阶偏微分方程组,然后对该方程组进行DG空间离散,得到半离散化的常微分方程组.最后,对此常微分方程组,应用加权的Runge-Kutta格式进行时间推进计算.数值结果表明,DG方法可以有效地求解Biot双相介质弹性波方程和D'Alembert介质波动方程,并能很好地压制因离散求解波动方程而产生的数值频散,获得清晰的各种地震波震相.
关键词: 数值模拟      间断有限元      双相介质      黏弹性介质     
Discontinuous Galerkin method for solving wave equations in two-phase and viscoelastic media
ZHANG JinBo1, YANG DingHui1, HE XiJun2, MA Xiao3     
1. Department of Mathematical Sciences, Tsinghua University, Beijing 100084, China;
2. Department of Mathematics, Colledge of Information Science and Technology, Hainan University, Haikou, 570228, China;
3. Department of Applied Mathematics, Northwestern Polytechnical University, Xi'an 710072, China
Abstract: The Discontinuous Galerkin (DG) method has great advantages in suppressing numerical dispersion and dealing with complex structures. Therefore, in this paper, we apply a new DG method to numerical simulations in two-phase and viscoelastic media, and suggest a DG method to solve both Biot elastic wave equations and the D'Alembert wave equations. For this, we first transform the Biot equations and the D'Alembert wave equations into a system of first-order equations with respect to time-space by introducing auxiliary variables. Then we transform the first-order equations into a semi-discrete ordinary differential equation (ODE) system using the DG method. Finally, we use a weighted Runge-Kutta method to solve the ODE system. The numerical results show that the DG method works very well for solving the Biot elastic wave equations and D'Alembert wave equations, and can effectively suppress the numerical dispersion and provide accurate information on the wave-field.
Key words: Numerical simulation    Discontinuous Galerkin method    Two-phase medium    Viscoelastic medium    
0 引言

在实际问题中,地下介质通常是由固相、液相和(或)气相组成的双相或多相系统.在这种双相或多相复杂系统中,不同地下介质的岩石孔隙度、孔隙流体和固体骨架的弹性模量、流体密度和骨架密度等差别较大.在此情况下,基于单相纯固体或纯流体的弹性波理论已不能准确描述实际情况.为了适应复杂介质情况下的勘探要求,以期更准确地描述实际地层结构和性质,地球物理学家建立了双相介质中的弹性波理论.这方面最早的理论是Gassmann模型(Gassmann, 1951),最重要的是Biot(1956, 1962)建立了著名的双相介质弹性波传播模型,这一模型奠定了多孔隙双相介质中弹性波传播的理论基础.Plona等(1980)在实验中观察到了Biot理论预测到的慢P波的存在,从而证实了Biot理论的正确性.基于这一理论,许多研究者深入细致地研究了油气储层介质中的波传播规律(Dvorkin et al., 1980;Parra, 1997;等等).

另一方面,黏弹性介质是介于固体完全弹性和流体黏性介质之间的一种复杂介质(牛滨华和孙春岩,2007),黏弹性介质的非弹性性质(介质会发生永久性形变),比弹性介质更接近实际情况.典型的黏弹性介质有D′Alembert介质、Boltzmann介质、Kelvin介质等.其中的D′Alembert粘弹性介质模型具有相对简单,但又能更准确刻画地震传播规律的波动方程,它是直接在弹性波方程基础上加入了黏性项而获得的一类黏弹性介质模型(牛滨华和孙春岩, 2004),在描述地震波传播方面十分重要.

总之,双相介质和黏弹性介质中波传播问题受到国内外很多学者的关注,在求解这两类波动方程及其波场模拟方面开展了大量的研究工作,获得了各种不同的数值求解方法.其中常用的数值求解方法有有限差分法(Kelly et al., 1976; 董良国等, 2004)、有限元法(杨顶辉, 2002; 李伟华等, 2010)和伪谱法(马永旺等, 2001),等等.针对双相介质的特殊性,de la Puente等(2008)等对双相介质中波传播的数值模拟做了很多研究;在国内,杨顶辉(2002)基于固-流相对位移给出了求解双相各向异性介质中弹性波方程的有限元方法,杨德宽等(2002)给出了求解双相各向同性介质中弹性波方程的有限差分法,Wang等(2008)利用交错网格法对三维双相介质中波传播进行了数值模拟,Yang等(2007)利用近似解析离散化方法及其改进版模拟了双相介质中的波传播.针对粘弹性介质,Carcione等(Carcione et al., 1988Carcione,1993)利用伪谱法对黏弹性介质中的波动方程进行了数值求解;杨慧珠和杜启振(2003)利用有限元法对黏弹性介质的波传播进行了数值模拟;在这些数值方法中,有限差分法因编程简单.储存量小、计算速度快、稳定性好(王美霞, 2012)等优点而受到普遍欢迎,但是粗网格条件下很容易产生严重的数值频散(Yang et al., 2003, 2006; Zheng et al., 2006);有限元法易处理复杂边界、网格剖分灵活,但计算效率低,不利于大规模问题求解(佘德平, 2004);伪谱法易于提高精度,但不能灵活适应复杂地质模型(牟永光和裴正林, 2005),并行效率低.

间断有限元方法最早是由Reed和Hill(1973)在求解中子输运方程时提出,至今已发展了多种多样的间断有限元方法,其中Cockburn和Shu(1989, 2001)提出的求解双曲守恒方程的龙格库塔间断有限元(Runge-Kutta Discontinuous Galerkin,简称RKDG)方法应用最广泛.RKDG是一种基于通量形式的间断有限元方法,它从局部单元入手,把有限元方法中需要求解的大规模矩阵转化成了局部单元上的小规模矩阵.相比其它几类方法,RKDG方法在处理复杂边界问题、起伏地表问题、模拟间断信号,以及在计算中更容易提高精度等方面具有很大的优势,可以更有效地模拟复杂介质中的地震波场.贺茜君(2015)分析了RKDG方法的数值频散和稳定性,并对单相介质中的声波和弹性波传播进行了数值模拟.结果表明,RKDG方法具有低数值频散等优良性质.在此基础上,He等发展出了一种加权的Runge-Kutta DG(WRKDG)方法(He et al., 2015),它采用隐式对角Runge-Kutta方法作为时间离散,并通过级数截断方法将其显式化,以达到提高库朗数进而加快计算速度的目的;Yang等将ONADM(optimal nearly analytic method,Yang et al., 2006)和WRKDG (He et al., 2015)两种方法相结合,发展了一种新的杂交算法(Yang et al., 2016),进一步提高了计算效率.在国际上,Käser和Dumbser课题组及其成员发展了一种基于通量函数的任意阶导数的时间步进DG方法(arbitrary high-order derivatives time stepping discontinuous Galerkin method: ADER-DG), 并将其应用到了求解弹性波方程,以及孔隙介质等模型中(Käser et al., 2007; De La Puente et al., 2008).Dupuy等人也利用DG方法对频率域中多孔隙复杂介质中的波传播过程进行了模拟(Dupuy et al., 2011).本文将He等提出的加权Runge-Kutta DG(WRKDG)方法引入到数值求解双相和黏弹性介质波方程中,发展了求解双相Biot方程和粘弹性D′Alembert介质波动方程的加权Runge-Kutta DG方法.

1 求解双相介质弹性波方程的DG方法

在对地下地质结构进行实际勘探时,弹性介质理论是一种相对简化的理论.而实际的含油气地下储层大多都是双相(固相和流相)或是多相(固相、流相和气相)介质.通常情况下含流体双相介质地震波传播理论中包含了更多的储层物性参数,可以更好地适应实际的地质结构.因此通过数值模拟研究地震波在多孔隙介质中的传播规律、认识储层参数的地震响应特征具有重要的实际意义.

Biot于1956年建立的弹性波在饱和多孔隙岩石中的传播理论具有如下形式:

(1)

x-z二维情况下的弹性波传播方程为

(2)

其中u1, u3, U1, U3分别为固体xz方向、流体xz方向的位移;ρ11, ρ12, ρ22表示固体、流体和固-流耦合的等效密度.λμ为固体骨架等效弹性模量,QR为固-流耦合等效弹性模量和流体等效弹性模量;b表示耗散系数,可以通过介质渗透率k,粘滞系数η,孔隙度φ计算:

(3)

对于方程(2),我们令U=(u1, u3, U1, U3),引入变量P=(p1, p3, P1, P3),Q=(q1, q3, Q1, Q3),且满足

(4)

加入震源项后,则方程(2)可以转化为

(5)

其中 =(əx, əz)是散度算子.令W =(u1, U1, u3, U3, p1, P1, p3, P3, q1, Q1, q3, Q3)T,并定义F(W)为

(6)

再令

这里I4表示4阶单位矩阵,表示震源项.则方程(5)可以改写为如下形式:

(7)

亦即

(8)

其中,H-1表示对矩阵H求逆运算.

当线性算子B=0时,表示在无耗散时的理想状态弹性波方程.下面构造求解方程(8)的DG方法.

假设计算区域为Ω,将其划分为一些子区域{Ωi},并且子区域{Ωi}之间不重叠.本文只考虑二维情况下将Ω划分为矩形网格的情况.取试验函数空间为

(9)

这里v表示在Ωi上不超过k次,在Ωi之外为0的多项式.因此,与传统有限元不同的是,这里的试验函数可以在{Ωi}之间存在间断,不需要在整个Ω区域上连续.

对于方程(8),引入试验函数v之后在Ωi上进行积分得

(10)

利用Green公式,可将方程(10)变为如下形式:

(11)

其中右端两项分别与震源和耗散有关,左端项中n表示积分边界əΩi的外法向量.方程(11)中的WəΩi上可能存在间断,因此在əΩi上的被积函数项H-1F(Wn需要用数值通量(Wint, Wext, n)来替代,其中WintWəΩi内部的值,WextWəΩi外部的值.这样方程(11)就转化为

(12)

根据He等(2015)的分析,结合上述基函数的定义可以看出,DG方法的基函数构造不需要定义在节点上,而且其选取对于数值模拟结果影响较小.但是从方程(12)可以看到,在计算方程左端项时,取合适的基函数可以在计算左端项的过程中避免计算一个线性代数方程,因而较常见的是取一组正交基底作为DG方法的基函数,这样可以有效提高计算效率.

我们现在把方程(12)中的WΩi中按照基函数展开:

(13)

其中N为基函数个数.

按照方程(13)中W的基函数展开形式,将其代入方程(12)可得

(14)

至此,方程(2)转化成为了一个半离散化的ODEs.

下面我们对方程(14)左端的每一项做一些说明.对于其中的第一项,利用基函数wl的正交性,可以转化为

(15)

而方程(14)中的第二项可写为

(16)

这里A1, A2分别为

对于方程(14)中的第三项:

(17)

这里j表示与i相邻的矩形单元,Lj分别为单元i与单元j的边界,Aint, Aext分别为

(18)

其中n1, n2分别为法向量n的分量, (Wint, Wext, n)为

C的取值为maxmin(a, b)≤s≤max(a, b)|H-1F′(sn|.

至此,利用数值积分公式对方程(14)中的每一项进行数值计算,完成了DG空间离散格式构造.令

则方程(14)可以转化为如下方程:

(19)

其中L是与方程(2)中的FHB相关的线性算子,F是和震源相关的项.利用加权的Runge-Kutta格式对方程(19)进行时间推进计算,其格式为

(20)

其中

以迭代计算得到的K1(n)K2(n)加权逼近K(n),以T(n)作为初始值迭代计算K(n)

其中以迭代计算得到的K 1(n)K2(n)逼近K(n).从而利用该格式对方程(19)进行时间推进计算.

2 黏弹性介质波动方程的DG方法

黏弹性介质是介于黏性和弹性之间的一种介质,其中最常见的是D′Alembert介质.D′Alembert介质中的波动方程是直接在弹性介质波动方程中加入了耗散项,以此来描述地震波在D′Alembert介质中的衰减规律.二维均匀各向同性D′Alembert介质中波动方程为

(21)

这里u1, u3分别表示xz方向的位移,λμ是Lamé常数,ρ表示介质密度,η表示耗散系数.

2.1 D′Alembert纵波方程的DG方法

D′Alembert介质的矢量波场可以分解为纵波和横波波场,其中x-z二维介质中的纵波方程为

(22)

为了给出求解方程(22)的DG方法,首先将其转化为关于时间-空间一阶偏导数的方程组.于是,对u分别引入变量p, q,且满足

(23)

,则方程(22)可以转化为

(24)

W=(u, p, q)TF(W)和g(W)分别为

则方程(24)可以表示为

(25)

至此,我们将D′Alembert介质波方程转化成了关于时间和空间的一阶偏导数方程组.现在把上式中的W利用式子(13)中的基函数展开可以得到

(26)

至此即完成D′Alembert介质波方程的DG空间离散,利用加权的Runge-Kutta进行时间推进计算即可求解.

2.2 P-SV波方程的DG方法

在2.1小节中已经给出了D′Alembert介质中纵波方程的关于时间-空间一阶偏导的形式,本小节将给出黏弹性介质中P-SV波方程的关于时间-空间一阶偏导的形式.其思路与2.1小节思路一致.

对于方程(21),为了构造其DG方法,首先将其转化为关于时间-空间一阶偏导的形式.对U= (u1, u3)T分别引入变量P=(p1, p3)TQ=(q1, q3)T其满足:

(27)

,则其可以转化为

(28)

W=(u1, u3, p1, p3, q1, q3)TF(W)和g(W)分别为

则方程(28)可以同样表示为

(29)

对于方程(29),可以将其转为方程(26)的形式,并同样通过加权的Runge-Kutta格式(20)进行时间推进计算.

3 数值模拟 3.1 Biot双相介质弹性波方程

在第一个数值试验中,我们选取一个单层介质模型,计算区域大小为0≤x≤10.0 km, 0≤z≤10.0 km.在区域设置(200×200)矩形网格,选取震源为纵波震源,设置震源位置为选取区域的中心位置(5.0 km, 5.0 km),选取震源函数为

(30)

其中震源的频率为f0=30.0 Hz.我们利用前面给出的WRKDG方法对其进行求解,方法中的参数r=/5,η=0.5.选取密度为ρs=2.32 g·cm-3, ρf=1.08 g·cm-3, ρa=0.083 g·cm-3,孔隙度为φ=0.2,其它参数为λ=9.6 GPa, μ=15.0 GPa, Q=0.953 GPa, R=0.331 GPa.

图 1图 2给出了时间步长为Δt=0.0014 s,利用WRKDG方法模拟得到的T=1.0 s时刻在理想状态(b=0)条件下固体和和流体分别在xz方向的位移分量u1, u3, U1, U3的波场快照图.从图中可以看出,波场快照中的快P波、S波和慢P波很清晰,无明显可见数值频散,表明DG方法能很好地压制数值频散,能给出清晰的波场模拟结果.

图 1 T=1.0 s时刻固体位移x方向(a)和z方向(b)的波场快照 Fig. 1 Snapshots for solid displacement on the (a) x-direction and (b) z-direction at T=1.0 s
图 2 T=1.0 s时刻流体位移x方向(a)和z方向(b)的波场快照 Fig. 2 Snapshots for liquid displacement on the (a) x-direction and (b) z-direction at T=1.0 s

图 3图 4是取耗散系数b=4.0后计算得到的波场快照.由于在耗散状态条件下慢P波会快速衰减,故从图 34的波场快照中观察不到慢P波.

图 3 耗散条件下T=1.0 s时刻固体位移x(a)和z(b)分量的波场快照 Fig. 3 Snapshots for solid displacement on the (a) x-direction and (b) z-direction at T=1.0 s with dissipation
图 4 耗散条件下T=1.0 s时刻流体位移x(a)和z(b)分量的波场快照 Fig. 4 Snapshots for liquid displacement on the (a) x-direction and (b) z-direction at T=1.0 s with dissipation
3.2 D′Alembert介质纵波波场模拟

在这个数值试验中,我们选取均匀介质模型,利用WRKDG方法对D′Alembert介质中的纵波传播进行模拟.实验中选取计算区域大小为0≤x≤10.0 km, 0≤z≤10.0 km,200×200的矩形网格,密度参数为ρ=2.1g·cm-3.震源位于计算区域中心,震源频率为f0=20 Hz,震源函数为

(31)

图 5给出了耗散参数r=2.0和8.0时纵波在1.4 s时刻的波场快照,波场非常清晰.图 6给出了不同耗散参数,即不同r值情况下在接收器(7.2 km, 7.8 km)处的波形记录图.从图 6可以看出,在耗散参数r增大时,波形振幅会减小,从而验证了D′Alembert介质纵波在传播过程中的衰减特性.

图 5 D′Alembert介质中纵波在T=1.4 s时刻的波场快照((a) r=2;(b) r=8) Fig. 5 Snapshots of P-wave in the D′Alembert medium for (a) r=2 and (b) r=8
图 6 不同耗散参数下的波形记录图 Fig. 6 Waveforms for different dissipation parameters
3.3 Marmousi模型

在这个试验中,我们选取Marmousi速度模型,利用本文的WRKDG方法模拟非均匀D′Alembert介质中的纵波波场.Marmousi速度模型有很强的非均匀性,其速度变化范围为1500~5500 m·s-1.选取计算区域大小为0≤x≤9216 m,0≤z≤2928 m;耗散参数为r=2.0.选取384×122的矩形网格,震源位于(4609 m, 12 m),震源频率f0=25 Hz,震源函数为方程(31)的形式.

图 7是Marmousi速度模型示意图,图 8图 9分别是纵波在0.6 s和1.2 s时刻的波场快照.图 89表明,在非均匀介质情况下,DG方法也能给出无可见数值频散的波场结果.这说明WRKDG方法可以有效地模拟非均匀复杂D′Alembert介质中的纵波传播.

图 7 Marmousi速度模型 Fig. 7 Marmousi model
图 8 对于Marmousi模型,D′Alembert介质中纵波在0.6 s时刻的波场快照 Fig. 8 Snapshots of P-wave in the D′Alembert medium at T=0.6 s for the Marnousi model
图 9 对于Marmousi模型,D′Alembert介质中纵波在1.2 s时刻的波场快照 Fig. 9 Snapshots of P-wave in the D′Alembert medium at T=1.2 s for the Marnousi model
3.4 P-SV波波场模拟

在这个例子中,选取计算区域大小为0≤x≤ 10.0 km,0≤z≤10.0 km;模型参数为ρ=2.1 g·cm-3, η=6.3, λ=9.6 GPa,μ=15.0 GPa.选取网格点数为200×200的矩形网格,震源位于计算区域的中心,震源频率f0=20 Hz,震源函数为方程(30)的形式,在xz方向加相同力源.

图 10是取r=4.0,P-SV波在0.6 s时刻的波场快照.图 10的左图和右图分别给出了位移u1, u3的非常清晰的波场快照.这表明WRKDG方法可以同时有效地模拟P波和SV波波场.图 11给出了不同耗散参数情况下于接收器(5.4 km, 5.8 km)处的波形记录.从图 11可以看出,随耗散参数增大波振幅减小,这说明了黏弹性介质中P波和SV波在传播过程中的衰减特性.

图 10 P和SV波在0.6 s时刻的波场快照((a) u1;(b) u3) Fig. 10 Snapshots of P-wave and SV-wave at T=0.6 s: (a) u1 and (b) u3
图 11 P和SV波在不同耗散参数情况下的波形记录((a) u1;(b) u3) Fig. 11 Waveforms of P-wave and SV-wave at T=1.2 s for different dissipation parameters: (a) u1 and (b)u3
4 结论

本文发展了数值求解Biot双相介质弹性波方程和D′Alembert介质波动方程的DG方法.首先,针对Biot双相介质弹性波方程和D′Alembert介质波动方程,通过对原方程引入变量并且进行波动方程变换,将包含空间和时间二阶偏导数的原波动方程转化成关于时间和空间一阶偏导的波动方程组.然后选取合适的数值通量、基函数、网格单元和数值求积公式,应用DG方法对Biot方程和D′Alembert介质波动方程所对应的一阶方程组进行空间离散,将其转化成半离散化的常微分方程组,并应用加权的Runge-Kutta格式进行时间推进计算.最后针对不同双相介质和D′Alembert介质的地质模型,利用DG方法进行了地震波场模拟.数值结果表明,DG方法能有效地求解Biot双相介质弹性波方程和D′Alembert介质波动方程,不仅能很好地压制数值频散,获得包含不同波相(快、慢P波和S波)的清晰的地震波场,而且能充分体现耗散介质中波的衰减特征.

需要指出的是,尽管DG方法具有压制数值频散、可灵活处理间断面和起伏地表等诸多优势,但DG方法所需的计算量和存储量相对于其它数值方法还是比较大,因此,在未来的研究和实际应用中,应重点考虑如何提高DG方法的计算效率.例如,对于具有复杂地质结构的双相和黏弹性介质波传播模拟,可考虑发展多尺度、非一致网格的数值模拟方法,或考虑将DG方法与其它高效数值方法相结合,发展新型的混合方法等,这些都是我们今后可以继续深入研究的方向.

参考文献
Biot M A. 1956. Theory of propagation of elastic waves in a fluid-saturated porous solid, I Low frequency range. The Journal of the Acoustical Society of America, 28(2): 168-178. DOI:10.1121/1.1908239
Biot M A. 1962. Mechanics of deformation and acoustic propagation in porous media. Journal of Applied Physics, 33(4): 1482-1498. DOI:10.1063/1.1728759
Carcione J M, Kosloff D, Kosloff R. 1988. Viscoacoustic wave propagation simulation in the earth. Geophysics, 53(6): 769-777. DOI:10.1190/1.1442512
Carcione J M. 1993. Seismic modeling in viscoelastic media. Geophysics, 58(1): 110-120. DOI:10.1190/1.1443340
Cockburn B, Shu C W. 1989. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws Ⅱ:General framework. Mathematics of Computation, 52(186): 411-435.
Cockburn B, Shu C W. 2001. Runge-kutta discontinuous galerkin methods for convection-dominated problems. Journal of Scientific Computing, 16(3): 173-261. DOI:10.1023/A:1012873910884
de la Puente J, Dumbser M, Käser M, et al. 2008. Discontinuous galerkin methods for wave propagation in poroelastic media. Geophysics, 73(5): T77-T97. DOI:10.1190/1.2965027
Dong G L, Li P M. 2004. Dispersive problem in seismic wave propagation numerical modeling. Natural Gas Industry, 24(6): 53-56.
Dvorkin J, Nur A. 1993. Dynamic poroelasticity:A unified model with the squirt and the Biot mechanisms. Geophysics, 58(4): 524-533. DOI:10.1190/1.1443435
Dupuy B, De Barros L, Garambois S, et al. 2011. Wave propagation in heterogeneous porous media formulated in the frequency-space domain using a discontinuous Galerkin method. Geophysics, 76(4): N13. DOI:10.1190/1.3581361
Gassmann F. 1951. Elastic waves through a packing of spheres. Geophysics, 16(4): 673-685. DOI:10.1190/1.1437718
He X J, Yang D H, Wu H. 2015. A weighted Runge-Kutta discontinuous Galerkin method for wavefield modeling. Geophysical Journal International, 200(3): 1389-1410. DOI:10.1093/gji/ggu487
He X J. 2015. Discontinuous Galerkin method for solving wave equations[Ph. D. thesie] (in Chinese). Beijing: Tsinghua University. https://en.wikipedia.org/wiki/DG_method
Käser M, Dumbser M, De La Puente J, et al. 2007. An arbitrary high-order Discontinuous Galerkin method for elastic waves on unstructured meshes-Ⅲ. Viscoelastic attenuation. Geophysical Journal International, 168(1): 224-242. DOI:10.1111/gji.2007.168.issue-1
Kelly K R, Ward R W, Treitel S, et al. 1976. Synthetic seismograms; a finite-difference approach. Geophysics, 41: 2-27. DOI:10.1190/1.1440605
Li W H, Liu Q H, Zhao C G. 2010. Three-dimensional viscous-spring boundaries in domain and dynamic analysis using explicit finite element method of saturated porous medium. Chinese journal of Geophysics, 53(10): 2460-2469. DOI:10.3969/j.issn.0001-5733.2010.10.020
Ma Y W, Zhan Z B, Gu H M. 2001. Numerical modeling for seismic wave propagation via pseudo spectral method. Geophysical Prospecting for Petroleum, 40(2): 42-48.
Mou Y G, Pei Z L. 2005. Seismic Numerical Simulation in Three-dimensional Complex Media. Beijing: Petroleum Industry Press.
Niu B H, Sun C Y. 2004. Developing theory of propagation of seismic waves-medium model and propagation of seismic waves. Progress in Geophysics, 19(2): 255-263.
Niu B H, Sun C Y. 2007. Half-Space Homogeneous Isotropic. Beijing: Geological Publishing House.
Parra J O. 1997. The transversely isotropic poroelastic wave equation including the Biot and the squirt mechanisms:theory and application. Geophysics, 62(1): 309-318. DOI:10.1190/1.1444132
Plona T J. 1980. Observation of a second bulk compressional wave in a porous medium at ultrasonic frequencies. Applied Physics Letters, 36(4): 259-261. DOI:10.1063/1.91445
Reed W H, Hill T R. 1973. Triangular mesh methods for the Neutron transport equation. Los Alamos Report. LA-UR, 73-479.
She D P. 2004. Wavefield numerical simulation technology. Progress in Exploration Geophysics, 27(1): 16-21.
Wang M X. 2012. Symplectic stereomodelling method for wave equations in two-phase and viscoelastic media and its numerical simulations[Master's thesis] (in Chinese). Beijing: Tsinghua University.
Wang Z, He Q, Wang D. 2008. The numerical simulation for a 3D two-phase anisotropic medium based on BISQ model. Applied Geophysics, 5(1): 24-34. DOI:10.1007/s11770-008-0011-9
Yang D H. 2002. Finite element method of the elastic wave equation and wavefield simulation in two-phase anisotropic media. Chinese Journal of Geophysics, 45(4): 575-583.
Yang D H, He X J, Ma X, et al. 2016. An optimal nearly analytic discrete-weighted Runge-Kutta discontinuous Galerkin hybrid method for acoustic wavefield modeling. Geophysics, 81(5): T251-T263. DOI:10.1190/geo2015-0686.1
Yang D H, Peng J M, Lu M, et al. 2006. Optimal nearly-analytic discrete approximation to the scalar wave equation. Bulletin of the Seismological Society of America, 96(3): 1114-1130. DOI:10.1785/0120050080
Yang D H, Song G J, Chen S, et al. 2007. An improved nearly analytical discrete method:an efficient tool to simulate the seismic response of 2-D porous structures. Journal of Geophysics and Engineering, 4(1): 40-52. DOI:10.1088/1742-2132/4/1/006
Yang D H, Teng J W, Zhang Z J, et al. 2003. A nearly analytic discrete method for acoustic and elastic wave equations in anisotropic media. Bulletin of the Seismological Society of America, 93(2): 882-890. DOI:10.1785/0120020125
Yang K D, Yang D H, Wang S Q. 2002. Wave-fields simulation based on the biot-squirt equation. Chinese Journal of Geophysics, 45(6): 853-861.
Zheng H S, Zhang Z J, Liu E R. 2006. Non-linear seismic wave propagation in anisotropic media using the flux-corrected transport technique. Geophysical Journal International, 165(3): 943-956. DOI:10.1111/gji.2006.165.issue-3
董良国, 李培明. 2004. 地震波传播数值模拟中的频散问题. 天然气工业, 24(6): 53–56.
贺茜君. 2015. 求解波动方程的间断有限元方法及其波场模拟[博士论文]. 北京: 清华大学. http://cdmd.cnki.com.cn/Article/CDMD-10003-1016712358.htm
李伟华, 刘清华, 赵成刚. 2010. 饱和多孔介质三维时域黏弹性人工边界与动力反应分析的显式有限元法. 地球物理学报, 53(10): 2460–2469. DOI:10.3969/j.issn.0001-5733.2010.10.020
马永旺, 詹正彬, 顾汉明. 2001. 伪谱法地震波传播数值模拟. 石油物探, 40(2): 42–48.
牟永光, 裴正林. 2005. 三维复杂介质地震数值模拟. 北京: 石油工业出版社.
牛滨华, 孙春岩. 2004. 地震波理论研究进展-介质模型与地震波传播. 地球物理学进展, 19(2): 255–263.
牛滨华, 孙春岩. 2007. 半无限空间各向同性黏弹性介质与地震波传播. 北京: 地质出版社.
佘德平. 2004. 波场数值模拟技术. 油气藏评价与开发, 27(1): 16–21.
王美霞. 2012. 双相及黏弹性介质中波传播方程的保辛算法及其波场模拟[硕士论文]. 北京: 清华大学. http://cdmd.cnki.com.cn/Article/CDMD-10003-1013017068.htm
杨顶辉. 2002. 双相各向异性介质中弹性波方程的有限元解法及波场模拟. 地球物理学报, 45(4): 575–583.
杨德宽, 杨顶辉, 王书强. 2002. 基于Biot-Squirt方程的波场模拟. 地球物理学报, 45(6): 853–861.
杨慧珠, 杜启振. 2003. 方位各向异性黏弹性介质波场有限元模拟. 物理学报, 52(8): 2010–2014. DOI:10.7498/aps.52.2010