地球物理学报  2015, Vol. 58 Issue (4): 1355-1366   PDF    
求解弹性波方程的辛RKN格式
刘少林1, 李小凡1, 汪文帅2, 朱童3    
1. 中国科学院地质与地球物理研究所, 中国科学院地球与行星物理重点实验室, 北京 100029;
2. 宁夏大学数学与计算机学院, 银川 750021;
3. 中国石油化工股份有限公司石油物探技术研究院, 南京 210014
摘要:将弹性波方程变换至Hamilton体系, 构造适用于弹性波模拟的高效显式二阶辛Runge-Kutta-Nyström (RKN)格式, 运用根数理论得到此格式的阶条件方程组. 通过给定系数的限定条件, 得到方程的对称解. 为了使时间离散误差达到极小,提出数值频率与真实频率比较,通过Taylor展开,得到关于辛系数的限定方程,求解方程组得到最小频散辛RKN格式. 对比分析时间演进方程的稳定性,得到使库朗数达到极大值的限定方程,求解方程组得到最稳定辛RKN格式. 发现此两种格式为同一格式. 新得到的辛RKN格式不依赖于空间离散方法,为了对比的需要,选取有限差分法进行空间离散. 在频散、稳定性分析中,与常见辛格式对比,从理论上分析了本文提出的格式在数值频散压制、稳定性提升等方面的优势, 数值实验进一步证实了理论分析的正确性.
关键词弹性波方程     辛RKN格式     稳定性条件     频散关系     有限差分法    
A symplectic RKN scheme for solving elastic wave equations
LIU Shao-Lin1, LI Xiao-Fan1, WANG Wen-Shuai2, ZHU Tong3    
1. Key Laboratory of Earth and Planetary Physics, Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China;
2. School of Mathematics and Computer Science, Ningxia University, Yinchuan 750021, China;
3. SINOPEC Geophysical Research Institute, Nanjing 210014, China
Abstract: The construction of low dispersive and strong stable numerical schemes is critical for the investigation of wave propagation. Solving seismic wave equations in the time domain, the dispersion and stability are related to the temporal and spatial discretizations. Here we focus on the temporal discretization and develop a high-efficient time integration scheme for the elastic wave equation.Following the transformation of the elastic seismic wave equation into a Hamiltonian system, an explicit second-order symplectic Runge-Kutta-Nyström (RKN) scheme is proposed for elastic wave modeling. The order conditions are obtained by the rooted trees theory. To obtain symmetric solutions, we assign a constraint condition to the symplectic coefficients. In order to minimize the error in temporal discretization, we propose a method that compares numerical angular frequency with real frequency, and obtain a constraint equation based on Taylor expansion. A symplectic RKN scheme with minimized numerical dispersion is developed by solving the order conditions associating with constraint equations. Through analyzing the time advance equation, we construct another constraint equation to allow stability limit to attain its maximal value, and then another symplectic RKN scheme with the same symplectic coefficients as the first scheme is developed. The optimization process for determining the coefficients of the RKN scheme is independent of spatial discretizations. For the ease of comparison, we simply choose finite difference method to approximate spatial derivatives. Generally, the numerical dispersion amount is lower than the conventional second-order schemes, and is slightly larger than the high-order schemes. The maximal Courant number for the new scheme is the largest compared to those for the conventional schemes. The largest time increment for the new scheme is nearly two times as large as those for the conventional second-order schemes, and is about 1.5 times as large as those for the high-order ones. In a practical simulation, all these schemes discussed in the paper consume approximately the same computer memory, but the numerical accuracy and computational speed are quite different. The numerical accuracy and complexity of our scheme indicate that it may be a good candidate for the balance of numerical accuracy and complexity.Both theoretical analyses and numerical experiments show that our scheme is superior to conventional symplectic schemes in some aspects including suppressing dispersion and increasing stability range. The optimization procedure proposed in the paper may be used to define the coefficients in other symplectic or non-symplectic schemes. Though striking numerical results are obtained by our scheme with a finite difference operator, a low-dispersive spatial operator is still needed to be combined with our scheme to form a more powerful tool for wave simulations.
Key words: Elastic wave equation     RKN     Dispersion relation     Stability condition     Finite difference method    
 1 引言

大尺度高精度地震波模拟与成像迫切需要发展高效高精度的地震波正演模拟方法(Etgen et al., 2009; Virieux and Operto, 2009). 对于地震波方程的直接求解而言,空间离散已经发展了众多方法,如有限差分法(Alford et al., 1974; Virieux, 19841986; Levander, 1988; Heiner et al., 1995; 董良国等,2000; Moczo et al., 2000; 王秀明等,2003; Zhang and Yao, 2013; Liu,2014),伪谱法(Gazdag,1981; Wang et al., 2013; Long et al., 2013),有限元法(Marfurt,1984; Hughes,1987; Padovani et al., 1994; 杨顶辉等,2002; 张怀等,2009; 李伟华等,2010; 刘有山等,2013),谱元法(Patera 1984; Seriani and Priolo, 1994; Komatitsch and Vilotte, 1998; Komatitsch et al., 2002; De Basabe and Sen, 20072010; 汪文帅等,2012),以及Yang等发展的高效高精度近似解析离散方法(NAD,Nearly analytic discrete method)(Yang et al., 20032009201020122014)等. 但是适用于地震波方程的时间离散方法的研究相对较少.

