一、实验目的

通过编制和运行相关C或matlab程序，理解地震子波的概念和特点，掌握地震合成记录的制作流程。

1. 认识子波，对子波的波形有直观的认识。（零相位子波，混合相位子波，最小相位子波；了解子波的分辨率与频宽的关系；）
   
2. 利用褶积公式合成一维地震记录
   

二、主要内容和原理

1、理论子波的计算及显示    

    （1）零相位子波、最小相位子波、钟型子波的计算

零相位子波的表达式

最小相位子波的表达式为；

n=m1/m2，为最大波峰m1和最大波谷m2之比，取2.0或1.5，本次试验采用2.0.；代表子波的中心频率， t=i*dt，dt为时间采样间隔，i为时间离散点序号;

本次试验通过零相位子波和最小相位子波的公式，利用0.002s的采样间隔，对零相位子波和最小相位子波进行离散采样，得到一组数据。利用excel绘制图像，得到中心频率分别为10,25,40,100hz的8张图像。

2、利用声波测井资料计算层速度，密度和反射系数序列

本实验通过C语言编程序处理获得的声波测井资料，把资料中的数据用于层速的计算，公式为

和密度计算，当速度与密度呈线性关系时，计算公式为：

本次试验假设速度和密度为线性关系，采用这个公式参与计算。

3、 制作合成地震记录

地震合成记录是在一定的假设条件下进行的，如果地下实际介质有n层面都是反射界面的话，地面可以接收到每一层界面的反射波，于是一个实际地震记录道上会记录n个反射波。忽略介质吸收、透射损失影响，使用地震子波和反射系数离散卷积的方法得到。因为离散波形卷积需要采样间隔相同，所以使用0.002s的间隔对第i个反射面到地表的垂直往返时间进行离散采样.计算公式为：

得到反射系数的序号，然后求取每个序号求取反射系数，当然并不是所有的序号都有反射系数，对没有反射系数的用0取代。得到反射系数的离散数据，在得到35hz和70hz的零相位子波的离散数据，根据离散卷积公式，

由此可得一系列随时间变化的卷积结果即为对应时刻的波形振幅值，绘制散点图即可观察到所记录的理论合成波记录。

三、实验结果分析（含图）

1、理论子波的计算及显示

（1）、零相位子波的波形

中心频率为10HZ

中心频率为25HZ

中心频率为40HZ

中心频率为100HZ

结果分析：根据零相位子波的计算结果和所绘制的图像，我们可以发现当n一定时（n=m1/m2,为最大波峰m1和最大波谷m2之比），零相位子波图像随着子波的中心频率fm的增加，而呈现有规律的变化。对于零相位子波来说，子波的中心频率越大，最大幅度值越小。

（2）、最小相位子波的波形

中心频率为10HZ

中心频率为25HZ

中心频率为40HZ

中心频率为100HZ

结果分析：根据最小相位子波的计算结果和所绘制的图像，我们可以发现当n一定时（n=m1/m2,为最大波峰m1和最大波谷m2之比），对于最小相位子波来说，随着子波中心频率的增加，振幅变化的速度加快了，但是最大幅值却没有明显的变化。

2、利用声波测井资料计算层速度，密度和反射系数序列

（1）、层速曲线

结果分析：通过编程计算得到层速度结果并绘制图件，分析结果和图件可以发现，层速度和深度并没有明显的线性关系，层速度应该与每层介质的密度有关。

（2）、密度曲线

结果分析：密度和深度并没有明显的线性关系，但是所有地层除了个别层密度较大外，基本变化不大。密度应该和层速度有密切的联系。

（3）、反射系数

结果分析：反射系数是根据离散取样而来，所以许多标号下，反射系数为零，因为相邻两层的的介质性质不同，所以，反射系数也不相同，有大有小，与相邻底层的波阻抗性质有关。

3、 制作合成地震记录

（1）、中心频率为35HZ

结果分析：

（2）、中心频率为70HZ

结果分析：观察分析两幅地震合成记录图件地震合成记录70Hz相比35Hz其最大幅值要大并且相同水平距离70Hz的波数大于35Hz，原因在于同一介质频率越高波长越短使得相同水平距离70Hz的波数大于35Hz，同样频率越高在相同时间间隔内叠加次数越多使得卷结果绝对值增大。

四、源程序代码

1、理论子波的计算及显示

#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;int fm1=10,fm2=25,fm3=40,fm4=100,i,j;double r1,r2,r3,r4,w1,w2,w3,w4;

fp1=fopen(dz1,"w+");//文件新建

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

fclose(fp1);

fclose(fp2);//关闭文件

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

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

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

