数值模拟实验一

数值模型模拟实验报告实验名称：地震记录数值模拟的褶积模型法实验学院：地球物理学院学号：2010050603xx教师：熊高军姓名：Blackheart--Mike日期：2010.6.13实验一一、实验题目地震记录数值模拟的褶积模型法二、实验目的掌握褶积模型基本理论、实现方法与程序编制，由褶积模型初步分析地震信号的分辨率问题。

三、原理公式1、褶积原理地震勘探的震源往往是带宽很宽的脉冲，在地下传播、反射、绕射到测线，传播经过中高频衰减，能量被吸收。

吸收过程可以看成滤波的过程，滤波可以用褶积完成。

在滤波中，反射系数与震源强弱关联，吸收作用与子波关联。

最简单的地震记录数值模拟，可以看成反射系数与子波的褶积。

通常，反射系数是脉冲，子波取雷克子波。

（1）雷克子波（2）反射系数：（3）褶积公式：数值模拟地震记录trace(t):trace(t) =rflct(t)wave(t)反射系数的参数由z变成了t，怎么实现？在简单水平层介质，分垂直和非垂直入射两种实现，分别如图1和图2所示。

1）垂直入射：2）非垂直入射：2、褶积方法（1）离散化（数值化）计算机数值模拟要求首先必须针对连续信号离散化处理。

反射系数在空间模型中存在，不同深度反射系数不同，是深度的函数。

子波是在时间记录上一延续定时间的信号，是时间的概念。

在离散化时，通过深度采样完成反射系数的离散化，通过时间采样完成子波的离散化。

如果记录是Trace(t)，则记录是时间的函数，以时间采样离散化。

时间采样间距以 t表示，深度采样间距以 z表示。

在做多道的数值模拟时，还有横向x的概念，横向采样间隔以 x表示。

离散化的实现：t=It× t；x=Ix× x；z=Iz× z或：It=t/ t; Ix=x/ x; Iz=z/ z（2）离散序列的褶积四、实验内容1、垂直入射地震记录数值模拟的褶积模型；2、非垂直入射地震记录数值模拟的褶积模型；3、点绕射的地震记录数值模拟的褶积模型；五、实验步骤（一）实验一：1、根据垂直入射褶积模型理论算法，填充程序（附后）的下划线部分，使程序完整，调试程序，算出结果，用“Fimage”显示软件显示褶积结果；2、根据非零偏移距算法，编制非零偏移距褶积模型程序，算出结果，用“Fimage”显示软件显示褶积结果。

（参考垂直入射褶积模型理论算法和程序，子波与反射层不变）；3、变换子波的长度Nw和主频fm，重复1和2；（1）Nw=80的情况下：1）fm=5，2）fm=10，3）fm=15，3）fm=20，3）fm=55；（2）Nw=160的情况下：1）fm=5，2）fm=10，3）fm=15，3）fm=20，3）fm=55；六、实验结果1、结果显示2、对比分析对于垂直入射而言，当子波长度保持不变而频率逐渐升高时画出的图像都在0.5下方而且一次变得细起来，而对于子波长度不同频率相同时子波长度打的图像位于0.5上方，同样子波长度不变频率增加时图像变细七、讨论建议1、实验收获：获得了垂直入射地震记录数值模拟的褶积模型和非垂直入射地震记录数值模拟的褶积模型2、存在问题：开始的时候Nt取值较小让其范围超出了最大的容纳范围，最后经过改正把Nt改的比原来大一些才得出了正确的结果，开始非垂直绘图出来都是直线，经检查是程序出现咯一点点小问题，但却花了很长时间才找出来错误3、心得体会：在做程序的时候一定要细心，一个小小的错误都会耽误自己做程序的效率，而且写程序的时候一定要对其这样发现错误检查起来也方便得多.实验二：两层反射界如图3所示为两层反射界面的模型，对比不同间距H12情况下，不同长子波长度Nw和主频fm的褶积模型。

1、H12=20m时：（1）Nw=80的情况下：1）fm=5，2）fm=10，3）fm=15，4）fm=20，5）fm=55；（2）Nw=160的情况下：1）fm=5，2）fm=10，3）fm=15，4）fm=20，5）fm=55；2、H12=40m 时：（1）Nw=80的情况下：1）fm=5，2）fm=10，3）fm=15，3）fm=20，3）fm=55；（2）Nw=160的情况下：1）fm=5，2）fm=10，3）fm=15，3）fm=20，3）fm=55；3、H12=60m时：（1）Nw=80的情况下：1）fm=5，2）fm=10，3）fm=15，3）fm=20，3）fm=55；（2）Nw=160的情况下：1）fm=5，2）fm=10，3）fm=15，3）fm=20，3）fm=55；4、H12=80m时：（1）Nw=80的情况下：1）fm=5，2）fm=10，3）fm=15，3）fm=20，3）fm=55；（2）Nw=160的情况下：1）fm=5，2）fm=10，3）fm=15，3）fm=20，3）fm=55；六、实验结果2、结果显示H=20H=40H=60H= 802、对比分析：从图中可以看出档Nw不变时频率一次增加时获得的图像变得细起来而且清晰度分辨率更好一些，当频率相同而Nw不同时图像分别位于0.5的上下方，变换子波的长度为64的位于下方，而变换子波长度为128的位于0.5上方七、讨论建议1、实验收获：获得了不同长子波长度Nw和主频fm的褶积模型。

