for(i=0;i<I;i++) for(j=0;j<J;j++) u2[i][j]=u3[i][j];
}
main() {
FILE *fp; int i,j,k; float v[I][J],Dt[I][J];
//定义文件
double **u1,**u2,**u3,**u4,**u5,w[KMAX];
通过以上的计算，可以将 k=130 和 k=240 时波场的数据进行导出，用 Matlab 绘出图像。其中，在 Matlab 里用到的成像程序如下（这里只列举第一张波前快照的程序，第二张的与之相同）： clear all; load wavefront1.dat;
m=wavefront1(1); n=wavefront1(2); len=length(wavefront1)-2; nn=len/m; if n~=nn; disp('The data file length error!'); return; end
u k 1 i, j

2uik, j

uik,
1 j

v2t 2 h 2
{
1 12
[
uk i2,
j

uk i2,
j
]

4 3
[uik1,
j

uk i 1,
j
]

5 2
uik,
j
}

v 2t 2 h 2
{
1 12
[
uik,
j2

uik, j2 ]
4 3
[uik,
j
1
主要集中在前部，大致上 t>35 的时候 S(t)都为零值， 故在绘制 S(t)图像时 k 取 100，来比较清晰得呈现雷 克子波。本实验 不涉及吸 收边界条件 的使用。这样
声波方程数值模 拟所需的（ 1）震源函 数（2）地层
速度（波速）（3）边界条件均定义完全。
为求（1）式的数值解，必须将此式离散化（包 括时间离散和空间离散，这里 dt 取 0.002，dh 取 4， 即用有限差分来逼近导数，用差商代替微商。为此，
图4
图5
通过调节增益（x(j)=(i-1)*dx+wavefront((i-1)*n+j+2)*12*dx，dx 前面的 12），对波场
值放大一定倍数，可以使能量小的增大，使成像清晰。从下图中可清楚看到波前的图像，数层
的地震波随时间的推移由震源开始向外扩散。k=130 时，大约是时间为 260ms 时的图像，此时
关键词：声波方程 数值模拟 地震波 有限差分 引言：地震波场模拟即地震正演，是指已知模型结构，通过物理或数值计算的
方法模拟该地质结构下的地震波的传播，最终合成地震记录，也可以认为其是野 外数据采集过程的室内再现。物理模拟花费昂贵，人们一般采用比较经济的数值 模拟技术。地震波场数值模拟是在给定数学模型（如弹性波方程，声波方程等）、 震源和地下几何界面、物性参数（岩层密度、速度等）情况下，研究弹性波或声 波的传播规律。地震波波动方程数值模拟方法主要包括克希霍夫积分法、傅里叶 变换法、有限元法和有限差分法等。相对于上述几种方法，有限差分法是一种更 为快速有效的方法。虽然其精度比不上有限元法，但因其具有计算速度快，占用 内存较小的优点，在地震学界受到广泛的重视与应用。本实验给定大小两个模型， 对于小模型，整个区域的速度值可设为常数，即只有一种介质，将震源点放在模 型中间，分别记录两个时刻的波前快照（即该时刻区域内所有网格点的波场值）。
由上面的程序可得震源离散后的图像，如下图：
图2
小模型波前快照
选用一个小模型，即 200*200 的范围，视为完全弹性模型，为一层介质，在介质中取速 度为定值，本实验中取为 2500，将震源放在模型中央，分别记录两个时刻的波前快照，（即 区域内所有网格点的波场值）。第一时刻为地震波还未传播到边界上的某时刻，第二时刻为 地震波已经传播到边界上的某时刻，体会其人工边界反射。
别对应哪些波。
实验内容：
对于二维速度-深度模型，地下介质中地震波的传播规律可以近似地用声波方程描述：
2u t 2

v
2
(2u x 2Fra bibliotek2u z 2
)

S (t )
（1）
其中 v(x,z)是介质在点(x , z)处的纵波速度，本实验中 v 为常数 2500，u 为描述速度位或者压力的波场，
s(t) 为震源函数，本实验选用雷克子波（零相位子波）
e s(t) (2f / )2t2 cos 2ft
上式中，t 为时间，fm 为中心频率，一般取为 20-40HZ，本实验取 20，γ为控制频带宽度的参数， 一般取 3-5,这里取 4，在实际计算过程中，需把此震 源函数离散，参 与波场计 算。由图像 可知子波能量
图3
根据网格设计可知，网格间距为 10m，则 200*200 的矩阵实际上大致是个 2000*2000 平方米的区域， 这样震源在（100，100）即中心点，即在地质剖面上，宽度为 1000m，深度为 1000m，假设介质均匀各向 同性，且边界为自由界面，所以在介质中速度 v 以一常数 2500m/s 的传播速度传播，应该在 400ms 左右时 到达边界。

uik, j1 ]
5 2
uik,
j
}
s(t) * (i i0 ) * ( j j0 )
（2）
一般来讲，差分时差分 格式阶数越高， 得到的波场值精 度越高，但稳定性 会受影响，差分 格式的稳定 性不仅与差分格式本身有 关，也随着差分 精度的不同而不同 ，而且与网格步 长之比的大小有 关，值得指出
u1[i][j]=0;u2[i][j]=0;u3[i][j]=0; //对速度，u1，u2，u3 赋初值
}
for(k=0;k<KMAX;k++)
{
w[k]=exp(-pow(2*PI*FM/R,2)*pow(k*DT,2))*cos(2*PI*FM*k*DT);
}
//波场值计算
for(i=0;i<I;i++)
������������ 为 Nyquist 频率，一般取为主频的两倍，G 为每个波长所占的网格点数，时间、空间为两阶差分的 情况 G 取 8，而时间、空间为四阶差分的情况 G 取 4。
则由上面两式得 v∈[1000，3061.8]
由上面定义的数值，用 Matlab 将震源函数 s(t) 成图，下面是成图的源程序：
的地震波的振幅愈来愈小。
以下为小模型波前快照的源程序： /********************************** *声波方程正演模拟——小模型波前快照* ***********************************/ #include <stdio.h> #include <math.h> #include <stdlib.h> #define PI 3.14159265 #define FM 20 #define R 4 #define KMAX 1000 #define DT 0.002 #define I 200 #define J 200 #define SI 100 #define SJ 100 #define DH 10
第一时刻为地震波还未传播到边界上的某时刻，第二时刻为地震波已经传播到边
界上的某时刻，体会其人工边界反射；对于大模型，定义为水平层状速度模型（至
少两层）；做两个实验，一是将震源点放在区域表层任一点，记录下某些时刻的
波前快照，体会地震波在两种介质的分界面上传播规律；二是合成一个地震记录，
即记录下与震源同一深度点的各点所有时刻的波场值，并指出记录上的同向轴分
m[i]=(double*)malloc(y*sizeof(double));
return m;
}
void exchange(double **u1,double **u2,double **u3) //循环交换赋值 {
int i,j; for(i=0;i<I;i++)
for(j=0;j<J;j++) u1[i][j]=u2[i][j];
3
vmax * t / h
的是差分格式的稳定性与微分方程无关，本实验所采用的（2）式为条件稳定，稳定条件：
8，
（2）式由差分方程近似替代微分方程所得，所以如果空间和时间采样间隔不当，就会产生频散现象，导致
波形畸变，甚至派生出多个同相轴，为减少频散，Δh 需满足 Dablain 的经验公式： h vmin /(GfN )
声波方程有限差分法数值模拟
学 院：海洋地球科学学院 专业：地球信息科学与技术 年级：2012 级 姓 名：王昊 指 导教师：宋鹏
声波方程有限差分法数值模拟
2012 级 地球信息科学与技术 12040032037 王昊
摘要：本实验应用声波方程作为正演模拟的波动方程，将所提供震源函数离散
后绘图，并将给定两个二维速度-深度模型（一个小模型；一个大模型），绘出图 形。通过模拟地震波在介质中的传播，理解实际勘探中地震波在地层中的传播规 律，在模拟水平层状速度模型中，体会地震波在两种介质分界面的传播规律，并 能够从地震记录中识别出反射波，透射波，多次波，折射波和绕射波。并通过模 拟人工合成的地震记录，体会地震勘探基本原理和方法，验证地震波传播能量波 形变化趋势。
在小模型中，可将介质看做是单一介质，即在整个介质中，波的传播速度是相同的，在 这里设地震波速度为 2500m/s，以下是 Matlab 成图的源程序：
clf for x=1:200
f(x)=2500; end plot(f);axis([1,200,1,3000]);title('小模型的速度-深度模型'); figure; x=[0 0 200 200]; y=[0 200 200 0]; patch(x,y,'g'); title('小模型'); axis ij; 用 matlab 编译后得图像，见下图：