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

成都理工大学正反演实验报告严光明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#include#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.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;i12;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;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+)fprintf(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; /*/*写出检验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=(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=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,y2);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的正确性*/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 t430,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 *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)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;elsed4j1-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); elsed1j1-21=2.5+5.0*(j1-21);/*左接收点*/K2ij1=-(d1j1-21/c2i); /printf(K2%d%d=%fn,i,j1,K2ij1);/printf(n);/*右方激发*/ for(i=0;i4;i+) /printf(第%d炮的斜率n,i); c3i=12.5+10*i;/*右激发点*/ for(j1=0;j130;j1+) if(j19) d4j1=1.5+(j1)*3.0;/*下接收点*/ K3ij1=(60-c3i)/(d4j1-27); else if(j1=9) d1j1-9=2.5+5.0*(j1-9);/*左接收点*/ K3ij1=(c3i-d1j1-9)/27; else d2j1-21=1.5+(j1-21)*3.0;/*上接收点*/ K3ij1=c3i/(27-d2j1-21); / printf(K3%d%d=%fn,i,j1,K3ij1); /printf(n); /*下方激发* for(i=0;i4;i+) /printf(第%d炮的斜率n,i); c4i=4.5+6*i;/*下激发点*/ for(j1=0;j133;j1+) if(j112) d1j1=2.5+5.0*(j1);/*左接收点*/ K4ij1=(60-d1j1)/c4i; else if(j1=12) d2j1-12=1.5+(j1-12)*3.0;/*上接收点*/ K4ij1=60/(c4i-d2j1-12); else d3j1-21=2.5+5.0*(j1-21);/*右接收点*/ K4ij1=(60-d3j1-21)/(c4i-27); /printf(K4%d%d=%fn,i,j1,K4ij1); /printf(n); /*/*/*/*/*计算R1的值和t1*/for(i=0;i4;i+)/第i炮 for(j=0;j30;j+)/第j接收点 Kij=K1ij; if(Kij=0)/平行于x轴,该行所在的每一个网格均经过,路程都是3 m=j; fun1(i,R,t,m); /printf(t%d%d=%fn,i,j,tij); else if(Kij!=0) k=Kij; o=ci=c1i; m=i;n=j; fun2(k,o,t,R,i,j); /printf(t%d%d=%lfn,i,j,tij); fp=fopen(time1.txt,w); for(i=0;i4;i+) for(j=0;j30;j+) fprintf(fp,%lfn,tij); fclose(fp); fp=fopen(系数矩阵R1的值.txt,w); for(i=0;i120;i+) for(j=0;j108;j+) fprintf(fp,%ft,Rij); fprintf(fp,n); fclose(fp);/*计算R3和t3*/ for(i=0;i4;i+)/第i炮 for(j=0;j30;j+)/第j接收点 Kij=K3ij; if(Kij=0)/平行于x轴,该行所在的每一个网格均经过,路程都是3 m=j; fun1(i,R,t,m); /printf(t%d%d=%fn,i,j,tij); else if(Kij!=0) k=Kij; o=ci=c2i; m=i;n=j; fun2(k,o,t,R,i,j); / printf(t%d%d=%lfn,i,j,tij); fp=fopen(time3.txt,w); for(i=0;i4;i+) for(j=0;j30;j+) fprintf(fp,%lfn,tij); fclose(fp); fp=fopen(系数矩阵R3的值.txt,w); for(i=0;i120;i+) for(j=0;j108;j+) fprintf(fp,%ft,Rij); fprintf(fp,n); fclose(fp);/*计算K2和t2*/ for(i=0;i4;i+)/第i炮 for(j=0;j33;j+)/第j接收点 Kij=K2ij; if(Kij=0)/平行于x轴,该行所在的每一个网格均经过,路程都是3 m=j; fun3(i,R1,t1,m); /printf(t%d%d=%fn,i,j,t1ij); else if(Kij!=0) k=Kij; o=c5i=c2i; m=i;n=j; fun4(k,o,t1,R1,i,j); /printf(t%d%d=%lfn,i,j,t1ij); fp=fopen(time2.txt,w); for(i=0;i4;i+) for(j=0;j33;j+) fprintf(fp,%lfn,t1ij); fclose(fp); fp=fopen(系数矩阵R2的值.txt,w); for(i=0;i120;i+) for(j=0;j108;j+) fprintf(fp,%ft,R1ij); fprintf(fp,n); fclose(fp);/*计算K4和t4*/ for(i=0;i4;i+)/第i炮 for(j=0;j33;j+)/第j接收点 Kij=K4ij; if(Kij=0)/平行于x轴,该行所在的每一个网格均经过,路程都是3 m=j; fun3(i,R1,t1,m); /printf(t%d%d=%fn,i