2、存在问题：开始得出的图像始终一样，最终检查出来是程序的错误，并且开始的时候Nt取值较小让其范围超出了最大的容纳范围，最后经过改正把Nt改的比原来大一些才得出了正确的结果，开始非垂直绘图出来都是单直线，经检查是程序出现咯一点点小问题，但却花了很长时间才找出来错误3、心得体会：在做程序的时候一定要细心，一个小小的错误都会耽误自己做程序的效率，而且写程序的时候一定要对齐，这样发现错误检查起来也方便得多，实验的结果和理论结果有一定的误差，但在允许的范围，总的来说达到了想要的结果.实验三点绕射的褶积模型（1）地质模型与参数：如图4所示为点绕射模型（下图）及其地震记录（上图）。

参数为：Nx=128，Nz=200，V=3500m/s,h=1000m。

如图4图下图所示，地下绕射点位置为：(Ix0=NX/2,Iz0=1000/20)。

其他参数与前述同。

如图4所示，绕射点正上方检波器的自激自收世间为t0：t0=2h/v任意位置检波器的自激自收世间为ti：（3）褶积模型爆炸反射界面传播到检波点的脉冲值游侠是计算：与子波w(It)褶积,求得褶积模型。

计算公式如下：),(It Ix r),(It Ix p六、实验结果七、讨论建议1、实验收获：获得咯较理想的点绕射的褶积模型的图像2、存在的问题以及其它问题：点绕射它既是发射点也是接收点所以其中有很多其他波的影响3、心得体会：在做程序的时候一定要细心，注意自己的格式，一个小小的错误都会耽误自己做程序的效率，而且写程序的时候一定要对齐，这样发现错误检查起来也方便得多，实验的结果和理论结果有一定的误差，但在允许的范围，总的来说达到了想要的结果.附：程序编制与分析//1. 预处理部分#include<math.h>#include<stdio.h>#include<string.h>int Cnltn(float,float);int Rflct(float,float,float);int Wave(float,float);#define Nx 128#define Nt 256#define Nw 32//////////////////////////////////////////////////////////////////////////2. 主程序波分int main(){float dt=0.004,dx=20,fm=25,h=1000,v=3000;int iflag_Co,iflag_Re,iflag_Wv;if(iflag_Wv=Wave(fm,dt)!=1)printf("Wave is error");if(iflag_Re=Rflct(dt,h,v)!=1)printf("Reflection is error");if(iflag_Co=Cnltn(dt,dx)!=1)printf("Convosion is error"); }//////////////////////////////////////////////////////////////////////// 3.函数实现部分// 3.1 Wave Formaing functionint Wave(float fm,float dt){FILE *fpw;int It;float Wa[Nw],t;double pai=3.1415926;if((fpw=fopen("wave.dat","wb"))==NULL)printf("Connot open file ""wave""");for(It=0;It<Nw;It++){t=(float)(It+1)*dt;Wa[It]=(float)(1.0-2*pai*pai*fm*fm*t*t)*(float)exp(-pai*pai*fm*fm*t); //雷克子波的计算？fwrite(&Wa[It],sizeof(Wa[It]),1,fpw);}fclose(fpw);return(1);}/////////////////////////////////////// 3.2 Reflect Formaing functionint Rflct(float dt,float h,float v){FILE *fpr;int It,Ix,Ltdpth;float t;float Re[Nt];if((fpr=fopen("Reflect.dat","wb"))==NULL)printf("Connot open file ""Reflect""");for(Ix=0;Ix<Nx;Ix++){for(It=0;It<Nt;It++){Re[It]=0.;}t=(float)2.*h/v;Ltdpth=(int)(t/dt);Re[Ltdpth]=1.;for(It=0;It<Nt;It++){fwrite(&Re[It],sizeof(Re[It]),1,fpr);}}fclose(fpr);return(1);}/////////////////////////////////////// 3.2 Convolution functionint Cnltn(float dt,float dx){FILE *fpc,*fpw,*fpr;int It,Ix,Itao;float Wa1[Nw],Wa[Nw],Re[Nt+Nw+Nw],Re1[Nt];float Con[Nt+Nw];if((fpc=fopen("Convosion.dat","wb"))==NULL)printf("Connot open file ""Convosion""");if((fpw =fopen("wave.dat","rb"))==NULL)printf("Connot open file ""wave""");if((fpr =fopen("Reflect.dat","rb"))==NULL)printf("Connot open file ""Reflect""");for(Ix=1;Ix<2;Ix++){for(It=0;It<Nw;It++){fread(&Wa1[It],sizeof(Wa1[It]),1,fpw);Wa[Nw-1-It]=Wa1[It];}}fclose(fpw);for(Ix=0;Ix<Nx;Ix++)for(It=0;It<Nt;It++){fread(&Re1[It],sizeof(&Re1[It]),1,fpr);}for(It=0;It<Nt+2*Nw;It++){Re[It]=0.;}for(It=0;It<Nt;It++){Re[Nw+Nt]=Re1[Nt];//反射系数数据移动: "0 到Nt" 移为"Nw 到Nw+Nt"}for(It=0;It<Nt+Nw;It++){Con[It]=0;for(Itao=0;Itao<Nw;Itao++){Con[It]=Con[It]+Wa[Itao]*Re[It+Itao];//褶积运算//printf("%f\n",Con[It]);}for(It=Nw/2;It<Nt+Nw/2;It++){fwrite(&Con[It],sizeof(Con[It]),1,fpc); }}fclose(fpw);fclose(fpr);fclose(fpc);return(1);。