Velocity analysis for transversely isotropic media

Summary Numerical and physical modelling studies have been used to investigate the accuracy and limitations of the small-offset P-wave NMO velocity estimation in transversely isotropic media. An expression for the P-wave small-offset NMO velocity for a single horizontal reflector in a vertically transversely isotropic layer has been previously obtained: V nmo =+ αδ 0 12 where αο and δ respectively represent the vertical P-wave velocity and the P-wave critical anisotropy at oblique angles. The above equation is presented as being exact and valid for any degree of anisotropy and has been used in estimating the small-off set NMO velocity in the offshore Zaire seismic data acquired by Unocal in West Africa. The accuracy of this expression is subject t o practical testing in this paper. A one-layer horizontal model was used to carry out the numerical modelling studies. A computer program which calculates Pwave NMO velocities in transversely isotropic media was used. This program uses the exact equations for the phase and ray velocities in generating velocity functions to be used in ray tracing; it is known to give accurate results. Using the publishe d elastic parameters for anisotropic rocks, the small-offset NMO velocities are computed and results are compared with values obtained using the above expression. The degree of anisotropy varied from very weak to strong and contrasting velocity function s (positive and negative δ) were used in numerical modelling simulations. Phenolite materials with contrasting velocity functions were used to simulate transversely isotropic media with vertical axis of symmetry in physical modelling tests. Experimental reflection data were collected and velocity analysis was conducted for small-offset tr aces, and the NMO velocity determined by the maximum stack response. The elastic parameters αο and δ of the Phenolite materials were recovered by an anisotropic inversion technique. These parameters were used to estimate the small-offset P-wave NMO velocity in Phenolite with the above expression and also with our computer program. The NMO velocities obtained from these two methods were compared to that obtained from velocity analysis of the experimental seismic data. These comparisons enabled the validity of the above expressio n to yield accurate P-wave NMO velocity, in the limit of small offset, to be tested. Numerical and physical modelling results indicate that the accuracy of the small-offset P-wave NMO equation depends on the nature and degree of the anisotropy encountered in a given area, and is not valid for all degrees of P-wave anisotropy. The sma lloffset NMO (stacking) velocity from horizontal reflectors is known to influence the effects of vertical transverse isotropy on normal moveout (NMO), dip moveout (DMO) corrections and time migration. Determination of the correct small-offset P-wave NMO velocity is essential for obtaining accurate seismic images, and subsequently correct interpretations.

Anisotropic P-wave velocity analysis and seismic imaging in onshore Kutch sedimentary basin of India

We consider a case where a 3D depth migration has been performed in the local angle domain (LAD) using rich-azimuth seismic data (e.g., conventional land surveys). The subsurface geologic model is characterized by considerable azimuthally anisotropic velocity variations. The background velocity field used for the migration can consist of azimuthally independent, e.g., vertical transverse isotropy, and/or azimuthally dependent (e.g., orthorhombic), velocity layers. The resulting 3D full-azimuth reflection angle gathers generated by the LAD migration represent in situ high-resolution amplitude preserved reflectivities associated with opening angles between incident and reflected slowness vectors in the specular directions. Residual moveouts (RMOs) automatically picked on these 3D image gathers along major horizons can indicate considerable residual periodic azimuthal variations. This situation is typical in depth imaging applied to unconventional shale plays, where the background velocity model doesn’t yet account for the aligned stress/fracture systems that exist in some of the target layers. We use the azimuthally dependent, phase-angle RMOs to update the interval parameters of the background model, accounting for the azimuthal anisotropy effect. Until now, this problem was mainly treated in the unmigrated time-offset domain, which is limited in describing the actual in situ changes of the velocity field with azimuths. The subsurface full-azimuth phase-angle domain RMOs provide better physical parameters to analyze the in situ azimuthal variations of the anisotropic media. Our method is grounded in a newly derived generalized Dix-based theory, where locally the background and updated models are assumed to be 1D anisotropic velocity models. At each lateral location, the orthorhombic axis x3 points in the vertical direction across all layers, but the azimuthal orientations of the orthorhombic layers change from layer to layer. An effective model for such a layered structure (background or updated) is represented by a single layer with a vertical time identical to that of the whole package, effective fast and slow normal moveout (NMO) velocities, and an effective azimuthal orientation of the slow NMO velocity. Our approach begins with computation of these effective parameters for the background model and conversion of the high-resolution RMOs into a dense set of updated, effective, azimuthally dependent NMO velocities, which are then converted into three effective parameters of the updated model. Next, we apply a generalized Dix-based inversion approach to estimate the local NMO parameters for each updated layer. Finally, we convert the local parameters into interval azimuthally varying anisotropic model parameters (e.g., TTI, orthorhombic, or tilted orthorhombic) within each layer. The 1D Dix-based approach presented in this work should not be considered an alternative to more accurate 3D global inversion approaches, such as global anisotropic tomography. However, the proposed method can be effectively used for moderately laterally varying models, and some of the principal physical rules derived for the 1D model can be further used to improve the formulation and geophysical constraints applied to 3D global inversion methods.

