NAME
radon3d - filter and interpolate 3D common offset or
stacked data using a 3d @ ( tau ,p,q)@ discrete Radon
transform.
SYNOPSIS
radon3d [ -Nfile_in ] [ -Ofile_out ] [ ildmdx ] [ -cldmdy ]
[ -dzdz ] [ -ilhwwx ] [ -clhwwy ] [ -ilazimxazim ] [ -clazi-
myazim ] [ -ilinterpnxinterp ] [ -clinterpnyinterp ] [ -f1f1
] [ -f2f2 ] [ -f3f3 ] [ -f4f4 ] [ -t1t1 ] [ -t2t2 ] [ -t3t3
] [ -t4t4 ] [ -startlinestartline ] [ -endlineendline ] [
-starttracestarttrace ] [ -endtraceendtrace ] [ -prewprew ]
[ -pminpmin ] [ -pmaxpmax ] [ -qminqmin ] [ -qmaxqmax ] [
-p1p1 ] [ -p2p2 ] [ -p3p3 ] [ -p4p4 ] [ -q1q1 ] [ -q2q2 ] [
-q3q3 ] [ -q4q4 ] [ -smaxsmax ] [ -s1s1 ] [ -s2s2 ] [ -s3s3
] [ -s4s4 ] [ -a1a1 ] [ -a2a2 ] [ -a3a3 ] [ -a4a4 ] [ -F ] [
-R ] [ -reject ] [ -rho ] [ -V ] [ -? ]
DESCRIPTION
radon3d reads in 3D common offset seismic data (including
zero offset, 3d stacked data) and rejects coherent seismic
noise according to user specified ranges of linear moveout
and azimuth using a 3D @ ( tau ,p,q)@ discrete Radon
transform.
ALGORITHM IMPLEMENTATION.
We begin by defining the inverse discrete Radon transform from @ ( omega , p,q) @ space to @ ( omega ,x,y) @ space:
@u = R under U @
where:
@U sub m = U( omega ,p sub m ,q sub m )@,
@u sub j = u( omega , x sub j , y sub j )@
@R sub {mj} = exp [ +i omega (p sub m x sub j +q sub m y sub j ) ]@,
@(x sub j, y sub j ) = @ the relative distance of the jth trace in the analysis window from the reference trace at @(x sub 0, y sub 0 )@,
@omega =@ the temporal frequency in radians/sec,
@p sub m = @ the mth apparent moveout in the inline (x) direction, and
@q sub m = @ the mth apparent moveout in the crossline (y) direction.
This allows us to define the forward Discrete Radon Transform using least squares:
@U = ( R under sup H R under + prew * I under ) sup {-1} R under sup H u@,
where
@R sub {mj} sup H = exp [ -i omega (p sub m x sub j +q sub m y sub j ) ]@,
@I under @ = the identity matrix, and
prew = the prewhitening factor entered on the command line.
We can efficiently mute the data in a space and time
invarient manner directly in the @( omega ,p,q)@ domain. We
define the diagonal weight matrix @ W under sub pq @ as:
@W sub {pq} = 0. ~ if ~ sqrt {p sup 2 + q sup 2 } <= s sub 1@ ,
@W sub {pq} = 1. ~ if ~ s sub 2 <= sqrt {p sup 2 + q sup 2 } <= s sub 3@,
@W sub {pq} = 0. ~ if ~ sqrt {p sup 2 + q sup 2 } >= s sub 4@ , and
@W sub {pq} = ~ tapered @, otherwise.
The weighted (muted) data in the @ ( omega ,p,q)@ domain are then given by:
@V = W under sub pq U sub pq @.
Finally, the weighted (muted) data in the @ ( omega ,x,y)@ domain are given by:
@ v = R under V = R under W under sub pq U @.
This least square modeled data, @ v @ can be efficiently
calculated from the original input data, u , using a
precomputed matrix, @ A under @:
@ v = A under u @
where
@ A = R under W under sub pq ( R under sup H R under + prew I under ) sup {-1} R sup H @.
In the reject mode, we subtract the least square modeled component of the data @ v @ from the original data set, u:
@ v sub r = u - v = (I under - A under ) u .@
Since the dimensions of @R under @ can be several times
larger than that of @A under @ for spatially aliased data,
it is significantly more efficient to run a single pass of
radon3d than to run the piped flow radon3d -F | radon3d -R.
However, if time or space variant @( tau ,p,q) @ filtering
or deconvolution are desired, then the intermediate @( tau
,p,q) @ data set needs to be generated using the -F option
described below. In this case, only those (p,q) components
for which @ W under sub pq @ are nonzero are generated.
COMMAND LINE ARGUMENTS:
radon3d gets all its parameters from command line arguments.
These arguments specify the input, output, spatial analysis
window, and dip discretization parameters.
-N file_in
Enter the input data set name or file immediately after
typing -N. This input file should include the complete
path name if the file resides in a different directory.
Example -N/export/data2/san_juan/time_stack tells the
program to look for file 'time_stack' in directory
'/export/data2/san_juan'. For this program, the data
must be stored as a rectangular grid of regularly
binned data. The number of traces denoted by lineheader
word 'NumTrc' defines the number of traces in the 'x'
direction. The number of records (seismic lines)
denoted by lineheader word 'NumRec' defines the number
of traces in the 'y' direction. Missing data should be
padded in with dead traces flagged by the dead trace
header flag 'StaCor'=30000.
-O file_out
Enter the output filtered data set name or file immedi-
ately 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).
-ildm dx
Enter the inline trace separation distance in m (ft).
(Default: read in dx from seismic line header word
'ILClIn')
-cldm dy
Enter the crossline line separation distance in m (ft).
(Default: read in dy from seismic line header word
'CLClIn')
-dz dz
Enter -dz to designate a 3d depth section (x,y,z) fol-
lowed by the vertical sample spacing in m (ft).
(Default, if -dz is entered but no value supplied, dz
is taken from line header word Dz1000; otherwise, time
domain data are assumed).
-ilhw wx
Enter the half width of the rectangular analysis window
in the inline (x) direction in m (ft). (Default,
wx=5.*dx, such that the x dimension of the analysis
window ranges from (-5:+5) trace about the analysis
point).
-clhw my
Enter the half width of the rectangular analysis window
in the cross line (y) direction in m (ft). (Default,
wy=5.*dy, such that the y dimension of the analysis
window ranges from (-5:+5) trace about the analysis
point).
-f1 f1
Frequency in cycles/sec (Hz) for time data, or in
cycles/km (cycles/kft) for depth data at which we begin
roll in of a Hamming zero phase band pass filter
(default f1=df=1./trace length).
-f2 f2
Frequency in cycles/sec (Hz) for time data, or in
cycles/km (cycles/kft) for depth data at which we end
roll in of a Hamming zero phase band pass filter
(default = f2=f1).
-f3 f3
Frequency in cycles/sec (Hz) for time data, or in
cycles/km (cycles/kft) for depth data at which we begin
roll out of a Hamming zero phase band pass filter
(default = f3=f4).
-f4 f4
Frequency in cycles/sec (Hz) for time data, or in
cycles/km (cycles/kft) for depth data at which we end
roll out of a Hamming zero phase band pass filter. The
cost increases with the value of f4. (default @f4=f
sub Nyquist@ )
-t1 t1
Enter the beginning of the analysis window roll in
taper in ms (m,ft). (Default = t1= first sample of the
trace, or 0 ms).
-t2 t2
Enter the end of the analysis window roll in taper in
ms (m,ft). (Default = t1=t2).
-t3 t3
Enter the beginning of the analysis window roll out
taper in ms (m,ft). (Default = t3=t4).
-t4 t4
Enter the end of the analysis window roll out taper in
ms (m,ft). (Default t1= last sample of the trace).
-startline startline
Enter the first output line to be generated (Default:
startline=1, the first line of the data set).
-endline endline
Enter the last output line to be generated (Default:
endline= the last line of the data set).
-starttrace starttrace
Enter the first output trace to be generated (Default:
starttrace=1, the first trace of the data set).
-endtrace endtrace
Enter the last output trace to be generated (Default:
endtrace= the last trace of the data set).
-ilinterp nxinterp
Enter the number of output traces to be output (inter-
polated) in the inline direction. (Default: nxinterp=1,
or no inline interpolation). (Warning: the interpola-
tion and reject options cannot be used together, since
only the 'noise' have been interpolated to the new,
intermediate locations.
-clinterp nyinterp
Enter the number of output traces to be output (inter-
polated) in the crossline direction. (Default: nyin-
terp=1, or no crossline interpolation) (Warning: the
interpolation and reject options cannot be used
together, since only the 'noise' have been interpolated
to the new, intermediate locations.
Definition of mutes using apparent dips
-pmin pmin
Minimum apparent dip to be modeled in the inline (x)
direction, measured in ms/m or ms/ft for time volumes
(in m/m or ft/ft for depth volumes). (Default: pmin=0.)
-pmax pmax
Maximum apparent dip to be modeled in the inline (x)
direction, measured in ms/m or ms/ft for time volumes
(in m/m or ft/ft for depth volumes). (Default: pmax=0.)
-qmin qmin
Minimum apparent dip to be modeled in the crossline (y)
direction, measured in ms/m or ms/ft for time volumes
(in m/m or ft/ft for depth volumes). (Default: qmin=0.)
-qmax qmax
Maximum apparent dip to be modeled in the crossline (y)
direction, measured in ms/m or ms/ft for time volumes
(in m/m or ft/ft for depth volumes). (Default: qmax=0.)
-p1 p1
Apparent dip in the inline (x) direction (in ms/m or
ms/ft for time volumes; in m/m or ft/ft for depth
volumes) at which we begin the roll in of the Hamming
weights on ray parameter p. (Default, p1=pmin).
-p2 p2
Apparent dip in the inline (x) direction (in ms/m or
ms/ft for time volumes; in m/m or ft/ft for depth
volumes) at which we end the roll in of the Hamming
weights on ray parameter p (Default, p2=p1)
-p3 p3
Apparent dip in the inline (x) direction (in ms/m or
ms/ft for time volumes; in m/m or ft/ft for depth
volumes) at which we begin the roll out of the Hamming
weights on ray parameter p (Default, p3=p4)
-p4 p4
Apparent dip in the inline (x) direction (in ms/m or
ms/ft for time volumes; in m/m or ft/ft for depth
volumes) at which we end the roll out of the Hamming
weights.on ray parameter p (Default, p4=pmax)
-q1 q1
Apparent dip in the crossline (y) direction (in ms/m or
ms/ft for time volumes; in m/m or ft/ft for depth
volumes) at which we begin the roll in of the Hamming
weights on ray parameter q. (Default, q1=qmin).
-q2 q2
Apparent dip in the crossline (y) direction (in ms/m
or ms/ft for time volumes; in m/m or ft/ft for depth
volumes) at which we end the roll in of the Hamming
weights on ray parameter q (Default, q2=q1)
-q3 q3
Apparent dip in the crossline (y) direction (in ms/m
or ms/ft for time volumes; in m/m or ft/ft for depth
volumes) at which we begin the roll out of the Hamming
weights on ray parameter q (Default, q3=q4)
-q4 q4
Apparent dip in the crossline (y) direction (in ms/m
or ms/ft for time volumes; in m/m or ft/ft for depth
volumes) at which we end the roll out of the Hamming
weights.on ray parameter q (Default, q4=qmax)
Definition of mutes using true dip
-smax smax
Maximum true dip to be modeled measured in ms/m or
ms/ft for time volumes (in m/m or ft/ft for depth
volumes). (Default: smax=0.)
-s1 s1
True dip (in ms/m or ms/ft for time volumes; in m/m or
ft/ft for depth volumes) at which we begin the roll in
of the Hamming Apparent dip in the inline (x) direction
filter weights. (Default, s1=0.).
-s2 s2
True dip (in ms/m or ms/ft for time volumes; in m/m or
ft/ft for depth volumes) at which we end the roll in of
the Hamming linear moveout filter weights. (Default,
s2=s1)
-s3 s3
True dip (in ms/m or ms/ft for time volumes; in m/m or
ft/ft for depth volumes) at which we begin the roll out
of the Hamming linear moveout filter weights. (Default,
s3=s4)
-s4 s4
True dip (in ms/m or ms/ft for time volumes; in m/m or
ft/ft for depth volumes) at which we end the roll out
of the Hamming linear moveout filter weights. (Default,
s4=smax)
-a1 a1
Azimuthal angle at which we begin the roll in of the
Hamming azimuth filter weights. If the -ilazim and
-clazim options are properly defined, then an output
azimuth of 0 degrees corresponds to North, while an
output azimuth of 90 degrees corresponds to East.
(Default: a1 = -181 degrees)
-a2 a2
Azimuthal angle at which we end the roll in of the Ham-
ming azimuth filter weights. (Default: a2=a1).
-a3 a3
Azimuthal angle at which we begin the roll out of the
Hamming azimuth filter weights. (Default: a3=a4).
-a4 a4
Azimuthal angle at which we end the roll out of the
Hamming azimuth filter weights. (Default: a4 = +181
degrees)
-ilazim xazim
Enter the azimuth (0 degrees being North, 90 degrees
being East) of the in line (x) axis, where the positive
x axis is the direction of subsequent traces as
presented to the program. This value is used to cali-
brate the azimuth rejection angles a1,a2,a3,a4 above.
(Default: Read xazim from the seismic line header word
'AziIln').
-clazim yazim
Enter the azimuth (0 degrees being North, 90 degrees
being East) of the cross line (y) axis, where the posi-
tive y axis is the direction of subsequent seismic
lines as presented to the program. This value is used
to calibrate the azimuth rejection angles a1,a2,a3,a4
above. (Default: Read yazim from the seismic line
header word 'AziXln').
-prew prew
Enter the prewhitening factor in percent to be used in
the least square matrix inversion associated with the
discrete Radon transform (default = 5.0 %).
-F Enter the command line argument '-F' if you only wish
to perform the 3D forward Discrete Radon Transform.
Input data will be in the (t,x,y) domain and output
data in the @ ( tau ,p,q)@ domain. Warning: be sure to
pipe the output data as they can be many times larger
than the input data!(Default: perform forward
transform, mute in @ ( tau ,p,q)@ domain, and reverse
transform back to (t,x,y) within a single run of
radon3d)).
-R Enter the command line argument '-R' if you only wish
to perform the 3D reverse Discrete Radon Transform.
Input @ ( tau ,p,q)@ data originally generated by the
flow radon3d -F will be muted using the -p1,-p2,-p3,-
p4,-q1,-q2,-q3,-q4 or -s1,-s3,-s3,-s4,-a1,-a2,-a3 and
-a4 options and transformed back to the (t,x,y) domain.
(Default: perform forward transform, mute in @ ( tau
,p,q)@ domain, and reverse transform back to (t,x,y)
within a single run of radon3d).
-reject
If present, reject the least square modeled data lieing
between (p2,p3) and (q2,q3) or (s2,s3) and (a2,a3)
by subtracting it from the input dataset. (Default:
pass data lieing between (p2,p3) and (q2,q3) or
(s2,s3) and (a2,a3)
-rho If present, use the less accurate rho filter convention
@( tau ,p,q)@ transform instead of the discrete Radon
transform.
-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.
NOTE 1:
To suppress the acquisition footprint for subsequent
coherency analysis or post stack depth migration, a 2D by 2D
pass is quite effective. In this example, we have modelled
our strong coherent energy up to a dip of smax=1.0 ms/m. We
don't see dipping reflectors having an apparent dip beyond
.4 ms/m, after which point we begin rejecting coherent
noise. The roll-in is complete beyond .6 ms/m beyond which
all (unaliased) dipping events are rejected. It is extremely
important to keep the values of f1 and f2 as low as possible
when filtering 3d time or depth migrated data, since images
of vertical faults will be mapped to the low frequency
range. The first pass of radon3d is a degenerate 2D
(2*mx+1)=21 point in line running window dip filter. The
second pass is a degenerate 2D (2*my+1)=21 point cross line
running window dip filter. Since the cost is proportional to
the number of points in the operator, 2 passes of a 21 by 1
filter is 20 times faster than a single pass of a 21 by 21
filter. Both will reject the acquisition footprint equally
well.
radon3d -N input_file ildm 25 -cldm 25 -ilhw 250 -clhw 0 -smax 1.0 -s1 .4 -s2 .6 -f1 2 -f2 5 -f3 60 -f4 80 |
radon3d ildm 25 -cldm 25 -ilwx 0 -clhw 250 -smax 1.0 -s1 .4 -s2 .6 -f1 2 -f2 5 -f3 60 -f4 80 -O output_file
NOTE 2:
To perform predicitive deconvolution in the @ ( tau ,p,q)@
domain one might use the following piped flow:
radon3d -N input_file -ildm 25 -cldm 25 -ilhw 50 -clhw 50 -smax 1.0 -f1 5 -f2 10 -f3 60 -f4 80 -F |
pred -p 40 -ol 125 |
radon3d -smax 1.0 -f1 5 -f2 10 -f3 60 -f4 80 -s1 .5 -s2 .6 -R -O output_file
NOTE 3:
To run under IKP, the user needs to supply each processor
with enough records to contain the computational grid. This
is acheived by invoking the roll-along mode in routine
splitr:
splitr -R -nxpadmx -Osubcube1 -Osubcube2 -Osubcube3 -Osubcube4
where @mx >=@ where wx/dx above.
See Also:
IKP, splitr, gather, vfilt3d, swfilt3d
REFERENCES:
Marfurt, K. J. (1995) Filtering and interpolation of 3D com-
mon offset data volumes using spherical discrete Radon
transforms. Geoscience Research Bulletin F-94-xxx. (in
press).
CONTRACT AGREEMENT
This product is a 4th Quarter, 1994, deliverable provided
under Research Contract Agreement D94-2518 (Seismic Signal
Analysis). IKP options added, 2/95. Interpolation and
reject options added, 6/95.
AUTHOR
Kurt. J. Marfurt (1994, EPTG, Tulsa).
COPYRIGHT
copyright 2001, Amoco Production Company
All Rights Reserved
an affiliate of BP America Inc.
Man(1) output converted with
man2html