时间离散使用最多的是二阶中心差分(Virieux, 19841986),在所有时间离散方法中此格式效率最高. 但中心差分精度较低,由时间离散引入的数值频散极其明显,为了减弱时间频散,需采用远小于稳定限(库朗数能取到的最大值)的时间步长(Liu et al., 2014; Tan and Huang, 2014),这必然造成计算量的成倍增加. 为了在时间精度和计算效率之 间优化平衡,Dablain(1986)首次将Lax与Wendroff(1964)提出的空间微分修正时间精度的Lax-Wendroff方法运用于声波的模拟中,并获得四阶时间精度. Lax-Wendroff方法需求解高阶空间微分,而空间网格有限个节点的波场值用差分的方式近似高阶微分精度往往较低,利用伪谱法近似(Chen,2009; Long et al,2013)精度较高,但伪谱法计算量较大. 最近Tong等(2013)利用NAD算子并考虑对角线上节点信息,发展了时间四阶空间八阶 的两步STEM方法(Stereo-modeling method). Liu和Sen(2013)选取空间有限差分模板(Finite difference stencil),借助声波方程的频率波数域平面波解,用Taylor展开对比系数,构造了时空任意偶数阶的有限差分模拟方法,Tan和Huang(2014)将此方法用于交错网格的声波模拟中,并证明基于频散关系的有限差分系数确定方法与Lax-Wendroff方法等价. Lax-Wendroff方法亦可以结合有限元法或谱元法获得任意偶数阶时间精度( De Basabe and Sen, 2010).

计算效率的提高除了提高时间精度以外,还可以提高方法的稳定性,使时间步长可以取较大值. 虽然隐式格式不受库朗数的限制,可以使用相对较大的时间步长,但隐式格式需要求解方程组,这必然增加计算量. 罗明秋等(20012003a2003b)对Pade近似的隐式二阶辛算法做了研究. 隐式辛算法涉及到空间差分算子的求逆,若运用LU分解精度高但计算量较大且消耗巨大计算机内存,谱因式分解计算量小却精度低,LU分解与谱因式分解相结合的混合方法虽然吸收两者的优势,但隐式格式仍面临着边界不易处理、边界或速度间断面处误差过大等弊端. 为了获得隐式格式的高稳定性,同时兼顾显式格式的计算效率,Yang等(200920102012)做了许多有益的研究,发展了基于隐式Runge-Kutta(RK)与Adams方法的显式化方法,以及预估校正的隐式RK方法.

Chen(2007)讨论了多种显式高阶时间离散格式,包括Lax-Wendroff法、三级四阶辛RKN方法以及三阶的PRK方法(Partitioned Runge-Kutta method),并通过大量的数值实验证实了后两种辛方法的时间误差压制能力远优于Lax-Wendroff方法. 马啸等(2010)Ma等(2011)时间上采用二阶的PRK格式,空间上采用NAD算子,得到了声波和弹性波模拟的NSPKR方法,从数值频散和计算精度等方面说明了该方法的优势. Li等(20112012)注重数值格式长时间的计算能力,时间离散上采用三级三阶的辛PRK格式,发展了三阶辛褶 积微分算子地震波模拟方法. Nissen-Meyer等(2008)在空间离散使用谱元法,时间离散上用四阶的PRK格式,模拟了固液2D地球中地震波的传播.

虽然高阶辛格式具有较高精度和较强的长时间计算能力,但高阶格式在每一个时间步涉及到较多的子迭代步,计算量较大. 二阶辛格式计算量较小,但精度相对较低. 如何最大限度的提升二阶格式的计算精度和稳定性是本文研究重点. 本文采用显式的二阶辛RKN格式离散弹性波方程,用根数理论得到阶条件方程,通过给定限定方程得到辛系数的对称解,在频散分析中,为了使时间离散引入的数值频散最小,提出数值频率与不含时间离散的真实频率比较,得到一种最小频散辛RKN格式. 在稳定性分析中,设定辛系数使库朗数最大化,得到一种最稳定的辛RKN格式,发现在对称限定条件下,这两种辛RKN格式为同一种格式. 数值频散分析表明,新得到的辛RKN格式,频散特征与高阶格式接近; 稳定性分析表明,该格式稳定限大约为常规二阶格式的2倍,大约为三阶和四阶格式的1.5倍. 在数值实验中通过对其精度、效率、稳定性等方面测试,数值结果进一步证实了理论分析的 正确性和本文提出的时间积分格式求解弹性波的能力. 2 辛RKN格式

真实地球介质极其复杂,除了各向异性以外,可能还有多相多孔性以及黏性(Ma and Liu, 2006; Magnoni et al., 2014; Yang et al., 2014). 在区域尺度地震勘探,地震地面运动中,地球介质常常被简化为弹性介质. 本文主要研究弹性波方程的数值求解问题. 2.1 二级二阶显式辛RKN格式

在不考虑外力的情形下,3D非均匀各向异性介质中弹性波方程可以表示为

