synthetic.c

su 的源代码库

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


#include "par.h"
#include "Reflect/reflpsvsh.h"
/******************************Self-Documentation*****************************/
/******************************************************************************
SYNTHETICS - Subroutine to compute a synthetic seismogram from the output
				of the reflectivity modeling code

*******************************************************************************
Function Prototypes:
void compute_synthetics (int verbose, int nt, int ntc, int nx, int nor, int nw,
	int nlayers, int lsource, int nf, int flt, int vsp, int win, int wtype, 
	float tlag, float red_vel, float decay, float tsec, float bx, float dx, 
	float w1, float w2, int *fil_phase, int nfilters, int *fil_type, 
	float *dbpo, float *f1, float *f2, int *lobs, float *cl, float *t, 
	complex ***response, float **reflectivity, FILE *outfp);

void red_vel_factor (float x, float red_vel, float tlag, complex wpie, 
	complex wsq, complex *rvfac);

void construct_tx_trace (int nw, int nt, int ntfft, int nfilters,
	int *filters_phase, float tsec, float unexp, float sphrd, complex *refw,
	int *fil_type, float *f1, float *f2, float *dbpo, float *reft);

void compute_Hanning_window (int iwin, int iw, int if1, int if2, int nf, 
	complex *win);

void apply_filters (int *min_phase, int nfilters, int nw, float tsec, float *f1,
	float *f2, float *dbpo, int *filter_type, complex *refw);
*******************************************************************************
compute_synthetics:
Input:
nt			number of time samples in the computed synthetics
tlag		time lag in seconds to be used in computing the synthetics
nf			number of frequencies to be used in computing the synthetics.
			Set to zero if all frequencies are to be used
flt			=1 to apply the earth flattening correction
vsp			=1 to compute synthetic vsp
win			=1 ray parameter windows ???
wtype		=1 for displacement
			=2 for velocity
			=3 for acceleration
w1			frequency to start lo end hanning taper. If w1=0, deafult=0.15*nf
			frequency to start hi end hanning taper. If w2=0, deafult=0.85*nf
red_vel		reducing velocity in km/s. If set to zero, use the highest
			model velocity
nfilters	number of filters to be applied to the synthetics
fil_phase	=1 for minimum phase filters
			=0 for zero phase filters
fil_type	array[nfilters] of filter flags: =1 for high cut, =2 for low cut,
			=3 for notch filter
dbpo		array[nfilters] of slopes in db/octave
f1			frequency 1 in hertz
f2 			frequency 2 in hertz (only in the case of notch filter)
response	array[nw]nx][nor] of computed wavefield amplitudes 

Output:
reflectivity	array[nx][nt] of computed reflectivities	
*******************************************************************************
compute_red_vel_factor
Input:
x			range (offset) of current trace 
red_vel		reducing velocity
tlag		time lag to be applied to the seismogram
wpie
wsq			complex factor

Output
rvfac		computed reducing velocity factor
*******************************************************************************
construct_tx_trace:
Input:
nw			number of frequencies
nt			number of time samples
ix			trace index
nfilt		number of filters to be applied
tsec		length of computed trace
unexp		parameter to undo decay factor	
sphrd		flattening earth correction factor
refw		array[nw+1] of complex trace

Output:
reft		array[nt] current trace in t-x domain
*******************************************************************************
compute_Hanning_window:
Input:
iwin		=1 for Hanning window
if1			lower f=window frequency
if2			higher window frequency
nw			number of frequencies in output data
iw			current frequency

Output:
win			computed Hanning window for given frequency
*******************************************************************************
apply_filters:
Input:
min_phase		=1 for minimum phase filters
				=0 for zero phase filters
nfilters		number of filters to apply
nw				number of frequencies
tsec
f1				array[nfilters] of low frequencies
f2				array[nfilters] of high frequencies
dbpo			array[nfilters] of filter slopes in db/oct
filtype			array[nfilters] of filter types:
				=1 low cut filter
				=2 high cut filter
				=3 notch filter
refw			array[nw+1] of complex samples to filter

Output:
refw			array[nw+1] of complex filtered samples 

Note:
refw contains only the positive frequencies, the negatives are simply their
complex conjugates and can be computed when necessary
*******************************************************************************/
/***********************************End Self_documentation********************/