Despite the complexity of wave propagation in anisotropic media, reflection moveout on conventional common‐midpoint (CMP) spreads is usually well described by the normal‐moveout (NMO) velocity defined in the zero‐offset limit. In their recent work, Grechka and Tsvankin showed that the azimuthal variation of NMO velocity around a fixed CMP location generally has an elliptical form (i.e. plotting the NMO velocity in each azimuthal direction produces an ellipse) and is determined by the spatial derivatives of the slowness vector evaluated at the CMP location. This formalism is used here to develop exact solutions for the NMO velocity in anisotropic media of arbitrary symmetry.For the model of a single homogeneous layer above a dipping reflector, we obtain an explicit NMO expression valid for all pure modes and any orientation of the CMP line with respect to the reflector strike. The contribution of anisotropy to NMO velocity is contained in the slowness components of the zero‐offset ray (along with the derivatives of the vertical slowness with respect to the horizontal slownesses) — quantities that can be found in a straightforward way from the Christoffel equation. If the medium above a dipping reflector is horizontally stratified, the effective NMO velocity is determined through a Dix‐type average of the matrices responsible for the ‘interval’ NMO ellipses in the individual layers. This generalized Dix equation provides an analytic basis for moveout inversion in vertically inhomogeneous, arbitrarily anisotropic media. For models with a throughgoing vertical symmetry plane (i.e. if the dip plane of the reflector coincides with a symmetry plane of the overburden), the semi‐axes of the NMO ellipse are found by the more conventional rms averaging of the interval NMO velocities in the dip and strike directions.Modelling of normal moveout in general heterogeneous anisotropic media requires dynamic ray tracing of only one (zero‐offset) ray. Remarkably, the expressions for geometrical spreading along the zero‐offset ray contain all the components necessary to build the NMO ellipse. This method is orders of magnitude faster than multi‐azimuth, multi‐offset ray tracing and, therefore, can be used efficiently in traveltime inversion and in devising fast dip‐moveout (DMO) processing algorithms for anisotropic media. This technique becomes especially efficient if the model consists of homogeneous layers or blocks separated by smooth interfaces.The high accuracy of our NMO expressions is illustrated by comparison with ray‐traced reflection traveltimes in piecewise‐homogeneous, azimuthally anisotropic models. We also apply the generalized Dix equation to field data collected over a fractured reservoir and show that P‐wave moveout can be used to find the depth‐dependent fracture orientation and to evaluate the magnitude of azimuthal anisotropy.

The paper presents estimates of the S-wave velocity and the crack density at which fractured reservoirs begin to play an important role in oil exploration. Transverse isotropy with a horizontal axis of symmetry (HTI) is the simplest azimuthally anisotropic model used to describe fractured reservoirs that contain parallel vertical cracks. A double profile concept is used to develop an equation for the P-S wave normal-moveout (NMO) velocity. The azimuthal NMO velocities of the P- and P-S waves can then be used to estimate the velocities of the S-waves and Thomsen's coefficient, γ. For multilayered media, a recursive equation is developed for the NMO velocity in each layer. The numerical results indicate that the S-wave NMO velocity can be accurately estimated using the P- and P-S wave NMO velocities in HTI media. An important parameter of fracture systems that can be measured from seismic data is the crack density which can be estimated using the NMO velocities of the P- and S-waves from horizontal reflectors. Therefore, fractures can be completely characterized by the joint inversion of P-waves and converted P-S waves in HTI media.