其中,ρ=ρ(x1,x2,x3)为介质密度,ui为第i个方向的位移分量,σij和εij分别为应力和应变张量,cijkl为四阶弹性张量. 引入变量vi=,位移和速度向量表示为 U =(u1,u2,u3)T和V =(v1,v2,v3)T,(1)式的第一个等式可以降阶为

考虑(2)式为Hamilton系统的情形,(2)式可以表示为

其中,H(U,V)为一标量函数,其具体形式为H(U,V)= .由(2)式和(3)式可知,f(U)可以表示为一标量函数的梯度,因此可以考虑(3)式的辛算法. 对于一般的s级RKN格式可以写为(Okunbor and Skeel, 19921994)

其中,Δt为时间步长,ci,aijj,bj为辛系数,下标n表示第n个时间层. 由 V n∧ U n= V n+1∧ U n+1(∧为外积,Wedge product)联合(4)式,容易证明(4)式为显式辛的,辛系数应满足如下等式(Okunbor and Skeel, 1992; 冯康和秦孟兆,2003)

在地震波模拟时,期望格式(4)的计算效率与二阶的PRK格式相近(马啸等,2010; Ma et al. 2011),当s=2时即满足高效地震波模拟的要求. 在精度上,期望格式(4)尽可能提高,对于非约化的辛RKN格式(即格式中无冗余级,bj≠0),满足(5)式的格式(4)是p阶的,可以根据用图形表示微分关系的根数理论(Hairer,2006; 冯康和秦孟兆,2003),即任何s-树,都有一个有根的s-树ρsτ由sτ的一胖节点提升得到,ρsτ的权与密度满足等式

其中,φ为ρsτ的权,γ为ρsτ的密度,权和密度的定义见文献(冯康和秦孟兆,2003). 如果二级的辛RKN格式是三阶的,(6)式中有4个非多余有根s-树需满足(6)式,将辛条件(5)式带入(6)式化简得到:

(7)式并无实数解,即二级的辛RKN格式只能达到二阶,(7)式的前两个等式需要满足. 求解方程(3)时常常期望获得满足(7)式的对称解(Okunbor and Skeel, 1994),即c1+c2=1,对于不可逆的动力过程,期望辛系数均为正实数(Li and Wu, 2011),二级二阶显式辛RKN格式辛系数满足的条件为

2.2 强稳定辛RKN格式

二级二阶显示辛RKN格式运用于弹性波方程的时间积分,时间与空间步长应满足一定的关系,否则计算失稳. 为了保持分析过程中量纲的一致,引入变量= V Δt,考虑如下的简谐解:

其中,ω为角频率,三个方向上的波数k1=ksinφcosφ,k2=ksinφsinφ,k3=kcosφ.φ为波数矢量 k 与x3轴之间的夹角,φ为 k 在x1x2坐标平面的投影与x1轴之间的夹角. 将(9)式带入(4)式可以得到t=(n+1)Δt和t=nΔt两个时刻的相满足的等式关系,联立(8)式即得到时间演进方程:

其中,M 为6×6映射矩阵(Jacobi矩阵),具体表示如下

其中,h为空间网格间距,e1=b1c21+b1c22,e2=b1b2(c2-c1),I为3×3的单位矩阵,Σ 是关于空间离散的3×3的矩阵.(11)式中 M 为6阶的矩阵,在空间离散未知的情况下,不可能直接求解映射矩阵得到具体特征值. 在这里引入相似变换(Igel et al., 1995; 董良国等,2000; Ma et al., 2014),假设存在可逆矩阵P Σ对角化,即有以下等式

其中,λi为 Σ 的特征值,且均为负数. M 的相似矩阵为

其中,的形式是 M 中包含的 Σ 矩阵换成其相似对角矩阵 Ω .求解的特征值相对容易,令θi=.假设 M 的特征值为ξ,由(13)式容易得到ξ满足的方程为

为了使(10)式稳定的充要条件为 ξ ≤1(De Basabe and Sen, 20072010). 当e2(1+c1-c2)= 时,θi能取到极小值,且为-16,联立(8)式,一种最稳定辛RKN格式系数为

2.3 弱频散辛RKN格式

由于用时间和空间的离散点值近似替代时间空间连续变化函数,必然会引起真实信号中的部分信号无法拾取而导致数值频散. 在构造数值格式时,时间离散对频散的讨论较少(Liu and Sen, 2013),因此构造弱频散的时间格式压制数值频散极其重要. 由(10)式和(14)式知,数值频率可以表示为

将(9)式带入(1)式知,仅有空间离散,无时间离散的真实角频率为

由(16)式和(17))式,将数值频率与真实频率比较,并对反正弦函数进行Taylor展开,令w=e2(1+c1-c2),有以下等式:

对(18)式中的再次利用Taylor展开得到等式

由(19)式知,当w= 时,数值频率与真实频率的比值为θi的高阶无穷小. 联立(8)式知方程无解,但可以得到限定条件|w- |→min的最优解,一种最小频散辛RKN格式为

由(15)式与(20)式知,最稳定格式和最小频散格式在对称限定条件下为同一种格式. 3 辛RKN格式特性分析

