SU初学习之有限差分法地震数据正演模拟(一)

（一）预备知识：

（1）多维数组的一维表达

二维数组变化为一维数组
编辑
二维变一维的原理：把每一行的行首接到上一行的行尾，然后把剩下的元素按原来的顺序继续接在原来行首的后面
同理，三维变为一维的方法为：data[iy][iz][ix]变为data[iynznx+iz*nx+ix]

（2）数组的初始化

memset(data,0,nsize*sizeof(float));

（3）地震数据正演流程

设置地质构造模型——波动方程（常选用声波方程）——数学方法求解——计算机编程实现——获得模拟地震数据

（二）地震数据正演

（1）设计速度场

示例：设置一个两层模型，上层速度为2500m/s，下层速度为3000m/s，观测系统横坐标为检波点位置，纵坐标为深度。

计算机编程实现：

#include<stdio.h>
#include<par.h>
#include<su.h>
#include<math.h>
#include<segy.h>
#include"main_lib.c"//main_lib.c为创建的一个包含附属函数的文件
int main()//mian函数为设置速度场的主函数
{
    int nx=401,nz=301,ix,iz;//nx表示观测系统横坐标范围即网格列数，同理nz为观测网格列坐标
    float *data;//定义一个速度数组
    data=alloc1float(nz*nx);//为数组分配空间
    memset(data,0,nx*nz*sizeof(float));//安全起见将数组中元素初始化为0
    for(ix=0;ix<nx;ix++)//两个for循环定义为data数组写入值
    {
        for(iz=0;iz<nz;iz++)
        {
            if(iz<150)
                data[iz*nx+ix]=2500.0;
            else
                data[iz*nx+ix]=3000.0;
        }
    }
    output("velocity.bin",data,nx,nz);//output为main_lib中的函数，作用：创建一个文件将数组存入文件中
    free1float(data);//设置一个数组就必须要在结束时释放数组
}

void output(char *filename,float *data,int nx,int nz)
{
	FILE *fp;//定义一个文件
	if((fp=fopen(filename,"wr+"))==NULL)//该if-else组合是为方便检验错误，当文件出错时可以报错
	{	printf("can not find '*filename'\n");
		exit(0);
	}
	else
	{
		fwrite(data,nx*nz*sizeof(float),1,fp);//将数组data写入filename文件中
	}
	fclose(fp);//用完一个文件后一定要关闭文件！
}

写一个Makefile文件并生成可执行文件main.sh，执行后输出得到如下速度模型：

为方便观看，并且通常情况下应设置为上下分层而非左右分层，因此进行转置：
将data数组进行转置：改变一维数组的排列。此时要先输出行，再输出列，因此改变两个for循环为：

for(ix=0;ix<nx;ix++)//两个for循环定义为data数组写入值
    {
        for(iz=0;iz<nz;iz++)
        {
            if(iz<150)
                data[ix*nz+iz]=2500.0;
            else
                data[ix*nz+iz]=3000.0;
        }
    }

得到输出结果为：

（2）设置一个雷克子波

#include<stdio.h>
#include<par.h>
#include<su.h>
#include<math.h>
#include<segy.h>
#include"main_lib.c"

int main()
{
	float *data;
	int nt=501;//采样点个数
	float dt=0.001;//采样点间隔
	float fm=30.0;//雷克子波主频
	data=alloc1float(nt);//为数组分配空间
	memset(data,0,nt*sizeof(float));//初始化数组
	generate_wavelet(data,nt,fm,dt);//创建一个雷克子波
	output_wavelet("ricker.bin",data,nt);//输出雷克子波
	free1float(data);//释放数组

}

void generate_wavelet(float *data,int nt,float fm,float dt)
{
	float temp;
	int it;
	for(it=0;it<nt;it++)
	{
		temp=PI*PI*fm*fm*(it*1.0)*dt*(it*1.0)*dt;
		data[it]=(1.0-2.0*temp)*exp(-1.0*temp);//雷克子波计算公式
	}
}

void output_wavelet(char *filename,float *data,int nt)
{
	FILE *fp;
	if((fp=fopen(filename,"wr+"))==NULL)
	{	printf("can not find '*filename'\n");
		exit(0);
	}
	else
	{
		fwrite(data,nt*sizeof(float),1,fp);
	}
	fclose(fp);
}

在终端输入命令,生成雷克子波图像。

make
./wavelet
xgraph<ricker.bin n1=501 d1=0.01&

生成图像为雷克子波的一个周期，为了能够更清楚看清雷克子波的形状，可以通过调整输出周期对图像进行偏移

void generate_wavelet(float *data,int nt,float fm,float dt)
{
	float temp;
	int it,itshift;
	for(it=0;it<nt;it++)
	{
		itshift=it-(10.0/fm/dt);//左加右减，向右偏移10个周期
		temp=PI*PI*fm*fm*(itshift*1.0)*dt*(itshift*1.0)*dt;
		data[it]=(1.0-2.0*temp)*exp(-1.0*temp);
	}

偏移后结果如图所示：