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

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

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

3、方法计算走时的正演；分块均匀模型1.1.1.1单边激发正演程序：#include <stdio.h>#include <math.h>void main()int v129;int m,n，i，j;FILE *fp0;fp0=fopen("速度.txt","r");for(i=0;i<12;i+)for(j=0;j<9;j+)fscanf(fp0,"%d",&vij); double b12; /截距; double xl1212; /斜率;double y_jf12,y_js12;/激发点

4、与接收点的纵坐标for(i=0;i<12;i+)y_jfi=1.5+3.0*i;/激发点点坐标的方程for(j=0;j<12;j+)y_jsj=1.5+3.0*j;/接收点坐标的方程xlij=(y_jsj-y_jfi)/(45-0);/斜率printf("%fn",xlij);for(i=0;i<12;i+)bi=1.5+i*3;/每条射线截距/以上在求射线的斜率和射线在纵轴上的截距 /double ft_t=0.0; /每一格的时间;double fl1212129; /每一格射线的长度;double Time1212; /每条射线的时间;double

5、X0,Y0; /第一个点坐标;double X1,Y1; /第二个点坐标;double x_0,x_1,y_0,y_1; /判定的 x，y;double x0,x1,y0,y1; /小格的边界;FILE *fp_ds;fp_ds=fopen("每一小格的距离.dat","w");for(i=0;i<12;i+)for(j=0;j<12;j+)Timeij=0.0; for(n=0;n<12;n+)for(m=0;m<9;m+)flijnm=0;x0=5*m;x1=x0+5;y0=3*n;y1=y0+3;y_0=xlij*x0+bi

6、;if(y_0<=y1&&y_0>=y0)X0=x0;Y0=y_0;y_1=xlij*x1+bi;if(y_1<=y1)&&(y_1>=y0)Y1=y_1;X1=x1;elsex_1=(y0-bi)/xlij;if(x_1<=x1&&x_1>=x0)Y1=y0;X1=x_1;elsex_1=(y1-bi)/xlij;X1=x_1;Y1=y1;flijnm=sqrt(X1-X0)*(X1-X0)+(Y1-Y0)*(Y1-Y0);elsey_1=xlij*x1+bi;if(y_1<=y1&&y

7、_1>=y0)X1=x1;Y1=y_1;x_0=(y0-bi)/xlij;if(x_0<=x1&&x_0>=x0)X0=x_0;Y0=y0;elsex_0=(y1-bi)/xlij;X0=x_0;Y0=y1;flijnm=sqrt(X1-X0)*(X1-X0)+(Y1-Y0)*(Y1-Y0);elsex_0=(y0-bi)/xlij;if(x_0<=x1&&x_0>=x0)X0=x_0;Y0=y0;x_1=(y1-bi)/xlij;X1=x_1;Y1=y1; flijnm=sqrt(X1-X0)*(X1-X0)+(Y1-Y0)*(Y

8、1-Y0);else flijnm=0.0;fprintf(fp_ds,"%f ",flijnm);Timeij+=flijnm/vnm;/每个单元格的走时；printf("%d ",vnm);/射线传播的总时间 fprintf(fp_ds,"n");/以上判定射线存在的单元格情况/FILE *fp1; fp1=fopen("走时.txt","w"); for(i=0;i<12;i+)for(j=0;j<12;j+)fprintf(fp1,"%fn",Timeij

9、); /fprintf(fp,"/n");/以上得出走时并写入txt文件中输出/ 1.1.1.2四边激发正演程序；#include<stdio.h>#include<math.h>#include<limits.h>double funfpo1(double xl,double b,FILE *fp,FILE *fp1)int i,j,m,n; double fl129; double X0,Y0; double X1,Y1; double x0,x1,y0,y1; double x_0,x_1,y_0,y_1; double v129,

10、Time=0.0; for(i=0;i<12;i+) for(j=0;j<9;j+) vij=3.0;v22=5.0;v32=5.0;v85=2.0;v86=2.0; for(n=0;n<12;n+) y0=3*n;y1=y0+3; for(m=0;m<9;m+) flnm=0.0;x0=5*m;x1=x0+5;y_0=xl*x0+b;if(y_0<=y1&&y_0>=y0) X0=x0;Y0=y_0; y_1=xl*x1+b; if(y_1<=y1)&&(y_1>=y0) Y1=y_1;X1=x1; else x

