Interpretation of shallow refraction seismic data by reflection/refraction tomography

A new algorithm for tomographic inversion of traveltimes of reflected and refracted seismic waves is developed. The inversion gives interface configurations and velocity distributions in layers. The important features of the algorithm are: (a) the inclusion of shot time delays in the list of unknown parameters; (b) the regularization is applied in such a way that the most probable model is characterized by the similarity of neighbouring interfaces. As the problem under consideration is non-linear, several iterations are necessary in order to obtain the final model. In the case of a very inexact initial model, a 'layer-by-layer' inversion strategy is recommended as a first inversion step. The inversion program is supplied with a user interface, thanks to which one can: (a) pick interactively and identify seismic traveltimes; (b) build and edit depth/velocity models; and (c) display calculated traveltime curves and compare them with picked traveltimes as well as with the original seismic sections. The efficiency of the inversion software developed is illustrated by a numerical example and a field example in which shallow seismic data are considered. Application to wide-aperture reflection/refraction profiling (WARRP) data is also possible.     

Depth and morphology of reflectors from the non-linear inversion of arrival times and waveform semblance data. Part II: modelling and interpretation of real data acquired in Southern Apennines, Italy

In order to retrieve a 2D background velocity model and to retrieve the geometry and depth of shallow crustal reflectors in the Southern Apennines thrust belt a separate inversion of first arrival traveltimes and reflected waveforms was performed. Data were collected during an active seismic experiment in 1999 by Enterprise Oil Italiana and Eni-Agip using a global offset acquisition geometry. A total of 284 on-land shots were recorded by 201 receivers deployed on an 18 km line oriented SW–NE in the Val D'Agri region (Southern Apennines, Italy).
The two-step procedure allows for the retrieval of a reliable velocity model by using a non-linear tomographic inversion and reflected waveform semblance data inversion. The tomographic model shows that the P wave velocity field varies vertically from approximately 3 km/s to 6 km/s within 4 km from the Earth's surface. Moreover, at a distance of approximately 11 km along the profile, there is an abrupt increase in the velocity field. In this zone indeed, an ascent from 2 km depth to 0 km above sea level of the 5.2 km/s iso-velocity contour can be noted. The retrieved velocity can be associated with Plio-Pleistocene clastic deposits outcropping in the basin zone and with Mesozoic limestone deposits. The inversion of waveform semblance data shows that a P-to-P reflector is retrieved at a depth of approximately 2 km. This interface is deeper in the north-eastern part of the profile, where it reaches 3 km depth and can be associated with a limestone horizon.     

A SIMULATED ANNEALING APPROACH TO SEISMIC MODEL OPTIMIZATION WITH SPARSE PRIOR INFORMATION1

It is well known that seismic inversion based on local model optimization methods, such as iterative use of linear optimization, may fail when prior information is sparse. Where the seismic events corresponding to reflectors of interest remain to be identified, a global optimization technique is required. We investigate the use of a global, stochastic optimization method, that of simulated annealing, to solve the seismic trace inversion problem, in which the two-way traveltimes and reflection coefficients are to be determined. The simulated annealing method is based on an analogy between the model-algorithm system and a statistical mechanical system. We exploit this analogy to produce improved annealing schedules. It is shown that even in cases of virtually no prior information about two-way traveltimes and reflection coefficients, the method is capable of producing reliable results.     

一种改进的地震反射层析成像方法

针对复杂介质的地震反射走时层析成像存在数据拾取困难问题，本文提出了一种新的地震反射层析成像速度模型建立方法，该方法用速度和地震射线走时描述模型，用地震反射波走时、地震波在源点和接收点处的传播方向信息反演模型.为提高反演的稳定性和计算效率，引入了Hamilton函数描述射线，在相空间计算反演所需的射线路径和目标函数对模型参数的导数，对理论模型和实际地震资料进行了试算，试算表明该方法对复杂介质具有较强的适应能力.     

Estimation of petrophysical parameters by linearized inversion of angle domain pre-stack data