在空间离散时,(11)式并未指定空间离散格式的具体形式,因此新得到的辛RKN格式具有一般意义,即不依赖于空间离散格式. 3D空间离散采用有限差分、伪谱法时得到3个不同特征值,如果采用有限元或谱元法离散时将得到众多特征值(刘少林等,2014). 为了分析的需要,设定介质为各向同性弹性介质,采用有限差分或伪谱法离散空间微分算子. 为了对比,记本文格式为M1,二阶的辛PRK格式(马啸等,2010; Ma et al., 2011)为M2,二阶优化的辛PRK格式(刘少林等,2013)为M3,三级三阶的PRK格式(Li et al., 20112012)为M4,三级四阶的RKN格式(Okumbor and Skeel, 1992; Chen,2007)为M5. 3.1 稳定性分析

得到空间算子各阶有限差分的稳定性解析表达式较困难,但可以借助数值手段得到其近似显示表达式(Ma et al., 2011),将稳定性条件设定为

其中,VP和VS分别为P波和S波速度,a和b为待定常数.让波速比r=VP/VS的变化范围为,步长为0.01,借助软件(如Matlab)由(11—14)式计算得到每一个r对应的库朗数取到的最大值,对这859个点用曲线拟合确定(21)式中的常数. 在1D-3D空间中,不同差分阶的稳定性参数如表 1所示,2D空间差分阶数为8阶时M1-M5稳定性参数如表 2所示.
表 1 空间为不同差分阶数的有限差分算子和时间为M1求解1D-3D弹性波方程的稳定性参数 Table 1 Stability parameters of M1 with various order central difference operators for solving 1D to 3D elastic wave equations

表 2 空间为八阶的有限差分算子和时间为M1-M5求解2D弹性波方程的稳定性参数 Table 2 Stability parameters of M1-M5 with an eighth-order central difference operator for solving 2D elastic wave equation

表 1可以看出,随着空间精度的增加,稳定限逐渐变小,波速比对稳定的影响逐渐减弱(a逐渐减小),稳定限大致与(m为空间维数)成反比. 由表 2可知,M1-M5的稳定限对波速比的依赖大致相同,a的取值在0.44左右. M1的稳定限分别为M2-M5的2、1.77、1.50、1.55倍,这意味着在空间离散相同时,M1能取最大的时间步长. 3.2 频散分析

在2D均匀各向同性弹性情形下(11)式中 M 有两个不同特征值,分别对应P波和S波,在相同频 率下,P波的波长要大于S波,因此在相同空间网格的情形下,一个波长内P波的采样点数要大于S波,同时值得注意的是,在相同采样率条件下S波的频散往往大于P波(De Basabe and Sen, 20072010; 刘少林等,2014). 因此频散分析中,对S波的频散分析更有意义.

由(10)式和(14)式知S的频散关系可以表示为

其中,s为采样率(一个波长采样点数的倒数),θS=max1,λ2).图 1给出了空间为八阶中心差分,波数比为2,库朗数为0.5时,M1-M5的φ∈[0,π/4]四个方向上S波的频散曲线.
图 1 M1-M5时间离散弹性波方程的S波频散曲线
(a),(c),(e)和(g)分别为S波沿φ=0,π/12,π/6和π/4传播时的频散曲线. 右边图为左边对应图的局部放大图. 蓝色、红色、青色、品红色和黄色线分别对应于M1-M5的频散曲线.
Fig. 1 Dispersion curves of S wave using an eighth-order central difference operator and different temporal schemes M1-M5
The left plots(a),(c),(e) and (g)are dispersion curves of wave propagation with propagation angles 0,π/12,π/6 and π/4,respectively. Blue,red,cyan,magenta, and yellow curves are generated by M1,M2,M3,M4, and M5,respectively. The right plots are the local views of the corresponding plots in the left panel.

图 1可以看出,沿不同传播方向,对于5种时间离散方法,数值频散变化趋势相同,在φ=π/12方向上频散较小,在对角线方向上频散最大. 除了图 1a,当s>0.3时M1频散略大于M2和M3以外,其他情形M1频散均小于M2和M3,而略大于M4和M5,由于M4与M5频散足够接近,两者几乎完全覆盖,图中只能看到4条频散曲线,图 1(b、d、f、h)显示当一个波长的采样点数大于4时M1-M5 频散误差分别在0.0064,0.011,0.0098,0.0034和0.0034以内. 总体而言M1频散误差小于同阶格式M2和M3,而大于高阶格式M4和M5. 4 数值实验 4.1 2D均匀介质模型

为了验证这本文格式对弹性波的模拟精度和效率,设计了如图 2所示的2D均匀介质模型,模型大小为2000 m×2000 m,集中力源位于模型中心,由主频为25 Hz的Ricker子波激发. 两个接收点R1和R2分别位于(1300m,1100m)和(1300 m,800 m)处,空间网格间距为5 m,时间步长为0.001 s. P波速度为2000 m·s-1,S波速度为1154.7 m·s-1. 定义接收点的总体误差E=|ui-ur|,其中,ui为t=t时刻的数值解,ur为解析解,Nt为总时间步. 为了防止人工边界反射,模型周围为20个网格的PML吸收层. 模拟均在Intel Xeon X5650(六核,2.66G Hz)CPU上进行. 计算结果如图 3,效率与总体误差如表 3所示.

