地震层析成像的正演实验报告.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]）;
升级会员 