当前位置:文档之家› 数值模拟实验报告

数值模拟实验报告

一、实验题目地震记录数值模拟的这几模型法 二、实验目的掌握褶积模型基本理论、实现方法与程序编制,由褶积模型初步分析地震 信号的分辨率问题 三、实验原理 1、褶积原理地震勘探的震源往往是带宽很宽的脉冲,在地下传播、反射、绕射到测线,传播经过中 高频衰减,能量被吸收。

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

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

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

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

(1) 雷克子波Wave(t) = (1−2π2f 2t 2)e −2π2f 2t 2(2) 反射系数:rflct(z)={1 z =z 反射界面0 z =others(3) 褶积公式:数值模拟地震记录trace(t): trace(t) =rflct(t)*wave(t)反射系数的参数由z 变成了t ,怎么实现?在简单水平层介质,分垂直和非垂直入射两 种实现,分别如图1 和图2 所示。

1) 垂直入射: t =2hv图一 垂直入射 2) 非垂直入射:t =2√h 2+x 2v图二非垂直入射 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) 离散序列的褶积 trace (It )=∑rflct(Itao) ×wave(It −Itao)∞Itao=−∞ 四、实验内容1、 垂直入射地震记录数值模拟的褶积模型;2、 非垂直入射地震记录数值模拟的褶积模型;3、点绕射的地震记录数值模拟的褶积模型;五、方法路线根据褶积模型的实验原理编写C++程序,完成对于垂直入射波的褶积。

改变子波的长度与主频的大小,关注其对于实验结果的影响。

通过增加一个地层来模拟地下两层界面的反射情况,通过改变界面的高度来说明其对于实验结果的影响,同时改变子波长度与主频。

非垂直入射改变时间t来改变褶积结果显示地面情况。

点绕射模型通过时间的改变,任意位置检波器的自激自搜时间来改变褶积结果,同时改变子波长度与主频来分析影响。

六、实验结果:1.垂直入射:单层界面1>h=1000,fm=10,Nw=80v=2000v=3250 v=4500 v=5750v=70002>h=1000,v=4000,fm=10Nw=80Nw=100Nw=120Nw=140Nw=1603>h=1000,v=4000,Nw=80fm=5fm=10fm=15fm=20fm-464>h=1000,v=4000,Nw=160fm=5fm=10fm=15fm=20fm=465>v=4000,fm=10,Nw=80h=800h=1000h=1200 h=1400h=1600双层界面1>H12=20h=1000,v=4000,Nw=80fm=5fm=10fm=15fm=20fm=46h=1000,v=4000,Nw=160fm=5fm=10fm=15fm=20fm=462>H12=40h=1000,v=4000,Nw=80fm=5fm=10fm=15fm=20fm=46h=1000,v=4000,Nw=160fm=5fm=10fm=15fm=20fm=463>H12=60h=1000,v=4000,Nw=80fm=5fm=10fm=15fm=20fm=46h=1000,v=4000,Nw=160 fm=5fm=10fm=15fm=20fm=46H12=80h=1000,v=4000,Nw=80 fm=5fm=10fm=15fm=20fm=46h=1000,v=4000,Nw=160fm=5fm=10fm=15fm=20fm=462.非垂直入射1> h=1000,fm=10,Nw=80v=3000v=4000 v=5000 v=6000 v=7000 2> h=1000,v=4000,fm=10 Nw=80Nw=100Nw=120Nw=140Nw=1604>h=1000,v=4000,Nw=80fm=5fm=10fm=15fm=20fm=55 h=1000,v=4000,Nw=160 fm=5fm=10fm=15fm=20fm=55v=4000,fm=10,Nw=80 h=800h=1000 h=1200 h=1400 h=16003.点绕射卷积模型图件Nw=80 fm=15相关程序:#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 32void main(){floatdt=0.004,dx=20.0,fm=25.0,h=1000.0,v=3000.0;int iflag_Co,iflag_Re,iflag_Wv;if(iflag_Wv=Wave(fm,dt)!=1)printf("Waveiserror"); if(iflag_Re=Rflct(dt,h,v)!=1)printf("Reflectioniserror ");if(iflag_Co=Cnltn(dt,dx)!=1)printf("Convosioniserror ");}//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("Connotopenfile""wave""");for(It=0;It<Nw;It++){t=(float)(It+1)*dt;Wa[It]=(float)(1.-2*pai*pai*fm*fm*t*t)*(float) exp(-2*pai*pai*fm*fm*t); //雷克子波的计算?fwrite(&Wa[It],sizeof(Wa[It]),1,fpw);}fclose(fpw);return(1);}//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("C onnot 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);}//Convolution functionint Cnltn(float dt,float dx){FILE*fpc,*fpw,*fpr;int It,Ix,Itao;floatWa1[Nw],Wa[Nw],Re[Nt+Nw+Nw],Re1[Nt];float Con[Nt+Nw];if((fpc=fopen("Convosion.dat","wb"))==NULL) printf("Connotopenfile""Convosion""");if((fpw=fopen("wave.dat","rb"))==NULL)printf("Con notopenfile""wave""");if((fpr=fopen("Reflect.dat","rb"))==NULL)printf("Co nnotopenfile""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[It]=Re1[It+Nw];//反射系数数据移动:"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[-Itao+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);}对比分析:1.垂直入射的结果是一条直线,反映了地下为水平界面;非垂直入射的结果是半条的双曲线,其反映了地下的非水平界面。

相关主题