速度分析_二进制文件输出.cpp

AVO正演模型实验

//---------------------------------------------------------------//
//                 地震速度分析（二层模型）                      //
//---------------------------------------------------------------//
#include <stdio.h>
#include<math.h>
#define  PI   3.1415926
#define  M    2.0        //最大波峰与最大波谷绝对值比值
#define  dt   0.002      //采样间隔 
#define  t1_0 0.74       //第一层在第一炮回声时间
#define  t2_0 1.1        //第二层在第一炮回声时间
#define  tl   700        //采样点数
#define  ln   31         //雷克子波长度
#define  A    1.0       //雷克子波“衰减因子”
//——————————观测系统参数————————————//
#define  Nd   6          //迭加次数
#define  dx   80         //道间距
#define  x_0  160        //最小偏移距
#define  L    24         //接收道数
#define  ds   160        //炮间距
#define  Ns   6          //炮数
#define  a1  10          //测线和倾向的夹角
#define  ks   5          //真倾角
//——————————滤波因子参数————————————//
#define  hn   50         //滤波因子半长度
#define  hf   38         //滤波因子主频
#define  hdf  15         //滤波因子半门宽

#define  ncdp  4          //cdp道数
#define  nvc   4          //速度个数
#define  nv    5          //速度曲线条数
#define  daoz  2.0        //绘图的道间隔

