4波动方程有限差分的几个问题
4.1初始条件
在进行波动方程有限差分数值模拟时，初始时刻的波场值是不知道的，而在计算差分方程时，只有知道时间间隔前一时刻的值才能递推计算出后面时刻的值。因此在计算开始时就必须要考虑到初始条件。
一般假设在震源发生震动前，各个质点均处于静止状态：
初始时刻波场的速度也同样为零，即
4.2 震源函数
震源函数的给出通常有两种方式：一种是初值法，就是将理论结果当作初始值给出；另外一种是力源法，就是以力源的方式给出。初值法的震源除了自由表面或内界面附近，可以在模型的其他任意处进行定义；而力源法可以将震源定义在自由表面的附近，但必须是在网格点上。两种方法各有其优缺点[1,26]。
在进行波动方程数值模拟计算的过程中，震源问题的处理对模拟的结果有着重要的影响。子波是描述震源的时间延续特征的时间函数，对于地震子波而言，子波延续的时间越短，频带越宽，子波的垂直分辨率也就越高。由于做有限差分计算时会出现数值频散，尤其当空间采样不足时，子波的高频成分频散就会更严重。所以要根据模型的速度及网格间距合理选择子波主频。本文采用了雷克子波和高斯子波作为震源函数[27]。
12
4.3 差分方程的稳定性及收敛性
一个有实际应用意义的数值模拟算法必须是稳定的，也就是说对算法的数值解与解析解的差值是有限定的。对于有限差分法，一般来说如果差分方程的解的误差不随时间的推进而增加，那么就说该差分方程的解是稳定的，这也就是说差分方程的解是有界的。差分的稳定条件根据差分阶数的不同而有所不同，差分的阶数越高，稳定性越差。对于不同的差分格式，参数的选取一般也是不同的[5,28,29]。
对于均匀介质，差分方程的解的稳定性是确定的，稳定性条件如下式所示：
其中V是均匀介质中波的传播速度，应当注意，这里的Δx是差分格点中Δx与Δz的最小值。而对于非均匀介质应满足：
其中Vmax是各种均匀介质中的最大波速值。
4.4 频散问题
在进行地震波数值模拟时，我们希望尽可能同时达到较高的模拟精度和较快的计算速度。但是在用微分方程数值模拟的方法求解波动方程时，就会产生数值频散或网格频散。这种频散是用微分方程求解波动方程时固有的本质特征，是避免不了的，只能尽量减弱这种频散[2]。
通过对有限差分法数值频散的理论分析研究可以知道，影响频散的因素有很多。在有限差分的计算过程中，不仅需要对空间进行网格化离散，同时还需要对时间进行离散。波动方程经过空间离散化之后就引入了差分算子误差；而时间间隔的选择也对数值模拟的精度和数值频散有一定影响。当子波频率越高，则一个波场内的网格点数就越少，则频散现象越严重；对于一定频率的子波，介质的速度越低，则一个波场内的网格点数也少，数值频散现象也越严重。因此，从理论的角度来说，在进行数值模拟时，对于给定的子波频率，提高时间和空间的差分阶数，减小网格
13
步长和时间间隔都可以提高数值计算的精度和改善数值频散情况。但是过细的网格剖分会增大计算所需的存储量，从而增加计算的时间，所以就需要选取合适的参数，在保证精度和计算速度的同时尽量减少频散[29-33]。
5 有限差分算法中边界条件的处理
我们在计算机上进行有限差分数值模拟计算时，由于受到计算机内存和计算速度等因素的制约，只能在一个有限的区域内进行计算。这样就会在两侧与底界上产生边界，这个边界就是人工边界。当波传播到人工边界时就会产生反射，反射必然会造成一定的干扰，产生边界效应。从图5-1（左）中可以看出，在没有加入边界条件时，当波传播到边界时，在边界处产生了很强的反射，这是我们在进行数值模拟计算时所不期望出现的。
为了消除或减弱这种边界反射效应，得到地质地层真实的反射信息，就需要对人工边界进行处理，从而得到更接近于实际空间中波的传播规律。在处理边界反射时常用的边界条件方法有高衰减区法、傍轴近似法、吸收边界法和反周期扩展法[1,2,25]。
高衰减区法是在靠近边界处引入一个高吸收区，通过耗散因子使入射波仔这个吸收区域内逐渐衰减，从而抑制边界附近的人工反射。旁轴近似法是将边界区的波动方程用单向的傍轴近似波动方程代替，只模拟波能量由差分网络向外边界的单向传播，从而消除边界反射的问题，不过这个方法只适用于小角度入射的情况。吸收边界法是先导出吸收边界条件方程，然后让方程与计算区域内的波动方程联立求解，从而使从人工边界向计算区域反射回来的波能够全部或部分被吸收。反周期扩展法，即利用正反周期函数极性相反的特点消除回绕波场，这个方法在理论上是能够完美的消除回绕波场而不产生任何的反射，但是在实际应用时，由于这种方法需要进行反周期波场的扩展，大大的增加了计算量，而且为了得到最后叠加的平均波场，正反周期波场都必须存储，这样同时也会占用很大内存空间[25,34,35]。
下面将就本文将要用到的透明边界条件进行讨论。
14
5.1 一维透明边界条件
对于一维均匀介质波动方程：
式中 V 表示声波的传播速度，U 表示声压。我们引入一个向右行波将这个平面波的方程代入上面的均匀介质波动方程（5-1）中可以解得：
（k =ω/V），
其中r表示反射系数。在右边界x=a上反射波与入射波的振幅相等，即|r|=1。当引入一个左行波
时，左边界x=a上同样也会产生反射，且|r|=1。由此可以看出，如果边界条件选
取不当会使左右两侧的边界上产生强反射，因此就需要在左右边界处选取的边界条件能使r=0。
则右行波就满足方程
左行波就满足方程
由此可以看出，如果在右行波上加上一个左行波项
，则只有当r = 0时方程（5-5）才
成立，即当选式（5-5）作为右边界条件时，则在右边界上不会产生反射，因此可以将式上式作为右边界条件。同理可得左边界条件为
通常用右行波和左行波之和来表示波函数，因此也只有当 r = 0时，波函数才满足式（5-4），所以在引入边界条件以后左边界和右边界上都不会产生反射。
15
5.2 二维透明边界条件
设在平面区域为维波动方程
，对于已知某初始条件下的二
如果给定一个边界条件
则平面波在边界x = a、 x = a及 z = b处会产生反射。考虑一个向右传播的平面波
其中θ表示入射角，即平面波前与 x 轴的夹角。那么这时反射波为
其中r为边界产生的反射波的反射系数。则总波场可以表示为
以x=a处为例，把式（5-11）代入U(a,z,t)=0可得|r|=1，即在右边界上产了反射。由于对于波动方程（5-6）
当波只在水平方向传播时，即θ=0时
在右边界x = a，如果
16