velpertan.c

su 的源代码库

/* Copyright (c) Colorado School of Mines, 2006.*/
/* All rights reserved.                       */

#include "par.h"
#define EPS_smooth FLT_MIN

/*********************** self documentation **********************/
char *sdoc[]={
"									",
" VELPERTAN - Velocity PERTerbation analysis in ANisotropic media to    ",
"             determine the model update required to flatten image gathers",  
"									",
" velpertan boundary= par=cig.par refl1= refl2= npicks1= 		",
"	npicks2= cdp1= cdp2= vfile= efile= dfile= nx= dx= fx= 		",
"	ncdp= dcdp= fcdp= off0= noff= doff= >outfile [parameters]	",
"									",
" Required Parameters:							",
" refl1=	file with picks on the 1st reflector	  		",
" refl2=	file with picks on the 2nd reflector  			",
" vfile=	file defining VP0 at all grid points from prev. iter.	", 
" efile=	file defining eps at all grid points from prev. iter.	", 
" dfile=	file defining del at all grid points from prev. iter.	", 
" boundary=	file defining the boundary above which 			",
"        	parameters are known; update is done below this		", 
"		boundary						", 
" npicks1=	number of picks on the 1st reflector			",
" npicks2=	number of picks on the 2nd reflector			",
" ncdp=		number of cdp's		 				",
" dcdp=		cdp spacing			 			",
" fcdp=		first cdp			 			",
" off0=		first offset in common image gathers 			",
" noff=		number of offsets in common image gathers  		",
" doff=		offset increment in common image gathers  		",
" cip1=x1,r1,r2,..., cip=xn,r1n,r2n         description of input CIGS	",
" cip2=x2,r1,r2,..., cip=xn,r1n,r2n         description of input CIGS	",
"	x	x-value of a common image point				",
"	r1	hyperbolic component of the residual moveout		",
"	r2	non-hyperbolic component of residual moveout		",
" Optional Parameters:							",
" method=akima         for linear interpolation of the interface       ",
"                       =mono for monotonic cubic interpolation of interface",
"                       =akima for Akima's cubic interpolation of interface ",
"                       =spline for cubic spline interpolation of interface ",
" VP0=2000	Starting value for vertical velocity at a point in the   ",
"							target layer	",
" x00=0.0	x-coordinate at which VP0 is defined			",
" z00=0.0	z-coordinate at which VP0 is defined			",
" eps=0.0	Starting value for Thomsen's parameter epsilon		",
" del=0.0	Starting value for Thomsen's parameter delta		",
" kz=0.0	Starting value for the vertical gradient in VP0		",
" kx=0.0	Starting value for the lateral gradient in VP0		",
" nx=100	number of nodes in the horizontal direction for the     ",
"							velocity grid 	",
" nz=100	number of nodes in the vertical direction for the	",
"							velocity grid	",
" dx=10	horizontal grid increment				",
" dz=10	vertical grid increment					",
" fx=0		first horizontal grid point				",
" fz=0		first vertical grid point				",
" dt=0.008	traveltime increment					",
" nt=500	no. of points on the ray				",
" amax=360	max. angle of emergence					",
" amin=0	min. angle of emergence					",
"									",
" Smoothing parameters:							",
" r1=0                  smoothing parameter in the 1 direction          ",
" r2=0                  smoothing parameter in the 2 direction          ",
" win=0,n1,0,n2         array for window range                          ",
" rw=0                  smoothing parameter for window function         ",
" nbound=2	number of points picked on the boundary			",
" tol=0.1	tolerance in computing the offset (m)			",
" Notes:								",
" This program is used as part of the velocity analysis technique developed",
" by Debashish Sarkar, CWP:2003.					",
"									",
" Notes:								",
" The output par file contains the coefficients describing the residual ",
" moveout. This program is used in conjunction with surelanan.		",
"									",
NULL};
/*
 * Author: CSM: Debashish Sarkar, December 2003 
 * based on program: velpert.c written by Zhenuye Liu.
 */

/**************** end self doc ***********************************/

/* functions defined and used internally */
/* one step along ray */
typedef struct RayStepStruct {
        float t;                /* time */
        float x,z;              /* x,z coordinates */
        float q1,p1,q2,p2;      /* Cerveny's dynamic ray tracing solution */
        int kmah;               /* KMAH index */
        float c,s;              /* cos(angle) and sin(angle) */
        float v,dvdx,dvdz;      /* velocity and its derivatives */
} RayStep;

/* one ray */
typedef struct RayStruct {
        int nrs;                /* number of ray steps */
        RayStep *rs;            /* array[nrs] of ray steps */
        int nc;                 /* number of circles */
        int ic;                 /* index of circle containing nearest step */
        void *c;                /* array[nc] of circles */
} Ray;