void main()
{
	int    i,j,k,l,m,n,flag,r1,r2,l1,l2=0;
	float  tc1[Ns][L],tc2[Ns][L],lz[ln][2],N[tl],tp[tl],tp1[tl]={0.0};
	float  xxx[tl],x[tl],xx[tl][Ns*L],v[nvc][nv];
	float  f_0[2]={40.0,35.0},t0[nvc]={0.5,0.74,1.1,1.5};
	float  fi,dv=50.0,x0,v0[nvc]={2400.0,2600.0,3000.0,3600.0};  
	float  cdp0[tl][ncdp*Nd*(nv+1)];
	float  vv[nvc],cdp[tl][(ncdp+1)*nvc-1]={0.0},Cdp[tl][(ncdp+1)*nvc-1],cdp1[tl][ncdp]={0.0};//,v2[nvc],v1[nvc],v1[nvc],v1[nvc];
//==================  子程序说明  =========================//
	void flz(float f_0[2],float lz[ln][2]);
	void fzy(float N[]);
	void fv(float v0[],float v[][nv],float dv);
	void ftc(float tc[Ns][L],float fi,float v0,float t0);
	void flljl(float xxx[tl],float tc_1,float tc_2,float N[tl],float lz[ln][2]);
	void fxdjz(float xxx[tl],float x[tl],float vv[],float fi,
		       float x0,float tc1[][L],float tc2[][L],int i1,int j1);
	void flb (float cdp[tl][(ncdp+1)*nvc-1],int l1,float Cdp[tl][(ncdp+1)*nvc-1]);
//==================  子程序说明  =========================//
	FILE *fp1;
	FILE *fp2;
	FILE *fp3;
	FILE *fp4;
	FILE *fp5;
    FILE *fp6;
	FILE *fp7;
	FILE *fp8;
	FILE *fp9;
	fi=(float)asin(cos(a1*PI/180.0)*sin(ks*PI/180.0));

	printf("程序运行中........");
//——————————计算雷克子波及噪音————————————//
	fp1=fopen("雷克子波.txt","w+");
	fprintf(fp1,"第一层雷克子波\t第二层雷克子波\n");
    flz(f_0,lz);
	for(i=0;i<ln;i++)
	fprintf(fp1,"%f\t%f\n",lz[i][0],lz[i][1]);
	
	fzy(N);
	fp2=fopen("噪音.txt","w+");
	fprintf(fp2,"噪音\n");
	for(i=0;i<tl;i++)
	fprintf(fp2,"%f\n",N[i]);

//———————  求取共反射道记录道的层反射时间  ———————//
	fv(v0,v,dv);
	fp3=fopen("速度曲线.txt","w+");
	fprintf(fp3,"t0\tv1\tv2\tv3\tv4\tv5\n");
	for(i=0;i<nvc;i++)
	fprintf(fp3,"%f\t%f\t%f\t%f\t%f\t%f\n",t0[i],v[i][0],v[i][1],v[i][2],v[i][3],v[i][4]);
	ftc(tc1,fi,v0[1],t0[1]);//第一层
	ftc(tc2,fi,v0[2],t0[2]);//第二层
    fp4=fopen("第一层反射走时.dat","w+");
    fp5=fopen("第二层反射走时.dat","w+");
	for(i=0;i<Ns;i++)
	{
		for(j=0;j<L;j++)
			fprintf(fp4,"%f\t",tc1[i][j]/dt);
		fprintf(fp4,"\n");
	}
	for(i=0;i<Ns;i++)
	{
		for(j=0;j<L;j++)
			fprintf(fp5,"%f\t",tc2[i][j]/dt);
		fprintf(fp5,"\n");
	}
//———————   形成理论记录  ———————//	
	fp6=fopen("地震理论记录.dat","wb+");
	l=0;
	for(i=0;i<Ns;i++)
		for(j=0;j<L;j++)
		{	
			flljl(xxx,tc1[i][j],tc2[i][j],N,lz);
			for(k=0;k<tl;k++)
				xx[k][l]=xxx[k];
			l++;			
		}
    for(j=0;j<Ns*L;j++)
	{
		for(i=0;i<tl;i++)
			tp[i]=xx[i][j];
		fwrite(tp,4,tl,fp6);	
	}
//—————————————   动校正  —————————————//
	l1=0;
	flag=0;
	for(n=0;n<nv;n++)   
	{
		for(m=0;m<nvc;m++)			
		{
			vv[m]=v[m][n];//printf("%d\t%f\n",n,vv[m]);
		}
	    l=0;
		for(r1=0;r1<ncdp;r1++)
			for(r2=0;r2<tl;r2++)
				cdp1[r2][r1]=0.0;
//—————————————  抽道、迭加求和  ——————————//
	    while(l<Ns*L)
		{
			for(k=0;k<tl;k++)
				xxx[k]=xx[k][l];
		    i=l/L;//
		    j=l%L;
		    x0=float(j*dx+x_0);
			fxdjz(xxx,x,vv,fi,x0,tc1,tc2,i,j);	
				for(r1=0;r1<ncdp;r1++)//选道是flag=1,否flag=0
				{
					if(j+L/Nd*i-L+L/Nd==r1)
					{
						
						flag=1;
						l2=r1*(nv+1)*Nd+i;
					    for(r2=0;r2<tl;r2++) 
						{
							cdp1[r2][r1]+=x[r2];
							cdp0[r2][l2+(n+1)*Nd]=x[r2];	
							if(n==0)cdp0[r2][l2]=xxx[r2];
						}
					}
				    else 
						flag=0;
				}
			l++;
		}
			for(j=0;j<ncdp;j++) 
			{
				for(i=0;i<tl;i++) 
					cdp[i][l1]=cdp1[i][j]/Nd;//振幅平均
				l1++;
			}
	}
//—————————————  迭加前  ————————————— //
	fp7=fopen("动校正前后1-4cdp地震道记录.dat","wb+");
	for(j=0;j<ncdp*Nd*(nv+1);j++)
	{
		for(i=0;i<tl;i++)
			tp[i]=cdp0[i][j];
		fwrite(tp,4,tl,fp7);
	}
	fp8=fopen("迭加地震道记录.dat","wb+");
	flag=-1;
	l=0;
	for(j=0;j<(ncdp+1)*nv;j++)
	{
		for(i=0;i<tl;i++)
			tp[i]=cdp[i][l];
		if(j-flag==ncdp+1&&j!=0)
		{
			fwrite(tp1,4,tl,fp8);
		    flag=j;
		}
		else 
		{
			fwrite(tp,4,tl,fp8);
			l++;
		}	
	}

//—————————————  滤波  —————————————//
	flb(cdp,l1,Cdp);
	
	fp9=fopen("滤波后cdp道集.dat","wb+");
	l=0;	
	flag=0;
	for(j=0;j<(ncdp+1)*nv;j++)
	{
		for(i=0;i<tl;i++)
			tp[i]=Cdp[i][l];
	    if(j-flag==ncdp+1&&j!=0)
		{
			fwrite(tp1,4,tl,fp9);
		    flag=j;
		}
		else 
		{
			fwrite(tp,4,tl,fp9);
			l++;
		}	
	}
	fclose(fp1);
	fclose(fp2);
	fclose(fp3);
	fclose(fp4);
	fclose(fp5);
	fclose(fp6);
	fclose(fp7);
	fclose(fp8);
	fclose(fp9);
//………………………程序运行结束…………………………………//
	printf("\n=======   程序运行完毕！ =======\n");
}
//========================================================//
//=============         子程序         ===================//
//========================================================//
//——————————雷克子波子程序———————————//
void flz(float f_0[],float lz[][2])
{
	int i,j;
	float la;
	for(j=0;j<2;j++)
	{
		la=float(2/3.0*pow(f_0[j],2)*A*log10(M));
		for(i=0;i<ln;i++)
		lz[i][j]=float(sin(2*PI*f_0[j]*i*dt)*exp(-la*i*dt*i*dt));	
	}
}
//——————————噪音子程序———————————//
void fzy(float N[])
{
	int i;
	double x[tl],c=2045.0,a=-0.5,b=0.5;
	x[0]=8388607.0/8388608.0;
	for(i=1;i<tl;i++)
		x[i]=c*x[i-1]-int(c*x[i-1]);
    for(i=0;i<tl;i++)
		N[i]=float(a+(b-a)*x[i]);
}
//——————————速度曲线子程序———————————//
void fv(float v0[],float v[][nv],float dv)
{
	int i,j;
	for(j=0;j<nvc;j++)//4
		for(i=0;i<nv;i++)//5
		{
			v[j][i]=v0[j]+(i-nv/2)*dv;//5条速度曲线
		}
}
//————求取公反射道记录道的层反射时间子程序————//
void ftc(float tc[Ns][L],float fi,float v0,float t0)
{
	int i,j;
	float h0,tom,x;
	for(i=0;i<Ns;i++)
		for(j=0;j<L;j++)
		{
			x=float(j*dx+x_0);
			h0=float(v0*t0/2+(ds*i+x/2)*sin(fi));
			tom=2*h0/v0;
			tc[i][j]=float(sqrt(pow(tom,2)+pow(x*cos(fi)/v0,2)));
		}
}
//—————————— 形成理论记录子程序   ——————————//
void flljl(float xxx[tl],float tc_1,float tc_2,float N[tl],float lz[ln][2])
{
	int i,t1,t2,k=1;
	float sum[tl]={0};
	t1=(int)(tc_1/dt);
	t2=(int)(tc_2/dt);
	for(i=0;i<ln;i++)
	{
		if(i+t1<tl)		sum[i+t1]=sum[i+t1]+lz[i][0];
		if(i+t2<tl)     sum[i+t2]=sum[i+t2]+lz[i][1];
	}
	for(i=0;i<tl;i++)
		xxx[i]=sum[i]+N[i];	
}
//————————————   动校正子程序    ———————————//
void fxdjz(float xxx[tl],float x[tl],float vv[],float fi,float x0,float tc1[][L],float tc2[][L],int i1,int j1)
{
	int i,k;
	float t_0,v,ddt;
	for(i=0;i<tl;i++)
	{
		t_0=float(i*dt);
		if(t_0<=tc1[i1][j1-1])                       v=vv[1];
		if(t_0>tc1[i1][j1-1]&&t_0<=tc2[i1][j1-1])    v=vv[2];
		if(t_0>tc2[i1][j1-1])                        v=vv[3];
		ddt=float(sqrt(pow(t_0,2)+pow(x0*cos(fi)/v,2))-t_0);
		if((t_0+ddt)<=tl*dt)
		{
			k=(int)((t_0+ddt)/dt);
			x[i]=xxx[k];
		}
		else 
			x[i]=0.0;
	}
}
//————————————   滤波子程序    ———————————//
void flb (float cdp[tl][(ncdp+1)*nvc-1],int l1,float Cdp[][(ncdp+1)*nvc-1])
{
   int   i,j,k,l;
   float h[hn],y[tl+hn]={0.0};
   FILE *fp;
//......................计算滤波因子(半)...........//
   fp=fopen("滤波因子.dat","w+");
   h[0]=1.0;
   for(i=1;i<hn;i++) 
   {
	   h[i]=(float)(sin(2*PI*hdf*i*dt)*cos(2*PI*hf*i*dt)*cos((i*PI)/(2.0*hn))/(2*PI*hdf*i*dt));      
	   fprintf(fp,"%f\n",h[i]);
   }
   for(k=0;k<l1;k++)
   {
	   for(l=0;l<tl+hn;l++)
		   y[l]=0.0;
	   for(i=0;i<hn;i++)
		   for(j=0;j<tl;j++)
			   y[i+j]+=cdp[j][k]*h[i];
	   for(i=0;i<tl;i++)
		   Cdp[i][k]=y[i];		   
  }
  fclose(fp);
}