地震记录数值模拟的褶积模型法实验报告.docx
《地震记录数值模拟的褶积模型法实验报告.docx》由会员分享，可在线阅读，更多相关《地震记录数值模拟的褶积模型法实验报告.docx（40页珍藏版）》请在冰豆网上搜索。
地震记录数值模拟的褶积模型法实验报告
地震记录数值模拟的褶积模型法
实验报告
专业：
勘查技术以工程
学号：
XXXXXXXXXXXXXXX
姓名：
XXXXXXXXXXXXX
实验题目：
地震记录数值模拟的褶积模型法实验报告
一、实验目的
掌握褶积模型基本理论;
实现方法与程序编制;
由褶积模型初步分析地震信号的分辨率问题。
二、实验原理和方法：
1、褶积原理
地震勘探的震源往往是带宽很宽的脉冲，在地下传播、反射、绕射到测线，传播经过中高频衰减，能量被吸收。
吸收过程可以看成滤波的过程，滤波可以用褶积完成。
在滤波中，反射系数与震源强弱关联，吸收作用与子波关联。
最简单的地震记录数值模拟，可以看成反射系数与子波的褶积。
通常，反射系数是脉冲，子波取雷克子波。
（1）雷克子波:
（2）反射系数:
（3）褶积公式:
数值模拟地震记录trace（t）:
trace（t）=rflct（t）*wave（t）
反射系数的参数由z变成了t，怎么实现？
在简单水平层介质，分垂直和非垂直入射两种实现，分别如图1和图2所示。
1）垂直入射：
t=2*h/v
2）非垂直入射：
2、褶积方法
（1）离散化（数值化）
计算机数值模拟要求首先必须针对连续信号离散化处理。
反射系数在空间模型中存在，不同深度反射系数不同，是深度的函数。
子波是在时间记录上一延续定时间的信号，是时间的概念。
在离散化时，通过深度采样完成反射系数的离散化，通过时间采样完成子波的离散化。
如果记录是Trace（t），则记录是时间的函数，以时间采样离散化。
时间采样间距以t表示，深度采样间距以z表示。
在做多道的数值模拟时，还有横向x的概念，横向采样间隔以x表示。
离散化的实现：
或：
（2）离散序列的褶积
三、实验内容
1、垂直入射地震记录数值模拟的褶积模型；
2、非垂直入射地震记录数值模拟的褶积模型；
3、点绕射的地震记录数值模拟的褶积模型；
三、实验步骤
（一）实验一：
1、根据垂直入射褶积模型理论算法，填充程序（附后）的下划线部分，使程序完整，调试程序，算出结果，用“Fimage”显示软件显示褶积结果；
2、根据非零偏移距算法，编制非零偏移距褶积模型程序，算出结果，用“Fimage”显示软件显示褶积结果。
（参考垂直入射褶积模型理论算法和程序，子波与反射层不变）；
3、变换子波的长度Nw和主频fm，重复1和2；
（1）Nw=80的情况下：
1）fm=5，2）fm=10，3）fm=15，3）fm=20，3）fm=55；
（2）Nw=160的情况下：
1）fm=5，2）fm=10，3）fm=15，3）fm=20，3）fm=55；
（二）实验二：
两层反射界面
如图3所示为两层反射界面的模型，对比不同间距H12情况下，不同长子波长度Nw和主频fm的褶积模型。
1、H12=20m时：
（1）Nw=80的情况下：
1）fm=5，2）fm=10，3）fm=15，4）fm=20，5）fm=55；
（2）Nw=160的情况下：
1）fm=5，2）fm=10，3）fm=15，4）fm=20，5）fm=55；
2、H12=40m时：
（1）Nw=80的情况下：
1）fm=5，2）fm=10，3）fm=15，3）fm=20，3）fm=55；
（2）Nw=160的情况下：
1）fm=5，2）fm=10，3）fm=15，3）fm=20，3）fm=55；
3、H12=60m时：
（1）Nw=80的情况下：
1）fm=5，2）fm=10，3）fm=15，3）fm=20，3）fm=55；
（2）Nw=160的情况下：
1）fm=5，2）fm=10，3）fm=15，3）fm=20，3）fm=55；
4、H12=80m时：
（1）Nw=80的情况下：
1）fm=5，2）fm=10，3）fm=15，3）fm=20，3）fm=55；
（2）Nw=160的情况下：
1）fm=5，2）fm=10，3）fm=15，3）fm=20，3）fm=55；
（三）点绕射的褶积模型
（1）地质模型与参数：
如图4所示为点绕射模型（下图）及其地震记录（上图）。
参数为：
Nx=128，Nz=200，V=3500m/s,h=1000m。
如图4图下图所示，地下绕射点位置为：
（Ix0=NX/2,Iz0=1000/20）。
其他参数与前述同。
（2）走时计算
如图4所示，绕射点正上方检波器的自激自收世间为t0：
t0=2h/v
任意位置检波器的自激自收世间为ti：
（3）褶积模型
爆炸反射界面传播到检波点的脉冲值游侠是计算：
（四）观察并记下实验结果图
四、实验结果分析
1、结果显示
（1）源程序段：
实验一.垂直入射地震记录数值模拟的褶积模型
单层反射界面
/////预处理部分/////
#include
#include
#include
intCnltn（float,float）;
intRflct（float,float,float）;
intWave（float,float）;
#defineNx128
#defineNt256
#defineNw160
/////主程序波分/////
voidmain（）
{
floatdt=0.004,dx=20,fm=20,h=1000,v=3000;
intiflag_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"）;
}
/////函数实现部分/////
//WaveFormaingfunction///
intWave（floatfm,floatdt）
{
FILE*fpw;
intIt;
floatWa[Nw],t;
doublepai=3.1415926;
if（（fpw=fopen（"wave.dat","wb"））==NULL）printf（"Connotopenfile""wave"""）;
for（It=0;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）;
}
/////ReflectFormaingfunction/////
intRflct（floatdt,floath,floatv）
{
FILE*fpr;
intIt,Ix,Ltdpth;
floatt;
floatRe[Nt];
if（（fpr=fopen（"Reflect.dat","wb"））==NULL）printf（"Connotopenfile""Reflect"""）;
for（Ix=0;Ix{
for（It=0;It{
Re[It]=0;
t=（float）2*h/v;
Ltdpth=（int）（t/dt）;
Re[Ltdpth]=1;
}
for（It=0;It{
fwrite（&Re[It],sizeof（Re[It]）,1,fpr）;
}
}
fclose（fpr）;
return
（1）;
}
/////Convolutionfunction/////
intCnltn（floatdt,floatdx）
{
FILE*fpc,*fpw,*fpr;
intIt,Ix,Itao;
floatWa1[Nw],Wa[Nw],Re[Nt+Nw+Nw],Re1[Nt];
floatCon[Nt+Nw];
if（（fpc=fopen（"Convosion.dat","wb"））==NULL）printf（"Connotopenfile""Convosion"""）;
if（（fpw=fopen（"wave.dat","rb"））==NULL）printf（"Connotopenfile""wave"""）;
if（（fpr=fopen（"Reflect.dat","rb"））==NULL）printf（"Connotopenfile""Reflect"""）;
for（Ix=1;Ix<2;Ix++）
{
for（It=0;It{
fread（&Wa1[It],sizeof（Wa1[It]）,1,fpw）;
/////子波数据反褶/////
Wa[Nw-1-It]=Wa1[It];
}
}
fclose（fpw）;
for（Ix=0;Ix{
for（It=0;It{
fread（&Re1[It],sizeof（&Re1[It]）,1,fpr）;
}
for（It=0;It{
Re[It]=0.;
}
for（It=0;It{
Re[It+Nw]=Re1[It];/////反射系数数据移动:
"0到Nt"移为"Nw到Nw+Nt"/////
}
for（It=0;It{
Con[It]=0;
for（Itao=0;Itao{
Con[It]=Con[It]+Wa[Itao]*Re[It+Itao];///褶积运算///
}
}
for（It=Nw/2;It{
fwrite（&Con[It],sizeof（Con[It]）,1,fpc）;
}
}
fclose（fpw）;
fclose（fpr）;
fclose（fpc）;
return
（1）;
}
双层反射界面
/////预处理部分/////
#include
#include
#include
intCnltn（float,float）;
intRflct（float,float,float,float）;
intWave（float,float）;
#defineNx128
#defineNt256
#defineNw160
/////主程序波分/////
voidmain（）
{
floatdt=0.004,dx=20,fm=15,h=1000,dh=80,v=4000;
intiflag_Co,iflag_Re,iflag_Wv;
if（iflag_Wv=Wave（fm,dt）!
=1）printf（"Waveiserror"）;
if（iflag_Re=Rflct（dt,h,dh,v）!
=1）printf（"Reflectioniserror"）;
if（iflag_Co=Cnltn（dt,dx）!
=1）printf（"Convosioniserror"）;
}
/////函数实现部分/////
/////WaveFormaingfunction/////
intWave（floatfm,floatdt）
{
FILE*fpw;
intIt;
floatWa[Nw],t;
doublepai=3.1415926;
if（（fpw=fopen（"wave.dat","wb"））==NULL）printf（"Connotopenfile""wave"""）;
for（It=0;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）;
}
/////ReflectFormaingfunction/////
intRflct（floatdt,floath,floatdh,floatv）
{
FILE*fpr;
intIt,Ix,Ltdpth,Ltdpth1;
floatt,t1;
floatRe[Nt];
if（（fpr=fopen（"Reflect.dat","wb"））==NULL）printf（"Connotopenfile""Reflect"""）;
for（Ix=0;Ix{
for（It=0;It{
Re[It]=0;
t=（float）2*h/v;
Ltdpth=（int）（t/dt）;
Re[Ltdpth]=1;
t1=（float）2*（h+dh）/v;
Ltdpth1=（int）（t1/dt）;
Re[Ltdpth1]=1;
}
for（It=0;It{
fwrite（&Re[It],sizeof（Re[It]）,1,fpr）;
}
}
fclose（fpr）;
return
（1）;
}
/////Convolutionfunction/////
intCnltn（floatdt,floatdx）
{
FILE*fpc,*fpw,*fpr;
intIt,Ix,Itao;
floatWa1[Nw],Wa[Nw],Re[Nt+Nw+Nw],Re1[Nt];
floatCon[Nt+Nw];
if（（fpc=fopen（"Convosion.dat","wb"））==NULL）printf（"Connotopenfile""Convosion"""）;
if（（fpw=fopen（"wave.dat","rb"））==NULL）printf（"Connotopenfile""wave"""）;
if（（fpr=fopen（"Reflect.dat","rb"））==NULL）printf（"Connotopenfile""Reflect"""）;
for（Ix=1;Ix<2;Ix++）
{
for（It=0;It{
fread（&Wa1[It],sizeof（Wa1[It]）,1,fpw）;
///子波数据反褶///
Wa[Nw-1-It]=Wa1[It];
}
}
fclose（fpw）;
for（Ix=0;Ix{
for（It=0;It{
fread（&Re1[It],sizeof（&Re1[It]）,1,fpr）;
}
for（It=0;It{
Re[It]=0.;
}
for（It=0;It{
Re[It+Nw]=Re1[It];
/////反射系数数据移动:
"0到Nt"移为"Nw到Nw+Nt"/////
}
for（It=0;It{
Con[It]=0;
for（Itao=0;Itao{
Con[It]=Con[It]+Wa[Itao]*Re[It+Itao];/////褶积运算/////
}
}
for（It=Nw/2;It{
fwrite（&Con[It],sizeof（Con[It]）,1,fpc）;
}
}
fclose（fpw）;
fclose（fpr）;
fclose（fpc）;
return
（1）;
}
实验二.非垂直入射地震记录数值模拟的褶积模型
/////预处理部分/////
#include
#include
#include
intCnltn（float,float）;
intRflct（float,float,float,float）;
intWave（float,float）;
#defineNx128
#defineNt256
#defineNw60
/////主程序波分/////
voidmain（）
{
floatdt=0.004,dx=20,fm=5,h=1000,v=3000;
intiflag_Co,iflag_Re,iflag_Wv;
if（iflag_Wv=Wave（fm,dt）!
=1）printf（"Waveiserror"）;
if（iflag_Re=Rflct（dt,h,v,dx）!
=1）printf（"Reflectioniserror"）;
if（iflag_Co=Cnltn（dt,dx）!
=1）printf（"Convosioniserror"）;
}
/////3.函数实现部分/////
/////3.1WaveFormaingfunction/////
intWave（floatfm,floatdt）
{
FILE*fpw;
intIt;
floatWa[Nw],t;
doublepai=3.1415926;
if（（fpw=fopen（"wave.dat","wb"））==NULL）printf（"Connotopenfile""wave"""）;
for（It=0;It{
t=（float）（It+1）*dt;
Wa[It]=（float）（1.-2*pow（（pai*fm*t）,2））*（float）exp（-2*t*pow（（pai*fm）,2））;/////雷克子波的计算/////
fwrite（&Wa[It],sizeof（Wa[It]）,1,fpw）;
}
fclose（fpw）;
return
（1）;
}
/////3.2ReflectFormaingfunction/////
intRflct（floatdt,floath,floatv,floatdx）
{
FILE*fpr;
intIt,Ix,Ltdpth;
floatt;
floatRe[Nt];
if（（fpr=fopen（"Reflect.dat","wb"））==NULL）printf（"Connotopenfile""Reflect"""）;
for（Ix=0;Ix{
for（It=0;It{
Re[It]=0.;
}
t=（float）sqrt（4*h*h+Ix*Ix*dx*dx）/v;
Ltdpth=（int）（t/dt）;
Re[Ltdpth]=1.;
for（It=0;It{
fwrite（&Re[It],sizeof（Re[It]）,1,fpr）;
}
}
fclose（fpr）;
return
（1）;
}
/////3.2Convolutionfunction/////
intCnltn（floatdt,floatdx）
{
FILE*fpc,*fpw,*fpr;
intIt,Ix,Itao;
floatWa1[Nw],Wa[Nw],Re[Nt+Nw+Nw],Re1[Nt];
floatCon[Nt+Nw];
if（（fpc=fopen（"Convosion.dat","wb"））==NULL）printf（"Connotopenfile""Convosion"""）;
if（（fpw=fopen（"wave.dat","rb"））==NULL）printf（"Connotopenfile""wave"""）;
if（（fpr=fopen（"Reflect.dat","rb"））==NULL）printf（"Connotopenfile""Reflect"""）;
for（Ix=1;Ix<2;Ix++）
{
for（It=0;It{
fread（&Wa1[It],sizeof（Wa1[It]）,1,fpw）;
Wa[Nw-1-It]=Wa1[It];/////子波数据反褶/////
}
}
fclose（fpw）;
for（Ix=0;Ix{
for（It=0;It{
fread（&Re1[It],sizeof（&Re1[It]）,1,fpr）;
}
for（It=0;It{
Re[It]=0.;
}
for（It=0;It{
Re[Nw+It]=Re1[It];/////反射系数数据移动:
"0到Nt"移为"Nw到Nw+Nt"/////
}
for（It=0;It{
Con[It]=0;
for（Itao=0;Itao{
Con[It]=Con[It]+Wa[Itao]*Re[It+Itao];/////褶积运算/////
}
}
for（It=Nw/2;It{
fwrite（&Con[It],sizeof（Con[It]）,1,fpc）;
}
}
fclose（fpw）;
fclose（fpr）;
fclose（fpc）;
return
（1）;
}
实验三.点绕射的褶积模型
/////预处理部分/////
#include
#include
#include
intCnltn（float,float）;
intRflct（float,float,float,float）;
intWave（float,float）
升级会员 