地震记录数值模拟的褶积模型法
实验课程数值模型模拟
学院名称地球物理学院
专业名称勘测技术与工程
学生姓名
学生学号
指导教师熊高君
实验地点 5417
实验成绩
2015年5月
成都理工大学
地震数值模拟》实验报告
实验报告
实验题目:
地震记录数值模拟的褶积模型法
二、实验目的:
掌握褶积模型基本理论、实现方法与程序编制,由褶积模型初步分析地震信号的分辨率问题。
三、原理公式
1、褶积原理
地震勘探的震源往往是带宽很宽的脉冲,在地下传播、反射、绕射到测线,传播经过中高频衰减,能量被吸收。吸收过程可以看成滤波的过程,滤波可以用褶积完成。在滤波中,反射系数与震源强弱关联,吸收作用与子波关联。最简单的地震记录数值模拟,可以看成反射系数与子波的褶积。通常,反射系数是脉冲,子波取雷克子波。
(1)雷克子波:
wave(t)=cos(2 ft)*
(2)反射系数:
(3)褶积公式:
数值模拟地震记录trace(t): trace(t) =rflct(t)*wave(t);
反射系数的参数由 z 变成了 t,怎么实现?在简单水平层介质,分垂直和非垂直
图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;
或:I t=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—1 Nw=32,h=1000,v=3000,fm=100地震记录数值模拟的褶积模型(左)和子波(右)
图3—2 Nw=32,h=1000,v=3000,fm=200地震记录数值模拟的褶积模型(左)和子波(右)
图3—3 Nw=32,h=1000,v=3000,fm=300地震记录数值模拟的褶积模型(左)和子波(右)
图3—4 h=1000,v=3000,fm=25,Nw=20地震记录数值模拟的褶积模型(左)和子波(右)
图3—5 h=1000,v=3000,fm=25,Nw=30地震记录数值模拟的褶积模型(左)和子波(右)
图3—7 Nw=32,v=3000,fm=25,h=1000地震记录数值模拟的褶积模型
图3—10 Nw=32, h=1000,fm=25,v=2000
地震记录数值模拟的褶积模型
2)非垂直入射
图4—1 Nw=32,h=1000,v=3000,fm=100地震记录数值模拟的褶积模型(左)和子波(右)
图4—4 h=1000,v=3000,fm=25,Nw=20地震记录数值模拟的褶积模型(左)和子波(右)
图4—7 Nw=32,v=3000,fm=25,h=1000地震记录数值模拟的褶积模型
图4—10 Nw=32, h=1000,fm=25,v=2000
地震记录数值模拟的褶积模型
子波振幅谱:
图5—1 fm=100的子波振幅谱
图5—2 fm=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、心得体会
在此次试验中,应特别注意褶积模型的理论实现过程的理解,以及程序调试时,要特别仔细地去检查每一个错误,每修改一处错误,就重新运行一次程序。
附程序代码:
#include
#include
float Rflct(float,float,float);
float Wave(float,float);
#define Nx 128
#define Nt 256
#define Nw 32
#define PI 3.1415926
//==================2. 主程序波分============== void main()
{
float dt=0.004,dx=20,fm=25,h=1000,v=3000; int
iflag_Co,iflag_Re,iflag_Wv; if(iflag_Wv=Wave(fm,dt)!=1) printf("Wave is error"); if(iflag_Re=Rflct(dt,h,v)!=1)
printf("Reflection is error"); if(iflag_Co=Cnltn(dt,dx)!=1)
printf("Convosion is error");
}
// =================3.函数实现部分=========
// ==============3.1 Wave Formaing function float Wave(float fm,float dt)
FILE *fpw;
int It;
float Wa[Nw],t; if((fpw=fopen("wave.dat","wb"))==NULL)
printf("Connot open file ""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.2 Reflect Formaing function==============// float Rflct(float dt,float h,float v) { FILE *fpr; int It,Ix,J,Ltdpth; float t,dx=20,x; float Re[Nt]; printf(" 请输入 J : \n"); scanf("%d",&J); if((fpr=fopen("Reflect.dat","wb"))==NULL) printf("Connot open file ""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); for(It=0;It Re[Ltdpth]=1; fwrite(&Re[It],sizeof(Re[It]),1,fpr); } fclose(fpr); return(1); } // =============3.3 Convolution function============= float Cnltn(float dt,float dx) { FILE *fpc,*fpw,*fpr; int It,Ix,Itao; float Wa1[Nw],Wa[Nw],Re[Nt+Nw+Nw],Re1[Nt],t; float Con[Nt+Nw]; if((fpc=fopen("Convosion.dat","wb"))==NULL) printf("Connot open file ""Convosion"""); if((fpw =fopen("wave.dat","rb"))==NULL) printf("Connot open file ""wave"""); if((fpr =fopen("Reflect.dat","rb"))==NULL) printf("Connot open file ""Reflect"""); for(Ix=1;Ix<2;Ix++) { for(It=0;It for(It=0;It { fread(&Wa1[It],sizeof(Wa1[It]),1,fpw); for(It=0;It Wa[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];// 反射系数准备 }