基于辛格式离散奇异褶积微分算子的弹性波场模拟

李一琼，李小凡，张美根

1 中国科学院地质与地球物理研究所，北京 100029

2 中国地震局地球物理研究所，北京 100081

基于辛格式离散奇异褶积微分算子的弹性波场模拟

李一琼1，2，李小凡1，张美根1

1 中国科学院地质与地球物理研究所，北京 100029

2 中国地震局地球物理研究所，北京 100081

本文发展了基于辛格式离散奇异褶积微分算子（SDSCD）的保结构方法模拟弹性波场，求解弹性波动方程时，引入辛差分格式进行时间离散，采用离散奇异褶积微分算子进行空间离散.相比于传统的伪谱方法，该方法提高了计算精度和稳定性.数值结果表明SDSCD方法可以有效地抑制数值频散，为解决大尺度、长时程地震波场模拟问题提供了合适的数值方法.

辛算法，离散奇异褶积微分算子，弹性波场模拟

1 引 言

研究地震波在地球介质中的传播规律，是地球物理学中认识地球内部结构的一项重要手段.在地震波场数值模拟方法中，波动方程数值模拟方法一直占有重要的地位，波动方程模拟包含了地震波丰富的动力学信息，为研究地震波的传播机理提供了更多的依据.描述地球内部的地震波场的波动方程多采用声波方程和弹性波方程.声波方程应用于地震波研究，虽已获得诸多成果，但声波模型多适用于勘探地震学.因地球介质一般被视为非均匀固体介质，在理论地震学研究中通常须考虑其弹性性质及地震波的矢量性质，简单、机械地套用声波近似，其精度及物理合法性均将受到质疑.弹性波近似较为繁复，更贴近实际、更广泛地应用于地震波研究.本质上，声波近似可视为弹性波近似的一个特例和简化.因此，多采用弹性波动方程描述地震波.

地震波波动方程的定解问题包括微分算子、震源、初始条件和边界条件等.常见的正演数值模拟方法依据处理微分项的数值方式而不同，如有限差分法［1］、伪谱法［2］和有限元法［3］等都已被运用多年，较为成熟.有限差分法是基于Taylor展开，用差分代替微分，将微分方程转化为代数方程，在计算上易于实现，但有些差分格式是根据微分算子凭经验凑算得到的，缺乏数学理论保证［4］.伪谱法采用了快速傅里叶变换计算空间微分，使其具有较快的计算速度，计算精度要高于经典有限差分法.但计算空间微分时动用了区域内的所有变量，有悖于微分的局部属性，对于复杂介质的地震波模拟结果并不理想［5－6］.以上两种数值方法都不可避免地存在差分近似和网格离散化近似，以及无穷级数的有限截断近似，将破坏波动方程及解的结构，造成数值解相对于精确解的偏离.这一数值方法的固有问题，在一定范围内，显现得不是很突出，从而被认为是有效的数值方法，其计算结果被接纳.与前两种算法相比，有限元法有较强的数学理论保障，其理论基础是变分原理，离散化构建在Sobolev函数空间框架下［7］，通过将单位元上的分片插值函数和最小位能理论相结合求解整个区域的波场信息，但其计算过程复杂，计算量大，应用到复杂介质中还有一定的困难.

在研究大尺度的地震波场模拟时，我们关注的是较长时程的地震波计算问题，其波动方程在时间域处理时，具有很强的时空继承性，其积累误差较为明显，需要选择具有高精度的计算方法.辛算法这一保结构算法正是在诸如这样的需求下提出的，算法构建在哈密尔顿体系下，其构造之初，就要求在方程进行离散化和有限化近似过程中，尽可能保持运动方程的结构及其相应的守恒量不变，或者，把数值解对于精确解的偏离控制在设定的精度以内，使得辛算法能从根本上提高解的精度和可信度.计算实验显示，辛算法优异的稳定性和长时间跟踪能力有着重要的应用前景.