图 2 2D均匀介质模型
震源位于模型中心,两个接收点R1和R2分别位于(1300 m,1100 m)和(1300 m,800 m)处.
Fig. 2 A 2D homogeneous model
A vertical force source is located at the central of the model; two receivers R1 and R2 are located at(1300 m,1100 m) and (1300 m,800 m),respectively

图 3 两个接收点R1(a)和R2(b)处位移垂直分量的波形对比及放大图(c,d)
黑色线为解析解,蓝色,红色,青色,品红和黄色线分别由M1-M5计算得到.(c)和(d)分别为(a)和(b)矩形框中波形放大图.
Fig. 3 Waveforms of vertical components of displacement at R1(a) and R2(b)generated by analytic solution and numerical methods
(c—d)are local views of(a—b),respectively. The black lines are calculated by analytic solution,which is obtained by Cagniard-De Hoop method(De Hoop,1960), and the blue,red,cyan,magenta, and yellow lines are calculated by M1-M5,respectively.

表 3 5种数值算法的计算效率与总体误差对比 Table 3 Computational efficiency of the five numerical schemes and the total error at R1 and R2

图 3为数值解与解析解(De Hoop,1960)的位移垂直分量波形对比,图 3a图 3b显示数值解与解析解总体上接近,即M1-M5均能较好模拟弹性波传播,但局部放大图,图 3c图 3d显示M4和M5与解析解最接近,M1与解析解的接近程度好于M2和M3而略逊于M4和M5. 表 3给出了5种格式的计算效率与总体误差的对比,就内存要求而言,5种方法内存要求基本相同,就计算时间而言,M1比最快的M2计算时间多出了5.46%,M4和M5的计算时间接近,M1比M5的计算时间减少了29.79%; 从R1和R2的总体误差E1和E2上看,M2的误差几乎为M1的2倍,M3的误差为M1的1.21倍,而M1的误差为M4的1.35倍. 综合考虑计算效率和精度发现,M1是一种较为高效方法. 值得注意的是,四阶的M5的精度略低于三阶的 M4,这可能由于M5本身的频散特征或其他未知原因造成(Liu et al., 2015). 通过分配额外的位移 中间数组可以减少M1,M3,M4和M5的计算时间(Liu et al., 2014),但实际计算表明,计算时间并不明显减少,反而内存显著增加,因此在利用显式辛算法模拟地震波传播时,建议广义坐标与动量的更新单独进行. 4.2 非均匀介质模型-Marmousi模型

为了检验M1在强烈非均匀介质中模型弹性波传播的稳定性,选取Marmousi做测试. 图 4为Marmousi模型的P波速度结构,速度变化范围为1500~5500 m·s-1,S波速度为VS=VP/3,密度由Cardner公式ρ=310×V0.25P得到. 模型由384×122个网格点组成,网格间距为10 m. 爆炸源位于模型中心,由主频为20 Hz的Ricker激发. 一个接收点位于(2060 m,750 m)处. 除了选取M1计算外,另外选取M2和M4做对比. 由表 2可知,对于此模型,M1,M2和M4能取到的最大时间步长分别为0.00247 s、0.00123 s和0.00165 s. 模拟结果如图 5图 6所示.

图 4 Marmousi模型P波速度结构Fig. 4 The velocity profile of P-wave in Marmousi model

图 5为M1,M2和M4三种方法模拟Marmousi 模型中弹性波传播位移水平分量0.2 s时的波场快照,图 5(a—c)中三种方法采用的时间步长均为 0.001 s,图 5d中M1采用的时间步长为0.002 s. 总体上看波场模拟结果相似. 当M2采用0.0013 s的时间步长,M4采用0.0017 s的时间步长时,M2和M4的数值结果溢出,计算无法进行,说明M2和M4已经失稳. 从数值稳定性上看,M1稳定性 远好于其他方法. 图 6为M1采用0.002 s的时间 步长与M4采用0.001 s的时间步长得到的合成地震记录对比,由于M4为高阶格式且采用较小的时间步长,因此M4结果可做为参考解. 从图 6对比可以看出,M1在主要震相处与参考解拟合较好,只在尾波部分有微弱偏差,说明M1采用大时间步长的模拟结果仍是可信的.

图 5 M1、M2和M4三种方法模型Marmousi模型中弹性波传播0.2 s时的位移水平分量波场快照
(a—c)三种方法采用0.001 s的时间步长;(d)M1采用0.002 s的时间步长.
Fig. 5 Snapshots of wavefield in Marmousi model at t=0.2 s
(a—c)are respectively computed by M1,M2 and M4 with the same time step Δt=0.001 s;(d)is computed by M1 with a doubled time increment.

图 6 M1用0.002 s的时间步长(点线)和M4用0.001 s的时间步长(实线)得到的地震记录对比
(a)和(b)分别为位移的水平分量和垂直分量.
Fig. 6 Synthetic seismic records obtained by M1 with a time step 0.002 s(dotted line) and M4 with time step 0.001 s(solid line)
(a)is the horizontal component of displacement and (b)is the vertical component of displacement. Two curves in each plots are well superposed with each other except very small deviations in the part of coda wave.
5 讨论