We describe a linear Bayesian inversion method to estimate the relevant petrophysical properties of the media forming a reflecting interface from the observations of amplitude variation with incidence angle. Three main steps characterize the proposed approach:
– information from borehole logs are statistically analysed to estimate the empirical models that describe the functional relationship between petrophysical (e.g. porosity, saturation, pressure or depth) and seismic variable(P and S velocities and density);
– the pure-mode (PP) reflection coefficient is parameterized in terms of the relevant petrophysical variables and is linearized in order to implement the linear inversion;
– the sought petrophysical parameters are estimated from the seismic reflected amplitudes by applying the linearized inversion where a priori information, data and model errors and solutions are described by probability density functions.
We test the method on synthetic and real data relative to reflections from a shale/gas-sand interface where the amplitude versus angle response, besides the lithological contrast, is mainly controlled by the saturation and porosity of the sand layer. The outcomes of the linearized inversion are almost identical to those obtained by a previously developed non-linear inversion method demonstrating the applicability of the linear inversion. It turns out that the gas-sand saturation in the range 0%–95% is a poorly resolved parameter while the porosity is the best resolved parameter. The issues of robustness and resolution of the inversion are discussed either through singular value decomposition analysis or the observation of the a posteriori probability density functions.
The linear inversion algorithm, compared with the previously developed non-linear method, reduces significantly the computation time allowing for more extensive applications.     

Resolution assessment of the three-dimensional acquisition geometry of the Anhui experiment in "Geoscience Yangtze Plan"

According to the actual observation conditions of the Yangtze River valley from Anqing city to Maanshan city, we designed the 3D acquisition geometry, and applied the multi-scale checkerboard semblance analysis to assess the preliminary resolution of the designed observation. The checkerboard semblance tests use the refraction and reflection travel-time simultaneous inversion algorithm to quantitatively provide both resolution assessment of velocity structure and Moho topography. The multi-scale checkerboard semblance recovery results show that while the checkerboard semblance threshold value is 0.5, the preliminary resolution of the designed acquisition geometry is better than 10 km in the upper crust (the depth is less than 10 km), around 15 km in the mid-crust (the depth is 10?25 km), and better than 20 km in the lower crust (the depth is 25?33 km). The preliminary tomographic resolution for the Moho topography is about 20 km in the ray-path coverage area beneath the acquisition geometry. While the checkerboard semblance threshold value is 0.75, the preliminary resolution is 20 km in the upper crust, around 20?25 km in the mid-crust and 25 km in the lower crust. And the preliminary tomographic resolution for the Moho topography is better than 30 km in the ray-path coverage area beneath the acquisition geometry. These non-linear checkerboard tests reveal that the designed acquisition geometry is suitable to image the crustal velocity structure of the Yangtze River valley in the Anhui province.     

INVERSION OF POST-STACK SEISMIC DATA USING SIMULATED ANNEALING1

Model-based inversion of seismic reflection data is a global optimization problem when prior information is sparse. We investigate the use of an efficient, global, stochastic optimization method, that of simulated annealing, for determining the two-way traveltimes and the reflection coefficients. We exploit the advantage of an ensemble approach to the inversion of full-scale target zones on 2D seismic sections. In our ensemble approach, several copies of the model-algorithm system are run in parallel. In this way, estimation of true ensemble statistics for the process is made possible, and improved annealing schedules can be produced. It is shown that the method can produce reliable results efficiently in the 2D case, even when prior information is sparse.     

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

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

Migration velocity analysis and waveform inversion

Least‐squares inversion of seismic reflection waveform data can reconstruct remarkably detailed models of subsurface structure and take into account essentially any physics of seismic wave propagation that can be modelled. However, the waveform inversion objective has many spurious local minima, hence convergence of descent methods (mandatory because of problem size) to useful Earth models requires accurate initial estimates of long‐scale velocity structure. Migration velocity analysis, on the other hand, is capable of correcting substantially erroneous initial estimates of velocity at long scales. Migration velocity analysis is based on prestack depth migration, which is in turn based on linearized acoustic modelling (Born or single‐scattering approximation). Two major variants of prestack depth migration, using binning of surface data and Claerbout's survey‐sinking concept respectively, are in widespread use. Each type of prestack migration produces an image volume depending on redundant parameters and supplies a condition on the image volume, which expresses consistency between data and velocity model and is hence a basis for velocity analysis. The survey‐sinking (depth‐oriented) approach to prestack migration is less subject to kinematic artefacts than is the binning‐based (surface‐oriented) approach. Because kinematic artefacts strongly violate the consistency or semblance conditions, this observation suggests that velocity analysis based on depth‐oriented prestack migration may be more appropriate in kinematically complex areas. Appropriate choice of objective (differential semblance) turns either form of migration velocity analysis into an optimization problem, for which Newton‐like methods exhibit little tendency to stagnate at nonglobal minima. The extended modelling concept links migration velocity analysis to the apparently unrelated waveform inversion approach to estimation of Earth structure: from this point of view, migration velocity analysis is a solution method for the linearized waveform inversion problem. Extended modelling also provides a basis for a nonlinear generalization of migration velocity analysis. Preliminary numerical evidence suggests a new approach to nonlinear waveform inversion, which may combine the global convergence of velocity analysis with the physical fidelity of model‐based data fitting.     

2D多尺度非线性地震速度成像

将遗传算法和单纯形算法相结合,得到了一种高效、健全的2D混合地震走时反演方法.把速度场划分为不同的空间尺度,定义网格节点上的速度作为待反演参数,采用双三次样条函数模型参数化,正问题采用有限差分走时计算方法,反问题采用多尺度混合反演方法.首先在较大的空间尺度内反演,然后减小空间尺度,将大尺度的反演结果作为次一级尺度反问题的初始模型,再进行混合反演,如此类推逐次逼近全局最优解.一个低速度异常体的数值模拟试验和抗走时扰动试验表明该方法是有效和健全的.我们将该方法应用到青藏高原东北缘阿尼玛卿缝合带东段上部地壳速度结构研究中,并与前人的成果进行了对比.     

Determination of a shallow velocity–depth model from seismic refraction data by coherence inversion1

Seismic refractions have different applications in seismic prospecting. The travel- times of refracted waves can be observed as first breaks on shot records and used for field static calculation. A new method for constructing a near-surface model from refraction events is described. It does not require event picking on prestack records and is not based on any approximation of arrival times. It consists of the maximization of the semblance coherence measure computed using shot gathers in a time window along refraction traveltimes. Time curves are generated by ray tracing through the model. The initial model for the inversion was constructed by the intercept-time method. Apparent velocities and intercept times were taken from a refraction stacked section. Such a section can be obtained by appling linea moveout corrections to common-shot records. The technique is tested successfully on synthetic and real data. An important application of the proposed method for solving the statics problem is demonstrated.     

Automatic cross-well tomography: an application of the differential semblance optimization to two real examples

An automatic tomography algorithm, based on differential semblance optimization (DSO), has been used to invert real cross-well seismic data for the background velocity. The method relies on the first-arrival transmitted waves. Given a background velocity model, the traveltimes between the sources and the receivers are computed, then semblance panels are created by back-propagating the data traces. If the velocity model is correct all the first-arrival transmitted waves will be aligned in the semblance panels. The DSO method consists of finding the background velocity by minimizing the L ₂-norm of the difference between adjacent back-propagated traces. Thanks to the good behaviour of this DSO cost function about the solution, a local (gradient) optimization can be performed. This provides a relatively fast algorithm when ray tracing and analytic computation of the gradient are used.
Unfortunately the method fails in the presence of caustics in the data. However, this difficulty can be circumvented by applying suitable masks to the data. This approach is first applied to a synthetic example then to two real data sets: the McElroy data set recorded in West Texas and the NIMR data set recorded in Oman. The results are quite encouraging and similar to those obtained with classical tomography.     

The effectiveness of a pseudo-inverse extended Born operator to handle lateral heterogeneity for imaging and velocity analysis applications

Wave equation–based migration velocity analysis techniques aim to construct a kinematically accurate velocity model for imaging or as an initial model for full waveform inversion applications. The most popular wave equation–based migration velocity analysis method is differential semblance optimization, where the velocity model is iteratively updated by minimizing the unfocused energy in an extended image volume. However, differential semblance optimization suffers from artefacts, courtesy of the adjoint operator used in imaging, leading to poor convergence. Recent findings show that true amplitude imaging plays a significant role in enhancing the differential semblance optimization's gradient and reducing the artefacts. Here, we focus on a pseudo-inverse operator to the horizontally extended Born as a true amplitude imaging operator. For laterally inhomogeneous models, the operator required a derivative with respect to a vertical shift. Extending the image vertically to evaluate such a derivative is costly and impractical. The inverse operator can be simplified in laterally homogeneous models. We derive an extension of the approach to apply the full inverse formula and evaluate the derivative efficiently. We simplified the implementation by applying the derivative to the imaging condition and utilize the relationship between the source and receiver wavefields and the vertical shift. Specifically, we verify the effectiveness of the approach using the Marmousi model and show that the term required for the lateral inhomogeneity treatment has a relatively small impact on the results for many cases. We then apply the operator in differential semblance optimization and invert for an accurate macro-velocity model, which can serve as an initial velocity model for full waveform inversion.     

Background velocity estimation using non-linear optimization for reflection tomography and migration misfit

We show that it is possible to estimate the background velocity for prestack depth migration in 2D laterally varying media using a non-linear optimization technique called very fast simulated annealing (VFSA). We use cubic splines in the velocity model parametrization and make use of either successive pairs of shot gathers or several constant-offset sections as input data for the inversion. A Kirchhoff summation scheme based on first-arrival traveltimes is used to migrate/model the input data during the velocity analysis. We evaluate and compare two different measures of error. The first is defined in the recorded data or (x,t) domain and is based on a reflection-tomography criterion. The second is defined in the migrated data or (x,z) domain and is based on a migration-misfit criterion. Depth relaxation is used to improve the convergence and quality of the velocity analysis while simultaneously reducing the computational cost. Further, we show that by coarse sampling in the offset domain the method is still robust. Our non-linear optimization approach to migration velocity analysis is evaluated for both synthetic and real seismic data. For the velocity-analysis method based on the reflection-tomography criterion, traveltimes do not have to be picked. Similarly, the migration-misfit criterion does not require that depth images be manually compared. Interpreter intervention is required only to restrict the search space used in the velocity-analysis problem. Extension of the proposed schemes to 3D models is straightforward but practical only for the fastest available computers.     

Stacking velocities in the presence of overburden velocity anomalies

Lateral velocity changes (velocity anomalies) in the overburden may cause significant oscillations in normal moveout velocities. Explicit analytical moveout formulas are presented and provide a direct explanation of these lateral fluctuations and other phenomena for a subsurface with gentle deep structures and shallow overburden anomalies. The analytical conditions for this have been derived for a depth-velocity model with gentle structures with dips not exceeding 12°. The influence of lateral interval velocity changes and curvilinear overburden velocity boundaries can be estimated and analysed using these formulas. An analytical approach to normal moveout velocity analysis in a laterally inhomogeneous medium provides an understanding of the connection between lateral interval velocity changes and normal moveout velocities. In the presence of uncorrected shallow velocity anomalies, the difference between root-mean-square and stacking velocity can be arbitrarily large to the extent of reversing the normal moveout function around normal incidence traveltimes. The main reason for anomalous stacking velocity behaviour is non-linear lateral variations in the shallow overburden interval velocities or the velocity boundaries.
A special technique has been developed to determine and remove shallow velocity anomaly effects. This technique includes automatic continuous velocity picking, an inversion method for the determination of shallow velocity anomalies, improving the depth-velocity model by an optimization approach to traveltime inversion (layered reflection tomography) and shallow velocity anomaly replacement. Model and field data examples are used to illustrate this technique.     

Enhancement of prestack diffraction data and attributes using a traveltime decomposition approach

Diffractions not only carry important information about small-scale subsurface structures, they also possess unique properties, which make them a powerful tool for seismic processing and imaging. Since a point diffractor scatters an incoming wave to all directions, a diffraction event implies better illumination than a reflection, because the rays travel through larger parts of the subsurface. Furthermore, unlike the reflection case, in which the emergence location of the reflected wave depends on the source position, in the case of non-Snell scattering, up-going and down-going raypaths are decoupled. Based on this decoupling, we introduce a diffraction traveltime decomposition principle, which establishes a direct connection between zero-offset and finite-offset diffraction wavefield attributes. By making use of this approach, we are able to enhance diffractions and obtain high-quality diffraction wavefield attributes at arbitrary offsets in the prestack domain solely based on zero-offset processing without any further optimization of attributes. We show the accuracy of the method by fitting diffraction traveltimes, and on simple waveform data. Application to complex synthetic data shows the ability of the proposed approach to enhance diffractions and provide high-quality wavefield attributes even in sparsely illuminated regions such as subsalt areas. The promising results reveal a high potential for improved prestack data enhancement and further applications such as efficient diffraction-based finite-offset tomography.     

杨氏模量和泊松比反射系数近似方程及叠前地震反演

杨氏模量和泊松比等岩石弹性参数是表征页岩气储集体岩石脆性、评价储层含气特征的重要特征参数,而叠前地震反演是从地震资料中获取此类参数的有效途径.地震波反射系数方程是叠前反演的基础.首先,在平面波入射等假设条件下推导了基于杨氏模量(Y)、泊松比(σ)和密度(D)的纵波反射系数线性近似方程(YPD反射系数近似方程),该方程建立了地震纵波反射系数与杨氏模量反射系数、泊松比反射系数和密度反射系数的线性关系;其次,对该方程的精度和适用条件进行了分析,模型分析表明,在入射角为40°时,该方程具有较高的计算精度;最后,建立了一种稳定获取杨氏模量和泊松比的叠前地震直接反演方法.模型试算和实际资料试处理表明,基于新方程的反演方法能够稳定合理的直接从叠前地震资料中获取杨氏模量和泊松比参数,提供了一种高可靠性的页岩气"甜点"地震识别方法.     

Efficient and accurate computation of seismic traveltimes and amplitudes

We describe two practicable approaches for an efficient computation of seismic traveltimes and amplitudes. The first approach is based on a combined finite‐difference solution of the eikonal equation and the transport equation (the ‘FD approach’). These equations are formulated as hyperbolic conservation laws; the eikonal equation is solved numerically by a third‐order ENO–Godunov scheme for the traveltimes whereas the transport equation is solved by a first‐order upwind scheme for the amplitudes. The schemes are implemented in 2D using polar coordinates. The results are first‐arrival traveltimes and the corresponding amplitudes. The second approach uses ray tracing (the ‘ray approach’) and employs a wavefront construction (WFC) method to calculate the traveltimes. Geometrical spreading factors are then computed from these traveltimes via the ray propagator without the need for dynamic ray tracing or numerical differentiation. With this procedure it is also possible to obtain multivalued traveltimes and the corresponding geometrical spreading factors. Both methods are compared using the Marmousi model. The results show that the FD eikonal traveltimes are highly accurate and perfectly match the WFC traveltimes. The resulting FD amplitudes are smooth and consistent with the geometrical spreading factors obtained from the ray approach. Hence, both approaches can be used for fast and reliable computation of seismic first‐arrival traveltimes and amplitudes in complex models. In addition, the capabilities of the ray approach for computing traveltimes and spreading factors of later arrivals are demonstrated with the help of the Shell benchmark model.     

Inversion for elliptically anisotropic velocity using VSP reflection traveltimes

This paper presents a traveltime inversion approach, using the reflection traveltimes from offset VSP data, to reconstruct the horizontal and vertical velocities for stratified anisotropic media. The inverse problem is reduced to a set of linear equations, and solved by the singular value decomposition (SVD) technique. The validity of this inversion scheme is verified using two sets of synthetic data simulated using a finite‐difference method, one for an isotropic model and the other for an elliptically anisotropic model. The inversion result demonstrates that our anisotropic velocity inversion scheme may be applied to both isotropic and anisotropic media. The method is finally applied to a real offset VSP data set, acquired in an oilfield in northwestern China.     

平均入射角道集PP波与PS波联合反演

在界面两侧地层的弹性参数弱反差的假设难以成立的情况下,本文提出用平均入射角道集进行PP波与PS波的联合反演.首先,在PP波与PS波AVA(amplitude versus angle,振幅随入射角变化)道集的基础上,分别选择小入射角范围与大入射角范围的AVA道集进行局部加权叠加,以获得由两个角度组成的平均入射角道集,并作为后续反演的输入数据.然后,再通过最小二乘原理建立了PP波与PS波联合反演目标函数,推导了模型修改量的向量公式,建立了平均入射角道集联合反演的流程.模型数据与实际数据的测试结果表明:在信噪比较低、地层弹性参数反差较大、层厚较薄的情况下,该反演方法的精度在很大程度上超过了基于近似反射系数的反演方法,为复杂油气藏勘探提供了新的思路.