Frequency domain reverse time migration for acoustic-elastic coupled media using the wavefield separation method

© 2025 by the Authors. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution -Noncommercial 4.0 International License (CC-by the license) ( https://creativecommons.org/licenses/by-nc/4.0/ )

Download PDF

Cite

XML

Abstract

Kang, S.-G., Ha, W., Shin, C., Hong, J.K. and Jin, Y.K., 2016. Frequency domain reverse time migration for acoustic-elastic coupled media using the wavefield separation method. Journal of Seismic Exploration, 25: 57-85. Recent research results concerning frequency domain reverse time migration based on the adjoint-state of the acoustic wave equation have highlighted several limitations imposed by the use of an acoustic-based algorithm. In marine seismic exploration, targeted areas are located within elastic media. Elastic wave components, such as S-waves, surface waves and mode converted waves can remain obscured by reverse time migration based on the acoustic wave equation in the targeted media. Several research papers addressing the topic of acoustic-elastic coupled media with full waveform inversion have shown that this method can generate more accurate inversion results for P-wave velocity models than acoustic-based algorithms. This paper formulates the frequency domain reverse time migration for acoustic-elastic coupled media based on the adjoint-state of the wave equation. It goes on to adopt the wavefield separation method to reduce the effects of crosstalk artifacts on migrated images. The validity of the proposed algorithm is demonstrated using a synthetic dataset generated by elastic staggered grid time modeling. The image of the frequency domain reverse time migration for acoustic-elastic coupled media calculated using the wavefield separation method is then compared to the results of the acoustic reverse time migration and reverse time migration for acoustic-elastic coupled media using a conventional zero-lag cross-correlation approach. Comparison of the migration images revealed that the images of acoustic-elastic coupled media from the wavefield separation method resolved geological structures with greater accuracy and exhibited fewer noise-contaminated components than those obtained from the acoustic and conventional acoustic-elastic coupled media imaging methods. We also analyze the reverse time migration method’s sensitivity to correct P- and S-wave inputs and density model assumptions.

Keywords

frequency domain

reverse time migration

acoustic-elastic coupled media

adjoint-state wave equation

marine seismic data

wavefield separation

References

Bae, H.S., Shin, C., Cha, Y.H., Choi, Y. and Min, D., 2010. 2D acoustic-elastic coupled
waveform inversion in the LaPlace domain. Geophys. Prosp., 58: 997-1010.
Baysal, E., Kosloff, D.D. and Sherwood, J.W.C., 1983. Reverse time migration. Geophysics, 48:
1514-1524.
Choi, Y., Min, D.-J. and Shin, C., 2008. Two-dimensional waveform in- version of multicomponent
data in acoustic-elastic coupled media. Geophysics, 56: 863-881.
Chung, W., Pyun, S., Bae, H.S., Shin, C. and Marfurt, K.J., 2012. Implementation of elastic
reverse-time migration using wavefield separation in the frequency domain. Geophys. J. Int.,
189: 1611-1625.
Cohen, G.C., 2002. High-Order Numerical Methods for Transient Wave Equations. Springer
Verlag, Berlin.
Kang, S.-G., Bae H.S. and Shin, C., 2012. Laplace-fourier-domain waveform inversion for
fluid-solid media. Pure Appl. Geophys., 169: 2165-2179.
REVERSE TIME MIGRATION 85
Kim, M.H., Choi, Y., Cha, Y.H. and Shin, C., 2009. 2-D frequency-domain waveform inversion
of coupled acoustic-elastic media with an irregular interface. Pure Appl. Geophys., 166:
1967-1985.
Kim, Y., Min, D. and Shin, C., 2011. Frequency-domain reverse-time migration with source
estimation. Geophysics, 76: $41-S49.
Komatitsch, D., Barnes, C. and Tromp, J., 2000. Wave propagation near a fluid-solid interface: A
spectral-element approach. Geophysics, 65: 623-631.
Lailly, P., 1983. The seismic inverse problem as a sequence of before stack migrations. Expanded
Abstr., Conf. Inverse Scattering: Theory and Application Mathematics: 206-220.
Lee, H., Lim, S., Min, D., Kwon, B. and Park, M., 2009. 2D time-domain acoustic-elastic coupled
modeling: a cell-based finite-difference method. Geosci. J., 13: 407-414.
Lee, J., Kim, Y. and Shin, C., 2012. Frequency-domain reverse time migration using the L,-norm.
J. Seismic Explor., 21: 281-300.
Loewenthal, D. and Mufti, I.R., 1983. Reversed time migration in spatial frequency domain.
Geophysics, 48: 627-635.
Martin, G.S., Marfurt, K.J. and Larsen, S. 2002. Marmousi-2: an updated model for the
investigation of AVO in structurally complex areas. Expanded Abstr., 72nd Ann. Internat.
SEG Mtg., Salt Lake City: 1979-1982.
McMechan, G.A., 1983. Migration by extrapolation of time-dependent boundary values. Geophys.
Prosp., 31: 413-420.
Plessix, R.-E. and Mulder, W.A., 2004. Frequency-domain finite-difference amplitude preserving
migration. Geophys. J. Internat., 157: 975-987.
Pratt, R.G., Shin, C. and Hicks, G.J., 1998. Gauss-Newton and full Newton method in frequency
domain seismic waveform inversion. Geophys. J. Internat., 133: 341-362.
Shin, C., Jang, S. and Min, D., 2001. Improved amplitude preservation for prestack depth migration
by inverse scattering theory. Geophys. Prosp., 49: 592-606.
Shin, C., Min, D.J., Yang, D. and Lee, S.K., 2003. Evaluation of poststack migration in terms of
virtual source and partial derivative wavefields. J. Seismic Explor., 12: 17-37.
Symons, N.P., Aldridge, D.F. and Haney, M.M., 2006. 3D acoustic and elastic modelling with
Marmousi-2. Expanded Abstr., 76th Ann. Internat. SEG Mtg., New Orleans: 2171-2175.
Tarantola, A., 1984. Inversion of seismic reflection data in the acoustic approximation. Geophysics,
49: 1259-1266.
Virieux, J., 1986. P-SV wave propagation in heterogeneous media: velocity-stress finite-difference
method. Geophysics, 51: 889-901.
Whitmore, N.D., 1983. Iterative depth migration by backward time propagation. Expanded Abstr.,
53rd Ann. Internat. SEG Mtg., Las Vegas: 382-385.
Xu, K., Zhou, B. and McMechan, G.A., 2010. Implementation of prestack reverse time migration
using frequency-domain extrapolation. Geophysics, 75: S61-S72.
Zienkiewicz, O.C., Taylor, R.L. and Zhu, J.Z., 2005. The Finite Element Method: Its Basis and
Fundamentals. Elsevier, Butterworth-Heinemann, New York.
JOURNAL OF SEISMIC EXPLORATION 25, 87-101 (2016) 87
A LINEARIZED GROUP VELOCITY APPROACH FOR
TWO-POINT qP RAY TRACING IN A LAYERED
ORTHORHOMBIC MEDIUM
P.F. DALEY and E.S. KREBES
Department of Geoscience, University of Calgary, 2500 University Drive, Northwest Calgary,
Alberta, Canada T2N 1N4. pfdaley@gmail.com, krebes@ucalgary.ca
(Received August 20, 2015; revised version accepted November 27, 2015)
ABSTRACT
Daley, P.F. and Krebes, E.S., 2016. A linearized group velocity approach for two-point qP ray
tracing in a layered orthorhombic medium. Journal of Seismic Exploration, 25: 87-101.
Using a linearized approximation for the quasi-compressional phase velocity, vip in an
orthorhombic anisotropic medium, which is a subset of the related quasi-compressional (qP) wave
propagation in a general 21 parameter anisotropic medium, a linearized compressional group velocity
may be derived as a function of group angles only. In addition, linearized analytic expressions for
the components of the slowness vector in terms of group velocities and angles are also obtained.
These expressions are used to define two nonlinear equations which are a generalization of Snell’s
Law. The solutions of these are used to determine the propagation directions of the reflected and
transmitted rays due to an incident ray at an interface between two orthorhombic media. The axes
of anisotropy, in both media are, in general, not aligned with the interface separating them.
Computer code has been written to consider ray tracing in media defined by a type of large scale
3D finite element blocking (blocky structures). However, a plane parallel layered medium will be
used in preliminary investigations. Additionally, in each of these elements (layers) the anisotropic
parameters are initially assumed to be constant, although provisions for at least minimal spatial
variations of the anisotropic parameters have been considered.
KEY WORDS: orthorhombic media, two-point ray tracing, group velocity,
layered anisotropic media.
INTRODUCTION
In the geophysical literature related to wave propagation in anisotropic
media, specifically quasi-compressional qP waves in a medium displaying
orthorhombic symmetry, a linearized approach to derive an approximate phase
velocity expression for quasi-compressional qP wave propagation has been
presented in Backus (1965) and more recently by Pseneik and Farra (2005).
88 DALEY & KREBES
These and other methods, such a perturbation theory, have been employed
to formulate expressions for quantities related to the propagation of qP waves
as well as for the two shear wave modes, qS, and qS,, in a general 21 parameter
anisotropic medium (Every, 1980; Every and Sachse, 1992; Jech and PSenéik,
1989; Pseneik and Gajewski, 1998; Song and Every, 2000; Song et al., 2001;
Pseneik and Farra, 2008). Once linearized phase velocity approximations have
been obtained, eikonal equations with comparable accuracy may be written.
From these, using the method of characteristics, (Courant and Hilbert, 1962;
Cerveny, 2002) the formulae for the vector components of group velocity may
be obtained.
Employing an extension of this method of approximation, Daley and
Krebes (2006) obtained expressions for the qP group velocities in an
orthorhombic anisotropic medium, expressed in terms of group angles. It was
further shown that for a weakly anisotropic medium the average deviation of this
approximation from the exact expression for the qP group velocity, in an
orthorhombic medium, over a range of polar angles for a number of azimuths,
®, was of the order of a few percent, for media that could be classified as
weakly anisotropic. The expression for the qP group velocity was also compared
with the linearized qP group velocity expressions obtained using the First Order
Ray Tracing (FORT) method presented in Pseneik and Farra (2005), with
similar accuracy, for an orthorhombic medium. For this reason it was chosen
for use in preliminary ray tracing procedures where speed, together with
reasonable accuracy, was required.
It will be convenient to first consider the case where the anisotropic
parameters and orientation of anisotropy are assumed constant within any finite
element. What is required to be determined is an analogue of Snell’s Law at the
plane interface between two adjacent elements. Within the context of
'two-point' ray tracing this topic will be treated in what follows for an
orthorhombic anisotropic medium.
THEORETICAL BACKGROUND
The linearized quasi-compressional (qP) phase (wavefront normal)
velocity, v,p(n,), in an orthorhombic medium may be written in Voigt notation,
where the density (po) normalized stiffness coefficients, A, j = Gio, Gj =
1,2,...,6) have the dimensions of velocity squared, as
vo) = Anni + A22n3 + Am + E,ning + Einin3 + E, nny , (1)
Nn = (N,,M),N;) = (sindcos¢,sindsing,cosé) , (2)
LINEARIZED GROUP VELOCITY APPROACH 89
(Backus, 1965; Daley and Krebes, 2006) where (6,) are the polar and
azimuthal angles associated with the phase or wavefront normal velocity, vp.
The E;, are the anellipsoidal terms, specifying the deviation of the slowness or
ray surface from the ellipsoidal in the 'i j' symmetry plane, (i j = 12,13,23).
Their definitions are as follows
Ey = (An + 2Ag6) 一 (An + An) ， @)
E,3 in 2(Ai 還 2Ass) ot (Ay aie A33) 2 (4)
E = 2(A2 + 2Ay) 一 (An + Ay) . 6)
These expressions may be compared to those given for an arbitrary anisotropic
orthorhombic medium and presented in Gassmann (1964) or Schoenberg and
Helbig (1996) as an indication of the how linearization simplifies the phase
velocity expression.
The 3D phase velocity propagation direction vector, n, is defined by
eq.(2). The related slowness vector in terms of these quantities is of the form
P 三 (pi,pz,p3) = tni/vup(no9],[n2/vuop(ng],[ns/vop(ng]) (6)
It was shown in Daley and Krebes (2006) that the slowness vector p, for
the qP case may be further approximated for the linearized case and written in
terms of group velocity quantities, those being the azimuth and polar ray angles
P 三 (pi,p2,p3) = {IN, Vae(Nd/Ari]s[N2Vgp(N,)/Aoo] 5IN3Vgp(N,)/Aga] } : (7)
In the above N is the orthonormal group (ray) vector defined by
N = (N,,N,N3) = (sinOcos®,sinOsin®,cos@) . (8)
It is along the ray that energy propagates from one point in a medium to
another. The linearized group velocity expression for a qP ray in an
orthorhombic medium is given by the formula
[LVap(N] = (NVAD + (N3/A22) + (N3/Ag3)
一 (EnN3N3/AyAg9) 一 (ENIN3ANAs) 一 (ExsNiN3/An As), (9)
with the E;; being previously defined in eqs. (3)-(5). Eq. (9) may be expressed
in terms of (8,®) using eq. (8).
90 DALEY & KREBES
As the slowness vector p = (p;,P2,P3) is assumed to be known, the sines
of the group angles (fi = sinO, & = sin®) will be used as an alternate
parameterization of group velocity and related quantities. In the ellipsoidal case
(E,. = EE = E,; = 0) fi and &, may be obtained analytically as
£, = (p,/ps)[(cos*&/Aj,) + (sin?&/A3,)]-“/A3,
x {1 + (p,/p3)[(cos?@/Aj,) + (sin?/A5,)] V/A}, (10)
and
&) = (A22pz/AiipD/[1 + (A2p2/Anp0 和 - (11)
The radial component of the slowness vector, p, ,is defined as
P, = (pt + ps)” . (12)
The relationship between the phase angles (0,¢) and the group angles
(8,®) in the ellipsoidal case are calculated by equating individual components
in eqs. (6) and (8) which are used to derive eqs. (10)-(12).
SNELL’S LAW ANALOGUE
Given a three dimensional slowness surface that may be rotated at any
orientation with respect to the model axis system, it would be expected that
Euler Angles be introduced. However, as one of the motivations for this work
is to introduce 3D ray tracing in an orthorhombic anisotropic medium in the
simplest possible manner, this will not be done. When more complex ray tracing
methods (comparatively) are addressed, such as FORT, this will be necessary.
As the slowness vector in terms of phase angles and velocities is not used here,
one level of theoretical and numerical complexity has been removed, allowing
for the subsequent formulation.
Two orthorhombic media in welded contact are separated by a plane
boundary with the axes of anisotropy in both the upper (1) and lower (2)
medium being aligned with some Cartesian model coordinates system within
which the plane boundary is parallel to the horizontal (x,,x,) model coordinate
plane. As previously mentioned, this is a very simplistic problem type, but once
solved may be extended to more complex geometries with minimal effort.
determine the ray angle and magnitude of the ray velocity for either the
reflected ray in the upper medium or the transmitted ray in the lower medium.
Specifically, the two horizontal components of the slowness vector, ph = (pi,p2)
LINEARIZED GROUP VELOCITY APPROACH 91
at the point of incidence at an interface, measured with respect to the model
coordinates, are required to be determined. Once this is done, the third
component of the slowness vector, p3, follows. It has been initially assumed that
the plane interface is aligned with the model coordinate axes as are the qP
slowness surfaces in both media. As a consequence of these assumptions, the
projections of the slowness surface onto the pip and pyp planes may be
realized. From eq. (7), the continuity of the horizontal components of the
slowness vector may be stated in terms of (£,,£) by the two equations
Fi(E1,8) = prAn — £1 — &)'Vep(E.£) = 0 , (13)
B,(E,,8) = pA» 一 £,£,Vap(E1.62) = 0 . (14)
These coupled nonlinear equations are required to be solved numerically.
For coupled equations of this type Newton’s Method, or some variant, would
be a reasonable choice as a solution method, the standard formulation of which
is given here for completeness,
Ea = Ee = FYE E19 (E.£)17' > (15)
where & = (go FE.) = ([F(E..£2),.Fo(E.,€)] and the Jacobian
$.(E1,£2), is defined as
aF,/dé, 9F,/dé,
$(E,E) = 1 (16)
aF,/dé, OF,/de, |,
Its inverse is given by
aF,/dE, —OF,/dé,
[SG = (1/D,) (17)
—dF,/aé, dF,/aé, |,
D, = [(0F,/0&,)(6F,/0&,) 一 (OF,/0E,)(OF,/0E,)], « (18)
The superscripts 'T' and '—1' indicate transpose and inverse, respectively.
The partial derivatives of Fi(g,,g,),G = 1,2) with respect to fi and #
analytic expressions for the partial derivatives of Vap(E 1,2) with respect to &,
and £, are given in the Appendix. The initial 'guesses' for the two components
of the quantity £) may be found using the ellipsoidal (degenerate) eqs. (10) and
(11). Once = (&,,£,) has been determined to within a specified tolerance, the
group velocity and as a result, the three components of the slowness vector in
the medium of reflection or transmission, may be determined using eqs. (7).
92 DALEY & KREBES
SNELL’S LAW ANALOGUE AT A PLANE INTERFACE WHERE THE
AXES OF ANISOTROPY ARE NOT ALIGNED WITH THE INTERFACE
For this problem, it is assumed that the slowness vector p, relative to the
model coordinates is known. To solve the problem, that is to determine the
polar ray angles (0,®) relative to the model coordinates, the angles in the
rotated system (0' = 9 十 Y, 6’ = 6+) must be determined, under the
assumption that the polar and azimuthal rotation angles (y,{ are also known.
The primed angles define the slowness vector, p’, in the rotated coordinate
system and the angle pair (y,{) specifies the relationship between (0,®) and
(0',8').
With pi = (pi,p2) known, a system of nonlinear equations in terms of :|
and &5 is required to be determined where
Ei
总
sin(O+y) = sin9' , (19)
and
sin(@+) = sin’ . (20)
The following nonlinear system is, apart from the 'primed' notation, the
same as that for the unrotated (unprimed) case and may be solved in a similar
manner.
Fi = PiAn 一 Van — &°)* = 0 , (21)
FA(E1.£2) = prAn 一 Vgr(EEdelEs 一 0 (22)
What is required to solve this system using numerical methods is initial
estimates for Ki and 3. This, as previously done, may be accomplished by
assuming the ellipsoidal case and solving for (gap_and (£%).ui, using the
sequence of eqs. (10)-(12) in the 'primed' system. As the system of nonlinear
equations above is essentially the same as the unrotated problem, the solution
for this case follows.
Once (£;,&3) has been determined, the three components of the slowness
vector in the primed system may be obtained. Given that the vector p’ is
known, p may be obtained using the inverse of the rotational transforms given
by eqs. (19) and (20).
It should be noted that the orientation of the local vertical axis in both
slowness space and ray space must be taken in the proper sense as a reflected
and transmitted ray and slowness vertical component of each individual vector
have different signs. It is often better to do away with the signs and measure the
angles as acute with respect to either the positive or negative direction of the
vertical axis and one of the horizontal axes in model coordinates.
LINEARIZED GROUP VELOCITY APPROACH 93
RAY TRACING METHOD
Assume a medium composed of N plane parallel homogeneous anisotropic
layers where both the source and receiver are located at the surface. The source
is situated at the origin of a Cartesian system aligned with the plane layered
medium, with the vertical Cartesian axis, x,, chosen to be positive downwards.
A description of unconverted qP ray propagation is given by the two coupled
nonlinear equations in terms of fi and &, as
M
Riog) = 1 -@DI' — Y (nthe jl - ED '1 - 41}
j=l
M
~ Y faye) EI- ET = 0 ， (23)
R2(Gnog2) = rt 一刻 Ef} E01 — EDI}
j=l
M
= » {njh(E)j(E)/-EDi4} = 0 , (24)
where M(M < N) is the deepest layer traversed by the ray. The integers n, are
the number of P-ray segments in the j-th layer, with the superscripts '个 ' and
')' indicating whether the ray segment(s) are upward or downward
propagating. As upward and downward propagating ray segments within a layer
do not generally make the same angle with vertical or horizontal axes they are
considered separately. The thickness of the j-th layer is taken to be h;. What
should be noted in eqs. (23)-(24) is that there are no 'primes' on the gi (i =
1,2). However, they are implied by the use of the 't' and '4' notation. The
angle ®..(Es)ures(1 — (ED) fue) “] is the angle measured in a positive (clockwise)
direction relative to the x, axis from the source to a receiver on the surface at
r = (xi,x§), with r = |r|.
Embedded in these two equations are the two coupled nonlinear equations
for the analogue of Snell’s Law discussed in an earlier section. Thus, this
'two-point' ray tracing scheme is quite computationally intensive. The method
or methods used in the solution of these numerical problems include those of
Newton among others. The technique chosen is usually based on previous
experience. It should however be capable of dealing with more complex
situations, a detail to be kept in mind when devising the ray tracing code. In the
next section examples will be presented.
Finally, the traveltime along the ray may also be written in terms of the
sines [(é;)''(i = 1,2), j=1, M) of the angles (Of', 1', j=1, M) as
94 DALEY & KREBES
M
EDGED] = Y nthy{l 一 LEDPPP} VIED} ED
j=l
M
+ Vinth/{1 — [€ PP} VIELE . 5)
j=l
COMMENTS ON THE SOLUTION AND NUMERICAL RESULTS
The problem discussed here may be approached in a number of ways,
depending on what the ray tracing code is being used for. It will be assumed
that ray tracing for the purposes of synthetic seismogram production,
Born-Kirchhoff migration or some related task is being considered.
The 3D ray tracing algorithm for qP rays in an orthorhombic anisotropic
medium is reasonably complex. However, for it to be used in a productive
manner, a relatively sophisticated model building program is an essential
addition. In 2D, a grid with all of the anisotropic parameters, together with
other required quantities, could be specified at each node of this grid. At
present, a realistic 3D grid is not an option unless top level computer hardware
is available. Few persons would have access to this current generation of 'super
computers'. A reasonably fast laptop with 2 - 4 gigabytes of memory should be
adequate if the model is specified as a 'blocky type structure', composed of
slabs or wedges or other similar types of constituents, separated by plane
interfaces. For some applications these interfaces may be smoothed, but the ray
tracing program should be written to handle jump discontinuities in the
anisotropic parameters defining each of the blocks in the model. Although these
types of model building programs exist, obtaining a freeware type is not as easy
as might be expected. Those which are in the public domain usually require
more than a minimum of effort to be of use for an existing ray tracing program.
It is highly advisable to choose and become acquainted with one of these and
then proceed to write the ray tracing code. It should not have to be mentioned
that the possibility of errors in these types of codes is high. Additionally, the
known 3D model building software requires at least one, usually two, licensed
numerical software packages, which introduces more inconveniences. If
alternatives in the public domain are used to replace the equivalent of the
licensed packages, integrating these into the model building code can produce
another source of error as well as an increased development time. This leaves
the preferable option of spending the time to write one’s own (not an
insignificant undertaking).
The ray tracing described here is computationally fast. Most of the time
required is iterating to a solution, which increases if bad initial guesses are used.
LINEARIZED GROUP VELOCITY APPROACH 95
Consequently, shooting the initial ray of a linear sequence at near vertical
incidence is done followed by considering a line of receivers, increasing in
distance at polar angle increments from the origin, at a constant azimuth. This
is applicable to both offset shooting and shooting to some reference depth.
In the first example considered here a plane layered three layer model is
used. The surface layer is chosen to be isotropic. Layer two has an
orthorhombic type structure similar to olivine, with the third layer being an
orthorhombic sandstone type medium. In layer 2 the symmetry axes are not
aligned with 'model coordinates' but rather rotated about the vertical axis by
an angle '®,,, = 30°', measured from the positive x,-axis. In the third layer,
a rotation of the ray (slowness) surface(s) of “Or = 15°' about the x,-axis is
introduced.
Two-point ray tracing which varies somewhat from that described in the
text is first considered. Rather than solve for Of and @}', (Gj = 1,2,3), the
parameters used are @}' (j = 1,2,3) and offset. This was done so that the &!*
(j = 1,2,3) could 'float' to provide an indication of how the azimuthal rotation
in the second layer affects the arrival at the surface receivers, which are a fixed
distance apart. Put another way, what is done here is to keep the azimuthal
phase angle constant for a given ray, maintain a constant receiver spacing and
solve for the polar ray angle.
Three views of this is given in Fig. (1). The scale varies along all axes
in the figure to produce plots that are easily viewed. It should be noted that
proper ray 3D plotting software could significantly enhance the presentation of
results. A similar plot for a number of azimuths is given in Fig. (2). Again, the
ray azimuthal angles, in each layer, along each ray, are not used as solution
parameters. The anisotropic parameters used are given in Tables 1 and 2. Next,
a ray tracing problem is considered in which the polar and azimuthal ray angles
in each layer and along each ray, together with the distance between the surface
located source and receivers are required to satisfy certain conditions. This is
the type of ray tracing which was discussed in this paper. Upon viewing Fig.3,
where results are presented for this case, it may be seen that the progression of
the receiver location describes a straight line from the source to the furthest
offset.
An alternate example uses the same model as in the previous case and the
'two-point' ray tracing described in the text is used. Here a number of azimuths
are chosen and the rays are computed from a source located at the surface down
to some reference depth (Fig. 4). The motivation for this is that there are
occasions in seismic processing when a 'time to depth' map is useful. If the
coverage is sufficient the computed data may be used in migration programs.
Usually, the coverage must be enhanced which may be done using any number
of methods or strategies. The pursuit of these is beyond the scope of this work.
96 DALEY & KREBES
Fig. 1. Three layer model. The first layer is isotropic, the second layer is orthorhombic with a
rotation in the x - y plane about the z-axis and the third layer is also orthorhombic with a rotation
in the x - z plane. The line shot is very close to two-point ray tracing, requiring some minor
tweaking. The analogue of Snell’s Law at an interface is a system of two coupled nonlinear equation
in two unknown angles (sines of the transmitted or reflected angles). The x- and y-axes in the three
panels have different scaling to enhance the effects of the rotations of the orthorhombic ray surfaces.
The line is shot at an angle of 70 degrees with respect to the x-axis. The modified two-point ray
tracing algorithm consists of two coupled nonlinear equations in two unknown angles per layer in
which is embedded the other nonlinear equation set for Snell’s Law. Different scaling parameters
are used on three panels to improve visibility.
LINEARIZED GROUP VELOCITY APPROACH 97
Table 1. The anisotropic parameters (density normalized stiffnesses) in the three layers in Voigt
notation. The A,, have the dimensions of velocity squared (km’/s’).
An An As Aw As Ac An An As
00 7.00 7.00 2.50 2.50 2.50 2.00 2.00 2.00
900 6.023 7.093 1.964 2.448 2.438 1.926 2.074 2.225
00 9.84 5.94 2.00 1.60 2.18 3.60 2.25 2.40
Table 2. The anellipsoidal parameters in the three layers. The E;; have the dimensions of velocity
squared (km’/s*) and are defined in terms of the Aij using eqs. (3)-(5) in the text.
Eu En Ey
0 0.0 0.0
一 2.319 一 3.053 一 0.81
—3.92 一 5.04 一 2.98
4 3
2 3
aren x (®=0)
Fig. 2. Three layer model. The model used is the same as that in Fig. 1, except that a number of
a fixed distance between adjacent receivers are again used as the dependent variables.
98 DALEY & KREBES
*x(@=0)
Fig. 3. Three layer model. The model used is the same as that in Fig. 2, except that the polar angle
azimuthal line are now in a straight line at the surface.
CONCLUSIONS
The basic formulae and solution method for tracing quasi-compressional
(qP) rays in a plane layered orthorhombic anisotropic medium where the axes
of anisotropy (polar and azimuthal) in any of the layers may be rotated with
respect to the general Cartesian coordinate system have been presented. A
linearized form of the qP ray velocity is used for this. Proceeding in this manner
removes one 'level' of nonlinear equations. Included in this discussion is a
presentation of an analogue of Snell’s Law for reflected and transmitted rays at
an interface separating two layers. In general, the 'two-point' ray tracing
problem consists of two nonlinear equations (Snell’s Law) embedded in two
nonlinear (qP) ray tracing equations. Writing the shear ray velocities as is done
for the qP group velocity does not present itself in a manner as that for the qP
case. Obtaining accurate group velocity expressions for the two possible shear
wave modes is also a fairly complex undertaking.
LINEARIZED GROUP VELOCITY APPROACH 99
= 7/2)
y(®
x(D =0)
Fig. 4. Three layer model. In this case what is wanted is time to some reference depth. A sampling
of three azimuths are shown. The anisotropic parameters describing the layers are given in Tables
1 and 2.
100 DALEY & KREBES
ACKNOWLEDGEMENTS
The first author wishes to thank the sponsors of CREWES and NSERC
(Professor G.F. Margrave, CRDPJ 379744-08) for financial support in
undertaking this work.
REFERENCES
Backus, G.E., 1965. Possible forms of seismic anisotropy of the uppermost mantle under oceans.
J. Geophys. Res., 70: 3429-3439.
Courant, R. and Hilbert, D., 1962. Methods of Mathematical Physics, Vol. II: Partial Differential
Equations. Interscience Publications, New York.
Cerveny, V., 2002. Seismic Ray Theory., Cambridge University Press, Cambridge.
Daley, P.F. and Krebes, E.S., 2006. Quasi-compressional group velocity approximations in a
weakly anisotropic orthorhombic medium. J. Seismic Explor., 14: 319-334.
Every, A.G., 1980. General closed-form expressions for acoustic waves in elastically anisotropic
solids. Phys. Rev. B, 22: 1746-1760.
Every, A.G. and Sachse, W., 1992. Sensitivity of inversion algorithms for recovering elastic
constants of anisotropic solids from longitudinal wavespeed data. Ultrasonics, 30: 43-48.
Farra, V. and Pgen¢ik, I., 2008. First-order ray computations of coupled S waves in inhomogeneous
weakly anisotropic media. Geophys. J. Internat., 173: 979-989.
Gassmann, F., 1964. Introduction to seismic travel time methods in anisotropic media. Pure Appl.
Geophys., 58: 63-112.
Jech, J. and P&enéik, I., 1989, First order perturbation method for anisotropic media. Geophys. J.
Internat., 99: 369-376.
Pseneik, I. and Gajewski, D., 1998. Polarization, phase velocity and NMO velocity of qP waves
in arbitrary weakly anisotropic media. Geophys., 63: 1754-1766.
P&encik, I. and Farra, V., 2005. First-order ray tracing for qP waves in inhomogeneous weakly
anisotropic media. Geophysics, 70: D65-D75.
Schoenberg, M. and Helbig, K., 1996. Orthorhombic media: Modeling elastic wave behavior in a
vertically fractured earth. Geophysics, 62: 1954-1974.
Song, L.-P. and Every, A.G., 2000. Approximate formulae for acoustic wave group slowness in
weakly orthorhombic media. J. Phys. D: Appl. Phys., 33: L81-L85.
Song, L.-P., Every, A.G. and Wright, C., 2001. Linearized approximations for phase velocities of
elastic waves in weakly anisotropic media. J. Phys. D: Appl. Phys., 34: 2052-2062.
LINEARIZED GROUP VELOCITY APPROACH 101

Previous article in this issue

Next article in this issue

Journal of Seismic Exploration, Electronic ISSN: 0963-0651 Print ISSN: 0963-0651, Published by AccScience Publishing

Publisher's Core Philosophy

We are committed to support the scientific community by publishing impactful research and enhancing communication among scientists. At AccScience Publishing, we are continuously looking for ways to accelerate scientific progress and to strive for transparency and open communication, making knowledge freely accessible without barrier.

9 Raffles Place, Republic Plaza 1 #06-00 Singapore 048619

+65 8182 1586

editorial@accscience.com