地震波波动方程数值模拟方法

地震波波动方程数值模拟方法主要包括克希霍夫积分法、傅里叶变换法、有限元法和有限差分法等。

克希霍夫积分法引入射线追踪过程，本质上是波动方程积分解的一个数值计算，在某种程度上相当于绕射叠加。该方法计算速度较快，但由于射线追踪中存在着诸如焦散、多重路径等问题，故其一般只能适合于较简单的模型，难以模拟复杂地层的波场信息。

傅里叶变换法是利用空间的全部信息对波场函数进行三角函数插值，能更加精确地模拟地震波的传播规律，同时，利用快速傅里叶变换(FFT)进行计算，还可以提高运算效率，其主要优点是精度高，占用内存小，但缺点是计算速度较慢，对模型的适用性差，尤其是不适应于速度横向变化剧烈的模型.

波动方程有限元法的做法是：将变分法用于单元分析，得到单元矩阵，然后将单元矩阵总体求和得到总体矩阵，最后求解总体矩阵得到波动方程的数值解；其主要优点是理论上可适宜于任意地质体形态的模型，保证复杂地层形态模拟的逼真性，达到很高的计算精度，但有限元法的主要问题是占用内存和运算量均较大，不适用于大规模模拟，因此该方法在地震波勘探中尚未得到广泛地应用。。

相对于上述几种方法，有限差分法是一种更为快速有效的方法。虽然其精度比不上有限元法，但因其具有计算速度快，占用内存较小的优点，在地震学界受到广泛的重视与应用。

声波方程的有限差分法数值模拟

对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述：

22u2u2uv(22)S(t) （4-1） t2xzv(x,z)是介质在点（x , z）处的纵波速度，u为描述速度位或者压力的波场，s(t)为震

源函数。

为求式（4-1）的数值解，必须将此式离散化，即用有限差分来逼近导数，用差商代替微商。为此，先把空间模型网格化（如图4-1所示）。

i1,ji,ji1,ji2,jz 设x、z方向的网格间隔长度为h，t为时间采样步长，则有：

xih （i为正整数）

zjh (j为正整数)tnt (n为正整数)

uik,j 表示在(i,j)点，k时刻的波场值。

1将uik,j在(i,j)点k时刻用Taylor展式展开：

uk1i,juki,jut12utk*t*t2t2tk*t*t2o(t2)

(4-2)

1将uik,j在(i,j)点k时刻用Taylor展式展开：

uk1i,juki,juttk*t12u*t2t2tk*t*t2o(t2)

(4-3)

1）将上两式相加，略去高阶小量，整理得(i,j)点k时刻的二阶时间微商为：

k1kk12uui,j2ui,jui,j t2t2 (4-4)

同理可得(i,j)点k时刻的二阶空间微商分别为：

u2 x2uik1,j2uik,juik1,jh2 (4-5)

(4-6)

这就实现了用网个点波场值的差商代替了偏微分方程的微商，将上三个式子代入(4-1)式中得：

k1(v*t/h)2*(ukukukuk)ui,ji1,ji1,ji,j1i,j1

2h24v2t2kk1(4-7) ()uui,ji,js(k)*(ii0)*(jj0)2h

1k1k1k12）采用四阶精度差分格式，（以X方向为例）即将uik2,j、ui1,j、ui1,j、ui1,j分别在(i,j)

kkk2uui,j12ui,jui,j12zh2点k时刻展开到四阶小量，消除四阶小量并解出二阶微分得：

2u11k4k5kkk{[uu][uu]ui,j} (4-8) i2,ji2,ji1,ji1,j32x2x2122u11k4k5kkk{[uu][uu]ui,j} i,j2i,j2i,j1i,j1221232zz

(4-9)

这就实现了用网个点波场值的差商代替了偏微分方程的微商，代入(4-1)式中得可得时间

二阶、空间四阶精度的声波方程差分格式为：：

uk1i,j2uki,juk1i,jv2t21k4k5kkk{[uu][uu]ui,j} i2,ji2,ji1,ji1,j1232h2

v2t21k4k5kkk{[uu][uu]ui,j} i2,ji2,ji1,ji1,j1232h2s(t)*(ii0)**(jj0) （4-10）

式中v(i,j)为介质速度的空间离散值，h是空间离散步长，t为时间离散步长，s(k)为震源函数，关于s(k)一般使用一个理论的雷克型子波代替，即：

s(t)e(2f/)tcos2ft （4-11）

上式中，t为时间， f为中心频率，一般取为20-40HZ，为控制频带宽度的参数，一般取3-5。在实际计算过程中，需把此震源函数离散，参与波场计算。

22(ii0)*(jj0)确定震源位置。

稳定性条件：对于特定的偏微分方程只有特定的几种有限差分格式是无条件或有条件稳定的，（4-7）、（4-10）式即是已被证明的有条件稳定格式，其稳定性条件分别为：

vmax*t/h2/23vmax*t/h

8 (4-12)

这里vmax表示的是地下介质的最大波速；若地下介质网格间隔、最小速度、及时间采样间隔不符合(4-12)式时，递推求解（4-10）式，波场值会出现误差（高阶小量）累积，出现不稳定现象。 频散关系式：

hvmin/(GfN) (4-13)

同时，在差分计算过程中，如果空间和时间采样间隔不当，就会导致波形畸变，甚至派生出多个同相轴，这种现象称为频散现象。

偏微分方程本身没有频散，网格频散是由于差分方程近似替代微分方程引起的。当波场按照波动方程所表示的微分方程传播时，波场的传播速度就是波动方程中的速度，但当波场按照波动方程离散化后的差分方程传播时，波场的传播速度就不再是波动方程中的速度了，而是与波的频率和波数有关的函数，具有不同频率和波数的波有不同的传播速度，因而在传播过程中会出现频散，发生畸变，且随走时的增加而增加。

式中vmin为最小速度，fN为Nyquist频率。一般取震源子波中的主频f的2倍值参与计算，G为每个波长所占的网格点数，对于空间二阶差分、时间二阶差分G取8，而对于空间为四阶差分的情况则G取4方能有效减少频散。

边界条件：在地震波场正演模拟中，必须引入人工边界来界定计算区域。人工边界若不做特殊处理，就会随着波场的递推计算在边界上产生虚假反射波从而扰乱波场，人工边界的处理是地震波场正演数值模拟的一个重要课题。