NAME
radonf - generalized Discrete Radon Transform from (t,x) to
$( tau ,p)$, $ ( omega ,p)$ or $ ( omega ,k) $ space.
SYNOPSIS
radonf [ -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 ] [ -K ] [ -taup ] [ -time ] [
-omega ] [ -livenlive ] [ -nxtapernxtaper ] [ -ntta-
pernttaper ] [ -Mamaxmem ] [ -V ] [ -? ]
DESCRIPTION
radonf takes irregularly sampled seismic (t,x) gathers and
forward transforms them into the generalized $( tau ,p), (
omega , p)$ or $( omega , k)$ domain. The user models the
seismic data with a series of np linear, parabolic, hyper-
bolic curves (or for the least squares discrete Fourier
transform, with sines and cosines) which are fit to the
data in a least squares sense.
When coupled with other programs, radonf can be used to
enhance coherent signal or attenuate coherent noise such as
linear noise trains and multiple reflections. The flow
radonf|polymute|radonr produces results identical to the
modeled output of routine rmmult 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|radonf|pred|polymute|radonr|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.
radonf 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 allow-
ing 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 (or sines and
cosines for the -K option) are used to model the data.
The moveout range for these curves is specified 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 and -K).
-P Enter the command line argument '-P' to use parabolic
curves in the data model. (No default. Other options
are -L, -H and -K).
-H Enter the command line argument '-H' to use time
invariant hyperbolic curves in the data model (No
default. Other options are -L, -P and -K).
-K Enter the command line argument '-K' to perform a least
square $ ( omega ,k) $ transform. (No default. Other
options are -L, -P and -H).
-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. This option is forced
when using the -K option described above).
-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 radonr) 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 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 least square Fourier (omega,k) Transform:
The conventional Fast Fourier Transform (FFT) assumes
equally spaced seismic traces, in which case the discrete
Fourier transform (DFT) as implemented in program fft2da
becomes orthogonal. In addition, the equal spacing allows
one to exploit a certain spatial invariance of the data
resulting in the 'fast' part of the FFT. For irregularly
sampled data, the inverse DFT as implemented in program fk
is non orthogonal. The least square DFT is closely related
to the least square DRT:
We define the reverse DFT (implemented in program radonr)
for irregularly spaced data to be :
${u}=[F] {U}$,
or
$u ( omega , x sub j ) = sum from {m= k sub min } to {m = k
sub max } F sub mj U sub m $,
where
$F sub {mj} = e sup { +i omega [k sub m x sub j ] }$,
$k sub m = $ the mth wavenumber component in the x direc-
tion,
$k sub min = 2 pi f sub 4 m sub min / x sub max $, and
$k sub max = 2 pi f sub 4 m sub max / x sub max $
We define the forward DFT such that when when it is followed by the above reverse transform we can reconstruct the original data, $u( omega , x sub j )$ in a least squares sense:
${U}=([F] sup * [F] + epsilon [I] ) sup {-1} [F] sup *
{u}$,
where
$F sub {mj} sup * = e sup { -i omega [k sub m x sub j ]
}$.
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 radonf|polymute|radonr, radonf|taupred|radonr and
radonf|polymute|taupred|radonr can be quite effective in
eliminating coherent noise and/or multiples. The flow
radonf|radonr 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 |radonf -semb |radonr|rwindow
-R can be used to significantly increase the coherency of
post stack data.
NOTE 5:
The flow rwindow -ntr200 -npad21
|radonf|taupred|radonr|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 radonf and radonr 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 radonf and radonr 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 radonf and radonr 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 radonf and radonr 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: Geophy-
sics, 46, 255-267.
BUGS
No bugs known at present.
SEE ALSO
radonr,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.
COPYRIGHT
copyright 2001, Amoco Production Company
All Rights Reserved
an affiliate of BP America Inc.
Man(1) output converted with
man2html