成都理工大学层析成像实验报告
本科生实验报告
实验课程地球物理层析成像
学院名称地球物理学院
专业名称勘查技术与工程
学生姓名
学生学号
指导教师曹俊兴
实验地点5417
实验成绩
二〇一五年三月二〇一五年四月
地震走时层析成像实验
——《地球物理正反演概论》课程结业报告
学号：201205060423
姓名：马力衡
专业：勘查技术与工程
手机：180********
摘要
运用C语言程序，正演得到地震走时和射线在传播过程中经过离散化处理单元格内的距离。通过反演程序反演出地下异常速度值，将反演所得速度值成图与原始速度成图进行比较，得出结论。离散化处理模型建立，单边激发,四边激发直射线正演，单边激发，四边激发反演异常值，用代数重建算法迭代慢度矩阵对单边，四边激发进行迭代。
关键词：离散化反演迭代
第1章地震走时层析成像实验1.1实验内容
1.1.1直射线正演：
使用直射线追踪方法计算走时的正演；分块均匀模型
1.1.1.1单边激发正演
程序：
#include
#include
void main()
{
int v[12][9];
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",&v[i][j]);
}
double b[12]; //截距;
double xl[12][12]; //斜率;
double y_jf[12],y_js[12];//激发点与接收点的纵坐标for(i=0;i<12;i++)
{
y_jf[i]=1.5+3.0*i;//激发点点坐标的方程
for(j=0;j<12;j++)
xl[i][j]=(y_js[j]-y_jf[i])/(45-0);//斜率
printf("%f\n",xl[i][j]);
}
}
for(i=0;i<12;i++)
{
b[i]=1.5+i*3;//每条射线截距
}
//以上在求射线的斜率和射线在纵轴上的截距 //
double ft_t=0.0; //每一格的时间; double fl[12][12][12][9]; //每一格射线的长度; double Time[12][12]; //每条射线的时间; double 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++)
{
Time[i][j]=0.0;
for(n=0;n<12;n++)
{
for(m=0;m<9;m++)
{
fl[i][j][n][m]=0;
x0=5*m;
x1=x0+5;
y0=3*n;
y1=y0+3;
y_0=xl[i][j]*x0+b[i];
if(y_0<=y1&&y_0>=y0)
{
X0=x0;
Y0=y_0;
y_1=xl[i][j]*x1+b[i];
if((y_1<=y1)&&(y_1>=y0))
{
Y1=y_1;
X1=x1;
{
x_1=(y0-b[i])/xl[i][j];
if(x_1<=x1&&x_1>=x0)
{
Y1=y0;
X1=x_1;
}
else
{
x_1=(y1-b[i])/xl[i][j];
X1=x_1;
Y1=y1;
}
}
fl[i][j][n][m]=sqrt((X1-X0)*(X1-X0)+(Y1-Y0)*(Y1-Y0));
}
else
{
y_1=xl[i][j]*x1+b[i];
if(y_1<=y1&&y_1>=y0)
{
X1=x1;
Y1=y_1;
x_0=(y0-b[i])/xl[i][j];
if(x_0<=x1&&x_0>=x0)
{
X0=x_0;
Y0=y0;
}
else
{
x_0=(y1-b[i])/xl[i][j];
X0=x_0;
Y0=y1;
}
fl[i][j][n][m]=sqrt((X1-X0)*(X1-X0)+(Y1-Y0)*(Y1-Y0));
}
else
if(x_0<=x1&&x_0>=x0)
{
X0=x_0;
Y0=y0;
x_1=(y1-b[i])/xl[i][j];
X1=x_1;
Y1=y1;
fl[i][j][n][m]=sqrt((X1-X0)*(X1-X0)+(Y1-Y0)*(Y1-Y0));
}
else fl[i][j][n][m]=0.0;
}
}
fprintf(fp_ds,"%f ",fl[i][j][n][m]);
Time[i][j]+=fl[i][j][n][m]/v[n][m];//每个单元格的走时；
printf("%d ",v[n][m]);//射线传播的总时间
}
}
fprintf(fp_ds,"\n");
}
}
//以上判定射线存在的单元格情况//
FILE *fp1;
fp1=fopen("走时.txt","w");
for(i=0;i<12;i++)
{
for(j=0;j<12;j++)
{
fprintf(fp1,"%f\n",Time[i][j]);
}
//fprintf(fp,"/n");
}
//以上得出走时并写入txt文件中输出//
1.1.1.2四边激发正演
程序；
#include
#include
#include
double funfpo1(double xl,double b,FILE *fp,FILE *fp1)
double fl[12][9];
double X0,Y0;
double X1,Y1;
double x0,x1,y0,y1;
double x_0,x_1,y_0,y_1;
double v[12][9],Time=0.0; for(i=0;i<12;i++)
{
for(j=0;j<9;j++)
{
v[i][j]=3.0;
v[2][2]=5.0;
v[3][2]=5.0;
v[8][5]=2.0;
v[8][6]=2.0;
}
}
for(n=0;n<12;n++)
{
y0=3*n;y1=y0+3;
for(m=0;m<9;m++)
{
fl[n][m]=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;
if((y_1<=y1)&&(y_1>=y0))
{
Y1=y_1;X1=x1;
}
else
{
x_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;
}
}
fl[n][m]=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; Y0=y1;
}
fl[n][m]=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;
fl[n][m]=sqrt((X1-X0)*(X1-X0)+(Y1-Y0)*(Y1-Y0));
}
}
}
fprintf(fp,"%f ",fl[n][m]);
Time+=fl[n][m]/v[n][m];
}
fprintf(fp,"\n");
}
fprintf(fp1,"%f\n ",Time);
return 0;
}
void main()
{
int i,j;
FILE *fp=fopen("每个单元格的距离.txt","w"),
*fp1=fopen("时间.txt","w");
double x_jf[4]={0.0},y_jf[4];//左侧激发点坐标
double x_js[4]={0.0},y_js[4];//右侧激发点的坐标 double X_jf[4],Y_begin[4]={0.0};//上下激发点坐标 double X_end[4],Y_js[4]={0.0};//下顶下激发点坐标
for(i=0;i<4;i++)
{
for(j=0;j<4;j++)
{
y_jf[i]=7.5+6*i;
x_js[j]=45.0;
y_js[j]=7.5+6*j;
X_jf[i]=7.5+10.0*i;
X_end[j]=7.5+10.0*j;
Y_js[j]=36.0;
}
}
double x_zuo[12]={0.0},y_zuo[12];//左侧接收点坐标double x_you[12]={0.0},y_you[12];//右侧接收点坐标double x_shang[9],y_shang[9]={0.0};//上侧接收点坐标double x_xia[9],y_xia[9]={0.0};//下侧接收点坐标
for(i=0;i<12;i++)
{
for(j=0;j<9;j++)
{
x_you[i]=45.0;
y_zuo[i]=1.5+3*i;
y_xia[j]=36.0;
x_shang[j]=2.5+5*j;
x_xia[j]=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_shang[j]-y_jf[i])/(x_shang[j]-x_jf[i]);
b_zuo1=10.5+6*i;
funfpo1( xl_zuo1,b_zuo1,fp,fp1);
xl_zuo3=(y_xia[j]-y_jf[i])/(x_xia[j]-x_jf[i]);
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_you[j]-y_jf[i])/(x_you[j]-x_jf[i]);
funfpo1( xl_zuo2,b_zuo2,fp,fp1);
}
}
//右侧激发//
double xl_you1;
double b_you1;//右侧激发时截距
double xl_you2;
double b_you2;//右侧激发时截距
double xl_you3;
double b_you3;//右侧激发时截距
for(i=0;i<4;i++)
{
for(j=0;j<9;j++)
{
xl_you1=(y_shang[j]-y_js[i])/(x_shang[j]-x_js[i]);
b_you1=y_js[i]-xl_you1*x_js[i];
funfpo1(xl_you1,b_you1,fp,fp1);
xl_you3=(y_xia[j]-y_js[i])/(x_xia[j]-x_js[i]);
b_you3=y_js[i]-xl_you3*x_js[i];
funfpo1(xl_you3,b_you3,fp,fp1);
}
}
for(i=0;i<4;i++)
{
for(j=0;j<12;j++)
{
xl_you2=(y_zuo[j]-y_js[i])/(x_zuo[j]-x_js[i]);
b_you2=y_js[i]-xl_you2*x_js[i];
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=(y_zuo[j]-Y_begin[i])/(x_zuo[j]-X_jf[i]);
b_shang1=Y_begin[i]-xl_shang1*X_jf[i];
funfpo1(xl_shang1,b_shang1,fp,fp1);
xl_shang2=(y_you[j]-Y_begin[i])/(x_you[j]-X_jf[i]);
b_shang2=Y_begin[i]-xl_shang2*X_jf[i];
funfpo1(xl_shang2,b_shang2,fp,fp1);
}
}
for(i=0;i<4;i++)
{
for(j=0;j<9;j++)
{
if(x_xia[j]==X_jf[i])
{
xl_shang3=INT_MAX;
b_shang3=Y_begin[i]-xl_shang3*X_jf[i];
}
{
xl_shang3=(y_xia[j]-Y_begin[i])/(x_xia[j]-X_jf[i]);
b_shang3=Y_begin[i]-xl_shang3*X_jf[i];
}
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;j++)
{
xl_xia1=(y_zuo[j]-Y_js[i])/(x_zuo[j]-X_end[i]);
b_xia1=Y_js[i]-xl_xia1*X_end[i];
funfpo1( xl_xia1, b_xia1,fp,fp1);
xl_xia2=(y_you[j]-Y_js[i])/(x_you[j]-X_end[i]);
b_xia2=Y_js[i]-xl_xia2*X_end[i];
funfpo1(xl_xia2,b_xia2,fp,fp1);
}
}
for(i=0;i<4;i++)
{
for(j=0;j<9;j++)
if(x_shang[j]==X_end[i])
{
xl_xia3=INT_MAX;
b_xia3=Y_js[i]-xl_xia3*X_end[i];
}
else
{
xl_xia3=(y_shang[j]-Y_js[i])/(x_shang[j]-X_end[i]);
b_xia3=Y_js[i]-xl_xia3*X_end[i];
}
funfpo1(xl_xia3,b_xia3,fp,fp1);
}
}
1.1.2 反演（矩阵方程求解）：
1.1.
2.1 单边激发反演
程序；
#include
#include
void main()
{
double r0, d0[12][12][12][9], a1,a2,Time12][12]，md[12][9];
int N，i,j,n,m;
double M[12][9];
double wucha[12][12]={0.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 ",&d0[i][j][n][m]);
}
/*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.2lf\n",d0[i][j][n][m]);
}
}
}*/
for(i=0;i<12;i++)
{
for(j=0;j<12;j++)
fscanf(fp8,"%lf",&Timei][j]);
}
/*for(i=0;i<12;i++)
{
for(j=0;j<12;j++)
{
printf("%5.6lf\n",Timei][j]);
}
}*/
for(n=0;n<12;n++)
{
for(m=0;m<9;m++)
{
M[n][m]=0.0;
md[n][m]=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(d0[i][j][n][m]!=0)
M[n][m]+=1; //算出系数矩阵中每一列不为零值的数的个数 }
}
}
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<12;n++)
{
for(m=0;m<9;m++)
r0+=d0[i][j][n][m]*md[n][m];
wucha[i][j]=(Timei][j]-r0); //printf("%f",wucha[i][j]);
}
}
}
//重建算法迭代慢度矩阵 //
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+=d0[i][j][m][n]*d0[i][j][m][n]; a2+=(wucha[i][j]*d0[i][j][m][n]);
}
}
md[m][n]+=a2/(a1*M[m][n]);
}
}
N++;
}
//以txt格式输出最后数据//
fp9=fopen("最终结果.txt","w");
for(i=0;i<12;i++)
{
fprintf(fp9,"%lf ",1.0/md[i][j]); }
fprintf(fp9,"\n");
}
}
成图；