地球物理学报  2013, Vol. 56 Issue (3): 961-970   PDF    
引用本文
孙伟家, 符力耘, 管西竹 等. 页岩气地震勘探中页岩各向异性的地震模拟研究. 地球物理学报, 2013, 56(3): 961-970, doi: 10.6038/cjg20130324.  
Sun W J, Fu L Y, Guan X Z, et al. A study on anisotropy of shale using seismic forward modeling in shale gas exploration. Chinese J. Geophys. (in Chinese), 2013, 56(3): 961-970, doi: 10.6038/cjg20130324.  
页岩气地震勘探中页岩各向异性的地震模拟研究
孙伟家 , 符力耘 , 管西竹 , 魏伟     
中国科学院地质与地球物理研究所, 北京 100029;中国科学院地球深部研究重点实验室, 北京 100029
收稿日期 2012-05-15, 2013-02-01 收修定稿
基金项目: 国家自然科学基金项目(41104079)资助
作者简介: 孙伟家, 男, 1984年生, 2010年毕业于中国科学院地质与地球物理研究所获得博士学位, 副研究员, 主要从事地震正演模拟和偏移成像的方法研究.E-mail:swj@mail.iggcas.ac.cn
摘要: 页岩的强各向异性特征挑战地震波传播数值模拟方法的精度极限, 特别是易引起频散的高频波(>100 Hz)传播的数值模拟.鉴于目前我国页岩气地震勘探主要以常规地震声波资料为主, 本文首先介绍了一种VTI介质声波方程的任意偶数阶有限差分数值模拟方法, 并讨论其稳定性条件和吸收边界条件.任意偶数阶的差分解可有效提高计算精度, 压制数值频散噪声.针对页岩较强的各向异性特征, 本文比较了不同模型的声波方程和VTI介质声波方程计算得到的地震响应.数值结果表明, 各向异性对地震波的运动学(相位)和动力学(振幅)特性影响作用明显.因此, 在页岩气地震勘探资料处理的各个环节必须充分考虑各向异性的影响, 采取有别于常规油气勘探的处理流程和技术.
关键词: 页岩      各向异性      有限差分      声波方程     
A study on anisotropy of shale using seismic forward modeling in shale gas exploration
SUN Wei-Jia, FU Li-Yun, GUAN Xi-Zhu, WEI Wei     
Key Laboratory of the Earth's Deep Interior, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China
Abstract: The accuracy of numerical methods is challenged by strong anisotropies shown by shale, especially for high-frequency components (>100 Hz). Since the shale gas exploration still employs the acoustic seismic data, we first introduce an artificial even-order differential solution to VTI acoustic wave equation in the paper. The high-order solution can improve the computational accuracy and suppress numerical dispersion effectively. To show the effect of strong anisotropy of shale, we compare the seismic response of VTI acoustic wave equation with that of acoustic wave equation for different models. Numerical results show that anisotropy has large influences on the kinetics (phases) and the dynamics (amplitudes) of seismic waves. We conclude that the anisotropy must be considered when processing seismic data..
Key words: Shale      Anisotropy      Finite-difference      Acoustic wave equation     
1 引言

页岩气是一种重要的非常规天然气, 其勘探开发已在北美洲(尤其是美国)获得了成功.页岩广泛分布在我国的盆地, 页岩气富集地质条件优越, 页岩气资源开发具有良好的前景及潜力[1-2].我国已将页岩气列为第172种独立矿种.页岩(shale)是由粒径小于0.0039mm的细粒碎屑、黏土、有机质等组成, 具页状或薄片状层理分布、易碎裂的一类沉积岩.岩石物理学研究表明:页岩表现出很强的各向异性特征[3-6].鉴于我国页岩气地震勘探目前主要以低成本的常规地震声波资料为主, 在页岩气勘探地震数据处理中, 如何针对页岩强各向异性而采取有别于常规油气勘探的处理流程和技术是目前需要解决的关键问题.通过数值模拟页岩的声波响应特征, 结合实际页岩地震反射资料分析, 对建立适合于页岩气勘探的地震数据处理流程和技术具有重要的实际意义.

