地震层析成像的正演实验报告.docx
- 文档编号:3516784
- 上传时间:2022-11-23
- 格式:DOCX
- 页数:41
- 大小:678.95KB
地震层析成像的正演实验报告.docx
《地震层析成像的正演实验报告.docx》由会员分享,可在线阅读,更多相关《地震层析成像的正演实验报告.docx(41页珍藏版)》请在冰豆网上搜索。
地震层析成像的正演实验报告
1-1-1正演的速度模型图
图1-1-2分块均匀的模型
1-2正演的后的走时图
图1-3反演前后速度对比图
图1-5-a第0炮,第5接收点的数据
图1-5-b-1正演第1炮,第8接收点(0为开始的激发点,0开始的接收点)
图1-5-b-2与1-5-b-2对应的验证图形(附注:
由于本人u盘被病毒入侵,导致本人做得CAD图丢失,此图引用廖松杰同学的CAD图像,但是1-5-b-2由本人程序自己得出,特此说明。
)
图1-5-c四边放炮,四边接收左方第2激发,图1-5-b单边接收第0炮,
第25接收的r图像第5接收点的r数据图
正反演的程序
单边放炮单边接收:
#include
#include
#include
voidfun1(intn,doubleR[144][108],doublet[12][12]);
voidfun2(doublek,doubleo,doublet[12][12],doubleR[144][108],intm,intn);//k为斜率,o为炮点坐标,相当于截距;j;
doublefun(doublex1,doubley1,doublex2,doubley2);
voidmain()
{
FILE*fp;
inti,j,m,n,l,f;
doublec[12],d[12],K[12][12],r[12][9]={0},v[12][9]={0.0},t[12][12]={0.0},u,o,w,R[144][108]={0.0},k,v2[144]={0};
floatv1;
//*******************************************************************************************************//
for(i=0;i<12;i++)//第i行
{
for(j=0;j<9;j++)//第j列
{v[i][j]=3000;
}
}
v[2][2]=5000.0;
v[3][2]=5000.0;
v[8][5]=2000.0;
v[8][6]=2000.0;
for(i=0;i<12;i++)//第i行
{
for(j=0;j<9;j++)//第j列
{
v2[j*12+i]=v[i][j];
}
}
fp=fopen("速度","wb");
for(i=0;i<12;i++)
{
for(j=0;j<9;j++)
{
v1=v[i][j];
fwrite(&v1,sizeof(float),1,fp);
}
}fclose(fp);
//*******************************************************************************************************//
//*************************************计算各点的斜率***************************************************//
for(i=0;i<12;i++)
{
//printf("第%d炮的斜率\n",i);
c[i]=(i+0.5)*5.0;//***************************************************************激发点********//
for(j=0;j<12;j++)
{
d[j]=(j+0.5)*5;//***********************************************************接收点********//
K[i][j]=(d[j]-c[i])/(9.0*3);//**********************************************K斜率*********//
//printf("K[%d][%d]=%f\n",i,j,K[i][j]);
}//printf("\n");
}
//**********************************************************************************************************//
for(i=0;i<12;i++)//第i炮
{
for(j=0;j<12;j++)//第j接收点
{
if(K[i][j]==0)//平行于x轴,该行所在的每一个网格均经过,路程都是3
{
fun1(i,R,t);
printf("t[%d][%d]=%f\n",i,j,t[i][j]);
}
elseif(K[i][j]!
=0)
{
k=K[i][j];
o=c[i];
m=i;n=j;
fun2(k,o,t,R,i,j);
printf("t[%d][%d]=%lf\n",i,j,t[i][j]);
}
}
}
fp=fopen("time.txt","w");
for(i=0;i<12;i++)
{
for(j=0;j<12;j++)
{
fprintf(fp,"%lf\n",t[i][j]);
}
}fclose(fp);
fp=fopen("系数矩阵R的值.txt","w");
for(i=0;i<144;i++)
{
for(j=0;j<108;j++)
{
fprintf(fp,"%f\t",R[i][j]);
}fprintf(fp,"\n");
}fclose(fp);
fp=fopen("原来的速度值.txt","w");
for(j=0;j<108;j++)
{
fprintf(fp,"%f\t",v2[j]);
}fclose(fp);
}
//**********************************************************************************************************//
//**********************************************************************************************************//
//*******************************当斜率k为0的时候,计算走时t的值********************************************//
//**********************************************************************************************************//
voidfun1(intn,doubleR[144][108],doublet[12][12])
{
FILE*fp1;
doubleb=0.0;
inti=0,j=0,q=0;//循环变量
doubler[12][9]={0.0},v[12][9];
//***************************************************************
for(i=0;i<12;i++)//第i行
{
for(j=0;j<9;j++)//第j列
{v[i][j]=3000;
}
}
v[2][2]=5000.0;
v[3][2]=5000.0;
v[8][5]=2000.0;
v[8][6]=2000.0;
for(j=0;j<9;j++)
{
r[n][j]=3.0;
}
//*********************************************************************//
//**********************写出检验r的值**********************************//
/*fp1=fopen("r的值.txt","w");
for(i=0;i<12;i++)
{
for(j=0;j<9;j++)
{
fprintf(fp1,"%f\t",R[i][j]);
}fprintf(fp1,"\n");
}fclose(fp);*/
//***********************************************************************//
//***********************************************************************//
for(i=0;i<12;i++)//第i行
{
for(j=0;j<9;j++)//第j列
{
b+=r[i][j]*(1/v[i][j]);
}
}
t[n][n]=b;
for(i=0;i<12;i++)//第i行
{
for(j=0;j<9;j++)//第j列
{
R[n*12+n][q++]=r[i][j];
}
}
}
doublefun(doublex1,doubley1,doublex2,doubley2)
{
doubles;
{s=(y2-y1)*(y2-y1)+(x2-x1)*(x2-x1);
returnsqrt(s);
}
}
//*****************************************************//
//**********************************************************************************************************//
voidfun2(doublek,doubleo,doublet[12][12],doubleR[144][108],intm,intn)//k为斜率,o为炮点坐标,相当于截距;
{
FILE*fp2;
inti=0,j=0,q=0;;//循环变量
intw1,w2,w3,w4;//中间变量,用来判断点在分块均匀上的位置
doublep=0,v[12][9]={0.0},r[12][9]={0.0};
doublex1,y1,x2,y2,x3,y3,x4,y4;
floatr1;
//***************************************************************
v[2][2]=5000.0;
v[3][2]=5000.0;
v[8][5]=2000.0;
v[8][6]=2000.0;
for(i=0;i<12;i++)//第i行
{
for(j=0;j<9;j++)//第j列
{v[i][j]=3000;
r[i][j]=0;
}
}
for(i=0;i<12;i++)
{
for(j=0;j<9;j++)
{
y1=i*5.0;
x1=(y1-o)/k;//计算交点1,由y计算x,交点一位于上边
y2=(i+1)*5.0;
x2=(y2-o)/k;//计算交点1,由y计算x,交点2位于下边
x3=j*3.0;
y3=k*x3+o;//计算交点3,由x计算y,交点3位于左边
x4=(j+1)*3.0;
y4=k*x4+o;//计算交点3,由x计算y,交点四位于右边
//**********************************************************************
//***判断射线是否经过分块均匀的网格点上,四个交点是否在网格的四条边上*****
//**********************************************************************
//***注意:
网格的上下两条边y值相等,网格的左右两边x的值相等***************
//**********************************************************************
w1=(x1>=(j*3.0))&&(x1<=((j+1)*3.0))&&(y1==(i*5.0));//上方
w2=(x2>=(j*3.0))&&(x2<=((j+1)*3.0))&&(y2==((i+1)*5.0));//下方
w3=(y3>=(i*5.0))&&(y3<=((i+1)*5.0))&&(x3==(j*3.0));//左方
w4=(y4>=(i*5.0))&&(y4<=((i+1)*5.0))&&(x4==((j+1)*3.0));//右方
//**********************************************************************************
//计算路径长度r,当有两个点存在时,有下面的六种情况。
if(w1!
=0&w3!
=0)
{
r[i][j]=fun(x1,y1,x3,y3);
}
elseif(w1!
=0&&w2!
=0)
{
r[i][j]=fun(x1,y1,x2,y2);
}
elseif(w2!
=0&&w4!
=0)
{
r[i][j]=fun(x2,y2,x4,y4);
}
elseif(w4!
=0&&w1!
=0)
{
r[i][j]=fun(x1,y1,x4,y4);
}
elseif(w3!
=0&&w2!
=0)
{
r[i][j]=fun(x2,y2,x3,y3);
}
elseif(w3==1&&w4==1)
{
r[i][j]=fun(x3,y3,x4,y4);
}
//**************************************************
//走时等于射线经度每一个单元格的时间与慢速和(1/v)的累加
p+=(r[i][j]*(1/v[i][j]));
}
}
t[m][n]=p;
for(i=0;i<12;i++)
{
for(j=0;j<9;j++)
{
R[m*12+n][q++]=r[i][j];
}
}
//*******************************************************************************************************//
//***********************************检验第0炮,第五接收点r的正确性*******************************************************//
if(m==0&&n==5)
{
fp2=fopen("r2的值.txt","w");
for(i=0;i<12;i++)
{
for(j=0;j<9;j++)
{
fprintf(fp2,"%f\t",r[i][j]);
}fprintf(fp2,"\n");
}fclose(fp2);
fp2=fopen("r2的值","wb");
for(i=0;i<12;i++)
{
for(j=0;j<9;j++)
{
r1=r[i][j];
fwrite(&r1,sizeof(float),1,fp2);
}
}fclose(fp2);
}
//******************************************************************************************************//}
四边放炮四边接收:
#include
#include
#include
voidfun1(intn,doubleR[120][108],doublet[4][30],intm);
voidfun2(doublek,doubleo,doublet[4][30],doubleR[120][108],intm,intn);//k为斜率,o为炮点坐标,相当于截距;
doublefun(doublex1,doubley1,doublex2,doubley2);
voidfun3(intn,doubleR1[133][108],doublet[4][33],intm);
voidfun4(doublek,doubleo,doublet1[4][33],doubleR1[133][108],intm,intn);//k为斜率,o为炮点坐标,相当于截距
main()
{
FILE*fp;
inti,j,j1,m,n,l,e;
doublec1[12],c2[9],c3[12],c4[9],c[12],c5[9];
doubled1[12],d2[9],d3[12],d4[9];
doubleK1[4][30]={0.0},K2[4][33]={0.0},K3[4][30]={0.0},K4[4][33]={0.0};
doubleK[12][42],r[12][9]={0},v[12][9],t[4][30]={0.0},t1[4][33]={0.0},u,o,w,R[120][108]={0.0},R1[133][108]={0.0},k;
//*******************************************************************************************************//
//*******************************************************************************************************//
//*************************************计算各点的斜率***************************************************//
//************************************************************左方激发***********************************//
for(i=0;i<4;i++)
{
c1[i]=12.5+10*i;//***************************************************************左激发点********//
for(j1=0;j1<30;j1++)
{
if(j1<9)
{
d2[j1]=1.5+j1*3.0;//***********************************************************上接收点********//
K1[i][j1]=-c1[i]/d2[j1];
}
elseif(j1<21&&j1>=9)
{
d3[j1-9]=2.5+(j1-9)*5.0;//***********************右接收点**************************************//
K1[i][j1]=(d3[j1-9]-c1[i])/27.0;
}
else
{
d4[j1-21]=1.5+(j1-21)*3.0;//***********************下接收点*************************************//
K1[i][j1]=(60-c1[i])/d3[j1-21];
}
//printf("K1[%d][%d]=%f\n",i,j1,K1[i][j1]);
}
}
//***********************************************************上方激发*************************************
for(i=0;i<4;i++)
{
//printf("第%d炮的斜率\n",i);
c2[i]=4.5+6.0*i;//***************************************************************上激发点********//
for(j1=0;j1<33;j1++)
{
if(j1<12)
{
d3[j1]=2.5+(j1)*5.0;//***********************************************************右接收点********//
K2[i][j1]=d3[j1]/(27-c2[i]);
}
elseif(j1<21&&j1>=12)
{
d4[j1-12]=1.5+(j1-12)*3.0;//***********************下接收点**************************************//
K2[i][j1]=60/(d4[j1-12]-c2[i]);
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 地震 层析 成像 实验 报告