Seismic velocity modeling is a crucial step in seismic processing that enables the use of velocity information from both seismic and wells to map the depth and thickness of subsurface layers interpreted from seismic images. The velocity can be obtained in the form of normal moveout (NMO) velocity or by an inversion (optimization) process such as in full-waveform inversion (FWI). These methods have several limitations. These limitations include enormous time consumption in the case of NMO due to manual and heavy human involvement in the picking. As an optimization problem, it incurs high cost and suffers from nonlinearity issues. Researchers have proposed various machine learning (ML) techniques including unsupervised, supervised, and semi-supervised learning methods to model the velocity more efficiently. The focus of the studies is mostly to automate the NMO velocity picking, improve the convergence in FWI, and apply FWI using ML directly from the data. In the purview of the digital transformation roadmap of the petroleum industry, this paper presents a chronologic review of these studies, appraises the progress made so far, and concludes with a set of recommendations to overcome the prevailing challenges through the implementation of more advanced ML methodologies. We hope that this work will benefit experts, young professionals, and ML enthusiasts to help push forward their research efforts to achieving complete automation of the NMO velocity and further enhancing the performance of ML applications used in the FWI framework.

Progress in seismic inversion and processing in anisotropic media depends on our ability to relate different seismic signatures to the anisotropic parameters. While the conventional notation (stiffness coefficients) is suitable for forward modeling, it is inconvenient in developing analytic insight into the influence of anisotropy on wave propagation. Here, a consistent description of P-wave signatures in transversely isotropic (TI) media with arbitrary strength of the anisotropy is given in terms of Thomsen notation. The influence of transverse isotropy on P-wave propagation is shown to be practically independent of the vertical S-wave velocity VS0, even in models with strong velocity variations. Therefore, the contribution of transverse isotropy to P-wave kinematic and dynamic signatures is controlled by just two anisotropic parameters, ε and δ, with the vertical velocity VP0 being a scaling coefficient in homogeneous models. The distortions of reflection moveouts and amplitudes are not necessarily correlated with the magnitude of velocity anisotropy. The influence of transverse isotropy on P-wave normal-moveout (NMO) velocity in a horizontally layered medium, on small-angle reflection coefficient, and on point-force radiation in the symmetry direction is entirely determined by the parameter δ. Another group of signatures of interest in reflection seisimology–the dip-dependence of NMO velocity, magnitude of nonhyperbolic moveout, time-migration impulse response, and the radiation pattern near vertical–is dependent on both anisotropic parameters (ε and δ) and is primarily governed by the difference between ε and δ. Since P-wave signatures are so sensitive to the value of ε − δ, application of the elliptical-anisotropy approximation (ε = δ) in P-wave processing may lead to significant errors. Many analytic expressions given in the paper remain valid in transversely isotropic media with a tilted symmetry axis. Moreover, the equation for NMO velocity from dipping reflectors, as well as the nonhyperbolic moveout equation, can be used in symmetry planes of any anisotropic media (e.g., orthorhombic).

The conventional normal‐moveout (NMO) and common‐midpoint (CMP) stacking process enhances reflections having a particular moveout velocity, while attenuating events (such as multiple reflections) having different moveout velocities. Unfortunately, this process also acts as a dip filter applied to the CMP stack. In other words, NMO and stacking enhances reflections having a particular slope in the CMP stack, while attenuating reflections having different slopes. NMO and stacking, like any dip filter, degrades lateral resolution. Fortunately, this dip‐filtering action can be suppressed by applying, in addition to NMO, a prestack process known variously as DEVILISH, prestack partial migration, and dip‐moveout. As the latter term implies, this process is a dip‐dependent moveout correction that enables reflections from both horizontal and dipping reflectors to be stacked with the same NMO velocity. Stated another way, NMO velocities estimated from dip‐moveout‐corrected seismograms are independent of the dips of subsurface reflectors. Dip‐moveout by Fourier transform is a method for performing dip‐moveout (DMO) correction in the frequency‐wavenumber domain. The implementation of this method, which resembles the implementation of a discrete Fourier transform, is quite different from and compares favorably with previously published finite‐difference DMO algorithms. DMO by Fourier transform, unlike DMO by finite‐differences, is accurate for all offsets and all dips, provided that velocity is constant. Because velocity is never constant, some accuracy is inevitably lost; but the application of DMO by Fourier transform to recorded seismograms demonstrates the ability of this process to enhance (1) the dip bandwidth of CMP stacks and (2) the accuracy of velocity estimates.

