实验要求

实验一（四课时）： A 、地震信号频谱分析； 给定一个地震记录（随机或计算得出），编 程分析该地震记录的振幅和相位情况，绘制其振 幅谱，相位谱。 注：可用 C 语言编程实现或软件实现。振幅 谱，相位谱用 Excel 绘制。 B 、褶积合成地震反射记录 ； 选择适合自己的编程语言（C、C++等）或 软件合成地震反射记录。结果用 Excel 绘制。

实验报告

 一、实验目的

运用卷积方法合成地震记录，并对地震记录进行频谱分析，获取地震记录的振幅谱和相位谱，并进行解释。

二、实验原理

A 、地震信号频谱分析

获取地震记录的振幅谱和相位谱，进行频谱分析。

因为信号是复数，所以把实部和虚部分开来计算，并加以合成，便可以得到振幅谱和相位谱。

B 、褶积合成地震反射记录 ；

卷积方法可以得到地震反射记录。原始数据为地层的反射系数系列，通过把反射系数序列和雷克子波卷积既可以合成地震记录。

雷克子波的计算公式为。

反射系数序列为自己制定，本次实验取四层地层。共取200个点，给四个点赋值，其他均为零。

把得到的反射系数序列和计算出来的雷克子波进行卷积即可以合成地震记录。、

三、实验结果分析

结果分析：地震记录为反射系数序列和雷克子波卷积而成，应该保证反射系数序列比雷克子波序列足够的长，否者就会出现地震记录不圆滑等很多问题，而且雷克子波应该把频率选择高一些，这样周期才比较短，在一定长度的雷克子波序列比较完整。这样卷积的到的结果更加真实。

四、源代码

A、频谱分析

#include<stdio.h>

#include<stdlib.h>

#include<math.h>

#define gcwt "file2.dat"

#define gcwt2 "file1.dat"

float pinpu(float*x,float *Xr,float *Xi);

int main()

{

float Xr[200],Xi[200],x[200],X[200],xiangwei[200];int i,j,k;

FILE *fp,*fp1;

fp=fopen(gcwt,"r+");

fp1=fopen(gcwt2,"w+");

for(i=0;i<200;i++)

{

fscanf(fp,"%f",&x[i]);

}//读取数据

fclose(fp);

pinpu(x,Xr,Xi);//计算实部虚部

for (i=0;i<200;i++)

{

X[i]=sqrt(Xr[i]*Xr[i]+Xi[i]*Xi[i]);//振幅谱

xiangwei[i]=atan(Xi[i]/Xr[i]);//相位谱

printf("%d\t%f\t%f\n",i,X[i],xiangwei[i]);

fprintf(fp1,"%d\t%f\t%f\n",i,X[i],xiangwei[i]);

}

fclose(fp1);

}

float pinpu(float*x,float *Xr,float *Xi)

{

int i,j;

for(i=0;i<200;i++)

{

Xr[i]=0.0;

Xi[i]=0.0;

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

{

Xr[i]=x[j]*cos(2*3.14/200*j*i)+Xr[i];//实部

Xi[i]=-x[j]*sin(2*3.14/200*j*i)+Xi[i];//虚部

}

}

return 0;

}

B、地震记录合成

#include<stdio.h>

#include<math.h>

#include<stdlib.h>

#define dz1 "file1.dat"

#define dz2 "file2.dat"

main ()

{ FILE *fp1,*fp2;

double dt=0.002,t[71];int fm1=10,i,j;double zibo[71],xishu[200],dzjl[300];

fp1=fopen(dz1,"w+");

fp2=fopen(dz2,"w+");

fclose(fp1);

fclose(fp2);

fp1=fopen(dz1,"r+");

fp2=fopen(dz2,"r+");

for(i=0;i<=35;i++)

{

t[i+35]=i*dt;

zibo[i+35]=exp(-pow(3.141592653*fm1*t[i+35],2))*(1-2*pow(3.141592653*fm1*t[i+35],2));

}

for(i=0;i<71;i++)

// printf("%f\n",zibo[i]);

for(i=0;i<=34;i++)

{

t[i]=-t[70-i];

zibo[i]=zibo[70-i];

}

for(i=0;i<=70;i++)

{printf("%f\t%f\n",t[i],zibo[i]);

fprintf(fp1," %f\t%f\n",t[i],zibo[i]);

}

for(i=0;i<=199;i++)

{

xishu[i]=0;

}

xishu[10]=0.2;xishu[65]=0.4;xishu[110]=-0.5;xishu[155]=0.7;xishu[190]=-0.36;

for(i=0;i<200;i++)

{ dzjl[i]=0.0;

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

{ if((i-j)>0)

dzjl[i]=dzjl[i]+zibo[j]*xishu[i-j] ; }

}

for(i=0;i<200;i++)

{

// printf("%d\t%f\n",i,dzjl[i]);

fprintf(fp2,"%d\t%f\n",i,dzjl[i]);

}

fclose(fp1);

fclose(fp2);

}

备注:代码可能不能脱离源码包执行，如果不能执行，请在下方留言，我把代码包发给你。当然如果有其他问题，都可以给我留言。