地震记录数值模拟的褶积模型法.docx
《地震记录数值模拟的褶积模型法.docx》由会员分享，可在线阅读，更多相关《地震记录数值模拟的褶积模型法.docx（24页珍藏版）》请在冰豆网上搜索。
地震记录数值模拟的褶积模型法
本科生实验报告
实验课程数值模型模拟
学院名称地球物理学院
专业名称勘测技术与工程
学生姓名
学生学号
指导教师熊高君
实验地点5417
实验成绩
2015年5月
成都理工大学
《地震数值模拟》实验报告
实验时间
2015年5月
开课单位
地球物理学院
指导教师
熊高君
实验题目：
地震记录数值模拟的褶积模型法
姓名
学号
班级
专业
勘测技术与工程（石油物探）
院（系）
地球物理学院地球探测与信息技术系
单项成绩
内容理解
写作结构
程序设计
模型设计
计算结果
结果分析
总成绩
实验报告
一、实验题目：
地震记录数值模拟的褶积模型法
二、实验目的：
掌握褶积模型基本理论、实现方法与程序编制，由褶积模型初步分析地震信号的分辨率问题。
三、原理公式
1、褶积原理
地震勘探的震源往往是带宽很宽的脉冲，在地下传播、反射、绕射到测线，传播经过中高频衰减，能量被吸收。
吸收过程可以看成滤波的过程，滤波可以用褶积完成。
在滤波中，反射系数与震源强弱关联，吸收作用与子波关联。
最简单的地震记录数值模拟，可以看成反射系数与子波的褶积。
通常，反射系数是脉冲，子波取雷克子波。
（1）雷克子波：
wave（t）=cos（2
ft）*
（2）反射系数：
（3）褶积公式：
数值模拟地震记录trace（t）:
trace（t）=rflct（t）*wave（t）；
反射系数的参数由z变成了t，怎么实现？
在简单水平层介质，分垂直和非垂直入射两种实现，分别如图1和图2所示。
图1图2
1）垂直入射：
t=2h/v；
2）非垂直入射：
t=
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）=
四、实验内容
1、垂直入射地震记录数值模拟的褶积模型；
2、非垂直入射地震记录数值模拟的褶积模型。
五、方法路线
1、根据垂直入射褶积模型理论算法，填充程序（附后）的下划线部分，使程序完整，调试程序，算出结果，用“Fimage”显示软件显示褶积结果；
2、根据非零偏移距算法，编制非零偏移距褶积模型程序，算出结果，用“Fimage”显示软件显示褶积结果。
（参考垂直入射褶积模型理论算法和程序，子波与反射层不变）；
3、变换子波的主频：
fm（10hz到300hz范围），重复1和2；
4、变换子波的长度：
Nw（80ms到160ms范围），重复1和2；
5、改变反射层深度：
h（800m到1600m范围），重复1和2；
6、改变介质速度：
v（2000m/s到7000m/s范围），重复1和2。
六、实验结果
1、结果显示
1）垂直入射
图3—1Nw=32,h=1000,v=3000,fm=100地震记录数值模拟的褶积模型（左）和子波（右）
图3—2Nw=32,h=1000,v=3000,fm=200地震记录数值模拟的褶积模型（左）和子波（右）
图3—3Nw=32,h=1000,v=3000,fm=300地震记录数值模拟的褶积模型（左）和子波（右）
图3—4h=1000,v=3000,fm=25，Nw=20地震记录数值模拟的褶积模型（左）和子波（右）
图3—5h=1000,v=3000,fm=25，Nw=30地震记录数值模拟的褶积模型（左）和子波（右）
图3—6h=1000,v=3000,fm=25，Nw=40地震记录数值模拟的褶积模型（左）和子波（右）
图3—7Nw=32,v=3000,fm=25，h=1000地震记录数值模拟的褶积模型
图3—8Nw=32,v=3000,fm=25，h=1200地震记录数值模拟的褶积模型
图3—9Nw=32,v=3000,fm=25，h=1400地震记录数值模拟的褶积模型
图3—10Nw=32,h=1000,fm=25，v=2000地震记录数值模拟的褶积模型
图3—11Nw=32,h=1000,fm=25，v=4000地震记录数值模拟的褶积模型
图3—12Nw=32,h=1000,fm=25，v=6000地震记录数值模拟的褶积模型
2）非垂直入射
图4—1Nw=32,h=1000,v=3000,fm=100地震记录数值模拟的褶积模型（左）和子波（右）
图4—2Nw=32,h=1000,v=3000,fm=200地震记录数值模拟的褶积模型（左）和子波（右）
图4—3Nw=32,h=1000,v=3000,fm=300地震记录数值模拟的褶积模型（左）和子波（右）
图4—4h=1000,v=3000,fm=25，Nw=20地震记录数值模拟的褶积模型（左）和子波（右）
图4—5h=1000,v=3000,fm=25，Nw=30地震记录数值模拟的褶积模型（左）和子波（右）
图4—6h=1000,v=3000,fm=25，Nw=40地震记录数值模拟的褶积模型（左）和子波（右）
图4—7Nw=32,v=3000,fm=25，h=1000地震记录数值模拟的褶积模型
图4—8Nw=32,v=3000,fm=25，h=1200地震记录数值模拟的褶积模型
图4—9Nw=32,v=3000,fm=25，h=1400地震记录数值模拟的褶积模型
图4—10Nw=32,h=1000,fm=25，v=2000地震记录数值模拟的褶积模型
图4—11Nw=32,h=1000,fm=25，v=4000地震记录数值模拟的褶积模型
图4—12Nw=32,h=1000,fm=25，v=6000地震记录数值模拟的褶积模型
子波振幅谱：
图5—1fm=100的子波振幅谱
图5—2fm=100的子波振幅谱
2、对比分析
a）由图3—1、图3—2、图3—3（或图4—1、图4—2、图4—3）可知，当子波长度（Nw=32）、深度（h=1000）、速度（v=3000）不变，子波频率变化时，褶积模型不变，且均在0.5到1s之间；
b）由图3—4、图3—5、图3—6可知，当深度（h=1000）、速度（v=3000）、子波频率（fm=25）不变,子波长度变化，且垂直入射时，褶积模型为直线模型，是因为垂直入射时，时间与深度为线性关系；
c）由图4—4、图4—5、图4—6可知，当深度（h=1000）、速度（v=3000）、子波频率（fm=25）不变,子波长度变化，且非垂直入射时，褶积模型前半部分为双曲线模型，后半部分为直线模型，是由于计算的褶积结果的实际长度小于所取的长度，计算机赋的随机数所致，所以，非垂直入射时的褶积模型只有图件上显示的前半部分的双曲线，是由于非垂直入射时，时间与深度为双曲线关系；
d）由图3—7、图3—8、图3—9，当子波长度（Nw=32）、速度（v=3000）、子波频率（fm=25）不变,深度变大时，垂直入射时，图件上显示的褶积模型的位置逐渐向下移，即地震波的旅行路程变大，旅行时变长；
e）由图4—7、图4—8、图4—9可知，当子波长度（Nw=32）、速度（v=3000）、子波频率（fm=25）不变,深度变大时，非垂直入射时，图件上显示的褶积模型的位置不变，这是由于非垂直入射时，相当于地下有一半圆形界面，在圆心处自激自收；
f）由图3—10、图3—11、图3—12（或图4—10、图4—11、图4—12）可知，当子波长度（Nw=32）、深度（h=1000）、子波频率（fm=25）不变,速度变大时，图件上显示的褶积模型的位置逐渐向上移，这是由于随着速度变大，地震波的旅行时变小；
g）由图5—1与图5—2可知，当子波频率变化时，其振幅谱不变。
七、讨论建议
1、实验收获
通过此次试验，初步掌握了褶积模型基本理论、实现方法与程序编制，由褶积模型初步分析地震信号的分辨率问题，掌握了褶积模型与子波主频、子波长度、界面深度、介质中地震波速度的关系。
2、存在问题
对褶积模型的理论实现过程不是十分清楚，对结果的物理含义理解不够深入。
3、其他问题
由于不同的计算机，计算精度不一样，计算得到的数据结果可能会有部分差异，但总体趋势不变。
4、心得体会
在此次试验中，应特别注意褶积模型的理论实现过程的理解，以及程序调试时，要特别仔细地去检查每一个错误，每修改一处错误，就重新运行一次程序。
附程序代码：
//=================1.预处理部分==============//
#include
#include
#include
floatCnltn（float,float）;
floatRflct（float,float,float）;
floatWave（float,float）;
#defineNx128
#defineNt256
#defineNw32
#definePI3.1415926
//==================2.主程序波分==============//
voidmain（）
{
floatdt=0.004,dx=20,fm=25,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"）;
}
//=================3.函数实现部分===============//
//==============3.1WaveFormaingfunction=============//
floatWave（floatfm,floatdt）
{
FILE*fpw;
intIt;
floatWa[Nw],t;
if（（fpw=fopen（"wave.dat","wb"））==NULL）
printf（"Connotopenfile""wave"""）;
for（It=0;It{
t=It*dt;
Wa[It]=cos（2*PI*fm*t）*exp（-2*PI*PI*fm*fm*t）;//形成子波
fwrite（&Wa[It],sizeof（Wa[It]）,1,fpw）;
}
fclose（fpw）;
return
（1）;
}
//============3.2ReflectFormaingfunction==============//
floatRflct（floatdt,floath,floatv）
{
FILE*fpr;
intIt,Ix,J,Ltdpth;
floatt,dx=20,x;
floatRe[Nt];
printf（"请输入J：
\n"）;
scanf（"%d",&J）;
if（（fpr=fopen（"Reflect.dat","wb"））==NULL）
printf（"Connotopenfile""Reflect"""）;
for（Ix=0;Ix{
for（It=0;It{
Re[It]=0.;
}
if（J==1）
t=2*h/v;//垂直入射反射界面由深度转换为自激自收时间
if（J==2）
{
x=Ix*dx;
t=2*sqrt（h*h+x*x）/v;//非垂直入射反射界面由深度转换为自激自收时间
}
if（（J!
=1）&&（J!
=2））
printf（"输入错误\n"）;
Ltdpth=（int）（t/dt）;
Re[Ltdpth]=1;
for（It=0;It{
fwrite（&Re[It],sizeof（Re[It]）,1,fpr）;
}
}
fclose（fpr）;
return
（1）;
}
//=============3.3Convolutionfunction=============//
floatCnltn（floatdt,floatdx）
{
FILE*fpc,*fpw,*fpr;
intIt,Ix,Itao;
floatWa1[Nw],Wa[Nw],Re[Nt+Nw+Nw],Re1[Nt],t;
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）;
}
for（It=0;ItWa[It]=Wa1[Nw-It-1];//褶积前子波准备
}
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];//反射系数准备
}
for（It=0;It{
Con[It]=0;
t=0;
for（Itao=0;Itao{
t+=Re[It+Itao]*Wa[Itao];//褶积运算
}
Con[It]=t;
}
for（It=Nw/2;It{
fwrite（&Con[It],sizeof（Con[It]）,1,fpc）;
}
}
fclose（fpw）;
fclose（fpr）;
fclose（fpc）;
return
（1）;
}
升级会员 