Asymptotic Ray Theory Research Articles

Representation of discontinuous seismic velocity fields by sigmoidal functions for ray tracing and traveltime modelling

SUMMARYWave-modelling methods based on asymptotic ray theory have a lower computational cost than full wave-equation methods but require a smooth velocity field, though discontinuities may be handled by imposing interface conditions between adjacent blocks. We propose to approximate discontinuous velocity fields with model parametrizations based on smooth, rapidly varying functions known as sigmoidal functions. We have implemented the proposed technique on Cartesian grids using the wavelet theory formalism. Numerical experiments with 2-D and 3-D initial-value and two-point ray tracing in heterogeneous media show that the ray paths and traveltimes produced with the sigmoidal representation are consistent with the results produced by conventional ray tracing in block structures, broadening the scope of classical algorithms based on smooth velocity fields.

Seismic imaging of subsurface geological structures by Kirchhoff’s migration based on extended Born approximation

Complex seismic signatures are generated because of the multifaceted nature of the subsurface. These features make the interpretation very complex. To understand the seismic behaviour, different numerical tools are available. In this present study, an attempt has been made to demonstrate both the modelling and imaging aspects of these complex subsurface features commonly encountered in petroleum exploration. The present work is an extended form of the Born approximation by using Green’s function based asymptotic ray theory. Subsequently, Kirchhoff’s depth migration has been applied to generate seismic shot gathers over structural as well as stratigraphic traps. From this analysis, it is observed that the technique is able to efficiently migrate both the structural and stratigraphic traps. The proposed technique also intends to handle strong velocity variation and amplitude restoration. However, some noise in terms of over-critical reflection has been observed in the depth migrated section corresponding to pinch-out and unconformity respectively.

A trajectory mechanics approach for the study of wave propagation in an anisotropic elastic medium

SUMMARYWe derive equations describing the path and traveltime of a coherent elastic wave propagating in an anisotropic medium, generalizing expressions from conventional high-frequency asymptotic ray theory. The methodology is valid across a broad range of frequencies and allows for subwavelength variations in the material properties of the medium. The primary difference from current ray methods is the retention of a term that is neglected in the derivation of the eikonal equation. The additional term contains spatial derivatives of the properties of the medium and of the amplitude field, and its presence couples the equations governing the evolution of the amplitude and phase along the trajectory. The magnitude of this term provides a measure of the validity of expressions based upon high-frequency asymptotic methods, such as the eikonal equation, when modelling wave propagation dominated by a band of frequencies. In calculations involving a layer with gradational boundaries, we find that asymptotic estimates do deviate from those of our frequency-dependent approach when the width of the layer boundaries become sufficiently narrow. For example, for a layer with boundaries that vary over tens of meters, the term neglected by a high-frequency asymptotic approximation is significant for frequencies around 10 Hz. The visible differences in the paths of the rays that traverse the layer substantiate this conclusion. For a velocity model derived from an observed well log, the majority of the trajectories calculated using the extended approach, accounting for the frequency-dependence of the rays, are noticeably different from those produced by the eikonal equation. A suite of paths from a source to a specified receiver, calculated for a range of frequencies between 10 and 100 Hz, define a region of sensitivity to velocity variations and may be used for an augmented form of tomographic imaging.

A Crustal Velocity Model for the Peninsular Ranges of Baja California and Southwestern Laguna Salada, Mexico

To see any change in seismic velocities that may be associated with an abrupt change in the regional geology (granitic rock in contact with sediments), we conducted a refraction seismic study in the Peninsular Ranges of Baja California, which is in the Mexico–southwestern Laguna Salada (LS) region. We installed 30 three‐component portable seismic stations, supplemented with two permanent six‐component stations of the Northwest Mexico Seismic Network (RESNOM). The stations, spaced ∼6 km along a refraction profile, recorded two blasts; these were the direct shot located to the south of the city of Ensenada and the reverse shot in the southwestern LS (southwest–northeast direction). Record sections show seismograms with impulsive P arrivals at nearby stations. Rays from the two blasts were modeled (using asymptotic ray theory) to obtain a P‐wave velocity model from 0 to ∼15 km depth along the refraction profile. Our modeling results are as follows: in the southwestern part of the profile (0–25 km distance), a low‐velocity zone of ∼2 km/s exists between the depths of 0 and 3.5 km; in Sierra Juarez, the mean P‐wave velocity is ∼5.6 km/s between the depths of 0 and 5 km; and in southwestern LS, a low‐velocity layer of ∼2.5 km/s exists between the depths of 0 and ∼3 km. We also modeled a layer of ∼6.5 km/s between 4 and 12 km in the Ensenada–Ojos Negros region, and between the depths of 4 and 8 km below the southwestern LS. From a profile distance of 0 to 50 km, a velocity zone of ∼6.7 km/s appears between the depths of 12 and 15 km.

