Dynamic Ray Tracing System Research Articles

Gaussian beam prestack depth migration in heterogeneous transversely isotropic media

The transversely isotropic (TI) media approximation is commonly applied to assist in the processing of seismic data acquired in sedimentary environments. Based on anisotropic kinematic and dynamic ray tracing systems, a P-wave Gaussian beam prestack depth migration (GB-PSDM) method for TI media is introduced in this paper. The imaging principle of anisotropic GB-PSDM and the corresponding migration parameters are presented on the basis of the GB-PSDM method in isotropic media. Tests of synthetic and field seismic data show that the method is an accurate and efficient anisotropic prestack depth migration method in TI media.

Transformation relations for second-order derivatives of travel time in anisotropic media

In the computation of paraxial travel times and Gaussian beams, the basic role is played by the second-order derivatives of the travel-time field at the reference ray. These derivatives can be determined by dynamic ray tracing (DRT) along the ray. Two basic DRT systems have been broadly used in applications: the DRT system in Cartesian coordinates and the DRT system in ray-centred coordinates. In this paper, the transformation relations between the second-order derivatives of the travel-time field in Cartesian and ray-centred coordinates are derived. These transformation relations can be used both in isotropic and anisotropic media, including computations of complex-valued travel times necessary for the evaluation of Gaussian beams.

Gaussian beams in inhomogeneous anisotropic layered structures

Summary Gaussian beams concentrated close to rays of high-frequency seismic body waves propagating in an inhomogeneous anisotropic layered structure are studied. The amplitude proflles of the Gaussian beam along the plane perpendicular to the ray and along the plane perpendicular to the slowness vector are Gaussian. The Gaussian proflle is controlled by the 2 £ 2 complex-valued matrix M of the second derivatives of the travel-time fleld at any point of the ray. The matrix M can be simply determined at any point of the ray if the ray-propagator matrix along the ray is known and if the value of M is specifled at a selected point of the ray. The ray-propagator matrix can be determined by dynamic ray tracing along the ray. In inhomogeneous anisotropic medium, the dynamic ray tracing can be performed alternatively in several coordinate systems: in global Cartesian coordinates, in ray-centred coordinates and in wavefront orthonormal coordinates. In addition, also simplifled dynamic ray tracing in global Cartesian coordinates can be used, which reduces the number of equations of the dynamic ray tracing system. The derived expressions for the Gaussian beams are applicable to general 3-D inhomogeneous layered structures of arbitrary anisotropy (specifled by upto 21 independent position-dependent elastic moduli). Possible simpliflcation of the procedure are outlined.

A note on dynamic ray tracing in ray-centered coordinates in anisotropic inhomogeneous media

Dynamic ray tracing plays an important role in paraxial ray methods. In this paper, dynamic ray tracing systems for inhomogeneous anisotropic media, consisting of four linear ordinary differential equations of the first order along the reference ray, are studied. The main attention is devoted to systems expressed in a particularly simple choice of ray-centered coordinates, here referred to as the standard ray-centered coordinates, and in wavefront orthonormal coordinates. These two systems, known from the literature, were derived independently and were given in different forms. In this paper it is proved that both systems are fully equivalent. Consequently, the dynamic ray tracing system, consisting of four equations in wavefront orthonormal coordinates, can also be used if we work in ray-centered coordinates, and vice versa.

Prestack Gaussian-beam depth migration in anisotropic media

Gaussian-beam depth migration is a useful alternative to Kirchhoff and wave-equation migrations. It overcomes the limitations of Kirchhoff migration in imaging multipathing arrivals, while retaining its efficiency and its capability of imaging steep dips with turning waves. Extension of this migration method to anisotropic media has, however, been hampered by the difficulties in traditional kinematic and dynamic ray-tracing systems in inhomogeneous, anisotropic media. Formulated in terms of elastic parameters, the traditional anisotropic ray-tracing systems aredifficult to implement and inefficient for computation, especially for the dynamic ray-tracing system. They may also result inambiguity in specifying elastic parameters for a given medium.To overcome these difficulties, we have reformulated the ray-tracing systems in terms of phase velocity.These reformulated systems are simple and especially useful for general transversely isotropic and weak orthorhombic media, because the phase velocities for these two types of media can be computed with simple analytic expressions. These two types of media also represent the majority of anisotropy observed in sedimentary rocks. Based on these newly developed ray-tracing systems, we have extended prestack Gaussian-beam depth migration to general transversely isotropic media. Test results with synthetic data show that our anisotropic, prestack Gaussian-beam migration is accurate and efficient. It produces images superior to those generated by anisotropic, prestack Kirchhoff migration.

