2012, Vol. 55
2. 宁夏大学数学计算机学院, 银川 750021
2. School of Mathematics and Computer Science, Ningxia University, Yinchuan 750021, China
地震波数值模拟是研究地震波传播特征的有效手段，有助于认识复杂介质中地震波传播规律，是地震勘探和地震学研究的重要组成部分.20世纪60、70年代以来，随着地震波理论和计算机技术的不断发展，各种数值模拟方法应运而生.就波动方程模拟而言，目前常见的正演数值模拟方法有有限差分法、有限元法、伪谱法和谱元法等.
经典有限差分法是偏微分方程的主要数值解法之一，在各种地震波数值模拟方法中，也是最早出现、现在运用最广泛的方法，其特点是计算速度快，计算效率高.Alterman&Karal[1]首先将经典有限差分法应用到层状介质弹性波传播的数值模拟中.此后，通过Boore[2]，Madariaga[3]和Virieux[4]等的努力，在精度、应用于复杂介质等方面做了深入推广；Mora[5]，Frankel[6]，Frankel[7]，Pitarka[8]，Olsen[9]和Zingg[10]等在有限差分法的三维模拟方面做了大量的工作，尽管这些工作对有限差分法在精度、紧致性和优化等方面做了有益的改进，但由于方法本身的理论基础是Taylor展开的，不可避免地存在如模拟面波时精度相对较低[11]、对不规则几何体(尤其是自由起伏表面)的处理等方面显得不足等局限.
有限元法(FEM)的出现是数值分析方法研究领域内的重大突破性进展，其理论基础是变分原理，构建在Sobolev空间上，通过将单元上的分片插值函数和最小位能理论相结合来求解整个区域的波场信息.由于网格剖分的任意性及它所依据的变分原理，可以模拟任意地质体的形态，因此有限元自然被运用到地震波模拟中；虽然有限单元法成功地模拟了很多波的传播问题，但同样存在许多不足之处.如为使结果达到可接受的精度，一般说来低阶的有限单元每个波长需要至少布置10个节点，由此导致计算时需要的内存较大，耗时较多，计算效率较低；此外，低阶有限元法传播特性较差，频散误差较大[12]，而高阶有限元法则会产生一些问题，如伪波的出现[13].同时有限元离散后需要求解大型的线性稀疏矩阵方程组，以致有限元法在大规模弹性波方程模拟领域一直没有得到很好的应用.
伪谱法可以看作是有限差分法的推广，该算法在求解方程时，通过对空间坐标的偏微分实施快速傅氏变换，避免了对空间坐标的差分运算，从而大大加快了计算速度并节省了大量的计算机内存.该方法在原问题的解是充分光滑的情况下，可以得到无穷阶收敛.谱方法和伪谱方法已经运用到了局部[14-15]和全球[16-17]地震波模拟中，可以在每个波长中用少量的网格点，且数值频散小，其表现优于有限差分法.但是和有限差分法一样，模拟面波时精度较低[11]，且由于伪谱法是全局算子的特性，并行性较差，又网格点的位置随着基函数的选择是固定的.所以，该方法限于光滑规则的几何介质，不能运用到复杂几何形状的介质中或者变化剧烈的非均匀介质中[18].
谱元法(SEM)[19-30]建立在波动方程的变分或弱形式理论基础上，它兼顾了伪谱方法具有全局高精度和有限元方法在处理区域分解方面的灵活性，被广泛运用到波动模拟中.在传统的FEM 中，一般选择低阶精度的基函数和高斯积分，而在SEM 中，选择高阶精度的Lagrange 插值基函数和Gauss-Lobatto-Legendre(GLL)积分可以得到对角化质量矩阵的显式时间域微分方程组，这给时间域的求解带来了便利，节约了耗时，并使得算法便于并行化.
因此鉴于谱元法在地震波模拟中的优势本文选用谱元法进行空间离散.但是随着模拟技术的发展和计算机处理数据能力的增加，对于地震波的模拟不再局限于短时的波场值了，对于某些波场(如地球自由振荡等)的长时程传播的计算及模拟问题成为热点，这就需要选择较为高效、稳定和具有较长时间模拟能力的数值算法.在“数值算法应尽可能多地保持原问题的本质特征"的指导原则下，冯康首先提出了保结构算法的思想[31]，并且引起了国内外学者的极大兴趣.保结构算法的基础是哈密尔顿体系，一切真实的、耗散可忽略不计的物理过程都可表示成哈密尔顿体系.研究表明离散型哈氏算法的体系结构与守恒律完整并行[31]，高度逼近于哈氏原型，而且拥有理论上无限长期的跟综能力，这就从根本上保证并说明了哈氏算法的独特性和优越性.怎样离散才能高度逼近哈氏原型呢?一条正确计算牛顿运动方程的途径就是先把方程哈密尔顿化，然后运用哈密尔顿算法求解使得问题原型的基本特征在离散后应该尽可能得到保持[31].因此离散化应尽可能在问题原型的同一形式框架中进行，则要求离散所用的差分格式应是辛的，即从上一时刻到下一时刻的映射(步进算子)是辛映射，这类方法称之为辛格式，相应的数值离散算法称为辛算法(见附录A).
有限元方法正是把离散纳入原型的解所组成的Sobolev空间的同一框架中进行，使得对称性、正定性和守恒性等基本特征得到保持[31]；钟万勰[32]证明了线性有限元质量矩阵对称即保辛，即有限元是自动保辛的，而谱元法是有限元的一种特殊情形[11]，也建立在变分原理的基础上，通过了特殊的插值基函数得到的质量矩阵为对角化的，所以谱元法也是自动保辛的；既然在空间域的离散格式保辛，正如王雨顺[33]指出:对哈密尔顿系统进行保结构离散，手段之一是首先进行空间离散，使得离散得到的常微分方程组是Hamilton 方程组，然后用辛算法求积.因此在时间域求解微分方程时，应该用保结构的辛算法.Iwatsu[34]基于分部的Runge-Kutta格式推导了一种新的三级三阶辛格式；李小凡等[35-36]将此格式结合褶积微分运用到声波及弹性波场的模拟中，取得了满意的效果，长时程运算时表现出其优势.我们的工作是寻求一类时间域离散的辛算法和谱元法结合，这样既能够保证时间域和空间域都具有保辛性，也能在计算耗时、存储量、稳定性等方面综合表现良好，即建立一种地震波时-空保结构模拟方法.本文推导了一种新的三级三阶辛算法(NTSTO)，首先通过与几种经典的三级三阶辛算法比较，说明新推导的三级三阶辛算法在长时程计算时其保结构能力优势明显；然后通过数值试验用建立的时-空保结构算法与将谱元法结合非辛的三阶Runge-Katta方法、Newmark 算法[37-38]和HHT-α算法[39-40]进行比较，试验表明，本文的保结构算法在计算时间、内存要求、稳定性和长时间跟踪能力方面都有上佳的表现；给出的起伏地表多层介质模型的数值算例验证了该算法处理复杂几何形状和复杂介质的有效性.该算法的推出，为谱元法的地震波模拟，尤其是如地球自由振荡等长时程模拟提供了新的较可靠的研究工具.
2 基本理论 2.1 弹性波方程通常可将地震波近似为弹性波.在非均匀各向异性介质空间Ω⊂R3 内的线弹性波动方程可表示为[24]
|
(1) |
初始条件
|
这里u为位移矢量，u0，$\dot{u}$0 分别为位移和速度的初值，σ 为对称的二阶应力张量，ε 为对称的二阶应变张量，ρ 表示为密度，f为震源项.符号上面的点代表微分，上标T 表示转置.对于三维各向异性介质，刚度张量有21个独立分量；对于二维完全各向异性介质，C中完全独立的分量个数为6 个，根据Hooke定律，
|
且有cij=cji，i，j= 1，2，…，6.如果退化为各向同性情况，则令c11 =c22 =c33 =λ+2μ，c12 =c13 =c23 =λ及c44 =c55 =c66 =μ，其中λ，μ为Lame常数.
2.2 弹性波动方程的弱积分表达式在谱元方法中，强形式的波动方程(1)首先要改写成变分或弱积分公式，通过对方程(1)两端点乘以任意一个测试向量函数w，然后在物理空间进行积分，即可得到波动方程的弱积分表达式:
|
(2) |
这里Γabs =Γa+Γf 为计算域Ω 的边界，Γa 与Γf 分别表示为人工边界和自由边界.在自由边界条件Γf上，σ·
|
(3) |
与有限元方法类似，谱元法的离散过程中首先将计算域Ω 剖分为Nel个互不重叠的单元Ωe，即有
|
(4) |
一般来说，谱元法对区域的剖分要求剖分为结构网格，如四边形(2D)或六面体(3D)单元，可以很方便地通过张量积来构造基函数，但对于模型较复杂时，不能够剖分为结构网格时，可剖分为非结构网格，如三角形(2D)[41-42]或四面体(3D)[43-44]，但需要特殊处理，如通过集合投影方式构建Fekete内插节点[45]，还需要特殊的投影算子等.这样一系列处理以后，事实上其原有谱方法高阶精度就很难保证了.研究也验证了二维三角形网格上的谱元法精度比四边形网格上的谱元法精度要低[20].
本文主要研究二维谱元法结合辛算法的地震波模拟，下面主要讨论二维的情况.二维区域剖分后，每个单元通过可逆的局部变换Fe:Λ ￫Ωe，将标准单元变换到一般单元，其中Λ = [-1,1]2 为单位正方形.在参照单元Λ 上引入Lagrange插值基函数，为了使最后所得的质量矩阵是对角的，选择插值节点与GLL 积分节点相同.GLL 积分节点为方程(1-ξ2)P′N(ξ)=0的(N+1)个根，记作ξi∈ [-1,1]，i=0，1，…，N，这里P′N(ξ)为N阶Legendre多项式的导数.用uNe，wNe分别表示为u和w的分量在参照单元Λ 上以ξi为插值节点的多项式逼近函数ueN，则
|
(5) |
其中hp(ξ)表示以ξi∈ [-1,1]，i=0，1，…，N为插值节点的一维Lagrange插值基函数，且具有性质hp(ξq)=δpq.在参照单元上，偏微分可以表示为
|
(6) |
这里h′ 为一维Lagrange插值函数的微分.
注意到(5)式，并利用GLL 积分技术，则有积分
|
(7) |
其中，J e为单元上的变换Fe 的Jacobian矩阵；GLL积分的权系数ωi>0是与每个单元无关的常数.则对问题(2)的离散变分问题可以表述为:
对任意的时间t，寻找uN，使得对任意的测试函数wN下式成立
|
(8) |
其中已引入记号ρeij=ρe(ξi，ηj)，we，Nij=we，N(ξi，ηj)等，且定义
|
|
利用公式(5)和(6)，将以上方程组进行整理，可得到时间域的二阶线性常微分方程组[24]:
|
(9) |
这里U表示为待求解的全局位移向量，Ütt表示关于时间t的二阶微分，M为质量矩阵，K为刚度矩阵，F为震源项.注意到上述方程组与有限元法、Chebyshev基函数SEM[30]不同的是，由于选择了特殊的插值基函数，且GLL 积分节点作为Lagrange插值节点，所以得到的质量矩阵M为对角化的[20].又谱元法的保辛性保证了问题原型的解与谱元离散的解在Sobolev空间的同一框架中，从而也保持了原有问题的结构.
2.4 时间离散先进行空间离散，使得离散得到的常微分方程组是Hamilton方程组，然后用辛算法求积，才可以保证整体离散结构的保辛性[33]，所以我们希望时间离散运用辛算法求积.自从冯康先生提出辛算法以来，发展了各种离散的辛格式，但是针对于谱元法的地震波模拟，并不是各类辛算法在地震波模拟中都能取得满意的效果，比如高阶精度格式一般都是多级格式，由于在一个时间间隔上进行多步运算，当步数较多时，必然降低了模拟效率；Iwatsu[34]建议，在模拟波的长时程传播时，考虑到计算耗时和精度问题，三阶的方法是一有效的方法.邢誉峰等[46]研究了显式辛算法的相位精度，分析得出辛算法相位精度与算法阶次具有不协调性，即辛算法的阶次高并不意味着其相位精度也高，指出三阶显式辛算法具有比较高的相位精度.所以针对谱元法和三阶辛算法的特点，我们选择显式的三阶辛算法作为时间域离散算法.下面我们推导一种新的三级三阶辛算法(NTSTO).
对于微分问题
|
式中点表示对时间的一阶微分，Vp(p)，Au(u)分别表示关于变量p，u的微分算子，则其三级三阶辛格式有
|
(10) |
其中ci，di，i=1，2，3 为满足辛格式条件[31, 47]的常数，Δt为时间间隔，根据相位误差最小原理[48]，选择系数如下:
|
(11) |
为了验证新推导的上述系数的优越性，选择常用的测试方程Fermi-Pasta-Ulam(FPU)ütt=-ω2u，初始条件为u(0)= 1/ω，$\dot{u}$ (0)= 1，求解步长Δt=0.001，ω=50，计算时间总长200s(20 万个时间步)，与经典的Ruth[49]、Iwatsu[34]、Suzuki[50]以及McLachlan和Atela (MLA)[51]的三级三阶辛格式进行比较；事实上，Suzuki格式与MLA 格式是从不同的途径推导的相同格式，所以比较时用其中一种格式即可.
图 1给出了四种三级三阶辛算法在20万个时间步内各周期内最大误差比较，可以看出，新推导的三级三阶辛算法(NTSTO)误差积累最小、稳定性好，适合长时程跟踪.
|
图 1 四种三级三阶辛算法求解FPU问题 200 s(20万个时间步）内的误差比较 Fig. 1 The residual comparison for FPU problem after 200 s (200 thousand time steps) generated by the four methods |
|
图 2 二维均匀介质模型 Fig. 2 2-D homogeneous medium model |
将辛格式(10)与谱元法进行结合，代入到空间域谱元法离散后的二阶线性常微分方程组(9)，进行时间域离散，引入变量P，使得MU·t=P，P·t=F-KU，则方程(1)的一个时间步的推进格式为
|
(12) |
为了验证方法的有效性，我们将本文提出的三级三阶辛算法结合谱元法(NTSOS-SEM)与以下几类算法进行比较:(1)Newmark 算法[37-38](见附录B)结合谱元法(Newmark-SEM):由于Newmark离散格式也是一种保能量算法，因此谱元法进行离散后运用Newmark 时间离散格式进行求解，在地震波模拟方面[11, 52-53]有了广泛的应用；(2)HHT-α 算法[39-40](见附录C)结合谱元法(HHT-α-SEM):谱元离散后时间域用HHT-α 算法进行求解，此算法是非辛算法；(3)非辛的三阶Runge-Kutta算法结合谱元法(RK3-SEM).
首先我们研究这几种算法与解析解的短时程对比.设计一均匀介质模型，模型大小为2550 m ×2550m，纵波波速Vp=3000 m/s，横波波速Vs=1732m/s，密度ρ=2400kg/m3，模型剖分为40×40个单元，每个单元上采用9个GLL 插值节点；震源(五角星所在位置)为位于(1280m，-1280m)处的集中力源，由主频为f0 =25 Hz的Ricker子波f(t)= [1-2(πf0(t-t0))2]exp(-(πf0(t-t0))2)激发产生，其中t0=0.03；时间步长为1ms；接收点(三角形所在的位置)为(1600m，-2000 m).上边界为自由边界条件，下边界、左右边界都用吸收边界条件[54].模型的精确解由DeHoop[55]给出.
图 3给出了四种算法在0.6s内的误差，可以看出四种算法中，NTSTO-SEM 算法和Newmark-SEM 算法表现比RK3-SEM 算法和HHT-α-SEM算法稍好，但四种算法差别很小，说明短时程计算时各类算法都能够有好的逼近.
|
图 3 四种算法与精确解在0.6 s内的误差比较 Fig. 3 The residual of elastic waveform comparisons of the four methods within time of 0.6 s |
为了测试以上四种算法对于长时程运算的有效性，我们考虑上述模型的四周边界都用自由边界条件，进行多步计算后观察波场的频散情况[13, 35-36].
图 4-7给出了四种算法在第16s(16000个时间步)的波场快照，从快照可以明显看出在长时程跟踪能力方面，NTSTO-SEM 算法优于Newmark-SEM、非辛的RK3-SEM、HHT-α-SEM 三种算法，说明NTSTO-SEM 算法对误差的积累能够很好地抑制；在上述模型几何参数下，我们又选用算法考量的通用参数指标:内存消耗、运算速度和稳定性几个方面分别对四种算法进行对比，对比结果见表 1.
|
图 4 均匀介质中用HHT-α-SEM算法求解的第16 s (16000时间步）的波场快照 Fig. 4 Snapshots of elastic wavefields in the homogeneous model at time 16 s (16000 time steps) generated by HHT-α-SEM |
|
图 5 均匀介质中用Newmark-SEM算法求解的第16 s (16000时间步）的波场快照 Fig. 5 Snapshots of elastic wavefields n the homogeneous model at time 16 s (16000 time steps) generated by Newmark-SEM |
|
图 6 均匀介质中用NTSTO-SEM算法求解的第16 s (16000时间步）的波场快照 Fig. 6 Snapshots of elastic wavefields in the homogeneous model at time 16 s (16000 time steps) |
|
图 7 均匀介质中用非辛RK3-SEM算法求解的第16 s (16000时间步）的波场快照 Fig. 7 Snapshots of elastic wavefields in the homogeneous model at time 16 s (16000 time steps) generated by non-symplectic RK3-SEM |
|
表 1 四种算法的计算指标对比 Table 1 The comparisons of computation index of the four methods |
从表中可以看出，内存的消耗上，由于HHT-α-SEM 及RK3-SEM 算法都属于预估-校正方法，所以占用内存较多，而NTSTO-SEM 与显式Newmark-SEM 算法只利用上一迭代步信息进行迭代，因此是四种算法中最少消耗内存的两种算法.在本试验参数条件下，NTSTO-SEM 算法占用内存量是HHT-α-SEM 算法的78.3%，是RK3-SEM算法的87.7%；从稳定的条件下可以选取的最大步长来看，NTSTO-SEM 算法的稳定性远优于其它三种算法，也表明了辛算法可用大步长数值积分的优势；如果以每种算法允许的最大稳定步长模拟1s时长的耗时来看，NTSTO-SEM 算法的耗时多于Newmark-SEM 和HHT-α-SEM 算法的，但仅是速度最快的Newmark-SEM 算法耗时的1.32倍，而跟踪能力远胜于Newmark-SEM、HHT-α-SEM算法；因此，综合内存消耗、稳定性、运算速度和长时程跟踪能力各考量指标来看，NTSTO-SEM 算法不乏为优越的算法，进行地震波场模拟时既能保结构，又能达到满意的模拟效果.
3.2 试验2为了进一步体现算法对于复杂几何模型的有效性，我们考虑如图 8的起伏地表的三层介质模型，模型大小横向为6000m，纵向平均大小为6300m，地表有较大的起伏，从上向下的三层介质分别记作Ⅰ、Ⅱ、Ⅲ，物性参数如表 2所示，其中介质Ⅰ和Ⅱ之间的界面为一起伏界面，介质Ⅱ和Ⅲ之间的界面为直线界面，四边形网格剖分如图 8所示，横向剖分的单元数为127个，介质Ⅰ的纵向单元数为47 个，介质Ⅱ的纵向单元数为35个，介质Ⅲ的纵向单元数45个.震源是位于(3000 m，3000 m)处的集中力源，由主频为f0=20Hz的Ricker子波f(t)= [1-2(πf0(t-t0))2]exp(- (πf0(t-t0))2)激发产生，其中t0=0.05s；时间步长为0.3ms，每个单元边上都用7个GLL 插值节点.模型的左、右以及下边界都用吸收边界条件.
|
图 8 起伏地表的多层介质模型 Fig. 8 Multilayered model with irregular topography |
|
|
表 2 各层介质物性参数 Table 2 The model parameters of the three layers |
为了说明NTSTO-SEM 算法对复杂模型的有效性，我们与Newmark-SEM 算法的特殊情形---Leapfrog(中心差分)结合谱元法(Leapfrog-SEM)(见附录B)进行比较.图 9 给出了Leapfrog-SEM算法和NTSTO-SEM 算法求解起伏地表多层介质模型所得的波形记录的对比情况，总时长1.5s，检波点设在(2160m，4500m)处；图 10给出了该模型用两种算法求解的水平分量在1.2s时刻的波场快照.图 9和图 10 的比较结果显示两种算法的波形和快照几乎都没有差别.同时从快照可以看出，波前对层间和对起伏表面的响应都非常清晰地得到了模拟，说明NTSTO-SEM 算法对于复杂几何形状和介质的模拟都非常有效.
|
图 9 Leapfrog-SEM算法和NTSTO-SEM算法的波形对比 Fig. 9 Waveform comparisons generated by the Leapfrog-SEM algorithm and NTSTO-SEM algorithm |
|
图 10 分别由Leapfrog -SEM算法(a)和NTSTOSEM算法(b)求解的起伏地表的多层介质模型在1. 2 s时刻水平分量的波场快照 Fig. 10 Snapshots of horizontal component for multilayered model with irregular topography at time 1.2 s generated by the Leapfrog-SEM algorithm (a) and the NTSTO-SEM algorithm (b) |
为了更进一步观察面波在地表起伏时的传播情况，将震源置于图 8 所示模型中的地表(3120 m，6050m)处，选择震源为主频f0=20 Hz的Ricker子波的爆炸源，时间间隔为0.2 ms，每个单元上的GLL插值节点选取为10个，在地表设置一个检波器排列，横向平均间距15 m，其它参数选取与试验2相同.
图 11、12分别给出了0.8s和1.2s时刻模型上半部分的波场快照；图 13给出了地表单个检波器(2950m，6220m)处接收的水平分量和垂直分量合成地震记录；图 14给出了地表排列检波器接收的水平位移和垂直位移分量的单炮地面记录.从图 11-14中可以明显地观察到沿自由表面及近自由表面区域传播的比较强的Rayleigh面波(R)，它跟在横波的后面，能量比较集中，而随着传播时间的增加，与纵横波逐渐分离，非常清晰，且Rayleigh 面波沿界面传播，随深度增加迅速衰减，这与Rayleigh 波属性完全吻合.注意到本试验采用平均大小为47m×50m 的四边形单元，最短波长中包含的单元数小于1.6 个，而模拟结果清晰可辨，说明NTSTO-SEM 算法精度较高，而且其对表面波的捕捉也非常清晰，这也说明NTSTO-SEM 算法对模拟起伏地表多层复杂几何模型是非常有效的.
|
图 11 NTSTO-SEM算法求解的爆炸源在地表(3120 m，6050 m)处产生的水平分量在0. 8 s时刻的部分波场快照 Fig. 11 Snapshots of horizontal component for multilayered model with irregular topography at time 0.8 s generated by the NTSTO-SEM algorithm and the explosion source is located at free surface point (3120 m，6050 m) |
|
图 12 TSTO-SEM算法求解的爆炸源在地表(3120 m，6050 m)处产生的水平分量在1. 2 s时刻的部分波场快照 Fig. 12 Snapshots of horizontal component for multilayered model with irregular topography at time 1.2 s generated by the NTSTO-SEM algorithm and the explosion source is located at free surface point (3120 m, 6050 m) |
|
图 13 地表接收点(950 m，6220 m)处水平分量和垂直分量(仄）的波形记录 Fig. 13 Waveforms of horizontal (upper) and vertical (below) displacement components at the receiver located at free surface point (2950 m, 6220 m) |
|
图 14 自由表面水平位移分量(a)和垂直位移分量(b)合成地面记录 Fig. 14 Synthetic Seismograms of horizontal (a) and vertical (b) displacement components along free surface computed by the NTSTO-SEM method |
本文构造了一种基于新的三阶辛和谱元法结合的多辛结构的时间-空间保结构算法，通过与解析解对比，说明该算法对地震波短时程模拟是非常有效的；并与目前广泛采用的Newmark-SEM 等算法的比较，得出了该算法具有内存消耗少、稳定性好、长时程跟踪能力强和运算速度较快的特点；最后通谱元法对于复杂几何形状和多层介质模型的有效性，对表面波的传播也得到清晰的模拟.各类试验表明，在地震波场模拟中，NTSTO-SEM 算法既能对复杂几何模型中地震波的传播特性进行高精度模拟，又能为地震波的长时程计算及模拟(如地球自由振荡数值模拟)研究提供了新的较可靠的研究工具.
附录A辛算法是求解Hamilton 系统的一种数值方法.为了深入了解这种数值方法，我们先从辛几何的几个概念讲起.
设向量x，y∈ R2n，定义辛内积(x，y)=
该内积具有性质
(1) 双线性:(x+y，v+w)= (x，v)+(x，w)+(y，v)+ (y，w).
(2) 反对称性:(x，y)=- (y，x).
(3) 非退化性:∀y≠0，∃x，(x，y)=0⇒x=0.
定义1 V为R2n上的向量空间，ω 为V×V上的双线性映射，如果该映射满足反对称性(2)和非退化性(3)，则(V，ω)为辛空间，ω 为辛结构.
定义2 ∀x，y∈V，S为(V，ω)上的线性变换，如果满足(Sx，Sy)= (x，y)，则称S为辛变换.
定义3 矩阵S∈ R2n×2n，如果STJ S=J ，或SJ ST =J ，则S称为辛矩阵，其中J 为单位辛矩阵.由所有辛矩阵构成的群，称为辛群，记为Sp(2n).
Hamilton力学的基本定理:任何一个Hamilton系统的解是一个单参数辛群.
求解Hamilton方程保持辛结构的离散数值算法称为辛算法.任何一个差分格式不论是显式或隐式的，它都能看成从这一时刻到下一时刻的映射，如果这个映射是辛的，称该差分格式是辛几何格式，简称辛格式.如，对于单步差分格式vk=Skvk-1 ，且Sk为辛矩阵，即由vk-1 到vk的变换为辛变换，则此格式为辛格式.
附录BNewmark 格式，由Newmark在文献[36]中引入.考虑一个标准的连续Lagrangian系统，已证明经典的Newmark格式实际上是变分的[53]，
|
用隐式关系由(qk，$\dot{q}$k)更新到(qk+1，$\dot{q}$k+1)，即
|
式中M为质量矩阵，V为势函数.
对于方程MÜ =F-KU，记G(t)=F-KU，则方程可写为:MÜ =G(t)，在时间[0，T]区间内，将区间分为[ts，ts+1]，已有t0 =0，tS=T，ts+1 =ts+Δt，s=0，1，…，(S-1).引入变量β，γ 为实数，则方程的Newmark[38]离散格式为
|
当β ≠
|
其稳定的Courant数条件为:nC = max[cΔt/Δx].当γ=0，β=
|
对于方程MÜ =F-KU的Leapfrog格式为
所以，Us+1 = (Δt)2[M-1(F-KUs)]+2Us-Us-1 .
附录CHHT-α 算法[39]是Hilber，Hughes 和Taylor构造的α 方法，通过引入α，β，γ 三个控制变量来实现.对于二阶微分方程y″ =f(t，y，y′)可以写为:y′ =z，z′ =a，a=f(t，y，z)，则算法为
|
方程MÜ +KU=F经过离散可写为:
|
式中
|
这里α ∈ [-1/3，0]，β = (1-α)2/4，γ = (1-2α)/2，则微分方程可改写为
|
式中已记
|
| [1] | Alterman Z, Karal F C Jr. Propagation of elastic waves in layered media by finite difference methods. Bull. Seisml. Soc. Am. , 1968, 58(1): 367-398. |
| [2] | Boore D M. Finite-difference methods for seismic wave propagation in heterogeneous materials. //Bolt B A ed. Methods in Computational Physics 11, New York: Academic Press. Inc., 1972. |
| [3] | Madariaga R. Dynamics of an expanding circular fault. Bull. Seismol. Soc. Am. , 1976, 66(3): 639-669. |
| [4] | Virieus J. P-SV wave propagation in heterogeneous media: velocity-stress finite-difference method. Geophysics , 1986, 51(4): 889-901. DOI:10.1190/1.1442147 |
| [5] | Mora P. Modeling anisotropic seismic waves in 3-D. 59th Ann. Internat. Mtg. Soc. Expl. Geophys., Expanded Abstracts, 1989, 2: 1039-1043. |
| [6] | Frankel A, Vidale J. A three-dimensional simulation of seismic waves in the Santa Clara valley, California, from the Loma Prieta aftershock. Bull. Seismol. Soc. Am. , 1992, 82(5): 2045-2074. |
| [7] | Frankel A. Three-dimensional simulations of ground motions in the San Bernardino valley, California, for hypothetical earthquakes on the San Andreas fault. Bull. Seismol. Soc. Am. , 1993, 83(4): 1020-1041. |
| [8] | Pitarka A. 3D elastic finite-difference modeling of seismic motion using staggered grids with nonuniform spacing. Bull. Seisml. Soc. Am. , 1999, 89(1): 54-68. |
| [9] | Olsen K B. Site amplification in the Los Angeles basin from three-dimensional modeling of ground motion. Bull. Seismol. Soc. Am. , 2000, 90(6B). |
| [10] | Zingg D W. Comparison of high-accuracy finite-difference methods for linear wave propagation. SIAM Journal on Scientific Computing , 2000, 22(2): 476-502. DOI:10.1137/S1064827599350320 |
| [11] | Chaljub E, Komatitsch D, Vilotte J P, et al. Spectral-element analysis in seismology. Advances in Geophysics , 2007, 48: 365-419. DOI:10.1016/S0065-2687(06)48007-9 |
| [12] | Marfurt K J. Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations. Geophysics , 1984, 49(5): 533-549. DOI:10.1190/1.1441689 |
| [13] | Komatitsch D, Vilotte J P. The spectral element method: An efficient tool to simulate the seismic response of 2D and 3D geological structures. Bull Seisml. Soc. Am. , 1998, 88(2): 368-392. |
| [14] | Tessmer E, Kosloff D. 3-D elastic modeling with surface topography by a Chebyshev spectral method. Geophysics , 1994, 59(3): 464-473. DOI:10.1190/1.1443608 |
| [15] | Carcione J M. The wave equation in generalized coordinates. Geophysics , 1994, 59(12): 1911-1919. DOI:10.1190/1.1443578 |
| [16] | Furumura M, Kennett B L N, Furumura T. Seismic wavefield calculation for laterally heterogeneous Earth models-II: the influence of upper mantle heterogeneity. Geophys.J. Int. , 1999, 139(3): 623-644. DOI:10.1046/j.1365-246x.1999.00962.x |
| [17] | Furumura T, Kennett B L N, Furumura M. Seismic wavefield calculation for laterally heterogeneous whole Earth models using the pseudospectral method. Geophys. J. Int , 1998, 135(3): 845-860. DOI:10.1046/j.1365-246X.1998.00682.x |
| [18] | Carcione J M, Wang P J. A Chebyshev collocation method for the wave equation in generalized coordinates. Comp. Fluid Dyn. J. , 1993, 2: 269-290. |
| [19] | Shen J, Tang T. Spectral and High-Order Methods with Applications. Beijing: Science Press, 2006 . |
| [20] | Komatitsch D, Martin R, Tromp J, et al. Wave propagation in 2-D elastic media using a spectral element method with triangles and quadrangles. J. Comput. Acoust. , 2001, 9(2): 703-718. DOI:10.1142/S0218396X01000796 |
| [21] | Chaljub E, Valette B. Spectral element modelling of three-dimensional wave propagation in a self-gravitating Earth with an arbitrarily stratified outer core. Geophys.J. Int. , 2004, 158(1): 131-141. DOI:10.1111/gji.2004.158.issue-1 |
| [22] | Chaljub E, Capdeville Y, Vilotte J P. Solving elastodynamics in a fluid-solid heterogeneous sphere: a parallel spectral element approximation on non-conforming grids. J. Comput. Phys. , 2003, 187(2): 457-491. DOI:10.1016/S0021-9991(03)00119-0 |
| [23] | Komatitsch D, Ritsema J, Tromp J. The spectral-element method, Beowulf computing, and global seismology. Science , 2002, 298(5599): 1737-1742. DOI:10.1126/science.1076024 |
| [24] | Komatitsch D, Barnes C, Tromp J. Simulation of anisotropic wave propagation based upon a spectral element method. Geophysics , 2000, 65(4): 1251-1260. DOI:10.1190/1.1444816 |
| [25] | Komatitsch D, Tromp J. Spectral-element simulations of global seismic wave propagation-I. Validation. Geophys.J. Int. , 2002, 149(2): 390-412. DOI:10.1046/j.1365-246X.2002.01653.x |
| [26] | Komatitsch D, Tromp J. Spectral-element simulations of global seismic wave propagation-II. Three-dimensional models, oceans, rotation, and self-gravitation. Geophys.J. Int. , 2002, 150(1): 303-318. DOI:10.1046/j.1365-246X.2002.01716.x |
| [27] | Komatitsch D, Tromp J. Introduction to the spectral-element method for three-dimensional seismic wave propagation. Geophys.J. Int. , 1999, 139(3): 806-822. DOI:10.1046/j.1365-246x.1999.00967.x |
| [28] | Komatitsch D, Barnes C, Tromp J. Wave propagation near a fluid-solid interface: a spectral-element approach. Geophysics , 2000, 65(2): 623-631. DOI:10.1190/1.1444758 |
| [29] | Liu H P, Anderson D L, Kanamori H. Velocity dispersion due to anelasticity: implications for seismology and mantle composition. Geophys.J. R. Astron. Soc. , 1976, 47(1): 41-58. DOI:10.1111/j.1365-246X.1976.tb01261.x |
| [30] | Priolo E, Carcione J M, Seriani G. Numerical simulation of interface waves by high-order spectral modeling techniques. J. Acoust. Soc. Am. , 1994, 95(2): 681-693. DOI:10.1121/1.408428 |
| [31] | 冯康, 秦孟兆. 哈密尔顿系统的辛几何算法. 杭州: 浙江科技出版社, 2003 . Feng K, Qin M Z. Symplectic Geometric Algorithm for Hamiltonian Systems (in Chinese). Hangzhou: Zhejiang Science and Technology Press, 2003 . |
| [32] | 钟万勰. 分析结构力学与有限元. 动力学与控制学报 , 2004, 2(4): 1–8. Zhong W X. Analytical structural mechanics and finite element. Journal of Dynamics and Control (in Chinese) , 2004, 2(4): 1-8. |
| [33] | 王雨顺, 王斌, 秦孟兆. 偏微分方程的局部保结构算法. 中国科学A辑: 数学 , 2008, 38(4): 377–397. Wang Y S, Wang B, Qin M Z. Local structure-preserving algorithms for partial differential equations. Science in China (Series A: Mathematics) (in Chinese) , 2008, 38(4): 377-397. |
| [34] | Iwatsu R. Two new solutions to the third-order symplectic integration method. Phys. Lett. A. , 2009, 373(34): 3056-3060. DOI:10.1016/j.physleta.2009.06.048 |
| [35] | Li X F, Li Y Q, Zhang M G, et al. Scalar seismic-wave equation modeling by a multisymplectic discrete singular convolution differentiator method. Bulletin of the Seismological Society of America , 2011, 101(4): 1701-1718. |
| [36] | Li X F, Wang W S, Lu M W, et al. Structure-preserving modelling of elastic waves: a symplectic discrete singular convolution differentiator method. Geophys.J. Int. , 2012, 188(3): 1382-1392. DOI:10.1111/gji.2012.188.issue-3 |
| [37] | Newmark N M. A method of computation for structural dynamics. Journal of the Engineering Mechanics Division, ASCE , 1959, 85(3): 67-94. |
| [38] | Kane C, Marsden J E, Ortiz M, et al. Variational integrators and the Newmark algorithm for conservative and dissipative mechanical systems. International Journal for Numerical Methods in Engineering , 2000, 49(10): 1295-1325. DOI:10.1002/(ISSN)1097-0207 |
| [39] | Negrut D, Rampalli R, Ottarsson G, et al. On the use of the HHT method in the context of index 3 differential algebraic equations of multibody dynamics. ASME Conf. Proc. , 2005, 6: 207-218. |
| [40] | Jay L O, Negrut D. Extensions of the HHT- α method to differential-algebraic equations in mechanics. Electronic Transactions on Numerical Analysis , 2007, 26: 190-208. |
| [41] | Cohen G, Fauqueux S. Mixed finite elements with mass-lumping for the transient wave equation. J. Comput. Acoust. , 2000, 8(1): 171-188. DOI:10.1142/S0218396X0000011X |
| [42] | Mercerat E D, Vilotte J P, Sánchez-Sesma F J. Triangular spectral element simulation of two-dimensional elastic wave propagation using unstructured triangular grids. Geophys.J. Int. , 2006, 166(2): 679-698. DOI:10.1111/gji.2006.166.issue-2 |
| [43] | Giraldo F X, Warburton T. A nodal triangle-based spectral element method for the shallow water equations on the sphere. J. Comput. Phys. , 2005, 207(1): 129-150. DOI:10.1016/j.jcp.2005.01.004 |
| [44] | Bittencourt M L. Fully tensorial nodal and modal shape functions for triangles and tetrahedra. Int. J. Numer. Meth. Engng. , 2005, 63(11): 1530-1558. DOI:10.1002/(ISSN)1097-0207 |
| [45] | Taylor M A, Wingate B A, Vincent R E. An algorithm for computing Fekete points in the triangle. SIAM Journal on Numerical Analysis , 2000, 38(5): 1707-1720. DOI:10.1137/S0036142998337247 |
| [46] | 邢誉峰, 冯伟. 李级数算法和显式辛算法的相位分析. 计算力学学报 , 2009, 26(2): 167–171. Xing Y F, Feng W. Phase analysis of Lie Series algorithm and explicit symplectic algorithm. Chinese Journal of Computational Mechanics (in Chinese) , 2009, 26(2): 167-171. |
| [47] | Hairer E, Lubich C, Wanner G. Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations. Berlin: Springer-Verlag, 2006 . |
| [48] | Van Der Houwen P J, Sommeijer B P. Runge-Kutta (-Nystrǒm) methods with reduced phase errors for computing oscillating solutions. SIAM Journal on Numerical Analysis , 1987, 24(3): 595-617. DOI:10.1137/0724041 |
| [49] | Ruth R D. A canonical integration technique. IEEE Transactions on Nuclear Science , 1983, 30(4): 2669-2671. DOI:10.1109/TNS.1983.4332919 |
| [50] | Suzuki M. General theory of higher-order decomposition of exponential operators and symplectic integrators. Physics Letters A , 1992, 165(5-6): 387-395. DOI:10.1016/0375-9601(92)90335-J |
| [51] | McLachlan R I, Atela P. The accuracy of symplectic integrators. Nonlinearity , 1992, 5(2): 541-562. DOI:10.1088/0951-7715/5/2/011 |
| [52] | Festa G, Vilotte J P. The Newmark scheme as velocity-stress time-staggering: an efficient PML implementation for spectral element simulations of elastodynamics. Geophys. J. Int. , 2005, 161(3): 789-812. DOI:10.1111/gji.2005.161.issue-3 |
| [53] | Zampieri E, Pavarino L F. Approximation of acoustic waves by explicit Newmark's schemesand spectral element methods. Journal of Computational and Applied Mathematics , 2006, 185(2): 308-325. DOI:10.1016/j.cam.2005.03.013 |
| [54] | Stacey R. Improved transparent boundary formulations for the elastic-wave equation. Bull. Seisml. Soc. Am. , 1988, 78(6): 2089-2097. |
| [55] | De Hoop A. A modification of Cagniard's method for solving seismic pulse problems. Appl. Sci. Res. , 1960, 8(1): 349-356. DOI:10.1007/BF02920068 |