波传播模拟强各向异性页岩要求数值方法具有较高的模拟精度, 特别是对易引起频散的高频波(>100Hz)数值模拟.有限差分法由于数值实施简单灵活而得到快速发展和广泛应用, 特别是基于弹性波动方程的交错网格有限差分数值实施被广泛用于研究裂缝介质的各向异性特征[7-13].Alkhalifah[14]从横向各向同性(TI)介质频散关系出发, 将沿对称轴方向的横波速度置为零获得垂直TI(VTI)介质的伪声波方程.直接将横波速度置为零会导致震源处产生横波, 在资料处理过程中被视为噪声.该问题是因为Thomsen各向异性参数[15] ε和δ在震源附近不等造成的.因此, 该问题可以较为简单的解决, 即在震源附近通过光滑使二者相等, 但这会导致相位有微小的差异而并不会影响对波传播的理解[16].本文中, 首先介绍了VTI介质的声波方程, 该方程为将横波速度置为零的伪声波方程.然后, 将任意偶数2M阶的展开系数用于求解VTI介质声波方程, 并讨论了稳定性条件和吸收边界条件.最后, 以模型为例计算了地面地震剖面和垂直地震剖面(VSP), 并分析了页岩的各向异性对地震波的运动学和动力学特征的影响.

2 VTI介质的声波方程

三维VTI介质常密度声波方程[17]表示如下:

(1)

其中, Q和P分别为地震波场的水平应力分量和垂直应力分量, VP为声波速度, ε和δ为Thomsen各向异性参数[15].二维VTI介质常密度声波方程表示如下:

(2)

当各向异性参数δ为常数时, 方程(1)-(2)可退化成Zhou等[18]提出的VTI介质声波方程.为了书写简单, 令, 则方程(2)可重新写作

(3)

各向异性参数ε和δ的计算公式[15]如下

(4)

其中, Ci, j(i, j=1, 2, 3, 4)为弹性刚度矩阵的系数.

2.1 任意偶数2M阶有限差分解

本文中, 时间导数采用如下二阶精度的差分格式:

(5)

空间导数和采取任意2M阶的精度, 具体差分格式如下

(6)

其中, b0=和bm=Bm为常系数, 不同阶数精度的常系数bm的值可参见文献[19-20].本文采取空间差分计算精度取12阶(M=6), 常系数分别为b0=-5369/1800, b1=12/7, b2=-15/56, b3=10/189, b4=-1/112, b5=2/1925和b6=-1/16632.同样地, 地震波场垂直应力分量的空间差分格式如下

(7)

其中, (i, k, n)分别表示空间x, z方向和时间的离散网格, (Δx, Δz, Δt)分别为空间x, z方向和时间的采样间隔.

将方程(6)~(7)代入方程(5), 整理得

(8)

和

(9)

2.2 稳定性条件

为了保证有限差分格式(见方程(8)和(9))的计算稳定, 空间最大采样间隔Δdmax=max(Δx, Δz)须满足如下条件:

(10)

其中, λmin为地震波场的最小波长, Vmin为介质的最小速度, fmax为地震波场的最大频率.方程(10)说明在一个地震波长内要至少有5个空间采样点.同样, 时间采样间隔要足够小, 才能保证方程(8)和(9)的计算稳定.Jin等[21]从各向同性介质的稳定性条件出发, 给出了VTI介质的稳定性条件为

(11)

其中,

(12)

式中, Δdmin=min(Δx, Δz)为空间最小采样间隔, Vmax为介质的最大速度, b0和bm(m=1, 2, …, M)为有限差分系数, εmax和δmax为各向异性参数的最大值.本文采用了12阶空间差分的计算精度, 将有限差分系数代入到方程(12)中, 有:

(13)

2.3 吸收边界条件

由于地下介质被视为无限半空间, 而计算区域是有限的.因此, 在人为截取的有限空间内求解波动方程会导致很强的人工边界反射.目前, 完美匹配层(PML)[22-23]类的吸收边界被认为是效果最好的吸收边界, 然而其数值程序实现过程较为繁琐.Cerjan等[24]给出了一种十分简单的衰减吸收边界条件, 即在边界区域施加如下衰减函数:

(14)

其中, N为衰减区域的网格点数, i为衰减区域内的网格点序号, a为衰减系数.Cerjan等建议N=20和a=0.015, 但是吸收边界条件不能有效的吸收人工边界反射.Bording [25]指出Cerjan吸收边界效果之所以不好是因为参数未进行优化导致的, 并给出了优化后的参数为N=45和a=0.0053.图 1中比较了均匀介质情况下两种吸收边界的吸收效果.在使用优化后的参数后, 人工边界反射的吸收效果得到明显提高.因此, 本文采用Bording吸收边界条件.

图 1 Cerjan吸收边界(a)和Bording吸收边界(b)吸收效果比较 Fig. 1 Comparison df absorbing boundary conditions (a) Cerjan type; (b) Bording type.
3 数值算例 3.1 VTI介质的均匀模型

本节中, 我们计算了一个理论背景均匀模型的地震响应.该模型的网格点为201×201, 空间采样间隔为10m×10m.震源函数为Ricker子波, 主频为20Hz, 时间采样间隔为1ms, 震源位于模型中央(1000m, 1000m).图 2给出了不同情况下0.25s的波场快照.图 2a为声波波动方程(即各向异性参数ε=0且δ=0)计算的波场快照, 图 2(b-d)为VTI介质声波波动方程计算得到的波场快照, 对应的各向异性参数(ε, δ)分别为(0.25, 0), (0.25, 0)和(0.25, 0.25).图 2b和图 2c采用了相同的各向异性参数, 但图 2b未做震源校正而产生了剪切波的数值噪声, 图 2c做了震源校正(即在震源附近通过光滑平均, 令ε=δ)而避免了数值噪声.尽管图 2c在震源附近对各向异性参数做了近似, 但并未改变或很小的改变准P的相位.图 2d中, 由于各向异性参数相等而未出现剪切波的数值噪声.需要注意的是, 图 2采用了相同的增益, 而各向异性声波方程计算得到的波场快照的振幅明显不同于各向同性声波方程计算的波场快照的振幅.

图 2 不同各向异性参数计算的在t=0. 25 s时的波场快照 (a) ε=0, δ=0.此时, VTI介质声波方程等同于各向同性声波方程; (b)ε=0.25, δ=0, 未做震源校正; (c)ε=0.25, δ=0, 带震源校正; (d)ε=0.25, δ0.25.四幅图采用相同的增益控制. Fig. 2 Snapshots at t=0. 25 s with different anisotropy parameters (a) ε=0, δ=0. In this case, the VTI acoustic wave equation is equivalent to the isotropic acoustic wave equation; (b)ε=0.25, δ=0 without source correction; (c)ε=0.25, δ=0 with source correction; (d)ε=0.25, δ=0.25. These four figures are plotted at the same scale.
3.2 多层页岩模型

页岩广泛的分布在我国的各个盆地[1-2], 其埋深大约为1500~2500 m.页岩的声波速度为4000~5200m/s, 密度为1.9~2.4g/cm3(参见文献[26]).图 3为一个多层的页岩模型, 水平方向和深度方向的空间采样间隔均为8m, 空间采样点分别为1001×626.模型的主要参数如表 1所示, 其中模型的第三层和第四层分别代表页岩层和砂岩层.震源函数为Ricker子波, 主频为15 Hz, 时间采样间隔为1ms, 计算时间长度为4s.

