声波波动方程正演模拟程序
程序介绍：
第一部分：加载震源，此处选用雷克子波当作震源。

编写震源程序后，我将输出的数据复制，然后我用excel做成了图片，以检验程序编写是否正确。

以下为雷克子波公式部分的程序：
for(it=0;it<Nt;it++)
{
t1=it*dt;
t2=t1-t0;
source[it]=(1.0-2.0*pow(PI*fm*t1,2.0))*exp(-pow(PI*fm*t1,2.0));
fprintf(fp,"%.8f %.8f\n",t1,source[it]);
}
此处，为了成图完整，我用的是t2，而不是t1，也就是把雷克子波向右移动了一段距离，使主要部分都显示出来。

（频率采用的是30hz)
从图中可以看出程序是正确的，符合理论上雷克子波的波形。

第二部分：主程序，编写声波正演模拟算子。

首先定义了各种变量，然后指定震源位置，选择权系数，给速度赋值，然后是差分算子的编写，这是主要部分，最后再进行时间转换，即把n-1时刻的值给n时刻，把n时刻的值给n+1时刻。

此处，我编写的是均匀介质声波方程规则网格的正演模拟程序，时间导数采用二阶中心差分、空间导数为2N阶差分精度，网格大小为200*200，总时间为400。

第三部分：这一部分就是记录文件。

首先记录Un文件，然后记录record文件。

模型构建与试算：
1、我首先建立了一个均匀介质模型，首先利用不同时间，进行了数值模拟，得到波场快照如图所示：
100ms 200ms 300ms
此处，纵波速度为v=3000m/s。

模型大小为200×200，空间采样间隔为dx=dz=10m。

采用30Hz的雷克子波作为震源子波，时间采样间隔为1ms，图中可以看出，波场快照中的同相轴是圆形的，说明在均匀各向同性介质中，点源激发的波前面是一个圆，这与理论也是吻合的。

并且随着时间的增大，波前面的面积逐渐增大，说明地震波从震源中心向外传播。

2、我在建立的均匀模型的基础上，改变差分算子的精度，分别采用2阶、6阶、12阶精度进行试算。

时间统一采用300ms的时候。

得到的波长快照如下：
2阶精度6阶精度12阶精度
图中可以看出，在阶数较低时，出现很多同相轴，说明数值频散现象严重；随着算子阶数的增加，对于高阶差分算子来说，算子阶数越高，压制数值频散效果越好，精度越高。

3、最后，我又建立了一个层状介质模型，上层介质速度为v=2000m/s，下层介质速度为v=4000m/s。

模型大小为200×200，空间采样间隔为dx=dz=10m。

采用30Hz的雷克子波作为震源子波，震源位于模型(70,100)处，时间采样间隔为1ms。

采用12阶差分算子进行数值模拟。

结果如下：
100ms
200ms 300ms
图中可以看出，在未遇到界面前，地震波在均匀介质中的波前面一个圆。

当遇到地层界面之后，在界面处发生了反射、透射和折射现象.沿测线方向的单炮记录如下图所示。

记录中存在两条直线状的同相轴和一条近似双曲线的同相轴。

由于直达波的时距曲线是直线，因此两条直线同相轴对应直达波；由于反射波的时距曲线是近似双曲线，因此近似双曲线同相轴对应的是反射波。