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

1、/才声Z乂庶本 科生实验报告实验课程地球物理层析成像学院名称地球物理学院专业名称勘查技术与工程学生姓名学生学号指导教师曹俊兴实验地点5417实验成绩二。一五年三月二。一五年四月在学习了地球物理层析成像之后，收获了很多专业知识，比如学 生 实 验 心学会了利用层析成像的手段反演出地下地质体的异常，同时也学会 了利用我们的专业知识解决不同的地质问题。程序语言作为一种工 具一方面起到了辅助作用，另一方面我们也学会了一种思维方式， 如何设计程序，如何用程序解决我们的复杂问题。在今后的学习工 作当中，进一步拓宽思路，勇于创新，能够获得更多的知识。得学生（签名）：指 导 教 师 评 语成绩评定：指导教师（

2、签名）：年 月 日2015年4月 28日地震走时层析成像实验地球物理正反演概论课程结业报告学号：201205060423姓名：马力衡专业：勘查技术与工程手机要运用c语言程序，正演得到地震走时和射线在传播过程中经过离散化处理单元格 内的距离。通过反演程序反演出地下异常速度值，将反演所得速度值成图与原始 速度成图进行比较，得出结论。离散化处理模型建立，单边激发，四边激发直射 线正演，单边激发，四边激发反演异常值，用代数重建算法迭代慢度矩阵对单边, 四边激发进行迭代。关键词：离散化反演迭代第1章地震走时层析成像实验*3.0*3.03.0夜)3.0-%3.0*3.03.03.

3、03.03.03.03.03.03.0 ；*3.03.05.03.03.03.03.03.03-0 i*303.05.03.03.03.03.03.03.0*30303.03.03.03.03.03.03.0,3.03.03.03.03.03.03.03Q*3.03.03.03.03.03.03.03.03.0*3.03.03.03.03.03.03.03.01 3.00*3.03.03.03.03.02.02.03.03-013.03.03.03.03.03.03.03.03.0*3.03.03.03.03.03.03.03.03.03.0 3.03.0 3.0L3.03.0 3.0 3.0