Bayesian seismic inversion for estimating fluid content and fracture parameters in a gas-saturated fractured porous reservoir

Understanding the effects of in situ fluid content and fracture parameters on seismic characteristics is important for the subsurface exploration and production of fractured porous rocks. The ratio of normal-to-shear fracture compliance is typically utilized as a fluid indicator to evaluate anisotropy and identify fluids filling the fractures, but it represents an underdetermined problem because this fluid indicator varies as a function of both fracture geometry and fluid content. On the bases of anisotropic Gassmann’s equation and linear-slip model, we suggest an anisotropic poroelasticity model for fractured porous reservoirs. By combining a perturbed stiffness matrix and asymptotic ray theory, we then construct a direct relationship between the PP-wave reflection coefficients and characteristic parameters of fluids (P- and S-wave moduli) and fractures (fracture quasi-weaknesses), thereby decoupling the effects of fluid and fracture properties on seismic reflection characterization. By incorporating fracture quasi-weakness parameters, we propose a novel parameterization method for elastic impedance variation with offset and azimuth (EIVOA). By incorporating wide-azimuth observable seismic reflection data with regularization constraints, we utilize Bayesian seismic inversion to estimate the fluid content and fracture parameters of fractured porous rocks. Tests on synthetic and real data demonstrate that fluid and fracture properties can be reasonably estimated directly from azimuthal seismic data and the proposed approach provides a reliable method for fluid identification and fracture characterization in a gas-saturated fractured porous reservoir.

Radiation properties of a uniaxial chiral quadratic inhomogeneous slab under oblique incidence

In this paper, the focussing of the transmitted electromagnetic field through a quadric inhomogeneous slab of lossless uniaxial chiral medium is derived using transmission coefficient under oblique incidence. The inhomogeneity in the uniaxial chiral slab has been incorporated through permittivity parameter. Asymptotic ray theory provides valid field everywhere except at focal point where it gives infinite value. Singularity of the field at focal point is addressed using Maslov's method. The derived analytical field expressions at caustic or focal point of uniaxial quadratic inhomogeneous slab have been solved numerically using MATHEMATICA. The effects of chirality parameter, axial permittivity, transvers permittivity, angle of incidence on the refracted field are discussed and the effects of Brewster angle on the focussed field are also discussed. The results obtained using Maslov's method are compared with Huygens–Kirchhoff's integral which are in good agreement.

Simulation of selected strong motion records of the 2003 MW = 6.6 Bam earthquake (SE Iran), the modal summation-ray tracing methods in the WKBJ approximation

A new analytical methodology for computing synthetic seismograms in 3-D anelastic media is applied to the modelling of the local records of the 2003 Bam Mw 6.6 earthquake, Southeastern Iran. The method is based on the combination of the modal summation technique with the asymptotic ray theory. The 3-D models are determined by a set of vertically heterogeneous sections (1-D structures) that are juxtaposed on a regular grid. The distribution of these sections in the grid is done in such a way to satisfy the condition of applicability of the WKBJ approximation, that is, the lateral variation of all the elastic parameters has to be small with respect to the prevailing wavelength. The procedure, described very briefly, has been applied for the area of the destructive 2003 December 26 Bam earthquake Mw 6.6, with a maximum cut-off frequency of 6 Hz. Synthetic seismograms show fine agreement with recorded signals when comparing shapes, amplitudes and spectra and evidence the relevance of 3-D effects.

Finite-frequency structural sensitivities of short-period compressional body waves

