地球物理学报  2017, Vol. 60 Issue (12): 4671-4680   PDF    
球坐标系下谱元法三维地震波场模拟
董兴朋, 杨顶辉     
清华大学数学科学系, 北京 100084
摘要:谱元法已成为区域性乃至大陆性尺度地震波场模拟的重要工具.对于区域或大陆尺度层析成像而言,地球曲率不可忽略,此时模拟地震波传播采用球坐标系更为合适.本文从球坐标系下弹性波动方程弱形式出发,基于球坐标系变分原理给出了球坐标系下求解三维地震波方程的谱元法.另一方面,计算Fréchet敏感核是进行全波形反演的关键,本文借助伴随原理,推导了全波走时层析成像三维Fréchet敏感核表达式.为了验证球坐标系下谱元法的精度,我们将数值模拟结果与normal mode方法得到的解析解在1-D PREM模型下进行了对比.同时,我们将此方法应用到华北克拉通区域,以期获得地球内部结构精确成像.基于3-D全球径向各向异性地幔模型S362ANI和3-D地壳模型Crust1.0,我们建立了华北克拉通初始3-D背景模型,并将数值模拟结果与实际观测台站记录波形资料进行对比分析,利用互相关方法提取走时残差,最后给出了Fréchet敏感核在3-D空间中的分布,这些工作为下一步开展球坐标系下三维大尺度全波形反演奠定了基础.
关键词: 球坐标系      谱元法      Fréchet敏感核      地震波模拟     
Numerical modeling of the 3-D seismic wavefield with the spectral element method in spherical coordinates
DONG Xing-Peng, YANG Ding-Hui     
Department of Mathematical Sciences, Tsinghua University, Beijing 100084, China
Abstract: The spectral element method has become an important tool for seismic wavefield simulation of regional and even continental scales. For regional or continental tomography, the curvature of the earth cannot be neglected, so it is more appropriate to use the spherical coordinate system to simulate the propagation of seismic waves. Starting from the weak form of elastic wave equation in spherical coordinates, the basic theory of the Legendre spectral element method in a spherical coordinate system is expounded based on variational principle of the spherical coordinate system. On the other hand, calculation of Fréchet kernel is critical to the full waveform inversion. With the help of the adjoint method, we derive the expression of 3-D Fréchet kernel.The numerical results obtained by the forward modeling, which are compared with the analytical solutions obtained by the normal mode method to verify the accuracy of the 1-D PREM model. Simultaneously, we apply this method to the North China Craton to create accurate imaging of the earth's interior structure. Based on the 3-D global radially anisotropic mantle model S362ANI with the 3-D crustal model Crust1.0, the numerical simulation of seismic wave propagation in North China is carried out, and compared with the data recorded by observation stations. We obtain traveltime residuals from cross correlation and the spatial distribution of full wave traveltime 3-D Fréchet kernel, which lays a foundation for the next large full-waveform inversion in the spherical coordinate system.
Key words: Spherical coordinates    Spectral element method    Fréchet kernel    Seismic modeling    
1 引言

