导数，用差分方法来计算时间的导数，可以看成是一种无限阶的有限差分法，是传统有限差分法的一个推广。伪谱法在粗网格上也能实现高精度的计算，相对有限元法实现起来较容易，在非线性波动问题及气象预测等领域有着广泛的应用[8-10]。
有限差分法是一种最常用的数值模拟方法，它是将波动方程中波场函数的空间导数和时间导数用相应的空间、时间的差分代替。有限差分法具有计算速度快、占用内存小等优点，该方法对于近远场及复杂边界都有广泛的适用性，能够准确地模拟波在各种介质及复杂结构地层中的传播规律。有限差分法方法简单、高效的优点是其他方法难以比拟的，因此有限差分法目前仍然是勘探地震学中应用最广泛的数值计算方法[11]。
2 波动方程有限差分法数值模拟
2.1 有限差分原理
有限差分法就是把求解区域划分为差分网格，然后用有限的网格节点代替连续的求解域，利用微商与差商的近似关系将描述介质传播的微分方程转化为差分方程进行求解。差分离散的方法有两种，一种是将单变量的二阶波动方程直接转化为时间空间的二阶中心差分进行离散求解；第二种方法是把用位移表示的二阶波动方程转化为用应力及质点速度表示的一阶方程组，然后用应力和速度的交错网格求解[,8]。
构造差分的方法也有很多种，一般采用泰勒级数展开法。常见的差分格式有显示差分、隐式差分和显隐交替格式；按照差分的精度又可以分为一阶差分、二阶差分和高阶差分等，目前通常见到的差分格式主要是几种差分格式的组合。我们知道，差分方程的建立首先选择网格布局和差分形式，然后以有限差分代替无限微分，以差商代替微商，并以差分方程代替微分方程及其边界条件，最后建立差分方程[12]。
在建立差分方程的时候要注意到两个方面。一是合理选择网格布局与步长。我们将离散后各相邻离散点之间的距离或者离散化单元的长度称为步长。在所选定的区域内进行网格划分是差分方程建立的第一步，其方法比较灵活，但是实际应用中往往遵循误差最小原则。常用的典型的网格剖分方法有矩形网格、三角形网格等，网格样式的选择和区域的形状有关。其次是将微分方程
2
转化为差分方程。这个过程就是选择一种差分格式代替微分形式的过程。构造差分的方法有多种形式，本文采用的是泰勒级数展开法[2,5,13]。
地震波场数值模拟以地震波动理论为基础，用有限差分法解波动方程时，对变量离散化，也即对连续的物理量只考虑其在离散的空间位置和离散时刻的值，然后把方程中的导数用这些离散的采样值表示[14]。对于一个单变量的函数f(x)，将其离散化，那么在采样点x =lΔx的采样值就是f(lΔx),其中Δx表示步长，l为整数。则有限差分法中f(x)在采样点x=lΔx的导数就可以近似表示为：
其中an是系数，N表示差分格式的长度，差分的格式是由这两个系数来决定。在实际应用中常用的差分格式有向前差分、向后差分以及中心差分： （1）向前差分：
（2）向后差分：
（3）中心差分：
2. 2 波动方程的建立
建立波动方程是在已知物体形状、位置、弹性常数及外力分布等参数的情况下求出物体的位移、应力与应变的分布。这个问题具体包括以下内容： （1）应力部分的平衡方程
3
（2）应变部分
（3）应力与应变的关系
将式（2-7）代入（2-6），推导可以得下式：
方程（2-8）是物体在平衡状态下的平衡方程。当物体处于不平衡状态时，式（2-8）则变为：
4
在二维情况下，我们只需要考虑x，z方向的位移分量，这时式（2-9）就可转换为关于x,z的方程：
当fx=fz=0时，上面的这两个方程就变为：
（2-11）就是二维非均匀介质弹性波方程。而在均匀介质中，λ、μ和ρ均为常数，则二维均匀介
质弹性波方程就为：
把（2-12）中的x和z分量合并就得到了二维均匀介质弹性波方程的矢量形 式表达式
5
在没有弹性横波只有弹性纵波存在时，对（2-13）两边取散度：
利用旋度与散度的关系，交换上式的微分次序并化简可得：
其中
表示纵波速度。根据赫姆霍兹在他著名的有关涡流运动的著作中证
明了下属定理：任何向量点函数，若它的散度和旋度具有位，则它可以表示为一个无旋部分和一个旋转部分之和。亦即对于任何一种向量场，如果在其定义域内有散度和旋度，则该向量场可表示为一个标量位的梯度场和一个向量位的旋度场之和，而空间传播的波动正是无旋运动和旋转运动这两种运动之和[15-20]
令
代入式（2-8）中，取其中的无旋部分则有：
然后再消去式中的梯度,就得到了用位移位表示的纵波方程：
3二维声波方程有限差分格式的建立
3.1 声波方程
一般地，二维非均匀介质的声波波动方程可表示为：
式中 U = U ( x , z , t)表示声压，V表示声波在介质中的传播速度，ρ是密度，它是随空间各点而变
6