SUMMARY We present an extension of the method recently introduced by Zhao & Chevrot for calculating Frechet kernels from a precomputed database of strain Green's tensors by normal mode summation. The extension involves two aspects: (1) we compute the strain Green's tensors using the Direct Solution Method, which allows us to go up to frequencies as high as 1 Hz; and (2) we develop a spatial interpolation scheme so that the Green's tensors can be computed with a relatively coarse grid, thus improving the efficiency in the computation of the sensitivity kernels. The only requirement is that the Green's tensors be computed with a fine enough spatial sampling rate to avoid spatial aliasing. The Green's tensors can then be interpolated to any location inside the Earth, avoiding the need to store and retrieve strain Green's tensors for a fine sampling grid. The interpolation scheme not only significantly reduces the CPU time required to calculate the Green's tensor database and the disk space to store it, but also enhances the efficiency in computing the kernels by reducing the number of I/O operations needed to retrieve the Green's tensors. Our new implementation allows us to calculate sensitivity kernels for high-frequency teleseismic body waves with very modest computational resources such as a laptop. We illustrate the potential of our approach for seismic tomography by computing traveltime and amplitude sensitivity kernels for high frequency P, PKP and Pdiff phases. A comparison of our PKP kernels with those computed by asymptotic ray theory clearly shows the limits of the latter. With ray theory, it is not possible to model waves diffracted by internal discontinuities such as the core–mantle boundary, and it is also difficult to compute amplitudes for paths close to the B-caustic of the PKP phase. We also compute waveform partial derivatives for different parts of the seismic wavefield, a key ingredient for high resolution imaging by waveform inversion. Our computations of partial derivatives in the time window where PcP precursors are commonly observed show that the distribution of sensitivity is complex and counter-intuitive, with a large contribution from the mid-mantle region. This clearly emphasizes the need to use accurate and complete partial derivatives in waveform inversion.

A Numerical Ray Generation and its Application to the Computation of Synthetic Seismograms for Complex Layered Media

Asymptotic ray theory, which represents a high frequency approximation to the solution of wave equation, is beginning to find its way into geophysical literature dealing with computation of synthetic seismograms. Because in this theory-as well as in any ray theory in general-the total wave field is decomposed into contributions attributed to the individual rays, a rational method of ray generation must be used in order to select from the infinite number of rays, which exist between the source and a receiver in any layered medium, only those with the significant amplitudes. If the medium consists of homogeneous layers separated by parallel plane interfaces the phases which travel along different paths but arrive simultaneously at each epicentral distance can be called kinematic analogues. Some of these have even the same amplitudes and will be called dynamic analogues. Given the number of layers and specifying the maximum number of path elements in a ray, number of both kinematic and dynamic analogues can be determined analytically (in case of few conversions) or with the aid of a computer (if a large number of conversions is required). Unfortunately, the concepts of kinematic and dynamic analogues cannot be used for media with curvilinear, dipping interfaces or in case of generally inhomogeneous layers. A general program, which computes codes of all possible rays for specified numbers of layers, ray paths and conversions (in case of media with curvilinear, dipping interfaces or for inhomogeneous layers) or creates the groups of kinematic and dynamic analogues (for media consisting of parallel homogeneous layers), was written and its results presented. It was shown that with the exception of the most trivial media the correct selection of all significant rays for construction of synthetic seismograms is too difficult for any human being and has to be left to a computer if serious inaccuracies are to avoided. This conclusion was inferred from synthetic seismograms computed for a twelvelayered crust model proposed by I. P. Kosminskaya (see Fig. 1) as well as from the results (see Fig. 4) obtained for a simplified three-layered model suggested for southern Alberta by E. R. Kanasewich and R. M. Clowes (Fig. 2). The first model, which consisted of homogeneous layers separated by parallel plane interfaces, also served as an example of the usefulness of splitting the rays into groups of kinematic and dynamic analogues: although more than 150 OOO refracted, reflected and head phases were examined and their contributions to the seismograms evaluated, it took less than 10 min to an IBM 360/67 computer to evaluate and select the most significant phases in the seismogram at each epicentral distance. On the basis of synthetic seismograms computed for the Alberta model the following conclusions were drawn: (1) A special attention has to be paid to the exact location of cusps on travel-time curves and to the end points of individual branches related to the same ray. As the number of these points and branches is often surprisingly high (see Fig. 3), a rather sophisticated algorithm has to be used in order to make computation reasonably efficient.