辛算法在众多科学计算领域得到了应用和发展.地震波场在高频近似条件下，其射线解刻画了地震波场的射线走向和波阵面形态，这在地震波场传播的渐近解研究、模型研究及波场成像的运用中具有广泛的实用价值.然而，算法本身的一个重大缺陷是在焦散区失效［8］，1965年前苏联学者Maslov将射线问题转化为辛空间中的Lagrange子流形问题，成功地解决了焦散问题，从而大大拓广了几何光学应用范围.陈景波［9］和高亮［7］等用辛几何算法解决地震射线追踪问题得到了较好的结果，和四阶Runge－Kutta方法相比，精度相当，但辛算法的计算速度要快.在基于波动方程的地震波数值模拟方面，Chen［10］将辛格式和伪谱法相结合求解了标量波波动方程，和非辛格式的算法相比，辛格式显示了较小的误差增长.运用辛算法能很好地解决时间域波动方程的时间继承性，因此，我们将为波动方程中时间微分项引入辛差分格式的保结构算子.

波动方程中空间微分项的近似计算同样至关重要.20世纪70年代早期，就有学者提出褶积微分算子的概念，但将其应用到地震波场模拟中比较晚.Holberg［11］运用Fourier变换，对褶积微分算子进行了设计，并在三维弹性波数值模拟中取得了较好的结果.Zhou等［12］、张中杰等［13］和戴志阳等［14］也对这一方法进行了深入研究和发展.龙桂华等［15］在前人的工作基础上，根据广义函数理论，引入奇异核褶积微分算子，经过加窗处理的微分算子与波场变量进行褶积，来求取波场变量关于空间微分.褶积微分算子和伪谱法相比，其优势在于褶积微分算子可以优化为短算子，更能突出空间微分本身的局部属性；而伪谱法在计算微分时须利用全局信息.褶积微分算子和有限差分算子相比，精度更高，能够弥补有限差分算子描述复杂介质的不足.因此褶积微分算子应用于波动方程的空间微分近似是非常有效的.

本文对弹性波动方程引入显辛差分格式进行时间离散，采用离散奇异褶积微分算子进行空间离散，建立一个新的弹性波场模拟方法，辛格式离散奇异褶积微分算子方法（SDSCD）.该方法尽管会增加一些计算量，但将提高计算精度和稳定性.与传统的伪谱方法相比，SDSCD方法是全局保结构的，具有高精度和较强的长时间跟踪能力，将为解决大尺度、长时程地震波场模拟问题提供有效的数值方法.

2 弹性波动方程的辛差分格式

从弹性动力学的本构方程出发，可以推导出各向同性介质中一阶速度应力方程：其中，vx和vz分别为沿x和z方向的速度场，σxx、σzz和σxz为应力分量，λ和μ为表征介质属性的弹性参数，ρ为介质密度.

令v＝（vx，vz），σ＝（σxx，σzz，σxz，σzx），则（v，σ）＝（vx，vz，σxx，σzz，σxz，σzx），得到方程（1）的矩阵形式：

注意到，矩阵P，Q中含有关于空间变量的微分算子，处理这些空间微分项有赖于第3节的奇异核褶积微分算子，暂不详述.

由我们的先前工作［16－17］，可得到（2）式的k级n阶显式辛格式

其中，波场变量的上标（t）表示时间序列中的当前时间点，（t＋1）为下一时间点，Δt表示时间步长，矩阵I为与P，Q同阶的单位矩阵.

当k＝3时，（3）式写为弹性波动方程（2）的单步三级三阶显辛格式解［18］

式（4）中的空间导数可采用第3节的奇异核褶积微分算子计算.

3 离散奇异褶积微分算子

根据广义函数理论，当已知某函数的离散序列｛u（xk）｝，该函数的微分形式可由奇异核函数的微分来近似［19］

其中，δ（x）是奇异核Delta函数，（q）表示关于x的q阶广义导数，W为算子半带宽.W和高斯包络σ是可变的，适当选取可提高微分的精度.

奇异核Delta函数常被用于偏微分方程的数值求解方面，在数值计算中，需要考虑用近似函数来逼近物理上很难实现的理论Delta函数.这里选用了归一化的Shannon奇异核函数

其中，Δ为空间步长，σ为高斯函数包络线宽度.对于任意给定的σ≠0，当Δ→0时，（6）式趋近于Delta函数.

当q＝1时，一阶微分算子表示为

这个褶积微分算子系数关于k＝0是中心对称的.

在实际计算过程中，对上述算子进行了加窗处理，从而有效地抑制由于算子截断而引起的Gibbs效应.选用Hanning窗函数