数值求解地震波方程是实施全波形反演的先决条件,在区域乃至全球地震学中,精确计算实际三维地球模型的理论地震图已经成为一种现实需要(Komatitsch and Tromp, 1999, 2002a).目前,地震波数值模拟方法主要有:有限差分法(Alford et al., 1974; Kelly et al., 1976)、伪谱法(Fornberg, 1988)、有限元法(Lysmer and Drake, 1972; Liu et al., 2014a)和谱元法(Priolo and Seriani, 1991; Komatitsch and Vilotte, 1998; Komatitsch and Tromp, 2002b)等.有限差分法以差分近似微分,具有编程简单,内存消耗小等优点,但也存在难以处理起伏地表和自由界面,并存在明显数值频散等问题(宋国杰等,2010Liu et al., 2014b刘有山,2015).针对有限差分数值频散问题,Yang等发展了近似解析离散化(Nearly Analytic Discrete, NAD)方法(Yang et al., 2003, 2004, 2006),该方法的思想是利用波场值及其梯度值构造局部插值,与传统的有限差分法相比,NAD方法考虑了更多的波场信息,有利于高精度地近似空间高阶导数,但其一个缺点在于网格剖分仍不灵活,对起伏地形的描述只能采用阶梯状网格(宋国杰,2011).伪谱法使用快速Fourier变换计算空间导数值,是一种全局方法,具有数值频散小和收敛速度快的特点.Gazdag(1981)Kosloff和Baysal(1982)分别将伪谱法用于声波和弹性波的模拟(Gazdag, 1981; Kosloff et al., 1984),随后,Reshef等(1988)使用伪谱法模拟了地震波在三维介质中的传播,由于伪谱法为全局长算子,因而只能处理简单模型,并存在边界不易处理的问题.有限元法是基于积分形式(弱形式)的数值模拟方法,基本理论是利用变分原理得到地震波方程的弱形式.有限元法克服了上述几种方法在处理边界问题上的缺点,能够适应于任意地质体的模型并且自动满足自由边界条件(刘有山等,2013; 刘少林等,2014; Liu et al., 2014a),但低阶有限元法精度较低,频散严重(Marfurt, 1984),而高阶有限元法则会因为插值导致Runge现象而出现伪波问题(Komatitsch and Vilotte, 1998).

近20年来发展起来的谱元法综合了有限元法的网格剖分灵活性和并具有谱方法的高精度指数收敛等优点(汪文帅等,2012刘有山等,2014),是进行大尺度地震波模拟一种较为理想的方法(Komatitsch and Tromp, 2002a, 2002b, 2003).它的基本思想是在单元内做谱展开,常选取Chebyshve或Legendre正交多项式作为插值基函数,通过提高插值阶数来提高解的精度.1984年,Petera最早将Chebyshev型谱元法引入计算流体之中,随后,Priolo和Seriani(1991)使用该方法模拟声波在二维介质中的传播,并发现当插值阶数为8阶即能精确模拟声波传播,但由于Chebyshev正交多项式为带权正交,无法获得对角的质量矩阵.Maday和Patera(1989)将以Legendre多项式为基的Lagrange插值引入到谱元中,并采用Gauss-Lobatto-Legendre(GLL)积分,从而产生对角质量矩阵,降低了存储量,提高了计算效率,目前大尺度地震波模拟与成像中通常采用Legendre型谱元法.

对于区域性乃至全球性尺度而言,所考虑的模型为球体或球体的一部分,因而需要在球坐标进行地震数据的表达与处理,并能够避免传统的直角坐标投影所引入的误差.另外,球坐标系谱元法还有如下优点:能自然拟合分层地层曲率,大尺度地震波模拟中不需要做坐标旋转变换,而只需要做尺度伸缩变换;网格剖分容易,沿经纬度均匀剖分即可.本文基于球坐标系变分原理,详细推导了球坐标系下Legendre型谱元法基本理论,并在一维PREM地球模型中,与normal mode得到的参考解对比验证了球坐标系下SEM方法的精度.之后,以华北克拉通地区为例,基于3-D背景模型,将正演模拟得到的三分量波形数据与实际观测台站记录的波形资料对比,并利用互相关方法求得走时残差,最后给出了全波走时层析成像三维Fréchet敏感核表达式和相关数值算例.

2 球坐标系下的谱元法 2.1 三维球坐标系下弹性波方程及其弱形式

地震波在地球介质中的传播常用弹性波方程描述,方程的求解可以基于强形式或者弱形式(Komatitsch and Tromp, 1999; Liu et al., 2017a).在强形式下,直接将弹性波方程和相应的边界条件写成微分形式,例如有限差分法和伪谱法;弱形式即采用积分法求解弹性波方程,如有限元法和谱元法等.

球坐标系(r, θ, φ)下,弹性波方程基本形式如下(Lapwood and Usami, 1981):

(1)

初始条件为:

(2)

其中,u为位移场,v表示速度场,ρ为质量密度,f为震源项,c为四阶弹性张量,共有21个独立分量,方程(1)中的▽为球坐标系下梯度算子,具有如下形式:

对方程(1)两端乘以任意一个测试向量函数ψ并在求解区域进行积分,然后利用分部积分可以得到波动方程的弱形式:

(3)

