地震记录数值模拟的褶积模型法

实验课程数值模型模拟学院名称地球物理学院专业名称勘测技术与工程学生姓名学生学号指导教师熊高君实验地点5417实验成绩2015年5月成都理工大学地震数值模拟》实验报告实验时间2015年5月开课单位地球物理学院指导教师熊高君实验题目：地震记录数值模拟的褶积模型法姓名学号班级专业勘测技术与工程（石油物探）院（系）地球物理学院地球探测与信息技术系单项成绩内容理解写作结构程序设计模型设计计算结果结果分析总成绩实验报告实验题目：地震记录数值模拟的褶积模型法二'实验目的：掌握褶积模型基本理论、实现方法与程序编制，由褶积模型初步分析地震信号的分辨率问题。三、原理公式1、褶积原理地震勘探的震源往往是带宽很宽的脉冲，在地下传播、反射、绕射到测线，传播经过中高频衰减，能量被吸收。吸收过程可以看成滤波的过程，滤波可以用褶积完成。在滤波中，反射系数与震源强弱关联，吸收作用与子波关联。最简单的地震记录数值模拟，可以看成反射系数与子波的褶积。通常，反射系数是脉冲，子波取雷克子波。雷克子波：wave(t)=cos(2ft)*反射系数：⑶褶积公式：数值模拟地震记录trace(t):trace(t)=rflct(t)*wave(t)；反射系数的参数由z变成了t,怎么实现？在简单水平层介质，分垂直和非垂直入射两种实现，分别如图1和图2所示。I图1图?1）垂直入射：t=2h/v；2）非垂直入射:2、褶积方法（1）离散化（数值化）计算机数值模拟要求首先必须针对连续信号离散化处理。反射系数在空间模型中存在，不同深度反射系数不同，是深度的函数。子波是在时间记录上一延续定时间的信号，是时间的概念。在离散化时，通过深度采样完成反射系数的离散化，通过时间采样完成子波的离散化。如果记录是Trace（t）,则记录是时间的函数，以时间采样离散化。时间采样间距以At表示，深度采样间距以Az表示。在做多道的数值模拟时，还有横向Ax的概念，横向采样间隔以Ax表示。离散化的实现：t=ItXAt；x=IxXAx；z=IzXAz；或：It=t/At:Ix=x/Ax;Iz=z/Az（2）离散序列的褶积trace(It)=四、实验内容1、垂直入射地震记录数值模拟的褶积模型；2、非垂直入射地震记录数值模拟的褶积模型。五、方法路线1、根据垂直入射褶积模型理论算法，填充程序（附后）的下划线部分，使程序完整，调试程序，算出结果，用“Fimage”显示软件显示褶积结果；2、根据非零偏移距算法，编制非零偏移距褶积模型程序，算出结果，用“Fimage”显示软件显示褶积结果。（参考垂直入射褶积模型理论算法和程序，子波与反射层不变）；3、变换子波的主频：fm（10hz到300hz范围），重复1和2；4、变换子波的长度：Nw（80ms到160ms范围），重复1和2；5、改变反射层深度：h（800m到1600m范围），重复1和2；6、改变介质速度：v（2000m/s到7000m/s范围），重复1和2。六、实验结果图3—1Nw=32.h=1000.v=3000,fm=100地震记录数值模拟的褶积模型（左）和子波（右）图3—2Nw=32.h=1000.v=3000,fm=200地震记录数值模拟的褶积模型（左）和子波（右）图3—3Nw=32.h=1000.v=3000,fm=300地震记录数值模拟的褶积模型（左）和子波（右）图3—4h=1000.v=3000,fm=25,Nw=20地震记录数值模拟的褶积模型（左）和子波（右）图3—621000,v二3000,和P25.Nw二40地震记录数值模拟的褶积模型（左）和子波（右）图3—8Nw=32,v二3000,和P25,h=1200地震记录数值模拟的褶积模型图3—9Nw二32,v二3000,和P25,h二1400地震记录数值模拟的褶积模型图3—11Nw=32,h=1000,和二25,v=4000地震记录数值模拟的褶积模型图3—12Nw二32,柘1000,和二25,v二6000地震记录数值模拟的褶积模型图4—1Nw=32,h=1000,v=3000.fm=100地震记录数值模拟的褶积模型（左）和子波（右）图4—2Nw二32.hH000,v=3000,fm=200地雇记录数值模拟的褶积模型（左）和子波（右）图4—3Nw二32,h-1000,v-3000.fm-300地蔑记录数值模拟的褶积模型（左）和子波（右）图4—4h=1000,v=3000.fm=25,Nw=20地震记录数值模拟的褶积模型（左）和子波（右）图4—5h=1000,v=3000,吊二25,Nw-30地装记录数值模拟的褶积模型（左）和子波（右）图4—6h=1000.v=3000.fm-25,Nw=40地震记录数值模拟的褶积模型（左）和子波（右）图4—7Nw=32.v=3000,fm=25,h=1000地震记录数值模拟的褶积模型

