Search

SUMMARY The inverse problem of finding the slowness vector from a known ray direction in general anisotropic elastic media is a challenging task, needed in many wave/ray-based methods, in particular, solving two-point ray bending problems. The conventional resolving equation set for general (triclinic) anisotropy consists of two fifth-degree polynomials and a sixth-degree polynomial, resulting in a single physical solution for quasi-compressional (qP) waves and up to 18 physical solutions for quasi-shear waves (qS). For polar anisotropy (transverse isotropy with a tilted symmetry axis), the resolving equations are formulated for the slowness vectors of the coupled qP and qSV waves (quasi-shear waves polarized in the axial symmetry plane), and independently for the decoupled pure shear waves polarized in the normal (to the axis) isotropic plane (SH). The novelty of our approach is the introduction of the geometric constraint that holds for any wave mode in polar anisotropic media: The three vectors—the slowness, ray velocity and medium symmetry axis—are coplanar. Thus, the slowness vector (to be found) can be presented as a linear combination of two unit-length vectors: the polar axis and the ray velocity directions, with two unknown scalar coefficients. The axial energy propagation is considered as a limit case. The problem is formulated as a set of two polynomial equations describing: (i) the collinearity of the slowness-related Hamiltonian gradient and the ray velocity direction (third-order polynomial equation) and (ii) the vanishing Hamiltonian (fourth-order polynomial equation). Such a system has up to twelve real and complex-conjugate solutions, which appear in pairs of the opposite slowness directions. The common additional constraint, that the angle between the slowness and ray directions does not exceed ${90^{\rm{o}}}$, cuts off one half of the solutions. We rearrange the two bivariate polynomial equations and the above-mentioned constraint as a single univariate polynomial equation of degree six for qP and qSV waves, where the unknown parameter is the phase angle between the slowness vector and the medium symmetry axis. The slowness magnitude is then computed from the quadratic Christoffel equation, with a clear separation of compressional and shear roots. The final set of slowness solutions consists of a unique real solution for qP wave and one or three real solutions for qSV (due to possible triplications). The indication for a qSV triplication is a negative discriminant of the sixth-order polynomial equation, and this discriminant is computed and analysed directly in the ray-angle domain. The roots of the governing univariate sixth-order polynomial are computed as eigenvalues of its companion matrix. The slowness of the SH wave is obtained from a separate equation with a unique analytic solution. We first present the resolving equation using the stiffness components, and then show its equivalent forms with the well-known parametrizations: Thomsen, Alkhalifah and ‘weak-anisotropy’. For the Thomsen and Alkhalifah forms, we also consider the (essentially simplified) acoustic approximation for qP waves governed by the quartic polynomials. The proposed method is coordinate-free and can be applied directly in the global Cartesian frame. Numerical examples demonstrate the advantages of the method.

