地震数值模拟实验报告Word下载.docx
- 文档编号:5819421
- 上传时间:2023-05-05
- 格式:DOCX
- 页数:32
- 大小:213.82KB
地震数值模拟实验报告Word下载.docx
《地震数值模拟实验报告Word下载.docx》由会员分享,可在线阅读,更多相关《地震数值模拟实验报告Word下载.docx(32页珍藏版)》请在冰点文库上搜索。
利用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<
stdio.h>
math.h>
stdlib.h>
conio.h>
string.h>
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"
}
if(Po2Judge(Nx)!
Nx=%disnotthePowerof2\n"
Nx);
Nx=%disthePowerof2\n"
if(Absorb(Labs)!
Absorbiserror\n"
);
exit(0);
Absorbisokay\n"
if(Rflct()!
Reflectioniserror\n"
Reflectionisokay\n"
};
if(Vlcty()!
Vlctyiserror\n"
Vlctyisokay\n"
if(WvFld0()!
WvFldiserror\n"
WvFldisokay\n"
if(PsFrwd()!
PsFurdiserror\n"
PsFurdisokay\n"
remove("
Wfld0r.dat"
Wfld0i.dat"
PsFrwrdMdl.ncb"
PsFrwrdMdl.dsp"
PsFrwrdMdl"
Abs"
//return();
//////////////////////////////////////////////////
intPo2Judge(intN)//JudgeNtorNxisPowerof2
doublea=0;
for(inti=0;
a<
N;
i++)
{
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"
Absbisokay\n"
if((fp_Abs=fopen("
Abs.txt"
w"
Absisokay\n"
for(Ix=0;
Ix<
Nx;
Ix++)
Abs[Ix]=1.;
Labs;
Abs[Ix]=sqrt(sin(PI*Ix/(2*(Labs-1))));
Abs[Nx-Ix-1]=Abs[Ix];
fwrite(&
Abs[Ix],sizeof(Abs[Ix]),Nx,fp_Absorb);
//printf("
%f\n"
Abs[Ix]);
fprintf(fp_Abs,"
fclose(fp_Absorb);
fclose(fp_Abs);
return
(1);
intRflct()//FormingReflectStructureModel
FILE*fp_Rflct;
intIx,Iz;
floatRflct[Nz];
if((fp_Rflct=fopen("
Rflct.dat"
Cannotopenfile"
Reflection"
for(Iz=0;
Iz<
Nz;
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);
intVlcty()//FormingVelociyModel
FILE*fp_Vlcty;
floatVlcty[Nz];
if((fp_Vlcty=fopen("
Vlcty.dat"
))==NULL)printf("
Vlcty"
for(Iz=0;
2*Nz/3;
{
Vlcty[Iz]=3000.;
//fwrite(&
Vlcty[Iz],sizeof(Vlcty[Iz]),1,fp_Vlcty);
}
for(Iz=(2*Nz/3);
Vlcty[Iz]=6000.;
Vlcty[0],sizeof(Vlcty[Iz]),Nz,fp_Vlcty);
fclose(fp_Vlcty);
intWvFld0()//WaveFieldInitialization
FILE*fp_Rflct,*fp_Wfld0r,*fp_Wfld0i,*fp_Wfld0isee,*fp_Wfld0rsee,*fp_Wfld0r0see;
intIx,Iz,It;
floatWfld0r[Nt],Wfld0i[Nt];
if((fp_Wfld0r=fopen("
CannotopenWfld0r.dat"
if((fp_Wfld0i=fopen("
CannotopenWfld0i.dat"
rb"
CannotopenRflct.dat"
//if((fp_Wfld0isee=fopen("
fp_Wfld0isee.txt"
Cannotopenfp_Wfld0isee.dat"
//if((fp_Wfld0rsee=fopen("
fp_Wfld0rsee.txt"
Cannotopenfp_Wfld0rsee.dat"
//if((fp_Wfld0r0see=fopen("
fp_Wfld0r0see.txt"
Cannotopenfp_Wfld0r0see.dat"
Wavefield0_FFT:
Ix=%d\n"
Ix);
{
fread(&
for(It=0;
It<
Nt;
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)!
FFTiserror"
}//原为时间域函数,变后为频率域。
(int)(Nt/2)+1;
fwrite(&
Wfld0r[It],sizeof(Wfld0r[It]),1,fp_Wfld0r);
Wfld0i[It],sizeof(Wfld0i[It]),1,fp_Wfld0i);
//fprintf(fp_Wfld0isee,"
Wfld0i[It]);
//fprintf(fp_Wfld0rsee,"
Wfld0r[It]);
//fprintf(fp_Wfld0isee,"
\n"
//fprintf(fp_Wfld0rsee,"
//fprintf(fp_Wfld0r0see,"
fclose(fp_Wfld0r);
fclose(fp_Wfld0i);
intPsFrwd()//WaveFieldPhaseShift
intPhaseShift();
intFrqcy2Time();
if(PhaseShift()!
=1){printf("
PhaseShiftiserror\n"
if(Frqcy2Time()!
Frqcy2Timeiserror\n"
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;
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速度与削波数据读入
))==NULL)printf("
for(Iz=0;
fread(&
if((fp_Absb=fopen("
CannotopenAbsb.dat"
Absb[Ix],sizeof(Absb[Ix]),1,fp_Absb);
fclose(fp_Absb);
//remove("
//3.OpenInitialWaveFieldFileandCurrentWaveFieldFileusingInWaveFieldExtrapolating波场文件打开
if((fp_Wfldr=fopen("
Wfldr.dat"
CannotopenWfldr.dat"
if((fp_Wfldi=fopen("
Wfldi.dat"
CannotopenWfldi.dat"
//4.WaveFiedExtrapolateWithEveryFrequency每个频率的波场延拓
for(Iw=0;
Iw<
Nw/2+1;
Iw++)
PhaseShift:
Iw=%d\n"
Iw);
//4.1InitializingWavefieldofEveryFrequency:
Allget0初始化当前波场
for(Ix=0;
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;
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);
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文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 地震 数值 模拟 实验 报告
![提示](https://static.bingdoc.com/images/bang_tan.gif)