11、_1=(y0-b)/xl; if(x_1<=x1&&x_1>=x0) Y1=y0; X1=x_1; else x_1=(y1-b)/xl; X1=x_1;Y1=y1; flnm=sqrt(X1-X0)*(X1-X0)+(Y1-Y0)*(Y1-Y0);else y_1=xl*x1+b;if(y_1<=y1&&y_1>=y0) X1=x1; Y1=y_1; x_0=(y0-b)/xl; if(x_0<=x1&&x_0>=x0) X0=x_0; Y0=y0; else x_0=(y1-b)/xl; X0=x_0; Y

12、0=y1; flnm=sqrt(X1-X0)*(X1-X0)+(Y1-Y0)*(Y1-Y0);else x_0=(y0-b)/xl;if(x_0<=x1&&x_0>=x0) X0=x_0; Y0=y0;x_1=(y1-b)/xl;X1=x_1; Y1=y1; flnm=sqrt(X1-X0)*(X1-X0)+(Y1-Y0)*(Y1-Y0);fprintf(fp,"%f ",flnm);Time+=flnm/vnm; fprintf(fp,"n"); fprintf(fp1,"%fn ",Time); ret

13、urn 0;void main() int i,j;/以下是读取原始速度值/FILE *fp=fopen("每个单元格的距离.txt","w"), *fp1=fopen("时间.txt","w"); double x_jf4=0.0,y_jf4;/左侧激发点坐标 double x_js4=0.0,y_js4;/右侧激发点的坐标 double X_jf4,Y_begin4=0.0;/上下激发点坐标 double X_end4,Y_js4=0.0;/下顶下激发点坐标for(i=0;i<4;i+)for(j=0;j

14、<4;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_endj=7.5+10.0*j; Y_jsj=36.0;double 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;i<12;i+)for(j=0;j<9;j+) x_youi=45.0; y_zuoi=1.

15、5+3*i;y_youi=1.5+3*i;y_xiaj=36.0; x_shangj=2.5+5*j; x_xiaj=2.5+5*j;/左侧激发/double xl_zuo1;double b_zuo1;/左侧激发时截距double xl_zuo2;double b_zuo2;/左侧激发时截距double xl_zuo3;double b_zuo3;/左侧激发时截距for(i=0;i<4;i+) for(j=0;j<9;j+) xl_zuo1=(y_shangj-y_jfi)/(x_shangj-x_jfi); b_zuo1=10.5+6*i; funfpo1( xl_zuo1,b

16、_zuo1,fp,fp1); xl_zuo3=(y_xiaj-y_jfi)/(x_xiaj-x_jfi); b_zuo3=10.5+6*i; funfpo1( xl_zuo3,b_zuo3,fp,fp1); for(i=0;i<4;i+) for(j=0;j<12;j+) xl_zuo2=(y_youj-y_jfi)/(x_youj-x_jfi); b_zuo2=10.5+6*i; funfpo1( xl_zuo2,b_zuo2,fp,fp1); /右侧激发/double xl_you1;double b_you1;/右侧激发时截距double xl_you2;double b_y

17、ou2;/右侧激发时截距double xl_you3;double b_you3;/右侧激发时截距for(i=0;i<4;i+) for(j=0;j<9;j+) xl_you1=(y_shangj-y_jsi)/(x_shangj-x_jsi);b_you1=y_jsi-xl_you1*x_jsi;funfpo1(xl_you1,b_you1,fp,fp1); xl_you3=(y_xiaj-y_jsi)/(x_xiaj-x_jsi);b_you3=y_jsi-xl_you3*x_jsi;funfpo1(xl_you3,b_you3,fp,fp1); for(i=0;i<4;

18、i+) for(j=0;j<12;j+) xl_you2=(y_zuoj-y_jsi)/(x_zuoj-x_jsi); b_you2=y_jsi-xl_you2*x_jsi; funfpo1(xl_you2,b_you2,fp,fp1); /上侧激发/double xl_shang1;double b_shang1;/上侧激发时截距double xl_shang2;double b_shang2;/上侧激发时截距double xl_shang3;double b_shang3;/上侧激发时截距for(i=0;i<4;i+) for(j=0;j<12;j+) xl_shang1

19、=(y_zuoj-Y_begini)/(x_zuoj-X_jfi);b_shang1=Y_begini-xl_shang1*X_jfi; funfpo1(xl_shang1,b_shang1,fp,fp1); xl_shang2=(y_youj-Y_begini)/(x_youj-X_jfi);b_shang2=Y_begini-xl_shang2*X_jfi; funfpo1(xl_shang2,b_shang2,fp,fp1); for(i=0;i<4;i+) for(j=0;j<9;j+) if(x_xiaj=X_jfi) xl_shang3=INT_MAX; b_shang

