Calculation of Multivalued Ray Theory Travel Times in 3-D Structures.

Petr Bulant

Introduction

The thesis is composed of three author's papers written during the three-year period of his PhD studies. The research was devoted mostly to solving the problems of the ray-theory based computation of seismic wavefields propagating in complex three-dimensional laterally varying block structures. The author continued in the development of the sophisticated two-point ray tracing algorithm, suggested in his master thesis. The robust algorithm was designed to provide all the two-point rays of given wave for given source-receiver configuration, including multiple arrivals. The main attention was paid to the reliability of the algorithm, and also to its proper and efficient coding and to the numerical tests of its computer implementation. Resulting code was incorporated into already existing Fortran 77 subroutine package CRT. Considerable attention has then been paid to the development of a new method for computation of multivalued arrivals at nodes of 3-D grids in complex 3-D structures. The method is in principle equivalent to the wavefront tracing method. It is based on construction of ray tubes, composed of ray cells, through the model volume, followed by interpolation inside individual ray cells. The discussion of the results is classified according to these two topics.

1 Two-point ray tracing

Recent ray-theory based numerical methods enable rays of various seismic waves to be computed in 3-D laterally varying inhomogeneous media with a block structure. Given the initial conditions for the ray, they can be used to compute the ray trajectory, the time of propagation, and the slowness vector along any ray. They also allow the ray-theory amplitudes to be computed and synthetic seismograms to be constructed. Thus, the problem of initial-value ray tracing is, to some extents, solved. On the other hand, we are usually more interested in computation of travel times and other quantities at given receivers (the points of given coordinates). The task of computation of the ray connecting a source with the receiver, two-point ray tracing, is much more complex. There are basically two approaches to the solution of the two-point ray tracing problem: bending and shooting. The bending method is dependent on the approximate ray tracer estimating the initial rays to be bent. In a case of multiples it thus need not provide us with complete solution (i.e. with all the arrivals). Thus it is more reliable to use any of shooting methods for solving the two-point ray tracing problem. In the shooting method, it is necessary to find the boundaries between individual regions of the same history of rays of the elementary wave under consideration. The ray history is an integer-valued function, which assigns the rays to various groups according to the structural blocks and interfaces through which the ray has propagated, as well as according to the position of its endpoint (at the bottom of the model, at the sides of the model, at the surface, overcritical incidence at an interface, ...), and the caustics encountered. The procedure of identification of individual regions must be done very carefully in order to avoid omitting any of the regions, which may contain two-point rays. Following these ideas, the two-point ray tracing program being described is based on the procedures to demarcate the boundaries between individual regions in the ray-parameter domain, and on the proper triangularization of these regions. The boundaries are demarcated by pairs of so called boundary rays, each of them being situated at the opposite side of the boundary. The mutual distance of these two rays is kept under a priori prescribed (very small) limit, because all the two-point rays situated inside the demarcation belts will be lost. The distance between individual pairs of boundary rays along the boundary is also controlled and kept under prescribed limit, depending on the curvature of the boundary being demarcated. This very accurate demarcation of the boundaries should ensure that any of the regions, which is greater than prescribed quantities (i.e. distances of boundary rays and of pairs of rays), will not be omitted. After demarcation of the boundaries, individual regions are triangularized by approximately equilateral triangles with the maximum length of the sides given. Two-point rays are then iteratively identified within individual triangles using paraxial ray approximation. The described two-point ray tracing method was suggested in [1], [2] and [4] and described in detail in [5], the importance of careful demarcation of the boundaries between individual regions was discussed in [3], differences between various triangularizations of individual regions were shown in [6] and [11].

2 Calculation of multivalued traveltimes at gridpoints of 3-D grids

