Dynamic Ray Tracing Research Articles

2.5-D true‐amplitude Kirchhoff migration to zero offset in laterally inhomogeneous media

The proposed new Kirchhoff‐type true‐amplitude migration to zero offset (MZO) for 2.5-D common‐offset reflections in 2-D laterally inhomogeneous layered isotropic earth models does not depend on the reflector curvature. It provides a transformation of a common‐offset seismic section to a simulated zero‐offset section in which both the kinematic and main dynamic effects are accounted for correctly. The process transforms primary common‐offset reflections from arbitrary curved interfaces into their corresponding zero‐offset reflections automatically replacing the geometrical‐spreading factor. In analogy to a weighted Kirchhoff migration scheme, the stacking curve and weight function can be computed by dynamic ray tracing in the macro‐velocity model that is supposed to be available. In addition, we show that an MZO stretches the seismic source pulse by the cosine of the reflection angle of the original offset reflections. The proposed approach quantitatively extends the previous MZO or dip moveout (DMO) schemes to the 2.5-D situation.

Comparison of weights in prestack amplitude‐preserving Kirchhoff depth migration

Three different theoretical approaches to amplitude‐preserving Kirchhoff depth migration are compared. Each of them suggests applying weights in the diffraction stack migration to correct for amplitude loss resulting from geometric spreading. The weight functions are given in different notations, but as is shown, all of these expressions are similar. A notation that is well suited for implementation is suggested: entirely in terms of Green's function quantities (amplitudes or point‐source propagators). For the most common prestack configurations (common‐shot and common‐offset) and 3-D, 2.5-D, and 2-D migrations, expressions of the weights are given in this notation. The quantities needed for calculation of the weights can be computed easily, e.g., by dynamic ray tracing.

True‐amplitude common‐shot migration revisited

