地震層析成像的正演實驗報告

1、成都理工大學正反演實驗報告嚴光明1-1-1正演的速度模型圖圖1-1-2分塊均勻的模型1-2正演的后的走時圖圖1-3 反演前后速度對比圖圖1-5-a第0炮，第5接收點的數(shù)據(jù) 圖1-5-b-1正演第1炮，第8接收點（0為開始的激發(fā)點，0開始的接收點）圖1-5-b-2 與1-5-b-2對應的驗證圖形（附注：由于本人u盤被病毒入侵，導致本人做得CAD圖丟失，此圖引用廖松杰同學的CAD圖像，但是1-5-b-2由本人程序自己得出，特此說明。）圖 1-5-c四邊放炮，四邊接收左方第2激發(fā)， 圖1-5-b單邊接收第0炮，第25接收的r圖像 第5接收點的r數(shù)據(jù)圖正反演的程序單邊放炮單邊接收：#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;/*激發(fā)點*/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軸,該行所在的每一個網(wǎng)格均經(jīng)過,路程都是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("系數(shù)矩陣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;/循環(huán)變量 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;/循環(huán)變量 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，交點四位于右邊 /*/*判斷射線是否經(jīng)過分塊均勻的網(wǎng)格點上,四個交點是否在網(wǎng)格的四條邊上*/*/*注意：網(wǎng)格的上下兩條邊y值相等,網(wǎng)格的左右兩邊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);/* /走時等于射線經(jīng)度每一個單元格的時間與慢速和（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;/*/*/*計算各點的斜率*/*左方激發(fā)*/for(i=0;i<4;i+)c1i=12.5+10*i;/*左激發(fā)點*/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);/*上方激發(fā)*for(i=0;i<4;i+)/printf("第%d炮的斜率n",i); c2i=4.5+6.0*i;/*上激發(fā)點*/ 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");/*右方激發(fā)*/ for(i=0;i<4;i+) /printf("第%d炮的斜率n",i); c3i=12.5+10*i;/*右激發(fā)點*/ 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"); /*下方激發(fā)* for(i=0;i<4;i+) /printf("第%d炮的斜率n",i); c4i=4.5+6

21、*i;/*下激發(fā)點*/ 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軸,該行所在的每一個網(wǎng)格均經(jīng)過,路程都是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("系數(shù)矩陣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軸,該行所在的每一個網(wǎng)格均經(jīng)過,路程都是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("系數(shù)矩陣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軸,該行所在的每一個網(wǎng)格均經(jīng)過,路程都是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("系數(shù)矩陣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軸,該行所在的每一個網(wǎng)格均經(jīng)過,路程都是3 m=j; fun3(i,R1,t1,m); /printf("t%d%d=%