taup.c

su 的源代码库

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


#include "taup.h"

/*********************** self documentation **********************/
/******************************************************************************
TAUP - Functions to perform forward and inverse taup transforms (radon or
	slant stacks) in t-x, F-K or F-X domains

fwd_FK_sstack		Performs forward taup tranaform in F-K domain, via
			8-point sinc interpolator. Only linear transform is
			possible. Very fast for taup transform of many traces.

fwd_tx_sstack		Performs forward taup transform in t-x domain.
			Only linear transform is implemented, but it is
			straight forward to implement other curves.

forward_p_transform	Performs forward tau-p transform in F-X domain. Uses
			Beylkin's approach. Not very fast but can compute
			linear parabolic or time-independent hyperbolic
			transform. Space coordinate does not need to be 
			uniformly sampled.

inv_FK_sstack		Performs inverse taup transform in F-K domain, via
			8-point sinc interpolation.

inv_tx_sstack		Performs inverse taup transform in t-x domain.

inverse_p_transform	Performs inverse taup transform in F-X domain (Beylkin's
			approach).

rho_filter		Computes rho filter in frequency domain for inverse
			t-x domain transform.

ga_xinterpolate        interpolate input data in space by placing a requested
			number of traces between each pair of input traces
*******************************************************************************
Function Prototypes:
void fwd_FK_sstack (float dt, int nt, int nx, float xmin, float dx, int np,
	float pmin, float dp, float fmin, float **traces, float **out_traces);

void fwd_tx_sstack (float dt, int nt, int nx, float xmin, float dx, int np,
	float pmin, float dp, float **traces, float **out_traces);

void forward_p_transform(int nx, int nt, float dt, float pmax, float pmin,
        float dp, float depthref, float f1, float f2, float freq1, float freq2,
        int lagc, int lent, int lenx, int xopt, int ninterp, int nwin,
        float prewhite, float interoff, float offref, int igopt, float dx,
        float fx, float pmula, float pmulb, float **in_traces,
        float **out_traces);

void inverse_p_transform(int nx, int nt, float dt, float pmax, float pmin,
	float dp, float depthref, float f1, float f2, float interoff,
	float offref, int igopt, float dx, float fx, float **in_traces,
	float **out_traces);

void inv_FK_sstack (float dt, int nt, int nx, float xmin, float dx, int np,
	 float pmin, float dp, float fmin, float **traces, float **out_traces); 

void inv_tx_sstack (float dt, int nt, int nx, int npoints, float xmin, 
	float dx, int np, float pmin, float dp, float **traces,
	float **out_traces); 

void rho_filter (int npoints, int nt, float dt, float *rho);

void ga_xinterpolate(float **in_traces, float **out_traces, int ninterp,
                int nt, int nx, float freq1, float freq2, int lagc,
                int lent, int lenx, int xopt, float dt, int iopt);

void runav(int n,int len,float *a,float *b);
*******************************************************************************
fwd_FK_sstack:
Input:
dt              time sampling interval
nt              number of time samples
nx              number of horizontal samples (traces)
np              number of slopes
pmin            minimum slope for tau-p transform
dp		slope sampling interval
fmin            minimum frequency of interest
traces          2-D array of input traces in t-x domain

Output:
traces          2-D array of output traces in tau-p domain
*******************************************************************************
Credits:
	Gabriel Alvarez (1994). Based on subroutine ..... in CWP program
	migbzo.c by Dave Hale.
*******************************************************************************
fwd_tx_sstack:
Input:
dt              time sampling interval
nt              number of time samples
nx              number of horizontal samples (traces)
np              number of slopes
pmin            minimum slope for tau-p transform
dp		slope sampling interval
traces          2-D array of input traces in t-x domain

Output:
out_traces      2-D array of output traces in tau-p domain
*******************************************************************************
Credits:
	Gabriel Alvarez (1994). 
*******************************************************************************
inv_FK_sstack:
Input:
dt              time sampling interval
nt              number of time samples
nx              number of horizontal samples
np              number of slopes
pmin            minimum slope for inverse tau-p transform
dp		slope sampling interval
fmin            minimum frequency of interest
traces          2-D array of input traces in tau-p domain

Output:
out_traces      2-D array of output traces in t-x domain
*******************************************************************************
Credits:
	Gabriel Alvarez (1994).  Based on subroutine ..... in CWP program
	migbzo.c by Dave Hale.
*******************************************************************************
inv_tx_sstack:
Input:
dt              time sampling interval
nt              number of time samples
nx              number of horizontal samples
np              number of slopes
pmin            minimum slope for inverse tau-p transform
dp		slope sampling interval
traces          2-D array of input traces in tau-p domain

Output:
out_traces      2-D array of output traces in t-x domain
*******************************************************************************
Credits:
	Gabriel Alvarez (1994).
*******************************************************************************
rho_filter:
Input:
npoints         number of point for the rho filter
nt              number of time samples
dt              time sampling interval

Output:
rho             1-D array of filter points
*******************************************************************************
Credits:
	Gabriel Alvarez (1994).
*******************************************************************************
forward_p_transform:
Input:
nx              number of input traces
nt              number of intput time samples
dt              sample rate in seconds
dx		offset sampling interval (distance between traces) (m).
fx		first offset (meters)
igopt           =1 parabolic transform g(x)=offset**2
                =2 Foster/Mosher pseudo hyperbolic transform
                        g(x)=sqrt(depth**2+offset**2)
                =3 linear tau-p g(x)=offset
                =4 abs linear taup g(x)=abs(offset)
offref          reference maximum offset to which maximum and minimum moveout
                times are associated
interoff        intercept offset to which tau-p times are associated
                (usually zero)
pmax            maximum moveout in ms on reference offset
pmin            minimum moveout in ms on reference offset
dp              moveout sampling interval (ms/m)
depthref        reference depth for Foster/Mosher hyperbolic transform
f1              high-end  frequency before taper off
f2              high-end frequency
prewhite        prewhitening factor in percent (usually between 0.01 and 0.1)
nwin            number of windows to use through the mute zone

Parameters with good suggested values:
freq1           low-end frequency for picking (usually 3 Hz)
freq1           high-end frequency for picking (usually 20 Hz)
lagc            length of AGC operator for picking (usually 400 ms)
lent            length of time smoother in samples for picking (usually 5)
lenx            length of space smoother in samples for picking (usually 1)
xopt            =1 use differences for spatial derivatives (works with
                irregular spacing)
                =0 use FFT derivative for spatial derivatives (more accurate
                but requires regular spacing and at least 16 input traces),
                will switch to differences automatically is this is not met
in_traces       2-D array of input t-x traces

Output:
out_traces      2-D array[np][nt] of output tau-p traces

Notes:
offsets are computed internally as offset[ix]=fx+ix*dx
************************************************************************
Credits:
	Adapted by Gabriel Alvarez (1995) from suradon.c written by John
        Anderson (1993)
************************************************************************
inverse_p_transform:
Input:
nx              number of output traces
nt              number of output time samples
dt              time sampling interval (seconds)
dx              spatial sampling interval (meters)
fx		first offset (meters)
igopt           =1 parabolic trransform g(x)=offset**2
                =2 Foster/Mosher pseudo hyperbolic transform
                        g(x)=sqrt(depth**2+offset**2)
                =3 linear tau-p g(x)=offset
                =4 abs linear taup g(x)=abs(offset)
		=5 new pseudo-hyperbolic ransform
			g(x)=1/ref_vel*sqrt((ref_time*ref_vel)**2+offset**2)
offref          reference maximum offset to which maximum and minimum moveout
                times are associated
interoff        intercept offset to which tau-p times are associated
                (usually zero)
pmax            maximum moveout in ms on reference offset
pmin            minimum moveout in ms on reference offset
dp              moveout sampling interval (ms/m)
depthref        reference depth for Foster/Mosher hyperbolic transform
f1              high-end  frequency before taper off (hz)
f2              high-end frequency (hz)
in_traces       2-D array[np][nt] of input taup traces

Output:
out_traces      2-D array[nx][nt] of output t-x traces
**************************************************************************
Credits:
	Adapted by Gabriel Alvarez (1995) from suradon.c written by John
        Anderson (1993)
**************************************************************************
ga_xinterpolate
Input:
int ninterp             number of traces to interpolate between each input trace
int nt                  number of time samples
int nx                  number of input traces
float freq1             low-end frequency for picking (good default: 3 Hz)
float freq2             high-end frequency for picking (good default: 20 Hz)
int lagc                length of AGC operator for picking(good default: 400 ms)
int lent                length of time smoother in samples for picker
                        (good default: 5 samples)
int lenx                length of space smoother in samples for picker
                        (good default: 1 sample)
int xopt                1 = use differences for spatial derivative
                            (works with irregular spacing)
                        0 = use FFT derivative for spatial derivatives
                            (more accurate but requires regular spacing and
                            at least 16 input tracs--will switch to differences
                            automatically if have less than 16 input traces)
float dt                sample rate in sec
int iopt		0 = interpolate: output 1+(nx-1)*(1+ninterp) traces
                           with ninterp traces between each pair of input traces
			1 = compute low-pass model: output nx traces
                            on original trace locations -- This is typically
                            used for Quality Control if the interpolator
                            is failing for any reason
			2 = compute dip picks in units of samples/trace: 
                            output nx traces on original trace locations
in_traces       	2-D array of input traces

Output:
out_traces      	2-D array of interpolated tau-p traces


Notes:
This routine outputs 'ninterp' interpolated traces between each pair of
input traces.  The values for lagc, freq1, and freq2 are only used for
event tracking. The output data will be full bandwidth with no agc.  The
suggested default parameters typically will do a satisfactory job of
interpolation for dips up to about 12 ms/trace.  Using a larger value for
freq2 causes the algorithm to do a better job on the shallow dips, but to
fail on the steep dips.  Only one dip is assumed at each time sample between
each pair of input traces.  The original input traces are passed through
this routine without modification.

The key assumption used here is that the low frequency data are unaliased
and can be used for event tracking.  Those dip picks are used to interpolate
the original full-bandwidth data, giving some measure of interpolation
at higher frequencies which otherwise would be aliased.  Using iopt equal
to 1 allows you to visually check whether the low-pass picking model is aliased.
If you can't visually pick dips correctly on the low-pass picking
model, this computer routine will fail.

The place this code is most likely to fail is on the first breaks.
************************************************************************
Credits:
	Adapted by Gabriel Alvarez (1995) from suradon.c written by John
        Anderson (1993)
************************************************************************
Notes:
Other subroutines, used internally, might be of interest:
gofx		computes offsets for linear, parabolic or hyperbolic transforms
freqweight	computes frequency dependent weigths
compute_r	computes top row of a Hermitian Toeplitz matrix
compute_rhs	computes hermitian matrix times data vector 
ctoep		complex hermitian Toeplitz solver
ctoephcg	Hestenes and Stiefel conjugate gradient algorithm especialized
		for solving Hermitian Toeplitz systems
rcdot		computes real part of a complex dot product where the first
		vector is the one complex conjugated
runav		computes a boxcar running average filter 


More documentation about these subroutines on their headings, below.
******************************************************************************/
/**************** end self doc ********************************/




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

        Compute global forward slant stack (tau-p transform) via
				FK  transform