漸近波線理論と差分法を用いた弾性波動のハイブリッドシミュレーション

現在，水平方向にも不均質な媒質の弾性波の波動場計算では差分法が頻繁に利用されている。差分法は，アルゴリズムが簡単であり複雑な地下構造にも対応できるという利点を持つ。しかし，浅部に伝播速度の低い構造を持つモデルや広大なモデルに対しての計算には，計算時間や計算機のメモリなどの問題が生じる。一方，波線理論の計算は複雑な構造に対しては適用が難しいが，簡単な構造に対しては計算時間も小さく，計算を効率的に行うことができる。 本研究では，先に上げたようなモデルに対しての差分法の波動伝播計算に関する問題を解決するために，浅部の簡単な構造に対して波線理論を弾性学に応用した漸近波線理論を利用し，それ以外の領域に差分法を用いるという二つの手法を用いたハイブリッドシミュレーションを行った。このシミュレーション手法は差分法と漸近波線理論の両手法の利点を生かすことができ，今後複雑化することが予想される資源探査にも応用できると考えられる。

FOCUSING OF ELECTROMAGNETIC PLANE WAVE INTO UNIAXIAL CRYSTAL BY A THREE DIMENSIONAL PLANO CONVEX LENS

A three dimensional plano-convex lens which is placed at a certain distance from a plane uniaxial interface has been considered. High frequency fields refracted by the geometry are derived. The treatment is based on Maslov's method. The method combines the simplicity of asymptotic ray theory and generality of the transform method to remedy the problem of geometrical optics around the caustic point of a focusing system. Field patterns are obtained which includes the observation points around the caustic region. The results are found in good agreement with obtained using Huygens-Kirchhoff Principle.

Migração Kirchhoff em profundidade de seções de ondas planas

We present a migration method for plane wave source data based on the Kirchhoff integral. Plane wave migration has been employed to increase the efficiency of the wavefield extrapolation. Point source records are synthesized into a plane wave and the resulting plane wave data can be migrated using the same scheme as the one for conventional shot-record migration. Plane wave migration reduces the migration computation time, because the extrapolation does not need to be done for all the individual shot records, but only for the constant plane wave sections. Here we have developed a version of the plane wave migration method based on the Kirchhoff integral. To calculate the traveltimes for the point source and also for plane wave we used a eikonal solver. Once all plane waves section are migrated independently of the others, the final image is obtained by stacking all common images gathers. Integral methods are easier than the wave equation methods to adapt to different recording geometries and also allow large and small models to be treated in the same way. Also, with the aid of the Asymptotic Ray Theory (ART), integral methods can be extended to handle arbitrary velocity structures without any dip limitation. In this paper the SEG-EAGE salt model example and a real dataset from Gulf of Mexico are used to test the theory and show that Kirchhoff plane wave depth migration can produce images comparables in quality to those obtained with Kirchhoff migration and also a reduced computation cost when comparable with the Kirchhoff migration of common shot gathers.

Generalized anelastic asymptotic ray theory

This paper deals with dispersion and attenuation of body waves in a wide range of materials representing realistic rock structures. The time-domain asymptotic ray theory is applied to a new generalized coordinate-free wave equation with an arbitrary tensor relaxation function. This function describes any general viscosity that can occur in 3D anisotropic and laterally inhomogeneous media of practical interest. Perfect elasticity and existing anelastic models (linear viscoelasticity, power-law models, Biot’s fluid flow, thermoelasticity, etc.) are considered as special cases. A ray-series solution is obtained for both the principal and additional components. This solution is suitable for the analysis of anelasticity problems involving propagation of evanescent waves provided that the causality principle is respected. Accuracy of the principal component is numerically tested.

Wave equation calculation of most energetic traveltimes and amplitudes for Kirchhoff prestack migration

