NAME
fradonf - generalized Discrete Radon Transform from (t,x)
to $( tau ,p)$, $ ( omega ,p)$ space.
SYNOPSIS
fradonf [ -Nfile_in ] [ -Ofile_out ] [ -mminmmin ] [
-mmaxmmax ] [ -npnp ] [ -sist ] [ -eiend ] [ -f1f1 ] [ -f2f2
] [ -f3f3 ] [ -f4f4 ] [ -tsemb1tsemb1 ] [ -tsemb2tsemb2 ] [
-tsemb3tsemb3 ] [ -tsemb4tsemb4 ] [ -xmaxxmax ] [ -zrefzref
] [ -ipwpw ] [ -alphaalpha ] [ -prewwhite ] [ -sigma1sigma1
] [ -sigma2sigma2 ] [ -sembwxsembwx ] [ -sembwtsembwt ] [
-nyquist ] [ -L ] [ -P ] [ -H ] [ -taup ] [ -time ] [ -omega
] [ -livenlive ] [ -nxtapernxtaper ] [ -nttapernttaper ] [
-Mamaxmem ] [ -V ] [ -? ]
DESCRIPTION
fradonf takes irregularly sampled seismic (t,x) gathers and
forward transforms them into the generalized $( tau ,p), (
omega , p)$ domain. The user models the seismic data with a
series of np linear, parabolic, hyperbolic curves which are
fit to the data in a least squares sense.
When coupled with other programs, fradonf can be used to
enhance coherent signal or attenuate coherent noise such as
linear noise trains and multiple reflections. The flow fra-
donf | polymute | fradonr produces results identical to the
modeled output of routine frmmult except that the user is
allowed to interact with the intermediate results in the $(
tau ,p) $ domain. Typically, input to the process is
unstacked common shot, common receiver, common midpoint or
common reflection point gathers. For multiple rejection, the
primary reflections should be flattened or overcorrected
using an appropriate NMO routine. Stacked or common offset
data may be processed by breaking the data into running win-
dows using the flow rwindow | fradonf | pred | polymute |
fradonr | rwindow -R.
Like conventional dip filtering, quadrant chopping, and
other $( omega ,k)$ techniques, this technique allows effec-
tive separation of coincident primaries and multiples in the
transform domain. With proper parameterization, one can
remove well separated multiples while minimizing any ampli-
tude effects on the primaries, including at the inside
traces.
By using semblance weighting, one can deal effectively with
spatially aliased data.
fradonf gets both its data and its parameters from command
line arguments. These arguments specify the input, output,
the temporal design window, the type of transformation to be
used, the minimum and maximum moveout to be modeled, the
number of parameters/curves to fit, tapers, and semblance
weighting if desired.
Command line arguments
-N file_in
Enter the input data set name or file immediately after
typing -N unless the input is from a pipe in which case
the -N entry must be omitted. This input file should
include the complete path name if the file resides in a
different directory. Example
-N/export/data2/china/cdp.gather tells the program to
look for file 'cdp.gather' in directory
'/export/data2/china'.
-O file_out
Enter the output data set name or file immediately
after typing -O. This output file is not required when
piping the output to another process. The output data
set also requires the full path name (see above).
-s ist
Window Start Time (msec). This window specifies the
region on the input data after which noise will be
suppressed. (default = first sample).
-e iend
Window End Time (msec). End of window for noise
suppression. (default = last sample)
-f1 f1
Frequency (Hz) at which we begin roll in of a Hamming
zero phase band pass filter (default = 5Hz)
-f2 f2
Frequency (Hz) at which we end roll in of a Hamming
zero phase band pass filter (default = f1).
-f3 f3
Frequency (Hz) at which we begin roll out of a Hamming
zero phase band pass filter (default = f4).
-f4 f4
Frequency (Hz) at which we end roll out of a Hamming
zero phase band pass filter . The cost increases with
the value of f4. (default = Nyquist)
-mmin mmin
Minimum value of differential moveout (msec) to be
modeled, measured at distance xmax. If the primary is
NMO-corrected to be flat, it is recommended to set
$mmin<=-200msec$ for a value of f1=5Hz , thereby
allowing the algorithm to model subtle AVO effects at
5Hz in the wavelet. No default.
-mmax mmax
Maximum value of differential moveout (msec) to be
modeled, measured at distance xmax. This value should
be slightly greater than the largest moveout the user
wishes to reject. No default.
-np np
Number of parameters or curves. About each record's
(actual or nonexisting) zero -offset trace a set of
linear, parabolic, hyperbolic curves are used to model
the data. The moveout range for these curves is speci-
fied in the -mmin, -mmax entries above. This parameter
determines how many curves are used to span that range.
The cost goes up linearly with the value of np.
(Default: np = $n sub trace$, where $n sub trace $ is
the number of traces in the gather, guaranteeing enough
locations to carry trace headers from forward back to
inverse transformation. see the -nyquist option below).
-nyquist
If present, np = 2*f4*.001*(mmax-mmin)+1 or two points
per wavelength at the farthest offset, as in routine
rmmult. (Default, use np=$ n sub trace $)
-xmax xmax
Maximum trace distance (ft or m) in the data. The
moveout entered after the -mmin, -mmax options are
referenced to this value of xmax. For linear moveout
(-L option), the minimum slowness $p sub min = m sub
min /x sub max = 1/v sub max $, while the maximum slow-
ness $p sub max = m sub max /x sub max = 1/v sub min $
. No default.
-prew white
Percent prewhitening used to stabilize the least-
squares inversion in the presence of noise. Default =
5%
-nttaper nttaper
Temporal taper length (msec). This taper is applied to
both the beginning and the end of the analysis window.
Smooth tapers are imperative in obtaining artifact free
Radon transforms! (default=100 msec)
-nxtaper nxtaper
Spatial taper length (traces). This taper is applied to
both the start and the end of the analysis window.
Smooth tapers are imperative in obtaining artifact free
Radon transforms! (default=5 traces)
-live minlive
Enter the minimum number of live traces accepted in a
gather. Gathers having less than this many live traces
will be zeroed out and flagged dead. (default = 3).
-M amaxmem
Enter the maximum memory allowed for the program in
Megawords. The program will try to store as many radon
transform matrices, $[R( omega ,x sub j ,p)]$, in
memory as will fit, thereby eliminating the need to
recompute them for the repeated source-receiver dis-
tances in subsequent gathers. (Default: amaxmem=16
Megawords on the Crays, amaxmem=4 Megawords on the
Sparcs and HPs).
-L Enter the command line argument '-L' to use linear
curves in the data model (No default. Other options are
-P, -H ).
-P Enter the command line argument '-P' to use parabolic
curves in the data model. (No default. Other options
are -L, -H ).
-H Enter the command line argument '-H' to use time
invariant hyperbolic curves in the data model (No
default. Other options are -L, -P ).
-zref zref
Reference depth for time invariant hyperbolic curves if
-H option selected. Hyperbolae defined as in Foster
and Mosher, Geophysics, 1992 . Default: zref=xmax).
-ipw pw
Weight the input (t,x) data using inverse power weight-
ing in the conventional $( tau ,p)$ transform before
normalization/orthogonalization. This weighting will
suppress high amplitude spurious events, including edge
effects. Maximum weight will be 1./pw times the mean
square energy of the current seismic gather to be pro-
cessed. (Default if ipw is entered, is 0.1, otherwise
the forward transform is performed without inverse
power weighting).
-alpha alpha
If present, form the forward $( tau ,p) $ transform
using an $alpha$ trim mean vs. conventional sum before
normalization/orthogonalization. This weighting will
suppress high amplitude spurious events, including edge
effects. $alpha$ = 1.0 will produce a conventional sum.
$alpha$ = .5 will sum half the data. (Default,
$alpha$=1.0, such that conventional summation is used).
-sigma1 sigma1
If present, reject data whose time average semblance
$sigma bar <= sigma sub 1 = sigma1$. See algorithm
description below. This weighting will accentuate
coherent/continuous events and suppress incoherent
events. (Default: NO semblance weighting. Suggestion:
sigma1 = 0.05)
-sigma2 sigma2
If present, pass data whose time average semblance
$sigma bar >= sigma sub 2 = sigma1$. See algorithm
description below. This weighting will accentuate
coherent/continuous events and suppress incoherent
events. (Default: NO semblance weighting.
Suggestion:sigma2 = 0.20)
-sembwx sembwx
Half width, $ n sub x $, in traces, of the running sem-
blance calculation window which is triggered by supply-
ing a non-zero value. See algorithm description below.
A reasonable number is sembwx = 5.
-sembwt sembwt
Half length, $w sub t =Kdt$ in msec, of the running
semblance calculation window. (Used only if -sigma1,-
sigma2 options invoked above. If so, default: sembwt =
24 msec)
-tsemb1 tsemb1
Time (msec) at which we begin roll in of semblance
weighting (Default = first sample).
-tsemb2 tsemb2
Time (msec) at which we end roll in of semblance
weighting (Default = tsemb2)
-tsemb3 tsemb3
Time (msec) at which we end roll off of semblance
weighting (Default = tsemb4)
-tsemb4 tsemb4
Frequency (Hz) at which we end roll out of a Hamming
zero phase band pass filter (Default = last sample).
-taup
Output the conventional $( tau ,p)$ transform ( vs. the
discrete Radon transform). This option allows one to
compare the subtle differences between discrete Radon
and coventional $( tau ,p) $ transforms. The conven-
tional $( tau ,p) $ transform is significantly cheaper
than the discrete Radon transformand sometimes accurate
results can be acheived for the -L option.
-time
If present, calculate the $( tau ,p) $ transform part
of the calculation in the time (vs. frequency) domain
by resampling the traces using a 6 point interpolator.
(Default if -sigma1, -sigma2, -ipw, or -alpha options
are used, otherwise the more efficient frequency domain
transform is used).
-omega
If present, output the results in the $( omega ,p)$ vs.
the $( tau ,p) $ domain (no matter how they were calcu-
lated). The results are packed with the first half of
the samples in each constant p trace containing the
amplitude, the second containing the phase. This option
is useful for examining aliasing effects and possible
match filtering applications.
-V Enter the command line argument '-V' to get additional
printout.
-? Enter the command line argument '-?' to get online
help. The program terminates after the help screen is
printed.
The Reverse 2D Radon Transform:
We define the reverse Radon transform (implemented in pro-
gram fradonr) for irregularly spaced data to be :
${u}=[R] {U}$,
where
$U sub m = U( omega ,p sub m )$,
$u sub j = u( omega , x sub j )$
$R sub {mj} = e sup { +i omega [p sub m theta (x sub j )
] }$,
$omega =$ the temporal frequency in radians/sec,
$p sub m = $ the mth apparent moveout in the x direction,
$x sub j = $ the position of the jth trace ,
$theta (x sub j ) = {x sub j } $ for a linear transform,
$theta (x sub j ) = {x sub j } sup 2 $ for a parabolic
transform,
$theta (x sub j ) = [{x sub j } sup 2 + {z sub ref } sup 2 ]
sup 1/2 $ for a hyperbolic trans form, and
$ z sub ref $ = an appropriate reference depth.
The Fast forward/inverse Radon transform using
The fast Unequally spaced FFT (USFFT) has been employed to
get a much faster forward and inverse Radon transform. The
forward and inverse discrete triginometric expansion is done
along for instance parabolic curves, with non equally spaced
points. The relative speed of the forward/inverse Radon
transform fradonf/fradonr as compared to the original
forward/inverse Radon transform radonf/radonr, increases
with the number of traces. For 800 traces, we get CPU
speedup in the order of 2.5-2.75/15-25 times for the
foward/inverse fast Radon transform respectively. The main
reason for the relative underperformance of the forward fast
Radon transform is the computationally demanding solution of
the Toeplitz matrix system of equations.
The Forward 2D Discrete Radon Transform:
We define the forward DRT such that when when it is followed
by the above reverse transform we can reconstruct the origi-
nal data, $u( omega , x sub j )$ in a least squares sense:
${U}=([R] sup * [R] + epsilon [I] ) sup {-1} [R] sup *
{u}$,
where
$R sub {mj} sup * = e sup { -i omega [p sub m theta (x sub
j ) ] }$,
$[I]$ = the identity matrix, and
$epsilon $ = a prewhitening factor.
The Conventional 2d Forward (tau,p) Transform:
The conventional $( tau , p )$ algorithm when followed by
above defined reverse transform will approximately recon-
struct the original data when using linear moveout (the -L
option below)
on equally sampled trace spacing. This algorithm should not
be used for nonlinear moveouts and sparse, irregularly sam-
pled data as encountered in 3D CMP gathers.
First define the unscaled $( tau ,p )$ transform in the time
(t,x) domain:
$U bar ( tau ,p) = sum from j=1 to J { u( tau -px sub j , x
sub j )} $
Next, transform from $( tau ,p)$ to $( omega ,p)$ and scale
with the so-called $ rho $ factor
:
$U ( omega , p) = rho ( omega ) U bar ( omega , p) $, where
$rho ( omega ) = sqrt {omega / 2 pi } $.
The Semblance weighted Forward (tau,p) Transform:
For gathers comprised of only a few ( < 20 ) traces, it is
useful to consider weighting/windo wing the output conven-
tional $( tau ,p) $ transform $U ( tau , p)$ by weights $w
( tau
,p) $ proportional to the time averaged semblance along the
same (linear) summation curves to obtain $ U hat ( tau ,p)
$:
$ U hat ( tau , p) = w( sigma bar ( tau ,p) U bar ( tau ,
p) $,
where
$w ( sigma bar ) = 0.$ if $ sigma bar <= sigma bar sub 1
$,
$w ( sigma bar ) = 1 over 2 [ 1 + cos({ sigma bar - sigma
bar sub 1 } over {sigma bar su b 2 - sigma bar sub 1} pi
)]$, if $ sigma bar sub 1 < sigma bar < sigma bar sub 2
$,
$w ( sigma bar ) = 1.$ if $ sigma bar >= sigma bar sub 2
$, and
$ sigma bar ( tau ,p)={ sum from k=-w/dt to k=+w/dt [ sum
from j=1 to j=J u ( tau -px sub j ,x sub j )] sup 2 } over
{J sum from k=-w/dt to k=+w/dt sum from j=1 to j=J [u ( tau
-px sub
j , x sub j )] sup 2 }$
The Running Window Semblance Weighted
Forward (tau,p) Transform:
For gathers comprised of many ( > 20 ) traces, it is useful
to consider weighting the input data u(t,x) by weights $w
( tau ,p, x sub j ) $ proportional to the time averaged semb
lance in a running window ranging from $-n sub x$ to $+n sub
x$ along the same (linear) summation curves used to calcu-
late $U ( tau ,p) $:
$ U tilde ( tau , p) = sum from j=1 to j=J w[ sigma bar (
tau-px sub j )] u( tau -px sub j )
$,
where
$w ( sigma bar )$ is defined as above, and
$ sigma bar ( tau - px sub j )={ sum from k=-w/dt to
k=+w/dt [ sum from {n=-n sub x } to {n= +n sub x } u ( tau
-px sub j+n ,x sub j+n )] sup 2 } over {J sum from {n=-n
sub x } to {n=+n
sub x } sum from j=1 to j=J [u ( tau -px sub j+n , x sub
j+n )] sup 2 }$
NOTE 1:
The flows fradonf | polymute | fradonr, fradonf | taupred |
fradonr and fradonf | polymute | taupred | fradonr can be
quite effective in eliminating coherent noise and/or multi-
ples. The flow fradonf | fradonr by itself can be used to
generate data and/or interpolate dead traces on a regular
surface grid.
NOTE 2:
Inverse velocity multiple rejection in common shot or
receiver domain can be difficult since the apices of the
parabolae/hyperbolae are not centered about x=0. Many more
parameters are needed to express these off centered
parabolae/hyperbolae than those than centered
parabolae/hyperbolae in the CMP domain. In general, filter-
ing of common shot, receiver of offset domain should be done
only with the -L or -F options. The -P and -H options show a
clean separation of signal and noise only in the common mid-
point and common reflection point (MBS) domains. (P. J.
Singer).
NOTE 3:
Considerable care must be taken in applying the radon
transform, yet preserving subtle AVO effects. Even if the
desired events are flat, the p=0 moveout coefficient is only
able to model a constant amplitude component of this event.
Other positive and/or negative p coefficients are needed to
model the observed AVO variation. Rejecting these essential
components of AVO in a multiple elimination process will
also wipe out the desired AVO effect! Thus, even if primary
events are flattened in CMP or CRP data, it is necessary to
model negative moveouts by setting $mmin<=-T$, where T
=1/f1. For example, if f1=5Hz, set $mmin<=-200msec$. (N. C.
Allegar).
NOTE 4:
The flow rwindow -ntr1 -npad5 | fradonf -semb | fradonr |
rwindow -R can be used to significantly increase the
coherency of post stack data.
NOTE 5:
The flow rwindow -ntr200 -npad21 | fradonf | taupred | fra-
donr | rwindow -R can be used in post stack suppression of
dipping multiples.
NOTE 6 STOLEN HEADERS:
In order to apply the discrete radon transform filtering to
3D post stack data, several trace headers need to be stolen.
Considerable care was taken not to steal commonly used
header words. Nevertheless, certain collisions may occur for
unforseen flows. Future plans call for using the dds system,
which will ameliorate the line header problem.
Line header words stolen by fradonf and fradonr for transform in the inline direction
(including normal shot gather, receiver gather and cmp gather transforms):
Header keyword variable
Crew00 't'
Crew01 'p'
Crew02 moveout ('L','P',or 'H')
MnUHTm ntr (number of input traces)
MxUHTm nfft (length of FFT)
TmMsFS f1
MutVel f4
NmSpMi df
MxRSEL zref
MnGrEl mmin
MxGrEl mmax
MxTrOf xmax
MxTrSt ist
Trace header words stolen by fradonf and fradonr for transform in the inline direction:
RedVel apparent velocity in ms/m (ms/ft).
DstUsg p*dt*1.e+07 (for use in routine taupred)
TVPT20 and TVPV20 p (as a floating point for use in routine mute3D)
TVPT19 x (inline distance as found in input word DstSgn)
MulSkw muteend (first non-zero sample of input data used to restore early mute)
StaCor dead trace flag 30000 changed to 30001
transform of dead trace flag 30001 changed to 30002
this allows pred and taupred to work
operation reversed on inverse transform
Line header words stolen by fradonf and fradonr for transform in the crossline direction
invoked by using the -Y option above:
Header keyword variable
Crew00 't'
Crew03 'q'
Crew04 moveout ('L','P',or 'H')
MnShDp ntr (number of input traces)
MxShDp nfft (length of FFT)
TmMsSl f1*1.e+06
TmSlIn f4*1.e+06
AERcPr df*1.e+06
IndAdj zref
MnSPEl mmin
MxSPEl mmax
MnTrOf xmax
MnTrSt ist
Trace header words stolen by fradonf and fradonr for transform in the crossline direction:
DstUsg q*dt*1.e+07 (for use in routine taupred)
TVPtsemb21 and TVPV21 q (as a floating point for use in routine mute3D)
TVPV19 y (crossline distance as found in input word DstSgn)
PREPIn muteend (first non-zero sample of input data used to restore early mute)
StaCor dead trace flag 30000 changed to 30001
transform of dead trace flag 30001 changed to 30002
this allows pred and taupred to work
operation reversed on inverse transform
Line header words stolen by rwindow in generating running window in the inline direction
(default mode in rwindow):
Header keyword variable
RATFld npad
OpGrFl nline (number of lines in the inline direction)
OrNREC nrec_inline (number of inline windowed records)
OrNTRC ntrpline (number of inline traces of original unwindowed line)
RATTrc ntrcumpads (ntr+2*npad)
Line header words stolen by rwindow in generating running window in the
Icrossline direction
(invoked by using the -Y option in rwindow):
Header keyword variable
FrstSP npad
DpN1SP nline (number of lines in the crossline direction)
NmDpIn nrec_inline (number of crossline windowed records)
NTrLnS ntrpline (number of crossline traces of original unwindowed line)
StWdFl ntrcumpads (ntr+2*npad)
NOTE 7:
Weighting/Muting the results of the transform according to
the semblance. To weight (mute) the transformed data by the
semblance along the entire moveout curves, simply supply the
-sigma1 and -sigma2 values. See Marfurt and Cottle (1994)
for examples.
NOTE 8:
Weighting the transform with running semblance To apply a
running semblance weight to the data as it is being
transformed, supply the -sembwx and -sembwt options as well
as the -sigma1 and -sigma2 values above. To avoid articacts
that may arise, the user should compare the results with
those obtained by the conventional $( tau ,p) $ transform
obtaining by invoking the -taup option. See Marfurt and Vas-
siliou (1994) for examples.
REFERENCES
Marfurt, K. J., Schneider, R. V. and Mueller, M. C., 1994,
Challenges in seismic processing of irregularly sampled,
finte aperture aliased data using discrete Radon $( tau ,p)
$ transforms: APR Geos. Res. Bul. F94-G-14.
Marfurt, K. J. and Cottle, D. A. 1994, A comparison of
coherency weighted $( tau ,p) $ filtering: Application to
poststack and common offset data: APR Geos. Res. Bul. F94-
G-16.
Marfurt, K. J. and Vassiliou, A. A., 1994, in press.
Yilmaz, O. and Tanner, T. 1994, Discrete plane wave decompo-
sition by least-squares method, Geophysics: 59, 973-982.
Stoffa, P. L., Buhl, P. Diebold, J. and Wenzel, F., 1981,
Direct mapping of seismic data to the domain of intercept
time and ray parameter - A plane-wave decomposition:
Geophysics, 46, 255-267.
BUGS
No bugs known at present.
SEE ALSO
fradonr,polymute,taupred,rmmult,rwindow,radon3d,swfilt3d,taupf,taupr,slstkf,slstkr
AUTHORS
Kurt J. Marfurt. Original parabolic moveout version built
upon earlier work in routine rmmult (1992). Alpha trim mean
option added 6/93 by K. D. Crawford. 2D by 2D post stack
processing features added 12/93. Semblance weighting
features added 3/94. Least square DFT option added 5/94.
Fast Unequal Spaced FFT (USFFT) implemented in fradonf in
10/97 by Anthony A. Vassiliou.
COPYRIGHT
copyright 2001, Amoco Production Company
All Rights Reserved
an affiliate of BP America Inc.
Man(1) output converted with
man2html