工程與環境物探教程實驗報告一個反射地震紀錄道的形成專業：勘查技術與工程學號：20100506XXXX姓名：XXXXXXX解決的問題假設在地面以下半空間內有n+1層彈性介質，則有n個波阻抗界面，在第I層的p波速度為v[i]，密度為p[i]，第i層反射界面上反射系數用R[i]表示，進行實驗一個反射地震紀錄道的形成。已知模型如圖所示，p1=2000kg/m3v1=800m/sp2=1600kg/m3v1=1000m/sp3=2600kg/m3v1=2000m/sp4=1800kg/m3v1=900m/sp5=3200kg/m3v1=2300m/sh1=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)=n=14.程序代碼#include<stdio.h>#include<math.h>#include<string.h>#definepi3.1416intWave(float,float);//定義子波intReflect(float,float);//定義反射系數intConvolution(float,float,float);//定義合成記錄intmain(){ floatd=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"); return1;}/***********子波的求取***********/intWave(floatf,floatd){ FILE*fpw; if((fpw=fopen("Wave.text","w"))==NULL)printf("error"); inti; floatWa[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);return1;}/*******反射序列的求取*******/intReflect(floath,floatd){ FILE*fpr; fpr=fopen("Rflct.text","w"); inti,c=0; floatv[5]={800,1000,2000,900,2300}; floatp[5]={2000,1600,2600,1800,3200}; floatRe[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);return1;}/**********褶積的計算*********/intConvolution(floatf,floatd,floath){ FILE*fpr,*fpw,*fpc; inti,iao; floatw[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);return1;}5.數據及解釋f=40Hz地震記錄f=25Hz地震記錄f=20Hz地震記錄由圖可知，當頻率f逐漸變小時，子波開始產生干涉現象，導致分辨率降低。當中間某層速度低于上下兩高速層時，