In Kirchhoff‐type migration, two dynamic ray‐tracing computations are usually needed for computing the complex weighting (Green's) functions necessary for recovering the source pulse with true amplitude. One computation is from the source point to the image point, the other is from the receiver point to the image point. Since it is a time‐consuming procedure, dynamic ray tracing is a main factor slowing down the performance speed of weighted diffraction stack migration. Here, the known weighting function for a common‐shot configuration is revisited and a new, alternative formula is developed. Because only the takeoff angles of rays are involved in this alternative formula, the module of the complex weighting function can be computed solely by kinematic ray tracing. Further, it is shown that the phase (caustic) correction is not essential for the stack process. As a consequence, the weighted diffraction stack migration can be implemented without using dynamic ray tracing at all. In other words, the subsurface structure can be imaged without using any Green's function. Therefore, using the new formula may accelerate the performance of the Kirchhoff‐type migration, especially in 3-D cases. In addition, the new formula may affect the model smoothing process necessary for using some traveltime computing methods based on ray tracing. As is known, kinematic quantities associated with a given ray are less sensitive to the velocity distribution than the dynamic ones. Thus, the new formula allows one to use a model less smoothed than that demanded for using dynamic ray tracing. As a result, less smoothing operation is needed by building a velocity model without interfaces. The latter points are vital for the accuracy and efficiency of the ray tracers for computing the traveltime of point‐diffracted rays.

Estimation of multivalued arrivals in 3D models using wavefront construction—Part I1

AbstractA new 3D wavefield modelling approach based on dynamic ray tracing is presented. This approach is called wavefront construction, and it can be used in 3D models with constant or smoothly varying material properties (S‐ and P‐velocity and density) separated by smooth interfaces. Wavefronts consisting of rays arranged in a triangular network are propagated stepwise through the model. At each time step, the differences in a number of parameters are checked between each pair of rays on the wavefront. New rays are interpolated whenever this difference between pairs of rays exceeds some predefined maximum value. A controlled sampling of the wavefront at all time steps is thus obtained. Receivers are given multiple‐event values by interpolation when the wavefronts pass them. The strength of the wavefront construction method is that it is robust and efficient.

Wavefield smoothing and the effect of rough velocity perturbations on arrival times and amplitudes

Geometric ray theory is an extremely efficient tool for modelling wave propagation through heterogeneous media. Its use is, however, only justified when the inhomogeneity satisfies certain smoothness criteria. These criteria are often not satisfied, for example in wave propagation through turbulent media. In this paper, the effect of velocity perturbations on the phase and amplitude of transient wavefields is investigated for the situation that the velocity perturbation is not necessarily smooth enough to justify the use of ray theory. It is shown that the phase and amplitude perturbations of transient arrivals can to first order be written as weighted averages of the velocity perturbation over the first Fresnel zone. The resulting averaging integrals are derived for a homogeneous reference medium as well as for inhomogeneous reference media where the equations of dynamic ray tracing need to be invoked. The use of the averaging integrals is illustrated with a numerical example. This example also shows that the derived averaging integrals form a useful starting point for further approximations. The fact that the delay time due to the velocity perturbation can be expressed as a weighted average over the first Fresnel zone explains the success of tomographic inversions schemes that are based on ray theory in situations where ray theory is strictly not justified; in that situation one merely collapses the true sensitivity function over the first Fresnel zone to a line integral along a geometric ray.

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.

Calculating received time-domain amplitude and phase from 3-D ray tracing results

For long-range sound propagation simulations, ray tracing has many advantages as well as drawbacks. Some of the advantages ray tracing offers are 3-D implementation with present-day computers and the capability of propagating only the acoustic field that will contribute to a single early arrival. Well-known drawbacks to ray tracing computations are the neglect of diffraction effects, infinite-amplitude results at caustics, and the often quoted ‘‘high-frequency approximation.’’ However, other simulation techniques have advantages and drawbacks as well, suggesting the community should use many of the available propagation techniques and build strong conclusions through assimilation of results. Here, a ray-endpoint density method is investigated for calculating the time-domain amplitude of a received early arrival. Also, computation of the received time-domain phase is studied. The methods are investigated using 2-D range-invariant, 2-D range variant, and 3-D time and space variant computational ocean sound-speed models. The ray-endpoint density method is compared against normal mode, PE, and dynamic ray tracing results. [Work supported by ONR-ASSERT.]

Time-domain modelling of PKIKP precursors for constraints on the heterogeneity in the lowermost mantle

SUMMARY Compact solutions for the scattering of plane elastic waves by spheres now exist that are valid for arbitrary levels of velocity and density perturbations and arbitrary ratios of wavelength/sphere radius. These solutions are easily incorporated into methods of seismogram synthesis that remain valid for the frequency dependence near caustics. Short-period precursors to PKIKP are synthesized by dynamic ray tracing and superposition of Gaussian beams in a representation of D″ heterogeneity by distributed spheres. Corrections to Rayleigh-Born scattering for arbitrary levels of perturbation, wavelengths on the order of or smaller than the scatterer dimension, and scattering into inner core branches of PKP tend to cancel one another, validating previous work on PKIKP precursors that assumed the validity of the Rayleigh-Born approximation. Results are similar to those studies: scalelengths on the order of 20–35 km and perturbations in P velocity and density of 10 per cent. Modelling in the time domain is helpful in interpreting the emergent arrival of the precursor wave train, constraining the depth of D″ heterogeneity. The length and excitation of precursor coda in the 120° to 130° range helps to constrain the maximum thickness of D″ heterogeneity. An apparent increase in precursor amplitude as distance increases toward the B caustic suggests a scattering mechanism that is more strongly concentrated in the forward direction than is predicted by models containing spherical scatterers or by isotropic Gaussian distributions having a single scalelength.

3-D crosswell transmissions: Paraxial ray solutions and reciprocity paradox

Transmission of seismic waves through a 3-D earth model is of fundamental importance in seismology. If the model consists of many layers separated by curved interfaces, the only feasible solution to the transmitted waves is the one given by the geometrical optics approximation. Transmitted rays, transmitted wavefield, and the first Fresnel zone associated with a transmission point can be expressed by four 2 × 2 constant matrices constituting the 4 × 4 linearized ray transformation matrix. Generally, the ray transformation matrix can be constructed by dynamic ray tracing. However, if the layers are homogeneous, it can be formulated in closed form by using elementary vector calculus and coordinate transformations. Using the symplecticity of the ray transformation matrix, transmissions in opposite directions can be formulated by the ray transformation matrix for only one direction, and a reciprocity relation can be established. After the decomposition theorem for the ray transformation matrix, the transmitted wavefield in a model with many curved interfaces can be computed in a cascaded way. Using the [Formula: see text]‐matrix decomposition theorem, the normalized geometrical spreading factor can be expressed by means of the area of the first Fresnel zone of a transmission point. If seismic waves propagate through a locally spherical interface, the reciprocity relation may not hold. Using ray theory, this fact is shown by formulating the transmitted wavefield with the two principal radii of curvature of the transmitted wavefront at the transmission point under consideration. Using wave theory, this fact is shown by analyzing the Debye integral with the method of stationary phase.

Transformations for dynamic ray tracing in anisotropic media

Six-dimensional dynamic ray tracing in (phase-space) Cartesian coordinates was introduced by Červený [Geophys. J.R. Astr. Soc. 29, 1–13 (1972)]. Hanyga [Tectonophysics 90, 243–251 (1982)] showed that it reduces to 4-dimensional dynamic ray tracing in (phase-space) ray-centered coordinates. This paper concentrates on the explicit transformation equations of dynamic ray tracing between Cartesian and ray-centred coordinates. Many of the transformation equations have not been published before even for isotropic medium. Also proposed is an efficient way of reducing the number of equations being solved when numerically evaluating the paraxial-ray propagator matrices, both in Cartesian and ray-centred coordinates.

Numerical simulation of high- and low-frequencyLg-wave propagation

SUMMARY Lg waves are sometimes extinguished when traversing geological structures such as mountain ranges or grabens. Previous investigations have applied numerical, full-waveform simulations at relatively low frequencies (on the order of 1 Hz) to investigate the relationship of different structural models and the amplitude anomalies of the Lg waves. These studies often fail to link specific known structures to the amplitude anomalies. However, the observed waveforms also often contain data of significantly higher frequency than the full waveform simulations. We compare such low-frequency simulations computed with the boundary-integral-equation (BIE) technique with synthetic seismograms generated by dynamic ray tracing (DRT). Beginning with a plane-layered model of the crust in central France, we show that the general characteristics of the Lg wave train over distances ranging from 100 km to around 500 km can be modelled by computing a relatively small number of rays reflected between the free surface and the Moho. Matching the details of the waveforms is more difficult, but can be improved by the addition of reflections between other interfaces in the model. Considering a one-layer model of the crust with a step-like change in Moho depth that leads to a rapid decrease in crustal thickness, we show again that the BIE and DRT solutions are very similar in character. Significantly, both solutions predict that the Lg wave will not be extinguished on traversing the Moho step. Finally, we consider a more complete and detailed model of the Pyrenees. Though it is more difficult to match exactly the ray-and integral-equation solutions, the general trends continue to be the same, and the Lg waves are still predicted to cross the tectonic structures associated with the mountains. Therefore, it appears that observed amplitude changes must be caused by something other than the general structure of the Pyrenees, such as scattering by smaller scale features within the lower crust.

The continent‐ocean crustal transition across the southwest Greenland margin

This paper examines the complete crustal transition across the nonvolcanic, southwest Greenland continental margin of the Labrador Sea using wide‐angle and coincident vertical‐incidence seismic profiles. Six ocean bottom seismometers and a sonobuoy record P and S wave first and multiple arrivals from the crust and upper mantle, which are analyzed by two‐dimensional dynamic ray tracing and one‐dimensional reflectivity modeling. The resulting seismic velocity model requires that the preexisting 30‐km thick continental crust is thinned abruptly to ∼3 km across the continental slope, primarily by removal of the lower crust. Farther seaward, the crust thickens to ∼6 km primarily through the addition of a high‐velocity (7.0–7.6 km/s) layer in the lower crust. This lower crustal layer is 4–5 km thick, extends for a horizontal distance of ∼80 km, and is interpreted as partially serpentinized upper mantle. It is overlain by a low‐velocity (4.0–5.0 km/s), upper layer which is interpreted as highly fractured upper continental crust. Our model suggests that seafloor spreading did not start until chrons 27–28, 13 Ma younger than previously suggested. This interpretation is supported by two‐dimensional modeling of gravity and magnetic data along the refraction line. Our results are consistent with a simple shear mechanism for the initial rifting, with the SW Greenland margin as the upper plate. However, a full characterization of the rifting mechanism must await comparison with a seismic model for the conjugate margin, east of Labrador.

Multiple weights in diffraction stack migration

Three‐dimensional (3-D) prestack diffraction‐stack migration methods (often called Kirchhoff migration/inversion) play a fundamental role in seismic imaging. In addition to estimating the location of arbitrarily curved reflectors and the angle‐dependent reflection coefficients upon them, they can also be used to provide useful kinematic and dynamic information about the specular reflection ray that connects the source and receiver via the unknown reflecting interface. This is achieved by performing a diffraction stack more than once upon the same seismic data set using identical stacking surfaces but different weights. Some of these weights can be applied simultaneously, i.e., as a weight‐vector. The approach offers the possibility of determining various useful quantities that help to compute and interpret migrated reflections. The vector‐weighted diffraction stack is principally intended to economize the amplitude‐preserving migration that normally would require a large amount of dynamic ray tracing. A simple 2-D synthetic example shows how the method works in principle.

3-D Salt and Sub-salt Imaging Strategy: A Case History from the Gulf of Mexico

The purpose of this paper is to demonstrate 3-D prestack depth migration as an ultimate technique for imaging structural targets associated with salt diapirism, and compare it with other migration strategies ? 2-D prestack, and 2-D/3-D poststack depth migrations. The structural problem is to image sedimentary bed terminations, bottom salt boundaries and undersalt layers around two large salt masses in the Gulf of Mexico. A large number of oil and gas reservoirs in the Gulf of Mexico is controlled by salt structures. Typical salt flank plays have been mapped quite successfully with time migration strategy. However, as exploration moves into more complicated targets around salt overhangs and subsalt traps, depth migration strategy becomes necessary. The first task is to build a detailed 3-D velocity-depth model. To accomplish this, the following robust approach has been chosen. Using a velocity-depth model for the overburden (sedimentary section) based on edited stacking velocities, 3-D poststack depth migration was performed to verify its accuracy. Iterative 3-D poststack depth migration then followed to define the salt masses. Top salt reflectors were interpreted in 3-D and constant salt velocity was assigned to the half-space below the top-salt. Again, 3-D poststack depth migration was performed and base-salt events were interpreted. Finally, the velocity model was updated to contain complete salt bodies and 3-D poststack migration was rerun to check results. This loop of three iterations was repeated until convergence was achieved ? the input model matched the result of the depth migration. Given the velocity model, a single output line was selected and 3-D dynamic raytracing was performed to link output image points with surface shot-geophone locations. Finally, 3-D Kirchhoff depth migration was run in 3-D prestack, 2-D prestack and 2-D poststack modes to produce comparison sections along the output line.

Application of dynamic ray tracing in the 3-D inversion of seismic-reflection data

Possible applications of dynamic ray tracing to the 3-D Born and Kirchhoff inversion of seismic-reflection data are discussed. It is shown that the most important quantities in the inversion integrals, such as the Beylkin determinant and the amplitudes of the ray-theory Green functions, can be determined using the dynamic ray-tracing procedure. Dynamic ray tracing can be simply performed along known rays in a very general, laterally varying, layered background medium. The algorithms do not require the determination of the derivatives of traveltimes which play an important role in some other methods. They will find valuable applications particularly in complex 3-D structures.

Three‐dimensional primary zero‐offset reflections

Zero‐offset reflections resulting from point sources are often computed on a large scale in three‐dimensional (3-D) laterally inhomogeneous isotropic media with the help of ray theory. The geometrical‐spreading factor and the number of caustics that determine the shape of the reflected pulse are then generally obtained by integrating the so‐called dynamic ray‐tracing system down and up to the two‐way normal incidence ray. Assuming that this ray is already known, we show that one integration of the dynamic ray‐tracing system in a downward direction with only the initial condition of a point source at the earth’s surface is in fact sufficient to obtain both results. To establish the Fresnel zone of the zero‐offset reflection upon the reflector requires the same single downward integration. By performing a second downward integration (using the initial conditions of a plane wave at the earth’s surface) the complete Fresnel volume around the two‐way normal ray can be found. This should be known to ascertain the validity of the computed zero‐offset event. A careful analysis of the problem as performed here shows that round‐trip integrations of the dynamic ray‐tracing system following the actually propagating wavefront along the two‐way normal ray need never be considered. In fact some useful quantities related to the two‐way normal ray (e.g., the normal‐moveout velocity) require only one single integration in one specific direction only. Finally, a two‐point ray tracing for normal rays can be derived from one‐way dynamic ray tracing.

SEISMIC REFLECTION AMPLITUDES1

AbstractSeismic amplitude variations with offset contain information about the elastic parameters. Prestack amplitude analysis seeks to extract this information by using the variations of the reflection coefficients as functions of angle of incidence. Normally, an approximate formula is used for the reflection coefficients, and variations with offset of the geometrical spreading and the anelastic attenuation are often ignored. Using angle of incidence as the dependent variable is also computationally inefficient since the data are recorded as a function of offset.Improved approximations have been derived for the elastic reflection and transmission coefficients, the geometrical spreading and the complex travel‐time (including anelastic attenuation). For a 1 D medium, these approximations are combined to produce seismic reflection amplitudes (P‐wave, S‐wave or converted wave) as a Taylor series in the offset coordinate. The coefficients of the Taylor series are computed directly from the parameters of the medium, without using the ray parameter.For primary reflected P‐waves, dynamic ray tracing has been used to compute the offset variations of the transmission coefficients, the reflection coefficient, the geometrical spreading and the anelastic attenuation. The offset variation of the transmission factor is small, while the variations in the geometrical spreading, absorption and reflection coefficient are all significant.The new approximations have been used for seismic modelling without ray tracing. The amplitude was approximated by a fourth‐order polynomial in offset, the traveltime by the normal square‐root approximation and the absorption factor by a similar expression. This approximate modelling was compared to dynamic ray tracing, and the results are the same for zero offset and very close for offsets less than the reflector depth.

A lithospheric velocity anomaly beneath the Shagan river Test Site. Part 2. Imaging and inversion with amplitude transmission tomography

AbstractA linear inversion procedure is introduced that images weak velocity anomalies using amplitudes of transmitted seismic waves. Using projection operators from geometrical ray theory, an image of an anomaly is constructed from amplitudes recorded at arrays of receivers using arrays of sources. The image is related to the velocity anomaly by a second-order partial-differential equation that is inverted using 2-D discrete Fourier transforms.As an example of the inversion procedure, magnitude residuals for European stations recording Shagan River explosions are used to image the deep lithospheric anomaly beneath the Shagan River test site described in Part 1. This formal inversion analysis confirms the existence of a small-scale lateral heterogeneity located 50 km west-northwest of the test site at a probable depth between 80 and 100 km and indicates that it is consistent with a deterministic 1.5% peak-to-peak (or 0.5% rms) velocity anomaly with a scale length of about 3 km. 3-D dynamic raytracing is then used to verify that the inferred laterally varying structure produces amplitude fluctuations consistent with observations.

The times and amplitudes of core phases for a variable core-mantle boundary layer

SUMMARY Synthetic pulse forms of core phases have been computed in models of an undulating core-mantle boundary or an undulating boundary layer at the base of the mantle. A simple detector has been applied to the pulse forms of PcP, PKP (BC and AB) and PKKP (BC). The detections are influenced by amplitude and waveform changes due to focusing, defocusing and interference. Numerical results are given for the increased bias and the change of variance of the traveltime residuals. The results for traveltime bias follow the observed trend in real data of core phases, although the computed bias is less than observations suggest. Special attention is paid to the underlying assumptions of seismic tomography. It is shown that these assumptions are not tenable for several of the phases and models considered here. In those cases, seismic tomography will produce models that are in poor agreement with the true model. The synthetics were computed by introducing the concepts of 3-D dynamic ray tracing with a Kirchhoff type of integral over a spherical reference surface. Simple expressions are given for propagating the ray propagators of 3-D dynamic ray tracing through spherically symmetric parts of the model.

Parameter estimation in a one‐dimensional anelastic medium

Surface seismic data with a known source pulse are used to estimate parameters in a stack of homogeneous layers. The unknown parameters in each layer are the P wave velocity, the S wave velocity, the density, the layer thickness, and the quality factor. The forward problem is solved by dynamic ray tracing. Primary and multiple P wave reflections are considered. The S wave velocities are estimated from the PP transmission and reflection coefficients. The model is fitted to the data using a stabilized least squares procedure. The Jacobian is computed analytically in each step. By introducing the layer thicknesses as unknown parameters, it is possible to keep the number of unknowns at a tractable level, and the Gauss‐Newton approximation to the Hessian matrix can be inverted exactly at each iteration. When only one event is considered, the scheme provides a nonlinear solution to the amplitude versus offset problem, where offset dependent corrections for absorption, geometrical spreading, and transmission losses are automatically included. All parameters are successfully estimated from synthetic data. We also perform a sensitivity analysis for the different parameters and discuss the reliability for the estimates in each parameter group.