Surface-consistent residual statics correction for land seismic data does not account for the source - receiver offset. Consequently, it requires normal moveout (NMO) corrected gathers to bring raypaths close to the normal incidence. When the NMO velocity is inaccurate or unavailable, the estimated statics suffer. Therefore, multiple passes of NMO velocity picking and residual statics estimation become essential, which are efforts and time consuming. To avoid this, we utilize a rank-based solution that is capable of estimating non-surface consistent residual statics. The method is based on the rank property of frequency slices in the midpoint-offset domain, where ideal seismic data is of low-rank nature, while data with residual statics exhibits higher rank. Accordingly, we estimate the statics that lead to the desired low-rank signal via means of low-rank approximation and cross-correlation in an iterative and multi-rank-scale approach. Since we estimate non surface-consistent statics by accounting for the offset of each trace, it is no longer required to have NMO corrected gathers. Consequently, the method does not require windowing over a noise-free section containing primaries or windowing to avoid the NMO stretch effect, which are required by conventional residual statics correction. Numerical results on simulated and field data suggest that the method has the potential of replacing existing residual statics correction techniques.

ABSTRACTWe study the azimuthally dependent hyperbolic moveout approximation for small angles (or offsets) for quasi‐compressional, quasi‐shear, and converted waves in one‐dimensional multi‐layer orthorhombic media. The vertical orthorhombic axis is the same for all layers, but the azimuthal orientation of the horizontal orthorhombic axes at each layer may be different. By starting with the known equation for normal moveout velocity with respect to the surface‐offset azimuth and applying our derived relationship between the surface‐offset azimuth and phase‐velocity azimuth, we obtain the normal moveout velocity versus the phase‐velocity azimuth. As the surface offset/azimuth moveout dependence is required for analysing azimuthally dependent moveout parameters directly from time‐domain rich azimuth gathers, our phase angle/azimuth formulas are required for analysing azimuthally dependent residual moveout along the migrated local‐angle‐domain common image gathers. The angle and azimuth parameters of the local‐angle‐domain gathers represent the opening angle between the incidence and reflection slowness vectors and the azimuth of the phase velocity ψphs at the image points in the specular direction. Our derivation of the effective velocity parameters for a multi‐layer structure is based on the fact that, for a one‐dimensional model assumption, the horizontal slowness and the azimuth of the phase velocity ψphs remain constant along the entire ray (wave) path. We introduce a special set of auxiliary parameters that allow us to establish equivalent effective model parameters in a simple summation manner. We then transform this set of parameters into three widely used effective parameters: fast and slow normal moveout velocities and azimuth of the slow one. For completeness, we show that these three effective normal moveout velocity parameters can be equivalently obtained in both surface‐offset azimuth and phase‐velocity azimuth domains.

Velocity anisotropy should be taken into account when analyzing the amplitude variation with offset (AVO) response of gas sands encased in shales. The anisotropic effects on the AVO of gas sands in transversely isotropic (TI) media are reviewed. Reflection coefficients in TI media are computed using a planewave formula based on ray theory. We present results of modeling special cases of exploration interest having positive reflectivity, near‐zero reflectivity, and negative reflectivity. The AVO reflectivity in anisotropic media can be decomposed into two parts; one for isotropy and the other for anisotropy. Zero‐offset reflectivity and Poisson’s ratio contrast are the most significant parameters for the isotropic component while the δ difference (Δδ) between shale and gas sand is the most important factor for the anisotropic component. For typical values of Tl anisotropy in shale (positive δ and ε), both δ difference (Δδ) and ε difference (Δε) amplify AVO effects. For small angles of incidence, Δδ plays an important role in AVO while Δε dominates for large angles of incidence. For typical values of δ and ε, the effects of anisotropy in shale are: (1) a more rapid increase in AVO for Class 3 and Class 2 gas sands, (2) a more rapid decrease in AVO for Class 1 gas sands, and (3) a shift in the offset of polarity reversal for some Class 1 and Class 2 gas sands.