4、 实验内容单元数：9x12单元边长：3.0 x5. 0m参数:速度 km/s)实验一：单边激发单边接收(左发右收)数据：12 x 12实验二：四边激发,四边接收(每边设置4个激发点.在其它 三边所有接收点接收：上下两边 的激发点位置参见紫色匿叉图标 所示；直立两边的激发点位置参 见蓝色星状图标所示)数据：4x(4x12)+2x(4 x 9)+ 2 x(4 x12)+ 4x(4 x 9)直射线正演：使用直射线追踪方法计算走时的正演；分块均匀模型单边激发正演程序：Sinclude Sinclude void main()(int v12 9;int m, n, i, j;FILE *fpO;fpO

5、=fopen(“速度.txt, r);for(i=0;i12;i+)(for(j=0;j9;j+)fscanf (fpO, %d”, &vi j);double b12; 截距；double xl 12 12; /斜率;double y_jf 12, y_js12;激发点与接收点的纵坐标 for(i=0;i12;i+)(y_jf i=L 5+3. o*i；激发点点坐标的方程for(j=0;j12;j+)y_jsj=1.5+3. 0*j;接收点坐标的方程 xli j = (y_jsj-y_jf i)/(45-0);斜率 printfC%fn, xlij);for(i=0;i12;i+)(bi=L

6、5+i*3;每条射线截距)以上在求射线的斜率和射线在纵轴上的截距double ft_t=O. 0;每一格的时间；double fl12 12 12 9; 每一格射线的长度;double Time12 12;每条射线的时间；double X0,Y0;第一个点坐标；double XI, Y1;第二个点坐标；double x_0, x_l,y_0, y_l; /判定的 x, y;double xO, xl, yO, yl;小格的边界；FILE *fp_ds;fp_ds=fopen(每一小格的距离.dat, w);for(i=0;i12;i+) (for(j=0;j12;j+)Timeij=0.0;f

7、or (n=0;n12;n+)(for(m=0;m9;m+)(fli j n m=0;x0=5*m;xl=x0+5;y0=3*n;yl=y0+3;y_O=Xlij*xO+bi;if(y_0=y0)X0=x0;Y0=y_0;y_l=xlij*xl+bi;if (y_K=yl)& (y_l=y0) (Yl=y_l;Xl=xl;elsex_l=(yO-bi)/xli j;if(x_l=xO)(Yl=yO;XI=x_l;else(x_l=(yl-bi)/xlij;Xl=x_l;Yl=yl;)flijnm=sqrt(Xl-XO)*(Xl-XO)+(Yl-YO)*(Yl-YO);) else (y_l=x

8、lij*xl+bi;if(y_l二yl&y_l=yO) (Xl=xl;Y1=y_l;x_O=(yO-bi)/xlij;if(x_0=x0)(X0=x_0;YO=yO; else (x_O=(yl-bi)/xlij;X0=x_0;YO=yl;flijnm=sqrt(Xl-XO)*(Xl-XO) + (Yl-YO)*(Yl-YO); )elsex_O=(yO-bi)/xlij;int i, j, m, n;if(x_O=xO)XO=x_O;YO=yO;x_l=(yl-bi)/xlij;XI=x_l;Yl=yl;flijnm=sqrt(Xl-XO)*(Xl-XO)+(Yl-YO)*(Yl-YO);)

9、else flijnm=O. 0;)fprintf(fp_ds, *%f , flijnm);Timei j+=fl i j nm;每个单元格的走时；printf(z，%d ”, vn m);射线传播的总时间)fprintf(fp_ds, n);)以上判定射线存在的单元格情况FILE *fpl;fpl二fopen(走时.txt, w);for(i=0;i12;i+)(for(j=0;j12;j+)(fprintf(fpl, Timeij);)/fprintf (fp, Vn);以上得出走时并写入txt文件中输出1.1.1.2四边激发正演程序；ftinclude#include#included

10、ouble funfpol(double xl,double b,FILE *fp,FILE *fpl)double fl129;double XO, YO;double XI, Yl；double xO,xl,yO,yl;double x_0, x_l, y_0, y_l;double v129, Time=0. 0;for(i=0;i12;i+)(for(j=0;j9;j+)(vij=3.0;v2 2 =5.0;v3 2 =5.0;v8 5 =2.0;v8 6 =2.0;for(n=0;n12;n+)(y0=3*n;yl=y0+3;for(m=0;m9;m+)(fl nm=0. 0;x0=

11、5*m；xl=x0+5;y_0=xl*x0+b;if(y_0=y0)X0=x0; Y0=y_0;y_l=xl*xl+b;if(y_l=yO)Yl=y_l;Xl=xl;)else(x_l=(yO-b)/xl;if(x_l=xO)(Yl=yO; XI=x_l;)else(x_l=(yl-b)/xl;XI=x_l;Yl=yl;)fln m =sqrt(Xl-XO)*(Xl-XO) + (Yl-YO)*(Yl-YO);)else(y_l=xl*xl+b;if(y_l=y0)(Xl=xl; Yl=y_l;x_0=(y0-b)/xl;if(x_0=x0)X0=x_0; YO=yO;elsex_0=(yl-

12、b)/xl;XO=x_O; YO=yl;)flnm=sqrt(Xl-XO)*(Xl-XO)+(Yl-YO)*(Yl-YO);else(x_0=(y0-b)/xl;if(x_0=x0)(XO=x_O; YO=yO;x_l=(yl-b)/xl;Xl=x_l; Yl=yl;flnm=sqrt(Xl-XO)*(Xl-XO)+(Yl-YO)*(Yl-YO);)fprintf (fp, z，%f ”, fl n m);Time+=flnm/vnm;fprintf(fp, Xn*);fprintf (fpl,，z%fn , Time);return 0;)void main()(int i, j;以下是读取

13、原始速度值FILE *fp=fopen(每个单元格的距离.txt，w)，*fpl=fopen (时间.txt，w)；double x_jf 4 = 0. 0, y_jf 4;左侧激发点坐标double x_js4 = 0. 0, y_js4;右侧激发点的坐标double X_jf 4, Y_begin4 = 0. 0;上下激发点坐标double X_end4, Y_js4 = 0. 0;下顶下激发点坐标 for(i=0;i4;i+)(for(j=0;j4;j+)(y_jfi=7. 5+6*i;x_jsj=45. 0;y_jsj=7. 5+6*j;X_jfi=7. 5+10. 0*i;X_end

14、j=7. 5+10. 0*j;Y_jsj=36.0;Jdouble x_zuo12 = 0. 0, y_zuo12;左侧接收点坐标double x_you12 = 0. 0, y_you12;右侧接收点坐标 double x_shang9, y_shang9 = 0. 0;上侧接收点坐标 double x_xia9, y_xia9 = 0. 0;下侧接收点坐标 for(i=0;i12;i+)(for(j=0;j9;j+)x_youi=45. 0;y_zuoi=l. 5+3*i;y_youi=l. 5+3*i;y_xiaj=36. 0;x_shangj=2. 5+5*j；x_xiaj=2. 5+

15、5*j;)左侧激发double xl_zuol;double b_zuol;左侧激发时截距double xl_zuo2;double b_zuo2;左侧激发时截距double xl_zuo3;double b_zuo3;左侧激发时截距for(i=0;i4;i+)(for(j=0;j9;j+)(xl_zuo1=(y_shangj-y_j fi)/(x_shangj-x_j fi);b_zuol=10. 5+6*i;funfpol( xl_zuol, b_zuol, fp, fpl);xl_zuo3=(y_xiaj-y_jfi)/(x_xiai);b_zuo3=10. 5+6*i;funfpol

16、( xl_zuo3, b_zuo3, fp, fpl);Jfor(i=0;i4;i+)(for(j=0;j12;j+)xl_zuo2=(y_youj-y_jfi)/(x_youj一x_j fi);b_zuo2=10. 5+6*i;funfpol( xl_zuo2, b_zuo2, fp, fpl);)右侧激发double xl_youl;doub 1 e b_you 1; 右侧激发时截距double xl_you2;double b_you2;右侧激发时截距double xl_you3;doub 1 e b_you3; 右侧激发时截距for(i=0;i4;i+)(for(j=0J9;j+)(x

17、l_youl=(y_shangj-y-jsi)/(x_shangj-x_j si);b_you1=y_jsi-xl_you1*x_jsi;funfpol (xl_youl, b_youl, fp, fpl);xl_you3= (y_xia j-y_js i)/(x_xiaj-x_jsi);b_you3=y_jsi-xl_you3*x_jsi;funfpol(xl_you3, b_you3, fp, fpl);)for(i=0;i4;i+)(for(j=0;j12;j+)(xl_you2=(y_zuoj-y_jsi)/(x_zuoj-x_jsi);b_you2=y_jsi-xl_you2*x_j

18、si;funfpol(xl_you2, b_you2, fp, fpl);上侧激发double xl_shangl;double b_shangl;上侧激发时截距double xl_shang2;double b_shang2;上侧激发时截距double xl_shang3;double b_shang3;上侧激发时截距for(i=0;i4;i+)(for(j=0;j12;j+)(xl_shangl=(y_zuoj-Y_begini)/(x_zuoi);b_shangl=Y_begini-xl_shangl*X_jfi;funfpol(xl_shangl, b_shangl, fp, fpl)

19、;xl_shang2=(y_youj-Y_begini)/(x_youj-X_jfi);b_shang2=Y_begini-xl_shang2*X_jfi;funfpol(xl_shang2, b_shang2, fp, fpl);)for(i=0;i4;i+)(for(j=0;j9;j+)(if(x_xiaj=X_jfi)(xl_shang3=INT_MAX;b_shang3=Y_begini-xl_shang3*X_jfi;elsefor(j=0;j9;j+)xl_shang3= (y_xiaj-Y_begini)/ (x_xiai);b_shang3=Y_begini-xl_shang3

20、*X_jfi;)funfpol(xl_shang3, b_shang3, fp, fpl);)下侧激发double xl_xial;double b_xial;下侧激发时截距double xl_xia2;double b_xia2;下侧激发时截距double xl_xia3;double b_xia3;下侧激发时截距for(i=0;i4;i+)(for(j=0;j12;j+)(xl_xial=(y_zuoj-Y_jsi)/(x_zuoj-X_endi);b_xial=Y_jsi-xl_xial*X_endi;funfpol( xl_xial, b_xial, fp, fpl);xl_xia2=

21、(y_youj-Y_jsi)/(x_youj-X_endi);b_xia2=Y_jsi-xl_xia2*X_endi;funfpol(xl_xia2, b_xia2, fp, fpl);)for(i=0;i4;i+)i f (x_shang j =X_end i)xl xia3=INT MAX;b_xia3=Y_jsi-xl_xia3*X_endi;else(xl_xia3=(y_shangj_Y_jsi)/(x_shangj-X_endi);b_xia3=Y_jsi-xl_xia3*X_endi;funfpol(xl_xia3, b_xia3, fp, fpl);)1.1.2 反演(矩阵方程

22、求解)：单边激发反演程序；ftinclude#include void main() (double rO, dO1212129, al, a2, Timel212, md129;int N, i, j, n, m;double M129;double wucha12 12 = 0. 0;反演误差FILE *fp7, *fp8, *fp9, *fp10;fp7=fopen(每一小格的距离.dat”, r)；fp8=fopen (走时.txt, r);fplO二fopen(不为零的个数.txt，w)；for(i=0;i12;i+)(for(j=0;j12;j+)(for(n=0;n12;n+)(

23、for(m=0;m9;m+)fscanf(fp7, %lf ”, &d0ijnm);for(j=0;j9;j+)/*for(i=0;i12;i+)(for (尸0;j12;j+)(for (n=0;n12;n+)(for(m=0;m9;m+)printf (/z%5. 21fn，dOi j n m);)水/for(i=0;i12;i+)(for(j=0;j12;j+)fscanf(fp8, STimeij);)/*for(i=0;i12;i+)(for(j=0;j12;j+)(printf (zz%5. 61fn, Timei j);水/for (n=0;n12;n+)(for (m:=0;m

24、9;m+)(Mn m=0.0;mdn m =0. 0;) for (n=0;n12;n+) (for(m=0;m9;m+)for(i=0;i12;i+)(for(j=0;j12;j+)(if(dOijn m!=0)Mn m+=l;算出系数矩阵中每一列不为零值的数的个数)fprintf(fplO, Mnm);)N=0; 迭代次数10000次计算误差while(N10000)(for(i=0;i12;i+)(for(j=0;j12;j+)r0=0. 0;for (n=0;n12;n+)(for (m:z0;m9;m+)r0+=d0ijnm*mdnm;wuchaEi j = (Timei j-rO)

25、;/printf wuchaEi j);)重建算法迭代慢度矩阵for(m=0;m12;m+)(for (n=0;n9;n+)(al=0. 0;a2=:0. 0for(i=0;i12;i+)(for(j=0;j12;j+)(al+=d0ijmn*d0ijmn;a2+=(wuchaij m n);)mdmn+=a2/(al*Mmn);)N+;以txt格式输出最后数据fp9=fopen (最终结果.txt，w);for(i=0;i Figure 1FileEditView InsertToolsDesktopWindowHe IpMBMBMB S。端|Q|%二等夏说县| 口目| 口也迭代伯次100次

26、2.935739 3.2812253. 111551 2. 336008 3. 658779 2. 551901 3. 316450 3. 458434 2. 8662933. 048761 3.0958812. 813202 2. 569809 3. 645845 2. 738465 3. 239735 3. 093387 2. 9846862. 950545 3.1712883.618861 3. 042182 3. 629389 2. 830167 3. 151999 3.019176 3.0103132. 930405 3. 1031313. 672115 3. 252944 3. 3

27、64247 2. 975754 3. 091184 2. 991957 3. 0063282. 937960 2. 9369752. 730594 3. 053914 3. 260941 3. 084954 3. 071011 3. 028878 2. 9228312. 879351 2. 9399442. 898927 2.913109 3. 212950 3. 122887 3. 105629 3. 108858 2. 8470062. 831008 2. 9853592. 967733 2. 898719 3. 103858 3. 145873 3. 206954 3. 138705 2

28、. 7978112. 798233 2. 9621103. 009536 2. 934635 2. 972843 3. 152150 3. 475847 2. 976618 2. 7961062. 809070 2.9235133.018924 2. 928858 2.918971 2. 133096 2.234360 2. 786387 2. 8552692. 843132 2. 8710103. 098449 2. 664298 3.431973 3. 169985 3. 347160 2. 794868 2. 8956092. 880881 2. 9707202.945991 2.611

29、897 3. 752789 2. 851127 3. 353616 2. 990678 2. 8859622. 806888 3.2828453. 145242 2. 354578 3. 746505 2. 659282 3.394118 3. 386711 2. 751432成图;1000 次2.9527533. 1723992.9482422.4681173.8310072. 646069 3. 163270 3. 359202 2. 8659572.9818863. 1287042.6494502.7545473. 7410132.715717 3.206425 3.149545 2.9

30、394712.9654693. 1104583. 8760022. 9724663. 5038222. 914696 3. 102020 3. 0673292.9685012.9526013. 0718773.9486603.0178513.4121273. 0044563. 0878823. 0248392. 9716942.9417153.0089302.6815413.0111923.2770463. 0939503.1042293.0110642.9396562.9105602.9712332.8044052.9713813.1650863. 1378013.1714213.01158

31、12. 8978472. 8723762.9600322.9109502.9404373. 0750383. 1578023.2677352.999200 2.8648342.8397902.9442683. 0002062.9220483.0031373. 1683153.3884512.9450332.8540082.8184982.9353303. 0434502.8914243. 0272522. 0701652.1874282.9088172.8503462.8116862.9398483. 0595172. 7719963.2891983. 0414193. 4458422.915

32、0242.8414822.8116652.9622583.0568542.6402313.6801222. 9053643.4022392.9592432.8288222.8013472.9952373. 1094192.5153854.0051262.8450643. 303550 3. 0295942.809051迭代10000次2.9996383. 1001652. 665187 2. 630206 4. 133789 2. 894126 2. 850578 3.217386 2. 9448792.9779283. 0977792.6200222.8165373. 6505672. 92

33、24043.036146 3.124486 2. 9488322.9570583. 0757244.1072782.8654033.4884362. 9996443. 0595193.0872552.9400882.9332463. 0527204.2102612. 9083633. 3558883. 0513293. 0999433.0565872.9299202.9130013. 0305842.7440812.9434303.2382453. 0963223. 1485893.0277002.9158262.8912263.0111132. 7938442. 9659993.151931

34、3. 1247473.2058333.0042232.8991072.8715762.9916312.8480552.9674833. 106060 3. 1294893. 2700392.9831812.8819322.8541682.9719852.9063222.9465233.1009673. 1101683. 3404442.962240 2. 8667602.8358812. 9561242.9639512.9105863.1312592.0321512.175356 2.945506 2.8500282.8206082. 9394463. 0200672.8636013.1960

35、793. 021726 3. 488907 2.926017 2. 8350352. 8080052.9166073.0786922.7962583.3353132. 939619 3.569501 2.906135 2.8226712. 8021352. 8763523.2067612.5673303.9850212.752465 3.596445 2. 889599 2.811198100000 次2.9796623. 1232972. 5793382. 7150484.1748832.8557512.915260 3.125875 2. 9802532.9615283. 1023442.6199582.9180083. 3973423.