As was discussed above, the described two-point ray tracing algorithm is based on demarcation of the boundaries between individual ray-parameter regions of the same history of the elementary wave under consideration, and on the triangularization of these regions. Once the individual regions of ray take-off parameters are triangularized, we can use the triangles to construct ray tubes through the model volume. Each ray tube is then bounded by three rays, the vertices of the corresponding triangle in the ray-parameter domain. Each ray tube is composed of ray cells. Each ray cell corresponds to the space between two consecutive wavefronts and/or structural interfaces. Travel times and other quantities may be then interpolated inside individual ray cells, like in the wavefront tracing method. The described method is a useful extension of two-point ray tracing programs based on triangularization. It provides interpolated values of ray-theory travel times and other quantities at arbitrary 3-D grids located in the seismic model, including multiple arrivals. Inside blocks a ray cell is determined by 6 points located on the corresponding three rays, and it is determined by 4, 5 or 6 points in the vicinity of structural interfaces. The described interpolation method then reduces into proper organization of the loop over the ray cells and in proper decision, which gridpoints of the target grid are located in the ray cell being examined. The interpolation formula for individual gridpoint is connected with the system of equations for decision of whether the gridpoint is located in the selected ray cell. Thus the interpolation of the travel time and other quantities is very straightforward once the gridpoint is identified as being located in the ray cell. Described method for computation of multivalued arrivals was suggested in [6] and [11], and the first numerical experiments were performed ([9],[10]). Recently the method is being coded in Fortran 77 as a post-processing program to the CRT program [8].

List of author's publications

[1] Bulant, P., 1994, Two-point Ray Tracing in 3-D. Diploma Thesis, Dep. Geophys., Charles Univ., Prague.

[2] Bulant, P., 1995, Two-point ray tracing in 3-D. In: Seismic Waves in Complex 3-D Structures, Report 3, pp. 37-64, Dep. Geophys., Charles Univ., Prague.

[3] Bulant, P., 1995, Numerical examples of two-point ray tracing in 3-D. In: Seismic Waves in Complex 3-D Structures, Report 3, pp. 65-92, Dep. Geophys., Charles Univ., Prague.

[4] Bulant, P., & Klimes, L., 1996, Two-point ray tracing in 3D heterogeneous structures. Extended Abstracts, C045, EAGE 58th Conference and Technical Exhibition, Amsterdam.

[5] Bulant, P., 1996, Two-point ray tracing in 3-D. Pure and Applied Geophysics 148, 421-446.

[6] Bulant, P., 1996, Two-point ray tracing and controlled initial-value ray tracing in 3-D heterogeneous block structures. In: Seismic Waves in Complex 3-D Structures, Report 4, pp. 61-76, Dep. Geophys., Charles Univ., Prague.

[7] Bulant, P., & Klimes, L., 1996, Examples of seismic models. Part 2. In: Seismic Waves in Complex 3-D Structures, Report 4, pp. 39-52, Dep. Geophys., Charles Univ., Prague.

[8] Bulant, P., Cerveny, V., Klimes, L. & Psencik, I., 1996, MODEL version 5.00 & CRT version 5.00 (2 floppy disks). In: Seismic Waves in Complex 3-D Structures, Report 4, p. 91, Dep. Geophys., Charles Univ., Prague.

[9] Bulant, P., 1997, Controlled initial-value ray tracing and two-point ray tracing in 3D. Extended Abstracts, P025, EAGE 59th Conference and Technical Exhibition, Geneva.

[10] Bulant, P., 1997, Two-point ray tracing and controlled initial-value ray tracing in 3-D heterogeneous block structures. Expanded Abstracts, 67th SEG International Exposition and Annual Meeting, Dallas, in press.

[11] Bulant, P., 1997, Two-point ray tracing and Controlled initial-value ray tracing in 3-D heterogeneous block structures. Geophysics, submitted.

[12] Bulant, P., 1997, Recent development of the shooting algorithm. In: Seismic Waves in Complex 3-D Structures, Report 6, submitted.


PhD Thesis, Department of Geophysics, Charles University, Prague, 1997.
SW3D - main page of consortium Seismic Waves in Complex 3-D Structures .