Because of its computational efficiency, Kirchhoff migration is the method of choice for 3D prestack depth migration, particularly in the initial velocity model-building stages where several iterations of migration are necessary. In the beginning, Kirchhoff migration used traveltimes calculated by either eikonal solvers or asymptotic ray theory, while amplitudes were calculated using relatively smooth geometrical spreading and obliquity considerations. More recently, Kirchhoff migration has been generalized to include more than one arrival time, with amplitudes calculated from geometrical optics, solution of the transport equations, or Gaussian beams (e.g., Hill, 2001).

Three‐dimensional modeling of Mount Vesuvius with sequential integrated inversion

A new image of Mount Vesuvius and the surrounding area is recovered from the tomographic inversion of 693 first P wave arrivals recorded by 314 receivers deployed along five profiles which intersect the crater, and gravity data collected in 17,598 stations on land and offshore. The final three‐dimensional (3‐D) velocity model presented here is determined by interpolation of five 2‐D velocity sections obtained from sequential integrated inversion (SII) of seismic and gravity data. The inversion procedure adopts the “maximum likelihood” scheme in order to jointly optimize seismic velocities and densities. In this way we recover velocity and density models both consistent with seismic and gravity data information. The model parameterization of these 2‐D models is chosen in order to keep the diagonal elements of the seismic resolution matrix in the order of 0.2–0.8. The highest values of resolution are detected under the volcano edifice. The imaged 6‐km‐thick crustal volume underlies a 25 × 45 km2 area. The interpolation is performed by choosing the right grid for a smoothing algorithm which prepares optimum models for asymptotic ray theory methods. Hence this model can be used as a reference model for a 3‐D tomographic inversion of seismic data. The 3‐D gravity modeling is straightforward. The results of this study clearly image the continuous structure of the Mesozoic carbonate basement top and the connection of the volcano conduit structure to two shallow depressions, which in terms of hazard prevention are the regions through which magma may more easily flow toward the surface and cause possible eruptions.

Asymptotic Ray Theory for Seismic Surface Waves in Laterally Inhomogeneous Media

A recurrence procedure is outlined for constructing asymptotic series for surface wave field in a half-space with weak lateral heterogeneity. Both horizontal variations of the elastic parameters and of the wave field are assumed small on the distances comparable with the wavelength. This is equivalent to the condition that the frequency ω is large. The Surface Wave Asymptotic Ray Theory (SWART) is an analog of the asymptotic ray theory (ART) for body waves. However the case of surface waves presents additional difficulty: the rate of amplitude variation is different in vertical and horizontal directions. In vertical direction it is proportional to the large parameter ω. To overcome this difficulty the transformation ‘equalizing’ vertical and horizontal coordinated is suggested, Z = ωz. In the coordinates x,y,Z the wave field is represented as an asymptotic series in inverse powers of ω. The amplitudes of successive terms of the series are determined from a recurrent system of equations. Attention is paid to similarity and difference of the procedures for constructing the ray series in SWART and ART. Applications of SWART to interpretation of seismological observations are discussed.

An approximate elastic two-dimensional Green’s function for a constant-gradient medium

Summary Approximate analytical formulae for elastic 2-D Green's functions for a constant-gradient propagation velocity medium are presented. These solutions correspond to unit line forces per unit length: the antiplane SH line source and the in-plane P–SV line sources, respectively. They are based on the asymptotic ray theory and account for both near-source effects and low frequencies. The analytical expressions obtained are tested by means of the pseudospectral method with a velocity–stress staggered grid. The agreement of the approximate analytical displacement and stress fields with their numerical counterparts is generally very good. A measure of relative error is proposed.

3-D traveltimes and amplitudes by gridded rays

Seismic imaging in three dimensions requires the calculation of traveltimes and amplitudes of a wave propagating through an elastic medium. They can be computed efficiently and accurately by integrating the eikonal equation on an elemental grid using finite‐difference methods. Unfortunately, this approach to solving the eikonal equation is potentially unstable unless the grid sampling steps satisfy stability conditions or wavefront tracking algorithms are used. We propose a new method for computing traveltimes and amplitudes in 3-D media that is simple, fast, unconditionally stable, and robust. Defining the slowness vector as [Formula: see text] and assuming an isotropic medium, the ray velocity v is related to the slowness vector by the relation [Formula: see text]. Rays emerging from gridpoints on a horizontal plane are propagated downward a single vertical grid step to a new horizontal plane. The components of the slowness vector are then interpolated to gridpoints on this next horizontal plane. This is termed regridding; the process of downward propagation of rays, one vertical grid step at a time, is continued until some prescribed depth is reached. Computation of amplitudes is achieved using a method similar to that for obtaining the zero‐order approximation in asymptotic ray theory. We show comparisons with a full‐wave method on readily accessible 3-D velocity models.

