有限差分法地震数据正演模拟（二）

（一）正演模拟

在（一）中已经学会了用C语言编程建立一个速度模型（velocity_model）和设定雷克子波（ricker），下面来学习正演，首先要明白正演的含义：

（1）概念：

正演模拟，是在地球物理勘探研究中，根据地质体的形状、产状和物性数据，通过构造数学模型计算得到其理论值（数学模拟），或通过构造实体模型来观测模型所产生的地球物理效应的数值（物理模拟）。

（2）意义：

在地球物理资料解释过程中，常常利用正演模拟结果与实际地球物理勘探资料进行比较，不断修正模型，使模拟结果与实际资料尽可能地接近，进而使解释结果更接近客观实际。

（3）实现过程：

在速度场中设置一个震源，根据波动方程来演绎波的传播情况。其结果也可以作为反演的验证。

（二）有限差分

声波波动方程：反应声波的传播规律，描述最简单的波现象的方程，p代表波场信息，v表示速度场信息因为速度只与速度场有关，所以只与x,z有关。

由导数定义可知，在点p(t,x,z)处对t进行求二阶偏导则需先求出该点处的左右临近点的一阶偏导：
左右临近点的一阶偏导分别为：


再由这两个一阶导数求出该点的二阶导数为：

p(t,x,z)对x，z求二阶偏导同理，整理得出该波动方程为：

整理后得到：

p(t+1,x,z)即为下一时刻的波长信息，该波动方程的意义就是通过前一时刻（p(t-1,x,z)）和现在时刻（p(t,x,z)）预测下一时刻波的传播情况。
为了计算结果更加精确因此引入了偶数阶高阶有限差分：差分阶数=参与的数据点数

（三）具体实现

下面来学习具体的实现步骤，为精简主函数的代码部分，将设置子函数在另一个c文件中用来调用。
Step1：建立两个c文件，mian.c用来存储主函数作为主程序，main_lib.c用来存储子函数为主函数提供服务
Step2：首先书写主函数main，需要包括：（1）设置速度场（2）设置震源，将波加入速度场中（3）用有限差分法根据波动方程预测波的传播情况。

#include<stdio.h>
#include<su.h>
#include<segy.h>
#include<math.h>
	float *wavelet,*velocity;//为震源和速度场信息设置数组
	float *wf1,*wf2,*wf3;//为波动方程设置数组，wf1,wf2,wf3分别是t-1，t，t+1时刻的波的位置
	float *coe;//为差分系数设置数组
	int nt,nx,nz,sx,sz,N;//sx，sz为震源的横纵坐标，N为差分阶数
	float fm,dt,dx,dz;
//需要注意，在主函数外定义的变量为全局变量，可以在main_lib.c中直接使用
#include"main_lib.c"//调用子函数所在文件
int main()
{
   
	dx=5.0;dz=5.0;
	nx=401;nz=301;nt=501;
	dt=0.001;
	sx=200;sz=150;
	fm=30.0;
	N=2;
	Alloc();//自定义函数，为数组分配空间并作初始化
	generate_velocity(velocity,nx,nz);//自定义函数，设置速度场
	generate_wavelet(wavelet,nt,dt,fm);//自定义函数，设置震源信息
	generate_coe(coe);//自定义函数，设定差分系数的值
	zhengyan_display(wf1,wf2,wf3,velocity,nx,nz,nt,dx,dz,dt,coe,N);//自定义方程，预测波的传播情况
	output_wavelet("ricker.bin",wavelet,nt);//自定义，将波的信息数组保存到文件中并输出
	output("velocity.bin",velocity,nx,nz);//自定义，将速度场的信息数组保存到文件中并输出