工程与环境物探教程实验报告一个反射地震纪录道的形成专业: 勘查技术与工程学号: ********XXXX姓名: XXXXXXX1.解决的问题假设在地面以下半空间内有n+1层弹性介质, 则有n个波阻抗界面, 在第I 层的p波速度为v[i], 密度为p[i], 第i层反射界面上反射系数用R[i]表示, 进行实验一个反射地震纪录道的形成。
已知模型如图所示,p1=2000kg/m3 v1=800m/s p2=1600kg/m3 v1=1000m/sp3=2600kg/m3 v1=2000m/s p4=1800kg/m3 v1=900m/sp5=3200kg/m3 v1=2300m/s h1=h2=h3=h3=50m2. 计算要求1.求每个反射界面的反射系数波阻抗: Z[i]=D[i]*v[i] — D[i]为第i层密度, v[i]为第i层速度反射系数:R[i]=(z[i+1]-z[i]/(z[i+1]+z[i]整理后反射系数表达式如下:R[i]=(V[i+1]*D[i+1]-V[i]*D[i])/(V[i+1]*D[i+1]+V[i]*D[i]) 实验选取选取Ricker （雷克）子波:z[i]=(1-π*f*i*π*f*i)*exp(-(π*f*i)*(π*f*i ）);2.思路方法主导思想: 地震记录道由振幅反射系数序列与子波的褶积构成多层反射界面, 在地面接收的反射波总和的一个地震到记录X （t ）可表示为地层反射系数序列R[i]和地震子波z[i]的褶积（卷积）的结果, 公式表示如下：X(t)=∑R [i ]∗z [i ]=R[itao]N n=1∗b[itao]4.程序代码# include<stdio.h># include<math.h># include<string.h># define pi 3.1416int Wave(float,float ); //定义子波int Reflect(float,float); //定义反射系数int Convolution(float,float,float ); //定义合成记录int main(){float d=0.004,f=30.,h=50; / /设定采样间隔 频率 厚度if(Wave(f,d)!=1) printf("error");if(Reflect(h,d)!=1) printf("error");if(Convolution(f,d,h)!=1) printf("error");return 1;}/***********子波的求取***********/ int Wave(float f,float d){FILE *fpw;if((fpw=fopen("Wave.text","w"))==NULL)printf("error");int i;float Wa[100],t;for(i=0;i<50;i++){t=(50-i)*d;Wa[i]=(float)(1-(pi*pi*f*f*t*t))*(float)exp(-2*pi*pi*f*f*t*t);}for(i=50;i<100;i++){t=(i-50)*d;Wa[i]=(float)(1-(pi*pi*f*f*t*t))*(float)exp(-2*pi*pi*f*f*t*t);}for(i=0;i<100;i++){fprintf(fpw,"%f\n",Wa[i]);}fclose(fpw);return 1;}/*******反射序列的求取*******/ int Reflect(float h,float d){FILE *fpr;fpr=fopen("Rflct.text","w");int i,c=0;float v[5]={800,1000,2000,900,2300};float p[5]={2000,1600,2600,1800,3200};float Re[50];for(i=0;i<50;i++)Re[i]=0;for(i=0;i<4;i++){c+=(int)(h/(v[i]*d));Re[c]=(float)(v[i+1]*p[i+1]-v[i]*p[i])/(v[i+1]*p[i+1]+v[i]*p[i]);}for(i=0;i<50;i++){fprintf(fpr,"%f\n",Re[i]);}fclose(fpr);return 1;}/**********褶积的计算*********/ int Convolution(float f,float d,float h){FILE *fpr,*fpw,*fpc;int i,iao;float w[100],w1[100],r[200],r1[50],Con[150];if((fpw=fopen("Wave.text","r"))==NULL) printf("error");if((fpr=fopen("Rflct.text","r"))==NULL) printf("error");if((fpc=fopen("Convolution.text","w"))==NULL) printf("error");for(i=0;i<100;i++){fscanf(fpw,"%f",&w1[i]);//printf("%f\n",w1[i]);w[100-1-i]=w1[i];//printf("%f\n",w[100-1-i]);}for(i=0;i<50;i++){fscanf(fpr,"%f",&r1[i]);}for(i=0;i<200;i++)r[i]=0;for(i=0;i<50;i++){r[100+i]=r1[i];}for(i=0;i<200;i++)printf("%f\n",r[i]);for(i=0;i<150;i++){Con[i]=0.;for(iao=0;iao<100;iao++)Con[i]=Con[i]+w[iao]*r[i+iao];}for(i=0;i<150;i++){//printf("%lf\n",Con[i]);fprintf(fpc,"%lf\n",Con[i]);}fclose(fpw);fclose(fpr);fclose(fpc);return 1;}5.数据及解释f=40Hz地震记录f=25Hz地震记录f=20Hz地震记录由图可知, 当频率f逐渐变小时, 子波开始产生干涉现象, 导致分辨率降低。
