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

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&l

2、t;stdio.h>#include<stdlib.h>#include<math.h>void 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,K12

3、12,r129=0,v129=0.0,t1212=0.0,u,o,w,R144108=0.0,k,v2144=0;float v1;/*/for(i=0;i<12;i+)/第i行 for(j=0;j<9;j+)/第j列 vij=3000; v22=5000.0; v32=5000.0; v85=2000.0; v86=2000.0; for(i=0;i<12;i+)/第i行 for(j=0;j<9;j+)/第j列 v2j*12+i=vij; fp=fopen("速度","wb");for(i=0;i<12;i+)for(j=

4、0;j<9;j+)v1=vij; fwrite(&v1,sizeof(float),1,fp);fclose(fp);/*/*计算各点的斜率*/for(i=0;i<12;i+)/printf("第%d炮的斜率n",i);ci=(i+0.5)*5.0;/*激发点*/for(j=0;j<12;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;i<12;i

5、+)/第i炮 for(j=0;j<12;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;i<12;i+) for(j=0;j<12;j+)fprintf

6、(fp,"%lfn",tij);fclose(fp);fp=fopen("系数矩阵R的值.txt","w"); for(i=0;i<144;i+) for(j=0;j<108;j+) fprintf(fp,"%ft",Rij); fprintf(fp,"n"); fclose(fp); fp=fopen("原来的速度值.txt","w"); for(j=0;j<108;j+)fprintf(fp,"%ft",v2j);

7、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;i<12;i+)/第i行 for(j=0;j<9;j+)/第j列 vij=3000; v22=5000.0; v32=5000.0; v85=2000.0; v86=2000.0; for(j=0;j<9;j+) rnj=3.0; /*/*写出检验r的值*/ /

8、*fp1=fopen("r的值.txt","w"); for(i=0;i<12;i+) for(j=0;j<9;j+) fprintf(fp1,"%ft",Rij); fprintf(fp1,"n"); fclose(fp);*/*/*/ for(i=0;i<12;i+)/第i行 for(j=0;j<9;j+)/第j列 b+=rij*(1/vij); tnn=b; for(i=0;i<12;i+)/第i行 for(j=0;j<9;j+)/第j列 Rn*12+nq+=rij; do

9、uble 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;

10、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;i<12;i+)/第i行 for(j=0;j<9;j+)/第j列 vij=3000; rij=0; for(i=0;i<12;i+) for(j=0;j<9;j+) y1=i*5.0; x1=(y1-o)/k;/计算交点1,由y计算x,交点一位于上边 y2=(i+1)*5.0; x2=(y2-o)/k;/计算交点1,由y计算x，交点2位于下边 x3=j*3.0;y3=k*x

11、3+o;/计算交点3,由x计算y，交点3位于左边x4=(j+1)*3.0;y4=k*x4+o;/计算交点3,由x计算y，交点四位于右边 /*/*判断射线是否经过分块均匀的网格点上,四个交点是否在网格的四条边上*/*/*注意：网格的上下两条边y值相等,网格的左右两边x的值相等*/*w1=(x1>=(j*3.0)&&(x1<=(j+1)*3.0)&&(y1=(i*5.0);/上方 w2=(x2>=(j*3.0)&&(x2<=(j+1)*3.0)&&(y2=(i+1)*5.0);/下方 w3=(y3>=(i

12、*5.0)&&(y3<=(i+1)*5.0)&&(x3=(j*3.0);/左方 w4=(y4>=(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

13、,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;i<12;i+)for(j=0;j<9;j+)Rm*12+nq+=rij;/*/*检验第0炮，第五接收点r的正确性*/if(m=0&&n=

14、5)fp2=fopen("r2的值.txt","w"); for(i=0;i<12;i+) for(j=0;j<9;j+) fprintf(fp2,"%ft",rij); fprintf(fp2,"n"); fclose(fp2); fp2=fopen("r2的值","wb"); for(i=0;i<12;i+) for(j=0;j<9;j+) r1=rij; fwrite(&r1,sizeof(float),1,fp2); fclose(fp

15、2); /*/ 四边放炮四边接收：#include<stdio.h>#include<stdlib.h>#include<math.h>void 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 R

16、1133108,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,

17、o,w,R120108=0.0,R1133108=0.0,k;/*/*/*计算各点的斜率*/*左方激发*/for(i=0;i<4;i+)c1i=12.5+10*i;/*左激发点*/for(j1=0;j1<30;j1+)if(j1<9)d2j1=1.5+j1*3.0;/*上接收点*/K1ij1=-c1i/d2j1;else if(j1<21&&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

18、-c1i)/d3j1-21;/printf("K1%d%d=%fn",i,j1,K1ij1);/*上方激发*for(i=0;i<4;i+)/printf("第%d炮的斜率n",i); c2i=4.5+6.0*i;/*上激发点*/ for(j1=0;j1<33;j1+) if(j1<12)d3j1=2.5+(j1)*5.0;/*右接收点*/K2ij1=d3j1/(27-c2i); else if(j1<21&&j1>=12) d4j1-12=1.5+(j1-12)*3.0;/*下接收点*/ K2ij1=60/(

19、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;i<4;i+) /printf("第%d炮的斜率n",i); c3i=12.5+10*i;/*右激发点*/ for(j1=0;j1<30;j1+) if(j1<9) d4j1=1.5+(j1)*3.0;/*下接收点*/ K3ij1=(60-c3i)/(

20、d4j1-27); else if(j1<21&&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;i<4;i+) /printf("第%d炮的斜率n",i); c4i=4.5+6

21、*i;/*下激发点*/ for(j1=0;j1<33;j1+) if(j1<12) d1j1=2.5+5.0*(j1);/*左接收点*/ K4ij1=(60-d1j1)/c4i; else if(j1<21&&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);

22、 /printf("n"); /*/*/*/*/*计算R1的值和t1*/for(i=0;i<4;i+)/第i炮 for(j=0;j<30;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

23、); fp=fopen("time1.txt","w"); for(i=0;i<4;i+) for(j=0;j<30;j+) fprintf(fp,"%lfn",tij); fclose(fp); fp=fopen("系数矩阵R1的值.txt","w"); for(i=0;i<120;i+) for(j=0;j<108;j+) fprintf(fp,"%ft",Rij); fprintf(fp,"n"); fclose(fp);/*

24、计算R3和t3*/ for(i=0;i<4;i+)/第i炮 for(j=0;j<30;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",&

25、quot;w"); for(i=0;i<4;i+) for(j=0;j<30;j+) fprintf(fp,"%lfn",tij); fclose(fp); fp=fopen("系数矩阵R3的值.txt","w"); for(i=0;i<120;i+) for(j=0;j<108;j+) fprintf(fp,"%ft",Rij); fprintf(fp,"n"); fclose(fp);/*计算K2和t2*/ for(i=0;i<4;i+)/第i炮 fo

26、r(j=0;j<33;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;i<4

27、;i+) for(j=0;j<33;j+) fprintf(fp,"%lfn",t1ij); fclose(fp); fp=fopen("系数矩阵R2的值.txt","w"); for(i=0;i<120;i+) for(j=0;j<108;j+) fprintf(fp,"%ft",R1ij); fprintf(fp,"n"); fclose(fp);/*计算K4和t4*/ for(i=0;i<4;i+)/第i炮 for(j=0;j<33;j+)/第j接收点 Kij=K4ij; if(Kij=0)/平行于x轴,该行所在的每一个网格均经过,路程都是3 m=j; fun3(i,R1,t1,m); /printf("t%d%d=%