The basic requirements to recover plane layers of constant interval velocity, arbitrary dip and strike from common depth point (CDP) recordings are the following four quantities related to the primary event of each reflector at the common midpoint of a CDP profile: a) Two‐way normal time b) Normal moveout velocity within one arbitrary CDP profile c) Time slope of normally reflected rays within the profile d) Time slope of normally reflected rays in some other direction. The solution of the inverse problem is obtained directly. The moveout velocity is expressed in terms of seismic parameters along the normal incidence path in three dimensions and the direction of the profile within the free surface. A formula connecting dip and strike of the emerging normal ray with the measured time gradients is given and discussed. The method includes, as a special case, the Dix formulas for plane parallel layers.

Surface reflection tomography is an inversion method that attempts to determine simultaneously the subsurface elastic parameters and reflector depths from pre-stack seismic data. To date, this technique has been applied to acoustic and isotropic elastic media. In this paper, we study surface reflection tomography in the anisotropic elastic case. Anisotropic media are described by many more elastic parameters than are needed to describe isotropic media; these additional parameters are useful for obtaining rock properties such as crack density, pore shape and fracture strike that cannot be found with isotropic methods. Similar to isotropic surface reflection tomography, many features of an anisotropic model are poorly resolved due to limited ray path coverage; the vertical smearing of isotropic tomography occurs in each elastic parameter in the anisotropic problem. Unlike the isotropic problem, however, there is additional indeterminacy in the solution of the anisotropic tomography problem because of ambiguity amongst the several elastic parameters needed to describe anisotropic media. Also unlike the isotropic problem, the reflector depths participate very strongly in this ambiguity. We investigate the nature of this indeterminacy by studying the null space for linearized tomography, that is, the class of model perturbations of a background medium which, to first order, cause no perturbation at all in the surface reflection traveltimes. Such model perturbations cannot be determined from the traveltime perturbations; describing these null space model perturbations gives insight into the indeterminacy in the anisotropic problem. Complementary to computational approaches towards identifying the null space for discrete formulations of tomography, we study a continuum formulation. Our first set of results concerns the ‘elastic’ null space: perturbations of the elastic parameters (with zero depth perturbation) which cause no perturbation in the traveltimes. The results here are very similar to previous results for cross-well transmission tomography. As expected, the ‘elastic’ null space is larger than the isotropic null space due to the ambiguity amongst the elastic parameters. We identify three categories of model perturbations in the ‘elastic’ null space. The first category consists of model perturbations for which the perturbation in each of the individual elastic parameters is itself either in the isotropic null space, or in a closely related set of model perturbations which we call the odd isotropic null space. Elements in the second category are anisotropic versions of the most well-known isotropic null space elements: perturbations which are polynomials in the horizontal variable with coefficients which are functions of the depth variable satisfying certain linear integral constraints; unlike the isotropic problem, the integral constraints in the anisotropic problem couple together the several elastic parameters. The third category consists of model perturbations satisfying zero boundary conditions at the surface and at the reflector for which a specific linear combination of integrals and derivatives of the several elastic parameters is in the isotropic null space. In particular, there are model perturbations in this third category which represent anomalies that are completely contained in the interior of the model and yet are in the ‘elastic’ null space; this behaviour is different from the isotropic problem. These categories are sufficient to describe the ‘elastic’ null space completely. We demonstrate that every model perturbation in the ‘elastic’ null space is the sum of an element in the first category (indicating an indeterminacy of the same nature as in the isotropic problem in each of the elastic parameters separately) and an element in the third category (indicating an ambiguity amongst the parameters). The second category gives a rich family of examples of sums of null space elements in the first and third categories, and thereby gives a sense of just how large the ‘elastic’ null space is. Moreover, we show that the traveltime perturbations caused by an elastic perturbation determine only a small number of features of the elastic perturbation which distinguish between the several elastic parameters. We identify these features pre-cisely: they are functions of the horizontal variable representing vertical averages of combinations of the elastic parameters and their derivatives. Elastic parameters influencing the vertical velocity appear more prominently in these features than those influencing the horizontal velocity, and these features are closely related to the zero-offset traveltimes and the normal moveout velocities. Our second set of results concerns the so-called velocity–depth ambiguity, more precisely here the elastic parameter versus the depth parameter ambiguity. We demonstrate that for any smooth depth perturbation, there is a perturbation in the elastic parameters for which the combined elastic and depth perturbations cause zero traveltime perturbations; moreover, this can be accomplished by an elliptically anisotropic elastic perturbation with a vertical axis of symmetry, perturbing the vertical velocity alone. This implies that the velocity–depth ambiguity cannot be resolved here as it can be in the isotropic case. We show an example of an anticlinal structure embedded in a homogeneous background whose reflection traveltimes exactly match those from a model with elliptically anisotropic elastic perturbations overlying a flat reflector. We also draw a number of conclusions from these results for the special case of transversely isotropic media with a vertical symmetry axis.