******************************************************************************/
void fwd_FK_sstack (float dt, int nt, int nx, float xmin, float dx, int np,
	float pmin, float dp, float fmin, float **traces, float **out_traces)
/******************************************************************************
Input:
dt              time sampling interval
nt              number of time samples
nx              number of horizontal samples (traces)
xmin		minimum horizontal offset (ignored)
np              number of slopes
pmin            minimum slope for tau-p transform
dp		slope sampling interval
fmin            minimum frequency of interest (ignored)
traces          2-D array of input traces in t-x domain

Output:
traces          2-D array of output traces in tau-p domain
******************************************************************************/
{
        int it,itau,ix,iw,ik,ip;/* loop counters */
        float fx,fw;		/* first frequency */
        float fp;            	/* first slope */
        float dw,dk;            /* frequency and wavenumber sampling interval */
        int nw,nk;              /* number of frequencies and wavenumbers */
        int ntpad;              /* number of time samples to pad */
        int ntfft;              /* number of padded time samples */
	int ntau;		/*  */
	/* int lk; */		/* last wavenumber */
	int nxfft;		/* number of x samples to pad */
	float fftscl;		/* scale factor for FFT */
        float w,p,k;            /* frequency, slope and wavenumber */
        float phase;            /* phase angle */
        float xshift;           /* shift to account for non-centered x-axis */
        float c,s;              /* auxiliary variables */
        float xmax;             /* maximum horizontal value */
        float fk;          	/* first wavenumber */
        float pmax;             /* maximum slope */
        int lwrap;              /* samples required to make k periodic */
        float fka;              /* first wavenumber with k periodic */
	float temp;		/* auxiliary variable */
        int nka;                /* number of wavenumbers with k perioic */
        float *tr_fft;          /* padded trace for FFT */
        complex **ctr;          /* F-K transformed trace */
        complex **ctr_p;        /*  */
        complex *tr_k;          /* K-transformed trace  */
        complex *tr_ka;         /*  */
        complex *tr_x;          /* W-transformed input tracescaled  by dx */
        complex *hp;            /* slant stacked single trace */
        float *kp;		/* K-transformed slant stacked single trace */
	complex czero;		/* complex number zero */

        /* compute slope sampling interval */
        pmax = pmin+(np-1)*dp;
	czero = cmplx(0.0*fmin,0.0*xmin);

	/* determine lengths and scale factors for FFT's */
	fx=0.0;
	fp=pmin+0.5*(pmax-pmin-(np-1)*dp);
	pmax = (dp<0.0)?fp:fp+(np-1)*dp;
        xmax = (dx<0.0)?fx:fx+(nx-1)*dx;
	ntau = nt;
	ntpad = ABS(pmax*xmax)/dt;
	ntfft = npfar(MAX(nt+ntpad,ntau));
	fftscl = 1.0/ntfft;
	lwrap = 8;      /* samples required to make k periodic */
	nxfft = npfa(2*(nx+lwrap));