地球物理反演与层析成像-结业论文与程序设计

********班 **** 学号:*********

（ 2 ) 反演程序 #include\#include\struct data {

float A[108]; float b; float r; }ray[144];

结构体 ray【i】表示 b 每条射线的走时

A 每条射线在每条单元格内的长度 r 反演每条射线计算出的误差

main()

{int m[108];

float x[108],R,a1,a2; x[i]反演的每个单元格的慢度 int n=0,j,i,j1;

FILE *fp1,*fp2,*fp3;/*,*fp3,*fp4*/

/*----------------------------------------------------------------*/读取正演的结果

fp1=fopen(\矩阵A2.txt\fp2=fopen(\走时.txt\

for(i=0;i<144;i++) for(j=0;j<108;j++)

fscanf(fp1,\

for(i=0;i<144;i++)

fscanf(fp2,\

for(i=0;i<108;i++)x[i]=0; 给第一次的慢度全部付零值 for(i=0;i<144;i++) for(j=0;j<108;j++)

printf(\

printf(\

for(j1=0;j1<108;j1++) 算出系数矩阵中每一列不为零值的数的个数 {m[j1]=0; x[j1]=0;

for(i=0;i<144;i++) if(ray[i].A[j1]!=0)m[j1]+=1; printf(\

********班 **** 学号:*********

}

n=20; 迭代次数20次 while(n) {

误差的计算

for(i=0;i<144;i++) { R=0;

for(j=0;j<108;j++) R+=ray[i].A[j]*x[j]; ray[i].r=ray[i].b-R; }

for(j1=0;j1<108;j1++) {

{ a2=0;

/*----------------------------------------------------------------*/ 利用代数重建算法迭代慢度矩阵

for(i=0;i<144;i++) {

a1=0.0;

for(j=0;j<108;j++) a1+=ray[i].A[j]*ray[i].A[j]; a2+=ray[i].r*ray[i].A[j1]/a1; }

x[j1]+=a2/m[j1]; }

//printf(\ 输出中间数据 }n--; }

// for(j=0;j<108;j++) //{ if(x[j]>2)x[j]=10; // if(x[j]<0)x[j]=0; // }

输出最后数据

fp3=fopen(\n=0;

for(i=0;i<12;i++){for(j=0;j<9;j++) {fprintf(fp3,\fprintf(fp3,\} }

********班 **** 学号:*********

可编译程序源代码及附带数据；

（3）实验结果，图及文字说明（至少包含五幅图）

原始数据454520次迭代4540200次迭代452000次迭代40354040353535303030302525252520202020151515151010101055101520253035551015202530355510152025303555101520253035通过上图正演和反演所得出的图可以看出迭代次数越多反演得出的数据和原始数据就越接近（如右上角的异常边缘可以看出2000次迭代的结果明显比二十次的要好，黑色的要少）但是无论迭代多少次图像的边缘都会产生很大的异常，这应该是sirt算法本身的缺陷导致，虽然如此，本程序的算法已经可以基本大致正确的反演出原始的慢度数据，并且异常值在迭代次数相对大的情况下比较明显。

20454540402004540200035353530303025252520202015151510101055101520253035551015202530355510152025上三图即是不同次数迭代与实际慢度做差值得到的误差，从图中更能看出边缘的异常与迭代的好坏比较。 给出求误差程序： #include \#include \main()

********班 **** 学号:*********

{

FILE *fp1,*fp2,*fp3,*fp4,*fp5,*fp6,*fp7; int i,j; float

in[12][9],in1[12][9],in2[12][9],in3[12][9],out1[12][9],out2[12][9],out3[12][9];

fp1=fopen(\慢度.txt\ fp2=fopen(\ fp3=fopen(\ fp4=fopen(\ fp5=fopen(\ fp6=fopen(\ fp7=fopen(\ for(i=0;i<12;i++) for(j=0;j<9;j++){ in[i][j]=0;

in1[i][j]=0; in2[i][j]=0; in3[i][j]=0; out1[i][j]=0; out2[i][j]=0; out3[i][j]=0;} for(i=0;i<12;i++)

for(j=0;j<9;j++)fscanf(fp1,\ for(i=0;i<12;i++)

for(j=0;j<9;j++)fscanf(fp2,\ for(i=0;i<12;i++)

for(j=0;j<9;j++)fscanf(fp3,\ for(i=0;i<12;i++)

for(j=0;j<9;j++)fscanf(fp4,\ for(i=0;i<12;i++)

for(j=0;j<12;j++)out1[i][j]=in[i][j]-in1[i][j]; for(i=0;i<12;i++)

for(j=0;j<12;j++)out2[i][j]=in[i][j]-in2[i][j]; for(i=0;i<12;i++)

for(j=0;j<12;j++)out3[i][j]=in[i][j]-in3[i][j]; for(i=0;i<12;i++) {

for(j=0;j<9;j++) {

fprintf(fp5,\ }

fprintf(fp5,\