表 1 层状页岩模型的介质参数 Table 1 Media parameters for layered shale model
图 3 层状页岩模型的介质参数 具体介质参数见表 1，第3层和第4层分别为页岩和砂岩. Fig. 3 The multi-layered shale model The media parameters are shown in Table 1, where the third line and the fourth line give the parameters of shale and sandstone, respectively.

图 4为在地表接收的地震记录, 震源位于模型的中央(4000, 0)m处.图 4a为声波方程计算的水平地震记录, 图 4b为VTI介质声波方程计算的水平地震记录.比较图 4a和4b可以发现, 各向异性对地震波的相位和振幅影响是十分明显的, 如图中白色箭头标注.以各向同性介质情况下的地震记录为参考, 利用Kristek等[27]和Lan和Zhang [28]的方法定量分析了垂直横向各向异性介质对地震波传播的影响, 如图 5所示.由于页岩表现出强的各向异性, 图 5中可以清晰的发现页岩各向异性对地震波的相位和振幅的影响十分明显.

图 4 声波方程(a)和VTI声波方程(b)计算的地表接受的地面地震记录 震源位于(4000, 0)m，震源主频为15 Hz.两图采用同样的绘图尺度. Horizontal seismic records calculated by (a) acoustic wave equation and (b) VTI acoustic wave equation The source is located at (4000, 0)m. The peak frequency of source is 15 Hz. These two figures are plotted at the same scale.
图 5 垂直横向各向同性介质对水平地震记录的振幅(a)和相位(b)的影响（模型为图 3所示) Fig. 5 The influence of VTI media to the amplitudes (a) and the phases (b) of horizontal seismograms (The model s shown in Fig. 3)

图 6为声波方程和VTI声波方程计算得到的垂直地震剖面, 震源位于4000m, 接收位置为X=5000 m.通过VSP地震记录, 可以更加清晰的了解地震波传播通过页岩层时的响应.当地震波通过页岩时, 地震波的振幅明显增强(图 6中A处)或衰减(图 6中B处).由于页岩各向异性的存在, 其对相位的影响则更加明显.此外, 图 6b中, 当地震波传播至页岩层时, 可以观察到有轻微的频散现象出现, 这可能是因为上下层的速度对比很强, 速度对比约为0.38.图 7以VSP地震记录定量分析了VTI介质对地震波传播的影响.图 7中可以明显看到各向异性对地震波走时和振幅的影响, 尤其是页岩层(2000~3000 m)对振幅的影响十分明显.因此, 对于针对页岩气的地震资料要十分谨慎地考虑各向异性的影响, 若采用各向同性的技术去处理页岩气的地震资料将会得到不准确的偏移结果.

图 6 声波方程(a)和VTI声波方程(b)计算的垂直地震剖面（VSP) 震源位于4000 m，接收位置为X=5000 m，震源主频为15 Hz.两图采用同样的绘图尺度. Fig. 6 Vertical seismic profiles calculated by (a) acoustic wave equation and (b) VTI acoustic wave equation The source is located at (4000, 0) m. The receivers is located at X=5000 m. The peak frequency of source is 15 Hz. These two figures are plotted at the same scale.
图 7 垂直横向各向同性介质对垂直地震记录的振幅(a)和相位(b)的影响（模型为图 2所示) Fig. 7 The influence of VTI media to the amplitudes and the phases of vertical seismic profiles (The model is shown in Fig.2)
3.3 页岩气勘探的震源主频

由于页岩埋深较浅, 且与砂岩以交互形式存在.常规地震勘探的主频大约为10~20 Hz, 难以满足埋深较浅的页岩气勘探对分辨率的需求.而页岩气开发过程中, 用于监测裂缝的微地震资料的主频已经达到数百赫兹甚至上千赫兹.此外, 页岩十分致密, 页岩气的开发需要水平井压裂, 清楚的知道页岩的顶底边界对于后期水平井开发是十分重要的.因此, 中浅层的页岩气勘探需要提高地震勘探的主频.煤层气及煤层勘探的地震主频目前为80~150Hz[29], 因此, 可以利用地震正演的方法进行初步的尝试, 将页岩气勘探的主频提高至100Hz.