Ray propagator matrices in three-dimensional anisotropic inhomogeneous layered media

SUMMARY Ray propagator matrices contain the complete solutions to the system of dynamic ray-tracing (DRT) equations connected with a given reference ray. They play an important role in studying the properties of complete four-parameteric systems of paraxial rays and offer many applications in both numerical modelling and practical interpretational problems of seismic ray fields in the high-frequency asymptotic approximation. Traditionally, ray propagator matrices have been expressed either in Cartesian or in ray-centred coordinates, connected with a reference ray. Both coordinate systems have certain advantages. For ray-centred coordinates, the dimensions of ray propagator matrices can be easily reduced from 6 × 6 to 4 × 4 (in a 3-D medium), thereby reflecting the strictly four-parameteric nature of a general paraxial ray field. On the other hand, in Cartesian coordinates, the computations are conceptually simpler and generally valid in isotropic and anisotropic media. In a Cartesian coordinate system, the DRT system and ray propagator matrices are well known both for isotropic and anisotropic layered media. In ray-centred coordinates, the DRT system and ray propagator matrices are known for isotropic layered structures and for anisotropic smooth media, but not for anisotropic media with structural interfaces. We propose a simple and invertible transformation between ray propagator matrices in both coordinate systems. It allows to perform conventional DRT in Cartesian coordinates, and to transform the resulting ray propagator matrix to ray-centred coordinates at any point of the ray where we need it. This avoids DRT in ray-centred coordinates altogether. Vice versa, we can compute, at any point of the reference ray, the ray propagator matrix in Cartesian coordinates by DRT in ray-centred coordinates. We propose several alternative versions of the transformation, each of them equally valid in isotropic and anisotropic media. For rays hitting a structural interface, the ray propagator matrix has to be transformed across the interface. The relevant transformation matrix is usually referred to as the interface propagator matrix. In Cartesian coordinates, the 6 × 6 interface propagator matrix has been published before, but in ray-centred coordinates only for the case of isotropic media. Based on the transformation from Cartesian to ray-centred coordinates, we present the 6 × 6 and 4 × 4 interface propagator matrices in ray-centred coordinates, valid for general isotropic and anisotropic media. The 4 × 4 interface propagator matrix in ray-centred coordinates can be factorized and used to derive the 4 × 4 surface-to-surface paraxial matrices. These matrices allow to relate the paraxial ray properties at different surfaces crossed by the reference ray and offer many important applications.

Reformulated Kinematic and Dynamic Ray Tracing Systems for Arbitrarily Anisotropic Media

I reformulate well-established systems for kinematic and dynamic ray tracing in 3D heterogeneous media with arbitrary anisotropy. Matrices of size 3 × 3, e.g. the Christoffel matrix, are substituted by six-component vectors, and the Christoffel matrix elements are expressed explicitly in terms of the elements of the 6 × 6 matrix of elastic coefficients, written in Voigt notation. Thereby, I find it easier to see the effects on the ray tracing systems of vanishing elastic coefficients and Christoffel matrix elements and of vanishing derivatives with respect to spatial coordinates and slowness vector components. The eigenvalue of the current wave and its derivatives with respect to the ray parameters are included explicitly, which may be favorable for ray tracing processes optimized with respect to speed. With the ANRAY program as a reference, I show that the new formulation requires less optimization than the conventional one.

Determination of the influence zone for surface wave paths

