成都理工大學(xué)層析成像實(shí)驗(yàn)報(bào)告

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

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

3、方法計(jì)算走時(shí)的正演；分塊均勻模型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ā)點(diǎn)

4、與接收點(diǎn)的縱坐標(biāo)for(i=0;i<12;i+)y_jfi=1.5+3.0*i;/激發(fā)點(diǎn)點(diǎn)坐標(biāo)的方程for(j=0;j<12;j+)y_jsj=1.5+3.0*j;/接收點(diǎn)坐標(biāo)的方程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; /每一格的時(shí)間;double fl1212129; /每一格射線的長度;double Time1212; /每條射線的時(shí)間;double

5、X0,Y0; /第一個(gè)點(diǎn)坐標(biāo);double X1,Y1; /第二個(gè)點(diǎn)坐標(biāo);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;/每個(gè)單元格的走時(shí)；printf("%d ",vnm);/射線傳播的總時(shí)間 fprintf(fp_ds,"n");/以上判定射線存在的單元格情況/FILE *fp1; fp1=fopen("走時(shí).txt","w"); for(i=0;i<12;i+)for(j=0;j<12;j+)fprintf(fp1,"%fn",Timeij

9、); /fprintf(fp,"/n");/以上得出走時(shí)并寫入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("每個(gè)單元格的距離.txt","w"), *fp1=fopen("時(shí)間.txt","w"); double x_jf4=0.0,y_jf4;/左側(cè)激發(fā)點(diǎn)坐標(biāo) double x_js4=0.0,y_js4;/右側(cè)激發(fā)點(diǎn)的坐標(biāo) double X_jf4,Y_begin4=0.0;/上下激發(fā)點(diǎn)坐標(biāo) double X_end4,Y_js4=0.0;/下頂下激發(fā)點(diǎn)坐標(biāo)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è)接收點(diǎn)坐標(biāo)double x_you12=0.0,y_you12;/右側(cè)接收點(diǎn)坐標(biāo)double x_shang9,y_shang9=0.0;/上側(cè)接收點(diǎn)坐標(biāo)double x_xia9,y_xia9=0.0;/下側(cè)接收點(diǎn)坐標(biāo)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ā)時(shí)截距double xl_zuo2;double b_zuo2;/左側(cè)激發(fā)時(shí)截距double xl_zuo3;double b_zuo3;/左側(cè)激發(fā)時(shí)截距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ā)時(shí)截距double xl_you2;double b_y

17、ou2;/右側(cè)激發(fā)時(shí)截距double xl_you3;double b_you3;/右側(cè)激發(fā)時(shí)截距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ā)時(shí)截距double xl_shang2;double b_shang2;/上側(cè)激發(fā)時(shí)截距double xl_shang3;double b_shang3;/上側(cè)激發(fā)時(shí)截距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ā)時(shí)截距double xl_xia2;double b_xia2;/下側(cè)激發(fā)時(shí)截距double xl_xia3;double b_xia3;/下側(cè)激發(fā)時(shí)截距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("走時(shí).txt","r"); fp10=fopen("不為零的個(gè)數(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ù)的個(gè)數(shù) fprintf(fp10,"%fn",Mnm); N=0; /迭代次數(shù)10000次 /計(jì)算誤差/ 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