3D P‐wave traveltime computation in transversely isotropic media using layer‐by‐layer wavefront marching

Subsurface rocks (e.g. shale) may induce seismic anisotropy, such as transverse isotropy. Traveltime computation is an essential component of depth imaging and tomography in transversely isotropic media. It is natural to compute the traveltime using the wavefront marching method. However, tracking the 3D wavefront is expensive, especially in anisotropic media. Besides, the wavefront marching method usually computes the traveltime using the eikonal equation. However, the anisotropic eikonal equation is highly non‐linear and it is challenging to solve. To address these issues, we present a layer‐by‐layer wavefront marching method to compute the P‐wave traveltime in 3D transversely isotropic media. To simplify the wavefront tracking, it uses the traveltime of the previous depth as the boundary condition to compute that of the next depth based on the wavefront marching. A strategy of traveltime computation is designed to guarantee the causality of wave propagation. To avoid solving the non‐linear eikonal equation, it updates traveltime along the expanding wavefront by Fermat's principle. To compute the traveltime using Fermat's principle, an approximate group velocity with high accuracy in transversely isotropic media is adopted to describe the ray propagation. Numerical examples on 3D vertical transverse isotropy and tilted transverse isotropy models show that the proposed method computes the traveltime with high accuracy. It can find applications in modelling and depth migration.     

2D inversion of refraction traveltime curves using homogeneous functions

A method using simple inversion of refraction traveltimes for the determination of 2D velocity and interface structure is presented. The method is applicable to data obtained from engineering seismics and from deep seismic investigations. The advantage of simple inversion, as opposed to ray‐tracing methods, is that it enables direct calculation of a 2D velocity distribution, including information about interfaces, thus eliminating the calculation of seismic rays at every step of the iteration process. The inversion method is based on a local approximation of the real velocity cross‐section by homogeneous functions of two coordinates. Homogeneous functions are very useful for the approximation of real geological media. Homogeneous velocity functions can include straight‐line seismic boundaries. The contour lines of homogeneous functions are arbitrary curves that are similar to one another. The traveltime curves recorded at the surface of media with homogeneous velocity functions are also similar to one another. This is true for both refraction and reflection traveltime curves. For two reverse traveltime curves, non‐linear transformations exist which continuously convert the direct traveltime curve to the reverse one and vice versa. This fact has enabled us to develop an automatic procedure for the identification of waves refracted at different seismic boundaries using reverse traveltime curves. Homogeneous functions of two coordinates can describe media where the velocity depends significantly on two coordinates. However, the rays and the traveltime fields corresponding to these velocity functions can be transformed to those for media where the velocity depends on one coordinate. The 2D inverse kinematic problem, i.e. the computation of an approximate homogeneous velocity function using the data from two reverse traveltime curves of the refracted first arrival, is thus resolved. Since the solution algorithm is stable, in the case of complex shooting geometry, the common‐velocity cross‐section can be constructed by applying a local approximation. This method enables the reconstruction of practically any arbitrary velocity function of two coordinates. The computer program, known as godograf , which is based on this theory, is a universal program for the interpretation of any system of refraction traveltime curves for any refraction method for both shallow and deep seismic studies of crust and mantle. Examples using synthetic data demonstrate the accuracy of the algorithm and its sensitivity to realistic noise levels. Inversions of the refraction traveltimes from the Salair ore deposit, the Moscow region and the Kamchatka volcano seismic profiles illustrate the methodology, practical considerations and capability of seismic imaging with the inversion method.     

VTI介质P波非双曲时差分析

具有垂直对称轴的横向各向同性介质模型(VTI)是目前各向异性理论研究和多波多分量地震资料叠前成像处理中最常用的一种各向异性模型.VTI介质中反射 P波时距曲线一般不再是双曲线.基于不同的相速度近似公式会得到不同的时距关系式.文中对几种典型的非双曲时距曲线与射线追踪得到的准确时距曲线在不同各向异性强度下进行了对比，结果表明Muir等和Stovas等提出的非双曲时距公式由于过高地考虑了横波垂直速度的影响与精确的时距曲线有很大偏差；Tsvankin等提出的弱各向异性非双曲时距公式在ε-δ<0时误差增大；Alkhalifah等提出的非双曲时距公式在大炮检距任意各向异性强度下都具有较高的精度，适于在实际资料处理中应用.     

横向各向同性介质中地震波走时模拟