SUMMARY An approximate description of the zone of influence around the propagation path for a surface wave is provided by investigating the Fresnel zones for the frequency range of interest. The influence zone about surface wave paths, over which surface waves are coherent in phase, is identified as approximately one-third of the width of the first Fresnel zone. A technique called Fresnel-area ray tracing (FRT) for surface waves has been used to estimate this region around the ray path for each frequency. The FRT technique is developed by combining two standard ray tracing methods, i.e. kinematic ray tracing (KRT) and dynamic ray tracing (DRT). To obtain the exact Fresnel area in a laterally heterogeneous structure would require the solution of a large number of KRT equations. In contrast, the FRT approach requires just a few ray tracing calculations. In the first step, the trajectory of the surface wave is computed by solving the KRT system for the phase-velocity distribution at the required frequency. In the next step, the behaviour of rays in the zone surrounding the KRT path is calculated by solving the DRT system twice; once from the source to the receiver and once more from the receiver to the source along the same trajectory. Finally, combining the solutions of these ray tracing systems using paraxial ray theory, the Fresnel area around a central ray can be estimated. Using FRT, stationary-phase fields can be constructed around a central ray path in a laterally heterogeneous structure. The influence zone around the ray path is then estimated from the stationary-phase function with simple assumptions concerning the perturbed wavefield. The estimate of the influence zone can be efficiently calculated in laterally heterogeneous structure by using the FRT technique, and allows an extension of current methods of surface wave analysis, which have commonly been based on geometrical ray theory and on the approximation of great-circle propagation. This approach allows the treatment of finite-width rays as well as deviations in propagation from the great circle induced by moderate lateral heterogeneity as revealed by recent tomography models. Such finite-width rays should be of major benefit in enhancing ray-based surface wave tomography.

Three-dimensional geometrical ray theory and modelling of transmitted seismic energy of data from the Nevada Test Site

This paper presents a geometrically based algorithm for computing synthetic seismograms for energy transmitted through a 3-D velocity distribution. 3-D ray tracing is performed to compute the traveltimes and geometrical spreading (amplitude). The formulations of both kinematic and dynamic ray-tracing systems are presented. The two-point ray-tracing problem is solved by systematically updating the initial conditions and adjusting the ray direction until the ray intersects the specified endpoint. The amount of adjustment required depends on the derivatives of the position with respect to the given starting angles between consecutive rays. The algorithm uses derivatives to define the steepest-descent direction and to update the initial directions. The convergence rate depends on the complexity of the model. Test seismograms compare favourably with those from a 2-D asymptotic ray theory algorithm and a 3-D Gaussian-beam algorithm. The algorithm is flexible in modelling arbitrary source and recorder geometries for various smoothly varying 3-D velocity distributions. The algorithm is further tested by simulating surface-to-tunnel vibroseis field data. Shear waves as well as compressional waves may be approximately included. Application of the algorithm to a data set from the Rainier Mesa of the Nevada Test Site produced a good fit to the transmitted (first arrival) traveltimes and amplitudes, with approximately 15 per cent variation in the local 3-D velocity.

Theory of anisotropic dynamic ray tracing in ray-centred coordinates

A 4×4-propagator matrix formalism is presented for anisotropic dynamic ray tracing, including the propagation across curved interfaces. The computations are organised in the same way as in Cervený's well-known isotropic propagator matrix formalism. Attention is paid to cases where double eigenvalues of the Christoffel matrix result in unstable expressions in the dynamic ray tracing system, but where geometrical spreading is well-behaved.

Vertical seismic profile synthetics by dynamic ray tracing in laterally varying layered anisotropic structures

Dynamic ray tracing (DRT) is important in evaluating high‐frequency seismic wave fields in complicated structures. Two formulations of the DRT equations for inhomogeneous anisotropic media are presented. One of them is represented by the classical DRT system, suggested by Červený 15 years ago. Both systems are specified in Cartesian coordinates. The DRT equations are supplemented with formulae for the transformation of DRT across a smoothly curved interface between two inhomogeneous anisotropic media. Červený's formulation of the DRT is applied to the computation of vertical seismic profile (VSP) synthetics. The results of DRT are used for the evaluation of the geometrical spreading and of the coefficients of the paraxial ray approximation for travel times and ray amplitudes. In addition, the DRT is also used in the interval ray tracing procedure, a procedure searching for rays starling from the source and terminating in a specified interval on a profile. Results of numerical modeling of VSP measurements in a three‐dimensional laterally varying structure consisting of Isotropic and anisotropic layers separated by curved interfaces are presented. Ray diagrams of selected elementary waves, time‐distance curves, and multisource three‐component VSP synthetics generated for two different source locations calculated for the anisotropic model are compared with the results determined for a reference Isotropic model. The latter is obtained by averaging phase velocities of the anisotropic model. This comparison clearly shows the effects of anisotropy and the lateral variation of the model on seismic wave fields. The reliability of the ray synthetics is briefly discussed.

