GEOPHYSICS, VOL. 76, NO. 3 (MAY-JUNE 2011); P. C41–C51, 12 FIGS., 1 TABLE.
10.1190/1.3565182
An accurate ray-based offset-to-angle transform from normal moveout
uncorrected multicomponent data in a transversely isotropic medium with
vertical symmetry axis
Pradip Kumar Mukhopadhyay1 and Subhashis Mallick1
ABSTRACT times. Repetition of this procedure for all zero-offset time
samples and for all desired angles produces an angle gather.
A new ray-based approach for converting multicomponent The procedure allows computation of both primary (P-P) and
normal moveout uncorrected offset domain prestack seismic mode-converted (P-SV) angle gathers. Also, by calculating the
data into angles is presented. From any two-way zero-offset exact raypath length, an accurate geometric spreading com-
time sample and a given angle of incidence, the exact ray and pensation can be computed and applied to these gathers. By
a few neighboring rays that contribute to the given angle of using the methodology for a horizontally stratified transversely
incidence are traced through the medium to their correspond- isotopic medium with a vertical symmetry axis (a VTI me-
ing source and receiver locations to compute their raypath dium), it is possible to extract a much more accurate angle-do-
lengths, traveltimes, and offsets. A polynomial then is fitted main response from the ray-based method than is possible
locally between these computed offsets and their correspond- from a traditional NMO-based method. Accuracy out to large
ing traveltimes. Next the traveltimes for the actual offsets on angles for both primary and mode-converted reflection sug-
data that lie within these traced rays are computed from this gests that this ray-based transform potentially could be used in
local polynomial fit, and their corresponding time samples are a multicomponent prestack waveform inversion scheme.
partially stacked and moved to their two-way zero-offset
INTRODUCTION direct estimates of subsurface reservoir properties (Sen and
Stoffa, 1991, 1992; Mallick, 1995, 1999; Sen and Roy, 2003).
The variations of primary and converted wave (P-P and In all these applications, an accurate offset-to-angle transform
P-SV) reflection amplitudes as functions of the angle of inci- plays one of the most vital roles.
dence can be calibrated to subsurface lithology and reservoir A current approach to an offset-to-angle transform is based pri-
fluids and are the fundamental basis for amplitude-variation- marily on the work by Todd and Backus (1985) and Resnick
with-offset-and-angle (AVO=AVA) analysis (Castagna, 1993). (1993), among others. In this procedure, the angle h for any two-
For isotropic elastic media, the angular variation of reflection way zero-offset time sample t0 is computed from the relation
amplitudes depends strongly on the P-wave-to-S-wave velocity- ( )
ratio (VP=VS) contrasts between the reflecting interfaces
Vint
sin h ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffixffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : (1)
(Ostrander, 1984). Reservoir parameters, such as P- and S-wave Vsmooth x2 þ ðV 2smootht0Þ
reflectivities, Poisson reflectivity, density reflectivity, and the
like, can be estimated from the AVO analysis using a set of In equation 1, x is the source-to-receiver offset, Vsmooth is a spa-
weighted stacks (Stolt and Weglein, 1985; Smith and Gidlow, tially varying velocity function derived by smoothing the stack-
1987) or from an AVO inversion (Kelly et al., 2001; Skidmore ing velocities over a cable length, and Vint is the interval veloc-
et al., 2001). Alternatively, synthetic data computed in the ity obtained from Vsmooth. Note that Vint and Vsmooth both are
angle=ray-parameter domain can be compared with the observed functions of t0. Input prestack seismic data first are compensated
prestack seismic data in a waveform inversion scheme to obtain for spherical divergence loss and corrected for normal moveout
Manuscript received by the Editor 29 October 2009; revised manuscript received 22 December 2010; published online 2 May 2011.
1University of Wyoming, Department of Geology and Geophysics, Laramie, Wyoming, U. S. A. E-mail: pmukhopa@uwyo.edu; smallick@uwyo.edu;
mallick_sub@yahoo.com.
VC 2011 Society of Exploration Geophysicists. All rights reserved.
C41
C42 Mukhopadhyay and Mallick
(Todd and Backus, 1985). Using equation 1, we then compute NMO-uncorrected prestack seismic data. We ray-trace dynami-
the source-to-receiver offset x that corresponds to a given angle cally with a given angle of incidence from each sample point of
h and a given zero-offset time sample t0¼ nDt, where Dt is the t0 to the source and receiver locations, thereby giving ourselves
time-sampling interval and n is the sample number. the flexibility of obtaining the angle gathers for both P-P and P-
For band-limited seismic data, the contribution to a given SV reflections and displaying them together in zero-offset P-
angle does not, however, come from this single offset but wave times.
instead comes from a range of offsets that lie within a quarter Conventional NMO-based angle gathers in the presence of an-
wavelength of the dominant seismic frequency. This concept of isotropy are group angle gathers, and for any quantitative inter-
the range of offsets contributing to the response to a single pretation for reservoir properties, they first must be converted to
angle is similar to the concept of a Fresnel zone, but for an off- their corresponding phase angles (Mallick, 2008). In a ray-based
set-to-angle transform it commonly is called an angle mute pat- offset-to-angle transform, we can compute either phase angle
tern or simply an angle mute (Resnick, 1993). In this paper, gathers or group angle gathers. Lastly, NMO correction is based
therefore, we will use the term “angle mute pattern” or “angle on a Taylor series approximation of the moveout equation (Yil-
mute” to describe this zone. maz, 1987). The ray-based offset-to-angle transform, on the
The seismic trace values for the time sample and for the range other hand, does not make any such approximation. Conse-
of offsets lying within the angle-mute pattern therefore are par- quently, we can compute an accurate traveltime, an anisotropic
tially stacked to obtain the response for a given angle and zero- geometric spreading compensation, and an angle mute pattern
offset time sample (Resnick, 1993). Repeating this operation for for the computed gathers.
all zero-offset time samples and for all desired angles then pro-
duces an angle gather. Although such an offset-to-angle trans- BACKGROUND THEORY
form is for horizontally stratified layers, it can be extended to
include dips, following the steps given by Resnick et al. (1987). For an anisotropic elastic medium, the quasi-P (qP)- and
The above procedures for an offset-to-angle transform for iso- quasi-SV (qSV)-wave velocities vary with the angle (Mallick,
tropic layers also can be extended to VTI layers by replacing 2008). The plane qP- or qSV-wavefields, excited from a point
Vsmooth with the corresponding VTI stacking velocity (Alkhali- source, propagate outward along the surfaces known as the
fah, 1997). However, it must be noted that for anisotropic media slowness surfaces, a 2D cross section of which is shown in
the angle h is the group angle, not the phase angle (Mallick, Figure 1. Any directed line from the source to the slowness sur-
2008). For P-P reflections, the offset-to-angle transform com- face represents a plane wave associated with the slowness vector
puted using the above method produces reasonable results for p and the phase angle h. The phase velocity Vph is also directed
isotropic media. Its accuracy rapidly diminishes with increasing along the phase angle h and has a magnitude equal to the recip-
angles, especially in the presence of anisotropy (Mallick, 2008). rocal of the slowness vector, that is,
Furthermore, for P-S reflections, the NMO equations are com-
1
plex (Thomsen, 1999), and an offset-to-angle transform does not Vph ¼ p̂: (2)
produce reasonable results even for small angles. The primary jpj
purpose of this investigation is to present a ray-based offset-to-
In equation 2, p̂ is the unit vector along the slowness direction
angle transform that can extract an accurate angle-domain
and jpj is the magnitude of the slowness vector. The phase ve-
response for VTI media for both primary (P-P) and mode-con-
locity corresponding to different phase angles also forms a sur-
verted (P-SV) reflections. In contrast to the current NMO-based
face known as the phase-velocity or inverse-slowness surface.
methodology, to which we refer in this paper as the conven-
Whereas the wavefields of constant phase with a given phase
tional approach, input data in our transform are offset-domain
angle h propagate outward along the slowness surface, the
energy propagates along the direction that is normal to the slow-
ness surface at any given point (see Figure 1). In a lossless elas-
tic medium, the group velocity is directed along this normal
with a magnitude equal to the reciprocal of the component of
the slowness vector (Auld, 1973); that is,
n̂
Vg ¼ ; (3)
n̂ p
where n̂ is the unit vector along the direction of the normal to
the slowness surface. Note from Figure 1 that for any given
slowness vector and the corresponding phase-velocity vector
with phase angle h, there is a corresponding group velocity vec-
tor associated with a group angle u. In a homogeneous isotropic
elastic medium, the propagation velocities of P- and SV-waves
are independent of the propagation direction. Consequently, the
slowness surfaces are spheres, the group and phase directions
Figure 1. Pictorial representations of phase slowness vector p, are identical, and the distinction between the phase velocity,
group velocity vector V , phase angle h, and group angle u in an group velocity, phase angle, and group angle is not necessary.g
anisotropic medium. However, this is different for anisotropic media.
Ray-based offset-to-angle transform C43
Assume a VTI medium in which six independent components of The parameters a0 and b0 in equations 7 are the vertical P-wave
the stress tensor s are related linearly to the six independent com- velocity and the SV-wave velocity, respectively. Note that the
p0onen1ts of0the strain tensor e as (Auld, 1973; Thomsen, 119086) 1 phase velocities given by equations 8 and 9 are exact phase velocities in Thomsen notation, without any weak anisotropicB sB xx C B
C11 C11 2C66 CC B 13
0 0 0 exx
s
BB yy CC BB
C CB C approximation.11 2C66 C11 0 0 0 0 CB eyy C
B szz C ¼ B C C C 0 0 0 CCCBBB e Czz CC Recall from Figure 1 that the group velocity vector is normalB C B 13 13 33@ syz A @ 0 0 0 C44 0 0 CAB
:
@ 2eyz CA to the slowness surface and is given by equation 3. For a givenszx 0 0 0 0 C44 0 2ezx phase angle h with slowness vector p (Figure 1), the tangent
sxy 0 0 0 0 0 C66 2exy vector t to the slowness surface is given as (Kreyszig, 1983)
(4) dp
t ¼ : (13)
Five independent elastic moduli for a VTI medium are, there- dh
fore, C11, C13, C33, C44, and C66. The phase velocity for a VTI From Figure 1, the slowness vector p can be written as
medium can be obtained from the solution of the Christoffel
equation (Auld, 1973) and has been given by Thomsen (1986). p ¼ jpjðx̂ sin hþ ẑ cos hÞ; (14)
Because a VTI medium is azimuthally isotropic, it is adequate where jpj is the magnitude of the slowness vector, and x̂ and ẑ
to consider propagation in a vertical (x-z) plane with a given are the unit vectors along the x- and z-axes, respectively. It fol-
phase angle h. The quasi-P (qP)- and quasi-SV (qSV)-wave lows from equations 13 and 14 that
phase velocities VP(h) and VSV(h) then can be expressed as (for
j j d p djpj
details, see Thomsen, 1986) t ¼ x̂ jpj cos hþ sin h ẑ jpj sin h cos h :
dh dh
V2ð Þ ¼ 1h ½C sin2 hþ C 2 (15)P q 11 33 cos hþ C44 þ DðhÞ;2
(5) By using equation 14, we also can write the tangent vector t
2 1V 2 2 and the unit normal vector n̂ in terms of phase velocity asSVðhÞ ¼ ½Cq 11 sin hþ C33 cos hþ C44 DðhÞ;2
1 djVphðhÞj
with qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi t ¼ j ð Þj jVphðhÞj cos h sin h x̂V 2ph h
dh
DðhÞ ¼ ½ðC11 C 244Þ sin hþ ðC C 2 2 djV ðhÞj44 33Þ cos h (6) jVphðhÞj sin hþ cos h ph ẑ ; (16)
þðC13 þ Þ dhC44 sin2 2h:
and
The parameter q in equations 5 is the density. Now, introducing
Thomsenspaffiffirffiffiaffimffiffi eters (Thosmffisffieffiffinffiffiffi, 1986), we have n̂ ¼ sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi1ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi d V ðhÞffiffiffiffiffiffi2
a ¼ C330 ; b0 ¼
C44
; e ¼ C11 C33 ; j
2 phVphðhÞj þq q 2C dh33 (7)
2 2 dj VC phðhÞj
c ¼ 66 C44 ¼ ðC13 þ C44Þ ðC33 Cd 44Þ; ; jVphðhÞj sin hþ cos h x̂dh
2C44 2C33ðC33 C44Þ
and we obtain from equations 5: þ jVphðhÞj d VphðhÞcos h sin h ẑ : (17)dh
V2PðhÞ ¼ a2 1þ e sin2 0 hþ D ðhÞ ; (8)
" # In equations 16 and 17, Vph is the phase velocity of either a qP-and or a qSV-wave. If for Vph we use the expression for qP-wave
2 2
ð Þ ¼ þ a a ð Þ phase velocity V from equation 8 and the expression for qSV-V2 h b2 0 2 0 PSV 0 1 e sin h D h ; (9)b2 b2 wave phase velocity VSV from equation 9, and we plug those0 0 terms into equations 16 and 17, respectively, we will provide
with ( ) tangent and normal directions for the respective wavefields. 1 Finally, from the definition of the group velocity vector given
2ð 1 4d 4ð1þ eÞehÞ ¼ 1 þ 2 h 2 hþ 4 h by equation 3, we obtain, after some algebra,D 1 sin cos sin 1 ;2 12 12
VgðhÞ ¼ j ðhÞj djVphðhÞjVph sin hþ cos h x̂(10) dh
þ j ð Þj djVphðhÞj
Vph h cos h sin h ẑ : (18)
d ¼ 1ð2d eÞ; (11) dh
and From equation 18 it follows that the magnitude of the group ve-
locity vector is given ass ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffib2 2
1 ¼ 1 0 : (12) j ð Þj ¼ j 2 djVphðhÞja2 Vg h VphðhÞj þ ; (19)0 dh
C44 Mukhopadhyay and Mallick
and the group angle u is given as RAY-BASED OFFSET-TO-ANGLE TRANSFORM
þ 1 d VphðhÞ First, consider a stack of horizontal layers for which the depthtan h V ðhÞ dh of each layer, the vertical velocities, and the anisotropic parame-tan½uðhÞ ¼ ph ð Þ : (20) ters are given. Because seismic gathers are expressed in time, tan h d Vph h1 we first convert this depth model to two-way P-P vertical time,VphðhÞ dh as shown in Figure 2. From any given zero-offset time sample
point (j) and for a given group angle (see Figure 2), we calcu-
Now, for a VTI medium, [dD*(h)]=(dh) from equation 10 can
late the group velocity immediately above the sample. We also
be derived as
calculate the ray parameter that corresponds to the given valuedDðhÞ ¼ sin 2h þ ½ð þ Þ 2 of the angle. We then calculate the raypath length, horizontal
dh 2DðhÞ þ d 2 1 e e d sin h : (21)1 offset, and traveltime for the layer above the sample. As the ray
hits an interface with different material properties, we then cal-
Using equation 21, we obtain from equation 8 culate the corresponding phase velocity and angle for the newd V ðhÞ 1a2 sin2h medium by refracting the ray across the interface. After refrac-p ¼ 0 ð Þ tion, we calculate the group velocity and group angle for thedh 2 Vp h
new layer and continue this procedure until we reach the source 1 and receiver locations. By summing the individual traveltimes,eþ 2 2
2Dð Þþ d þ2 1eþ e d sin h ;h 1 raypath lengths, and horizontal offsets from each segment, we
get the total traveltime, raypath length, and horizontal offset.
(22)
Therefore, the ray-based offset-to-angle transform for P-P reflec-
and from equation 9 we have tions for a horizontal stack of layers is
djVSVðhÞj¼1 a
2
0 sin2h 1) For any given zero-offset time sample point j, choose the
dh 2 jV ðhÞj
angle h (phase or group). If that angle is the phase angle,SV compute the phase velocity Vp from equation 8, group veloc- 1e 2ð Þþ d þ2 1eþe d sin2h : ity Vg from equation 19, and group angle u from equation 20.2D h 1 From the phase velocity Vp and the phase angle h, compute
(23) the ray parameter p, given as p¼ sin h=Vp. If the angle is the
j j group angle, then compute the corresponding phase angle rayPlugging the respective expressions for [d VP(h) ]=(dh) and
j j parameter p, and the group velocity V .[d VSV(h) ]=(dh) from equations 22 and 23 into equations 19 g
2) Compute the offset xj at the jth sample point asand 20, respectively, provides the analytical expressions of the
x ¼ ðDt ÞVjj =2 P tanu j, where Dt is the sampling interval and Vgroup velocity and group angle for the qP- and qSV-waves. P
is the vertical P-wave velocity at sample point j. Also calcu-
late the raypath length ‘j and traveltime tj using the group ve-
locity for this sample point (see Figure 2 for details).
3) Repeat steps 1 and 2 to ray trace upward, and calculate offsets
(xj–1, xj–2, …, x1), raypath lengths (‘j–1, ‘j–2, …, ‘1), and trav-
eltimes (tj–1, tj–2, …, t1) for each interface. To apply Snell’s
law at each interface, it is necessary first to calculate the
phase angle from the ray parameter and then to calculate the
phase velocities, the group velocities, and the group angle.
Once these parameters are obtained, it is straightforward to
compute the offsets, raypath lengths, and traveltimes for each
interface.
4) Repeat steps 1 through 3 for all samples and angles to obtain
the desired offset-to-angle transform.
5) If the source and receiver are located oPn the same horizontal1
plane, then the total offPset X is X ¼ 2 i¼j xi, the total ray-
path lePngth L 1is L ¼ 2 i¼j ‘i, and the total traveltime T is
T ¼ 12 i¼j ti. In the case in which the source and receiver
locations are at different depths, we calculate the offset, ray-
path length, and traveltime from the sample point to the
source and receiver locations separately, using steps 1
through 3, and then we add them together to get the total off-
Figure 2. Steps involved in an offset-to-angle transform. For a set, raypath length, and traveltime.
given sample Ppoint and a given angle u, we nupward to the source and receiver locations to obPeed to ray tracetain the total off- Note that for a P-SV reflection, we need to ray trace froml l
set value X ¼ Pi¼j xi, total raypath length L ¼ i¼j ‘i, and total each sample point to the source location by using P-wave prop-l
traveltime T ¼ i¼j ti. We need to carry on this process for all erties and to the receiver location by using SV-wave properties.
time samples and angles to form angle-domain traces. Also note that at every interface, we ray trace using the group
Ray-based offset-to-angle transform C45
velocity and angle, but when we move from one interface to the Because the ti and xi values in equation 24 are known through
next, we need to refract using the corresponding phase velocities ray tracing, the equation can be solved easily for the coefficients
and phase angles. For a given phase angle, it is straightforward c0, c2, and c4. Once these coefficients are known, the traveltimes
to compute the phase velocity, group velocity, and group angle for all the offsets that are within the angle mute pattern between
from equations 8, 9, 19, and 20. For a given group angle, how- un–1 and unþ1 for the input data can be found, partially stacked,
ever, the phase angle, phase velocity, and group velocity must and moved to their corresponding zero-offset times. Repeating
be computed iteratively even for a VTI medium. To do this, for this procedure for all zero-offset time samples and for all
a given group angle u we first assume the corresponding phase desired values of angles gives us the angle gather. Note that our
angle h¼u. Using this initial guess of phase angle, we then local polynomials are in general different for different angles
compute the corresponding group angle and keep modifying the and therefore will handle angle-dependent velocity variations
phase angle using a conjugate gradient procedure (Press et al., correctly. Also note that to demonstrate our procedure of com-
1989) until the computed group angle equals u within an ac- puting traveltimes, in Figure 3 we have shown only one side of
ceptable limit of tolerance. Once the phase angle is obtained, the ray tracing from a given sample point to the receiver loca-
we then can compute the phase and group velocities from equa- tions. Clearly, however, to obtain the traveltimes and offsets
tions 8, 9, and 19. shown in equation 24, it is necessary to perform ray tracing
For dipping layers, the dip angles must be taken into consid- from the sample point to the source location also.
eration in addition to the ray parameter to correctly refract rays A crucial step in our methodology is the accuracy of the
across each interface. In the presence of dipping isotropic layers, polynomial fitting that is used to get the exact time samples.
the same procedure can be used as long as we take the local The procedure described above in equation 24 is essentially a
dips into consideration and ray trace from the given zero-offset three-term polynomial fitting and could be achieved using a
sample point to the source and receiver locations separately. Lagrange polynomial (Press et al., 1989). For practical applica-
However, for dipping anisotropic layers that include VTI, the tions, however, we believe that a constrained cubic spline inter-
procedure is not straightforward. Not only do we need to take polation (Press et al., 1989) is more accurate than a Lagrange
the local dips into consideration, but we also need to consider polynomial fitting. In the examples shown in this paper, we
the apparent change in the elastic stiffness matrix by rotating therefore used a constrained cubic spline interpolation to calcu-
the medium by the dip angles using a Bond transformation ma- late the traveltimes for the range of offsets described as above.
trix (Auld, 1973). Extension of this methodology to dipping ani- For best results, it is necessary also that the given offset values
sotropic layers is being investigated currently. be spread uniformly within each angle mute. To achieve that,
Because we work with NMO-uncorrected gathers we not only we choose the neighboring angle values while ray tracing in
need to obtain the set of offsets within the angle mute pattern such a way that the calculated offsets fall evenly on either side
that contribute to a given angle but also their corresponding of the zone. In addition, we also have found that using five
traveltimes. For a given value of group angle un, we therefore closely spaced rays, two on either side of the desired angle,
use the above procedure to calculate the traveltime, horizontal gives the optimum result. An exact calculation of the angle
offset, and raypath length. To partially stack the traces in the mute pattern is extremely important for the accuracy of our
angle mute pattern, we also calculate horizontal offsets and trav-
eltimes for a few more angles that are close to un and that lie
within the angle mute pattern (e.g., un–1 and unþ1). As is shown
in Figure 3, the raypaths PB and PC are such neighboring ray-
paths of PA. As we ray-trace from a given zero-offset time sam-
ple, the computed offset values OA, OB, and OC for these
angles will not necessarily fall on any given offset sample on
the input data. This is also true for the calculated traveltimes
BB1, AA1, and CC1. However, we can use these calculated off-
sets (OA, OB, and OC) and corresponding traveltimes (BB1,
AA1, and CC1) and locally fit a polynomial to find the travel-
times for the range of offsets that are present in input data and
that lie within the angle mute pattern.
To be specific, assume that for a given angle un from a given
zero-offset time sample P in Figure 3, we ray trace and compute
the traveltimes, and also the offsets for angles un–1, un, and
unþ1. Assume too that the calculated traveltimes BB1, AA1, CC1
and the corresponding offsets OB, OA, and OC are given respec-
tively as tn–1, tn, tnþ1 and xn–1, xn, and xnþ1. We then assume that
the time-distance relation in the neighborhood of un is expressed Figure 3. Steps involved in computing an offset-to-angle trans-
as a polynomial of the type t2¼ c0þ c 2 42 x þ c4 x , and we write form by fitting a polynomial around the given angle un. For a
given two-way zero-offset traveltime (point P) and a given angle
t2 2 4n1 ¼ c0 þ c2xn1 þ c4xn1 (un), the corresponding raypath length (PA), traveltime (AA1), off-
2 ¼ þ 2 þ 4 set (OA), and the angle mute pattern (MN) are shown. For neigh-tn c0 c2xn c4xn (24) boring angle values (un–1 and unþ1) within this angle mute
2 ¼ þ 2 þ 4 pattern, the traveltimes and offsets that are required to fit a poly-tnþ1 c0 c2xnþ1 c4xnþ1: nomial locally are also shown.
C46 Mukhopadhyay and Mallick
method (Monk, 2010). Note that the angle mute pattern is a Figure 6 shows the P-P angle gathers. Figure 6a is the angle
function of offset, angle, depth, and dominant seismic fre- gather computed from the synthetic offset gather using a con-
quency, and an incorrect calculation of the angle mute pattern ventional method — that is, after NMO correction of the verti-
can cause destructive interference to the computed angle-domain cal component of the response, computation of the range of off-
responses. sets for the given angles from equation 1, and partial stacking
of the traces over the angle mute pattern, as described by Todd
EXAMPLES WITH DISCUSSION and Backus (1985) and Resnick (1993). Because our input
model is VTI, we replaced the isotropic stacking velocities in
Figure 4 shows a multilayer anisotropic model in two-way P- equation 1 with their corresponding anisotropic stacking veloc-
wave time and depth. This anisotropic model is derived from a ities (see Alkhalifah, 1997 for details). Figure 6b is the angle
real well-log and anisotropic velocity analysis of prestack seis- gather computed by convolving the source wavelet with an
mic data. Figure 5 shows the vertical and radial component of exact (Zoeppritz) P-P reflection series for different angles, and
the response in offset domain, computed for this model. These Figure 6c is the new ray-based angle gather. Note that the con-
responses were computed using a wave-equation-based modeling ventional and ray-based angle gathers are computed from the
algorithm for a stratified anisotropic medium (for details, see same input model shown in Figure 4, and the Zoeppritz gather
Mallick and Frazer, 1987, 1988, 1990, 1991). is shown for reference to common practice.
Figure 7 shows the converted wave angle gathers for the
same model. Only the ray-based gather is compared with the
Zoeppritz gather in this figure. For clarity, the converted wave
angle gathers of Figure 7 were magnified by a factor of ten
compared with the P-P gathers of Figure 6. Notice that P-SV
angle gathers are shown in two-way P-P time, not in P-SV time.
The example shown in Figures 6 and 7 demonstrates the su-
periority of the ray-based offset-to-angle transformation com-
pared with the results from a conventional approach. Notice that
the P-P events in the conventional gather are bent upward with
increasing angles (Figure 6a). Additionally, the reflection ampli-
tudes in the conventional angle gather tend to dim with increas-
ing traveltime. Although such bending of events is controlled by
the choice of the VTI parameters e and d, the main reason for
such bending and dimming of amplitudes is the increasing inac-
curacy of the NMO equation for large offsets. NMO equations
for the mode-converted reflections are complex even for an iso-
tropic medium. By accurately computing the mode-converted
Figure 4. Multilayer anisotropic earth model in P-P time and
depth, derived from a real well log.
Figure 6. Angle-domain primary (P-P) seismograms computed
using the model shown in Figure 4 and from the offset-domain
vertical-component response shown in Figure 5a. (a) P-P angle
Figure 5. (a) Vertical component and (b) radial component of the gather computed using a conventional approach, (b) P-P Zoeppritz
synthetic seismic response, computed from the model shown in angle gather, and (c) P-P angle gather using the new ray-based
Figure 4. method.
Ray-based offset-to-angle transform C47
traveltimes, not only have we computed accurate angle gathers from the gathers shown in Figures 6 and 7. Figure 8 shows a
but we also have shown a convenient way of referencing them simple five-layer VTI model and Figure 9 shows the computed
to P-wave time (Figure 7) so that they can be compared jointly vertical (9a) and radial component (9b) synthetic seismograms
with P-P reflections in an AVO=AVA analysis. By computing for that model. Figures 10 and 11 show, respectively, P-P and
the exact traveltimes through ray tracing, we are able to P-SV angle gathers extracted from those synthetic responses.
improve the quality of the computed angle gathers drastically Figures 12 and 13 show the extracted amplitudes respectively
(Figures 6c, 7b). of all four P-P and P-SV reflection events seen in Figures 10
In Figures 6 and 7, notice that even though we used exactly and 11.
the same wavelet, the Zoeppritz gathers (Figures 6b and 7a) The five-layer examples shown in Figures 10 and 11 allow
appear to have higher frequency content than the wave equation direct comparison of the extracted amplitudes with Zoeppritz
gathers have (Figures 6a, 6c, 7b), especially with increasing amplitudes. Note that for P-P reflections, shown in Figure 12,
time. The primary reason for this is the fact that the wave equa- the extracted reflection amplitudes from the conventional angle
tion angle gathers, irrespective of whether they are conventional gather deteriorate with increasing time as a result of the pres-
or ray-based, are obtained from the offset-domain synthetics. ence of anisotropy in the model. Even though Event 4 (Figure
They are computed using a wave equation methodology that 12d) appears to be correct for the conventional gather, note that
included all effects of wave propagation, such as interbed multi- this reflection is very weak for large angles, and the apparent
ple reflections and mode conversions. These propagation effects proximity of the conventional angle gather reflections to the
in a finely layered medium cause an apparent loss of frequency, Zoeppritz reflections for this event results from the fact that the
known as the O’Doherty-Anstey effect or the stratigraphic-filter-
ing effect (for details, see Banik et al., 1985a and Banik et al.,
1985b). Zoeppritz gathers, on the other hand, are computed by
convolving the angle-dependent reflection coefficients with the
given wavelet and therefore do not contain any such effect.
Consequently, the frequency content is higher in the Zoeppritz
gather than in the wave equation angle gathers.
Because of wavefield interference effects, demonstrating the
accuracy of the ray-based offset-to-angle transform is difficult
Figure 8. A simple five-layer VTI model in P-P time and depth.
Figure 7. Angle-domain mode-converted (P-SV) gathers com-
puted using the model shown in Figure 4 and from the offset-do-
main radial-component response shown in Figure 5b. (a)
Zoeppritz angle gather and (b) ray-based angle gather. Note that Figure 9. (a) Vertical-component and (b) radial-component syn-
both these angle gathers are displayed in P-wave time. thetic seismic responses for the model shown in Figure 8.
C48 Mukhopadhyay and Mallick
reflection amplitudes themselves are nearly zero. The ray-based tions for events 1 and 4. For events 2 and 3, however, the match
P-P gathers, on the other hand, are very close to the Zoeppritz is good only to 50. That is because in the offset gather (Figure
reflections for all four events shown in Figure 12. 9b), the radial component of response at large offsets is affected
For mode-converted reflections (Figure 13), notice that the by interference from the multiply reflected P-wave modes for
ray-based amplitudes match very closely with Zoeppritz reflec- these reflections. This interference is apparent in the P-SV angle
gather for these events beyond 50.
Figures 14 and 15 show the sensitivity of the ray-based off-
set-to-angle transform to the input anisotropic parameters. From
the five-layer anisotropic model of Figure 8, we generated two
additional models: (1) an isotropic model created by setting
e¼ d¼ 0 for the entire model, and (2) a 10% anisotropic model
created by setting e¼ 0.10 for maximum anisotropy in the
model of Figure 8 and leaving d values as they are. Figures 14
and 15 compare the conventional and ray-based angle gathers
Figure 10. Angle-domain primary (P-P) seismograms computed
using the model shown in Figure 8 and the vertical component
response shown in Figure 9a. (a) P-P angle gather computed using
a conventional approach, (b) P-P Zoeppritz angle gather, and (c)
P-P angle gather using the ray-based method.
Figure 12. Comparisons of P-P reflection amplitudes for all four
reflection events of Figure 10. Amplitudes for (a) event 1, (b)
event 2, (c) event 3, and (d) event 4.
Figure 11. Angle-domain mode-converted (P-SV) gather com-
puted using the model shown in Figure 8 and the radial-compo- Figure 13. Comparisons of P-SV reflection amplitudes for all
nent response shown in Figure 9b. (a) Zoeppritz angle gather and four reflection events of Figure 11. Amplitudes for (a) event 1, (b)
(b) ray-based angle gather. event 2, (c) event 3, and (d) event 4.
Ray-based offset-to-angle transform C49
for these two new models and the original model of Figure 8. In phase angle and the P-P and P-SV reflection coefficients for
subsequent discussion, we refer to these three models as (1) the model 1 are shown respectively in Figure 16a, b, and c. Note
isotropic model, (2) the 10% anisotropic model, and (3) the that these reflection coefficients are exact anisotropic reflection
20% anisotropic model. Figure 14a is the conventional angle-
gather for the isotropic model, Figure 14b is the conventional
angle-gather for the 10% anisotropic model, and Figure 14c is
the conventional angle-gather for the 20% anisotropic model. In
Figure 15a, b, and c we have the respective ray-based angle
gathers for the same set of models.
The conventional NMO-based offset-to-angle transform pro-
vides reasonable results for isotropic media, but its accuracy in
anisotropic media diminishes rapidly with increasing angle and
with increasing anisotropy. From Figure 14, note that the con-
ventional angle gathers are increasingly inaccurate as anisotropy
increases in the model. Ray-based angle gathers (Figure 15), on
the other hand, are not. The accuracy of ray-based angle gathers
out to large angles both for primary (P-P) and for mode-con-
verted (P-SV) reflections is the result of the accuracy of the
transform in computing the large-angle traveltimes. This accu-
racy out to large angles in fact places a strong constraint on
model estimation and potentially be used in a multicomponent
waveform inversion. The use of ray-based angle gathers for
such an inversion currently is being investigated and will be dis-
cussed in a separate paper.
Finally, note that the ray-based transform can properly
Figure 15. Same sensitivity analysis as Figure 14, but using the
account for the difference between the phase angles and the ray-based offset-to-angle transform.
group angles. It is important to be aware of the differences
between the phase and group velocities and angles when one is
dealing with anisotropic AVO analyses; otherwise, any quantita-
tive interpretation for reservoir properties will lead to erroneous
results. Figure 16 is a plot of group angle versus phase angle
and reflection coefficients for two sets of models. Both models
consist of a top shale layer and a bottom sand layer. Elastic con-
stants for these layers were taken from the real well data of Fig-
ure 4 and are given in Table 1. The group angle versus the
Figure 16. Relation between group angle and phase angle and P-P
and P-SV reflection coefficients for two anisotropic models
extracted from real well data. Both of the models consist of a top
shale layer and a bottom sand layer. The model parameters are given
in Table 1. (a) Dependence of phase angle on group angle for model
1. The solid curve is the phase angle as a function of group angle,
and the dashed curve is shown as a reference in which the phase
angle is assumed to be equal to the group angle. (b) P-P reflection
coefficient as a function of group angle for model 1. The solid curve
is computed by transforming the group angles to their corresponding
phase angles and then computing the reflection coefficients for those
Figure 14. Analysis of the sensitivity of the conventional offset- phase angles, whereas the dashed curve is computed by assuming
to-angle transform to the input anisotropic model. Conventional that the group angles are equal to the phase angles. (c) The same as
offset-to-angle transform using (a) the isotropic model, (b) the (b), but for the P-SV reflection coefficient as a function of group
10% anisotropic model, and (c) the 20% anisotropic model. angle for model 1. (d-f) The same as (a-c), but for model 2.
C50 Mukhopadhyay and Mallick
Table 1. Models used for computing the group angle and tion of angle gathers. It also allows computation of converted
phase angle and the P-P and P-SV reflection coefficients for wave angle gathers in P-wave times, thereby offering a conven-
two anisotropic models, shown in Figure 16. ient way to jointly analyze primary and mode-converted reflec-
tion responses in AVO=AVA applications. Because the ray-
Model 1 Model 2 based approach accurately models large-angle traveltime, in
principle it places a stronger constraint on model estimation.
Model Top Bottom Top Bottom Therefore, the ray-based approach potentially can be used to
parameters layer layer layer layer compute angle gathers effectively from observed seismic data
and to match them iteratively with synthetic angle gathers in a
VP0 (km=s) 2.581 2.372 2.170 2.295 multicomponent prestack waveform inversion scheme.
VS0 (km=s) 1.269 1.474 0.888 1.209
q (gm/cm3) 2.2 2.19 2.26 2.2
ACKNOWLEDGMENTS
e 0.12 0.002 0.2 0.002
d 0.01 0.0002 0.02 0.0002 We thank Phil Anno, the SEGAssociate Editor, and three anony-
mous reviewers, whose constructive comments greatly helped
improve the quality of this paper. This research is funded by the
coefficients computed using the procedure outlined by Mallick U.S. Department of Energy Grant DE-NT0004730 and DE-FE-
and Frazer (1991). Because group angles are the ones that deter- 0001160.
mine the traveltimes of seismic reflections, we chose to use the
group angle as the independent variable and to show the phase REFERENCES
angle and reflection coefficients as functions of the group angle.
Although the anisotropy for the first model is relatively weak Alkhalifah, T., 1997, Velocity analysis using nonhyperbolic moveout in
(12%), significant differences still exist between the group transversely isotropic media: Geophysics, 62, 1839–1854.
Auld, B. A., 1973, Acoustic waves in solids: 1, John Wiley and Sons.
angles and phase angles for this model, as seen in Figure 16a. Banik, N., I. Lerche, J. R. Resnick, and R. T. Shuey, 1985a, Stratigraphic
Whereas the solid curve in Figure 16a shows the group-angle filtering, Part II: Model spectra: Geophysics, 50, 2775–2783,
versus phase-angle relation for model 1, for convenience we doi:10.1190/1.1441898.
Banik, N., I. Lerche, and R. T. Shuey, 1985b, Stratigraphic filtering, Part
also have included in this figure a dashed line along which the I: Derivation of the O’Doherty-Anstey formula: Geophysics, 50, 2768–
phase angle is equal to the group angle. Notice that when the 2774, doi:10.1190/1.1441897.
Castagna, J. P., 1993, AVO analysis — Tutorial and review, in J. P. Casta-group angle is 60 , the corresponding phase angle is approxi-
gna and M. M. Backus, eds., Offset-dependent reflectivity: Theory andmately 52 for model 1. In the reflection-coefficient curves of practice of AVO analysis: SEG Investigations in Geophysics No. 8, 3–36.
Figure 16b and c, the curves with solid lines are the reflection Kelly, M., C. Skidmore, and D. Ford, 2001, AVO inversion, Part I: Isolat-
ing rock property contrasts: The Leading Edge, 20, 320–323,
coefficients computed by first converting the group angles to doi:10.1190/1.1438940.
their corresponding phase angles and then computing the reflec- Kreyszig, E., 1983, Advanced engineering mathematics: John Wiley and Sons.
tion coefficients for those phase angles. The dashed curves, on Mallick, S., 1995, Model-based inversion of amplitude-variation-with-off-
set data using a genetic algorithm: Geophysics, 60, 939–954,
the other hand, were computed by assuming that the group doi:10.1190/1.1443860.
angles are equal to the phase angles. The computed results for ——, 1999, Some practical aspects of prestack waveform inversion using
a genetic algorithm: An example from the east Texas Woodbine gas
the second model are shown in Figure 16d, e, and f. Notice that sand: Geophysics, 64, 326–336, doi:10.1190/1.1444538.
the anisotropy in the top shale layer in the second model is 20% ——, 2008, Interpretation of angle gathers for transversely isotropic me-
and significant differences exist between the phase angles and dium: 78th Annual International Meeting, SEG, Expanded Abstracts,
2937–2941.
the group angles. In real seismic data, although our extracted Mallick, S., and L. N. Frazer, 1987, Practical aspects of reflectivity model-
angles are group angles, the reflection amplitudes will be pro- ing: Geophysics, 52, 1355–1364, doi:10.1190/1.1442248.
portional to the reflection coefficients for their corresponding ——, 1988, Rapid computation of multioffset vertical seismic profile syn-
thetic seismograms for layered media: Geophysics, 53, 479–491,
phase angles. In other words, these reflection amplitudes will be doi:10.1190/1.1442479.
proportional to the solid curves shown in Figure 16b, c, e, and ——, 1990, Computation of synthetic seismograms for stratified azimu-
thally anisotropic media: Journal of Geophysical Research, 95, B6,
f. In AVO=AVA analysis, on the other hand, if we ignore this 8513–8526, doi:10.1029/JB095iB06p08513.
fact and assume the group angle to be equal to the phase angle, ——, 1991, Reflection=transmission coefficients and azimuthal anisot-
the reflection amplitudes for the same models will follow the ropy in marine seismic studies: Geophysical Journal International, 105,
no. 1, 241–252, doi:10.1111/j.1365-246X.1991.tb03459.x.
dashed curves. Consequently, an AVO=AVA-based inversion of Monk, D. J., 2010, Fresnel-zone binning: Fresnel-zone shape with offset and
the angle-domain data in which the difference between group velocity function: Geophysics, 75, no. 1, T9–T14, doi:10.1190/1.3294576.
angles and phase angles is ignored will result in incorrect reser- Ostrander, W. J., 1984, Plane-wave reflection coefficients for gas sands at
nonnormal angles of incidence: Geophysics, 49, 1637–1648,
voir properties. doi:10.1190/1.1441571.
Press, W. H., B. P. Flannery, S. A. Teukolsky, and W. T. Veterling, 1989,
Numerical recipes: Cambridge University Press.
Resnick, J. R., 1993, Seismic data processing for AVO and AVA analysis,
CONCLUSIONS in J. P. Castagna and M. M. Backus, eds., Offset-dependent reflectivity:
Theory and practice of AVO analysis: SEG Investigations in Geophy-
Computing an accurate offset-to-angle transform is crucial to sics No. 8, 175–189.
Resnick, J. R., P. Ng, and K. Larner, 1987, Amplitude versus offset analy-
predicting subsurface lithology and fluid properties from seismic sis in the presence of dip: 57th Annual International Meeting, SEG,
data. Whereas a conventional NMO-based transform is increas- Expanded Abstracts, 617–620.
Sen, M. K., and I. G. Roy, 2003, Computation of differential seismograms
ingly inaccurate with increasing angles, especially in the pres- and iteration adaptive regularization in prestack waveform inversion:
ence of anisotropy, a ray-based method allows accurate calcula- Geophysics, 68, 2026–2039, doi:10.1190/1.1635056.
Ray-based offset-to-angle transform C51
Sen, M. K., and P. L. Stoffa, 1991, Nonlinear one-dimensional seismic Stolt, R. H., and A. B. Weglein, 1985, Migration and inversion of seismic
waveform inversion using simulated annealing: Geophysics, 56, 1624– data: Geophysics, 50, 2458–2472, doi:10.1190/1.1441877.
1638, doi:10.1190/1.1442973. Thomsen, L., 1986, Weak elastic anisotropy: Geophysics, 51, 1954–1966,
— —, 1992, Genetic inversion of AVO: The Leading Edge, 11, 27–29, doi:10.1190/1.1442051.
doi:10.1190/1.1436845. ——, 1999, Converted-wave reflection seismology over inhomogeneous,
Skidmore, C., M. Kelly, and R. Cotton, 2001, AVO inversion, Part 2: Iso- anisotropic media: Geophysics, 64, 678–690, doi:10.1190/1.1444577.
lating rock property contrasts: The Leading Edge, 20, 425–428, Todd, C. P., and M. M. Backus, 1985, Offset-dependent reflectivity in a
doi:10.1190/1.1438966. structural context: 55th Annual International Meeting, SEG, Expanded
Smith, G. C., and P. M. Gidlow, 1987, Weighted stacking for rock prop- Abstracts, 586–588.
erty estimation and detection of gas: Geophysical Prospecting, 35, 993– Yilmaz, O., 1987, Seismic data processing: SEG Investigations in Geo-
1014, doi:10.1111/j.1365-2478.1987.tb00856.x. physics No. 2.