Nonhyperbolic (long‐spread) moveout provides essential information for a number of seismic inversion/processing applications, particularly for parameter estimation in anisotropic media. Here, we present an analytic expression for the quartic moveout coefficient A4that controls the magnitude of nonhyperbolic moveout of pure (nonconverted) modes. Our result takes into account reflection‐point dispersal on irregular interfaces and is valid for arbitrarily anisotropic, heterogeneous media. All quantities needed to compute A4can be evaluated during the tracing of the zero‐offset ray, so long‐spread moveout can be modeled without time‐consuming multioffset, multiazimuth ray tracing.The general equation for the quartic coefficient is then used to study azimuthally varying nonhyperbolic moveout of P‐waves in a dipping transversely isotropic (TI) layer with an arbitrary tilt ν of the symmetry axis. Assuming that the symmetry axis is confined to the dip plane, we employed the weak‐anisotropy approximation to analyze the dependence of A4on the anisotropic parameters. The linearized expression for A4is proportional to the anellipticity coefficient η ≈ ε − δ and does not depend on the individual values of the Thomsen parameters. Typically, the magnitude of nonhyperbolic moveout in tilted TI media above a dipping reflector is highest near the reflector strike, whereas deviations from hyperbolic moveout on the dip line are substantial only for mild dips.The azimuthal variation of the quartic coefficient is governed by the tilt ν and reflector dip φ and has a much more complicated character than the NMO–velocity ellipse. For example, if the symmetry axis is vertical (VTI media, ν = 0) and the dip φ < 30°, A4goes to zero on two lines with different azimuths where it changes sign. If the symmetry axis is orthogonal to the reflector (this model is typical for thrust‐and‐fold belts), the strike‐line quartic coefficient is defined by the well‐known expression for a horizontal VTI layer (i.e., it is independent of dip), while the dip‐line A4is proportional to cos4φ and rapidly decreases with dip. The high sensitivity of the quartic moveout coefficient to the parameter η and the tilt of the symmetry axis can be exploited in the inversion of wide‐azimuth, long‐spread P‐wave data for the parameters of TI media.

It is necessary to account for anisotropy in full waveform inversion (FWI) of wide-azimuth and wide-aperture seismic data in most geologic environments, for correct depth positioning of reflectors, and for reliable estimations of wave speeds as a function of the direction of propagation. In this framework, choosing a suitable anisotropic subsurface parameterization is a central issue in monoparameter and multiparameter FWI. This is because this parameterization defines the influence of each physical parameter class on the data as a function of the scattering angle, and hence the resolution of the parameter reconstruction, and on the potential trade-off between different parameter classes. We apply monoparameter and multiparameter frequency-domain acoustic vertical transverse isotropic FWI to synthetic and real wide-aperture data, representative of the Valhall oil field. We first show that reliable monoparameter FWI can be performed to build a high-resolution velocity model (for the vertical, the horizontal, or normal move-out velocity), provided that the background models of two Thomsen parameters describe the long wavelengths of the subsurface sufficiently accurately. Alternatively, we show the feasibility of the joint reconstruction of two wave speeds (e.g., the vertical and horizontal wave speeds) with limited trade-off effects, while Thomsen parameter [Formula: see text] is kept fixed during the inversion. The influence of the wave speeds on the data for a limited range of scattering angles when combined each other can, however, significantly hamper the resolution with which the two wave speeds are imaged. These conclusions inferred from the application to the real data are fully consistent with those inferred from the theoretical parameterization analysis of acoustic vertical transverse isotropic FWI performed in the companion report.