地震层析成像的正演实验报告.docx

1、地震层析成像的正演实验报告1-1-1正演的速度模型图图1-1-2分块均匀的模型1-2正演的后的走时图图1-3 反演前后速度对比图图1-5-a第0炮，第5接收点的数据 图1-5-b-1正演第1炮，第8接收点（0为开始的激发点，0开始的接收点）图1-5-b-2 与1-5-b-2对应的验证图形（附注：由于本人u盘被病毒入侵，导致本人做得CAD图丢失，此图引用廖松杰同学的CAD图像，但是1-5-b-2由本人程序自己得出，特此说明。）图 1-5-c四边放炮，四边接收左方第2激发， 图1-5-b单边接收第0炮，第25接收的r图像 第5接收点的r数据图正反演的程序单边放炮单边接收：#include#incl

2、ude#includevoid fun1(int n,double R144108,double t1212);void fun2(double k,double o,double t1212,double R144108,int m,int n);/k为斜率，o为炮点坐标，相当于截距；j；double fun(double x1,double y1,double x2,double y2);void main() FILE *fp; int i,j,m,n,l,f; double c12,d12,K1212,r129=0,v129=0.0,t1212=0.0,u,o,w,R144108=0.

3、0,k,v2144=0; float v1;/*/ for(i=0;i12;i+)/第i行 for(j=0;j9;j+)/第j列 vij=3000; v22=5000.0; v32=5000.0; v85=2000.0; v86=2000.0; for(i=0;i12;i+)/第i行 for(j=0;j9;j+)/第j列 v2j*12+i=vij; fp=fopen(速度,wb); for(i=0;i12;i+) for(j=0;j9;j+) v1=vij; fwrite(&v1,sizeof(float),1,fp); fclose(fp);/*/* 计算各点的斜率*/ for(i=0;i1

4、2;i+) /printf(第%d炮的斜率n,i); ci=(i+0.5)*5.0;/*激发点*/ for(j=0;j12;j+) dj=(j+0.5)*5;/*接收点*/ Kij=(dj-ci)/(9.0*3);/*K斜率*/ /printf(K%d%d=%fn,i,j,Kij); /printf(n); /*/ for(i=0;i12;i+)/第i炮 for(j=0;j12;j+)/第j接收点 if(Kij=0)/平行于x轴,该行所在的每一个网格均经过,路程都是3 fun1(i,R,t); printf(t%d%d=%fn,i,j,tij); else if(Kij!=0) k=Kij;

5、o=ci; m=i;n=j; fun2(k,o,t,R,i,j); printf(t%d%d=%lfn,i,j,tij); fp=fopen(time.txt,w); for(i=0;i12;i+) for(j=0;j12;j+) fprintf(fp,%lfn,tij); fclose(fp); fp=fopen(系数矩阵R的值.txt,w); for(i=0;i144;i+) for(j=0;j108;j+) fprintf(fp,%ft,Rij); fprintf(fp,n); fclose(fp); fp=fopen(原来的速度值.txt,w); for(j=0;j108;j+) fp

6、rintf(fp,%ft,v2j); fclose(fp); /*/*/*当斜率k为0的时候，计算走时t的值*/*/void fun1(int n,double R144108,double t1212) FILE *fp1; double b=0.0; int i=0,j=0,q=0;/循环变量 double r129=0.0,v129;/* for(i=0;i12;i+)/第i行 for(j=0;j9;j+)/第j列 vij=3000; v22=5000.0; v32=5000.0; v85=2000.0; v86=2000.0; for(j=0;j9;j+) rnj=3.0; /*/*写

7、出检验r的值*/ /*fp1=fopen(r的值.txt,w); for(i=0;i12;i+) for(j=0;j9;j+) fprintf(fp1,%ft,Rij); fprintf(fp1,n); fclose(fp);*/ /*/*/ for(i=0;i12;i+)/第i行 for(j=0;j9;j+)/第j列 b+=rij*(1/vij); tnn=b; for(i=0;i12;i+)/第i行 for(j=0;j9;j+)/第j列 Rn*12+nq+=rij; double fun(double x1,double y1,double x2,double y2)double s;s=

