avo正演模型实验_直接形成cdp.cpp

AVO正演模型实验

//---------------------------------------------------------------//
//                 AVO 技术正演模型实验                          //
//                 3 层模型                                      //
//                 AVO 技术正演模型实验                          //
//                 40炮24道                                       //
//                 地层倾角为0度                                 //
//---------------------------------------------------------------//
//---------------------------------------------------------------//
#include <stdio.h>
#include<math.h>
#define  PI   3.1415926
#define nc 2             //地层层数
#define np 4             //地层参数个数
#define  dt   0.002      //采样间隔 
#define  zbn   51        //子波长度
#define  zf   30         //子波主频
#define  zA   1.0        //子波最大振幅



//——————————观测系统参数————————————//
#define  Nd   24         //迭加次数
#define  dx   5         //道间距
#define  ds   5         //炮间距//
#define  x_0  10        //最小偏移距
#define  L    48         //接收道数
#define  Ns   40         //炮数
#define  tl   300        //采样点数
#define  ncdp  100      //cdp道数


void main()
{
	int   i,j,k;
	float par[nc][np];
	float x[Ns][L],zb[zbn];
	FILE *fp1;
	FILE *fp2;
	FILE *fp3;
	FILE *fp4;
	FILE *fp5;
    FILE *fp6;
    //====  子程序说明  ====//
	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");
	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%2==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);
}
//=======================  子波子程序  =====================//
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,j,k,k2,kn,r1,r2,l2,flag;
	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];//含油气层参数
	float sa2_j,sa_j,ti0,t00,vi;
	float cdp_1[tl]={0.0},cdp_2[tl]={0.0},cdp_3[tl]={0.0},tp1[tl]={0.0};//角度叠加约束
	float cdp[tl]={0.0},tp[tl];
	printf("*********\n");
	FILE *fp1;
	FILE *fp2;
	FILE *fp3;
	FILE *fp4;
	FILE *fp5;
	FILE *fp6;
    //===== 计算基本参数 =====//
	fp1=fopen("反射系数.txt","w+");
	fp2=fopen("地震记录.dat","wb+");
	fp3=fopen("地震记录(文本).dat","w+");
	fp4=fopen("含油气层参数.txt","r");
	fp5=fopen("角度叠加(3组角度范围＋全角度叠加共4个300CDP).dat","wb+");
	fp6=fopen("角度叠加叠加区域.txt","w+");
	fprintf(fp6,"道号\t采样点\t值\n");
	//===== 读取含油气参数 =====//
	fscanf(fp4,"%f\t%f\t%f\t%f\n",&par_0[0],&par_0[1],&par_0[2],&par_0[3]);
   //——————————考虑AVO的地震记录形成————————————//
	
	for(i=0;i<ncdp;i++)
	{
		x1=x_0+dx*i;
		if(x1/dx>=100&&x1/dx<200)//含油气参数替换
			for(k2=0;k2<np;k2++)
				par[1][k2]=par_0[k2];
	    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];
	       	t_0[i1]+=2*(t[i1]);
	       	sum+=t[i1]*par[i1][0]*par[i1][0];
		}
		vr=sqrt((sum*2)/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;
			cdp[i1]=0.0;
			cdp_1[i1]=0.0;
			cdp_2[i1]=0.0;
			cdp_3[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)));
			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];		
		}
		fprintf(fp1,"%f\t%f\n",asin(sqrt(sa2[0]))*180/PI,R1[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];
			fprintf(fp3,"%f\t",yy[i1]+0.005*i1);//形成理论CDP道集
		}
		fwrite(yy,4,tl,fp2);
		fprintf(fp3,"\n");
		flag=1;
        while(flag<5)
		{
			//———  角度叠加～～确定角度范围↓...———//
		    for(r2=0;r2<tl;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<12)
		    		cdp_1[r2]+=yy[r2];
		    	else
		    		cdp_1[r2]+=0.0;
	    		if(sa_j>=6&&sa_j<19)
				{
					cdp_2[r2]+=yy[r2];
					fprintf(fp6,"%d\t%d\t%d\n",r2,i,1);
				}
	    		else
				{
					cdp_2[r2]+=0.0;
					fprintf(fp6,"%d\t%d\t%d\n",r2,i,0);
				}
	    		if(sa_j>=22&&sa_j<35)
	    			cdp_3[r2]+=yy[r2];
		    	else
		    		cdp_3[r2]+=0.0;
	    		//———  角度叠加～～确定角度范围↑...———//
	    		cdp[r2]+=yy[r2];						
			}
			if(flag==1)  fwrite(cdp_1,4,tl,fp5);
			if(flag==2)  fwrite(cdp_2,4,tl,fp5);
			if(flag==3)  fwrite(cdp_3,4,tl,fp5);
		    if(flag==4)  fwrite(cdp,4,tl,fp5);
			fwrite(tp1,4,tl,fp5);
			//printf("%d\n",flag);
			flag++;
		}

	}
	printf("X总测线长度：%f\t总道数：%f\n",x1,x1/dx);
	fclose(fp1);
	fclose(fp2);
	fclose(fp3);
	fclose(fp4);
	fclose(fp5);
	fclose(fp6);
}