得到处理后的微分算子：

在这里，选取Hanning窗函数的参数α＝0.54和β＝8.

4 算法的稳定性条件

在时间域地震波场数值模拟中，计算的时间步长与空间步长的离散尺度必须满足一定的条件，才能保证迭代计算的可行.否则，可能产生随指数量级增长的计算误差，造成严重的数值频散，最终导致数值溢出使计算无法进行.对一种数值解法，需要知道计算稳定的离散参数区间，即分析解法的稳定性.

在第2节得到波场变量关于时间迭代的辛差分显格式（3）式，现仅考虑三级三阶的格式.令矩阵Mi＝（I＋ciΔtP）（I＋diΔt Q），i＝1，2，3，矩阵M＝M3M2M1，则（3）式写为

形如（10）式的差分格式算法稳定的充要条件是矩阵M的谱半径ρ（M）≤1.由于矩阵P，Q中含有关于空间变量的微分，不能直接考察矩阵M，需要采取Fourier分析法［20－21］.

利用第3节的褶积微分算子，可得函数u（x）微分的离散近似表达式

将（11）式两边进行Fourier变换，并由算子系数的中心对称，得到：

对（10）式进行关于空间变量的Fourier变换，P，Q中关于空间变量的微分项由（11）式的算子代替，变换后的矩阵M仍记为M.

考虑如下定理，设A∈Cn×n，则对Cn×n上的任一矩阵范数‖·‖，有ρ（A）≤‖A‖.

由此算法稳定的充要条件等价为矩阵M的任一范数‖M‖≤1.

在实际推导过程中发现，矩阵M不具有特殊性质，不能给出‖M‖一个比较简单的表达式，故将数值实验选取的参数代入矩阵M，采取试算的方式，来保证满足稳定性条件.在数值实验中，对于均匀介质模型，采用SDSCD方法，最大Courant数为0.684，而在同一模型中，常规Fourier伪谱算法的最大Courant数为0.64.

5 数值实验

为了验证方法的有效性，将本文提出的SDSCD算法结果与非辛算法的结果进行比较.这里所说的非辛算法是指空间离散采用Fourier伪谱算子近似，时间离散用向前差分近似的常规Fourier伪谱算法.下文中除特别说明外，非辛算法均指此伪谱算法.

在验证算例中，首先采用了一个具有粗糙倾斜界面的分层介质模型.模型的物性参数见图1，五角星★表示震源所在位置（1280，1080），三角形▲为检波点所在位置（1180，980）；震源为胀缩源，由主频为30Hz的Ricker子波激发产生；时间步长为1ms，模型沿横向和纵向各剖分256个网格点，空间间隔均为10m.

图2和图3分别给出了SDSCD算法和伪谱算法在0.3s时刻的波场快照.这些快照都非常清晰地刻画了该模型中的直达P波、反射P波、反射P－S波、透射P波、透射P－S波以及散射波.可以看到模型上下两层的波速对比度较大，散射波信息十分丰富.图4和图5给出了检波点处0.3s内两种算法分别在水平波场和垂直波场的地震记录.波场快照和地震记录都说明了SDSCD算法的有效性.

图1 具有粗糙倾斜界面的分层均匀介质模型Fig.1 Homogeneous layered model with rough inclined interface

为了检验SDSCD算法的保结构特性，我们对一个均匀介质模型的波场情况进行了长时程的跟踪，模型边缘没有加任何边界条件.此做法相当于增大了模型的尺度，通过进行较多次时间迭代步，考察算法是否保持稳定.Chen［22］采用同样的做法对标量波模型的辛Nyström算法进行了验证.

模型的物性参数为：纵波波速4000m／s，横波波速2309.3m／s；密度2400kg／m3，模型沿横向和纵向各剖分出256个网格点，空间间隔均为20m；震源位于模型中央（2560m，2560m）处，震源为集中力源，由主频为20Hz的Ricker子波激发产生；时间步长为2ms.图6和图7分别给出了SDSCD算法和伪谱算法在300次时间迭代步后，即0.6s时刻的水平波场快照，可以看到，伪谱算法的结果数值频散较为明显.图8和图9分别给出了SDSCD算法和伪谱算法在5000次时间迭代步后，即10s时刻的水平波场快照，SDSCD算法模拟的波前面十分清晰，而伪谱算法的结果，则很难辨清波前的形状.由此可知，SDSCD方法具有较明显的保结构性，能够高精度处理长时程波场模拟.