横向各向同性介质是地球内部广泛分布的一种各向异性介质.针对这种介质,我们对各向同性介质的最小走时树走时模拟方法进行了推广,推广后的方法可适用于非均匀、对称轴任意倾斜的横向各向同性介质模型.为保证计算效率,最小走时树的构建采用了一种子波传播区域随地震波传播动态变化的改进算法.对于弱各向异性介质,我们使用了一种新的地震波群速度近似表示方法,该方法基于用射线角近似表示相角的思想,对3种地震波（qP, qSV和qSH）均有较好的精度.应用本文地震波走时模拟方法对均匀介质、横向非均匀介质模型进行了计算,并将后者结果与弹性波方程有限元方法的模拟结果进行了对比,结果表明两者符合得很好.本文方法可用于横向各向同性介质的深度偏移及地震层析成像的深入研究.     

Post‐stack diffraction imaging in vertical transverse isotropy media using non‐hyperbolic moveout approximations

Diffractions carry valuable information about local discontinuities and small‐scale objects in the subsurface. They are still not commonly used in the process of geological interpretation. Many diffraction imaging techniques have been developed and applied for isotropic media, whereas relatively few techniques have been developed for anisotropic media. Ignoring anisotropy can result in low‐resolution images with wrongly positioned or spurious diffractors. In this article, we suggest taking anisotropy into account in two‐dimensional post‐stack domain by considering P‐wave non‐hyperbolic diffraction traveltime approximations for vertical transverse isotropy media, previously developed for reflection seismology. The accuracy of the final images is directly connected to the accuracy of the diffraction traveltime approximations. We quantified the accuracy of six different approximations, including hyperbolic moveout approximation, by the application of a post‐stack diffraction imaging technique on two‐dimensional synthetic data examples.     

A comparison of cross-hole electrical and seismic data in fractured rock

Cross‐hole anisotropic electrical and seismic tomograms of fractured metamorphic rock have been obtained at a test site where extensive hydrological data were available. A strong correlation between electrical resistivity anisotropy and seismic compressional‐wave velocity anisotropy has been observed. Analysis of core samples from the site reveal that the shale‐rich rocks have fabric‐related average velocity anisotropy of between 10% and 30%. The cross‐hole seismic data are consistent with these values, indicating that observed anisotropy might be principally due to the inherent rock fabric rather than to the aligned sets of open fractures. One region with velocity anisotropy greater than 30% has been modelled as aligned open fractures within an anisotropic rock matrix and this model is consistent with available fracture density and hydraulic transmissivity data from the boreholes and the cross‐hole resistivity tomography data. However, in general the study highlights the uncertainties that can arise, due to the relative influence of rock fabric and fluid‐filled fractures, when using geophysical techniques for hydrological investigations.     

Depth migration anisotropy analysis in the time domain

In areas of complex geology such as the Canadian Foothills, the effects of anisotropy are apparent in seismic data and estimation of anisotropic parameters for use in seismic imaging is not a trivial task. Here we explore the applicability of common‐focus point (CFP)‐based velocity analysis to estimate anisotropic parameters for the variably tilted shale thrust sheet in the Canadian Foothills model. To avoid the inherent velocity‐depth ambiguity, we assume that the elastic properties of thrust‐sheet with respect to transverse isotropy symmetry axis are homogeneous, the reflector below the thrust‐sheet is flat, and that the anisotropy is weak. In our CFP approach to velocity analysis, for a poorly imaged reflection point, a traveltime residual is obtained as the time difference between the focusing operator for an assumed subsurface velocity model and the corresponding CFP response obtained from the reflection data. We assume that this residual is due to unknown values for anisotropy, and we perform an iterative linear inversion to obtain new model parameters that minimize the residuals. Migration of the data using parameters obtained from our inversion results in a correctly positioned and better focused reflector below the thrust sheet. For traveltime computation we use a brute force mapping scheme that takes into account weakly tilted transverse isotropy media. For inversion, the problem is set up as a generalized Newton's equation where traveltime error (differential time shift) is linearly dependent on the parameter updates. The iterative updates of parameters are obtained by a least‐squares solution of Newton's equations. The significance of this work lies in its applicability to areas where transverse isotropy layers are heterogeneous laterally, and where transverse isotropy layers are overlain by complex structures that preclude a moveout curve fitting.     

Pseudo‐spectral method using rotated staggered grid for elastic wave propagation in 3D arbitrary anisotropic media

