地震子波波形显示及一维地震合成记录

地震勘探原理实验一

地震子波波形显示及一维地震合成记录

姓名： 学号：专业：地球物理勘察技术与工程 级 一、实验目的

1. 认识子波，对子波的波形有直观的认识。（名词：零相位子波，混合相位

子波，最小相位子波；了解子波的分辨率与频宽的关系；） 2. 利用褶积公式合成一维地震记录。 二、实验步骤 1. 雷克子波

()()

))(21(22

t f e

t r m t f m ππ-=- 零相位子波

())2sin()

ln(222

t f e

t w m n t f m π-= （最小相位子波） n= m1/m2为最大波

峰m1和最大波谷m2之比

()())2cos(log *22xw t f e

t w m m

t f m +=-π 钟型子波 xw 为初相m 为时间域主波峰与次波峰之比

w(t)=exp(-2*Fm^2*t^2*ln(n))*sin(T-2*pi*Fm*t) n=m1/m2 最大相位子波

（最大相位子波请同学们自己查找相关文献完成，非必须完成）

其中

f代表子波的中心频率， t =i*dt，dt为时间采样间隔，i为时间

m

离散点序号; 这里可以为

f = 10，25，40，100 Hz等，采样间隔dt=0.002

m

秒，i为0～256;

2.根据公式编程实现不同频率的零相位子波的波形显示；

不同中心频率的零相位子波图

f = 25:

m

f = 100:

m

3.其地质模型为：

设计反射系数)(n r (n=512)，n 为地层深度，其中0.1)100(=r ，为第一层介质深度；7.0)200(-=r ，为第二层介质深度；5.0)300(=r ，为第三层介质深度；4.0)400(=r ，为第四层介质深度；6.0)450(=r ，为第五层介质深度；其它为0。

地震波在介质中传播，当到达介质分界面时，发生反射和透射，反射波被检波器接受，生成地震记录。反射系数表示地震波在两层介质分界面的能量重新分配，如r(100)=1.0，表示地震波入射到分界面时，只有一种波，反射纵波（或反射横波）。反射系数不为1.0时，表示当地震波入射到分界面时，产生两种反射波。反射系数为正，表示反射波相位与入射波相位相差2π；反射系数为负，表示反射波相位与入射波相位相差π，存在半波损失。

4. 应用褶积公式∑=-=*=N

m m n w m r n w n r n f 1)()()()()(合成一维地震记录，并图

形显示；

应用褶积公式求f （n ）的程序为：

#include #include #define PI 3.1415926 #define FM 100 void main() {

double fac(double x[],double y[],double z[],int m,int n);

FILE *fp;

int i,j,x;

double W,dt=0.002,t,a[256];

double b[512]={0};

double r[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("please input x:\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);

else if (x==2)

W=(1-2*pow(PI*FM*t,2))*exp(-pow(PI*FM*t,2)); else if (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]);

}

}

double fac(double x[],double y[],double z[],int m,int n)

{

int i,j;

for(i=0;i<=m+n-1;i++)

{

double sum=0.0;

for(j=0;j<=m;j++)

{

if(i-j>0&&i-j<=256)

sum+=x[j]*y[i-j];

}

z[i]=sum;

}

}

三、实验结果

一维反射系数序列图形显示为：

零相位子波与反射系数褶积后的地震记录图形显示：f = 25:

m

f = 100:

m

最小相位子波与反射系数褶积后的地震记录图形显示：

f = 25:

m

f = 100:

m

混合相位子波与反射系数褶积后的地震记录图形显示：

f = 25:

m

f = 100:

m

最大相位子波与反射系数褶积后的地震记录图形显示：

f = 25:

m

f = 100:

m

零相位振幅图形显示：f = 25:

m

f = 100:

m

零相位幅角图形显示：f = 25:

m

f = 100:

m

最小相位振幅图形显示：f = 25:

m

f = 100:

m

最小相位幅角图形显示：f = 25:

m

f = 100:

m

混合相位幅角图形显示：

f = 25:

m

f = 100:

m

混合相位振幅图形显示：

f = 25:

m

f = 100:

m

最大相位幅角图形显示：

f = 25:

m

f = 100:

m

最大相位振幅图形显示：f = 25:

m

f = 100:

m

四、实验分析

根据所学知识对实验结果进行分析；

地震子波由震源激发，在地层中传播，因为在沉积地层中，每层介质的物理性质不相同，从而使得地震波的传播速度也不相同。当地震波传播到两层介质的分界面时，会发生反射，由于每层介质的反射系数不同，所以反射波的能量也不相同，检波器接收到不同时刻的、不同能量的反射波，形成一个地震记录。

由合成地震记录中可以看出，最小相位子波相对零相位子波来说是相位滞后的，能量延迟的，但两者为同一家族的子波。合成地震记录中横坐标为时间，纵坐标为振幅。每一时刻的值由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

#define PI 3.1415926

#define FM 100

void main()

{

double fac(double x[],double y[],double z[],int m,int n);

FILE *fp,*fpr,*fpre,*fpi,*fpamp,*fpha;

int i,j,x;

double

W,dt=0.002,t,a[256],pr[512],pi[512]={0.0},fr[512],fi[512],amp[512],p ha[512];

double b[512]={0};

double r[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("please input x:\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));

else if (x==2)

W=exp(-2*FM*FM*t*t*log(2))*sin(2*PI*FM*t);

else if (x==3)

W=exp(-2*pow(FM*t,2)*log(2))*sin(0.512-2*PI*FM*t);

else if(x==4)

W=exp(-FM*FM*t*t*log(2))*cos(2*PI*FM*t+PI/4); a[i]=W;

}