6 结论与讨论

本文给出SDSCD算法用于弹性波场模拟，该方法在计算过程中，存储较少，易于实现，虽然相比本文用于对比的经典伪谱方法，增加了一些计算迭代步骤和时间，但在抑制数值频散和提高精度方面都有明显改善，并且具有较强的稳定性.可以在一定范围内加大时间步长和空间间隔，作为研究大尺度、长时程的地震波场模拟的有效方法.

作为方法研究，本文还有一些工作尚待完成，如SDSCD算法的稳定性条件尚未给出解析表达式，仅采取了试算，以保证数值实验的可行.这种方式对于客观评价本算法的优劣略显不便.这需要进一步的理论研究，使之更为完善.

（References）

［1］ Dablain M A.The application of high－order differencing to the scalar wave equation.Geophysics，1986，51（1）：54－66.

［2］ Fornberg B.The pseudospectral method：Comparisons with finite differences for the elastic wave equation.Geophysics，1987，52（4）：483－501.

［3］ Smith W D.The application of finite element analysis to body wave propagation problems.Geophys.J.Int.，1975，42（2）：747－768.

［4］ 钟万勰.应用力学的辛数学方法.北京：高等教育出版社，2006.Zhong W X.Symplectic Solution Methodology in Applied Mechanics（in Chinese）.Beijing：Higher Education Press，2006.

［5］ Özdenvar T，McMechan G A.Causes and reduction of numerical artifacts in pseudo－spectral wavefield extrapolation.Geophys.J.Int.，1996，126（3）：819－828.

［6］ Li X F，Zhu T，Zhang M G，et al.Seismic scalar wave equation with variable coefficients modeling by a new convolutional differentiator.Computer Physics Communications，2010，181（11）：1850－1858.

［7］ 高亮，李幼铭，陈旭荣等.地震射线辛几何算法初探.地球物理学报，2000，43（3）：402－410.Gao L，Li Y M，Chen X R，et al.An attempt to seismic ray tracing with symplectic algorithm.Chinese J.Geophys.（in Chinese），2000，43（3）：402－410.

［8］ 李世雄.波动方程的高频近似与辛几何.北京：科学出版社，2001.Li S X.The High－Frequency Approximation of Wave Equations and Symplectic Algorithm（in Chinese）.Beijing：Science Press，2001.

［9］ 陈景波，秦孟兆.辛几何算法在射线追踪中的应用.数值计算与计算机应用，2000，21（4）：255－265.Chen J B，Qin M Z.Ray tracing by symplectic algorithm.Numerical Computation and Computer Application（in Chinese），2000，21（4）：255－265.

［10］ Chen J B.Lax－Wendroff and Nyström methods for seismic modelling.Geophysical Prospecting，2009，57（6）：931－941.

［11］ Holberg O.Computational aspects of the choice of operator and sampling interval for numerical differentiation in largescale simulation of wave phenomena.Geophysical Prospecting，1987，35（6）：629－655.

［12］ Zhou B，Greenhalgh S A.Seismic scalar wave equation modeling by a convolutional differentiator.Bull.Seism.Soc.Am.，1992，82（1）：289－303.

［13］ 张中杰，滕吉文，杨顶辉.声波与弹性波场数值模拟中的褶积微分算子法.地震学报，1996，18（1）：63－69.Zhang Z J，Teng J W，Yang D H.The convolutional differentiator method for numerical modeling of acoustic and elastic wave－field.Acta Seismologica Sinica（in Chinese），1996，18（1）：63－69.

［14］ 戴志阳，孙建国，查显杰.地震波场模拟中的褶积微分算子法.吉林大学学报（地球科学版），2005，35（4）：520－524.Dai Z Y，Sun J G，Zha X J.Seismic wave field modeling with convolutional differentiator algorithm.Journal of Jilin University（Earth Science Edition）（in Chinese），2005，35（4）：520－524.