A new approach to Gaussian beams on a sphere: theory and application to long-period surface wave propagation

Summary The subject of this paper is the application of the Gaussian beam method to the propagation of surface waves on a laterally heterogeneous, spherical earth taking fully into account the curvature of the Earth's surface. This is achieved by constructing a ray centred coordinate system on general curved surfaces. The application of this method to a spherical shell leads to a ray centred coordinate system adjusted to the spherical geometry. Solving the elastodynamic equation in this coordinate system by a perturbation ansatz yields two interesting results. (1) The well-known parabolic differential equation for the zeroth order displacement follows directly from a solvability condition imposed on the right-hand side of an inhomogeneous differential equation arising from the perturbation ansatz. (2) The dynamic ray-tracing system which is derived from the parabolic equation is a special case of a general dynamic ray-tracing system valid for arbitrary curved surfaces which can be derived by perturbing the ray equations. It is shown that the geometry of the surface enters only by its intrinsic curvature into the dynamic ray-tracing system. Numerical computations of rays and synthetic seismograms based on this new approach for the model M84C of Woodhouse & Dziewonski show that there can occur considerable deviations of the ray-paths from great circles leading to visible amplitude anomalies in the synthetic seismograms. Comparison of synthetic seismograms with selected data recorded at the Black Forest Observatory (BFO) indicate that the model M84C is still too smooth to yield satisfactory agreement of the synthetics with the data.

Applications of dynamic ray tracing

Dynamic ray tracing has recently found broad application in the numerical modelling of high-frequency seismic wave fields in complex 2-D and 3-D layered structures. It involves solving a system of ordinary differential equations (dynamic ray tracing system) along a ray. This paper briefly explains the main principles of dynamic ray tracing in ray centred coordinates, introduces the ray propagator matrix and describes its properties, summarizes the most important applications of dynamic ray tracing and of the ray propagator matrix. Together with known applications, several new applications are discussed in greater detail. Among these, a method to determine the geometrical spreading from the travel-time measurements is proposed.

An application of the propagator matrix of dynamic ray tracing: The focusing and defocusing of body waves by three-dimensional velocity structure in the source region

Summary. Since the dynamic ray tracing system can be written in the linear form dW/ds=SW, it permits the definition of fundamental matrix and propagator matrix solutions. The propagator matrix can be exploited to connect complex 3-D and 2-D regions to 1-D regions of a model. For example, if the three-dimensional variations of a model are confined only to the portion of the ray path between O’ and Oo, numerical integration of the linear system need only be performed between O’ and Oo to construct the propagator π(Oo, O'). Along the segment of the ray path between O’ and Os, traversing a 1-D or homogeneous portion of the model, analytic solutions exist for π(O', Os). Thus π(Oo, Os) can be simply obtained by multiplying the analytic and numerically obtained π matrices. With this procedure, quantities are calculated that are needed for the synthesis of teleseismic P waves by superposition of Gaussian beams. These synthetics are used to investigate the focusing and defocusing of teleseismic P-waves by 3-D structure in the vicinity of the source. A source region, having a fluctuation in velocity of 4 per cent over a characteristic scale length of 50 to 100 km, produces factors of two fluctuations in amplitude and several tenths of a second in travel time at teleseismic range in experiments in which either source location is varied at constant azimuth or in which azimuth is varied at constant source location. A model having a maximum fluctuation as small as 0.8 per cent is capable of producing caustics and multipaths at teleseismic range, depending on its distribution of scale lengths and its ratio of vertical with respect to horizontal scale length. The multipaths and caustics of such a model, however, occur over too small an area and are too closely spaced in arrival time to be resolved with seismograph systems in the 0.01 to 4 Hz band.

