Comparative study among implementations of several free-surface boundaries with perfectly matched layer conditions
-
摘要: 传统的完全匹配层技术是一种能够较为有效地消除边界反射的边界条件,但是当表层为泊松比较高的自由表面时,该技术可能会产生不稳定的现象.针对传统的完全匹配层技术固有的不稳定和掠射情况下吸收效果不佳等缺陷,发展了多轴完全匹配层、卷积完全匹配层以及将两者结合的多轴卷积完全匹配层等3种边界条件.本文介绍了水平自由表面的不同处理方法以及传统、多轴、卷积和多轴卷积等4种完全匹配层条件的原理,通过二维半无限空间模型的交错网格有限差分正演模拟对比,分析了几种自由边界实施方法在这几种完全匹配层条件下的稳定性,并通过提取单道波形与解析解进行对比,定性分析了水平自由表面几种不同处理方法的准确性以及各自的适用条件. 结果表明,泊松比和水平自由表面实施方法对波场模拟效果及其稳定性有重要影响.Abstract: Classical perfectly matched layer (PML) absorbing boundary is certified an efficient method to suppress spurious edge reflections in seismic modeling. However, when modeling Rayleigh waves with the existence of free surface, the classical PML algorithm may become unstable in the case that Poisson’s ratio of the medium is high. On the basis of defects of classical PML, such as its inherent instability and the effects of grazing incidence, multiaxial perfectly matched layer, convolutional perfectly matched layer and multiaxial convolutional perfectly matched layer are developed. Six implementations of planar free-surface boundary conditions and principles of four perfectly matched layer methods are introduced, and forward modeling of 2D half-infinite model using staggered grid finite-difference method are given to study the stability of different methods. Through the comparison of waveform between numerical solution and analytical solution, we can illustrate the accuracy and applicability of these methods in this paper. The results show that Poisson’s ratio and implementations of planar free-surface have an important influence on the effects and stability of wavefield simulation.
-
Keywords:
- PML /
- planar free-surface /
- Poisson’s ratio /
- stability /
- accuracy
-
引言
云南地区位于青藏高原东南缘,晚元古代以来,由于受到古大洋向其北部和东部地台的长期挤压和俯冲,以及各时期主应力方向不断地变化,形成了该地区地壳缩短、 高原隆升和块体滑移等复杂的构造格局,是中国大陆构造运动最强烈的地区之一(白志明,王椿镛,2004; 张智等,2006). 据统计,该地区20世纪发生M≥5.0破坏性地震333次,其中M≥7.0地震13次. 就地震频度而言,M5.0—5.9地震平均每年发生3次,M6.0—6.9地震平均每两年发生1次,M7.0—7.9地震平均每8年发生1次(毛玉平,韩新民,2003). 此外,由于印度板块与欧亚板块的相互作用,云南地区也是青藏高原“物质东流”的挤出通道(曾融生,孙为国,1992). 因此,深入研究该地区的地壳结构,对认识青藏高原构造动力学、 川滇地区地震活动性以及提高云南地区测震台网的定位精度等都具有积极意义.
云南地区一直是众多学者研究的重点地区,该地区曾开展过大量的人工地震测深工作,其中包括“滇深82”、 “滇深86-87”、 “腾深99”以及“昆深05”等(胡鸿翔等,1986; 阚荣举,林中洋,1986; 熊绍柏等,1986; 林中洋等,1993; 白志明,王椿镛, 2003,2004; Wang,Huangfu,2004; 张中杰等, 2005a,b; 张智等,2006). 接收函数、 体波和面波层析成像等方法也用于该地区的壳幔结构研究(吴建平等,2001; Huang et al,2002; 王椿镛等,2002; 胡家富等,2003; 何正勤等, 2004a,b; 张智等,2008; 张晓曼等,2011). 这些成果加深了我们对云南地区壳幔结构和构造动力学的认识,但由于受人工地震测深剖面位置、 台站分布以及反演方法非唯一性等影响,云南地区壳幔速度结构仍然存在诸多问题,仍需进一步研究. 例如: Royden等(1997)的中下地壳流模型能够较好地解释青藏高原的演化和地表形变特征,但低速层分布范围及形态仍存在诸多争议; 云南地区强震发生通常与规模较大的断裂带有关,但强震的孕震环境仍不完全清楚,如红河断裂和小江断裂带的速度结构与强震之间的关系.
背景噪声层析成像是近年来快速发展的一种成像方法. 数值试验和理论研究表明,在一个均匀散射场中,任意两点间的格林函数可以从这两点的位移互相关函数中提取(Lobkis,Weaver,2001; Wapenaar,2004). Shapiro等(2005)通过对地震台站长时间的地震噪声记录进行互相关计算,提取出台站间的经验格林函数,并且得到了美国南加州7.5 s和15 s周期的瑞雷面波群速度分布图. 此后,背景噪声成像技术得到飞速发展,无论是在大尺度区域(如美国、 欧洲、 澳大利亚、 南非、 中国大陆),还是在小尺度区域(如美国加州、 韩国、 新西兰)均能得到高信噪比的经验格林函数和高分辨率的面波层析成像结果(Bensen et al, 2007a,2009; Cho et al,2007; Lin et al,2007; Yang et al, 2007,2008; Zheng et al, 2008,2011; Saygin,Kennett,2010; Zhou et al,2012). 2006年以来,国内逐渐开始使用该方法. 华北、 四川、 西藏东南缘和黑龙江等地均已取得很多研究成果(Yao et al, 2006,2008; 房立华等,2009; Li et al,2009; 李昱等,2010; 高东辉等,2011). 简单来说,背景噪声层析成像方法被广泛采用的主要原因有: ① 该方法使用的数据不是天然地震资料,因此不依赖震源位置分布; ② 天然地震成像中震中位置和发震时刻存在误差,而背景噪声成像方法使用的是精确的台站经纬度,横向分辨率取决于台站之间的间距,只要台网内台站密集就可以得到高分辨率成像结果; ③ 背景噪声可以在不同时间段多次重复记录并计算其频散曲线,因此可以通过经验格林函数相互叠加来提高信噪比和对结果进行不确定性分析; ④ 背景噪声成像方法提取出的频散曲线周期主要集中在6—70 s,可以反映地壳结构方面的信息,而传统面波频散曲线的周期以长周期为主,远震体波成像方法在地壳中成像分辨率又比较差(Yang et al,2008). 综上所述,背景噪声层析成像方法是对传统研究方法很好的补充.
Yao等(2006,2008)利用中美合作的流动地震观测台网所记录的噪声数据及面波数据反演了青藏高原东南缘地壳上地幔三维速度结构,但是由于受台站覆盖区域和空间分辨率的限制,云南地区地壳速度结构的分辨率仍然不高. 本文利用云南地区固定台网的数据资料,从背景噪声互相关函数中提取瑞雷波经验格林函数,并采用面波层析成像技术建立云南地区的相速度分布图,为下一步反演该地区的三维剪切波速度结构奠定基础.
1. 资料及方法
本文收集了2007年11月—2009年10月云南区域台网46个固定台站(图 1)的连续背景噪声记录资料. 云南区域台网台站观测系统由EDAS2-24 IP数据采集器、 宽频带地震计(频宽60 s—50 Hz)、 IP协议转换器和电源构成. 整个台网数据字长为24位,采样率为100 Hz,系统动态范围可达140 dB,台网各子台在1—20 Hz频带内实际观测动态范围为95 dB左右,可测量的最小地动信号为10-8 m/s. 本研究中,我们使用地震台站垂直分量数据. 数据处理流程参照Bensen等(2007b)的方法,该方法主要由5部分组成: ① 单台数据预处理,主要包括1 Hz重采样、 去除直流分量及线性趋势、 5—50 s的带通滤波、 时域归一化和谱白化(时域归一化使用的是移动均值法; 谱白化即在频谱里将幅值拉平,这样可以拓宽地震背景噪声的频谱,抑制某一单频信号的干扰); ② 互相关计算和叠加; ③ 频散曲线的测量; ④ 质量控制和误差分析; ⑤ 面波层析成像. 在数据处理过程中,为了保证得到比较可靠的频散曲线,我们需要满足以下3个方面的条件: ① 选取两年数据进行叠加; ② 每个中心周期经验格林函数的信噪比SNR>20,不足的予以剔除; ③ 两个台站间的距离必须大于3个波长,一个简单的确定法则就是τmax=Δ/12,τmax为最大测量周期,Δ为两个台站间的距离(Bensen et al,2007b). 云南区域台网覆盖范围约为600 km×600 km,所以只选择40 s以下的周期进行研究. 图 2为昆明台(HLT)与其余各台站路径上的面波经验格林函数.
本文采用自适应时频分析技术计算群速度和相速度频散曲线(Levshin,Ritzwoller,2001; Bensen et al,2007b; Lin et al,2008),该方法给定一个参考相速度来避免相位模糊. 一般来说,垂向背景噪声提取出的互相关函数的相慢度可以表示为
图 1 云南地区区域构造略图(阚荣举,韩源,1992)及台站和地震分布图F1: 怒江断裂; F2: 北澜沧江—昌宁—双江断裂; F3: 金沙江断裂; F4: 红河断裂; F5: 剑川断裂;F6: 程海断裂; F7: 元谋—绿汁江断裂; F8: 普渡河断裂; F9: 小江断裂带; F10: 畹町断裂; F11: 南汀河断裂; F12: 弥勒断裂; F13: 南澜沧江断裂. ① 滇缅泰板块; ② 印支板块; ③ 扬子板块; ④ 华南亚板块. Tc: 腾冲块体; Ba: 保山块体; Ls: 兰坪—思茅弧后盆地Figure 1. Regional tectonic map of Yunnan area(after Kan,Han,1992) as well as locations of stations and earthquakesF1: Nujiang fault; F2: Northern Lancangjiang- -Changning- -Shuangjiang fault; F3: Jinshajiang fault; F4: Honghe fault; F5: Jianchuan fault; F6: Chenghai fault; F7: Yuanmou- -Lüzhijiang fault; F8: Puduhe fault; F9: Xiaojiang fault zone; F10: W and ing fault; F11: Nantinghe fault; F12: Mile fault; F13: Southern Lancangjiang fault. ① Yunnan- -Myanmar- -Thail and block; ② Indo-China block; ③ Yangtze block; ④ South China block. Tc: Tengchong block; Ba: Baoshan block; Ls: Lanping- -Simao back arc basin. The triangle in the right panel indicates the station and the circle indicates the earthquake M≥5.0 occurred since the year AD 1500式中: SC和SU分别为相慢度和群慢度; ω为频率; Δ为台站间距; tU为Δ/U,U为群速度,(tU)为相应频率对应的相位; N=0,±1,±2,…. 为了确定N值,需要参考相速度模型. 对于长周期面波,N值的改变对测量结果有很大影响,因而很容易确定. 实际操作中,我们首先根据频散曲线得到群速度,然后计算长周期面波的相速度,依次计算较短周期的相速度,并且保证相速度随周期的连续性. 计算时我们参考了Yao等(2006,2008)和何正勤等(2004b)的研究结果,并选取相应的初始参考模型(表 1). 图 3给出了由台站对FUN-TNC的背景噪声数据得到的面波经验格林函数和由自适应时频分析技术求取的瑞雷波群速度频散曲线.
表 1 相速度值的初始参考模型Table 1. The initial reference phase velocity model at different periods图 3 台站对FUN-TNC的经验格林函数(a)及群速度频散曲线测量示意图(b)(台站FUN和TNC的位置分布见图 1)Figure 3. Example of the empirical Green’s functions and dispersion analysis of the station pair FUN-TNC(locations of stations FUN and TNC,see Fig. 1)(a)The empirical Green’s functions of station pair FUN-TNC;(b)Group velocity dispersion curve measurements2. 面波层析成像
本文采用Ditmar和Yanovskaya(1987)以及Yanovskaya和Ditmar(1990)的面波层析成像反演方法分别对8—40 s范围内15个中心周期进行反演,得到各个周期的相速度分布图. 该面波层析成像反演方法是传统Backus-Gilbert一维方法在面波二维反演情况下的推广(Backus,Gilbert,1968),是面波层析成像中广泛应用的方法之一(房立华等,2009; Guo et al,2009; 潘佳铁等,2011). 模型的基函数使用群到时的积分形式表示,不需要设定初始参数和约束条件. 假设相速度的分布用函数Ce(θ,φ)表示,θ和φ分别表示经度和纬度,通过使目标函数ψ(Ce)最小化来获得每个周期的相速度分布,即
其中
式中: r=r(θ,φ)是位置矢量; C0是初始模型速度,一般取该周期所有路径上的平均相速度; ti是第i条路径的观测走时; ti0是根据初始模型计算的走时; l0i是第i条路径的长度; s是参与反演的路径; α是正则化参数. 正则化参数控制着反演结果的光滑程度,α越大,反演结果越光滑,分辨率越低; 反之,则反演结果的分辨率越高,误差也越大,反演结果中会经常出现高低速交替的情形. 在反演过程中,α分别选取0.1,0.15,0.2,0.25,0.3来尝试确定最终值. 当α选取为0.2时,所得的结果比较光滑,且误差较小,所以确定为最终值.
由于观测数量有限,层析成像的结果也存在非唯一性. 面波层析成像的分辨率主要取决于射线密度和方向分布. 本文通过计算每个格点空间平均分辨率核函数来评估层析成像结果的分辨率(Backus,Gilbert,1968; Yanovskaya,Ditmar,1990; Guo et al,2009). 研究区内任意一点x0(θ0,φ0)的分辨率核函数R(θ,φ; x0)可以表示为
式中: δĈ(x0)/C表示点x0(θ0,φ0)的速度变化量; δC(θ,φ)/C表示整个区域的真速度变化量; Ω表示包含所有射线的区域. 图 4给出了云南地区射线覆盖分布图,周期为8 s,18 s和40 s的射线数目分别为606,901和278. 图 5给出了周期为8 s,18 s和40 s的分辨率半径分布图. 分辨率评价表明,在研究区域中部的大部分地区分辩率可达到40—60 km,其边界附近的分辨率约为80 km. 受分辨率限制,研究区域边缘的速度异常可能并不准确.
3. 相速度分布特征
相速度反映约为1/3个波长范围内的S波平均速度,不同周期的相速度反映不同深度范围内的构造. 通过相速度图就可以了解某一深度范围内的速度结构横向变化特征,相速度随周期从短到长的变化反映了介质从浅到深的速度结构变化特征.
短周期(8—15 s)相速度分布图反映大约为上地壳的平均速度变化情况. 图 6a为8 s周期的相速度分布图. 该周期瑞雷波相速度横向变化与区域地质构造(沉积层厚度、 断层分布等)关系密切. 北北西向展布的兰坪—思茅弧后盆地,在相速度分布图中表现为低速异常,这反映出此凹陷中巨厚中生代盖层的影响. 红河断裂(北段)东侧的楚雄盆地由于巨厚的中生代沉积,呈现出低速异常. 小江断裂带是构造碎裂化比较显著的断裂带,介质的完整性比较差,显示为低速异常. 而弥勒断裂南侧存在非常明显的速度分界面,该断裂南侧的滇东南区域为稳定抬升区,呈现出高速异常的特性. 腾冲南侧地区呈现出低速异常,可能与区域内断裂分布、 火山构造和盆地分布对低速异常的控制作用有关. 楼海等(2002)及Wang和Huangfu(2004)分析了人工地震测深资料,提出腾冲以南上地壳的低速异常可归因于上地幔源区的热活动和岩浆分异作用. 该结果与Yao等(2006)工作相比较,相互重叠地区的速度异常形态大体一致. 总体上看,相速度的异常与沉积层厚度分布和断裂分布有一定的关系.
中周期(16—25 s)相速度分布图反映大约为中下地壳深度范围内平均速度的变化情况. 图 6b为18 s周期的相速度分布图. 从图中可以看出,川滇菱形块体内主要为低速异常. 从分布形态上看,与Yao等(2006)研究结果相似. 该低速异常可能是川滇菱形块体内中下地壳存在低速层所致,但该低速层是否连续存在还有待商榷. 红河断裂西侧和弥勒断裂南侧为明显的速度变化分界面. 滇东南区域仍为高速区,说明该区域中下地壳构造比较稳定. 腾冲块体和保山块体低速异常逐渐消失. 印支块体南部思茅—普洱低速异常逐渐消失,可能是思茅盆地内沉积层的影响消失所致. Yao等(2006)利用噪声数据和天然地震数据分别对该地区面波成像,其中周期(20 s)相速度分布图显示,研究区西南边缘(约24°N,100°E附近)呈现不同结果. 该位置在本研究中位于分辨率较高地区,本文支持Yao等(2006)噪声数据的成像结果,该地区附近确实出现高速异常.
长周期(26—40 s)相速度分布图反映大约为下地壳至上地幔范围内平均速度的变化情况,而且与地壳的厚度密切相关. 图 6c为40 s周期的相速度分布图. 从图中可以看出,红河断裂和小江断裂带(昆明附近)的低速异常仍然存在,说明它们可能是切穿地壳的超壳断裂. Lei等(2009)根据红河断裂下方15—70 km深度的低速异常认为,该断裂穿透地壳进入了上地幔. 胥颐等(2003)采用体波层析成像技术研究红河断裂认为,该断裂在莫霍附近的低速异常不会是地壳厚度增加所致,而是与壳-幔边界的热动力状况有关. 小江断裂带中部地区(昆明附近)热流值较高(周龙泉等,2009),是仅次于云南腾冲的高热流区,所以该地区出现低速异常很可能与深部的热作用有关. 滇东南区域高速异常的特性已不明显. 腾冲东南(龙陵附近)显示为高速异常,可能与早期火山作用中难以挥发的高密度物质残留体有关. 该结果与Wei等(2010)“在下地壳深度上,腾冲火山地区存在高速异常”观点基本吻合,但与Huang等(2002)的成像结果不一致. 在40 s周期时,台站射线路径减少且腾冲位于研究区边缘,此高速异常也可能是相速度层析成像分辨率变差所造成的假象.
云南的大多数地震震源深度为5—25 km,属于上、 中地壳范围. 从周期为8 s和18 s的相速度分布图来看,相速度低速异常区均属于云南地区的强地震活动区. 红河断裂北段、 小江断裂带附近、 腾冲地区和思茅地区的低速异常区与周龙泉等(2009)得出的云南地区低Q值区域基本一致. 因为低速异常体强度相对较弱,在横向挤压的构造应力场作用下,容易发生地震. 滇东南地区构造比较稳定,相速度分布图上呈现为高速异常,且该地区Q值偏高(周龙泉等,2009),是云南地区的弱地震活动带. 周期8 s和18 s相速度分布图中的低速区与低Q值区和高热流值区域(周龙泉等,2009)具有很好的一致性. 值得注意的是,滇东北地区虽为地震活动区,但在周期8 s的相速度分布图中表现为高速异常. 王椿镛等(2002)在研究川滇地区的地震分布时指出,地震发生与纵向上地壳深部孕震环境有关,而且川滇地区构造复杂,每次强烈地震的发生还应有其特定的背景,所以反演出S波速度结构后可进一步了解该地区的发震机制. 红河断裂作为一条重要的构造边界带,强震并没有沿整条断裂均匀分布,而仅多发生于北段(洱源—弥渡)上. 在短周期相速度分布图上,红河断裂南北段呈现差异,表明该断裂南北段周边介质存在物性差异,这可能是地震活动呈现南北差异的主要原因.
小江断裂带南端与红河断裂交汇部位是川滇活动块体的南部边界,从8—40 s相速度异常分布图中可以看出,该地区逐渐转变为高速异常. GPS观测资料显示,小江断裂带中段及北段的左旋滑动速率为8—10 mm/a,南段左旋走滑速率降低为4 mm/a(Shen et al,2005). Lev等(2006)利用剪切波分裂方法对川滇地区进行各向异性研究时得出,该区域高速异常体附近的快波方向转变为近南北向. 由此可见,高速异常体对川滇块体深部介质向南运动起到了阻挡的作用. 该结果与吴建平等(2013)P波层析成像结果相一致,进一步证实了川滇活动块体南部边界的形成不仅与红河断裂有关,而且可能与该地区中下地壳及上地幔存在高速异常体有关.
Royden等(1997)认为,青藏高原东缘中下地壳存在一个15 km厚的黏性流动层,向高原外流动. 由于鲜水河断裂以南的川滇地区地壳强度相对软弱,中下地壳物质流展的通道开阔,地壳在上千千米的尺度上逐渐增厚,并且导致地表大范围的隆升,形成了宽阔、 地貌反差不明显、 边界模糊的青藏高原东南缘. 通道流模型要求中下地壳存在能够发生流动的S波低速异常带,结合该地区存在高热流、 高泊松比、 高导层等现象以及本文中长周期相速度图所显示的低速区等,可能反映了该地区中下地壳物质的这种软弱特性. 若要详细了解该地区地壳物质的这种特性,则需要进一步反演出三维剪切波速度结构.
4. 结论
本文使用噪声层析成像方法获得了云南地区不同周期瑞雷面波相速度分布图像. 通过背景噪声数据获取的频散数据主要集中在短周期部分,因此主要揭示了云南地区的地壳结构特征. 研究结果显示,云南地区地壳内部存在明显的横向不均匀性,短周期相速度低速异常与地表断裂带分布和沉积层厚度密切相关. 在短周期相速度分布图上,红河断裂南北段呈现差异,表明红河断裂南北段周边介质存在物性差异,这可能是造成地震活动南北差异的主要原因. 在长周期相速度分布图上,红河断裂和小江断裂带的低速异常可能与深部热作用有关,它们可能是切穿地壳的超壳断裂. 相速度低速异常区(周期8—25 s)均属于云南地区的强地震活动区. 本文结果为下一步反演该地区的三维剪切波速度结构奠定了基础. 本研究中缺乏四川地区的资料,如果加入可取得更好的效果. 噪声层析成像方法采用实时记录的噪声资料且不依赖天然地震,这为获取精细地壳结构以及实时监测地下介质的物性变化提供了可能.
审稿专家为本文提出了建设性意见,澳大利亚麦考瑞大学杨英杰博士和俄罗斯圣彼得堡大学Yanovskaya教授提供了相关程序,云南大学胡家富教授指导分析解释过程,在此一并表示感谢.
-
图 1 二维交错网格速度-应力空间分布图(引自Levander,1988)
Figure 1. Spatial stencils of two-dimensional staggered grid for the velocity and stress update Dot and open circle represent horizontal and verical components of velocity,respectively; solid square and hollow square separately represent normal stress and tangential stress
图 3 地震记录 (a)直接法;(b)应力镜像-2阶格式法;(c)应力镜像-降阶法;(d)应力镜像-2阶速度展开法;(e)应力镜像-虚拟层赋零法;(f)应力镜像-紧致差分法
Figure 3. Seismic records (a)Direct method;(b)Stress image and second-order scheme method;(c)Stress image and reduced-order method;(d)Stress image and second-order velocity evolving method;(e)Stress image and fictitious layer assigns to zero method;(f)Stress image and compact difference method
图 4 传统PML(a)、M-PML(b)、C-PML(c)和MC-PML(d)条件下采用应力镜像-紧致差分法得到的t=139 ms时刻的波场快照(泊松比为0.48). 图中画圈处为不稳定机理
Figure 4. Wavefield snapshots computed by stress image and compact difference method when the Poisson’s ratio is 0.48 with PML(a),M-PML(b),C-PML(c) and MC-PML(d)at time instant t=139 ms. The red circle represents the instability
图 5 传统PML边界条件下泊松比为0.48时t=115,132,139和149 ms时刻的波场快照 (a)应力镜像-紧致差分法;(b)应力镜像-降阶法;(c)应力镜像-2阶速度展开法
Figure 5. Wavefield snapshots when the Poisson’s ratio is 0.48 with the classical PML at time instant t=115,132,139 and 149 ms (a)Stress image and compact difference method;(b)Stress image and reduced-order method;(c)Stress image and second-order velocity evolving method
图 6 单道波形(虚线)与解析解(实线)的对比 (a)直接法;(b)应力镜像-2阶格式法;(c)应力镜像-降阶法;(d)应力镜像-虚拟层赋零法;(e)应力镜像-2阶速度展开法;(f)应力镜像-紧致差分法
Figure 6. Comparison of waveforms between numerical(dashed line) and analytical solutions(solid line) (a)Direct method;(b)Stress image and second-order scheme method;(c)Stress image and reduced-order method;(d)Stress image and fictitious layer assigns to zero method;(e)Stress image and second-order velocity evolving method;(f)Stress image and compact difference method
表 1 传统PML在不同水平自由表面边界条件下的对比总结
Table 1 Comparison among implementations of planar free-surface boundary conditions with the classical PML
-
田坤, 黄建平, 李振春, 李娜, 孔雪, 李庆洋. 2013. 实现复频移完全匹配层吸收边界条件的递推卷积方法[J]. 吉林大学学报: 地球科学版, 43(3): 1022-1032. Tian K, Huang J P, Li Z C, Li N, Kong X, Li Q Y. 2013. Recursive convolution method for implementing complex frequency-shifted PML absorbing boundary condition[J]. Journal of Jilin University: Earth Science Edition, 43(3): 1022-1032 (in Chinese).
王守东. 2003. 声波方程完全匹配层吸收边界[J]. 石油地球物理勘探, 38(1): 31-34. Wang S D. 2003. Perfectly matched layer absorbing boundary of acoustic equation[J]. Oil Geophysical Prospecting, 38(1): 31-34 (in Chinese).
张伟. 2006. 含起伏地形的三维非均匀介质中地震波传播的有限差分算法及其在强地面震动模拟中的应用[D]. 北京: 北京大学地球与空间科学学院: 29-37. Zhang W. 2006. Finite Difference Seismic Wave Modeling in 3D Heterogeneous Media With Surface Topography and Its Implementation in Strong Ground Motion Study[D]. Beijing: School of Earth and Space science, Peking University: 29-37 (in Chinese).
Bayliss A, Jordan K E, Lemesurier B J, Turkel E. 1986. A fourth-order accurate finite-difference scheme for the computation of elastic-waves[J]. Bull Seism Soc Am, 76(4): 1115-1132.
Berenger J P. 1994. A perfectly matched layer for the absorption of electromagnetics waves[J]. Journal of Computational Physics, 114(2): 185-200.
Collino F, Tsohka C. 2001. Application of the perfectly matched layer model to the linear elastodynamic problem in anisotropic heterogeneous media[J]. Geophysics, 66(1): 294-307.
De Hoop A T. 1960. A modification of Cagniard's method for solving seismic pulse problems[J]. Appl Sci Res, 8(B): 349-356.
Gottschammer E, Olsen K B. 2001. Accuracy of the explicit planar free-surface boundary condition implemented in a fourth-order staggered-grid velocity-stress finite-difference scheme[J]. Bull Seism Soc Am, 91(3): 617-623.
Graves R W. 1996. Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences[J]. Bull Seism Soc Am, 86(4): 1091-1106.
Jastram C, Tessmer E.1994. Elastic modeling on a grid with vertically varying spacing[J]. Geophys Prosp, 42(4): 357-370.
Komatitsch D, Martin R. 2007. An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation[J]. Geophysics, 72(5): SM155-SM167.
Kristek J, Moczo P, Archuleta R J. 2002. Efficient methods to simulate planar free surface in the 3D 4(th)-order staggered-grid finite-difference schemes[J]. Studia Geophysica et Geodaetica, 46(2): 355-381.
Lan H Q, Zhang Z J. 2010. Comparative study of the free-surface boundary condition in two-dimensional finite-difference elastic wave field simulation[J]. J Geophys Eng, 8(2): 275-286.
Lele S K. 1992. Compact finite difference schemes with spectral-like resolution[J]. J Comput Phys, 103(1): 16-42.
Levander A R. 1988. Fourth-order finite-difference P-SV seismograms[J]. Geophysics, 53(11): 1425-1436.
Ohminato T, Chouet B A. 1997. A free-surface boundary condition for including 3D topography in the finite-difference method[J]. Bull Seism Soc Am, 87(2): 494-515.
Oprsal I, Zahradnik J. 1999. Elastic finite-difference method for irregular grids[J]. Geophysics, 64(1): 240-250.
Pitarka A, Irikura K. 1996. Modeling 3D surface topography by finite-difference method: Kobe-JMA station site, Japan, case study[J]. Geophys Res Lett, 23(20): 2729-2732.
Robertsson J. 1996. A numerical free-surface condition for elastic/viscoelastic finite-difference modeling in the presence of topography[J]. Geophysics, 61(6): 1921-1934.
Zahradnik J. 1995. Simple elastic finite-difference scheme[J]. Bull Seism Soc Am, 85(6): 1879-1887.
Zeng C, Xia J H, Miller R D, Tsoffias G P. 2011. Application of the multiaxial perfectly matched layer (M-PML) to near-surface seismic modeling with Rayleigh waves[J]. Geophysics, 76(3): T43-T52.