地震子波波形显示及一维地震合成记录.docx
《地震子波波形显示及一维地震合成记录.docx》由会员分享，可在线阅读，更多相关《地震子波波形显示及一维地震合成记录.docx（22页珍藏版）》请在冰点文库上搜索。
地震子波波形显示及一维地震合成记录
地震勘探原理实验一
地震子波波形显示及一维地震合成记录
姓名：
_学号：
专业：
地球物理勘察技术与工程级
一、实验目的
1.认识子波，对子波的波形有直观的认识。
（名词：
零相位子波，混合相位
子波，最小相位子波；了解子波的分辨率与频宽的关系；）
2.利用褶积公式合成一维地震记录。
二、实验步骤
1.
雷克子波
22
wte^mtln（n）sin（2：
f』）（最小相位子波）n=m1/m2为最大波峰ml和最大波谷m2之比
wt）=e一fm厂*"mcos（2nfmt+xw）钟型子波xw为初相m为
最大相位子波
时间域主波峰与次波峰之比
w（t）=exp（-2*FmA2*tA2*ln（n））*sin（T-2*pi*Fm*t））n=m1/m2
E10-6^1宣克子施
（最大相位子波请同学们自己查找相关文献完成，非必须完成）
其中fm代表子波的中心频率，t=i*dt,dt为时间采样间隔，i为时间离散点序号；这里可以为fm=10,25,40,100Hz等，采样间隔dt=0.002秒，i为0〜256;
2.根据公式编程实现不同频率的零相位子波的波形显示;
不同中心频率的零相位子波图
fm=25：
fm=100：
3.其地质模型为:
设计反射系数r（n）（n=512）,n为地层深度，其中r（100）=1.0,为第一层介质深度；r（200）=-0.7，为第二层介质深度；r（300）=0.5，为第三层介质深度；r（400）=0.4，为第四层介质深度；r（450）=0.6，为第五层介质深度；其它为0。
量重新分配，如r（100）=1.0，表示地震波入射到分界面时，只有一种波，反射纵波（或反射横波）。
反射系数不为1.0时，表示当地震波入射到分
界面时，产生两种反射波。
反射系数为正，表示反射波相位与入射波相位相差2n;反射系数为负，表示反射波相位与入射波相位相差n，存在半波
损失。
N
4.应用褶积公式f（n）=r（n）”w（n）—r（m）w（n-m）合成一维地震记录，并图md
形显示；
应用褶积公式求f（n）的程序为：
#include
#include
#definePI3.1415926
#defineFM100
voidmain（）
{
doublefac（doublex[],doubley[],doublez[],intm,intn）;
FILE*fp;
inti,j,x;
doubleW,dt=0.002,t,a[256];
doubleb[512]={0};
doubler[512]={0};
r[100]=1.0;
r[200]=-0.7;
r[300]=0.5;
r[400]=0.4;
r[450]=0.6;fp=fopen（"Date.txt","w+"）;printf（"pleaseinputx:
\n"）;scanf（"%d",&x）;
for（i=0;i<256;i++）
{
t=i*dt;
if（x==1）W=exp（-2*FM*FM*t*t*log
（2））*sin（2*PI*FM*t）;
elseif（x==2）
W=（1-2*pow（PI*FM*t,2））*exp（-pow（PI*FM*t,2））;elseif（x==3）
W=exp（-FM*FM*t*t*log
（2））*cos（2*PI*FM*t+PI/4）;a[i]=W;
}
fac（r,a,b,512,256）;for（j=0;j<512;k++）
{
fprintf（fp,"%f\n",b[j]）;
}
}
doublefac（doublex[],doubley[],doublez[],intm,intn）{
inti,j;
for（i=0;i<=m+n-1;i++）
{
doublesum=0.0;
for（j=0;j<=m;j++）
{
if（i-j>0&&i-j<=256）sum+=x[j]*y[i-j];
}z[i]=sum;
}
}
三、实验结果一维反射系数序列图形显示为：
Q.S
04
-04
-05
ttoiuo•・gigi
nszQrjloi
scNiriT^oriinri.oriLnir'O
41iiiiipi
零相位子波与反射系数褶积后的地震记录图形显示:
fm=25：
fm=100:
0.4
02
-0.4
-0.6
最小相位子波与反射系数褶积后的地震记录图形显示:
fm=25：
系列1
混合相位子波与反射系数褶积后的地震记录图形显示:
fm=25：
fm=100:
-1
最大相位子波与反射系数褶积后的地震记录图形显示:
fm=25：
O.b
-0.6
系列1
零相位振幅图形显示:
fm=25：
—枭列1
fm=100：
400E+61
3.50E+（J4
3OOE0
2.S0Etfi4
200E+64
150E+61
系列L
LOJt+M
5OQE465
□.OOE-rfJJ
nffitnihps.mtnh-fflminr+zmGiirTZEbunri、阿孑训戸〜
■n寸SODRICN寸9尽6口“*—旳理EPM屮CJEEREEF寸F寸寸寸〒寸
零相位幅角图形显示:
fm=25：
ICO
100
最小相位振幅图形显示:
fm=25：
最小相位幅角图形显示:
fm=25：
100
83
43
23
D
23
-43
-33
1OT
III1IIIU1IIIIII1IIIIII
mrnrnnnrmiirriiiiriiiririiriiriiniriirii
103
ao
60
4D
2D
-SO
-S3
-100
混合相位幅角图形显示:
fm=25：
1CQ
-ICO
ICO
混合相位振幅图形显示:
fm=25：
fm=100:
最大相位幅角图形显示:
乘列二
fm=25：
100
-ioa
100
最大相位振幅图形显示:
fm=25：
斗3E+£4
3.WEt64
250E+64
1.50E1-64
1.00E+64
j.OOE4-e3
□JME4-00
四、实验分析
根据所学知识对实验结果进行分析；
地震子波由震源激发，在地层中传播，因为在沉积地层中，每层介质的物理性质不相同，从而使得地震波的传播速度也不相同。
当地震波传播到两层介质的分界面时，会发生反射，由于每层介质的反射系数不同，所以反射波的能量也不相同，检波器接收到不同时刻的、不同能量的反射波，形成一个地震记录。
由合成地震记录中可以看出，最小相位子波相对零相位子波来说是相位滞后的，能量延迟的，但两者为同一家族的子波。
合成地震记录中横坐标为时间，纵坐标为振幅。
每一时刻的值由m个值的和组成，m为反射系数r（n）的长度，
整个地震记录由m+n-1个时刻的值组成。
对于零相位的地震记录来说，当
r（m）=1.0时，即j=100时，i=100时，w（i-j）=1.0，是能量最大的，即
w（0）=1.0。
同理，当n=200,、300、400、450时，能量也是最大的。
对于最小相位的地震记录来说，当r（m）=1.0时，即j=100时，但i=100时，w（i-j）不是最大能量的，即最大能量不是在w（0）出现，而是延迟出现。
同理，当n=200
300、400、450时，能量也不是最大的，而是要延迟出现。
由振幅图及幅角图可知，零相位子波能量聚集在首部，开始时就具有最大能量，无积累过程，当振幅最大时，相位为零，即此时波的振幅为实数，达到最大值；最小相位子波能量聚集在序列首部，是最小能量延迟的，信号随时间增大而减小，当振幅最大时，相位不为零，是非零相位的，相对零相位子波来说，最大能量是延迟的；混合相位子波的能量聚集在序列中部，是混合能量延迟的；最大相位子波能量聚集在后部。
最大相位子波和混合相位子波的信号信号不随时间增大而减小。
五、附：
源程序代码
#include#include"13KFFT.C"
#include
#definePI3.1415926
#defineFM100
voidmain（）
{
doublefac（doublex[],doubley[],doublez[],intm,intn）;
FILE*fp,*fpr,*fpre,*fpi,*fpamp,*fpha;
inti,j,x;
double
W,dt=0.002,t,a[256],pr[512],pi[512]={0.0},fr[512],fi[512],amp[512],pha[512];
doubleb[512]={0};
doubler[512]={0};
r[100]=1.0;
r[200]=-0.7;
r[300]=0.5;
r[400]=0.4;
r[450]=0.6;
fp=fopen（"褶积结果.txt","w+"）;
fpr=fopen（"反射系数.csv","w+"）;
fpre=fopen（"实部.txt","w+"）;
fpi=fopen（"虚部.txt","w+"）;
fpamp=fopen（"振幅.csv","w+"）;
fpha=fopen（"相位.csv","w+"）;
for（i=0;i<512;i++）
{
fprintf（fpr,"%f\n",r[i]）;
}
fclose（fpr）;
printf（"pleaseinputx:
\n"）;
scanf（"%d",&x）;
for（i=0;i<256;i++）
{
t=i*dt;
if（x==1）
W=（1-2*pow（PI*FM*t,2））*exp（-pow（PI*FM*t,2））;
elseif（x==2）
W=exp（-2*FM*FM*t*t*log
（2））*sin（2*PI*FM*t）;
elseif（x==3）
W=exp（-2*pow（FM*t,2）*log
（2））*sin（0.512-2*PI*FM*t）;
elseif（x==4）
W=exp（-FM*FM*t*t*log
（2））*cos（2*PI*FM*t+PI/4）;
a[i]=W;
fac（r,a,b,512,256）;
for（j=0;j<512;j++）
{
fprintf（fp,"%f\n",b[j]）;
}
for（i=0;i<512;i++）
{
pr[i]=b[j];
}
for（i=0;i<512;i++）
{
pr[i]=fr[i+127];
}
kfft（pr,pi,512,9,fr,fi,0,1）;
for（i=0;i<512;i++）
{
fprintf（fpre,"%e\n",fr[i]）;
fprintf（fpi,"%e\n",fi[i]）;
fprintf（fpamp,"%e\n",pr[i]）;
fprintf（fpha,"%f\n",pi[i]）;
fclose（fpre）;
fclose（fpi）;
fclose（fpamp）;
fclose（fpha）;
}
doublefac（doublex[],doubley[],doublez[],intm,intn）{
inti,j;
for（i=0;i<=m+n-1;i++）
{
doublesum=0.0;
for（j=0;j<=m;j++）
{
if（i-j>0&&i-j<=256）
sum+=x[j]*y[i-j];
}
z[i]=sum;
升级会员 