Staggering grid is a very effective way to reduce the Nyquist errors and to suppress the non‐causal ringing artefacts in the pseudo‐spectral solution of first‐order elastic wave equations. However, the straightforward use of a staggered‐grid pseudo‐spectral method is problematic for simulating wave propagation when the anisotropy level is greater than orthorhombic or when the anisotropic symmetries are not aligned with the computational grids. Inspired by the idea of rotated staggered‐grid finite‐difference method, we propose a modified pseudo‐spectral method for wave propagation in arbitrary anisotropic media. Compared with an existing remedy of staggered‐grid pseudo‐spectral method based on stiffness matrix decomposition and a possible alternative using the Lebedev grids, the rotated staggered‐grid‐based pseudo‐spectral method possesses the best balance between the mitigation of artefacts and efficiency. A 2D example on a transversely isotropic model with tilted symmetry axis verifies its effectiveness to suppress the ringing artefacts. Two 3D examples of increasing anisotropy levels demonstrate that the rotated staggered‐grid‐based pseudo‐spectral method can successfully simulate complex wavefields in such anisotropic formations.     

First‐arrival traveltime tomography with modified total‐variation regularization

First‐arrival traveltime tomography is a robust tool for near‐surface velocity estimation. A common approach to stabilizing the ill‐posed inverse problem is to apply Tikhonov regularization to the inversion. However, the Tikhonov regularization method recovers smooth local structures while blurring the sharp features in the model solution. We present a first‐arrival traveltime tomography method with modified total‐variation regularization to preserve sharp velocity contrasts and improve the accuracy of velocity inversion. To solve the minimization problem of the new traveltime tomography method, we decouple the original optimization problem into the two following subproblems: a standard traveltime tomography problem with the traditional Tikhonov regularization and a L₂ total‐variation problem. We apply the conjugate gradient method and split‐Bregman iterative method to solve these two subproblems, respectively. Our synthetic examples show that the new method produces higher resolution models than the conventional traveltime tomography with Tikhonov regularization, and creates less artefacts than the total variation regularization method for the models with sharp interfaces. For the field data, pre‐stack time migration sections show that the modified total‐variation traveltime tomography produces a near‐surface velocity model, which makes statics corrections more accurate.     

A practical review of geostatistical processing applied to geophysical data: methods and applications

Nowadays, geostatistics is commonly applied for numerous gridding or modelling tasks. However, it is still under used and unknown for classical geophysical applications. This paper highlights the main geostatistical methods relevant for geophysical issues, for instance to improve the quality of seismic data such as velocity cubes or interpreted horizons. These methods are then illustrated through four examples. The first example, based on a gravity survey presents how a geostatistical interpolation can be used to filter out a global trend, in order to better define real anomalies. In the second case study, dedicated to refraction surveying, geostatistical filtering is used to filter out acquisition artefacts and identify the main geological structures. The third one is an example of porosity being integrated geostatistically with a seismic acoustic impedance map. The last example deals with geostatistical time to depth conversion; the interest of performing geostatistical simulations is finally discussed.     

各向异性TI介质qP反射波走时层析成像

地震走时层析成像是反演地层各向异性参数分布的有效方法,但是关于地震各向异性介质走时层析成像的研究并不多,其技术远远没有达到成熟的阶段.在野外数据采集时,地表反射波观测方式相对井间和垂直地震剖面观测方式的成本更低,利用qP反射波走时反演各向异性参数具有更加广泛的实用价值.本文实现的TI介质地震走时层析成像方法结合了TI介质反射波射线追踪算法、走时扰动方程和非线性共轭梯度算法,它可以对任意强度的TI介质模型进行反演,文中尝试利用qP反射波走时重建TI介质模型的参数图像.利用qP反射波对层状介质模型和块状异常体模型进行走时反演,由于qP波相速度对弹性模量参数和Thomsen参数的偏微分不同,所以可以分别反演弹性模量参数和Thomsen参数.数值模拟结果表明:利用qP反射波可以反演出TI介质模型的弹性模量参数与Thomsen参数,不同模型的走时迭代反演达到了较好的收敛效果,与各向同性介质走时反演结果相比较,各向异性介质走时反演结果具有较好的识别能力.     

Pure acoustic wave propagation in transversely isotropic media by the pseudospectral method

Seismic wave propagation in transversely isotropic (TI) media is commonly described by a set of coupled partial differential equations, derived from the acoustic approximation. These equations produce pure P‐wave responses in elliptically anisotropic media but generate undesired shear‐wave components for more general TI anisotropy. Furthermore, these equations suffer from instabilities when the anisotropy parameter ε is less than δ. One solution to both problems is to use pure acoustic anisotropic wave equations, which can produce pure P‐waves without any shear‐wave contaminations in both elliptical and anelliptical TI media. In this paper, we propose a new pure acoustic transversely isotropic wave equation, which can be conveniently solved using the pseudospectral method. Like most other pure acoustic anisotropic wave equations, our equation involves complicated pseudo‐differential operators in space which are difficult to handle using the finite difference method. The advantage of our equation is that all of its model parameters are separable from the spatial differential and pseudo‐differential operators; therefore, the pseudospectral method can be directly applied. We use phase velocity analysis to show that our equation, expressed in a summation form, can be properly truncated to achieve the desired accuracy according to anisotropy strength. This flexibility allows us to save computational time by choosing the right number of summation terms for a given model. We use numerical examples to demonstrate that this new pure acoustic wave equation can produce highly accurate results, completely free from shear‐wave artefacts. This equation can be straightforwardly generalized to tilted TI media.     