其中Γ=Γf+Γa为计算域Ω的边界,ΓaΓf分别表示人工边界和自由边界.在自由边界Γf上,0,在人工边界条件上,需要特殊处理以压制虚假反射波.

在大尺度地震模拟中,因为研究区域比较大,震源的破裂往往是几公里或几十公里,相对于如此大的模型而言可以视为点源(Dahlen and Tromp, 1998),更复杂的震源的处理方式是我们以后研究的重点.地震震源可以由方程(1)中点力源f来表示,其形式如下(Komatitsch and Tromp, 1999):

(4)

其中,xs为震源位置,S(t)为震源时间函数,Mw为地震矩张量,δ(xxs)为震源位置xs处的狄拉克分布.

2.2 空间离散

三维谱元法空间离散需要将计算区域分解为N个互不重叠的六面体单元,利用自然边界条件和自由边界条件,方程(3)可以写成如下形式:

(5)

将每一个六面体单元都映射到标准立方体[-1, 1]3,我们选择的插值节点同Gauss-Legendre-Lobatto(GLL)积分节点一致(Cohen, 2002; Liu et al., 2017a).每个单元的插值基函数如下:

(6)

其中,αi表示i阶插值基函数,ξγη为等参坐标系下三方向坐标,l为拉格朗日插值多项式.基于方程(6)和插值点上的离散值ue, i, ve, i, σe, iψe, i,单元上位移、速度、应力与测试函数的连续值可以近似为

(7)

这里n为插值多项式的阶数,将方程(7)代入方程(5)中并使用GLL积分法则,可以得到时间域一阶线性常微分方程:

(8)

式(8)中U为全局位移向量,M为对角化的质量矩阵,K1K2为刚度矩阵,F为震源项.

2.3 时间离散

弹性波方程经空间离散后得到的常微分方程组,还需要进行时间离散化才能得到时间上的递推格式.常用的时间域算法主要有Newmark算法(Kane et al., 2000)、Runge-Katta方法、辛算法(刘少林等, 2013, 2015; Liu et al., 2014a, 2015)等,由于大尺度地震波模拟的时间长度可能超过150 s(Tape et al., 2009),需要高效、长时程能力的时间域算法,Liu等(2017b)在这方面做了大量出色的工作,提出了一种修正辛分块龙格-库塔方法,该方法在计算精度、稳定性和长时程追踪能力方面均有上佳的表现.但是高阶方法往往需要多次矩阵与向量相乘(Liu et al., 2017a, 2017c),计算量较大,在现有的计算资源条件下较难应用于大尺度地震波传播与成像(Tape et al., 2010).本文使用交错时间离散格式实现时间域递推,该方法虽然只有二阶精度,但需要中间变量少,易于处理边界条件.其递推形式如下:

(9)

在数值模拟地震波传播时,在有限的计算资源条件下模拟无限或半无限区域中地震波传播规律,这就需要引入人工边界对研究区域进行截断.如果处理不当,在人工边界上产生的虚假反射波会污染整个计算区域波场.针对如何处理人工边界,前人做了大量的工作,主要分为边界条件和吸收边界层两种策略.边界条件通常基于波动方程的旁轴近似(Clayton and Engquist, 1977; Komatitsch and Vilotte, 1998),但对于大入射角或面波不能有效地吸收,并且存在数值稳定性的问题;一种相对高效的边界处理技术是Bérenger于1994年提出的最佳匹配层(perfectly matched layer, 简称PML)方法,该方法可以较好的吸收来自各个方向和各个频率的入射波,且不会产生明显的边界虚假反射波(Bérenger, 1994).

经典的PML方法会导致一个新的且更大的边界方程系统,而使用各向异性完全匹配层(Anisotropic perfectly matched layer, 简称APML)方法则可以避免这种情况(Zheng and Huang, 2002).APML方法是在PML区域引入单轴各向异性有耗介质,通过适当选择单轴各向异性介质的本构参数使入射地震波无反射地进入PML并且被衰减吸收.APML的优点是:1)地震波动方程保持原有的物理形式,易于理解;2)在三维情况下,所占用的内存最小.因而本文选择使用APML方法处理人工边界问题.