［15］ 龙桂华，李小凡，张美根.基于Shannon奇异核理论的褶积微分算子在地震波场模拟中的应用.地球物理学报，2009，52（4）：1014－1024.Long G H，Li X F，Zhang M G.The application of convolutional differentiator in seismic modeling based on Shannon singular kernel theory.Chinese J.Geophys.（in Chinese），2009，52（4）：1014－1024.

［16］ 李一琼，李小凡，朱童.基于辛格式奇异核褶积微分算子的地震标量波场模拟.地球物理学报，2011，54（7）：1827－1834.Li Y Q，Li X F，Zhu T.The seismic scalar wave field modeling by symplectic scheme and singular kernel convolutional differentiator.Chinese J.Geophys.（in Chinese），2011，54（7）：1827－1834.

［17］ Li X F，Li Y Q，Zhang M G，et al.Scalar seismic－wave equation modeling by a multisymplectic discrete singular convolution differentiator method.Bull.Seism.Soc.Am.，2011，101（4）：1710－1718.

［18］ Iwatsu R.Two new solutions to the third－order symplectic integration method.Phys.Lett.A.，2009，373（34）：3056－3060.

［19］ Wei G W.Discrete singular convolution for the solution of the Fokker－Planck equations.J.Chem.Phys.，1999，110（18）：8930－8942.

［20］ 黄志祥，沙威，吴先良等.高阶辛算法的稳定性与数值色散性分析.计算物理，2010，27（1）：82－88.Huang Z X，Sha W，Wu X L，et al.Stability and numerical dispersion of high order symplectic schemes.Chinese Journal of Computational Physics（in Chinese），2010，27（1）：82－88.

［21］ 董良国，马在田，曹景忠.一阶弹性波方程交错网格高阶差分解法稳定性研究.地球物理学报，2000，43（6）：856－864.Dong L G，Ma Z T，Cao J Z.A study on stability of the staggered－grid high－order difference method of first－order elastic wave equation.Chinese J.Geophys.（in Chinese），2000，43（6）：856－864.

［22］ Chen J B.Modeling the scalar wave equation with Nyström methods.Geophysics，2006，71（5）：151－158.

The elastic wave fields modeling by symplectic discrete singular convolution differentiator method

LI Yi－Qiong1，2，LI Xiao－Fan1，ZHANG Mei－Gen1
1 Institute of Geology and Geophysics，Chinese Academy of Sciences，Beijing100029，China
2 Institute of Geophysics，China Earthquake Administration，Beijing100081，China

In this paper，we introduce a structure－preserving method based on symplectic discrete singular convolution differentiator（SDSCD）for simulating elastic wave fields.In the method presented for solving elastic wave equations，physical space is discretized by singular convolution differentiator，whereas a symplectic difference scheme is used for the time discretization.The computational accuracy and stability of this method have been greatly improved compared with traditional pseudo－spectral method.Numerical results suggest the SDSCD algorithm can suppress effectively numerical dispersion，and it is suitable for modeling the large－scale and long－term seismic wave propagation.

Symplectic algorithm，Discrete singular convolution differentiator，Elastic wave fields modeling

10.6038／j.issn.0001－5733.2012.05.029

P631

2011－11－03，2012－03－28收修定稿

国家自然科学基金重点项目（40437018）、国家自然科学基金项目（40874024，41174047）、国家重点基础研究发展计划“973”计划（2007CB209603）和中国地震局地球物理研究所基本科研业务费专项（DQJB11B10）联合资助.

李一琼，女，满族，1979年生，2011年于中国科学院地质与地球物理研究所获博士学位，现为中国地震局地球物理研究所在站博士后，主要从事地震波传播理论和数值模拟的研究工作.E－mail：liyq＠mail.igcas.ac.cn，liyq＠cea－igp.ac.cn

李一琼，李小凡，张美根.基于辛格式离散奇异褶积微分算子的弹性波场模拟.地球物理学报，2012，55（5）：1725－1731，

10.6038／j.issn.0001－5733.2012.05.029.

Li Y Q，Li X F，Zhang M G.The elastic wave fields modeling by symplectic discrete singular convolution differentiator method.Chinese J.Geophys.（in Chinese），2012，55（5）：1725－1731，doi：10.6038／j.issn.0001－5733.2012.05.029.

（本文编辑 何 燕）