成都理工大學層析成像實驗報告

1、本科生實驗報告實驗課程 地球物理層析成像 學院名稱 地球物理學院 專業(yè)名稱 勘查技術與工程 學生姓名 學生學號 指導教師 曹俊興 實驗地點 5417 實驗成績 二一五 年 三 月 二一五 年 四 月學生實驗 心得 在學習了地球物理層析成像之后，收獲了很多專業(yè)知識，比如學會了利用層析成像的手段反演出地下地質體的異常，同時也學會了利用我們的專業(yè)知識解決不同的地質問題。程序語言作為一種工具一方面起到了輔助作用，另一方面我們也學會了一種思維方式，如何設計程序，如何用程序解決我們的復雜問題。在今后的學習工作當中，進一步拓寬思路，勇于創(chuàng)新，能夠獲得更多的知識。學生（簽名）： 2015年 4 月 28 日指

2、導教師評語成績評定：指導教師（簽名）： 年 月 日地震走時層析成像實驗地球物理正反演概論課程結業(yè)報告學號：201205060423姓名：馬力衡專業(yè)：勘查技術與工程手機要運用C語言程序，正演得到地震走時和射線在傳播過程中經過離散化處理單元格內的距離。通過反演程序反演出地下異常速度值，將反演所得速度值成圖與原始速度成圖進行比較，得出結論。離散化處理模型建立，單邊激發(fā),四邊激發(fā)直射線正演，單邊激發(fā)，四邊激發(fā)反演異常值，用代數重建算法迭代慢度矩陣對單邊，四邊激發(fā)進行迭代。關鍵詞：離散化 反演 迭代 第1章 地震走時層析成像實驗1.1實驗內容1.1.1直射線正演：使用直射線追蹤

3、方法計算走時的正演；分塊均勻模型1.1.1.1單邊激發(fā)正演程序：#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;/激發(fā)點

4、與接收點的縱坐標for(i=0;i<12;i+)y_jfi=1.5+3.0*i;/激發(fā)點點坐標的方程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四邊激發(fā)正演程序；#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;/左側激發(fā)點坐標 double x_js4=0.0,y_js4;/右側激發(fā)點的坐標 double X_jf4,Y_begin4=0.0;/上下激發(fā)點坐標 double X_end4,Y_js4=0.0;/下頂下激發(fā)點坐標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;/左側激發(fā)/double xl_zuo1;double b_zuo1;/左側激發(fā)時截距double xl_zuo2;double b_zuo2;/左側激發(fā)時截距double xl_zuo3;double b_zuo3;/左側激發(fā)時截距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); /右側激發(fā)/double xl_you1;double b_you1;/右側激發(fā)時截距double xl_you2;double b_y

17、ou2;/右側激發(fā)時截距double xl_you3;double b_you3;/右側激發(fā)時截距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); /上側激發(fā)/double xl_shang1;double b_shang1;/上側激發(fā)時截距double xl_shang2;double b_shang2;/上側激發(fā)時截距double xl_shang3;double b_shang3;/上側激發(fā)時截距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); /下側激發(fā)/double xl_xia1;double b_xia1;/下側激發(fā)時截距double xl_xia2;double b_xia2;/下側激發(fā)時截距double xl_xia3;double b_xia3;/下側激發(fā)時截距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 單邊激發(fā)反演程序；#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