8、(y2-y1)*(y2-y1)+(x2-x1)*(x2-x1);return sqrt(s);/*/*/void fun2(double k,double o,double t1212,double R144108,int m,int n)/k为斜率，o为炮点坐标，相当于截距； FILE *fp2; int i=0,j=0,q=0;/循环变量 int w1,w2,w3,w4;/中间变量，用来判断点在分块均匀上的位置 double p=0,v129=0.0,r129=0.0; double x1,y1,x2,y2,x3,y3,x4,y4; float r1;/* v22=5000.0; v32=

9、5000.0; v85=2000.0; v86=2000.0; for(i=0;i12;i+)/第i行 for(j=0;j9;j+)/第j列 vij=3000; rij=0; for(i=0;i12;i+) for(j=0;j=(j*3.0)&(x1=(j*3.0)&(x2=(i*5.0)&(y3=(i*5.0)&(y4=(i+1)*5.0)&(x4=(j+1)*3.0);/右方 /* /计算路径长度r，当有两个点存在时，有下面的六种情况。 if(w1!=0&w3!=0) rij=fun(x1,y1,x3,y3); else if(w1!=0&w2!=0) rij=fun(x1,y1,x2,y

10、2); else if(w2!=0&w4!=0) rij=fun(x2,y2,x4,y4); else if(w4!=0&w1!=0) rij=fun(x1,y1,x4,y4); else if(w3!=0&w2!=0) rij=fun(x2,y2,x3,y3); else if(w3=1&w4=1) rij=fun(x3,y3,x4,y4); /* /走时等于射线经度每一个单元格的时间与慢速和（1/v）的累加 p+=(rij*(1/vij); tmn=p; for(i=0;i12;i+) for(j=0;j9;j+) Rm*12+nq+=rij; /*/*检验第0炮，第五接收点r的正确性*/

11、 if(m=0&n=5) fp2=fopen(r2的值.txt,w); for(i=0;i12;i+) for(j=0;j9;j+) fprintf(fp2,%ft,rij); fprintf(fp2,n); fclose(fp2); fp2=fopen(r2的值,wb); for(i=0;i12;i+) for(j=0;j9;j+) r1=rij; fwrite(&r1,sizeof(float),1,fp2); fclose(fp2); /*/ 四边放炮四边接收：#include#include#includevoid fun1(int n,double R120108,double t4

12、30,int m );void fun2(double k,double o,double t430,double R120108,int m,int n);/k为斜率，o为炮点坐标，相当于截距；double fun(double x1,double y1,double x2,double y2);void fun3(int n,double R1133108,double t433,int m);void fun4(double k,double o,double t1433,double R1133108,int m,int n);/k为斜率，o为炮点坐标，相当于截距main() FILE

13、 *fp; int i,j,j1,m,n,l,e; double c112,c29,c312,c49,c12,c59; double d112,d29,d312,d49; double K1430=0.0,K2433=0.0,K3430=0.0,K4433=0.0; double K1242,r129=0,v129,t430=0.0,t1433=0.0,u,o,w,R120108=0.0,R1133108=0.0,k; /*/*/* 计算各点的斜率*/*左方激发*/ for(i=0;i4;i+) c1i=12.5+10*i;/*左激发点*/ for(j1=0;j130;j1+) if(j19)

14、 d2j1=1.5+j1*3.0;/*上接收点*/ K1ij1=-c1i/d2j1; else if(j1=9) d3j1-9=2.5+(j1-9)*5.0;/*右接收点*/ K1ij1=(d3j1-9-c1i)/27.0; else d4j1-21=1.5+(j1-21)*3.0;/*下接收点*/ K1ij1=(60-c1i)/d3j1-21; / printf(K1%d%d=%fn,i,j1,K1ij1); /*上方激发* for(i=0;i4;i+) /printf(第%d炮的斜率n,i); c2i=4.5+6.0*i;/*上激发点*/ for(j1=0;j133;j1+) if(j112) d3j1=2.5+(j1)*5.0;/*右接收点*/ K2ij1=d3j1/(27-c2i); else if(j1=12) d4j1-12=1.5+(j1-12)*3.0;/*下接收点*/ K2ij1=60/(d4j1-12-c2i);