main.cpp

地震勘探纵波双曲时差速度分析程序

#include"Susgy.h"
#include<stdlib.h>
#include<malloc.h>
#include"math.h"
#include"stdio.h"
#include"fstream.h"


void main()
{	
	segy  tracehead;
	int   i,j,r,k,l;
	float t0,v,tx;
	float sc;
	int   vnum,t0num;
	int   tmin,tmax; 
	int   *x;
	float sum1,sum2,sum3,sum4;
	int m;
	
	
	////////////////可以改的参数///////////////////////////////
	
	int     ntrace=49,samnum=3001;
	float	cy=0.002,dt=0.02;
	float	vmin=1500,vmax=6000,dv=20;
	float     t0min=0,t0max=cy*samnum;
	
	//////////////////////////////////////////////////////////////////
	//  fp为要读的文件 fp1为得到的相似系数文件
	FILE  *fp,*fp1;   
	float **p,**p1; 
	vnum=1+int((vmax-vmin)/dv);//cout<<"vnum"<<vnum<<endl;
	t0num=1+int((t0max-t0min)/dt); //cout<<"tonum"<<t0num<<endl;
	//  p为道信息 p1为相似系数矩阵 x为炮检距
	p=(float**)malloc(ntrace*sizeof(float*));
	p1=(float**)malloc(vnum*sizeof(float*));
	x=(int*)malloc(ntrace*sizeof(int));
	for(i=0;i<ntrace;i++)
	{
		p[i]=(float*)malloc(samnum*sizeof(float));	
   	}
	for(i=0;i<vnum;i++)
	{
		p1[i]=(float*)malloc(t0num*sizeof(float));
	}
	fp=fopen("zz1100.segy","rb"); 
    //读入炮检距到数组x中
	fseek(fp,3600,0);
	for(i=0;i<ntrace;i++)  
	{
		
		fseek(fp,36,1);
		fread(&x[i],sizeof(int),1,fp);
		//	    printf("%d\n",x[i]);
		fseek(fp,200,1);
		fread(p[i],sizeof(float),samnum,fp);
	}// cout<<p[1][100]<<endl; 
	fclose(fp);
	//return;	
	//求相似系数
	
	for(k=0;k<t0num;k++)       //t0点数的变化
	{ //cout<<k<<endl;
		t0=k*dt;
		for(l=0;l<vnum;l++)
		{ 
			
			v=vmin+l*dv;
			tmin=int((t0-dt)/cy);       //cy为采样间隔，dt为时窗宽度的一半
			tmax=int((t0+dt)/cy);
			sum3=0.0;
			sum4=0.0;
			sc=0.0;
			for(i=tmin;i<tmax;i++)//时窗
			{   
				sum1=0.0;
				sum2=0.0;		          
				if(i>=0&&i<samnum)
				{					          
					
					for(j=0;j<ntrace;j++)
					{ 
						tx=sqrt(t0*t0+float(x[j]*x[j])/(v*v));
						r=int((tx-t0)/cy);
						m=i+r;
						if(m>=0&&m<samnum)
						{
							sum1=sum1+p[j][m];//printf("%f \n",sum1);
							sum2=sum2+pow(p[j][m],2);
						}
					}
				}
				sum3=sum3+sum1*sum1;
				sum4=sum4+sum2;
			}
			if(sum4!=0)
				sc=sum3/(ntrace*sum4);
			p1[l][k]=sc;
			//		cout<<"p1["<<l<<"]"<<"["<<k<<"]="<<p1[l][k]<<endl; 
		}
		cout<<t0num-k<<endl;
	} 
	//////////////////////////
	//为相似系数写道头
	fp1=fopen("zz1100.segy","wb"); 
	tracehead.dt=int(dt*1000000)+1; 
	tracehead.ntr=vnum; //横坐标点数 
	tracehead.d2=dv;
    tracehead.ns=t0num;//采样点数
	////////////////////////
	
	/////////////////////////////
	for(l=0;l<vnum;l++)
	{
		fwrite(&tracehead,sizeof(segy),1,fp1);
		for(k=0;k<t0num;k++)
		{ 
			fwrite(&p1[l][k],sizeof(float),1,fp1);
			//		 	printf("%f\n",**p1);
		} 
	}
	/////////////////////////////////////
	fclose(fp1);
	
}