本节建立了一个多层页岩-砂岩交互模型(如图 8所示), 模型物理参数如表 2所示.含气页岩的埋深较浅, 最浅埋深为800m.模型的水平和深度方向的网格点数分别为3001和2001, 采样间隔均为1m, 震源函数为Ricker子波, 主频为100 Hz, 频带范围为300 Hz.时间采样间隔为0.1ms, 时间采样长度为1.5s.

图 8 层状页岩-砂岩模型(该模型包含有一个含气页岩薄层) Fig. 8 The layered shale-sandstone model, which contains a thin shale gas layer
表 2 层状页岩砂岩模型的介质参数 Table 2 Media parameters for layered shale-sand-stone model

与常规油气相比, 页岩气的埋深较浅, 因此提高地震勘探的主频是可行的.为了说明高频勘探的必要性, 本文首先给出了主频为15 Hz的地面地震记录和垂直地震记录, 如图 9所示.震源函数为Ricker子波, 位于(1000, 0)m, 垂直地震记录的接收位置为X=1500m.图 9中, 由于地震主频较低导致无法分辨来自薄层顶底界面的反射, 且较深的反射同相轴振幅较弱.图 9b中发现, 来自薄层顶底界面的反射和来自层3的反射重叠在一起, 这主要是因为地震主频低引起的.

图 9 震源位于(1000, 0)m处的地面地震记录(a)和位于1500 m处的垂直地震剖面(b) 震源函数为Ricker子波，主频为15 Hz Fig. 9 The surface seismic records (a) and vertical seismic profile (b) at X=1500 m. The source is located at (1000, 0) m The source function is Ricker wavelet with peak frequency of 15 Hz.

图 10为震源位于(1000, 0)m和(2000, 0)m处的地面地震记录, 震源主频为100 Hz.图 10中发现, 可以明显分辨出薄层顶底界面的反射, 有利于地震偏移成像时分辨出顶底界面, 为后续的页岩气水平井压裂开发奠定了基础.图 11和图 12分别为震源位于(1000, 0)m和(2000, 0)m处在X=1500m记录的垂直地震剖面.图 11和12中, 由于地震主频的提高, 来自各个界面的反射以及界面之间的多次反射可以清晰的观察到.较为明显的现象是, 由于薄层的存在, 地震波“陷入”在薄层内形成“导波”(waveguides).

图 10 不同位置的地面地震记录 (a) 1000 m; (b) 2000 m.震源函数为Ricker子波，主频为100 Hz. Fig. 10 The surface seismograms fired at different locations (a) 1000 m; (b) 2000 m. The source function is Ricker wavelet with the dominant frequency of 100 Hz.
图 11 位于X=1500 m处的垂直地震记录 震源位于（1000, 0)m.震源函数为Ricker子波，主频为100Hz.图（b)为图（a)的部分窗口放大. Fig. 11 The vertical seismic profile at X=1500 The source location is (1000, 0) m. The source function is Ricker wavelet with the dominant frequency of 100 Hz. Figure 11b is a zoom window of Figure 11a.
图 12 位于X=1500 m处的垂直地震记录 震源位于(2000, 0)m.震源函数为Ricker子波，主频为100 Hz.图(b)为图(a)的部分窗口放大. Fig. 12 The vertical seismic profile at X=1500 The source location is (2000, 0) m. The source function is Ricker wavelet with the dominant frequency of 100 Hz. Figure 12b is a zoom window of Figure 12a.
4 结论