3 数值模拟与分析 3.1 计算结果与理论解拟合对比

本节中,正演模拟基于一维的PREM地球模型,模拟最短周期为5 s, 理论地震图计算时间为500 s, 时间步长0.05 s,研究区范围为110°E—126°E,34°N—44°N,深度方向为地表至地下271 km.沿纬度方向,经度方向和深度方向分别划分126,145和32个谱元单元,即共1150720个六面体单元,单元采用4阶Lagrange多项式插值.研究所使用的地震事件为201304220517560,发震位置为122.32°E,42.9°N,震源深度12 km,震源时间函数为带通滤波的Heaviside函数.震源矩张量参数由CMT(http://ds.iris.edu/spud/momenttensor)得到(单位:N·m),如下:

地震事件201304220517560垂直方向速度分量在地球表面传播的波场快照如图 1所示,图中可以清晰看到地震波场在地表传播形态,在t=50 s时,地震波场刚刚抵达研究区上边界,APML边界能有效吸收入射波;在t=100 s时,入射波几乎被完全吸收.

图 1 地震波垂向速度分量在地球表面的50 s (a)、100 s (b)波场快照 Fig. 1 50 seconds (a)、100 seconds (b) snapshots of the vertical velocity component on the surface of the earth

图 2显示了台站LN.YKO(40.68°N, 122.60°E)处球坐标系下谱元法与normal mode方法垂直分量的波形对比图,二者吻合较好,但在波峰(谷)处幅值存在小的偏差,出现这种现象的一个主要原因是:我们在使用normal mode方法求解理论地震图时,受限于计算时间,只选择了径向阶数n≤18的部分,损失了更高阶的成分.Normal mode位移空间分布可以写成如下形式(Morse and Feshbach, 1953; Dahlen and Tromp, 1998):

图 2 球坐标系下SEM方法(蓝)和normal mode法(红)在台站LN.YKO处垂直分量速度地震图对比 Fig. 2 Comparision of vertical velocity waveforms generated by spherical coordinate system SEM (blue) and normal mode method (red) at LN.YKO

其中m表示方位阶数,n表示径向阶数.对于频率低于120 mHz而言,n的取值范围为:0≤n≤341,一共有57710种环形振荡和94887种球形振荡需要计算,计算量巨大,因而本文只计算了径向阶数小于等于18的部分.

3.2 实际数据与理论波形对比

本节基于3-D全球径向各向异性模型S362ANI(Kustowski et al., 2008)和3-D地壳模型crust1.0(Laske et al., 2013)模型建立了华北克拉通三维初始速度模型,并在此基础上进行大尺度地震波场正演模拟.地震事件和所选用台站位置如图 3所示,球坐标系下SEM方法三分量位移正演模拟结果与实际观测值对比见图 4,可以看到,P波成分相位和振幅均吻合较好,而S波成分相位和振幅差别较大.这也从侧面说明了全波走时层析成像的必要性,即:它能够利用整个波场的信息(P波,S波和面波)反演模型参数.

图 3 地震事件201304220517560(红色五角星)与台站SD.RCH(蓝色圆圈)位置 Fig. 3 The location of earthquake event 201304220517560 (red pentagram) and SD.RCH station (blue cycle)
图 4 球坐标系下SEM在台站SD.RCH处三分量位移(蓝色)与实际观测值(红色)对比 Fig. 4 Comparison of waveforms generated by spherical coordinate system SEM (blue) and observed data (red) at SD.RCH
3.3 互相关走时Fréchet敏感核

计算Fréchet敏感核是进行全波形反演的关键,它表征了空间中各点对于模型扰动的敏感性(Huang et al., 2016).在全波形走时层析成像中,需要最小化如下的误差函数(Zhao et al., 2000; Dahlen et al., 2000; Hung et al., 2000; Tromp et al., 2005):

其中,Δτ为走时残差,m为模型参量.基于连续伴随方法(Fichtner, 2011),伴随方程可以写成如下形式:

上式中vσ分别为伴随速度场和伴随应力场,xr为地表台站位置,误差敏感核▽mχ满足如下方程:

走时残差Δτ一般通过理论波形与实际数据做互相关方法得到(详见附录B),我们给出了3.2节中台站SD.RCH处三分量位移的互相关结果(见图 5,其中CC,dT与dA分别表示互相关系数,走时延迟和振幅扰动).从图 5可以看到,一共有7个互相关走时窗口,得到7个互相关走时残差,其中垂直分量三个(1.6 s, 1.7s和9.2 s),南北分量两个(1.5 s, 3.3 s),东西分量两个(5.2 s, 3.4 s).图 6为走时Fréchet敏感核在地表 0 km、地下30 km和地下50 km分布图,可以看到:随着深度增加,敏感核强度逐渐减弱;敏感核水平切面轮廓为近似椭圆,之所以不是标准椭圆主要是因为初始模型是三维非均匀模型.

图 5 互相关方法提取台站SD.RCH处走时残差 Fig. 5 Traveltime residuals from cross correlation at SD.RCH
图 6 走时Fréchet敏感核在地表 0 km (a)、地下30 km (b)和50 km (c)水平切片图,五角星表示震源,三角形表示台站(色标单位s2·m-4) Fig. 6 Horizontal slice of traveltime Fréchet kernel at 0 km (a), 30 km (b) and 50 km (c) earthquake (red pentagram), station (triangle) (colorbar unit s2·m-4)
4 讨论与结论

谱元法是近年来在地震波数值模拟领域流行的一种数值方法,高效且高精度的特性使其成为大尺度地震波场模拟一种较为理想的工具.对于区域乃至全球层析成像而言,地球曲率不可忽略,因而在球坐标系下进行三维大尺度地震波传播数值模拟尤为重要.相较于传统的直角坐标系SEM方法,球坐标系SEM方法有着如下的优势:1)能自然拟合分层地层曲率;2)大尺度地震波模拟中不需要做旋转变换,只需要做尺度伸缩变换;3)容易剖分,沿经纬度均匀剖分即可;4)能够避免传统的直角坐标投影所引入的误差.

