ARTICLE

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

SEUNG-GOO KANG1 WANSOO HA2 CHANGSOO SHIN3 JONG KUK HONG1 YOUNG KEUN JIN1
Show Less
1 Division of Polar Earth-System Sciences, Korea Polar Research Institute, KIOST, Incheon 406-840, South Korea.,
2 Department of Energy Resources Engineering, Pukyong National University, Busan 608-737, South Korea. wansooha@gmail.com,
3 Department of Energy Systems Engineering, Seoul National University, Seoul 151-742, South Korea.,
JSE 2016, 25(1), 57–101;
Submitted: 13 July 2015 | Accepted: 20 November 2015 | Published: 1 February 2016
© 2016 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/ )
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
  1. Bae, H.S., Shin, C., Cha, Y.H., Choi, Y. and Min, D., 2010. 2D acoustic-elastic coupledwaveform inversion in the LaPlace domain. Geophys. Prosp., 58: 997-1010.
  2. Baysal, E., Kosloff, D.D. and Sherwood, J.W.C., 1983. Reverse time migration. Geophysics, 48:1514-1524.
  3. Choi, Y., Min, D.-J. and Shin, C., 2008. Two-dimensional waveform in- version of multicomponentdata in acoustic-elastic coupled media. Geophysics, 56: 863-881.
  4. Chung, W., Pyun, S., Bae, H.S., Shin, C. and Marfurt, K.J., 2012. Implementation of elasticreverse-time migration using wavefield separation in the frequency domain. Geophys. J. Int.,189: 1611-1625.
  5. Cohen, G.C., 2002. High-Order Numerical Methods for Transient Wave Equations. SpringerVerlag, Berlin.
  6. Kang, S.-G., Bae H.S. and Shin, C., 2012. Laplace-fourier-domain waveform inversion forfluid-solid media. Pure Appl. Geophys., 169: 2165-2179.REVERSE TIME MIGRATION 85
  7. Kim, M.H., Choi, Y., Cha, Y.H. and Shin, C., 2009. 2-D frequency-domain waveform inversionof coupled acoustic-elastic media with an irregular interface. Pure Appl. Geophys., 166:1967-1985.
  8. Kim, Y., Min, D. and Shin, C., 2011. Frequency-domain reverse-time migration with sourceestimation. Geophysics, 76: -S49.
  9. Komatitsch, D., Barnes, C. and Tromp, J., 2000. Wave propagation near a fluid-solid interface: Aspectral-element approach. Geophysics, 65: 623-631.
  10. Lailly, P., 1983. The seismic inverse problem as a sequence of before stack migrations. Expanded
  11. Abstr., Conf. Inverse Scattering: Theory and Application Mathematics: 206-220.
  12. Lee, H., Lim, S., Min, D., Kwon, B. and Park, M., 2009. 2D time-domain acoustic-elastic coupledmodeling: a cell-based finite-difference method. Geosci. J., 13: 407-414.
  13. Lee, J., Kim, Y. and Shin, C., 2012. Frequency-domain reverse time migration using the L,-norm.J. Seismic Explor., 21: 281-300.
  14. Loewenthal, D. and Mufti, I.R., 1983. Reversed time migration in spatial frequency domain.Geophysics, 48: 627-635.
  15. Martin, G.S., Marfurt, K.J. and Larsen, S. 2002. Marmousi-2: an updated model for theinvestigation of AVO in structurally complex areas. Expanded Abstr., 72nd Ann. Internat.SEG Mtg., Salt Lake City: 1979-1982.
  16. McMechan, G.A., 1983. Migration by extrapolation of time-dependent boundary values. Geophys.Prosp., 31: 413-420.
  17. Plessix, R.-E. and Mulder, W.A., 2004. Frequency-domain finite-difference amplitude preservingmigration. Geophys. J. Internat., 157: 975-987.
  18. Pratt, R.G., Shin, C. and Hicks, G.J., 1998. Gauss-Newton and full Newton method in frequencydomain seismic waveform inversion. Geophys. J. Internat., 133: 341-362.
  19. Shin, C., Jang, S. and Min, D., 2001. Improved amplitude preservation for prestack depth migrationby inverse scattering theory. Geophys. Prosp., 49: 592-606.
  20. Shin, C., Min, D.J., Yang, D. and Lee, S.K., 2003. Evaluation of poststack migration in terms ofvirtual source and partial derivative wavefields. J. Seismic Explor., 12: 17-37.
  21. Symons, N.P., Aldridge, D.F. and Haney, M.M., 2006. 3D acoustic and elastic modelling with
  22. Marmousi-2. Expanded Abstr., 76th Ann. Internat. SEG Mtg., New Orleans: 2171-2175.
  23. Tarantola, A., 1984. Inversion of seismic reflection data in the acoustic approximation. Geophysics,49: 1259-1266.
  24. Virieux, J., 1986. P-SV wave propagation in heterogeneous media: velocity-stress finite-differencemethod. Geophysics, 51: 889-901.
  25. Whitmore, N.D., 1983. Iterative depth migration by backward time propagation. Expanded Abstr.,53rd Ann. Internat. SEG Mtg., Las Vegas: 382-385.
  26. Xu, K., Zhou, B. and McMechan, G.A., 2010. Implementation of prestack reverse time migrationusing frequency-domain extrapolation. Geophysics, 75: S61-S72.
  27. Zienkiewicz, O.C., Taylor, R.L. and Zhu, J.Z., 2005. The Finite Element Method: Its Basis and
  28. Fundamentals. Elsevier, Butterworth-Heinemann, New York.
  29. JOURNAL OF SEISMIC EXPLORATION 25, 87-101 (2016) 87A LINEARIZED GROUP VELOCITY APPROACH FORTWO-POINT qP RAY TRACING IN A LAYEREDORTHORHOMBIC MEDIUMP.F. DALEY and E.S. KREBES
  30. Department of Geoscience, University of Calgary, 2500 University Drive, Northwest Calgary,
  31. Alberta, Canada T2N 1N4. pfdaley@gmail.com, krebes@ucalgary.ca(Received August 20, 2015; revised version accepted November 27, 2015)ABSTRACT
  32. Daley, P.F. and Krebes, E.S., 2016. A linearized group velocity approach for two-point qP raytracing in a layered orthorhombic medium. Journal of Seismic Exploration, 25: 87-101.
  33. Using a linearized approximation for the quasi-compressional phase velocity, vip in anorthorhombic anisotropic medium, which is a subset of the related quasi-compressional (qP) wavepropagation in a general 21 parameter anisotropic medium, a linearized compressional group velocitymay be derived as a function of group angles only. In addition, linearized analytic expressions forthe components of the slowness vector in terms of group velocities and angles are also obtained.
  34. These expressions are used to define two nonlinear equations which are a generalization of Snell’s
  35. Law. The solutions of these are used to determine the propagation directions of the reflected andtransmitted rays due to an incident ray at an interface between two orthorhombic media. The axesof anisotropy, in both media are, in general, not aligned with the interface separating them.
  36. Computer code has been written to consider ray tracing in media defined by a type of large scale3D finite element blocking (blocky structures). However, a plane parallel layered medium will beused in preliminary investigations. Additionally, in each of these elements (layers) the anisotropicparameters are initially assumed to be constant, although provisions for at least minimal spatialvariations of the anisotropic parameters have been considered.
  37. KEY WORDS: orthorhombic media, two-point ray tracing, group velocity,layered anisotropic media.INTRODUCTION
  38. In the geophysical literature related to wave propagation in anisotropicmedia, specifically quasi-compressional qP waves in a medium displayingorthorhombic symmetry, a linearized approach to derive an approximate phasevelocity expression for quasi-compressional qP wave propagation has beenpresented in Backus (1965) and more recently by Pseneik and Farra (2005).0963-065 1/16/.00 © 2016 Geophysical Press Ltd.88 DALEY & KREBES
  39. These and other methods, such a perturbation theory, have been employedto formulate expressions for quantities related to the propagation of qP wavesas well as for the two shear wave modes, qS, and qS,, in a general 21 parameteranisotropic medium (Every, 1980; Every and Sachse, 1992; Jech and PSenéik,1989; Pseneik and Gajewski, 1998; Song and Every, 2000; Song et al., 2001;
  40. Pseneik and Farra, 2008). Once linearized phase velocity approximations havebeen obtained, eikonal equations with comparable accuracy may be written.
  41. From these, using the method of characteristics, (Courant and Hilbert, 1962;
  42. Cerveny, 2002) the formulae for the vector components of group velocity maybe obtained.
  43. Employing an extension of this method of approximation, Daley and
  44. Krebes (2006) obtained expressions for the qP group velocities in anorthorhombic anisotropic medium, expressed in terms of group angles. It wasfurther shown that for a weakly anisotropic medium the average deviation of thisapproximation from the exact expression for the qP group velocity, in anorthorhombic 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 asweakly anisotropic. The expression for the qP group velocity was also comparedwith the linearized qP group velocity expressions obtained using the First Order
  45. Ray Tracing (FORT) method presented in Pseneik and Farra (2005), withsimilar accuracy, for an orthorhombic medium. For this reason it was chosenfor use in preliminary ray tracing procedures where speed, together withreasonable accuracy, was required.
  46. It will be convenient to first consider the case where the anisotropicparameters and orientation of anisotropy are assumed constant within any finiteelement. What is required to be determined is an analogue of Snell’s Law at theplane interface between two adjacent elements. Within the context of'two-point' ray tracing this topic will be treated in what follows for anorthorhombic anisotropic medium.THEORETICAL BACKGROUND
  47. 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, asvo) = Anni + A22n3 + Am + E,ning + Einin3 + E, nny , (1)
  48. Nn = (N,,M),N;) = (sindcos¢,sindsing,cosé) , (2)LINEARIZED GROUP VELOCITY APPROACH 89(Backus, 1965; Daley and Krebes, 2006) where (6,) are the polar andazimuthal angles associated with the phase or wavefront normal velocity, vp.
  49. The E;, are the anellipsoidal terms, specifying the deviation of the slowness orray surface from the ellipsoidal in the 'i j' symmetry plane, (i j = 12,13,23).Their definitions are as followsEy = (An + 2Ag6) 一 (An + An) , @)E,3 in 2(Ai 還 2Ass) ot (Ay aie A33) 2 (4)E = 2(A2 + 2Ay) 一 (An + Ay) . 6)
  50. These expressions may be compared to those given for an arbitrary anisotropicorthorhombic medium and presented in Gassmann (1964) or Schoenberg and
  51. Helbig (1996) as an indication of the how linearization simplifies the phasevelocity expression.
  52. The 3D phase velocity propagation direction vector, n, is defined byeq.(2). The related slowness vector in terms of these quantities is of the form
  53. P 三 (pi,pz,p3) = tni/vup(no9],[n2/vuop(ng],[ns/vop(ng]) (6)
  54. It was shown in Daley and Krebes (2006) that the slowness vector p, forthe qP case may be further approximated for the linearized case and written interms of group velocity quantities, those being the azimuth and polar ray angles# and ©, and the group/ray velocity, V,»(0,®), as
  55. P 三 (pi,p2,p3) = {IN, Vae(Nd/Ari]s[N2Vgp(N,)/Aoo] 5IN3Vgp(N,)/Aga] } : (7)
  56. In the above N is the orthonormal group (ray) vector defined byN = (N,,N,N3) = (sinOcos®,sinOsin®,cos@) . (8)
  57. It is along the ray that energy propagates from one point in a medium toanother. The linearized group velocity expression for a qP ray in anorthorhombic 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 expressedin terms of (8,®) using eq. (8).90 DALEY & KREBES
  58. As the slowness vector p = (p;,P2,P3) is assumed to be known, the sinesof the group angles (fi = sinO, & = sin®) will be used as an alternateparameterization 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)
  59. The radial component of the slowness vector, p, ,is defined asP, = (pt + ps)” . (12)
  60. The relationship between the phase angles (0,¢) and the group angles(8,®) in the ellipsoidal case are calculated by equating individual componentsin eqs. (6) and (8) which are used to derive eqs. (10)-(12).SNELL’S LAW ANALOGUE
  61. Given a three dimensional slowness surface that may be rotated at anyorientation with respect to the model axis system, it would be expected that
  62. Euler Angles be introduced. However, as one of the motivations for this workis to introduce 3D ray tracing in an orthorhombic anisotropic medium in thesimplest possible manner, this will not be done. When more complex ray tracingmethods (comparatively) are addressed, such as FORT, this will be necessary.
  63. 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, allowingfor the subsequent formulation.
  64. Two orthorhombic media in welded contact are separated by a planeboundary with the axes of anisotropy in both the upper (1) and lower (2)medium being aligned with some Cartesian model coordinates system withinwhich the plane boundary is parallel to the horizontal (x,,x,) model coordinateplane. As previously mentioned, this is a very simplistic problem type, but oncesolved may be extended to more complex geometries with minimal effort.
  65. Given that the incident ray angles (©,,6;) are known it is required todetermine the ray angle and magnitude of the ray velocity for either thereflected ray in the upper medium or the transmitted ray in the lower medium.
  66. Specifically, the two horizontal components of the slowness vector, ph = (pi,p2)LINEARIZED GROUP VELOCITY APPROACH 91at the point of incidence at an interface, measured with respect to the modelcoordinates, are required to be determined. Once this is done, the thirdcomponent of the slowness vector, p3, follows. It has been initially assumed thatthe plane interface is aligned with the model coordinate axes as are the qPslowness surfaces in both media. As a consequence of these assumptions, theprojections of the slowness surface onto the pip and pyp planes may berealized. From eq. (7), the continuity of the horizontal components of theslowness vector may be stated in terms of (£,,£) by the two equations
  67. Fi(E1,8) = prAn — £1 — &)'Vep(E.£) = 0 , (13)B,(E,,8) = pA» 一 £,£,Vap(E1.62) = 0 . (14)
  68. These coupled nonlinear equations are required to be solved numerically.
  69. For coupled equations of this type Newton’s Method, or some variant, wouldbe a reasonable choice as a solution method, the standard formulation of whichis 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 asaF,/dé, 9F,/dé,$(E,E) = 1 (16)aF,/dé, OF,/de, |,Its inverse is given byaF,/dE, —OF,/dé,[SG = (1/D,) (17)—dF,/aé, dF,/aé, |,
  70. D, = [(0F,/0&,)(6F,/0&,) 一 (OF,/0E,)(OF,/0E,)], « (18)
  71. The superscripts 'T' and '—1' indicate transpose and inverse, respectively.
  72. 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 componentsof the quantity £) may be found using the ellipsoidal (degenerate) eqs. (10) and(11). Once = (&,,£,) has been determined to within a specified tolerance, thegroup velocity and as a result, the three components of the slowness vector inthe medium of reflection or transmission, may be determined using eqs. (7).92 DALEY & KREBES
  73. SNELL’S LAW ANALOGUE AT A PLANE INTERFACE WHERE THE
  74. AXES OF ANISOTROPY ARE NOT ALIGNED WITH THE INTERFACE
  75. For this problem, it is assumed that the slowness vector p, relative to themodel coordinates is known. To solve the problem, that is to determine thepolar ray angles (0,®) relative to the model coordinates, the angles in therotated system (0' = 9 十 Y, 6’ = 6+) must be determined, under theassumption that the polar and azimuthal rotation angles (y,{ are also known.
  76. The primed angles define the slowness vector, p’, in the rotated coordinatesystem and the angle pair (y,{) specifies the relationship between (0,®) and(0',8').
  77. With pi = (pi,p2) known, a system of nonlinear equations in terms of :|and &5 is required to be determined whereEi总sin(O+y) = sin9' , (19)andsin(@+) = sin’ . (20)
  78. The following nonlinear system is, apart from the 'primed' notation, thesame as that for the unrotated (unprimed) case and may be solved in a similarmanner.Fi = PiAn 一 Van — &°)* = 0 , (21)FA(E1.£2) = prAn 一 Vgr(EEdelEs 一 0 (22)
  79. What is required to solve this system using numerical methods is initialestimates for Ki and 3. This, as previously done, may be accomplished byassuming the ellipsoidal case and solving for (gap_and (£%).ui, using thesequence of eqs. (10)-(12) in the 'primed' system. As the system of nonlinearequations above is essentially the same as the unrotated problem, the solutionfor this case follows.
  80. Once (£;,&3) has been determined, the three components of the slownessvector in the primed system may be obtained. Given that the vector p’ isknown, p may be obtained using the inverse of the rotational transforms givenby eqs. (19) and (20).
  81. It should be noted that the orientation of the local vertical axis in bothslowness space and ray space must be taken in the proper sense as a reflectedand transmitted ray and slowness vertical component of each individual vectorhave different signs. It is often better to do away with the signs and measure theangles as acute with respect to either the positive or negative direction of thevertical axis and one of the horizontal axes in model coordinates.LINEARIZED GROUP VELOCITY APPROACH 93RAY TRACING METHOD
  82. Assume a medium composed of N plane parallel homogeneous anisotropiclayers where both the source and receiver are located at the surface. The sourceis situated at the origin of a Cartesian system aligned with the plane layeredmedium, with the vertical Cartesian axis, x,, chosen to be positive downwards.
  83. A description of unconverted qP ray propagation is given by the two couplednonlinear equations in terms of fi and &, asMRiog) = 1 -@DI' — Y (nthe jl - ED '1 - 41}j=lM~ Y faye) EI- ET = 0 , (23)R2(Gnog2) = rt 一刻 Ef} E01 — EDI}j=lM= » {njh(E)j(E)/-EDi4} = 0 , (24)where M(M < N) is the deepest layer traversed by the ray. The integers n, arethe number of P-ray segments in the j-th layer, with the superscripts '个 ' and')' indicating whether the ray segment(s) are upward or downwardpropagating. As upward and downward propagating ray segments within a layerdo not generally make the same angle with vertical or horizontal axes they areconsidered separately. The thickness of the j-th layer is taken to be h;. Whatshould 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. Theangle ®..(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 atr = (xi,x§), with r = |r|.
  84. Embedded in these two equations are the two coupled nonlinear equationsfor the analogue of Snell’s Law discussed in an earlier section. Thus, this'two-point' ray tracing scheme is quite computationally intensive. The methodor methods used in the solution of these numerical problems include those of
  85. Newton among others. The technique chosen is usually based on previousexperience. It should however be capable of dealing with more complexsituations, a detail to be kept in mind when devising the ray tracing code. In thenext section examples will be presented.
  86. Finally, the traveltime along the ray may also be written in terms of thesines [(é;)''(i = 1,2), j=1, M) of the angles (Of', 1', j=1, M) as94 DALEY & KREBESMEDGED] = Y nthy{l 一 LEDPPP} VIED} EDj=lM+ Vinth/{1 — [€ PP} VIELE . 5)j=lCOMMENTS ON THE SOLUTION AND NUMERICAL RESULTS
  87. 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 assumedthat ray tracing for the purposes of synthetic seismogram production,
  88. Born-Kirchhoff migration or some related task is being considered.
  89. The 3D ray tracing algorithm for qP rays in an orthorhombic anisotropicmedium is reasonably complex. However, for it to be used in a productivemanner, a relatively sophisticated model building program is an essentialaddition. In 2D, a grid with all of the anisotropic parameters, together withother required quantities, could be specified at each node of this grid. Atpresent, a realistic 3D grid is not an option unless top level computer hardwareis available. Few persons would have access to this current generation of 'supercomputers'. A reasonably fast laptop with 2 - 4 gigabytes of memory should beadequate if the model is specified as a 'blocky type structure', composed ofslabs or wedges or other similar types of constituents, separated by planeinterfaces. For some applications these interfaces may be smoothed, but the raytracing program should be written to handle jump discontinuities in theanisotropic parameters defining each of the blocks in the model. Although thesetypes of model building programs exist, obtaining a freeware type is not as easyas might be expected. Those which are in the public domain usually requiremore than a minimum of effort to be of use for an existing ray tracing program.
  90. It is highly advisable to choose and become acquainted with one of these andthen proceed to write the ray tracing code. It should not have to be mentionedthat the possibility of errors in these types of codes is high. Additionally, theknown 3D model building software requires at least one, usually two, licensednumerical software packages, which introduces more inconveniences. Ifalternatives in the public domain are used to replace the equivalent of thelicensed packages, integrating these into the model building code can produceanother source of error as well as an increased development time. This leavesthe preferable option of spending the time to write one’s own (not aninsignificant undertaking).
  91. The ray tracing described here is computationally fast. Most of the timerequired is iterating to a solution, which increases if bad initial guesses are used.LINEARIZED GROUP VELOCITY APPROACH 95
  92. Consequently, shooting the initial ray of a linear sequence at near verticalincidence is done followed by considering a line of receivers, increasing indistance at polar angle increments from the origin, at a constant azimuth. Thisis applicable to both offset shooting and shooting to some reference depth.
  93. In the first example considered here a plane layered three layer model isused. The surface layer is chosen to be isotropic. Layer two has anorthorhombic type structure similar to olivine, with the third layer being anorthorhombic sandstone type medium. In layer 2 the symmetry axes are notaligned with 'model coordinates' but rather rotated about the vertical axis byan 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 isintroduced.
  94. Two-point ray tracing which varies somewhat from that described in thetext is first considered. Rather than solve for Of and @}', (Gj = 1,2,3), theparameters 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 rotationin the second layer affects the arrival at the surface receivers, which are a fixeddistance apart. Put another way, what is done here is to keep the azimuthalphase angle constant for a given ray, maintain a constant receiver spacing andsolve for the polar ray angle.
  95. Three views of this is given in Fig. (1). The scale varies along all axesin the figure to produce plots that are easily viewed. It should be noted thatproper ray 3D plotting software could significantly enhance the presentation ofresults. A similar plot for a number of azimuths is given in Fig. (2). Again, theray azimuthal angles, in each layer, along each ray, are not used as solutionparameters. 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 anglesin each layer and along each ray, together with the distance between the surfacelocated source and receivers are required to satisfy certain conditions. This isthe 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 ofthe receiver location describes a straight line from the source to the furthestoffset.
  96. 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 azimuthsare chosen and the rays are computed from a source located at the surface downto some reference depth (Fig. 4). The motivation for this is that there areoccasions in seismic processing when a 'time to depth' map is useful. If thecoverage is sufficient the computed data may be used in migration programs.
  97. Usually, the coverage must be enhanced which may be done using any numberof methods or strategies. The pursuit of these is beyond the scope of this work.96 DALEY & KREBES
  98. Fig. 1. Three layer model. The first layer is isotropic, the second layer is orthorhombic with arotation in the x - y plane about the z-axis and the third layer is also orthorhombic with a rotationin the x - z plane. The line shot is very close to two-point ray tracing, requiring some minortweaking. The analogue of Snell’s Law at an interface is a system of two coupled nonlinear equationin two unknown angles (sines of the transmitted or reflected angles). The x- and y-axes in the threepanels have different scaling to enhance the effects of the rotations of the orthorhombic ray surfaces.
  99. The line is shot at an angle of 70 degrees with respect to the x-axis. The modified two-point raytracing algorithm consists of two coupled nonlinear equations in two unknown angles per layer inwhich is embedded the other nonlinear equation set for Snell’s Law. Different scaling parametersare used on three panels to improve visibility.LINEARIZED GROUP VELOCITY APPROACH 97
  100. Table 1. The anisotropic parameters (density normalized stiffnesses) in the three layers in Voigtnotation. The A,, have the dimensions of velocity squared (km’/s’).An An As Aw As Ac An An As00 7.00 7.00 2.50 2.50 2.50 2.00 2.00 2.00900 6.023 7.093 1.964 2.448 2.438 1.926 2.074 2.22500 9.84 5.94 2.00 1.60 2.18 3.60 2.25 2.40
  101. Table 2. The anellipsoidal parameters in the three layers. The E;; have the dimensions of velocitysquared (km’/s*) and are defined in terms of the Aij using eqs. (3)-(5) in the text.Eu En Ey0 0.0 0.0一 2.319 一 3.053 一 0.81—3.92 一 5.04 一 2.984 32 3aren x (®=0)
  102. Fig. 2. Three layer model. The model used is the same as that in Fig. 1, except that a number ofquasi-azimuthal lines are displayed. The modified ray tracing procedure using the polar angle © anda fixed distance between adjacent receivers are again used as the dependent variables.98 DALEY & KREBES*x(@=0)
  103. Fig. 3. Three layer model. The model used is the same as that in Fig. 2, except that the polar angle© and azimuthal angle ® are used as the dependent variables. The surface receivers for someazimuthal line are now in a straight line at the surface.CONCLUSIONS
  104. The basic formulae and solution method for tracing quasi-compressional(qP) rays in a plane layered orthorhombic anisotropic medium where the axesof anisotropy (polar and azimuthal) in any of the layers may be rotated withrespect to the general Cartesian coordinate system have been presented. Alinearized form of the qP ray velocity is used for this. Proceeding in this mannerremoves one 'level' of nonlinear equations. Included in this discussion is apresentation of an analogue of Snell’s Law for reflected and transmitted rays atan interface separating two layers. In general, the 'two-point' ray tracingproblem consists of two nonlinear equations (Snell’s Law) embedded in twononlinear (qP) ray tracing equations. Writing the shear ray velocities as is donefor the qP group velocity does not present itself in a manner as that for the qPcase. Obtaining accurate group velocity expressions for the two possible shearwave modes is also a fairly complex undertaking.LINEARIZED GROUP VELOCITY APPROACH 99= 7/2)y(®x(D =0)
  105. Fig. 4. Three layer model. In this case what is wanted is time to some reference depth. A samplingof three azimuths are shown. The anisotropic parameters describing the layers are given in Tables1 and 2.100 DALEY & KREBESACKNOWLEDGEMENTS
  106. The first author wishes to thank the sponsors of CREWES and NSERC(Professor G.F. Margrave, CRDPJ 379744-08) for financial support inundertaking this work.REFERENCES
  107. Backus, G.E., 1965. Possible forms of seismic anisotropy of the uppermost mantle under oceans.J. Geophys. Res., 70: 3429-3439.
  108. Courant, R. and Hilbert, D., 1962. Methods of Mathematical Physics, Vol. II: Partial DifferentialEquations. Interscience Publications, New York.
  109. Cerveny, V., 2002. Seismic Ray Theory., Cambridge University Press, Cambridge.
  110. Daley, P.F. and Krebes, E.S., 2006. Quasi-compressional group velocity approximations in aweakly anisotropic orthorhombic medium. J. Seismic Explor., 14: 319-334.
  111. Every, A.G., 1980. General closed-form expressions for acoustic waves in elastically anisotropicsolids. Phys. Rev. B, 22: 1746-1760.
  112. Every, A.G. and Sachse, W., 1992. Sensitivity of inversion algorithms for recovering elasticconstants of anisotropic solids from longitudinal wavespeed data. Ultrasonics, 30: 43-48.
  113. Farra, V. and Pgen¢ik, I., 2008. First-order ray computations of coupled S waves in inhomogeneousweakly anisotropic media. Geophys. J. Internat., 173: 979-989.
  114. Gassmann, F., 1964. Introduction to seismic travel time methods in anisotropic media. Pure Appl.Geophys., 58: 63-112.
  115. Jech, J. and P&enéik, I., 1989, First order perturbation method for anisotropic media. Geophys. J.Internat., 99: 369-376.
  116. Pseneik, I. and Gajewski, D., 1998. Polarization, phase velocity and NMO velocity of qP wavesin arbitrary weakly anisotropic media. Geophys., 63: 1754-1766.
  117. P&encik, I. and Farra, V., 2005. First-order ray tracing for qP waves in inhomogeneous weaklyanisotropic media. Geophysics, 70: D65-D75.
  118. Schoenberg, M. and Helbig, K., 1996. Orthorhombic media: Modeling elastic wave behavior in avertically fractured earth. Geophysics, 62: 1954-1974.
  119. Song, L.-P. and Every, A.G., 2000. Approximate formulae for acoustic wave group slowness inweakly orthorhombic media. J. Phys. D: Appl. Phys., 33: L81-L85.
  120. Song, L.-P., Every, A.G. and Wright, C., 2001. Linearized approximations for phase velocities ofelastic waves in weakly anisotropic media. J. Phys. D: Appl. Phys., 34: 2052-2062.LINEARIZED GROUP VELOCITY APPROACH 101
Share
Back to top
Journal of Seismic Exploration, Electronic ISSN: 0963-0651 Print ISSN: 0963-0651, Published by AccScience Publishing