本文介绍了任意偶数2M阶有限差分求解VTI介质声波方程的离散格式, 其计算稳定, 几乎无频散现象发生.然后讨论了稳定性条件和吸收边界条件, 其数值实施简单有效.针对多层页岩模型, 比较了声波方程和VTI介质声波方程的地面地震记录和垂直地震剖面.结果发现, 由于页岩表现出较强的各向异性, 而导致地震波的走时和振幅较声波方程计算的有明显差异.这说明在处理有关页岩气的地震资料时, 必须充分考虑各向异性的影响, 否则将得不到准确的偏移结果.由于含气页岩层顶底界面的准确位置对于后期的水平井压裂和开发十分重要, 且页岩气通常埋深较浅, 提高勘探主频可以明显提高地震资料的分辨率, 有助于提高岩性反演和储层反演的精度.因此, 提高地震勘探主频对于页岩气勘探是十分必要的.

参考文献
[1] 张金川, 徐波, 聂海宽, 等. 中国页岩气资源勘探潜力. 天然气工业 , 2008, 28(6): 136–141. Zhang J C, Xu B, Nie H K, et al. Exploration potential of shale gas resources in China. Natural Gas Industry (in Chinese) , 2008, 28(6): 136-141.
[2] 张金川, 姜生玲, 唐玄, 等. 我国页岩气富集类型及资源特点. 天然气工业 , 2009, 29(12): 109–114. Zhang J C, Jiang S L, Tang X, et al. Accumulation types and resources characteristics of shale gas in China. Natural Gas Industry (in Chinese) , 2009, 29(12): 109-114.
[3] Prasad M, Pal-Bathija A, Johnston M, et al. Rock physics of the unconventional. The Leading Edge , 2009, 28(1): 34-38. DOI:10.1190/1.3064144
[4] Delle P C, Dewhurst D, Sarout J. TI or not TI? Stress effects on shale anisotropy. 14th International Workshop on Seismic Anisotropy, 2010:41-42.
[5] 史謌, 邓继新. 地层条件下泥、页岩衰减各向异性研究. 中国科学(D辑:地球科学) , 2005, 48(11): 1882–1890. Shi G, Deng J X. The attenuation anisotropy of mudstones and shales in subsurface formations. Science in China Series D:Earth Sciences (in Chinese) , 2005, 48(11): 1882-1890. DOI:10.1360/03yd0426
[6] Franquet J, Patterson D, Moos D, et al. Advanced dipole borehole acoustic processing-Rock Physics and Geomechanics Applications. 81st Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 2011:1876-1871.
[7] 杨顶辉, 滕吉文. 各向异性介质中三分量地震记录的FCT有限差分模拟. 石油地球物理勘探 , 1997, 32(2): 181–190. Yang D H, Teng J W. FCT finite difference modeling of three-component seismic records in anisotropic medium. Oil Geophysical Prospecting (in Chinese) , 1997, 32(2): 181-190.
[8] Virieux J. P-SV wave propagation in heterogeneous media:Velocity-stress finite-difference method. Geophysics , 1986, 51(4): 889-901. DOI:10.1190/1.1442147
[9] Saenger E H, Gold N, Shapiro S A. Modeling the propagation of elastic waves using a modified finite-difference grid. Wave Motion , 2000, 31(1): 77-92. DOI:10.1016/S0165-2125(99)00023-2
[10] McGarry R, Pasalic D, Ong C. Anisotropic elastic modeling on a Lebedev grid:Dispersion reduction and grid decoupling. 81st Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 2011:2829-2833.
[11] 祝贺君, 张伟, 陈晓非. 二维各向异性介质中地震波场的高阶同位网格有限差分模拟. 地球物理学报 , 2009, 52(6): 1536–1546. Zhu H J, Zhang W, Chen X F. Two-dimensional seismic wave simulation in anisotropic media by non-staggered finite difference method. Chinese J. Geophys. (in Chinese) , 2009, 52(6): 1536-1546.
[12] 兰海强, 刘佳, 白志明. VTI介质起伏地表地震波场模拟. 地球物理学报 , 2011, 54(8): 2072–2084. Lan H Q, Liu J, Bai Z M. Wave-filed simulation in VTI media with irregular free surface. Chinese J. Geophys. (in Chinese) , 2011, 54(8): 2072-2084.
[13] 李桂花, 冯建国, 朱光明. 黏弹性VTI介质频率空间域准P波正演模拟. 地球物理学报 , 2011, 54(1): 200–207. Li G H, Feng J G, Zhu G M. Quasi-P wave forward modeling in viscoelastic VTI meida in frequency-space domain. Chinese J. Geophys. (in Chinese) , 2011, 54(1): 200-207.
[14] Alkhalifah T. Acoustic approximations for processing in transversely isotropic media. Geophysics , 1998, 63(2): 623-631. DOI:10.1190/1.1444361
[15] Thomsen L. Weak elastic anisotropy. Geophysics , 1986, 51(10): 1954-1966. DOI:10.1190/1.1442051
[16] Grechka V, Zhang L B, Rector J W. Shear waves in acoustic anisotropic media. Geophysics , 2004, 69(2): 576-582. DOI:10.1190/1.1707077
[17] Duveneck E, Milcik P, Bakker P M, et al. Acoustic VTI wave equations and their application for anisotropic reverse-time migration. 78th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 2008:2186-2190.
[18] Zhou H, Zhang G, Bloor R. An anisotropic acoustic wave equation for VTI media. 68th EAGE Conference and Exhibition, 2006:194-197.
[19] Liu Y, Sen M K. A practical implicit finite-difference method:examples from seismic modelling. Journal of Geophysics and Engineering , 2009, 6(3): 231-249. DOI:10.1088/1742-2132/6/3/003
[20] Pei Z L, Fu L Y, Sun W J, et al. Anisotropic finite-difference algorithm for modeling elastic wave propagation in fractured coalbeds. Geophysics , 2012, 77(1): C13-C26. DOI:10.1190/geo2010-0240.1
[21] Jin S W, Jiang F, Ren Y Q, et al. Comparison of isotropic, VTI and TTI reverse time migration:an experiment on BP anisotropic benchmark dataset. 80th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 2010:3198-3203.
[22] Berenger J P. A perfectly matched layer for the absorption of electromagnetic waves. Journal of Computational Physics , 1994, 114(2): 185-200. DOI:10.1006/jcph.1994.1159
[23] Komatitsch D, Martin R. An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation. Geophysics , 2007, 72(5): SM155-SM167. DOI:10.1190/1.2757586
[24] Cerjan C, Kosloff D, Kosloff R, et al. A nonreflecting boundary condition for discrete acoustic and elastic wave equations. Geophysics , 1985, 50(4): 705-708. DOI:10.1190/1.1441945
[25] Bording R P. Finite difference modeling-nearly optimal sponge boundary conditions. 74th Ann. Internat Mtg., Soc. Expi. Geophys. Expanded Abstracts, 2004:1921-1924.
[26] 陈颙, 黄庭芳, 刘恩儒. 岩石物理学. 合肥: 中国科学技术大学出版社, 2009 . Chen Y, Huang T F, Liu E R. Rock Physics (in Chinese). Hefei: Press of University of Science and Technology of China, 2009 .
[27] Kristek J, Moczo P, Archuleta R. Efficient methods to simulate planar free surface in the 3D 4th-order staggered-grid finite-difference schemes.. Stud.Geophys. Geod. , 2002, 46(2): 355-381. DOI:10.1023/A:1019866422821
[28] Lan H Q, Zhang Z J. Comparative study of the free-surface boundary condition in two-dimensional finite-difference elastic wave field simulation. Journal of Geophysics and Engineering , 2011, 8(2): 275-286. DOI:10.1088/1742-2132/8/2/012
[29] Sun W J, Zhou B Z, Hatherly P, et al. Seismic wave propagation through surface basalts-implications for coal seismic surveys. Exploration Geophysics , 2010, 41(1): 1-8. DOI:10.1071/EG09015