Inferences on the upper crustal structure of Southern Apennines (Italy) from seismic refraction investigations and subsurface data

This paper presents an interpretation of crustal seismic refraction data from the northern sector of the Southern Apennines thrust belt, a region that in historical times experienced large destructive earthquakes. The data were acquired in 1992 along a seismic line 75 km long and parallel to the Apenninic chain, in order to determine a detailed 2-D P-wave velocity model of the upper crust in an area that had not been deeply investigated by geophysical methods previously. We have used a 2-D ray tracing technique based on asymptotic ray theory to model travel times of first and reflected P-wave arrivals. Synthetic seismograms have been produced by finite difference simulations in order to check the reliability of the velocity model inferred by ray-tracing modelling. The interpretation of the velocity model is constrained by stratigraphic and sonic velocity logs from wells for oil exploration located close to the seismic line. Gravity data modelling allows to check the velocity model and to extend the structural interpretation in 3-D. In the shallow crust, up to a depth of 3–4 km, strong lateral variations of the modelled velocities are produced by the overlapping of thrust sheets formed by: (1) Cenozoic flyschoid cover and basinal successions that underlie the seismic profile with P-wave velocities in the 2.8–4.1 km/s range and thicknesses varying between 0.5 and 4.5 km; (2) Mesozoic basinal sequences with a velocity of 4.8 km/s and a depth of 1.5–2.1 km in the northern part of the profile; (3) Mesozoic limestones of the Western Carbonate Platform with a velocity of 6.0 km/s and a depth of 0.1–0.8 km in the southern part of the profile. At a greater depth, the model becomes more homogeneous. A continuous seismic interface 3.0–4.5 km deep with a velocity of 6.0 km/s is interpreted as the top of the Meso-Cenozoic Carbonate Multilayer of the Apulia Platform, characterized by an increase in seismic velocity from 6.2 to 6.6 km/s at depths of 6–7 km. A lower P-wave velocity (about 5.0 km/s) is hypothesized at depths ranging between 9.5 and 11 km. As inferred by commercial seismic lines and data from two deep wells located in the Apulia foreland and Bradano foredeep, this low-velocity layer can be related to Permo-Triassic clastic deposits drilled at the bottom of the Apulia Platform. Seismic data do not allow us to identify possible deeper seismic interfaces that could correspond to the top of the Paleozoic crystalline basement; this is probably due to the low-velocity layer at the bottom of the Carbonate Multilayer that reflects and attenuates a great part of the seismic energy. The joint interpretation of seismic refraction and well data, in accordance with gravity data, provides the first detailed P-wave velocity model of the upper crust of the northern sector of the Southern Apennines, which differs considerably from previous 1-D velocity models used to study the seismicity of the region, and reveals new information about the structure of the thrust belt.

Asymptotic ray theory in poro- and viscoelastic media

Asymptotic ray methods are applied to two kinds of dispersive and attenuating media: viscoelasticity and Biot’s model of poroelastic media. Asymptotic solutions of the first kind are based on the assumption that the relaxation time and the period of the signal are small with respect to the microinhomogeneity scale (short memory asymptotics). Depending on the ratio of viscous stress to the elastic stress a viscoelastic medium is either strongly or weakly attenuating. A short memory asymptotic solution in a strongly attenuating case involves complex rays. A general theory of short memory asymptotics in hereditary media including poroelastic media is constructed. For poroelastic media the behavior of the solution near the wavefront is examined by generalized wavefront asymptotics. The signal at the wavefront vanishes with all its derivatives and its peak is delayed with respect to the wavefront. Comparison with numerically exact solutions shows excellent agreement close to the wavefront.