Kinematical characteristics of the factorized velocity model

The factorized velocity model that incorporates both vertical heterogeneity and constant anisotropy is one of the complicated analytical models used in seismic data processing and interpretation. In this paper, I derive the analytic equations for offset, traveltime and relative geometrical spreading for the quasi‐compressional (qP‐) waves that can be used for modelling and inversion of the traveltime parameters. I show that the presence of anelliptic anisotropy usually dominates over the vertical heterogeneity with respect to the non‐hyperbolicity of the factorized velocity model.     

Non‐uniqueness with refraction inversion – a syncline model study

Non‐uniqueness occurs with the 1D parametrization of refraction traveltime graphs in the vertical dimension and with the 2D lateral resolution of individual layers in the horizontal dimension. The most common source of non‐uniqueness is the inversion algorithm used to generate the starting model. This study applies 1D, 1.5D and 2D inversion algorithms to traveltime data for a syncline (2D) model, in order to generate starting models for wave path eikonal traveltime tomography. The 1D tau‐p algorithm produced a tomogram with an anticline rather than a syncline and an artefact with a high seismic velocity. The 2D generalized reciprocal method generated tomograms that accurately reproduced the syncline, together with narrow regions at the thalweg with seismic velocities that are less than and greater than the true seismic velocities as well as the true values. It is concluded that 2D inversion algorithms, which explicitly identify forward and reverse traveltime data, are required to generate useful starting models in the near‐surface where irregular refractors are common. The most likely tomogram can be selected as either the simplest model or with a priori information, such as head wave amplitudes. The determination of vertical velocity functions within individual layers is also subject to non‐uniqueness. Depths computed with vertical velocity gradients, which are the default with many tomography programs, are generally 50% greater than those computed with constant velocities for the same traveltime data. The average vertical velocity provides a more accurate measure of depth estimates, where it can be derived. Non‐uniqueness is a fundamental reality with the inversion of all near‐surface seismic refraction data. Unless specific measures are taken to explicitly address non‐uniqueness, then the production of a single refraction tomogram, which fits the traveltime data to sufficient accuracy, does not necessarily demonstrate that the result is either ‘correct’ or the most probable.     

Estimating azimuthal stress‐induced P‐wave anisotropy from S‐wave anisotropy using sonic log or vertical seismic profile data

Most sedimentary rocks are anisotropic, yet it is often difficult to accurately incorporate anisotropy into seismic workflows because analysis of anisotropy requires knowledge of a number of parameters that are difficult to estimate from standard seismic data. In this study, we provide a methodology to infer azimuthal P‐wave anisotropy from S‐wave anisotropy calculated from log or vertical seismic profile data. This methodology involves a number of steps. First, we compute the azimuthal P‐wave anisotropy in the dry medium as a function of the azimuthal S‐wave anisotropy using a rock physics model, which accounts for the stress dependency of seismic wave velocities in dry isotropic elastic media subjected to triaxial compression. Once the P‐wave anisotropy in the dry medium is known, we use the anisotropic Gassmann equations to estimate the anisotropy of the saturated medium. We test this workflow on the log data acquired in the North West Shelf of Australia, where azimuthal anisotropy is likely caused by large differences between minimum and maximum horizontal stresses. The obtained results are compared to azimuthal P‐wave anisotropy obtained via orthorhombic tomography in the same area. In the clean sandstone layers, anisotropy parameters obtained by both methods are fairly consistent. In the shale and shaly sandstone layers, however, there is a significant discrepancy between results since the stress‐induced anisotropy model we use is not applicable to rocks exhibiting intrinsic anisotropy. This methodology could be useful for building the initial anisotropic velocity model for imaging, which is to be refined through migration velocity analysis.     

Traveltime computation with the linearized eikonal equation for anisotropic media