/******************************************************************************

	Subroutine to generate the synthetic seismograms based on the
		information contained in the files wxp,wxr and wxz

******************************************************************************/
void compute_synthetics (int verbose, int nt, int ntc, int nx, int nor, int nw,
	int nlayers, int lsource, int nf, int flt, int vsp, int win, int wtype, 
	float tlag, float red_vel, float decay, float tsec, float bx, float dx, 
	float w1, float w2, int *fil_phase, int nfilters, int *fil_type, 
	float *dbpo, float *f1, float *f2, int *lobs, float *cl, float *t, 
	complex ***response, float **reflectivity, FILE *outfp)
/****************************************************************************** 
Input:
nt			number of time samples in the computed synthetics
tlag		time lag in seconds to be used in computing the synthetics
nf			number of frequencies to be used in computing the synthetics.
			Set to zero if all frequencies are to be used
flt			=1 to apply the earth flattening correction
vsp			=1 to compute synthetic vsp
win			=1 ray parameter windows ???
wtype		=1 for displacement
			=2 for velocity
			=3 for acceleration
w1			frequency to start lo end hanning taper. If w1=0, deafult=0.15*nf
			frequency to start hi end hanning taper. If w2=0, deafult=0.85*nf
red_vel		reducing velocity in km/s. If set to zero, use the highest
			model velocity
nfilters	number of filters to be applied to the synthetics
fil_phase	=1 for minimum phase filters
			=0 for zero phase filters
fil_type	array[nfilters] of filter flags: =1 for high cut, =2 for low cut,
			=3 for notch filter
dbpo			array[nfilters] of slopes in db/octave
f1			frequency 1 in hertz
f2 			frequency 2 in hertz (only in the case of notch filter)
response	array[nw]nx][nor] of computed wavefield amplitudes 

Output:
reflectivity	array[nx][nt] of computed reflectivities	
******************************************************************************/
{	
	int il,iw,ix,iz=0;	/* loop counters */
	int kz;			/* auxiliary index */
	int ntfft;		/* number of FFT samples */
	int nw1;		/* number of frequencies for FFT */
	float sum=0.0;		/* auxiliary variable */
	float *zr,zs=0.0;	/* depths of source and receivers */
	float sphrd=1.0;	/* factor for flattening earth correction */
	int if1=0,if2=0; 	/* freq indices to start and end Hanning window */
	float unexp;		/* parameter to undo decay factor */
	float x,z;		/* range and depth */
	float w;		/* frequency in rad/sec */
	complex wpie;		/* complex frequency */
	complex wsq;		/* complex frequency squared */
	complex rvfac;		/* reducing velocity factor */
	complex *refw;		/* scratch array for complex wavefield */
	complex window;		/* frequency domain Hanning window */

	/* if requested, check input parameters */
	if (verbose==1||verbose==3) {
		fprintf (stderr,"in compute reflectivity series, checking input ...\n");
		fprintf (stderr,"nt=%d nx=%d nor=%d nw=%d",nt,nx,nor,nw);
		fprintf (stderr," nlayers=%d lsource=%d nf=%d\n",nlayers,lsource,nf);
		fprintf (stderr,"flt=%d vsp=%d win=%d wtype=%d\n",flt,vsp,win,wtype);
		fprintf (stderr,"tlag=%g red_vel=%g decay=%g ",tlag,red_vel,decay);
		fprintf (stderr,"tsec=%g bx=%g dx=%g w1=%g w2=%g\n",tsec,bx,dx,w1,w2);
		fprintf (stderr,"nfilters=%d\n",nfilters);
	}
	if (verbose==2||verbose==3) {
		fprintf (outfp,"in compute reflectivity series, checking input ...\n");
		fprintf (outfp,"nt=%d nx=%d nor=%d nw=%d",nt,nx,nor,nw);
		fprintf (outfp," nlayers=%d lsource=%d nf=%d\n",nlayers,lsource,nf);
		fprintf (outfp,"flt=%d vsp=%d win=%d wtype=%d\n",flt,vsp,win,wtype);
		fprintf (outfp,"tlag=%g red_vel=%g decay=%g ",tlag,red_vel,decay);
		fprintf (outfp,"tsec=%g bx=%g dx=%g w1=%g w2=%g\n",tsec,bx,dx,w1,w2);
		fprintf (outfp,"nfilters=%d\n",nfilters);
	}

	/* compute number of fft samples and frequencies */
	ntfft=npfa(ntc);
	nw1=ntfft/2;
	if (nf>0 && nf<=nw) nw=nf;
	fprintf (stderr,"\n\n nt=%d ntfft=%d nw=%d nw1=%d\n",nt,ntfft,nw,nw1);

	/* allocate working space */
	zr=alloc1float(nor);
	refw=alloc1complex(ntfft);

	/* if not provided, set reducing velocity to highest pwave velocity */
	if (red_vel<0.0) {
		red_vel=1.51;	
		for (il=0; il<nlayers; il++) {
			if (cl[il]>red_vel) red_vel=cl[il];
		}
	}

	/* compute depths of source and receivers */
	iz=0;
	for (il=0; il<nlayers-1; il++) {
		kz=il+1;
		sum +=t[il];
		if (kz==lsource-1) zs=sum;	/* depth of the source */
		if ((iz<nor)&&(kz==lobs[iz]-1)) {
			zr[iz]=sum;		/* depth of the receivers */
			iz++;
		}
	}

	/* if requested, output processing information */
	if (verbose==1||verbose==3) {
		fprintf(stderr,"\nSOURCE AND RECEIVER INFORMATION\n");
		fprintf(stderr,"\nReducing velocity=%g source depth=%g\n",red_vel,zs);
		fprintf(stderr,"\nReceiver Number  Depth   \n");
		for (iz=0; iz<nor; iz++) {
			fprintf (stderr,"%8d%15.3f\n",iz,zr[iz]);
		}
	}
	if (verbose==2||verbose==3) {
		fprintf(outfp,"\nSOURCE AND RECEIVER INFORMATION\n");
		fprintf(outfp,"\nReducing velocity=%g source depth=%g\n",red_vel,zs);
		fprintf(outfp,"\nReceiver Number  Depth   \n");
		for (iz=0; iz<nor; iz++) {
			fprintf (outfp,"%8d%15.3f\n",iz,zr[iz]);
		}
	}
	
	/* compute windowed frequencies */
	if (win==1) {
		if (w1==0) if1=0.15*nf;
		else  if1=w1*tsec;
		if (w2==0) if2=0.85*nf;
		else if2=w2*tsec;
	}
	unexp=decay*tsec;					/* to undo exponential decay */

	/* main loop to compute a vsp or a seismogram */
	if (vsp==1) {
		for (ix=0; ix<nx; ix++) {

			/* compute range and flattening correction factor */
			x=bx+dx*ix;
			if (flt==1) sphrd=sqrt((x/RSO)/ABS(sin(x/RSO)));

			/* loop over the receivers */
			for (iz=0; iz<nor; iz++) {
				z=zr[iz];

				/* loop over frequencies */
				for (iw=0; iw<nw1; iw++) {

					/* initialize output array */
					if (iw<nw) {
						w=(iw+1)*PI*2/tsec;
						refw[iw+1]=response[iw][ix][iz];
					} else {
						w=0.0;
						refw[iw+1]=cmplx(0.0,0.0);
					}

					/* compute complex frequencies */
					wpie=cmplx(w,decay);
					wsq=cneg(cmul(cmplx(0.0,1.0),wpie));

					/* compute reducing velocity factor */
					red_vel_factor (x, red_vel, tlag, wpie,wsq, &rvfac);

					/* compute Hanning window */
					compute_Hanning_window (win, iw, if1, if2, nw, &window);

					/* apply Hanning window */
					refw[iw+1]=cmul(refw[iw+1],cmul(rvfac,window));
				}
			
				/* construct time image */
				construct_tx_trace (nw1, nt, ntc, ntfft, nfilters, fil_phase,
					tsec, unexp, sphrd, refw, fil_type, f1, f2, dbpo,
					reflectivity[ix]);
			}
		}

	} else {		/* compute "normal" seismogram */

		/* loop over recievers */	
		for (iz=0; iz<nor; iz++) {
			z=zr[iz];

			/* loop over ranges (seismogram traces) */
			for (ix=0; ix<nx; ix++) {
				x=bx+dx*ix;

				/* if earth flattening correction requested */
				if (flt==1) sphrd=sqrt((x/RSO)/ABS(sin(x/RSO)));
				
				/* loop over frequencies */
				for (iw=0; iw<nw1; iw++) {

					/* initialize output array */
					if (iw<nw) {