对于弹性波方程,利用变分原理,得到了弹性波动方程的弱形式,并详细阐述了球坐标系下Legendre型谱元法基本理论与数学公式.为了验证本文所给方法的精度,我们基于一维PREM地球模型,将球坐标系SEM方法计算结果同normal mode解析解相对比,两种方法得到的结果吻合较好,只是在波峰(谷)处存在微小差别,主要原因是我们在使用normal mode计算理论地震图时,限于计算能力,只计算了径向阶数n≤18部分.随后,以华北克拉通地区为例,基于3-D背景模型,将正演模拟得到的三分量波形数据与实际观测台站记录的波形资料对比,并利用互相关方法求的走时残差,最后给出了全波走时层析成像三维Fréchet敏感核表达式和相关数值算例.这些结果为下一步开展基于球坐标SEM方法进行大尺度全波形反演奠定了方法基础.

附录A  球坐标系下应力张量散度和应变张量

在球坐标系中,应力张量散度▽·σ和应变张量ε精确形式如下(Lapwood and Usami, 1981):

(A1)

(A2)

(A3)

(A4)

(A5)

(A6)

(A7)

(A8)

(A9)

附录B  互相关走时残差

由于L2范数下地震波形误差带来了极强的非线性,Luo和Schuster(1991)提出利用理论地震图与实际数据做互相关的方法来提取走时延迟.考虑接收器xr处观测波形的第i个分量,ui(xr, t; m)obs,和相应的理论波形ui(xr, t; m)calm为模型变量,互相关走时函数定义如下:

(B1)

我们寻找最佳的时间移位τ使得理论波形“最好”匹配观测波形,走时残差Δτ定义为互相关函数取得最大值时的时间移位:

(B2)

上式中T为预估的观测波形与理论波形走时差异的最大值,在Δτ处互相关函数满足必要条件:

(B3)

致谢

本文得到国家自然科学基金重点项目“基于全波方程的岩石圈结构成像方法及其应用研究”(资助号:41230210)的资助.