A linearized eikonal equation is developed for transversely isotropic (TI) media with a vertical symmetry axis (VTI). It is linear with respect to perturbations in the horizontal velocity or the anisotropy parameter η. An iterative linearization of the eikonal equation is used as the basis for an algorithm of finite-difference traveltime computations. A practical implementation of this iterative technique is to start with a background model that consists of an elliptically anisotropic, inhomogeneous medium, since traveltimes for this type of medium can be calculated efficiently using eikonal solvers, such as the fast marching method. This constrains the perturbation to changes in the anisotropy parameter η (the parameter most responsible for imaging improvements in anisotropic media). The iterative implementation includes repetitive calculation of η from traveltimes, which is then used to evaluate the perturbation needed for the next round of traveltime calculations using the linearized eikonal equation. Unlike isotropic media, interpolation is needed to estimate η in areas where the traveltime field is independent of η, such as areas where the wave propagates vertically.
Typically, two to three iterations can give sufficient accuracy in traveltimes for imaging applications. The cost of each iteration is slightly less than the cost of a typical eikonal solver. However, this method will ultimately provide traveltime solutions for VTI media. The main limitation of the method is that some smoothness of the medium is required for the iterative implementation to work, especially since we evaluate derivatives of the traveltime field as part of the iterative approach. If a single perturbation is sufficient for the traveltime calculation, which may be the case for weak anisotropy, no smoothness of the medium is necessary. Numerical tests demonstrate the robustness and efficiency of this approach.     

Quantitative imaging of complex structures from dense wide-aperture seismic data by multiscale traveltime and waveform inversions: a case study

An integrated multiscale seismic imaging flow is applied to dense onshore wide‐aperture seismic data recorded in a complex geological setting (thrust belt). An initial P‐wave velocity macromodel is first developed by first‐arrival traveltime tomography. This model is used as an initial guess for subsequent full‐waveform tomography, which leads to greatly improved spatial resolution of the P‐wave velocity model. However, the application of full‐waveform tomography to the high‐frequency part of the source bandwidth is difficult, due to the non‐linearity of this kind of method. Moreover, it is computationally expensive at high frequencies since a finite‐difference method is used to model the wave propagation. Hence, full‐waveform tomography was complemented by asymptotic prestack depth migration to process the full‐source bandwidth and develop a sharp image of the short wavelengths. The final traveltime tomography model and two smoothed versions of the final full‐waveform tomography model were used as a macromodel for the prestack depth migration. In this study, wide‐aperture multifold seismic data are used. After specific preprocessing of the data, 16 frequency components ranging from 5.4 Hz to 20 Hz were inverted in cascade by the full‐waveform tomography algorithm. The full‐waveform tomography successfully imaged SW‐dipping structures previously identified as high‐resistivity bodies. The relevance of the full‐waveform tomography models is demonstrated locally by comparison with a coincident vertical seismic profiling (VSP) log available on the profile. The prestack depth‐migrated images, inferred from the traveltime, and the smoothed full‐waveform tomography macromodels are shown to be, on the whole, consistent with the final full‐waveform tomography model. A more detailed analysis, based on common‐image gather computations, and local comparison with the VSP log revealed that the most accurate migrated sections are those obtained from the full‐waveform tomography macromodels. A resolution analysis suggests that the asymptotic prestack depth migration successfully migrated the wide‐aperture components of the data, allowing medium wavelengths in addition to the short wavelengths of the structure to be imaged. The processing flow that we applied to dense wide‐aperture seismic data is shown to provide a promising approach, complementary to more classical seismic reflection data processing, to quantitative imaging of complex geological structures.     

Traveltime computation for 3D anisotropic media by a finite-difference perturbation method

The first-order perturbation theory is used for fast 3D computation of quasi-compressional (qP)-wave traveltimes in arbitrarily anisotropic media. For efficiency we implement the perturbation approach using a finite-difference (FD) eikonal solver. Traveltimes in the unperturbed reference medium are computed with an FD eikonal solver, while perturbed traveltimes are obtained by adding a traveltime correction to the traveltimes of the reference medium. The traveltime correction must be computed along the raypath in the reference medium. Since the raypath is not determined in FD eikonal solvers, we approximate rays by linear segments corresponding to the direction of the phase normal of plane wavefronts in each cell. An isotropic medium as a reference medium works well for weak anisotropy. Using a medium with ellipsoidal anisotropy as a background medium in the perturbation approach allows us to consider stronger anisotropy without losing computational speed. The traveltime computation in media with ellipsoidal anisotropy using an FD eikonal solver is fast and accurate. The relative error is below 0.5% for the models investigated in this study. Numerical examples show that the reference model with ellipsoidal anisotropy allows us to compute the traveltime for models with strong anisotropy with an improved accuracy compared with the isotropic reference medium.