地震数值模拟实验报告.docx
- 文档编号:9985995
- 上传时间:2023-02-07
- 格式:DOCX
- 页数:32
- 大小:213.82KB
地震数值模拟实验报告.docx
《地震数值模拟实验报告.docx》由会员分享,可在线阅读,更多相关《地震数值模拟实验报告.docx(32页珍藏版)》请在冰豆网上搜索。
地震数值模拟实验报告
本科生实验报告
实验课程数值模型模拟
学院名称地球物理学院
专业名称勘查技术与工程
学生姓名ZRY
学生学号
指导教师
实验地点624
实验成绩
二〇一五年4月二〇一五年5月
成都理工大学
《地震数值模拟》实验报告
实验时间
2015年5月31日
开课单位
地球物理学院
指导教师
实验题目:
叠加地震记录的相移波动方程正演模拟实验
姓名
学号
班级
专业
勘查技术与工程
院(系)
地球物理学院
单
项
成
绩
内容理解
写作结构
程序设计
模型设计
计算结果
结果分析
总成绩
实验二叠加地震记录的相移波动模拟实方程正演验
摘要
利用C语言编制地质模型的相移波动方程正演模拟,改变绕射点位置、速度,再做正演模拟。
关键字:
地震模型;正演记录
1.1实验目的
掌握各向同性介质任意构造、水平层状速度结构地质模型的相移波动方程正演模拟基本理论、实现方法与程序编制,由正演记录初步分析地震信号的分辨率。
1.2实验内容
1、基本要求:
(1)点绕射构造和水平层状速度模型(参数如图1所示)的正演数值模拟;
1)削波的正演;
2)无削波的震正演;
(2)计算中点和两个边界的信号位置,分析实验结果的正确性;
(3)做同样模型的褶积模型数值模拟,对比分析分析两者的异同。
(4)改变绕射点位置、速度,再做正演模拟。
2、较高要求:
(1)使用雷克子波做爆炸源,对三个不同的主频:
25hz、50hz和75hz分别做点绕射
模型的正演模拟;
(2)设计复杂反射构造模型,再做正演模拟。
1.3实验原理
1、地震波传播的波动方程
设(x,z)为空间坐标,t为时间,地震波传播速度为v(x,z),则二位介质中任意位置、任意时刻的地震波场为p(z,x,t):
压缩波——纵波。
则二维各向同性均匀介质中地震波传播的遵循声波方程为
2、傅里叶变换的微分性质
p(t)与其傅里叶变换的P(ω)的关系:
则有时间微分性质
ω为频率,ω=2π/T,T为周期。
同理有空间微分性质:
k为频率,k=2π/λ,λ为波长。
3、地震波传播的相移外推公式
令速度v不随x变化,只随z变化,则利用傅里叶变换微分性质(3)和(4)式,把波
动方程
(1)式变换到频率-波数域,得:
或:
令:
则(5)式的解为:
包括上行波和下行波两项。
正演模拟取上行波:
若
和
间隔为△Z,速度v(z)在此间隔内不随Z变的常数,(7)式实现波场从
到
的延拓,即:
在深度Zj+1开始向上延拓到Zj,若延拓深度为零,即:
∆Z=Zj+1-Zj=0,则
对于任意深度Zj+1到Zj的延拓,可得正演模拟中地震波的传播方程(延拓公式)
4、初始条件和边界条件
按照爆炸界面理论,反射界面震源在t=0时刻同时起爆,此时刻的波场就是震源。
根据不同情况,可直接使用反射系数脉冲或子波作震源。
如果直接使用反射系数作震源脉冲,则初始条件可表示为:
对时间t和空间x做二维傅立叶变换,则得频率-波数域的初始波场
边界条件:
其他参数都是在
范围内定义的。
5、边界处理
(1)边界反射问题
把实际无穷空间区域中求解波场的问题化为有穷区域求解时,左右两边使用零边界条件。
物理上假设探区距
两个端点很远,在两个端点上收到的反射波很弱。
但是,上述条件在实际中不能成立,造成零边界条件反而成为绝对阻止波通过的强反射面。
在正演模拟的剖面上出现了边界假反射干涉正常界面的反射。
(2)边界强反射的处理
镶边法、削波法、吸收边界都能有效消除边界强反射。
削波法就是在波场延拓过程中,没延拓一次,在其两侧均匀衰减到零,从而消除边界强反射的影响。
假设横向总长度为NX,以两边Lx道吸波为例,有以下吸波公式:
6、数字化
根据数字信号处理的采样定理,把连续的信号变为计算机能处理的数字信号,使相移法正演模拟得以实现。
频域抽样定理:
一个频谱受限信号f(t),如果时间只占据-tm—+tm的范围,若在频域以不大于1/2tm频率间隔∆f≤1/2tm对信号f(t)的频谱F(ω)采样,则抽样到的离散信号F1(ω)可以唯一表示原信号。
时域抽样定理:
一个时间受限信号f(t),如果频谱只占据-ωm—+ωm的范围,则信号f(t)可以用等间隔的抽样值唯一表示出来,而时间∆t抽样间隔必须不大于1/2fm,ωm=2πfm,∆t≤1/2fm。
1.4实验过程
//#include"PsFrwrdMdlParameter.h"
#include
#include
#include
#include
#include
intkkfft(floatpr[],floatpi[],intn,intl);
intAbsorb(intLabs);
intPo2Judge(int);
intRflct();
intVlcty();
intWvFld0();
intPsFrwd();
intexp_ikzDz(floateikzdz[],intIx,floatVc,intIw,floatDw,floatDkx);
intPhaseShift();//相移
intFrqcy2Time();
#defineNx128//TraceNumber
#defineNt256//RecordNumber
#defineNz100//DepthNumber
#defineDx20.//TraceInterval
#defineDt0.004//RecordInterval
#defineDz20.//DepthInterval
#definePI3.1415
voidmain()
{
intLabs=10;//LengthOfBoundaryAbsorbing
if(Po2Judge(Nt)!
=1){printf("Nt=%disnotthePowerof2\n",Nt);exit(0);}
else{printf("Nt=%disthePowerof2\n",Nt);}
if(Po2Judge(Nx)!
=1){printf("Nx=%disnotthePowerof2\n",Nx);exit(0);}
else{printf("Nx=%disthePowerof2\n",Nx);}
if(Absorb(Labs)!
=1){printf("Absorbiserror\n");exit(0);}
else{printf("Absorbisokay\n");}
if(Rflct()!
=1){printf("Reflectioniserror\n");exit(0);}
else{printf("Reflectionisokay\n");};
if(Vlcty()!
=1){printf("Vlctyiserror\n");exit(0);}
else{printf("Vlctyisokay\n");}
if(WvFld0()!
=1){printf("WvFldiserror\n");exit(0);}
else{printf("WvFldisokay\n");}
if(PsFrwd()!
=1){printf("PsFurdiserror\n");exit(0);}
else{printf("PsFurdisokay\n");}
remove("Wfld0r.dat");
remove("Wfld0i.dat");
remove("PsFrwrdMdl.ncb");
remove("PsFrwrdMdl.dsp");
remove("PsFrwrdMdl");
remove("Abs");
//return();
}
//////////////////////////////////////////////////
intPo2Judge(intN)//JudgeNtorNxisPowerof2
{
doublea=0;
for(inti=0;a { a=pow(2,i); } if(fabs(N-a)>=1)return(0); elsereturn (1); } ////////////////////////////////////////////////// intAbsorb(intLabs)//FormingAbsorbingBoudary { FILE*fp_Absorb,*fp_Abs; intIx; floatAbs[Nx]; if((fp_Absorb=fopen("Absb.dat","wb"))==NULL) printf("Cannotopenfile""Absorb"""); else{printf("Absbisokay\n");} if((fp_Abs=fopen("Abs.txt","w"))==NULL) printf("Cannotopenfile""Abs"""); else{printf("Absisokay\n");} for(Ix=0;Ix { Abs[Ix]=1.; } for(Ix=0;Ix { Abs[Ix]=sqrt(sin(PI*Ix/(2*(Labs-1)))); Abs[Nx-Ix-1]=Abs[Ix]; } for(Ix=0;Ix { fwrite(&Abs[Ix],sizeof(Abs[Ix]),Nx,fp_Absorb); //printf("%f\n",Abs[Ix]); fprintf(fp_Abs,"%f\n",Abs[Ix]); } fclose(fp_Absorb); fclose(fp_Abs); return (1); } ////////////////////////////////////////////////// intRflct()//FormingReflectStructureModel { FILE*fp_Rflct; intIx,Iz; floatRflct[Nz]; if((fp_Rflct=fopen("Rflct.dat","wb"))==NULL) printf("Cannotopenfile""Reflection"""); for(Ix=0;Ix { for(Iz=0;Iz { Rflct[Iz]=0.0; if(Ix==(int)(Nx/2)&&Iz==(int)(Nz/2))Rflct[Iz]=1.; fwrite(&Rflct[Iz],sizeof(Rflct[Iz]),1,fp_Rflct); } } fclose(fp_Rflct); return (1); } ////////////////////////////////////////////////// intVlcty()//FormingVelociyModel { FILE*fp_Vlcty; intIx,Iz; floatVlcty[Nz]; if((fp_Vlcty=fopen("Vlcty.dat","wb"))==NULL)printf("Cannotopenfile""Vlcty"""); for(Ix=0;Ix { for(Iz=0;Iz<2*Nz/3;Iz++) { Vlcty[Iz]=3000.; //fwrite(&Vlcty[Iz],sizeof(Vlcty[Iz]),1,fp_Vlcty); } for(Iz=(2*Nz/3);Iz { Vlcty[Iz]=6000.; //fwrite(&Vlcty[Iz],sizeof(Vlcty[Iz]),1,fp_Vlcty); } fwrite(&Vlcty[0],sizeof(Vlcty[Iz]),Nz,fp_Vlcty); } fclose(fp_Vlcty); return (1); } ////////////////////////////////////////////////// intWvFld0()//WaveFieldInitialization { FILE*fp_Rflct,*fp_Wfld0r,*fp_Wfld0i,*fp_Wfld0isee,*fp_Wfld0rsee,*fp_Wfld0r0see; intIx,Iz,It; floatRflct[Nz]; floatWfld0r[Nt],Wfld0i[Nt]; if((fp_Wfld0r=fopen("Wfld0r.dat","wb"))==NULL)printf("CannotopenWfld0r.dat"); if((fp_Wfld0i=fopen("Wfld0i.dat","wb"))==NULL)printf("CannotopenWfld0i.dat"); if((fp_Rflct=fopen("Rflct.dat","rb"))==NULL)printf("CannotopenRflct.dat"); //if((fp_Wfld0isee=fopen("fp_Wfld0isee.txt","wb"))==NULL)printf("Cannotopenfp_Wfld0isee.dat"); //if((fp_Wfld0rsee=fopen("fp_Wfld0rsee.txt","wb"))==NULL)printf("Cannotopenfp_Wfld0rsee.dat"); //if((fp_Wfld0r0see=fopen("fp_Wfld0r0see.txt","wb"))==NULL)printf("Cannotopenfp_Wfld0r0see.dat"); for(Ix=0;Ix { //printf("Wavefield0_FFT: Ix=%d\n",Ix); for(Iz=0;Iz { fread(&Rflct[Iz],sizeof(Rflct[Iz]),1,fp_Rflct); for(It=0;It { Wfld0r[It]=0.; Wfld0i[It]=0.; if(It==0)Wfld0r[It]=Rflct[Iz]; } //fprintf(fp_Wfld0r0see,"%f\t",Wfld0r[Iz]); if(kkfft(Wfld0r,Wfld0i,Nt,0)! =1){printf("FFTiserror");exit(0);}//原为时间域函数,变后为频率域。 for(It=0;It<(int)(Nt/2)+1;It++) { fwrite(&Wfld0r[It],sizeof(Wfld0r[It]),1,fp_Wfld0r); fwrite(&Wfld0i[It],sizeof(Wfld0i[It]),1,fp_Wfld0i); //fprintf(fp_Wfld0isee,"%f\t",Wfld0i[It]); //fprintf(fp_Wfld0rsee,"%f\t",Wfld0r[It]); } } //fprintf(fp_Wfld0isee,"\n"); //fprintf(fp_Wfld0rsee,"\n"); //fprintf(fp_Wfld0r0see,"\n"); } fclose(fp_Rflct); fclose(fp_Wfld0r); fclose(fp_Wfld0i); return (1); } ////////////////////////////////////////////////// intPsFrwd()//WaveFieldPhaseShift { intPhaseShift(); intFrqcy2Time(); if(PhaseShift()! =1){printf("PhaseShiftiserror\n");exit(0);} if(Frqcy2Time()! =1){printf("Frqcy2Timeiserror\n");exit(0);} return (1); } ////////////////////////////////////////////////// intPhaseShift()//相移 { //1.Prepprocedure预处理 FILE*fp_Wfldr,*fp_Wfldi,*fp_Wfld0r,*fp_Wfld0i,*fp_Vlcty,*fp_Absb; floatkz[2]; intIx,Iz,Iw,Nw=Nt,Ikx,Nkx=Nx; longMgrtn; floatVlcty[Nz]; floatWfld_r,Wfld_i; floatAbsb[Nx],Wfld0r[Nx],Wfld0i[Nx],Wfldr[Nx],Wfldi[Nx]; floatKxmax,Dkx,Wmax,Dw; Wmax=PI/Dt; Dw=Wmax/(float)Nt; Kxmax=PI/Dx; Dkx=Kxmax/(float)Nx; //2.ReadinVelocityandAbsorbingBoundaryData速度与削波数据读入 if((fp_Vlcty=fopen("Vlcty.dat","rb"))==NULL)printf("CannotopenRflct.dat"); for(Iz=0;Iz { fread(&Vlcty[Iz],sizeof(Vlcty[Iz]),1,fp_Vlcty); } fclose(fp_Vlcty); if((fp_Absb=fopen("Absb.dat","rb"))==NULL)printf("CannotopenAbsb.dat"); for(Ix=0;Ix { fread(&Absb[Ix],sizeof(Absb[Ix]),1,fp_Absb); } fclose(fp_Absb); //remove("Absb.dat"); //3.OpenInitialWaveFieldFileandCurrentWaveFieldFileusingInWaveFieldExtrapolating波场文件打开 if((fp_Wfld0r=fopen("Wfld0r.dat","rb"))==NULL)printf("CannotopenWfld0r.dat"); if((fp_Wfld0i=fopen("Wfld0i.dat","rb"))==NULL)printf("CannotopenWfld0i.dat"); if((fp_Wfldr=fopen("Wfldr.dat","wb"))==NULL)printf("CannotopenWfldr.dat"); if((fp_Wfldi=fopen("Wfldi.dat","wb"))==NULL)printf("CannotopenWfldi.dat"); //4.WaveFiedExtrapolateWithEveryFrequency每个频率的波场延拓 for(Iw=0;Iw { printf("PhaseShift: Iw=%d\n",Iw); //4.1InitializingWavefieldofEveryFrequency: Allget0初始化当前波场 for(Ix=0;Ix { Wfldr[Ix]=0.; Wfldi[Ix]=0.; } //4.2WaveFiedExtrapolateFromZ=ZmaxtoZ=Zmin波场从Iz=Nz-1最深处开始,延拓到Iz=1测线深度 for(Iz=Nz-1;Iz>0;Iz--) { //4.2.1FormingNewWavefieldforExtrapolatingonSpaceDomain形成新波长 for(Ix=0;Ix { Mgrtn=(Ix*Nz+Iz)*(Nw/2+1)+Iw; //TakeoutInitialWaveFieldDataWitheveryDepth取出当前深度的初始波场 fseek(fp_Wfld0r,sizeof(Wfld0r[Ix])*Mgrtn,0); fread(&Wfld0r[Ix],sizeof(Wfld0r[Ix]),1,fp_Wfld0r); fseek(fp_Wfld0i,sizeof(Wfld0i[Ix])*Mgrtn,0); fread(&Wfld0i[Ix],sizeof(Wfld0i[Ix]),1,fp_Wfld0i); //NewWaveFiedEqualtotheInitialonPlusCurrentone新波场=初始波场+从下面延拓到此处的波场 Wfldr[Ix]=Wfldr[Ix]+Wfld0r[Ix]; Wfldi[Ix]=Wfldi[Ix]+Wfld0i[Ix]; //BoundaryabsorbingofNewWaveFied边界削波: 新波场=新波场×削波因子 Wfldr[Ix]=Wfldr[Ix]*Absb[Ix]; Wfldi[Ix]=Wfldi[Ix]*Absb[Ix]; } //4.2.2FFTofNewWavefieldFromSpacetoWaveNumber新波场FFT到波数
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 地震 数值 模拟 实验 报告