Gaussian beams in anisotropic elastic media

We derive the parabolic equation for high-frequency solutions of anisotropic elasticity propagating in the neighbourhood of a ray. The dynamic ray tracing system governing the evolution of a Gaussian beam along its central ray is derived. It is shown that a Gaussian beam is non-singular at the caustics. The decomposition of a point source in an anisotropic medium into Gaussian beams is derived and its application to the complex DRT of seismic waves in anisotropic media is discussed.

Gaussian beams for surface waves in laterally slowly-varying media

Summary. Asymptotic ray theory is applied to surface waves in a medium where the lateral variations of structure are very smooth. Using ray-centred coordinates, parabolic equations are obtained for lateral variations while vertical structural variations at a given point are specified by eigenfunctions of normal mode theory as for the laterally homogeneous case. Final results on wavefields close to a ray can be expressed by formulations similar to those for elastic body waves in 2-D laterally heterogeneous media, except that the vertical dependence is described by eigenfunctions of ‘local’ Love or Rayleigh waves. The transport equation is written in terms of geometrical-ray spreading, group velocity and an energy integral. For the horizontal components there are both principal and additional components to describe the curvature of rays along the surface, as in the case of elastic body waves. The vertical component is decoupled from the horizontal components. With complex parameters the solutions for the dynamic ray tracing system correspond to Gaussian beams: the amplitude distribution is bell-shaped along the direction perpendicular to the ray and the solution is regular everywhere, even at caustics. Most of the characteristics of Gaussian beams for 2-D elastic body waves are also applicable to the surface wave case. At each frequency the solution may be regarded as a set of eigenfunctions propagating over a 2-D surface according to the phase velocity mapping.

Dynamic ray tracing in an anisotropic medium

We derive a dynamic ray tracing system for seismic waves in an anisotropic medium. It consists of ten equations and can be solved numerically. For isotropic media it does not reduce to the DRT systems derived by other authors.

On dynamic ray tracing in three dimensional inhomogeneous media

Among the various dynamic ray tracing systems described by Cervený and Hron[1] is one particular linear system of second order that readily provides identical parameters (in the ray centred coordinate system) to those that fall out of the system by Popov and Psencik[3, 4]. Hence the initial conditions of the latter system for sources and interfaces can easily be used to provide those for the linear system of second order.

The ray series method and dynamic ray tracing system for three-dimensional inhomogeneous media

The ray series method is used to study propagation of seismic waves in the three-dimensional media consisting of generally inhomogeneous layers separated by curved interfaces. The investigation is carried out with the help of the so-called ray centered corrdinate system which was proposed by Popov and Psencik (1978a). It is shown that in this system the principal components of the amplitude coefficients in the ray series for S waves do not rotate about the ray with respect to the basis vectors when the wave progresses, even though they rotate with respect to the unit vectors n^ and b^ along the direction of the normal and binormal to the ray, respectively. This considerably simplifies the final expressions for the amplitude coefficients for S waves, whose two principal components are decoupled in the ray-centered coordinate system. The ray-centered coordinate system is also applied to the eikonal equation in order to produce a dynamic ray tracing system consisting of three nonlinear ordinary differential equations of the first order determining the second derivatives of the time field and, in this manner, even the basic geometrical properties of the wave fronts (e.g., principal curvatures and geometrical spreading) along the ray. Several different modifications of the dynamic ray tracing system are presented. It is demonstrated that in the case of generally inhomogeneous two-dimensional media the dynamic ray tracing system reduces, under certain not too restrictive conditions, to the single first order differential equation of the Riccati type. Finally, the phase matching method is used to determine discontinuities of individual quantities in the dynamic ray tracing system when the wave is impinging on a curved interface separating two generally inhomogeneous media. Since all basic equations are presented in a computationally convenient matrix formulation, they can be readily employed for any numberical evaluation of dynamic properties of seismic waves propagating through structurally complicated media. As the paper describes all basic features of asymptotic ray theory (the name under which the ray series method is known on this continent), it can serve as a starting point for anyone wishing to develop computer programs for the computation of synthetic seismograms.