成都理工大学层析成像实验报告

地球物理层析成像实验报告姓名：XXXXXX、学号：20100XXXXXX指导老师：XXXXX电话:XXXXXXXXXX2013年11月8日实习一:Matlab初识与概述部分及编程实验要求：计算下述函数在自变量采样点集上的值：编写脚本程序，分别创建采样向量，计算相应的函数值向量，并利用函数绘图函数fplot(x,y)绘制x-y曲线。Matlab源程序：clear,clcy(i)=-x(i)^2+1;x=-3:0.01:3;elseN=length(x);y(i)=(-x(i)^2+4*x(i)-3)/2;y=zeros(1,N);endfori=1:Nendifx(i)>=-3&x(i)<-1plot(x,y)y(i)=(-x(i)^2-4*x(i)-3)/2;elseifx(i)>=-1&x(i)<1程序运行结果：按要求完成下述练习编写函数maxi，返回给定向量中最大元素的下标值。换言之，imax=maxi(v)返回的imax满足v(imax)不小于任何v(i)。

将你的学号作为测试数据（每一位对应一个向量元素）进行测试functionimax=fun(v)t=v(1,1);imax=1;fori=1:11ifv(1,i+1)>tt=v(1,i+1);imax=i+1;%else%imax1=0;endendimaxend实习二、射线追踪实验内容：模型为一三层模型，范围100m×100m。该模型含有三个水平层，上下两个层速度为3000m/s，中间为一高速层，速度为4500m／s，高速层厚度范围是50～60m。激发点有19个，坐标依次为是El(0,5)、E2(0,10)、E3(0,15)，…，E10(0,95)接收点有21个，坐标分别为Rl(100，0)、R2(100,5)、R3(100,10)，…，R21(100,100)。将此模型网格化为10×10个单元格．根据LTI法追踪出19×21条射线的路径和旅行时。实验源程序：#include<stdio.h>#include<math.h>doubleconv(doubleWa[],doubleRe[],doubleCon[],intn_f,intm_f);#defineA5//定义道数#defineB900//定义角度个数#defineC61//定义子波长度#defineD70//定义采样点个数#definepi3.14159265//定义pi常量voidmain(){doubleh=80.0,v=2000.0,dx=15.0,dt=0.002,f=30.0;//定义第一层厚度,第一层速度，检波点间隔，采样间隔,子波频率//doublea=0.05*pi/180.0;//定义初始角度 doublels,lt,w;//各种数组的定义//intre[A];//定义反射系数序列doubleWa[C+D-1],Re[C+D-1],Con[C+D-1];//定义子波数组,反射系数数组,褶积结果数组doubleangle[B],receive[B],t[A];//定义角度数组,接收点数组，褶积数组doubledata[A][D],S[1],R[A];//定义地震记录数组,炮点数组,检波点数组FILE*fp,*fp1,*fp2,*fp3,*fp4;fp=fopen("子波.txt","w+");fp1=fopen("反射系数.segy","wb+");fp2=fopen("接收点.txt","w+");fp3=fopen("地震记录.segy","wb+");fp4=fopen("接受时间.txt","w+");for(inti=0;i<C;i++) {w=pi*pi*f*f*(i-30)*(i-30)*dt*dt;Wa[i]=(1-2.0*w)*exp((-1.0)*w);fprintf(fp,"%lf\n",Wa[i]);}fclose(fp); S[0]=0.0;h=80.0;for(i=0;i<A;i++){R[i]=95.0+i*dx;} //开始射线追踪 for(i=0;i<B;i++){angle[i]=a*(i+1); receive[i]=2.0*h*tan(angle[i]); printf("%lf\n",receive[i]); for(intj=0;j<A;j++) {if(fabs(receive[i]-R[j])<0.1) {fprintf(fp2,"%lf\n",receive[i]); t[j]=receive[i]/sin(angle[i])/v; re[j]=(int)(t[j]/dt); printf("%d%d\n",j+1,re[j]); fprintf(fp4,"%lf%d\n",t[j],re[j]);}}} fclose(fp2); fclose(fp4); for(i=0;i<A;i++) {for(intj=0;j<D;j++) {if(j==re[i]) Re[j]=0.2; elseRe[j]=0.0; lt=float(Re[j]);fwrite(<,sizeof(float),1,fp1);}//褶积函数的调用// conv(Wa,Re,Con,C,D); for(j=0;j<D;j++) {data[i][j]=Con[j+30]; ls=float(data[i][j]); fwrite(&ls,sizeof(float),1,fp3);}}fclose(fp1); fclose(fp3);}//定义褶积函数doubleconv(doubleWa[],doubleRe[],doubleCon[],intn_f,intm_f){inti,j;intN=n_f+m_f-1;doublesum;for(i=n_f;i<N;i++) Re[i]=0.0;for