数值模拟实验报告

一、实验题目地震记录数值模拟的这几模型法二、实验目的学握褶枳模型基本理论、实现方法与程序编制，由褶积模型初步分析地震信号的分辨率问题三、实验原理1、褶枳原理地震勘探的震源往往是带宽很宽的脉冲，在地卜传播、反射、绕射到测线，传播经过中高频衰减，能量被吸收。吸收过程町以看成滤波的过程，滤波町以用褶枳完成。在滤波中，反射系数与震源强弱关联，吸收作用与子波关联。最简单的地震记录数值模拟，可以看成反射系数与子波的褶积。通常，反射系数是脉冲，子波取雷克子波。(1)雷克子波Wave(t)=(1—2n2f2t2)e_21t2f2t2(2)反射系数：回arflctfzlJ1z=z反射杯(0z=others(3)褶积公式：数值模拟地廉记录trace(t):trace)t)=rflct(t)*wave(t)反射系数的参数由z变成了t,怎么实现？在简单水平层介质，分垂直和非垂直入射两种实现，分别如图1和图2所示。1)垂直入射：2)非垂直入射:2ht=Y图一垂直入射2、褶枳方法2Vh22ht=Y图一垂直入射2、褶枳方法2Vh2+x2t=V图二非垂直入射离散化(数值化)计算机数值模拟要求首先必须针对连续信石离散化处理。反射系数在空间模型中存在，不同深度反射系数不同，是深度的函数。子波是在时间记录上一延续定时间的信号，是时间的概念。在离散化时，通过深度采样完成反射系数的离散化，通过时间釆样完成子波的离散化。如果记录是Trace©则记录是时间的函数，以时间采样离散化。时间采样间距以At表示，深度采样间距以Az表示。在做多道的数值模拟时，还有横向x的概念，横向采样间隔以Ax表示。离散化的实现:t=ltxAt；x=lxxAx：z=lzxAz或：lx=x/Ax;lz=z/4z离散序列的褶积0trace(It)=Xitao=-oorflct(Itao)xwave(It—Itao)四、实验内容1、垂直入射地震记录数值模拟的褶积模型：2、非垂直入射地震记录数值模拟的褶积模型：3、点绕対的地震记录数值模拟的褶积模型五、方法路线根据褶枳模型的实验原理编写C++程序，完成对于垂直入射波的褶积。改变子波的长度与主频的大小，关注其对于实验结果的影响。通过增加一个地层来模拟地I、两层界面的反射情况，通过改变界面的高度来说明其对丁实验结果的影响，同时改变子波长度与主频。非垂直入射改变时间t來改变褶积结果显示地面情况。点绕射模型通过时间的改变，任意位置检波器的自激自搜时间來改变褶积结果，同时改变子波长度与主频来分析影响。六、实验结果：1.垂直入射：单层界而1>h=1000tfm=10»Nw=80v=2000v=3250v=70002>h=1000,v=4000,fm=10Nw=80v=4500v=5750Nw=100Nw=120Nw=140fm=20liiL—土」蔥丄ailkfa1Lk^*d.kiLU^LL*uLil.ttl'Nw=160fm-463>h=1000,v=40O0.Nw=80nmi«4>h=1000,v=4000,Nw=160fm=5fm二20■注―和BiMMM比■■■fm=465>v=4000tfm=10,Nw=80h=800h=1000h=1200上JLLi—1>」11・=讣||U..—k^iuULh^i如4iLLafauJJu.wklk-L«LmU■“L—亠4』山—ba.Fl~~WH~~h-1400CK<T<一-h"—“丄-厶」一・！■“■!h^LLiL.X*.bLl^・」丄ksl■亠Uu亠「丄山4HLhA~bd丄hufaahfUL——】LI4A-h=1600XI双层界而1>H12=20h=1000,v二4000,Nw=80fm=5fm=10h=1000,v=4000>Nw=160fm=5h=1000>v=4000.Nw=160fm=5fm=10fm=10fm=15fm=20fm=5h=1000»v=4000.Nw=80Imfm=20fm=15fm=46fm=20fm~46fm~46h=1000,v=4000>Nw=160fm=5fm=10h=1000.v-4000>Nw=160fm=5fm=10fm=20fm=20fm=20fm=46H12=80h=1000,v=4000,Nw=80fm=5fm=462非垂直入射l>h=1000・fm=10.Nw=80Nw=160Illfm=5h=10O0,v=4000,Nw=160fm=5fm=10fm=10fm=15fm=15fm=20fm=20~T11fm=553点绕射卷积模型图件:::相关程序：#includevmath・h>#includevstdio・h>#includevstring・h>intCnltnffloatJIoat);intRflctlfloatJIoatJIoat);intWave(floatfloat);//defineNx128林defineNt256#defineNw32voidmain(){floatdt=0・004/dx=20・O/fm=25・O/h=1OO0・0严=3000・0intiflag_Cojflag_Rejflag_Wv;if(iflag_Wv=Wave(fm/dt)!=l)printf(,,Waveiserror")if(iflag_Re=Rflct(dt/h/V)!=l)printf(,,Reflectioniserror“；if(iflag_Co=Cnltn(dt/dx)!=l)printf(HConvosioniserror//WaveFormaingfunctionintWave(floatfmjloatdt){FILE*fpw;intlt;floatWafNwJ.t;doublepai=3.1415926;if((fpw=fopen(Mwave.datl7，wb"))==NULL)printf("Connotopenfile""wave""")；for(lt=0;lt<Nw;lt++){t=(float)(lt+l)*dt;Wa[lt]=(float)(l.-2*pai*pai*fm*fm*t*t)*(float)exp(-2*pai*pai*fm*fm*t);//雷克子波的计算？fwrite(&Wa[lt],sizeof(Wa[lt]),l,fpw);}fclose(fpw);retum(l);}//ReflectFormaingfunctionintRflctffloatdtjloathjloatv){FILEYpr;intlt,lx,Ltdpth;floatt;floatRe[Nt];if((fpr=fopen("Reflectdat“,”wb'))=NULL)printf(“Connotopenfile""Reflect""")；for(lx=0;lx<Nx;lx++){for(lt=0;lt<Nt;lt++)for(lx=l;lx<2;lx++){for(lt=0;lt<Nw;lt++){fread(&Wal[lt],sizeof(Wal[lt]),l,fpw);Wa[Nw-l-lt]=Wal[lt];}//子波数据反褶?}fclose(fpw);tor(lx=0;lx<Nx;lx++){for(lt=0;lt<Nt;lt++){fread(&Rel[It],sizeof(&Rel[lt])几fpr);}for(lt=0;lt<Nt+2*Nw;lt++){Re[lt]=0.;}for(lt=0;lt<Nt;lt++)Re[lt]=Rel[lt+Nw];〃反射系数数据移Re[lt]=0.;}t=(float)2.*h/v;Ltdpth=(int)(t/dt);Re[Ltdpth]=l.;for(It=0;lt<Nt;lt++){fwrite(&Re[lt],sizeof(Re[lt]),l/fpr)；}}fclose(fpr);return(l);}//ConvolutionfunctionintCnltnffloatdt/loatdx){FILE*fpc/*fpw/*fpr;intltjxjtao;floatWa2[NwLWa[NwlRe[Nt+Nw+Nw],Rel[Nt];floatCon[Nt+Nw];if((fpc=fopen(HConvosion.datl7，wb,,))==NULL)printf(,,Connotopenfile,,,,Convosion,"");if((fpw=fopen(,iwave.dat'7，rb,,))==NULL)printf(,,Connotopenfile""wave""");if((fpr=fopen("Reflectdat","rb"))==NULL)printf("Connotopenfile""Reflect“””);动:”0到Nt“移为“Nw到Nw+Nt"}for(11=0;lt<Nt+Nw;lt++){Con[lt]=0;for(ltao=0;ltao<Nw;ltao++){Con[lt]=Con[lt]+Wa[ltao]*Re[-ltao+lt];//褶枳运算}}for(It=Nw/2;lt<Nt+Nw/2;lt++)fwrite(&Con[lt],sizeof(Con[lt]),Xfpc);}}}fclose(fpw);fclose(fpr);fclose(fpc);return(l);对比分析：垂立入射的结果是一条直线，反映了地下为水平界面：非垂II入射的结果是半条的双曲线，其反映了地下的非水平界面。子波的频率越高其反射结果面积越大，子波长度越大其时间越短。在两层反射界面垂直距离的情况卜，其界而间距越人则反射信号越容易识別，子波长度越长，主频越大垂直距离越大则识别极限越小。褶积结果与图4匕半部分基本吻合，I大1为任意检波器的自激自收情况与计算的褶枳过