成都理工大學(xué)層析成像實驗報告

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

2、導(dǎo)教師評語成績評定：指導(dǎo)教師（簽名）： 年 月 日地震走時層析成像實驗地球物理正反演概論課程結(jié)業(yè)報告學(xué)號：201205060423姓名：馬力衡專業(yè)：勘查技術(shù)與工程手機要運用C語言程序，正演得到地震走時和射線在傳播過程中經(jīng)過離散化處理單元格內(nèi)的距離。通過反演程序反演出地下異常速度值，將反演所得速度值成圖與原始速度成圖進行比較，得出結(jié)論。離散化處理模型建立，單邊激發(fā),四邊激發(fā)直射線正演，單邊激發(fā)，四邊激發(fā)反演異常值，用代數(shù)重建算法迭代慢度矩陣對單邊，四邊激發(fā)進行迭代。關(guān)鍵詞：離散化 反演 迭代 第1章 地震走時層析成像實驗1.1實驗內(nèi)容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;/左側(cè)激發(fā)點坐標 double x_js4=0.0,y_js4;/右側(cè)激發(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;/左側(cè)接收點坐標double x_you12=0.0,y_you12;/右側(cè)接收點坐標double x_shang9,y_shang9=0.0;/上側(cè)接收點坐標double x_xia9,y_xia9=0.0;/下側(cè)接收點坐標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;/左側(cè)激發(fā)/double xl_zuo1;double b_zuo1;/左側(cè)激發(fā)時截距double xl_zuo2;double b_zuo2;/左側(cè)激發(fā)時截距double xl_zuo3;double b_zuo3;/左側(cè)激發(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); /右側(cè)激發(fā)/double xl_you1;double b_you1;/右側(cè)激發(fā)時截距double xl_you2;double b_y

17、ou2;/右側(cè)激發(fā)時截距double xl_you3;double b_you3;/右側(cè)激發(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); /上側(cè)激發(fā)/double xl_shang1;double b_shang1;/上側(cè)激發(fā)時截距double xl_shang2;double b_shang2;/上側(cè)激發(fā)時截距double xl_shang3;double b_shang3;/上側(cè)激發(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); /下側(cè)激發(fā)/double xl_xia1;double b_xia1;/下側(cè)激發(fā)時截距double xl_xia2;double b_xia2;/下側(cè)激發(fā)時截距double xl_xia3;double b_xia3;/下側(cè)激發(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("不為零的個數(shù).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; /算出系數(shù)矩陣中每一列不為零值的數(shù)的個數(shù) fprintf(fp10,"%fn",Mnm); N=0; /迭代次數(shù)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格式輸出最后數(shù)據(jù)/ fp9=fopen("

27、最終結(jié)果.txt","w"); for(i=0;i<12;i+)for(j=0;j<9;j+) fprintf(fp9,"%lf ",1.0/mdij); fprintf(fp9,"n"); 數(shù)據(jù)；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