20、3=Y_begini-xl_shang3*X_jfi; else xl_shang3=(y_xiaj-Y_begini)/(x_xiaj-X_jfi); b_shang3=Y_begini-xl_shang3*X_jfi; funfpo1(xl_shang3,b_shang3,fp,fp1); /下侧激发/double xl_xia1;double b_xia1;/下侧激发时截距double xl_xia2;double b_xia2;/下侧激发时截距double xl_xia3;double b_xia3;/下侧激发时截距for(i=0;i<4;i+) for(j=0;j<12;

21、j+) xl_xia1=(y_zuoj-Y_jsi)/(x_zuoj-X_endi); b_xia1=Y_jsi-xl_xia1*X_endi; funfpo1( xl_xia1, b_xia1,fp,fp1); xl_xia2=(y_youj-Y_jsi)/(x_youj-X_endi); b_xia2=Y_jsi-xl_xia2*X_endi; funfpo1(xl_xia2,b_xia2,fp,fp1); for(i=0;i<4;i+) for(j=0;j<9;j+) if(x_shangj=X_endi)xl_xia3=INT_MAX;b_xia3=Y_jsi-xl_xia

22、3*X_endi;else xl_xia3=(y_shangj-Y_jsi)/(x_shangj-X_endi); b_xia3=Y_jsi-xl_xia3*X_endi; funfpo1(xl_xia3,b_xia3,fp,fp1); 1.1.2 反演（矩阵方程求解）： 1.1.2.1 单边激发反演程序；#include<stdio.h>#include<math.h>void main() double r0, d01212129, a1,a2,Time1212，md129; int N，i,j,n,m; double M129;double wucha1212=0

23、.0;/反演误差FILE *fp7,*fp8,*fp9,*fp10; fp7=fopen("每一小格的距离.dat","r"); fp8=fopen("走时.txt","r"); fp10=fopen("不为零的个数.txt","w");for(i=0;i<12;i+) for(j=0;j<12;j+) for(n=0;n<12;n+) for(m=0;m<9;m+) fscanf(fp7,"%lf ",&d0ijnm);