void interpolate(int ncdp, float *xcdp, float *xint, float *zint, float *zcdp, float *zdcdp, int np, char *method);
float velpertan_time(float x00,float z00,float *xxbound,float *zzbound,float x,float z,
	float zd,float off,float VP0,float eps,float del,float kz,float kx,
	int nx,int nz, float fx,float dx,float fz,float dz,int nt,float dt,
	float amin, float amax,float tol,float *p,float *vp,float *e,float *d,
	float r1,float r2,int *win,float rw);

/* functions (velocity interpolation)*/
void* vel2Alloc (int nx, float dx, float fx,
        int nz, float dz, float fz, float **v);
void vel2Free (void *vel2);
void vel2Interp (void *vel2, float x, float z,
        float *v, float *vx, float *vz, float *vxx, float *vxz, float *vzz);

Ray *makeRay (float x0, float z0, float a0, int nt, float dt,
	float **a1111xz, float **a3333xz, float **a1133xz, float **a1313xz, 
	float **a1113xz, float **a3313xz, 
	int nx, float dx, float fx, int nz, float dz, float fz, float amax, 
	float amin);
float zbrentou(float da,float da_2,float tol,float off,float x, float z,
	float a_normal,int nt,float dt,float **a1111,
	float **a3333,float **a1133, float **a1313,float **a1113,
	float **a3313,int nx,float dx,float fx,int nz,float dz,
	float fz,float amax,float amin);
void freeRay (Ray *ray);
void conjugate_gradient(float **A, float *x, float *b, int nrx, float tol);
void smooth2 (int n1, int n2, float r1, float r2, int *win, float rw, float **v);