This paper is the second in a sequel of two papers and dedicated to the computation of paraxial rays and dynamic characteristics along the stationary rays obtained in the first paper. We start by formulating the linear, second‐order, Jacobi dynamic ray tracing equation. We then apply a similar finite‐element solver, as used for the kinematic ray tracing, to compute the dynamic characteristics between the source and any point along the ray. The dynamic characteristics in our study include the relative geometric spreading and the phase correction due to caustics (i.e. the amplitude and the phase of the asymptotic form of the Green's function for waves propagating in 3D heterogeneous general anisotropic elastic media). The basic solution of the Jacobi equation is a shift vector of a paraxial ray in the plane normal to the ray direction at each point along the central ray. A general paraxial ray is defined by a linear combination of up to four basic vector solutions, each corresponds to specific initial conditions related to the ray coordinates at the source. We define the four basic solutions with two pairs of initial condition sets: point–source and plane‐wave. For the proposed point–source ray coordinates and initial conditions, we derive the ray Jacobian and relate it to the relative geometric spreading for general anisotropy. Finally, we introduce a new dynamic parameter, similar to the endpoint complexity factor, presented in the first paper, used to define the measure of complexity of the propagated wave/ray phenomena. The new weighted propagation complexity accounts for the normalized relative geometric spreading not only at the receiver point, but along the whole stationary ray path. We propose a criterion based on this parameter as a qualifying factor associated with the given ray solution. To demonstrate the implementation of the proposed method, we use several isotropic and anisotropic benchmark models. For all the examples, we first compute the stationary ray paths, and then compute the geometric spreading and analyse these trajectories for possible caustics. Our primary aim is to emphasize the advantages, transparency and simplicity of the proposed approach.

Kinematical characteristics of reflected waves in anisotropic elastic media play an important role in the seismic imaging workflow. Considering compressional and converted waves, we derive new, azimuthally dependent, slowness‐domain approximations for the kinematical characteristics of reflected waves (radial and transverse offsets, intercept time and traveltime) for layered orthorhombic media with varying azimuth of the vertical symmetry planes. The proposed method can be considered an extension of the well‐known ‘generalized moveout approximation’ in the slowness domain, from azimuthally isotropic to azimuthally anisotropic models. For each slowness azimuth, the approximations hold for a wide angle range, combining power series coefficients in the vicinity of both the normal‐incidence ray and an additional wide‐angle ray. We consider two cases for the wide‐angle ray: a ‘critical slowness match’ and a ‘pre‐critical slowness match’ studied in Parts I and II of this work, respectively. For the critical slowness match, the approximations are valid within the entire slowness range, up to the critical slowness. For the ‘pre‐critical slowness match’, the approximations are valid only within the bounded slowness range; however, the accuracy within the defined range is higher. The critical slowness match is particularly effective when the subsurface model includes a dominant high‐velocity layer where, for nearly critical slowness values, the propagation in this layer is almost horizontal. Comparing the approximated kinematical characteristics with those computed by numerical ray tracing, we demonstrate high accuracy.

Considering all types of pure-mode and converted waves, we derive the azimuthally dependent, fourth-order normal moveout (NMO) velocity functions, and hence the corresponding effective anellipticity functions, for horizontally layered orthorhombic media. We emphasize that this paper does not suggest a new nonhyperbolic traveltime approximation; rather, it provides exact expressions of the NMO series coefficients, computed for normal-incidence rays, which can then be further used within known azimuthally dependent traveltime approximations for short to moderate offsets. We do not assume weak anisotropy or acoustic approximation for P-waves. At each layer, the elastic parameters, thickness, and azimuth of the orthorhombic vertical symmetry planes are considered to be different. We distinguish between two different azimuths: slowness azimuth (part 1 of this paper) and offset azimuth (part 2 of this paper). In part 1, the slowness-azimuth domain NMO is approximated as a series of either infinitesimal horizontal slowness (slowness-azimuth/slowness domain) or infinitesimal offsets (slowness-azimuth/offset domain). Similarly, in part 2, we distinguish between two offset-azimuth domains: offset-azimuth/slowness and offset-azimuth/offset. Note that the azimuthally dependent NMO velocity functions of each of the four cases are different. The validity of the method is tested by introducing our derived azimuthally dependent, fourth-order effective anellipticity, into the well-known azimuthally dependent, asymptotic nonhyperbolic traveltime approximation, in which we compare the traveltime approximation versus exact numerical ray tracing for short to moderate offsets. It is clearly shown that for these types of azimuthally anisotropic layered models, the fourth-order terms are essential even for relatively small horizontal-slowness values or short offsets.

Based on the theory derived in part 1, in which we obtained the azimuthally dependent fourth-order normal-moveout (NMO) velocity functions for layered orthorhombic media in the slowness-azimuth/slowness and the slowness-azimuth/offset domains, in part 2, we extend the theory to the offset-azimuth/slowness and offset-azimuth/offset domains. We reemphasize that this paper does not suggest a new nonhyperbolic traveltime approximation; rather, it provides exact expressions of the NMO series coefficients, computed for normal-incidence rays, which can then be further used within known azimuthally dependent traveltime approximations for short to moderate offsets. The same type of models as in part 1 are considered, in which the layers share a common horizontal plane of symmetry, but the azimuths of their vertical symmetry planes are different. The same eight local (single-layer) and global (overburden multilayer) effective parameters are used. In addition, we have developed an alternative set of global effective parameters in which the “anisotropic” effective parameters are normalized, classified into two groups: two “azimuthally isotropic” parameters and six “azimuthally anisotropic” parameters. These parameters have a clearer physical interpretation and they are suitable for inversion purposes because they can be controlled and constrained. Next, we propose a special case, referred to as “weak azimuthal anisotropy,” in which only the azimuthally anisotropic effective parameters are assumed to be weak. The resulting NMO velocity functions are considerably simplified, reduced to the form of the slowness-azimuth/slowness formula. We verify the correctness of our method by applying it to a multilayer orthorhombic medium with strong anisotropy. We introduce our derived, fourth-order slowness-azimuth/offset domain NMO velocity function into the well-known nonhyperbolic asymptotic traveltime approximation, and we compare the approximate traveltimes with exact traveltimes obtained by two-point ray tracing. The comparison shows an accurate match up to moderate offsets. Although the accuracy with the weak azimuthal anisotropic formula is inferior, it can still be considered reasonable for practical use.

Asymptotically bounded velocity profiles describe the vertical velocity variations in compacted sediments in a more realistic way than unbounded velocity models, and allow presenting the subsurface by a smaller number of thicker layers. The first and the simplest asymptotically bounded model is the Hyperbolic velocity profile proposed by Muscatin 1937, and our paper is an extension of this early study. The Hyperbolic model has an advantage over other bounded models: The velocity increases with depth and approaches the limiting value with a more smooth and gradual rate. We derive the time-depth relationships, forward and backward transforms between the instantaneous velocity profile and the effective models (average, RMS and fourth order average velocities), study the trajectories for pre-critical and post-critical curved rays and derive the equations for traveltime, lateral propagation and arc length. We compare the ray paths obtained with the Hyperbolic model and with the other bounded velocity profiles.

We use residual moveouts measured along continuous full azimuth reflection angle gathers, in order to obtain effective horizontal transversely isotropic model parameters. The angle gathers are generated through a special angle domain imaging system, for a wide range of reflection angles and full range of phase velocity azimuths. The estimation of the effective model parameters is performed in two stages. First, the background horizontal transversely isotropic (HTI)/vertical transversely isotropic (VTI) layered model is used, along with the values of reflection angles, for converting the measured residual moveouts (or traveltime errors) into azimuthally dependent normal moveout (NMO) velocities. Then we apply a digital Fourier transform to convert the NMO velocities into azimuthal wavenumber domain, in order to obtain the effective HTI model parameters: vertical time, vertical compression velocity, Thomsen parameter delta and the azimuth of the medium axis of symmetry. The method also provides a reliability criterion of the HTI assumption. The criterion shows whether the medium possesses the HTI type of symmetry, or whether the azimuthal dependence of the residual traveltime indicates to a more complex azimuthal anisotropy. The effective model used in this approach is defined for a 1D structure with a set of HTI, VTI and isotropic layers (with at least one HTI layer). We describe and analyse the reduction of a multi-layer structure into an equivalent effective HTI model. The equivalent model yields the same NMO velocity and the same offset azimuth on the Earth's surface as the original layered structure, for any azimuth of the phase velocity. The effective model approximates the kinematics of an HTI/VTI layered structure using only a few parameters. Under the hyperbolic approximation, the proposed effective model is exact.

Summary The simplest assumption on an anisotropic medium with azimuthally dependent traveltime is an HTI. We use residual moveouts measured along continuous full azimuth reflection angle gathers, in order to obtain effective HTI model parameters. The angle gathers are generated through a special angle domain imaging system, for a wide range of reflection angles and full range of phase velocity azimuths. The estimation of the effective model parameters is performed in two stages. First, the background HTI/VTI layered model is used, along with the values of reflection angles, for converting the measured residual moveouts (or traveltime errors) into azimuthally dependent NMO velocities. Then, we apply a Fourier transform to convert the NMO velocities into azimuthal wavenumber domain, in order to obtain the effective HTI model parameters: vertical time, vertical compression velocity, Thomsen parameter delta, and the azimuth of the medium axis of symmetry. The method also provides a reliability criterion of the HTI assumption. The effective model used in this approach is defined for a 1D model with a set of HTI, VTI, and isotropic layers. We describe and analyze the reduction of a multi-layer structure into an equivalent effective HTI model. The equivalent model yields the same NMO velocity and the same offset azimuth on the earth surface as the original layered structure, for any azimuth of the phase velocity. The effective model approximates the kinematics of an HTI/VTI layered structure using only a few parameters. Under the hyperbolic approximation, the proposed effective model is exact.

The ray-tracing algorithm presented in this paper is based on formulae derived for the common reflecting element (CRE) stacking method. A 2D, smooth, laterally-varying media is assumed where offset rays and traveltimes are evaluated from normal-incidence (central) rays. The method uses a second-order asymmetrical approximation for rays and an additional oblique spherical approximation of the central wavefronts for calculating offset traveltimes. In order to solve the two-point ray-tracing problem for the common midpoint (CMP) configuration of source-receiver pairs located symmetrically around the CMP stations, the central rays are perturbed to satisfy the above-mentioned asymmetrical distribution. Although the accuracy of the calculations is limited for far offsets, it is still good for distances of the order of the reflecting depths. Since only a few normal-incidence rays are traced through the medium, the method is very fast and is found to be most attractive for iterative inversions in macromodel estimation.

Summary Curved Rays Tomography updates background anisotropy velocity parameters in the time-migrated domain. Th e tomography uses image gathers generated by Anisotropy Curved Rays Kirchhoff Time Migration. A locally varying 1D Vertical Transverse Isotropy (VTI) model is assu med. The background anisotropy parameters are the instantaneous (interval) vertical compression veloc ity V and the two Thomsen anisotropy parameters δ and e . Interval velocity (or alternatively δ ) is updated from short offsets reflection events, while e is updated from the available long offset data. Two complementary approaches are presented in this study: local and global. In t he local approach, the medium parameters are updated from top down, layer by layer, one parameter at a time. The residual anisotropy parameters, that best fit the residual m oveout curves, are picked. The residual moveout includes overburden and current layer components. In the glo bal approach, all parameters are inverted simultaneousl y. Due to a large number of offsets, the problem becomes o verdefined, and we solve it by a constrained least-squ ares minimization. The cost function accounts for data a nd model variances, which reflect the reliability of t he data and control parameter variations, respectively. The updated parameters are constrained to a feasible range.

Paraxial rays are used for 3D two point ray tracing in complex structures. In this work we use dynamic ray tracing which is formulated in cartessian coordinates. This formulation allows for more flexibility in perturbing the initial slowness and takeoff point of the rays, as compared to the more common dynamic ray tracing in ray centered coordinates. Furthermore, the paraxial rays are adjusted iteratively to conform to the surface and interface topography, and to have slowness vectors with the correct length. This extends considerably the range of validity of the paraxial approximation. We describe implementation of different types of two point ray tracing which include CRP ray tracing, normal incidence ray tracing and CMP ray tracing.