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

本科生实验报告实验课程 数值模型模拟 学院名称 地球物理学院 专业名称 勘测技术与工程 学生姓名 学生学号 指导教师 熊高君 实验地点 5417 实验成绩 2015年5月 成都理工大学地震数值模拟实验报告实验时间2015年 5月开课单位地球物理学院指导教师熊高君实验题目：地震记录数值模拟的褶积模型法姓名学号班级专业勘测技术与工程（石油物探）院（系）地球物理学院 地球探测与信息技术系单项成绩内容理解写作结构程序设计模型设计计算结果结果分析总成绩实验报告一、 实验题目：地震记录数值模拟的褶积模型法二、 实验目的：掌握褶积模型基本理论、实现方法与程序编制，由褶积模型初步分析地震信号的分辨率问题。三、原理公式1、褶积原理地震勘探的震源往往是带宽很宽的脉冲，在地下传播、反射、绕射到测线，传播经过中 高频衰减，能量被吸收。吸收过程可以看成滤波的过程，滤波可以用褶积完成。在滤波中， 反射系数与震源强弱关联，吸收作用与子波关联。最简单的地震记录数值模拟，可以看成反射系数与子波的褶积。通常，反射系数是脉冲，子波取雷克子波。（1） 雷克子波：wave(t)=cos(2ft)*e-22f2t（2） 反射系数：（3） 褶积公式：数值模拟地震记录trace(t): trace(t) =rflct(t)*wave(t)；反射系数的参数由 z 变成了 t，怎么实现？在简单水平层介质，分垂直和非垂直入射两种实现，分别如图 1 和图 2 所示。图1 图21） 垂直入射：t=2h/v；2）非垂直入射：t=2h2+x2v2、褶积方法（1） 离散化（数值化）计算机数值模拟要求首先必须针对连续信号离散化处理。反射系数在空间模型中存在，不同深度反射系数不同，是深度的函数。子波是在时间记录上一延续定时间的信号，是时间的概念。在离散化时，通过深度采样完成反射系数的离散化，通过时间采样完成子波的离散化。如果记录是 Trace(t)，则记录是时间的函数，以时间采样离散化。时间采样间距以t 表示，深度采样间距以z 表示。在做多道的数值模拟时，还有横向x 的概念，横向采样间隔以x 表示。离散化的实现：t=Itt；x=Ixx；z=Izz；或：It=t/t; Ix=x/x; Iz=z/z（2） 离散序列的褶积trace(It)= Itao-+rflct Itao* wave(It -Itao)四、实验内容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。六、实验结果1、结果显示1）垂直入射图31 Nw=32,h=1000,v=3000,fm=100地震记录数值模拟的褶积模型（左）和子波（右）图32 Nw=32,h=1000,v=3000,fm=200地震记录数值模拟的褶积模型（左）和子波（右）图33 Nw=32,h=1000,v=3000,fm=300地震记录数值模拟的褶积模型（左）和子波（右）图34 h=1000,v=3000,fm=25，Nw=20地震记录数值模拟的褶积模型（左）和子波（右）图35 h=1000,v=3000,fm=25，Nw=30地震记录数值模拟的褶积模型（左）和子波（右）图36 h=1000,v=3000,fm=25，Nw=40地震记录数值模拟的褶积模型（左）和子波（右）图37 Nw=32,v=3000,fm=25，h=1000地震记录数值模拟的褶积模型图38 Nw=32,v=3000,fm=25，h=1200地震记录数值模拟的褶积模型图39 Nw=32,v=3000,fm=25，h=1400地震记录数值模拟的褶积模型图310 Nw=32, h=1000,fm=25，v=2000地震记录数值模拟的褶积模型图311 Nw=32, h=1000,fm=25，v=4000地震记录数值模拟的褶积模型图312 Nw=32, h=1000,fm=25，v=6000地震记录数值模拟的褶积模型2）非垂直入射图41 Nw=32,h=1000,v=3000,fm=100地震记录数值模拟的褶积模型（左）和子波（右）图42 Nw=32,h=1000,v=3000,fm=200地震记录数值模拟的褶积模型（左）和子波（右）图43 Nw=32,h=1000,v=3000,fm=300地震记录数值模拟的褶积模型（左）和子波（右）图44 h=1000,v=3000,fm=25，Nw=20地震记录数值模拟的褶积模型（左）和子波（右）图45 h=1000,v=3000,fm=25，Nw=30地震记录数值模拟的褶积模型（左）和子波（右）图46 h=1000,v=3000,fm=25，Nw=40地震记录数值模拟的褶积模型（左）和子波（右）图47 Nw=32,v=3000,fm=25，h=1000地震记录数值模拟的褶积模型图48 Nw=32,v=3000,fm=25，h=1200地震记录数值模拟的褶积模型图49 Nw=32,v=3000,fm=25，h=1400地震记录数值模拟的褶积模型图410 Nw=32, h=1000,fm=25，v=2000地震记录数值模拟的褶积模型图411 Nw=32, h=1000,fm=25，v=4000地震记录数值模拟的褶积模型图412 Nw=32, h=1000,fm=25，v=6000地震记录数值模拟的褶积模型子波振幅谱：图51 fm=100的子波振幅谱图52 fm=100的子波振幅谱2、对比分析a) 由图31、图32、图33（或图41、图42、图43）可知，当子波长度（Nw=32）、深度（h=1000）、速度(v=3000)不变，子波频率变化时，褶积模型不变，且均在0.5到1s之间；b) 由图34、图35、图36可知，当深度（h=1000）、速度（v=3000）、子波频率（fm=25）不变, 子波长度变化，且垂直入射时，褶积模型为直线模型，是因为垂直入射时，时间与深度为线性关系；c) 由图44、图45、图46可知，当深度（h=1000）、速度（v=3000）、子波频率（fm=25）不变, 子波长度变化，且非垂直入射时，褶积模型前半部分为双曲线模型，后半部分为直线模型，是由于计算的褶积结果的实际长度小于所取的长度，计算机赋的随机数所致，所以，非垂直入射时的褶积模型只有图件上显示的前半部分的双曲线，是由于非垂直入射时，时间与深度为双曲线关系；d) 由图37、图38、图39，当子波长度（Nw=32）、速度（v=3000）、子波频率（fm=25）不变,深度变大时，垂直入射时，图件上显示的褶积模型的位置逐渐向下移，即地震波的旅行路程变大，旅行时变长；e) 由图47、图48、图49可知，当子波长度（Nw=32）、速度（v=3000）、子波频率（fm=25）不变,深度变大时，非垂直入射时，图件上显示的褶积模型的位置不变，这是由于非垂直入射时，相当于地下有一半圆形界面，在圆心处自激自收；f) 由图310、图311、图312（或图410、图411、图412）可知，当子波长度（Nw=32）、深度（h=1000）、子波频率（fm=25）不变,速度变大时，图件上显示的褶积模型的位置逐渐向上移，这是由于随着速度变大，地震波的旅行时变小；g) 由图51与图52可知，当子波频率变化时，其振幅谱不变。七、讨论建议1、实验收获通过此次试验，初步掌握了褶积模型基本理论、实现方法与程序编制，由褶积模型初步分析地震信号的分辨率问题，掌握了褶积模型与子波主频、子波长度、界面深度、介质中地震波速度的关系。2、存在问题对褶积模型的理论实现过程不是十分清楚，对结果的物理含义理解不够深入。3、其他问题由于不同的计算机，计算精度不一样，计算得到的数据结果可能会有部分差异，但总体趋势不变。4、心得体会在此次试验中，应特别注意褶积模型的理论实现过程的理解，以及程序调试时，要特别仔细地去检查每一个错误，每修改一处错误，就重新运行一次程序。附程序代码：/=1. 预处理部分=/#include#include#includefloat Cnltn(float,float);float Rflct(float,float,float);float Wave(float,float);#define Nx 128#define Nt 256#define Nw 32#define PI 3.1415926/=2. 主程序波分=/void 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 function=/float Wave(float fm,float dt)FILE *fpw; int It; float WaNw,t; if(fpw=fopen(wave.dat,wb)=NULL)printf(Connot open file wave); for(It=0;ItNw;It+)t=It*dt;WaIt=cos(2*PI*fm*t)*exp(-2*PI*PI*fm*fm*t);/形成子波fwrite(&WaIt,sizeof(WaIt),1,fpw); fclose(fpw); return(1);/ =3.2 Reflect Formaing function=/float Rflct(float dt,float h,float v)FILE *fpr; int It,Ix,J,Ltdpth; float t,dx=20,x; float ReNt;printf(请输入J：n);scanf(%d,&J); if(fpr=fopen(Reflect.dat,wb)=NULL)printf(Connot open file Reflect); for(Ix=0;IxNx;Ix+)for(It=0;ItNt;It+)ReIt=0.;if(J=1)t=2*h/v;/垂直入射反射界面由深度转换为自激自收时间if(J=2)x=Ix*dx;t=2*sqrt(h*h+x*x)/v;/非垂直入射反射界面由深度转换为自激自收时间if(J!=1)&(J!=2)printf(输入错误n); Ltdpth=(int)(t/dt); ReLtdpth=1; for(It=0;ItNt;It+)fwrite(&ReIt,sizeof(ReIt),1,fpr); fclose(fpr); return(1);/ =3.3 Convolution function=/float Cnltn(float dt,float dx)FILE *fpc,*fpw,*fpr;int It,Ix,Itao; float Wa1Nw,WaNw,ReNt+Nw+Nw,Re1Nt,t; float ConNt+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;Ix2;Ix+)for(It=0;ItNw;It+)fread(&Wa1It,sizeof(Wa1It),1,fpw);for(It=0;ItNw;It+)WaIt=Wa1Nw-It-1;/褶积前子波准备fclose(fpw);for(Ix=0;IxNx;Ix+)for(It=0;ItNt;It+)fread(&Re1It,sizeof(&Re1It