图4—8Nw=32,v二3000,fm二25,h二1200地震记录数值模拟的褶积模型图4—9Nw二32,v二3000,fm二25,h=1400地震记录数值模拟的褶积模型图4—10Nw=32.h=1000.fm=25,v=2000地震记录数值模拟的褶积模型

51图4—11Nw=32,h=1000,fm=25,v=4000地震记录数值模拟的褶积模型图4—12Nw二32,h=1000,fm二25,v二6000地震记录数值模拟的褶积模型子波振幅谱：图5—1fm=100的子波振幅谱

2、1.2图5—2fm=100的子波振幅谱对比分析a)由图3—1、图3—2、图3—3（或图4—1、图4—2、图4—3）可知，当不变，子波频子波长度（Nw=32）、深度（h=1000）、速度（v=3000）率变化时，褶积模型不变，且均在0.5到Is2、1.2图5—2fm=100的子波振幅谱对比分析a)不变，子波频b)由图3—4、图3—5、图3—6可知，当深度（h=1000）、速度且垂直入射（v=3000）、子波频率（fm=25）不变，子波长度变化,时，褶积模型为直线模型，是因为垂直入射时，时间与深度为线性关系;由图4一4、图4-5.图4-6可知，当深度（h=1000）、速度且垂直入射（V=3000）、子波频率（fm=25）不变，子波长度变化，且非垂直入射时，褶积模型前半部分为双曲线模型，后半部分为直线模型，是由于计算的褶积结果的实际长度小于所取的长度，计算机赋的随机数所致，所以，非垂直入射时的褶积模型只有图件上显示的前半部分的双曲线，是由于非垂直入射时，时间与深度为双曲线关系;时，褶积模型前半部分为双曲线模型，后半部分为直线模型，是由于计算的褶积结果的实际长度小于所取的长度，计算机赋的随机数所致，所以，非垂直入射时的褶积模型只有图件上显示的前半部分的双曲线，是由于非垂直入射时，时间与深度为双曲线关系;d)由图3—7、图3-8、图3-9,当子波长度（Nw=32）、速度（v-3000）、子波频率（fm=25）不变,深度变大时，垂直入射时，图件上显示的褶积模型的位置逐渐向下移，即地震波的旅行路程变大，旅行时变长;d)e）由图4—7、图4—8、图4—9可知，当子波长度（Nw=32）、速度（v=3000）、子波频率（fm=25）不变,深度变大时，非垂直入射时，图件上显示的褶积模型的位置不变，这是由于非垂直入射时，相当于地下有一半圆形界面，在圆心处自激自收；f）由图3—10、图3—11、图3—12（或图4—10、图4—11、图4—12）可知，当子波长度（Nw=32）、深度（h=1000）、子波频率（fm=25）不变，速度变大时，图件上显示的褶积模型的位置逐渐向上移，这是由于随着速度变大，地震波的旅行时变小；g）由图5—1与图5—2可知，当子波频率变化时，其振幅谱不变。七、讨论建议1、实验收获通过此次试验，初步掌握了褶积模型基本理论、实现方法与程序编制，由褶积模型初步分析地震信号的分辨率问题，掌握了褶积模型与子波主频、子波长度、界面深度、介质中地震波速度的关系。2、存在问题对褶积模型的理论实现过程不是十分清楚，对结果的物理含义理解不够深入。3、其他问题由于不同的计算机，计算精度不一样，计算得到的数据结果可能会有部分差异，但总体趋势不变。4、心得体会在此次试验中，应特别注意褶积模型的理论实现过程的理解，以及程序调试时，要特别仔细地去检查每一个错误，每修改一处错误，就重新运行一次程序。附程序代码://1.预处理部分#include<math.h>#include<stdio.h>#include<string.h>floatCnltn(float,float);floatRflet(float,float,float);floatWave(float,float);SdefineNx128SdefineNt256SdefineNw32SdefinePI3.1415926//==================2.主程序波分=======voidmain(){//floatdt=0.004,dx=20,fm=25,h=1000,v二3000;intiflag_Co,iflag_Re,iflag_Wv;if(iflag_Wv=Wave(fm,dt)!=1)printf("Waveiserror");if(iflag_Re=Rflet(dt,h,v)!=1)printf(，zReflectioniserror");if(iflag_Co=Cnltn(dt,dx)!=1)printf("Convosioniserror");}//=================3函数实现部分=========//=======3.1WaveFormaingfunctionfloatWave(floatfm,floatdt)FILE*fpw;intIt:floatWa[Nw],t;if((fpw=fopen("wave.dat〃，〃wb〃))=NULL)printf("Connotopenfile〃〃wave〃〃〃)；for(It=0:It<Nw;It++)t二It*dt;Wa[It]=cos(2*PI*fm*t)*exp(-2*PI*PI*fni*fni*t);//形成子波fwrite(&Wa[It]>sizeof(Wa[It]),1,fpw);}fclose(fpw);return(1);)//======3.2ReflectFormaingfunction=======〃floatRflet(floatdt,floath,floatv){FILE*fpr;intIt,lx,J,Ltdpth;floatt,dx=20,x;floatRe[Nt];printf(，z请输入J:\n〃);scanf(〃％d",&J);if((fpr=fopen(^Reflect.dat〃，〃wb〃))=NULL)printf(^Connotopenfile〃〃Reflect〃〃〃)；for(Ix=0;Ix<Nx;Ix++)(for(It=0;It<Nt;It++){Re[It]=0.;)if(J==Dt=2*h/v;//垂直入射反射界面由深度转换为自激自收时间if(J==2){x=Ix*dx;t=2*sqrt(h*h+x*x)/v;//非垂直入射反射界面由深度转换为自激自收时间)if((J!=l)&&(J!=2))printfC输入错误\n〃);Ltdpth=(int)(t/dt);Re[Ltdpth]=l;fwrite(&Re[It],sizeof(Re[It]),1,fpr);}fclose(fpr);return(1);}//=======3.3Convolutionfunction=======floatCnltn(floatdt,floatdx)(FILE*fpc,*fpw,*fpr;intIt,lx,Itao;floatWai[Nw],Wa[Nw],Re[Nt+Nw+Nw],Rel[Nt],t;floatCon[Nt+Nw];if((fpc=fopen("Convosion.dat”,〃wb〃))==NULL)printf(^Connotopenfile〃〃Convosion〃〃")；if((fpw=fopen("wave.dat〃,〃:rb"))==NULL)printf(^Connotopenfile〃〃wave〃〃〃)；if((fpr二fopen("Reflect.dat〃，〃:rb〃))=NULL)printf(^Connotopenfile〃〃Reflect”〃〃)；for(Ix=l;Ix<2;Ix++)fread(&Wal[It],sizeof(Wai[It]),1,fpw);for(It=0;It<Nw;It++)Wa[It]=Wal[Nw-It-l];//褶积前子波准备}fclose(fpw);for(Ix=0;Ix<Nx;Ix++)(for(It=0;It<Nt;It++){fread(&Rel[It],sizeof(&Rel[It]),1,fpr);)