陈 丽，朱兴文，张朝元*

（1.大理大学工程学院，云南大理 671003；2.大理大学数学与计算机学院，云南大理 671003）

开发精度高、效率高的模拟地震波计算方法是目前涉及地震波动方程开展反演的一个重要方向。过去的计算格式〔1-5〕在进行地震波模拟时会造成严重的数值频散现象而达不到当前规模大的地震波模拟效度，杨顶辉教授团队于2003年首次在地震波正演模拟中引入近似解析离散化（nearly analytic discrete，NAD）算子〔6〕，目前已获得系列NAD算子的数值算法〔7-10〕，这些算法的模拟频散效果均很不错，然而这些算法仅为四阶的空间精度。

为提高地震波计算模拟效度，本文针对具有哈密尔顿系统的二维弹性波方程，结合离散空间高阶偏导数的八阶NAD算子和离散时间导数的二阶辛分部Runge-Kutta算法〔4，11〕，得到了八阶NADSPRK算法。针对该方法，从理论和数值计算两方面研究了其稳定性条件、数值频散和计算效率。结果表明：同四阶NSPRK算法〔11〕、八阶Lax-Wendroff（LWC）算法〔12〕和八阶交错网格（SG）算法〔12〕相比，八阶NAD-SPRK算法压制数值频散的能力显著优于传统数值计算方法，且具有最小的数值误差和最高的计算效率。

1 方法推导

设二维弹性波方程为：

其中，j=1，3，fi和ui分别为i方向的力源分量和位移分量，ρ=ρ（x，z）为介质密度，σij为应力张量。

由应力与张力之间的物理关系，方程（1）式的向量方程为：

记vi=∂ui/∂t（i=1，2，3），则v=（v1，v2，v3）T，方程（2）式为：

其中D称为偏微分算子矩阵，且

O为三阶零方阵，I为三阶单位矩阵。

此时，表达式（4）拥有哈密尔顿系统模式，故采纳哈密尔顿系统模式的求解流程对二维弹性波方程表达式（1）的对应方程（4）式进行求解。……