本文讨论了弹性波模拟的显式二阶辛RKN格式,得到了一种不依赖于空间离散的优化格式. 在频散分析中,发现新得到的格式频散误差小于常见的二阶格式和优化二阶格式,稳定限约为二阶格式的2倍和三阶与四阶格式的1.5倍. 数值实验进一步证实了本文格式在计算精度、效率以及稳定性等方面的优越性,同时也说明本文提出的频率误差最小化思想用于时间离散优化的正确性. 在大尺度地震波模拟与成像中,可以使用本文格式采用较大的时间步长节省计算时间同时保持较高精度. 本文空间离散上仅采用有限差分,对于复杂起伏地表模型,有限差分法将受到局限,虽然空间离散采用有限元或谱元法将会得到类似的结论,但仍需要定量讨论. 空间上采用其他弱频散离散方法,以及高阶正系数辛算法的讨论,我们将另文给出.

致谢 感谢与西北工业大学应用数学系马啸老师的有益讨论.
参考文献
[1] 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.
[2] Chen J B. 2007. High-order time discretizations in seismic modeling. Geophysics, 72(5): SM115-SM122.
[3] Chen J B. 2009. Lax-Wendroff and Nystrm methods for seismic modelling. Geophys. Prospect., 57(6): 931-941.
[4] Dablain M A. 1986. The application of high-order differencing to the scalar wave equation. Geophysics, 51(1): 54-66.
[5] De Basabe J D, Sen M K. 2007. Grid dispersion and stability criteria of some common finite-element methods for acoustic and elastic wave equation. Geophysics, 72(6): T81-T95.
[6] De Basabe J D, Sen M K. 2010. Stability of the high-order finite elements for acoustic or elastic wave propagation with high-order time stepping. Geophys. J. Int., 181(1): 577-590.
[7] De Hoop A T. 1960. A modification of Cagniard's method for solving seismic pulse problems. Appl. Sci. Res., 8(1): 349-356.
[8] Dong L G, Ma Z T, Gao J Z. 2000. A study on stability of the staggered-grid high-order difference method of first-order elastic wave equation. Chinese J. Geophys. (in Chinese), 43(6): 856-864.
[9] Etgen J, Gray S H, Zhang Y. 2009. An overview of depth imaging in exploration geophysics. Geophysics, 74(6): WCA5-WCA17.
[10] Feng K, Qin M Z. 2003. Symplectic Geometric Algorithm for Hamiltonian Systems (in Chinese). Hangzhou: Zhejiang Science & Technology Press.
[11] Gazdag J. 1981. Modeling of the acoustic wave equation with transform methods. Geophysics, 46(6): 854-859.
[12] Hairer E. 2006. Geometric Numerical Integration Ⅰ (2nd ed). Berlin and New York: Springer-Verlag.
[13] Heiner I, Mora P, Riollet B. 1995. Anisotropic wave propagation through finite-difference grids. Geophysics, 60(4): 1203-1216.
[14] Hughes T J R. 1987. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. New Jersey: Prentice-Hall.
[15] Igel H, Mora P, Riollet B. 1995. Anisotropic wave propagation through finite-difference grids. Geophysics, 60(4): 1203-1216.
[16] 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. Seismiol. Soc. Am., 88(2): 368-392.
[17] Komatitsch D, Ritsema J, Tromp J. 2002. The spectral-element method, Beowulf computing, and global seismology. Science, 298(5599): 1737-1742.
[18] Lax P D, Wendroff B. 1964. Difference schemes for hyperbolic equations with high order of accuracy. Commun. Pure. Appl. Math., 17(3): 381-398.
[19] Levander A R. 1988. Fourth-order finite-difference P-SV seismograms. Geophysics, 53(11): 1425-1436.
[20] Li W H, Liu Q H, Zhao C G. 2010. Three-dimensional viscous-spring boundaries in time domain and dynamic analysis using explicit finite element method of saturated porous medium. Chinese J. Geophys. (in Chinese), 53(10): 2460-2469, doi: 10.3969/j.issn.0001-5733.2010.10.020.
[21] Li R, Wu X. 2011. Two new fourth-order three-stage sympletic integrators. Chinese Phys. Lett., 28(7): 070201, doi: 10.1088/0256-307X/28/7/070201.
[22] Li X F, Li Y Q, Zhang M G, et al. 2011. Scalar seismic-wave equation modeling by a multisymplectic discrete singular convolution differentiator method. Bull. Seism. Am. Soc., 101(4): 1710-1718.
[23] Li X F, Wang W S, Lu M W, et al. 2012. Structure-preserving modeling of elastic waves: a symplectic discrete singular convolution differentiator method. Geophys. J. Int., 188(3): 1382-1392.
[24] Liu H F, Dai N X, Niu F L, et al. 2014. An explicit time evolution method for acoustic wave propagation. Geophysics, 79(3): T117-T124.
[25] 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.
[26] Liu S L, Li X F, Wang W S, et al. 2014. A new kind of optimal second-order symplectic scheme for seismic wave simulations. Science China: Earth Sciences, 57(4): 751-758.
[27] Liu S L, Li X F, Liu Y S, et al. 2014. Dispersion analysis of triangle-based finite element method for acoustic and elastic wave simulations. Chinese J. Geophys. (in Chinese), 57 (8): 2620-2630, doi: 10.6038/cjg20140821.
[28] Liu S L, Li X F, Wang W S, et al. 2013. Optimal symplectic scheme and generalized discrete convolutional differentiator for seismic wave modeling. Chinese J. Geophys. (in Chinese), 56(7): 2452-2462, doi: 10.6038/cjg20130731.
[29] Liu Y, Sen M K. 2013. Time-space domain dispersion-relation-based finite-difference method with arbitrary even-order accuracy for the 2D acoustic wave equation. J. Comput. Phys., 232(1): 327-345.
[30] Liu Y. 2014. Optimal staggered-grid finite-difference schemes based on least-squares for wave equation modelling. Geophys. J. Int., 197(2): 1033-1047.
[31] 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., (in Chinese), 56(9): 3085-3099, doi: 10.6038/cjg20130921.
[32] Long G H, Zhao Y B, Zhou J. 2013. A temporal fourth-order scheme for the first-order acoustic wave equations. Geophys. J. Int., 194(3): 1473-1485.
[33] Luo M Q, Liu H, Li Y M. 2001. Seismic wave modeling with implicit symplectic method based on spectral factorization on helix. Chinese J. Geophys., (in Chinese), 44(3): 379-388.
[34] Luo M Q, Liu H, Li Y M. 2003a. LU decomposition with spectral factorization in seismic imaging. Chinese J. Geophys., (in Chinese), 46(3): 421-427.
[35] Luo M Q, Zhu G T, Liu H, et al. 2003b. A hybrid method of matrix inversion suited for 3D implicit prestack depth migration. Chinese J. Geophys., (in Chinese), 46(5): 684-689.
[36] Ma S, Liu P C. 2006. Modeling of the perfectly matched layer absorbing boundaries and intrinsic attenuation in explicit finite-element methods. Bull. Seismiol. Am. Soc., 96(5): 1779-1794.
[37] Ma X, Yang D H, Zhang J H. 2010. Symplectic partitioned Runge-Kutta method for solving the acoustic wave equation. Chinese J. Geophys.(in Chinese), 53(8): 1993-2003, doi: 10.3969/j.issn.0001-5733.2010.08.026.
[38] Ma X, Yang D H, Liu F Q. 2011. A nearly analytic symplectically partitioned Runge-Kutta method for 2-D seismic wave equations. Geophys. ,J. Int. 187(1): 480-496.
[39] Ma X, Yang D H, Song G J, et al. 2014. A low-dispersive symplectic partitioned Rung-Kutta method for solving seismic-wave equations: Ⅰ. scheme and theoretical analysis. Bull. Seismol. Soc. Am., 104(5), 2206-2225.
[40] Marfurt K. 1984. Accuracy of finite-difference and finite-element modeling of scalar and elastic wave equations. Geophysics, 49(5): 533-549.
[41] Magnoni F, Casarotti E, Michelini A, et al. 2014. Spectral-element simulations of seismic waves generated by the 2009 L′Aquila earthquake. Bull. Seismiol. Am. Soc., 104(1): 73-94.
[42] Moczo P, Kristek J, Halada L. 2000. 3D fourth-order staggered-grid finite-difference scheme: stability and grid dispersion. Bull. Seismiol. Soc. Am., 90(3): 587-603.
[43] Nissen-Meyer T, Fournier A, Dahlen F A. 2008. A 2-D spectral-element method for computing spherical-earth seismograms-Ⅱ. Waves in solid-fluid media. Geophys. J. Int., 174(3): 873-888.
[44] Okumbor D, Skeel R D. 1992. Explicit canonical methods for Hamiltonian systems. Math. Comp., 59(200): 439-455.
[45] Okunbor D, Skeel R D. 1994. Canonical Runge-Kutta-Nystrm methods of orders five and six. J. Comp. Appl. Math., 51(3): 375-382.
[46] Padovani E, Priolo E, Seriani G. 1994. Low and high order finite element method: experience in seismic modeling. J. Comp. Acoust., 2(4): 371-422.
[47] Patera A T. 1984. A spectral element method for fluid dynamics: laminar flow in a channel expansion. J. Comp. Phys., 54(3): 468-488.
[48] Seriani G, Priolo E. 1994. Spectral element method for acoustic wave simulation in heterogeneous media. Finite Elem. Anal. Des., 16(3-4): 337-348.
[49] Tan S R, Huang L J. 2014. An efficient finite-difference method with high-order accuracy in both time and space domains for modelling scalar-wave propagation. Geophys. J. Int., 197(2): 1250-1267.
[50] Tong P, Yang D H, Hua B L, et al. 2013. A high-order stereo-modeling method for solving wave equations. Bull. Seismiol. Am. Soc., 103(2): 811-833.
[51] Virieux J, Operto S. 2009. An overview of full-waveform inversion in exploration geophysics. Geophysics, 74(6): WCC1-WCC26.
[52] Virieux J. 1984. SH-wave propagation in heterogeneous media: velocity-stress finite-difference method. Geophysics, 49(11): 1933-1942.
[53] Virieux J. 1986. P-SV wave propagation in heterogeneous media: velocity-stress finite-difference method. Geophysics, 51(4): 889-901.
[54] 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. (in Chinese), 55(10): 3427-3439, doi:10.6038/j.issn.0001-5733.2012.10.026.
[55] Wang X M, Zhang H L, Wang D. 2003. Modelling of seismic wave propagation in heterogeneous poroelastic media using a high-order staggered finite-difference method. Chinese J. Geophys. (in Chinese), 46(6): 842-849.
[56] Wang Y B, Takenaka H, Jiang X H, et al. 2013. Modelling two-dimensional global seismic wave propagation in a laterally heterogeneous whole-Moon model. Geophys. J. Int., 192(3): 1271-1287.
[57] Yang D H. 2002. Finite element method of the elastic wave equation and wavefield simulation in two-phase anisotropic media. Chinese J. Geophys. (in Chinese), 45(4): 575-583.
[58] 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. Seismiol. ,Soc. Am. 93(2): 882-890.
[59] Yang D H, Wang N, Chen S, et al. 2009. An explicit method based on the implicit Runge-Kutta algorithm for solving wave equations. Bull. Seismiol. Soc. Am., 99(6): 3340-3354.
[60] Yang D H, Wang L, Deng X Y. 2010. An explicit split-step algorithm of the implicit Adams method for solving 2-D acoustic and elastic wave equations. Geophys. J. Int., 180(1): 291-310.
[61] Yang D H, Wang N, Liu E. 2012. A strong stability-preserving predictor-corrector method for the simulation of elastic wave propagation in anisotropic media. Commun. Comput. Phys., 12(4): 1006-1032.
[62] Yang D H, Wang M X, Ma X. 2014. Symplectic steromodelling method for solving elastic wave equations in porous media. Geophys. J. Int., 196(1): 560-579.
[63] Zhang H, Zhou Y Z, Wu Z L, et al. 2009. Finite element analysis of seismic wave propagation characteristics in Fuzhou basin. Chinese J. Geophys. (in Chinese), 52(5): 1270-1279, a href="http://dx.doi.org/10.3969/j.issn.0001-5733.2009.05.016 " target=_blank>doi: 10.3969/j.issn.0001-5733.2009.05.016.
[64] Zhang J H, Yao Z X. 2013. Optimized explicit finite-difference schemes for spatial derivatives using maximum norm. J. Comput. Phys., 250(1): 511-526.
[65] 董良国, 马在田, 曹景忠. 2000. 一阶弹性波方程交错网格高阶差分解法稳定性研究. 地球物理学报, 43(6): 856-864.
[66] 冯康, 秦孟兆. 2003. 哈密尔顿系统的辛几何算法. 杭州: 浙江科学出版社.
[67] 李伟华, 刘清华, 赵成刚. 2010. 饱和多孔介质三维时域黏弹性人工边界与动力学反应分析的显式有限元法. 地球物理学报,53(10): 2460-2469, doi: 10.3969/j.issn.0001-5733.2010.10.020.
[68] 刘少林, 李小凡, 刘有山等. 2014. 三角网格有限元法声波与弹性波模拟频散分析. 地球物理学报, 57(8): 2620-2630,doi: 10.6038/cjg20140821.
[69] 刘少林, 李小凡, 汪文帅等. 2013. 最优化辛格式广义离散奇异核褶积微分算子地震波场模拟. 地球物理学报, 56(7): 2452-2462, doi: 10.6038/cjg20130731.
[70] 刘有山, 滕吉文, 刘少林等. 2013. 稀疏存储的显式有限元三角网格地震波数值模拟及其PML吸收边界. 地球物理学报, 2013, 56(9): 3085-3099, doi: 10.6038/cjg20130921.
[71] 罗明秋, 刘洪, 李幼铭. 2001. 基于螺旋线上谱因式分解的地震波场隐式辛算法. 地球物理学报, 44(3): 379-388.
[72] 罗明秋, 刘洪, 李幼铭. 2003a. 用于波场成像的谱法LU分解. 地球物理学报, 46(3): 421-427.
[73] 罗明秋, 朱国同, 刘洪等. 2003b. 适用于三维隐式叠前深度偏移中矩阵求逆的混合算法. 地球物理学报, 46(5): 684-689.
[74] 马啸, 杨顶辉, 张锦华. 2010. 求解声波方程的辛可分Runge-Kutta方法. 地球物理学报, 53(8): 1993-2003, doi: 10.3969/j.issn.0001-5733.2010.08.026.
[75] 汪文帅, 李小凡, 鲁明文等. 2012. 基于多辛结构谱元法的保结构地震波场模拟. 地球物理学报, 55(10): 3427-3439, doi:10.6038/j.issn.0001-5733.2012.10.026.
[76] 王秀明, 张海澜, 王东. 2003. 利用高阶交错网格有限差分法模拟地震波在非均匀孔隙介质中的传播.地球物理学报, 46(6): 842-849.
[77] 杨顶辉. 2002. 双相各向异性介质中弹性波方程的有限元解法及波场模拟. 地球物理学报, 45(4): 575-583.
[78] 张怀, 周元泽, 吴忠良等. 2009. 福州盆地强地面运动特征的有限元法数值模拟. 地球物理学报, 52(2): 1270-1279, doi: 10.3969/j.issn.0001-5733.2009.05.016.