参考文献
Alford R M, Kelly K R, Boore D M. 1974. Accuracy of finite difference modeling of the acoustic wave equation. Geophysics, 39(6): 834-842. DOI:10.1190/1.1440470
Bérenger J P. 1994. A perfectly matched layer for the absorption of electromagnetic waves. Journal of Computational Physics, 114(2): 185-200. DOI:10.1006/jcph.1994.1159
Clayton R, Engquist B. 1977. Absorbing boundary conditions for acoustic and elastic wave equations. Bull. Seismol. Soc. Am., 67(6): 1529-1540.
Cohen G. 2002. Higher-Order Numerical Methods for Transient Wave Equations. Berlin Heidelberg: Springer-Verlag.
Dahlen F A, Tromp J. 1998. Theoretical Global Seismology. Princeton, NJ: Princeton University Press.
Dahlen F A, Hung S H, Nolet G. 2000. Fréchet kernels for finite-frequency traveltimes-Ⅰ. Theory. Geophys. J. Int., 141(1): 157-174. DOI:10.1046/j.1365-246X.2000.00070.x
Fichtner A, Bunge H P, Igel H. 2006a. The adjoint method in seismology:Ⅰ. Theory. Physics of the Earth and Planetary Interiors, 157(1-2): 86-104. DOI:10.1016/j.pepi.2006.03.016
Fichtner A, Bunge H P, Igel H. 2006b. The adjoint method in seismology:Ⅱ. Applications:traveltimes and sensitivity functionals. Physics of the Earth and Planetary Interiors, 157(1-2): 105-123. DOI:10.1016/j.pepi.2006.03.018
Fichtner A. 2011. Full Seismic Waveform Modelling and Inversion. Berlin Heidelberg: Springer-Verlag.
Fornberg B. 1988. The pseudospectral method:accurate representation of interfaces in elastic wave calculations. Geophysics, 53(5): 625-637. DOI:10.1190/1.1442497
Gazdag J. 1981. Modeling of the acoustic wave equation with transform methods. Geophysics, 46(6): 854-859. DOI:10.1190/1.1441223
Huang X Y, Yang D H, Tong P, et al. 2016. 3D nearly analytic central difference method for computation of sensitivity kernels of wave-equation-based seismic tomography. Bull. Seismol. Soc. Am., 106(6): 2877-2899. DOI:10.1785/0120150309
Hung S H, Dahlen F A, Nolet G. 2000. Fréchet kernels for finite-frequency traveltime-Ⅱ. Examples. Geophys. J. Int., 141(1): 175-203. DOI:10.1046/j.1365-246X.2000.00072.x
Kane C, Marsden J E, Ortiz M, et al. 2000. Variational integrators and the Newmark algorithm for conservative and dissipative mechanical systems. International Journal for Numerical Methods in Engineering, 49(10): 1295-1325. DOI:10.1002/(ISSN)1097-0207
Kelly K R, Ward R W, Treitel S, et al. 1976. Synthetic seismograms:a finite-difference approach. Geophysics, 41(1): 2-27. DOI:10.1190/1.1440605
Komatitsch D, Vilotte J P. 1998. The spectral element method:an efficient tool to simulate the seismic response of 2D and 3D geological structures. Bull. Seismol. Soc. Am., 88(2): 368-392.
Komatitsch D, Tromp J. 1999. Introduction to the spectral element method for three-dimensional seismic wave propagation. Geophys. J. Int., 139(3): 806-822. DOI:10.1046/j.1365-246x.1999.00967.x
Komatitsch D, Tromp J. 2002a. Spectral-element simulations of global seismic wave propagation-Ⅰ. Validation. Geophys. J. Int., 149(2): 390-412. DOI:10.1046/j.1365-246X.2002.01653.x
Komatitsch D, Tromp J. 2002b. Spectral-element simulations of global seismic wave propagation-Ⅱ. Three-dimensional models, oceans, rotation, and self-gravitation. Geophys. J. Int., 150(1): 303-318.
Komatitsch D, Tromp J. 2003. A perfectly matched layer absorbing boundary condition for the second-order seismic wave equation. Geophys. J. Int., 154(1): 146-153. DOI:10.1046/j.1365-246X.2003.01950.x
Kosloff D D, Baysal E. 1982. Forward modeling by a Fourier method. Geophysics, 47(10): 1402-1412. DOI:10.1190/1.1441288
Kustowski B, Ekström, Dziewoński A M. 2008. Anisotropic shear-wave velocity structure of the Earth's mantle:A global model. J. Geophys. Res., 113(B6): B06306.
Lapwood E R, Usami T. 1981. Free Oscillations of the Earth. Cambridge: Cambridge University Press.
Laske G, Masters G, Ma Z T, et al. 2013. Update on CRUST1.0-A 1-degree Global Model of Earth's Crust.//EGU General Assembly 2013. Vienna, Austria.
Liu S L, Li X F, Liu Y S, et al. 2013. Symplectic RKN schemes for seismic scalar wave simulations. Chinese Journal of Geophysics, 56(12): 4197-4205. DOI:10.6038/cjg20131222
Liu S L, Li S F, Liu Y S, et al. 2014. Dispersion analysis of triangle-based finite element method for acoustic and elastic wave simulations. Chinese Journal of Geophysics, 57(8): 2620-2630. DOI:10.6038/cjg20140821
Liu S L, Li X F, Wang W S, et al. 2014a. A mixed-grid finite element method with PML absorbing boundary conditions for seismic wave modelling. J. Geophys. Eng., 11(5): 055009. DOI:10.1088/1742-2132/11/5/055009
Liu S L, Li X F, Wang W S, et al. 2014b. A new kind of optimal second-order symplectic scheme for seismic wave simulations. Sci. China Earth Sci., 57(4): 751-758. DOI:10.1007/s11430-013-4805-0
Liu S L, Li X F, Wang W S, et al. 2015. A modified symplectic scheme for seismic wave modeling. J. Appl. Geophys., 116: 110-120. DOI:10.1016/j.jappgeo.2015.03.007
Liu S L, Li X F, Wang W S, et al. 2015. A symplectic RKN scheme for solving elastic wave equations. Chinese Journal of Geophysics, 58(4): 1355-1366. DOI:10.6038/cjg20150422
Liu S L, Yang D H, Dong X P, et al. 2017a. Element-by-element parallel spectral-element methods for 3-D teleseismic wave modeling. Solid Earth, 8: 969-986. DOI:10.5194/se-8-969-2017
Liu S L, Yang D H, Ma J. 2017b. A modified symplectic PRK scheme for seismic wave modeling. Computers & Geosciences, 99: 28-36.
Liu S L, Yang D H, Lang C, et al. 2017c. Modified symplectic schemes with nearly-analytic discrete operators for acoustic wave simulations. Computer Physics Communications, 213: 52-63. DOI:10.1016/j.cpc.2016.12.002
Liu Y S, Teng J W, Liu S L, et al. 2013. Explicit finite element method with triangle meshes stored by sparse format and its perfectly matched layers absorbing boundary condition. Chinese J. Geophys., 56(9): 3085-3099. DOI:10.6038/cjg20130921
Liu Y S, Teng J W, Xu T, et al. 2014. Numerical modeling of seismic wavefield with the SEM based on Triangles. Progress in Geophysics, 29(4): 1715-1726. DOI:10.6038/pg20140430
Liu Y S. 2015. Seismic wavefield modeling using the spectral element method and migration/inversion. Beijing: Institute of Geology and Geophysics, Chinese Academy of Sciences.
Luo Y, Schuster G T. 1991. Wave-equation traveltime inversion. Geophysics, 56(5): 645-653. DOI:10.1190/1.1443081
Maday Y, Patera A T. 1989. Spectral Element Methods for the Incompressible Navier-stokes Equations. York: ASME: 71-143.
Marfurt K J. 1984. Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations. Geophysics, 49(5): 533-549. DOI:10.1190/1.1441689
Morse P M, Feshbach H. 1953. Methods of Theoretical Physics. New York: McGraw-Hill.
Priolo E, Seriani G. 1991. A numerical investigation of chebyshev spectral element method for acoustic wave propagation.//Proceedings of the 13th IMACS Conference on Comparative Applied Mathematics. Dublin:Ireland, 551-556.
Reshef M, Kosloff D, Edwaeds M, et al. 1988. Three-dimensional elastic modeling by the Fourier method. Geophysics, 53(9): 1184-1193. DOI:10.1190/1.1442558
Song G J, Yang D H, Chen Y L, et al. 2010. Non-uniform grid algorithm based on the WNAD method and elastic wave-field simulations. Chinese J. Geophys., 53(8): 1985-1992. DOI:10.3969/j.issn.0001-5733.2010.08.025
Song G J. 2010. Parallel WNAD algorithm and its wave fields simulation[Ph. D] (in Chinese). Beijing:Tsinghua University.
Tape C, Liu Q Y, Maggi A, et al. 2009. Adjoint tomography of the southern California crust. Science, 325(5943): 988-992. DOI:10.1126/science.1175298
Tape C, Liu Q Y, Maggi A, et al. 2010. Seismic tomography of the southern California crust based on spectral-element and adjoint methods. Geophys. J. Int., 180(1): 433-462. DOI:10.1111/gji.2010.180.issue-1
Tromp J, Carl T, Liu Q Y. 2005. Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels. Geophys. J. Int., 160(1): 195-216.
Wang W S, Li X F, Lu M W, et al. 2012. Structure-preserving modeling for seismic wavefields based upon a multisymplectic spectral element method. Chinese J. Geophys., 55(10): 3427-3439. DOI:10.6038/j.issn.0001-5733.2012.10.026
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. Bull. Seismol. Soc. Am., 93(2): 882-890. DOI:10.1785/0120020125
Yang D H, Lu M, Wu R S, et al. 2004. An optimal nearly analytic discrete method for 2D acoustic and elastic wave equations. Bull. Seismol. Soc. Am., 94(5): 1982-1992. DOI:10.1785/012003155
Yang D H, Peng J M, Lu M, et al. 2006. Optimal nearly analytic discrete approximation to the scalar wave equation. Bull. Seismol. Soc. Am., 96(3): 1114-1130. DOI:10.1785/0120050080
Zhao L, Jordan T H, Chapman C H. 2000. Three-dimensional Fréchet differential kernels for seismicdelay times. Geophys. J. Int., 141(3): 558-576. DOI:10.1046/j.1365-246x.2000.00085.x
Zheng Y B, Huang X J. 2002. Anisotropic perfectly matched layers for elastic waves in Cartesian and curvilinear coordinates. Earth Resources Laboratory, Consortium Report. Cambridge:MIT.
刘少林, 李小凡, 刘有山, 等. 2013. 求解声波方程的辛RKN格式. 地球物理学报, 56(12): 4197–4205. DOI:10.6038/cjg20131222
刘少林, 李小凡, 刘有山, 等. 2014. 三角网格有限元法声波与弹性波模拟频散分析. 地球物理学报, 57(8): 2620–2630. DOI:10.6038/cjg20140821
刘少林, 李小凡, 汪文帅, 等. 2015. 求解弹性波方程的辛RKN格式. 地球物理学报, 58(4): 1355–1366. DOI:10.6038/cjg20150422
刘有山, 滕吉文, 刘少林, 等. 2013. 稀疏存储的显式有限元三角网格地震波数值模拟及其PML吸收边界条件. 地球物理学报, 56(9): 3085–3099. DOI:10.6038/cjg20130921
刘有山, 滕吉文, 徐涛, 等. 2014. 三角网格谱元法地震波场数值模拟. 地球物理学进展, 29(4): 1715–1726. DOI:10.6038/pg20140430
刘有山. 2015. 谱元法地震波场数值模拟与偏移成像. 北京: 中科院地质与地球物理研究所.
宋国杰, 杨顶辉, 陈亚丽, 等. 2010. 基于WNAD方法的非一致网格算法及其弹性波场模拟. 地球物理学报, 53(8): 1985–1992. DOI:10.3969/j.issn.0001-5733.2010.08.025
宋国杰. 2011. 并行WNAD算法及其波场模拟[博士论文]. 北京: 清华大学.
汪文帅, 李小凡, 鲁明文, 等. 2012. 基于多辛结构谱元法的保结构地震波场模拟. 地球物理学报, 55(10): 3427–3439. DOI:10.6038/j.issn.0001-5733.2012.10.026