24、/*for(i=0;i<12;i+) for(j=0;j<12;j+) for(n=0;n<12;n+) for(m=0;m<9;m+) printf("%5.2lfn",d0ijnm); */ for(i=0;i<12;i+) for(j=0;j<12;j+) fscanf(fp8,"%lf",&Timeij); /*for(i=0;i<12;i+) for(j=0;j<12;j+) printf("%5.6lfn",Timeij); */ for(n=0;n<12;n+

25、) for(m=0;m<9;m+) Mnm=0.0; mdnm=0.0; for(n=0;n<12;n+) for(m=0;m<9;m+) for(i=0;i<12;i+) for(j=0;j<12;j+) if(d0ijnm!=0)Mnm+=1; /算出系数矩阵中每一列不为零值的数的个数 fprintf(fp10,"%fn",Mnm); N=0; /迭代次数10000次 /计算误差/ while(N<10000) for(i=0;i<12;i+) for(j=0;j<12;j+) r0=0.0; for(n=0;n<1

26、2;n+) for(m=0;m<9;m+) r0+=d0ijnm*mdnm; wuchaij=(Timeij-r0); /printf("%f",wuchaij); /重建算法迭代慢度矩阵 / for(m=0;m<12;m+) for(n=0;n<9;n+) a1=0.0; a2=0.0; for(i=0;i<12;i+) for(j=0;j<12;j+) a1+=d0ijmn*d0ijmn; a2+=(wuchaij*d0ijmn); mdmn+=a2/(a1*Mmn); N+; /以txt格式输出最后数据/ fp9=fopen("

27、最终结果.txt","w"); for(i=0;i<12;i+)for(j=0;j<9;j+) fprintf(fp9,"%lf ",1.0/mdij); fprintf(fp9,"n"); 数据；10次2.784751 3.280478 3.261611 2.454797 3.711088 2.493625 3.307422 3.441376 2.771530 3.145613 3.042454 2.956347 2.564969 3.632866 2.761524 3.229051 3.030151 3.01

28、5298 3.024139 3.211870 3.507230 3.014290 3.717386 2.895093 3.202817 2.976226 2.986879 2.889704 3.105180 3.377056 3.284153 3.382822 2.999874 3.057579 3.020697 2.938980 2.886266 2.913581 2.910847 3.054787 3.280013 3.066888 3.096661 3.088917 2.967574 2.805710 2.954641 2.924597 2.918879 3.156665 3.08183

29、5 3.116152 3.125369 2.938478 2.843217 2.995171 2.971224 2.908500 3.053036 3.126632 3.201618 3.182704 2.910128 2.839169 2.927154 2.960580 2.883011 3.025253 2.984553 3.225567 3.049169 2.823326 2.884555 3.017693 2.989972 2.914996 2.886734 2.328265 2.426117 2.750940 2.851556 2.853319 2.855602 3.080046 2

30、.621081 3.484502 3.060188 3.186040 2.730749 2.849618 2.961355 2.888923 3.022542 2.756842 3.839824 2.798330 3.292418 2.982118 2.857067 2.660109 3.406152 3.320099 2.485385 3.758592 2.526407 3.350241 3.421877 2.677194成图；100次2.935739 3.281225 3.111551 2.336008 3.658779 2.551901 3.316450 3.458434 2.86629

31、3 3.048761 3.095881 2.813202 2.569809 3.645845 2.738465 3.239735 3.093387 2.984686 2.950545 3.171288 3.618861 3.042182 3.629389 2.830167 3.151999 3.019176 3.010313 2.930405 3.103131 3.672115 3.252944 3.364247 2.975754 3.091184 2.991957 3.006328 2.937960 2.936975 2.730594 3.053914 3.260941 3.084954 3

32、.071011 3.028878 2.922831 2.879351 2.939944 2.898927 2.913109 3.212950 3.122887 3.105629 3.108858 2.847006 2.831008 2.985359 2.967733 2.898719 3.103858 3.145873 3.206954 3.138705 2.797811 2.798233 2.962110 3.009536 2.934635 2.972843 3.152150 3.475847 2.976618 2.796106 2.809070 2.923513 3.018924 2.92

33、8858 2.918971 2.133096 2.234360 2.786387 2.855269 2.843132 2.871010 3.098449 2.664298 3.431973 3.169985 3.347160 2.794868 2.895609 2.880881 2.970720 2.945991 2.611897 3.752789 2.851127 3.353616 2.990678 2.885962 2.806888 3.282845 3.145242 2.354578 3.746505 2.659282 3.394118 3.386711 2.751432成图；1000次

34、2.952753 3.172399 2.948242 2.468117 3.831007 2.646069 3.163270 3.359202 2.865957 2.981886 3.128704 2.649450 2.754547 3.741013 2.715717 3.206425 3.149545 2.939471 2.965469 3.110458 3.876002 2.972466 3.503822 2.914696 3.102020 3.067329 2.968501 2.952601 3.071877 3.948660 3.017851 3.412127 3.004456 3.0

35、87882 3.024839 2.971694 2.941715 3.008930 2.681541 3.011192 3.277046 3.093950 3.104229 3.011064 2.939656 2.910560 2.971233 2.804405 2.971381 3.165086 3.137801 3.171421 3.011581 2.897847 2.872376 2.960032 2.910950 2.940437 3.075038 3.157802 3.267735 2.999200 2.864834 2.839790 2.944268 3.000206 2.9220

36、48 3.003137 3.168315 3.388451 2.945033 2.854008 2.818498 2.935330 3.043450 2.891424 3.027252 2.070165 2.187428 2.908817 2.850346 2.811686 2.939848 3.059517 2.771996 3.289198 3.041419 3.445842 2.915024 2.841482 2.811665 2.962258 3.056854 2.640231 3.680122 2.905364 3.402239 2.959243 2.828822 2.801347

37、2.995237 3.109419 2.515385 4.005126 2.845064 3.303550 3.029594 2.809051迭代10000次2.999638 3.100165 2.665187 2.630206 4.133789 2.894126 2.850578 3.217386 2.944879 2.977928 3.097779 2.620022 2.816537 3.650567 2.922404 3.036146 3.124486 2.948832 2.957058 3.075724 4.107278 2.865403 3.488436 2.999644 3.059

38、519 3.087255 2.940088 2.933246 3.052720 4.210261 2.908363 3.355888 3.051329 3.099943 3.056587 2.929920 2.913001 3.030584 2.744081 2.943430 3.238245 3.096322 3.148589 3.027700 2.915826 2.891226 3.011113 2.793844 2.965999 3.151931 3.124747 3.205833 3.004223 2.899107 2.871576 2.991631 2.848055 2.967483 3.106060 3.129489 3.270039 2.983181 2.881932 2.854168 2.971985 2.906322 2.946523 3.100967 3.110168 3.340444 2.962240 2.866760 2.835881 2.956124 2.963951 2.910586 3.131259 2.032151 2.175356 2.945506