📄 avo_2正演模型实验_加噪音_角度区域.cpp
字号:
//---------------------------------------------------------------//
// AVO 技术正演模型实验 //
// 3 层模型（一） //
// AVO 技术正演模型实验 //
// 40炮24道 //
// 地层倾角为0度 //
//---------------------------------------------------------------//
//---------------------------------------------------------------//
#include <stdio.h>
#include<math.h>
#define PI 3.1415926
#define nc 3 //地层层数
#define np 4 //地层参数个数
#define dt 0.002 //采样间隔
#define zbn 51 //子波长度
#define zf 30 //子波主频
#define zA 1.0 //子波最大振幅
//——————————观测系统参数————————————//
#define Nd 24 //迭加次数
#define dx 20 //道间距
#define ds 20 //炮间距//
#define x_0 40 //最小偏移距
#define L 48 //接收道数
#define Ns 40 //炮数
#define tl 300 //采样点数
#define ncdp 30 //cdp道数
void main()
{
int i,j,k;
float par[nc][np];
float x[Ns][L],zb[zbn];
float par0[np]={0.0};
FILE *fp1;
FILE *fp2;
FILE *fp3;
FILE *fp4;
FILE *fp5;
FILE *fp6;
FILE *fp7;
//==== 子程序说明 ====//
void fzb(float zb[zbn]);
void flljl(float x[Ns][L],float zb[zbn],float par[nc][np]);
printf("程序运行请稍后......\n\n");
//======================= 程序主体 =====================//
if((fp1=fopen("地层模型参数表2-1.txt","r"))==NULL)
printf("地层参数文件没有找到！\n");
else
printf("参数文件打开成功！\n");
fp3=fopen("含油气层参数.txt","r");
// for(i=0;i<np;i++)
//{
//fscanf(fp7,"%f\n",&par_0[i]);
fscanf(fp3,"%f\t%f\t%f\t%f\n",&par0[0],&par0[1],&par0[2],&par0[3]);
printf("%f\t%f\t%f\t%f\n",par0[0],par0[1],par0[2],par0[3]);
//}
for(i=0;i<nc;i++)
{
fscanf(fp1,"%f\t%f\t%f\t%f\n",&par[i][0],&par[i][1],&par[i][2],&par[i][3]);
printf("%f m/s\t%f m/s\t%f g/cm3\t%f m\n",par[i][0],par[i][1],par[i][2],par[i][3]);
if(i>=0)
printf("-R-----------------------------------------------------------R-\n");
}
//======== 子波 =====//
fzb(zb);
fp2=fopen("子波.txt","w+");
for(i=0;i<zbn;i++)
fprintf(fp2,"%f\t%f\n",dt*i,zb[i]);
//==== 计算各道反射系数，形成理论记录，抽道集 ====//
flljl(x,zb,par);
fclose(fp1);
fclose(fp2);
fclose(fp3);
}
//======================= 子波子程序 =====================//
void fzb(float zb[zbn])
{
int i,j;
float t;
for(i=0;i<=zbn/2;i++)
{
t=dt*i;
zb[zbn/2+i]=zA*cos(2*PI*zf*t)*exp(-0.5*pow(PI*zf*t,2));
zb[zbn/2-i]=zA*cos(2*PI*zf*t)*exp(-0.5*pow(PI*zf*t,2));
}
}
//======================= 计算各道反射系数，形成理论记录，抽道集子程序 =====================//
void flljl(float x[Ns][L],float zb[zbn],float par[nc][np])
{
int i,i1,i2,j,k,kn,r1,r2,l2,flag,flag2;
int j1;
float y[tl+zbn-1],yy[tl];
float t_0[nc-1]={0.0},t[nc-1],vr,sum=0.0;
float R1[nc-1],x1,sa2[nc-1],ta2,vp[2],vs[2],ro[2],K,da,a,A,B,C;//sina2->sa2和一些过渡参数
float R[tl]={0.0};
float par_0[np]={1896.0,50.0,1.50,100.0};//含油气层参数
int X,k2;
float sa2_j,sa_j,ti0,t00,vi;
float cdp_1[tl][ncdp]={0.0},cdp_2[tl][ncdp]={0.0},cdp_3[tl][ncdp]={0.0},tp1[tl]={0.0};//角度叠加约束
float cdp[tl][ncdp]={0.0},cdp0[tl][ncdp*Nd],tp[tl];
double zx[tl],zc=2045.0,zza=-0.06 ,zzb=0.06,zN[tl];
FILE *fp1;
FILE *fp2;
FILE *fp3;
FILE *fp4;
FILE *fp5;
FILE *fp6;
FILE *fp7;
FILE *fp8;
FILE *fp9;
//===== 计算基本参数 =====//
fp1=fopen("反射系数.txt","w+");
fp2=fopen("地震记录.dat","wb+");
fp3=fopen("地震记录(文本).dat","w+");
fp4=fopen("抽道得的道集.dat","wb+");
fp5=fopen("抽道叠加后30个cdp点.dat","wb+");
fp6=fopen("抽道得的道集（文本）.dat","w+");
fp7=fopen("含油气层参数.txt","r");
fp8=fopen("角度叠加(3组角度范围＋全角度叠加共4个30CDP).dat","wb+");
fp9=fopen("查_角度叠加(3组角度范围＋全角度叠加共4个30CDP).dat","w+");
//===== 读取含油气参数 =====//
for(i=0;i<np;i++)
{
fscanf(fp7,"%f\n",&par_0[i]);printf("%f\n",par_0[i]);
}
//fscanf(fp7,"%f\t%f\t%f\t%f ",&par_0[0],&par_0[1],&par_0[2],&par_0[3]);
//——————————噪音———————————//
zx[0]=8388607.0/8388608.0;
for(i=1;i<tl;i++)
zx[i]=zc*zx[i-1]-int(zc*zx[i-1]);
for(i=0;i<tl;i++)
zN[i]=float(zza+(zzb-zza)*zx[i]);
//——————————噪音———————————//
for(i=0;i<Ns;i++)
for(j=0;j<L;j++)
{
x1=x_0+dx*j;
X=ds*i+x1;
//if(X/dx>=30&&X/dx<40)//含油气参数替换
if((X+ds*i)/dx/2.0>=30&&(X+ds*i)/dx/2.0<32)
{
printf("炮：%d\n",i);
for(k2=0;k2<np;k2++)
{
par[1][k2]=par_0[k2];
// if(i==0)printf("****%d\t%f\n",j,par_0[k2]);
}
flag2=1;
}
else
flag2=0;
for(i1=0;i1<nc-1;i1++)
{
t_0[i1]=0.0;//注意清零！！
}
sum=0.0;
for(i1=0;i1<nc-1;i1++)
{
t[i1]=par[i1][3]/par[i1][0];
for(i2=0;i2<=i1;i2++)
t_0[i1]+=2*(t[i2]);
sum+=t[i1]*par[i1][0]*par[i1][0];
if(i==5&&j==0)printf("i=%d\tt[i1]=%f\tt_0[i1]=%f\n",i,t[i1],t_0[i1]);
}
vr=sqrt((sum*2)/t_0[nc-2]);
//if(i==0)printf("%d\t%f\n",j,t_0[nc-2]);
for(i1=0;i1<tl+zbn-1;i1++)//每道地震记录前归零！
y[i1]=0.0;
for(i1=0;i1<tl;i1++)
R[i1]=0.0;
for(k=0;k<nc-1;k++)
{
sa2[k]=pow(x1*par[k][0],2)/pow(vr,2)/(pow(vr*t_0[k],2)+pow(x1,2));
//sa2[k]=pow(par[k][0],2)*(x1-vr*t_0[k]*atan(x1/(vr*t_0[k])))/(x1*pow(vr,2));
//if(i==0&&k==1)printf("%d\t%f\n",j,vr);
if(sa2[k]<0)printf("%d\t,%d\t,%d\tsa2[]=%f\n",i,j,k,sa2[k]);
ta2=sa2[k]/(1-sa2[k]);
vp[0]=(par[k+1][0]+par[k][0])/2.0;
vp[1]=par[k+1][0]-par[k][0];
vs[0]=(par[k+1][1]+par[k][1])/2.0;
vs[1]=par[k+1][1]-par[k][1];
ro[0]=(par[k+1][2]+par[k][2])/2.0;
ro[1]=par[k+1][2]-par[k][2];
K=(pow(par[k+1][1],2)/pow(par[k+1][0],2)+pow(par[k][1],2)/pow(par[k][0],2));//K=vs2/vp2
//过渡参数计算↑↓……
A=(vp[1]/vp[0]+ro[1]/ro[0])/2.0;
B=vp[1]/(vp[0]*2.0)-4*K*vs[1]/vs[0]-2*K*ro[1]/ro[0];
C=vp[1]/(2*vp[0]);
R1[k]=A+B*sa2[k]+C*sa2[k]*ta2;
//printf("%d,%d,%d\n",i,j,k);
kn=t_0[k]/dt;
R[kn]=R1[k];
}
if(i==0)
{
fprintf(fp1,"%f\t%f\t%f\t%f\n",asin(sqrt(sa2[0]))*180/PI,R1[0],asin(sqrt(sa2[1]))*180/PI,R1[1]);
// printf("%d\t%f\t%f\tcha=%f\n",j,asin(sqrt(sa2[0]))*180/PI,R1[0],par[0][0]);
}
//查！！反射系数计算
for(i1=0;i1<tl;i1++)
for(j1=0;j1<zbn;j1++)
{
y[i1+j1]+=R[i1]*zb[j1];
}
for(i1=0;i1<tl;i1++)
{
yy[i1]=y[i1+zbn/2]+zN[i1];//;
fprintf(fp3,"%f\t",yy[i1]+0.005*i1);
}
fwrite(yy,4,tl,fp2);
fprintf(fp3,"\n");
//——— 抽道、迭加求和↓ ———//
for(r1=0;r1<ncdp;r1++)//选道是flag=1,否flag=0
{
flag=r1/2;
if(i>=flag)
if(j+L/Nd*i-L+L/Nd==r1)//?
{
l2=r1*Nd+i-flag;//-r1%Nd;//因该是除以叠加次数的余数！！
if(r1+1==1)printf("CDP道：i=%d\tj=%d\tcdp=%d\tl2=%d\n",i+1,j+1,r1+1,l2+1);
//if(r1==0)printf("X=%d\t%d\n",X,X/dx);
//if(flag2==1)printf("**&&&cdp=%d,X=%dpar=%f\n",r1+1,X,par[1][1]);
for(r2=0;r2<tl;r2++)
{
cdp0[r2][l2]=yy[r2];
//——— 角度叠加～～确定角度范围↓...———//
t00=r2*dt;
for(k=0;k<nc-1;k++)
{
if(t00<=t_0[k])
{
vi=par[k][0];
break;
}
else
vi=par[nc-2][0];
}
sa2_j=pow(vi,2)*(x1-vr*t00*atan(x1/(vr*t00)))/(x1*pow(vr,2));
sa_j=asin(sqrt(sa2_j))*180/PI;
//三组角度范围角度叠加//
if(sa_j>=0&&sa_j<5)
cdp_1[r2][r1]+=yy[r2];
else
cdp_1[r2][r1]+=0.0;
if(sa_j>=5&&sa_j<19)
cdp_2[r2][r1]+=yy[r2];
else
cdp_2[r2][r1]+=0.0;
if(sa_j>=19&&sa_j<35)
cdp_3[r2][r1]+=yy[r2];
else
cdp_3[r2][r1]+=0.0;
//——— 角度叠加～～确定角度范围↑...———//
cdp[r2][r1]+=yy[r2];
}
}
}
//——— 抽道、迭加求和↑ ———//
}//6炮24道循环结束
printf("X总测线长度：%d\t总炮数：%d\n",X,X/ds);
/*for(j=0;j<ncdp;j++)
for(i=0;i<tl;i++)
{
cdp[i][j]=cdp[i][j]/Nd;//振幅平均
cdp_1[i][j]=cdp[i][j]/Nd;
cdp_2[i][j]=cdp[i][j]/Nd;
cdp_3[i][j]=cdp[i][j]/Nd;
}*/
/* for(j=0;j<ncdp*Nd;j++)
{
for(i=0;i<tl;i++)
tp[i]=cdp0[i][j];
fwrite(tp,4,tl,fp4);
}*/
for(j=0;j<ncdp*Nd;j++)
{
for(i=0;i<tl;i++)
{
tp[i]=cdp0[i][j];
fprintf(fp6,"%f\t",tp[i]+0.005*i);
}
fprintf(fp6,"\n");
fwrite(tp,4,tl,fp4);
}
for(j=0;j<ncdp;j++)
{
for(i=0;i<tl;i++)
{
tp[i]=cdp[i][j];
}
fwrite(tp,4,tl,fp5);
}
//角度叠加写入数据
flag=1;
while(flag<5)
{
for(j=0;j<ncdp;j++)
{
for(i=0;i<tl;i++)
{
if(flag==1) { tp[i]=cdp_1[i][j];fprintf(fp9,"%f\t",tp[i]+0.005*i);}
if(flag==2) tp[i]=cdp_2[i][j];
if(flag==3) tp[i]=cdp_3[i][j];
if(flag==4) tp[i]=cdp[i][j];
}
if(flag==1)fprintf(fp9,"\n");
fwrite(tp,4,tl,fp8);
}
fwrite(tp1,4,tl,fp8);
printf("%d\n",flag);
flag++;
}
fclose(fp1);
fclose(fp2);
fclose(fp3);
fclose(fp4);
fclose(fp5);
fclose(fp6);
fclose(fp7);
fclose(fp8);
fclose(fp9);
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -