Finite-difference modelling of S-wave splitting in anisotropic media

We have implemented a 3D finite‐difference scheme to simulate wave propagation in arbitrary anisotropic media. The anisotropic media up to orthorhombic symmetry were modelled using a standard staggered grid scheme and beyond (monoclinic and triclinic) using a rotated staggered grid scheme. The rationale of not using rotated staggered grid for all types of anisotropic media is that the rotated staggered grid schemes are more expensive than standard staggered grid schemes. For a 1D azimuthally anistropic medium, we show a comparison between the seismic data generated by our finite‐difference code and by the reflectivity algorithm; they are in excellent agreement. We conducted a study on zero‐offset shear‐wave splitting using the finite‐difference modelling algorithm using the rotated staggered grid scheme. Our S‐wave splitting study is mainly focused on fractured media. On the scale of seismic wavelenghts, small aligned fractures behave as an equivalent anisotropic medium. We computed the equivalent elastic properties of the fractures and the background in which the fractures were embedded, using low‐frequency equivalent media theories. Wave propagation was simulated for both rotationally invariant and corrugated fractures embedded in an isotropic background for one, or more than one, set of fluid‐filled and dry fractures. S‐wave splitting was studied for dipping fractures, two vertical non‐orthogonal fractures and corrugated fractures. Our modelling results confirm that S‐wave splitting can reveal the fracture infill in the case of dipping fractures. S‐wave splitting has the potential to reveal the angle between the two vertical fractures. We also notice that in the case of vertical corrugated fractures, S‐wave splitting is sensitive to the fracture infill.     

Monte-Carlo Bayesian look-ahead inversion of walkaway vertical seismic profiles

We describe a method to invert a walkaway vertical seismic profile (VSP) and predict elastic properties (P‐wave velocity, S‐wave velocity and density) in a layered model looking ahead of the deepest receiver. Starting from Bayes's rule, we define a posterior distribution of layered models that combines prior information (on the overall variability of and correlations among the elastic properties observed in well logs) with information provided by the VSP data. This posterior distribution of layered models is sampled by a Monte‐Carlo method. The sampled layered models agree with prior information and fit the VSP data, and their overall variability defines the uncertainty in the predicted elastic properties. We apply this technique first to a zero‐offset VSP data set, and show that uncertainty in the long‐wavelength P‐wave velocity structure results in a sizable uncertainty in the predicted elastic properties. We then use walkaway VSP data, which contain information on the long‐wavelength P‐wave velocity (in the reflection moveout) and on S‐wave velocity and density contrasts (in the change of reflectivity with offset). The uncertainty of the look‐ahead prediction is considerably decreased compared with the zero‐offset VSP, and the predicted elastic properties are in good agreement with well‐log measurements.     

Extended reflectivity method for modelling the propagation of diffusive–viscous wave in dip‐layered media

The reflectivity method plays an important role in seismic modelling. It has been used to model different types of waves propagating in elastic and anelastic media. The diffusive–viscous wave equation was proposed to investigate the relationship between frequency dependence of reflections and fluid saturation. It is also used to describe the attenuation property of seismic wave in a fluid‐saturated medium. The attenuation of diffusive–viscous wave is mainly characterised by the effective attenuation parameters in the equation. Thus, it is essential to obtain those parameters and further characterise the features of the diffusive–viscous wave. In this work, we use inversion method to obtain the effective attenuation parameters through quality factor to investigate the characteristics of diffusive–viscous wave by comparing with those of the viscoacoustic wave. Then, the reflection/transmission coefficients in a dip plane‐layered medium are studied through coordinate transform and plane‐wave theory. Consequently, the reflectivity method is extended to compute seismograms of diffusive–viscous wave in a dip plane multi‐layered medium. Finally, we present two models to simulate the propagation of diffusive–viscous wave in a dip plane multi‐layered medium by comparing the results with those in a viscoacoustic medium. The numerical results demonstrate the validity of our extension of reflectivity method to the diffusive–viscous medium. The numerical examples in both time domain and time–frequency domain show that the reflections from a dip plane interface have significant phase shift and amplitude change compared with the results of horizontal plane interface due to the differences in reflection/transmission coefficients. Moreover, the modelling results show strong attenuation and phase shift in the diffusive–viscous wave compared to those of the viscoacoustic wave.     

Elastic wave mode separation for tilted transverse isotropy media

Seismic waves propagate through the earth as a superposition of different wave modes. Seismic imaging in areas characterized by complex geology requires techniques based on accurate reconstruction of the seismic wavefields. A crucial component of the methods in this category, collectively known as wave‐equation migration, is the imaging condition that extracts information about the discontinuities of physical properties from the reconstructed wavefields at every location in space. Conventional acoustic migration techniques image a scalar wavefield representing the P‐wave mode, in contrast to elastic migration techniques, which image a vector wavefield representing both the P‐ and S‐waves. For elastic imaging, it is desirable that the reconstructed vector fields are decomposed into pure wave modes, such that the imaging condition produces interpretable images, characterizing, for example, PP or PS reflectivity. In anisotropic media, wave mode separation can be achieved by projection of the reconstructed vector fields on the polarization vectors characterizing various wave modes. For heterogeneous media, because polarization directions change with position, wave mode separation needs to be implemented using space‐domain filters. For transversely isotropic media with a tilted symmetry axis, the polarization vectors depend on the elastic material parameters, including the tilt angles. Using these parameters, we separate the wave modes by constructing nine filters corresponding to the nine Cartesian components of the three polarization directions at every grid point. Since the S polarization vectors in transverse isotropic media are not defined in the singular directions, e.g., along the symmetry axes, we construct these vectors by exploiting the orthogonality between the SV and SH polarization vectors, as well as their orthogonality with the P polarization vector. This procedure allows one to separate all three modes, with better preserved P‐wave amplitudes than S‐wave amplitudes. Realistic synthetic examples show that this wave mode separation is effective for both 2D and 3D models with strong heterogeneity and anisotropy.     

平面S波在饱和多孔热弹性介质边界的反射问题研究

针对饱和多孔介质中热弹性波的传播特性问题,基于多孔介质理论和广义的热弹性模型,研究平面S波在饱和多孔热弹性介质边界上的反射问题.以考虑流-固耦合的饱和多孔介质波动方程和热-弹耦合的广义热弹性基本方程出发,建立饱和多孔介质的热-流-固耦合弹性波动模型.通过引入势函数并考虑自由透水和绝热的边界条件,经过理论推导最终给出在饱...     

Spectral sparse Bayesian learning reflectivity inversion

A spectral sparse Bayesian learning reflectivity inversion method, combining spectral reflectivity inversion with sparse Bayesian learning, is presented in this paper. The method retrieves a sparse reflectivity series by sequentially adding, deleting or re‐estimating hyper‐parameters, without pre‐setting the number of non‐zero reflectivity spikes. The spikes with the largest amplitude are usually the first to be resolved. The method is tested on a series of data sets, including synthetic data, physical modelling data and field data sets. The results show that the method can identify thin beds below tuning thickness and highlight stratigraphic boundaries. Moreover, the reflectivity series, which is inverted trace‐by‐trace, preserves the lateral continuity of layers.     

L1 norm inversion method for deconvolution in attenuating media

In order to perform a good pulse compression, the conventional spike deconvolution method requires that the wavelet is stationary. However, this requirement is never reached since the seismic wave always suffers high‐frequency attenuation and dispersion as it propagates in real materials. Due to this issue, the data need to pass through some kind of inverse‐Q filter. Most methods attempt to correct the attenuation effect by applying greater gains for high‐frequency components of the signal. The problem with this procedure is that it generally boosts high‐frequency noise. In order to deal with this problem, we present a new inversion method designed to estimate the reflectivity function in attenuating media. The key feature of the proposed method is the use of the least absolute error (L1 norm) to define both the data and model error in the objective functional. The L1 norm is more immune to noise when compared to the usual L2 one, especially when the data are contaminated by discrepant sample values. It also favours sparse reflectivity when used to define the model error in regularization of the inverse problem and also increases the resolution, since an efficient pulse compression is attained. Tests on synthetic and real data demonstrate the efficacy of the method in raising the resolution of the seismic signal without boosting its noise component.     

On the electrical reflectivity tensor in d.c. electric field modelling

We use the integral equation for a d.c. electric field, published in the literature, to introduce the concept of the electrical reflectivity tensor into d.c. electric field modelling. It is shown that in d.c. electric field modelling, the electric reflectivity tensor can be obtained in exactly the same way as in electromagnetic modelling. As a result, for a d.c. electric field, the quasi‐linear and the quasi‐analytical approximations, as well as the quasi‐analytical series, can be constructed in exactly the same way as in electromagnetic modelling. If the primary field is uniform, and if the anomalous body is a uniform circular cylinder or a uniform sphere, the reflectivity tensor is zero order (constant), relating to the free surface charge density. Thus, for some homogeneous bodies that have simple shapes and are embedded in a uniform primary field, the electrical reflectivity tensor is not only a mathematical mechanism for obtaining approximate solutions, but also a physical reality. Indeed, the free surface charge density is defined as the change of the electric displacement vector across the boundary surface under consideration. If the primary field is caused by a point source, and if the anomalous body is a uniform sphere, the reflectivity tensor is second order, varying slowly within the sphere. The relationship to the free surface charge density can be established only when both the reflectivity tensor and the free surface charge density are approximated by the first terms of their series solutions. If the point source is far from the centre of the sphere, the corresponding reflectivity tensor reduces to zero order, and is independent of the observation position within the sphere, i.e. it is a constant. Therefore, the basic idea of the quasi‐analytical approximation, i.e. taking the reflectivity tensor outside the integral operator, is justified in the case considered here.     

基于快速非支配排序遗传算法的VTI介质多分量叠前联合反演

在VTI介质中,由于引入了各向异性参数使得多分量多参数地震反演问题的非线性程度显著增加,因此采用传统的权重加权法构建单目标函数进行反演得到的反演结果往往并不理想.本文以反射率法为基础,结合快速非支配排序遗传算法研究了一种VTI介质的多分量叠前联合反演方法.该方法以反射率法为正演方程,应用互相关原理构建PP波和PSV波的多目标函数,进而采用快速非支配排序遗传算法全局寻优获得VTI介质的厚度、纵横波速度、密度和各向异性参数等多个参数.在正演的过程中,反射率法可以考虑几何扩散、吸收衰减、透射损失、多次波以及纵横波旅行时不匹配等地震波传播效应,更能精确地描述地震波在地下地层中的真实传播情况;在反演的过程中,快速非支配排序遗传算法可以在不引入权重系数的条件下同时优化多个目标函数,获得联合反演问题的Pareto最优解,既不添加权重系数影响又充分利用多分量地震数据.模型测试结果验证了该反演方法的有效性和可行性.     

TTI介质有限差分逆时偏移的稳定性探讨

在沉积学中,可假设在相同时期的沉积层具有相近的物理性质和演化过程.因此,沿层传播的地震波和垂直于地层传播的地震波具有各向异性的特点.在纵波资料的处理中,考虑各向异性对逆时偏移的影响,通常假设介质的横波速度为零,这样可以得到纵波在TTI介质中的传播方程,但是该方程在实际计算中仍存在数值稳定性问题.本文加入横波分量可有效解决数值稳定性问题,并选取适当的横波速度减小对纵波成像的影响,实现地震波在TTI介质中的逆时偏移.实际测算表明,P-SV波的方程中包含横波分量,若假设SV的速度为零,则会导致方程的差分格式不稳定;若加入SV波,选择合适的SV波速度可以使SV波的全区各向异性和反射系数达到极小,并可有效的抑制SV波对纵波勘探的影响.本文的方法是一种稳定的TTI介质中的逆时偏移方法.     

On the retrieval of moment tensors from borehole data

The complete moment tensors of seismic sources in homogeneous or vertically inhomogeneous isotropic structures cannot be retrieved using receivers deployed in one vertical borehole. The complete moment tensors can be retrieved from amplitudes of P‐waves, provided that receivers are deployed in at least three boreholes. Using amplitudes of P‐ and S‐waves, two boreholes are, in principle, sufficient. Similar rules also apply to transversely isotropic media with a vertical axis of symmetry. In the case of limited observations, the inversion can be stabilized by imposing the zero‐trace constraint on the moment tensors. However, this constraint is valid only if applied to observations of shear faulting on planar faults in isotropic media, which produces double‐couple mechanisms. For shear faulting on non‐planar faults, for tensile faulting, and for shear faulting in anisotropic media, the zero‐trace constraint is no longer valid and can distort the retrieved moment tensor and bias the fault‐plane solution. Numerical modelling simulating the inversion of the double‐couple mechanism from real data reveals that the errors in the double‐couple and non‐double‐couple percentages of the moment tensors rapidly decrease with increase in the number of boreholes used. For noisy P‐ and S‐wave amplitudes with noise of 15% of the top amplitude at each channel and for a velocity model biased by 10%, the errors in the double‐couple percentage attain 25, 13 and 6% when inverting for the double‐couple mechanism from one, two and three boreholes.     

2D acoustic‐elastic coupled waveform inversion in the Laplace domain

Although waveform inversion has been intensively studied in an effort to properly delineate the Earth's structures since the early 1980s, most of the time‐ and frequency‐domain waveform inversion algorithms still have critical limitations in their applications to field data. This may be attributed to the highly non‐linear objective function and the unreliable low‐frequency components. To overcome the weaknesses of conventional waveform inversion algorithms, the acoustic Laplace‐domain waveform inversion has been proposed. The Laplace‐domain waveform inversion has been known to provide a long‐wavelength velocity model even for field data, which may be because it employs the zero‐frequency component of the damped wavefield and a well‐behaved logarithmic objective function. However, its applications have been confined to 2D acoustic media. We extend the Laplace‐domain waveform inversion algorithm to a 2D acoustic‐elastic coupled medium, which is encountered in marine exploration environments. In 2D acoustic‐elastic coupled media, the Laplace‐domain pressures behave differently from those of 2D acoustic media, although the overall features are similar to each other. The main differences are that the pressure wavefields for acoustic‐elastic coupled media show negative values even for simple geological structures unlike in acoustic media, when the Laplace damping constant is small and the water depth is shallow. The negative values may result from more complicated wave propagation in elastic media and at fluid‐solid interfaces. Our Laplace‐domain waveform inversion algorithm is also based on the finite‐element method and logarithmic wavefields. To compute gradient direction, we apply the back‐propagation technique. Under the assumption that density is fixed, P‐ and S‐wave velocity models are inverted from the pressure data. We applied our inversion algorithm to the SEG/EAGE salt model and the numerical results showed that the Laplace‐domain waveform inversion successfully recovers the long‐wavelength structures of the P‐ and S‐wave velocity models from the noise‐free data. The models inverted by the Laplace‐domain waveform inversion were able to be successfully used as initial models in the subsequent frequency‐domain waveform inversion, which is performed to describe the short‐wavelength structures of the true models.     

Velocity model building from seismic reflection data by full‐waveform inversion

Full‐waveform inversion is re‐emerging as a powerful data‐fitting procedure for quantitative seismic imaging of the subsurface from wide‐azimuth seismic data. This method is suitable to build high‐resolution velocity models provided that the targeted area is sampled by both diving waves and reflected waves. However, the conventional formulation of full‐waveform inversion prevents the reconstruction of the small wavenumber components of the velocity model when the subsurface is sampled by reflected waves only. This typically occurs as the depth becomes significant with respect to the length of the receiver array. This study first aims to highlight the limits of the conventional form of full‐waveform inversion when applied to seismic reflection data, through a simple canonical example of seismic imaging and to propose a new inversion workflow that overcomes these limitations. The governing idea is to decompose the subsurface model as a background part, which we seek to update and a singular part that corresponds to some prior knowledge of the reflectivity. Forcing this scale uncoupling in the full‐waveform inversion formalism brings out the transmitted wavepaths that connect the sources and receivers to the reflectors in the sensitivity kernel of the full‐waveform inversion, which is otherwise dominated by the migration impulse responses formed by the correlation of the downgoing direct wavefields coming from the shot and receiver positions. This transmission regime makes full‐waveform inversion amenable to the update of the long‐to‐intermediate wavelengths of the background model from the wide scattering‐angle information. However, we show that this prior knowledge of the reflectivity does not prevent the use of a suitable misfit measurement based on cross‐correlation, to avoid cycle‐skipping issues as well as a suitable inversion domain as the pseudo‐depth domain that allows us to preserve the invariant property of the zero‐offset time. This latter feature is useful to avoid updating the reflectivity information at each non‐linear iteration of the full‐waveform inversion, hence considerably reducing the computational cost of the entire workflow. Prior information of the reflectivity in the full‐waveform inversion formalism, a robust misfit function that prevents cycle‐skipping issues and a suitable inversion domain that preserves the seismic invariant are the three key ingredients that should ensure well‐posedness and computational efficiency of full‐waveform inversion algorithms for seismic reflection data.     

S‐wave kinematics in acoustic transversely isotropic media with a vertical symmetry axis

Acoustic transversely isotropic models are widely used in seismic exploration for P‐wave processing and analysis. In isotropic acoustic media only P‐wave can propagate, while in an acoustic transversely isotropic medium both P and S waves propagate. In this paper, we focus on kinematic properties of S‐wave in acoustic transversely isotropic media. We define new parameters better suited for S‐wave kinematics analysis. We also establish the travel time and relative geometrical spreading equations and analyse their properties. To illustrate the behaviour of the S‐wave in multi‐layered acoustic transversely isotropic media, we define the Dix‐type equations that are different from the ones widely used for the P‐wave propagation.     

3D description and inversion of reflection moveout of PS‐waves in anisotropic media

Common‐midpoint moveout of converted waves is generally asymmetric with respect to zero offset and cannot be described by the traveltime series t²(x²) conventionally used for pure modes. Here, we present concise parametric expressions for both common‐midpoint (CMP) and common‐conversion‐point (CCP) gathers of PS‐waves for arbitrary anisotropic, horizontally layered media above a plane dipping reflector. This analytic representation can be used to model 3D (multi‐azimuth) CMP gathers without time‐consuming two‐point ray tracing and to compute attributes of PS moveout such as the slope of the traveltime surface at zero offset and the coordinates of the moveout minimum. In addition to providing an efficient tool for forward modelling, our formalism helps to carry out joint inversion of P and PS data for transverse isotropy with a vertical symmetry axis (VTI media). If the medium above the reflector is laterally homogeneous, P‐wave reflection moveout cannot constrain the depth scale of the model needed for depth migration. Extending our previous results for a single VTI layer, we show that the interval vertical velocities of the P‐ and S‐waves (V_P0 and V_S0) and the Thomsen parameters ε and δ can be found from surface data alone by combining P‐wave moveout with the traveltimes of the converted PS(PSV)‐wave. If the data are acquired only on the dip line (i.e. in 2D), stable parameter estimation requires including the moveout of P‐ and PS‐waves from both a horizontal and a dipping interface. At the first stage of the velocity‐analysis procedure, we build an initial anisotropic model by applying a layer‐stripping algorithm to CMP moveout of P‐ and PS‐waves. To overcome the distorting influence of conversion‐point dispersal on CMP gathers, the interval VTI parameters are refined by collecting the PS data into CCP gathers and repeating the inversion. For 3D surveys with a sufficiently wide range of source–receiver azimuths, it is possible to estimate all four relevant parameters (V_P0, V_S0, ε and δ) using reflections from a single mildly dipping interface. In this case, the P‐wave NMO ellipse determined by 3D (azimuthal) velocity analysis is combined with azimuthally dependent traveltimes of the PS‐wave. On the whole, the joint inversion of P and PS data yields a VTI model suitable for depth migration of P‐waves, as well as processing (e.g. transformation to zero offset) of converted waves.     

Anisotropic reverse-time migration for tilted TI media

Seismic anisotropy in dipping shales results in imaging and positioning problems for underlying structures. We develop an anisotropic reverse‐time depth migration approach for P‐wave and SV‐wave seismic data in transversely isotropic (TI) media with a tilted axis of symmetry normal to bedding. Based on an accurate phase velocity formula and dispersion relationships for weak anisotropy, we derive the wave equation for P‐wave and SV‐wave propagation in tilted transversely isotropic (TTI) media. The accuracy of the P‐wave equation and the SV‐wave equation is analyzed and compared with other acoustic wave equations for TTI media. Using this analysis and the pseudo‐spectral method, we apply reverse‐time migration to numerical and physical‐model data. According to the comparison between the isotropic and anisotropic migration results, the anisotropic reverse‐time depth migration offers significant improvements in positioning and reflector continuity over those obtained using isotropic algorithms.     

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.     

Double-square-root one-way wave equation prestack tau migration in heterogeneous media

In this paper, source‐receiver migration based on the double‐square‐root one‐way wave equation is modified to operate in the two‐way vertical traveltime (τ) domain. This tau migration method includes reasonable treatment for media with lateral inhomogeneity. It is implemented by recursive wavefield extrapolation with a frequency‐wavenumber domain phase shift in a constant background medium, followed by a phase correction in the frequency‐space domain, which accommodates moderate lateral velocity variations. More advanced τ‐domain double‐square‐root wave propagators have been conceptually discussed in this paper for migration in media with stronger lateral velocity variations. To address the problems that the full 3D double‐square‐root equation prestack tau migration could meet in practical applications, we present a method for downward continuing common‐azimuth data, which is based on a stationary‐phase approximation of the full 3D migration operator in the theoretical frame of prestack tau migration of cross‐line constant offset data. Migrations of synthetic data sets show that our tau migration approach has good performance in strong contrast media. The real data example demonstrates that common‐azimuth prestack tau migration has improved the delineation of the geological structures and stratigraphic configurations in a complex fault area. Prestack tau migration has some inherent robust characteristics usually associated with prestack time migration. It follows a velocity‐independent anti‐aliasing criterion that generally leads to reduction of the computation cost for typical vertical velocity variations. Moreover, this τ‐domain source‐receiver migration method has features that could be of help to speed up the convergence of the velocity estimation.