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

本科生实验报告实验课程数值模型模拟 ________________________________ 学院名称______________________________ 专业名称___________________________ 学生姓名__________________________________________________ 学生学号__________________________________________________ 指导教师熊高君______________________________________ 实验地点5417 ___________________________________ 实验成绩__________________________________________________2015年5月成都理工大学《地震数值模拟》实验报告实验报告实验题目:地震记录数值模拟的褶积模型法二、实验目的：掌握褶积模型基本理论、实现方法与程序编制，由褶积模型初步分析地震信号的分辨率问题。

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

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

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

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

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

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

图1 图21) 垂直入射：t=2h/v ;2) 非垂直入射：t=2、褶积方法(1)离散化(数值化)计算机数值模拟要求首先必须针对连续信号离散化处理。

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

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

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

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

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

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

离散化的实现：t=It X A t ; x=lx X A x; z=lz X A z;或：It=t/ △ t; Ix=x/ △x; lz=z/ A z(2)离散序列的褶积trace(lt)= 「今rflct (Itao) * waveflt - Itao i四、实验内容1、垂直入射地震记录数值模拟的褶积模型；2、非垂直入射地震记录数值模拟的褶积模型。

五、方法路线1、根据垂直入射褶积模型理论算法，填充程序（附后）的下划线部分，使程序完整，调试程序，算出结果，用“ Fimage”显示软件显示褶积结果;2、根据非零偏移距算法，编制非零偏移距褶积模型程序，算出结果，用“Fimage”显示软件显示褶积结果。

（参考垂直入射褶积模型理论算法和程序，子波与反射层不变）；3、变换子波的主频：fm（10hz到300hz范围），重复1和2；4、变换子波的长度：Nw（ 80m$到160ms范围），重复1 和2；5、改变反射层深度：h（800n到1600n范围），重复1 和2；6、改变介质速度：V （2000m/s 到7000m/s 范围），重复1 和2六、实验结果1、结果显示1）垂直入射31 ]IL I 质$1 111 I滋图3— 1 Nw=32,h=1000,v=3000,fm=100地震记录数值模拟的褶积模型（左）和子波（右）im i制图3— 2 Nw=32,h=1000,v=3000,fm=200地震记录数值模拟的褶积模型（左）和子波（右）国im izi口口O iu-图3— 3 Nw=32,h=1000,v=3000,fm=300地震记录数值模拟的褶积模型（左）和子波（右）图3— 4 h=1000,v=3000,fm=25 , Nw=20地震记录数值模拟的褶积模型（左）和子波（右）9 IQI图3—7 Nw=32,v=3000,fm=25 , h=1000地震记录数值模拟的褶积模型图3— 6 h=1000,v=3000,fm=25 , Nw=40地震记录数值模拟的褶积模型（左）和子波（右）5] ]01 曲5] ]0L 曲图3—8 Nw=32,v=3000,fm=25 , h=1200地震记录数值模拟的褶积模型图3—9 Nw=32,v=3000,fm=25 , h=1400地震记录数值模拟的褶积模型图3—10 Nw=32, h=1000,fm=25 , v=2000地震记录数值模拟的褶积模型5] ]0L 曲图3—11 Nw=32, h=1000,fm=25 , v=4000地震记录数值模拟的褶积模型图3—12 Nw=32, h=1000,fm=25 , v=6000地震记录数值模拟的褶积模型2）非垂直入射图4— 1 Nw=32,h=1000,v=3000,fm=100地震记录数值模拟的褶积模型（左）和子波（右）图4— 2 Nw=32,h=1000,v=3000,fm=200地震记录数值模拟的褶积模型（左）和子波（右）图4— 3 Nw=32,h=1000,v=3000,fm=300地震记录数值模拟的褶积模型（左）和子波（右）图4— 4 h=1000,v=3000,fm=25 , Nw=20地震记录数值模拟的褶积模型（左）和子波（右）*51 IM图4— 6 h=1000,v=3000,fm=25 , Nw=40地震记录数值模拟的褶积模型（左）和子波（右）图4—7 Nw=32,v=3000,fm=25 , h=1000地震记录数值模拟的褶积模型15】SOL LZ95] ]0L 曲图4—8 Nw=32,v=3000,fm=25 , h=1200地震记录数值模拟的褶积模型51 10L LZ3图4—9 Nw=32,v=3000,fm=25 , h=1400地震记录数值模拟的褶积模型图4—10 Nw=32, h=1000,fm=25 , v=2000地震记录数值模拟的褶积模型图4— 12 Nw=32, h=1000,fm=25 , v=6000地震记录数值模拟的褶积模型子波振幅谱:圏表区i--- 系列10.4 1 3 5 7 9 11 13 15 17 19 21 23 25 27 29 31图5— 1 fm=100的子波振幅谱图4— 11 Nw=32, h=1000,fm=25 , v=4000地震记录数值模拟的褶积模型 】叽 ■0.50 ]I>LOS图5—2 fm=100的子波振幅谱2、对比分析a）由图3—1、图3—2、图3—3 （或图4—1、图4—2、图4—3）可知，当子波长度（Nw=32、深度（h=1000）、速度（v=3000）不变，子波频率变化时，褶积模型不变，且均在0.5到1s之间；b）由图3—4、图3—5、图3—6可知，当深度（h=1000）、速度（v=3000）、子波频率（fm=25）不变,子波长度变化，且垂直入射时，褶积模型为直线模型，是因为垂直入射时，时间与深度为线性关系；c）由图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）不变，深度变大时，垂直入射时，图件上显示的褶积模型的位置逐渐向下移，即地震波的旅行路程变大，旅行时变长；e）由图4—7、图4—8、图4—9可知，当子波长度（Nw=32、速度（v=3000）、子波频率（fm=25）不变，深度变大时，非垂直入射时，图件上显示的褶积模型的位置不变，这是由于非垂直入射时，相当于地下有一半圆形界面，在圆心处自激自收；当子波长度（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> float Cnltn(float,float);float Rflct(float,float,float);float Wave(float,float);#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 Wa[Nw],t;if((fpw=fopen("wave.dat","wb"))==NULL) printf("Connot open file ""wave"""); for(It=0;It<Nw;It++){t=It*dt;Wa[It]=cos(2*PI*fm*t)*exp(-2*PI*PI*fm*fm*t);//fwrite(&Wa[It],sizeof(Wa[It]),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 Re[Nt];printf(" 请输入J：\n");形成子波//scanf("%d",&J);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.;}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);Re[Ltdpth]=1;fwrite(&Re[It],sizeof(Re[It]),1,fpr);}for(It=0;It<Nt;It++)fclose(fpr);return(1);}// =============3.3 Convolution function=============// float Cnltn(float dt,float dx){FILE *fpc,*fpw,*fpr;int It,Ix,Itao;float Wa1[Nw],Wa[Nw],Re[Nt+Nw+Nw],Re1[Nt],t;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<Nt;It++)fread(&Wa1[It],sizeof(Wa1[It]),1,fpw);for(It=0;It<Nw;It++)Wa[It]=Wa1[Nw-It-1];// 褶积前子波准备} 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[It+Nw]=Re1[It];// 反射系数准备}Con[It]=0;for(It=0;It<Nt+Nw;It++){t=0;for(Itao=0;Itao<Nw;Itao++){t+=Re[It+Itao]*Wa[Itao];//褶积运算}Con[It]=t;}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);}。