{ t=i*dt;

r1=exp(-pow(3.141592653*fm1*t,2))*(1-2*pow(3.141592653*fm1*t,2));//

     r2=exp(-pow(3.141592653*fm2*t,2))*(1-2*pow(3.141592653*fm2*t,2));

r3=exp(-pow(3.141592653*fm3*t,2))*(1-2*pow(3.141592653*fm3*t,2));

     r4=exp(-pow(3.141592653*fm4*t,2))*(1-2*pow(3.141592653*fm4*t,2));

     w1=exp(-2*fm1*fm1*t*t*log(2))*sin(2*3.14*fm1*t);

w2=exp(-2*fm2*fm2*t*t*log(2))*sin(2*3.14*fm2*t);

w3=exp(-2*fm3*fm3*t*t*log(2))*sin(2*3.14*fm3*t);

     w4=exp(-2*fm4*fm4*t*t*log(2))*sin(2*3.14*fm4*t);//计算子波

         printf("%d %f %f %f\n",i,t,r1,w1);

    fprintf(fp1,"%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n",t,r1,r2,r3,r4,w1,w2,w3,w4);

}//输出子波

    fclose(fp1);

fclose(fp2);

}

2、利用声波测井资料计算层速度，密度和反射系数序列

3、制作合成地震记录

#include<stdio.h>//速度密度，反射系数和制作地震记录程序在一块

#include<stdlib.h>//程序中有标明

#include<math.h>

main(){

FILE *fp1,*fp2,*fp3;double vi[200],midu1[200],midu2[200],ti[200],rn[2000];int n[200];int s1[200];int s2[200];float s3[200];

double dt=0.002,t;int fm1=35,fm2=70,i,j;double r1[200],w1[200],r2[200],w2[200],xn1[2000],xn2[2000];

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

{      t=(j+1)*dt;

     r1[j]=exp(-pow(3.141592653*fm1*t,2))*(1-2*pow(3.141592653*fm1*t,2));

r2[j]=exp(-pow(3.141592653*fm2*t,2))*(1-2*pow(3.141592653*fm2*t,2));

     w1[j]=exp(-2*fm2*fm2*t*t*log(2))*sin(2*3.141592653*fm1*t);

     w2[j]=exp(-2*fm2*fm2*t*t*log(2))*sin(2*3.141592653*fm2*t); }//计算子波

if((fp1=fopen("file1","rb+"))==NULL){//第二小问层速，密度程序开始

printf("Cannot open file strike any key exit!"); }

if((fp2=fopen("file2","r+"))==NULL){

printf("Cannot open file strike any key exit!");

}

if((fp3=fopen("file3","r+"))==NULL){

printf("Cannot open file strike any key exit!"); }

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

fscanf(fp1,"%d %d %f \n",&s1[i],&s2[i],&s3[i]);/读取文件中数据

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

{    vi[i]=1000000.0/s2[i];//层速计算

midu1[i]=1.62+0.00021*vi[i];

midu2[i]=0.31*pow(vi[i],0.25); }//密度计算，后边结果采用线性密度公式

ti[0]=2*s3[0]/vi[0];//初值

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

{ ti[i]=ti[i-1]+2*(s3[i]-s3[i-1])/vi[i]; } //时间

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

{ n[i]=int(ti[i]/0.002); }

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

     {          rn[i]=0.0;                  }

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

{rn[(n[i])]=(midu1[i+1]*vi[i+1]-midu1[i]*vi[i])/(midu1[i+1]*vi[i+1]+midu1[i]*vi[i] ;}//反射系数计算第二小问层速，密度，反射系数程序结束

for(i=0;i<2000;i++)//70hz子波卷积，地震记录合成

{    xn1[i]=0.0;

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

     { if((i-j)>0)

            xn1[i]=xn1[i]+r1[j]*rn[i-j] ;     }

}

for(i=0;i<2000;i++)//70hz子波卷积，地震记录合成

{    xn2[i]=0.0;

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

     { if((i-j)>0)

            xn2[i]=xn2[i]+r2[j]*rn[i-j] ;}

}

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

{

fprintf(fp2,"%d\t%f\t%f\t%f\t%f\t%d\n",i+1,vi[i],midu1[i],midu2[i],ti[i],n[i]);

}//输出层速，密度到文件

for(i=0;i<2000;i++)// 输出反射系数，卷积到文件

{ fprintf(fp3,"%d\t%f\t%f\t%f\n",i+1,rn[i],xn1[i],xn2[i]); }

printf("处理完毕请打开文件寻找数据\n");

fclose(fp1);

fclose(fp2);

fclose(fp3); }