int main (int argc, char **argv)
{
	int npicks1,npicks2,ipicks,ncip1,icdp,win[4];
	int ncdp,noff,ioff,nx,nz,nt,ipar1,ipar2,nbound;
	float dcdp,fcdp,doff,off0,off,VP0,kz,kx,eps,del,dt,p;
	float fx,fz,dx,dz,amin,amax,tol,x00,z00;
	float r1,r2,rw;
	float *xcdp, *zcdp1, *zcdp2, *zdcdp1, *zdcdp2;
	float *r11, *r12, *r21, *r22, *x, *xbound, *zbound;
	float *xxbound, *zzbound, *zdbound;
	float *xcip, *r11cip, *r12cip, *r21cip, *r22cip, temp[3];
	float ***t, ***tkz, ***tkx, ***td, *b1, *b2, *b, **A1, **A2, **A;
	float ***dt1, ***dt2, **zz1, **zz2, ***pp, ***tdV, ***tde;
	float **dt1avg,**dt2avg,*zz1avg,*zz2avg,*vp,*e,*d;
	float sum1,sum2,sum3,sum4,sum5,sum6,sum7,sum8,sum9,sum10,sum11,sum12;
	char *refl1,*refl2,*vfile,*efile,*dfile;
	float *zint1,*xint1,*zint2,*xint2,*norm;
	char *method="akima";
	char *boundary="boundary";
	FILE *rf1,*rf2,*bound,*vf,*ef,*df;

	/* hook up getpar to handle the parameters */
	initargs(argc,argv);
	requestdoc(0);

	/* get required parameters */
	if (!getparint("ncdp",&ncdp)) err("must specify ncdp!\n");
	if (!getparfloat("dcdp",&dcdp)) err("must specify dcdp!\n");
	if (!getparfloat("fcdp",&fcdp)) err("must specify fcdp!\n");
	if (!getparint("npicks1",&npicks1)) err("must specify npicks1!\n");
	if (!getparint("npicks2",&npicks2)) err("must specify npicks2!\n");
	if (!getparfloat("off0",&off0)) err("must specify off0!\n");
	if (!getparfloat("doff",&doff)) err("must specify doff!\n");
	if (!getparint("noff",&noff)) err("must specify noff!\n");
	if (!getparstring("refl1",&refl1)) err("must specify reflector 1!\n");
	if (!getparstring("refl2",&refl2)) err("must specify reflector 2!\n");
	if (!getparstring("vfile",&vfile)) err("must specify velocity file!\n");
	if (!getparstring("efile",&efile)) err("must specify epsilon file!\n");
	if (!getparstring("dfile",&dfile)) err("must specify delta file!\n");
	if (!getparstring("boundary",&boundary)) err("must specify boundary!\n");

	/* get optional parameters */
	getparstring("method",&method);
	if(!getparfloat("VP0",&VP0)) VP0=2000;
	if(!getparfloat("kz",&kz)) kz=0.0;
	if(!getparfloat("kx",&kx)) kx=0.0;
	if(!getparfloat("eps",&eps)) eps=0.0;
	if(!getparfloat("del",&del)) del=0.0;
	if(!getparint("nx",&nx)) nx=100;
	if(!getparint("nz",&nz)) nz=100;
	if(!getparfloat("dx",&dx)) dx=10;
	if(!getparfloat("dz",&dz)) dz=10;
	if(!getparfloat("fx",&fx)) fx=0;
	if(!getparfloat("fz",&fz)) fz=0;
	if(!getparfloat("dt",&dt)) dt=0.008;
	if(!getparint("nt",&nt)) nt=500;
	if(!getparfloat("amax",&amax)) amax=180.0;
	if(!getparfloat("amin",&amin)) amin=0.0;
	if(!getparfloat("tol",&tol)) tol=0.1;
	if(!getparint("nbound",&nbound)) nbound=2;
	if(!getparfloat("x00",&x00)) x00=0;
	if(!getparfloat("z00",&z00)) z00=0;
	if (!getparint("win",win)) {
                win[0] = 0;
                win[1] = nz;
                win[2] = 0;
                win[3] = nx;
                }
	if (!getparfloat("rw",&rw)) rw = 0;
	if (!getparfloat("r1",&r1)) r1 = 0;
	if (!getparfloat("r2",&r2)) r2 = 0;
	
	ncip1 = countparname("cip1");
	if (ncip1<1) err("Number of CIGS must be greater 0!\n");

	rf1=fopen(refl1,"r");
	rf2=fopen(refl2,"r");
	bound=fopen(boundary,"r");
	vf=fopen(vfile,"r");
	ef=fopen(efile,"r");
	df=fopen(dfile,"r");

	/* allocate space */
	xcip = alloc1float(ncip1);
	r11cip = alloc1float(ncip1);
	r21cip = alloc1float(ncip1);
	r22cip = alloc1float(ncip1);
	r12cip = alloc1float(ncip1);
	xcdp = alloc1float(ncdp);
	zcdp1 = alloc1float(ncdp);
	zcdp2 = alloc1float(ncdp);
	zdcdp1 = alloc1float(ncdp);
	zdcdp2 = alloc1float(ncdp);
	r11 = alloc1float(ncdp);
	r12 = alloc1float(ncdp);
	r21 = alloc1float(ncdp);
	r22 = alloc1float(ncdp);
	zint1 = alloc1float(npicks1);
	xint1 = alloc1float(npicks1);
	zint2 = alloc1float(npicks2);
	xint2 = alloc1float(npicks2);
	t     = alloc3float(noff,ncdp,2);
	tkx    = alloc3float(noff,ncdp,2);
	tkz    = alloc3float(noff,ncdp,2);
	td    = alloc3float(noff,ncdp,2);
	tdV    = alloc3float(noff,ncdp,2);
	tde    = alloc3float(noff,ncdp,2);
	b1 = alloc1float(5); 
	b2 = alloc1float(5); 
	b = alloc1float(5); 
	A1 = alloc2float(5,5); 
	A2 = alloc2float(5,5); 
	A = alloc2float(5,5); 
	dt1 = alloc3float(noff,ncdp,5);
	dt2 = alloc3float(noff,ncdp,5);
	zz1 = alloc2float(noff,ncdp);
	zz2 = alloc2float(noff,ncdp);
	dt1avg = alloc2float(ncdp,5);
	dt2avg = alloc2float(ncdp,5);
	zz1avg = alloc1float(ncdp);
	zz2avg = alloc1float(ncdp);
	xbound = alloc1float(nbound);
	zbound = alloc1float(nbound);
	xxbound = alloc1float(nx);
	zzbound = alloc1float(nx);
	zdbound = alloc1float(nx);
	x = alloc1float(5);
	pp = alloc3float(noff,ncdp,2);
	vp = alloc1float(nx*nz);
	e = alloc1float(nx*nz);
	d = alloc1float(nx*nz);
	norm = alloc1float(5);
	
	memset((void *) pp[0][0], 0, FSIZE*noff*ncdp*2);
        memset((void *) dt1[0][0], 0, FSIZE*noff*ncdp*5);
        memset((void *) dt2[0][0], 0, FSIZE*noff*ncdp*5);
        memset((void *) t[0][0], 0, FSIZE*noff*ncdp*2);
        memset((void *) tkx[0][0],0, FSIZE*noff*ncdp*2);
        memset((void *) tkz[0][0], 0, FSIZE*noff*ncdp*2);
        memset((void *) td[0][0], 0, FSIZE*noff*ncdp*2);
        memset((void *) tdV[0][0], 0, FSIZE*noff*ncdp*2);
        memset((void *) tde[0][0], 0, FSIZE*noff*ncdp*2);
        memset((void *) dt1avg[0], 0, FSIZE*5*ncdp);
        memset((void *) dt2avg[0], 0, FSIZE*5*ncdp);
        memset((void *) zz1[0], 0, FSIZE*noff*ncdp);
        memset((void *) A1[0], 0, FSIZE*5*5);
        memset((void *) A2[0], 0, FSIZE*5*5);
        memset((void *) A[0], 0, FSIZE*5*5);
        memset((void *) b1, 0, FSIZE*5);
        memset((void *) b2, 0, FSIZE*5);
        memset((void *) b, 0, FSIZE*5);
        memset((void *) xbound, 0, FSIZE*nbound);
        memset((void *) zbound, 0, FSIZE*nbound);
        memset((void *) norm, 0, FSIZE*5);

	/* input velocity, epsilon, and delta fields */
	fread(vp,sizeof(float),nz*nx,vf);
	fread(e,sizeof(float),nz*nx,ef);
	fread(d,sizeof(float),nz*nx,df);

	for(icdp=0; icdp<ncip1; ++icdp){
		getnparfloat(icdp+1,"cip1",temp);
		xcip[icdp] = temp[0];
		r12cip[icdp] = temp[1];
		r11cip[icdp] = temp[2];
	}

	for(icdp=0; icdp<ncip1; ++icdp){
		getnparfloat(icdp+1,"cip2",temp);
		r22cip[icdp] = temp[1];
		r21cip[icdp] = temp[2];
	}
	/* Input picked interface 1 */
           for(ipicks=0;ipicks<npicks1;ipicks++)
                fscanf(rf1,"%f %f\n", &zint1[ipicks], &xint1[ipicks]);
	/* Input picked interface 2 */
           for(ipicks=0;ipicks<npicks2;ipicks++)
                fscanf(rf2,"%f %f\n", &zint2[ipicks], &xint2[ipicks]);

	/* input the boundary to demarcate the region to update */
	   for(ipicks=0;ipicks<nbound;ipicks++)
		 fscanf(bound,"%f %f\n",&zbound[ipicks],&xbound[ipicks]);
	/* interpolate the boundary to all grid positions: spline, akima, monotonic*/
	for(ipicks=0;ipicks<nx;ipicks++) xxbound[ipicks]=fx+ipicks*dx;
	interpolate(nx,xxbound,xbound,zbound,zzbound,zdbound,nbound,method);

	/* compute uniformly sampled cdp */
                 for(icdp=0; icdp<ncdp; ++icdp)
                        xcdp[icdp] = fcdp+icdp*dcdp;

	/* reflector interpolation: spline, akima, monotonic */
	interpolate(ncdp,xcdp,xint1,zint1,zcdp1,zdcdp1,npicks1,method);
	interpolate(ncdp,xcdp,xint2,zint2,zcdp2,zdcdp2,npicks2,method);

	/* linear interpolation for the residual moveout components */
	intlin(ncip1,xcip,r11cip,r11cip[0],r11cip[ncip1-1],ncdp,xcdp,r11);
	intlin(ncip1,xcip,r12cip,r12cip[0],r12cip[ncip1-1],ncdp,xcdp,r12);
	intlin(ncip1,xcip,r21cip,r21cip[0],r21cip[ncip1-1],ncdp,xcdp,r21);
	intlin(ncip1,xcip,r22cip,r22cip[0],r22cip[ncip1-1],ncdp,xcdp,r22);

	/* find true and perturbed times for each cdp and each offset */
	for(icdp=0;icdp<ncdp;icdp++){
		xcdp[icdp]=fcdp+dcdp*icdp;
		if(xcdp[icdp]>(fx+(nx-1)*dx)){
			warn("cdp exceeds grid dimensions\n"); break;
			}
			for(ioff=0;ioff<noff;ioff++){
			off = off0+doff*ioff; p=0;
		t[0][icdp][ioff]=velpertan_time(x00,z00,xxbound,zzbound,xcdp[icdp],
		zcdp1[icdp],zdcdp1[icdp],off,VP0,eps,del,kz,kx,nx,nz,fx,
		dx,fz,dz,nt,dt,amin,amax,tol,&p,vp,e,d,r1,r2,win,rw);
		pp[0][icdp][ioff]=p; p=0;
		t[1][icdp][ioff]=velpertan_time(x00,z00,xxbound,zzbound,xcdp[icdp],
		zcdp2[icdp],zdcdp2[icdp],off,VP0,eps,del,kz,kx,nx,nz,fx,
		dx,fz,dz,nt,dt,amin,amax,tol,&p,vp,e,d,r1,r2,win,rw); 
		pp[1][icdp][ioff]=p; p=0;

	       	td[0][icdp][ioff]=velpertan_time(x00,z00,xxbound,zzbound,
		xcdp[icdp],zcdp1[icdp],zdcdp1[icdp],off,VP0,eps,del+0.05,