The problem of velocity inversion in refraction seismics: some observations from modelling results

The applicability of seismic refraction profiling for the detection of velocity inversion, which is also known as a low-velocity layer (LVL), is investigated with the aid of synthetic seismogram computations for a range of models. Our computational models focus on the inherent ambiguities in the interpretation of first-arrival time delays or 'skips' in terms of LVL model parameters. The present modelling results reveal that neither the measure nor even the existence of a shadow zone and/or a time shift (skip) in first arrivals is necessarily indicative of an LVL. Besides attenuation effects, the cap-layer velocity gradient is a critical parameter, determining the termination point of the cap-layer diving wave and thus the time skip.
We suggest that shallow LVLs can be delineated more reliably by traveltime and amplitude modelling of coherent phases reflected from their top and bottom boundaries, often clearly observed in the pre- and near-critical ranges in seismogram sections of refraction profiling experiments with a close receiver spacing. We demonstrate the applicability of this approach for a field data set of a refraction profile in the West Bengal Basin, India. The inferred LVL corresponds to the Gondwana sediments underlying the higher-velocity layer of the Rajmahal Traps. This interpretation is consistent with the data from a nearby well in the region.     

Assessing uncertainty in refraction seismic traveltime inversion using a global inversion strategy

To analyse and invert refraction seismic travel time data, different approaches and techniques have been proposed. One common approach is to invert first‐break travel times employing local optimization approaches. However, these approaches result in a single velocity model, and it is difficult to assess the quality and to quantify uncertainties and non‐uniqueness of the found solution. To address these problems, we propose an inversion strategy relying on a global optimization approach known as particle swarm optimization. With this approach we generate an ensemble of acceptable velocity models, i.e., models explaining our data equally well. We test and evaluate our approach using synthetic seismic travel times and field data collected across a creeping hillslope in the Austrian Alps. Our synthetic study mimics a layered near‐surface environment, including a sharp velocity increase with depth and complex refractor topography. Analysing the generated ensemble of acceptable solutions using different statistical measures demonstrates that our inversion strategy is able to reconstruct the input velocity model, including reasonable, quantitative estimates of uncertainty. Our field data set is inverted, employing the same strategy, and we further compare our results with the velocity model obtained by a standard local optimization approach and the information from a nearby borehole. This comparison shows that both inversion strategies result in geologically reasonable models (in agreement with the borehole information). However, analysing the model variability of the ensemble generated using our global approach indicates that the result of the local optimization approach is part of this model ensemble. Our results show the benefit of employing a global inversion strategy to generate near‐surface velocity models from refraction seismic data sets, especially in cases where no detailed a priori information regarding subsurface structures and velocity variations is available.     

一种基于高效迭代解法的频率域全波形反演

巨大的计算量是制约全波形反演（FWI）生产实用化的难题之一.为此,本文提出了一种高效的波场迭代解法,将其应用于频率域常密度声波方程FWI,并给出了详细的反演流程.通过建立用于波场迭代的目标函数,推导相应梯度、步长公式,新方法将反演中波场正传和残差波场反传过程转化为无约束优化问题,从理论上分析了新方法的计算效率显著高于常规FWI.在数值试验中,本文方法通过几次迭代便能获得高精度的正传、残差反传波场,收敛速度明显高于未经预处理的GMRES方法.进一步引入高效编码策略,新方法的计算时间约为常规编码FWI的1/8,与理论分析结果吻合（波场迭代次数为8,模型未知量个数约为7万）,且波场迭代次数为6时,反演效果已与常规编码FWI相近.

     

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.     

The detection of King Xerxes' Canal by the use of shallow reflection and refraction seismics — Preliminary results

Shallow reflection and refraction seismic studies were carried out in Greece in the eastern leg of the Chalkidiki peninsula, in order to test the validity of reports in history books which describe a legendary canal built by the engineers of King Xerxes during the major Persian invasion of Europe through Greece in the 5th century B.C.   In the narrowest part of the Athos peninsula, where it is 2 km wide, an 85 m profile was topographically defined almost centrally between the two coastlines. The position of this profile was based on palaeogeographical, geomorphological and topographic studies and observations. A sledgehammer was used as the seismic source for the shallow target. Despite the presence of significant urban and coherent noise, a final stacked section was produced by a suitable choice of acquisition and processing parameters.   Both the reflection and refraction seismic studies illustrated clearly the existence of a channel-like structure of trapezoid cross-section, almost midway between the two opposite sides of the peninsula.     

Simple inversion procedure for shallow seismic refraction in continuously nonhomogeneous soils

Seismic refraction analysis is presented for uniformly deposited shallow soil strata with wave velocity increasing continuously with depth due to differential compaction effects. Closed form solutions of the surface-to-surface travel time as well as for the depth to maximum penetration are derived for realistic velocity-depth functions. An inversion procedure based on simple formulas is presented and its applicability is demonstrated and discussed. The method takes into account the information along the entire refraction line in one step, in contrast to the discrete ray-tracing technique and is, therefore, less sensitive to the natural scatter of the data.     

折射波静校正与层析静校正技术适用性分析

对于任何一种类型的油气藏,构造都是十分重要的主控因素,而构造因素的影响往往与地震资料处理中的静校正密切相关,因此,长波长静校正问题越来越多地引起人们的关注.折射波静校正与层析静校正是目前资料处理中常用的静校正技术.前者适合近地表有稳定的折射层、表层速度和厚度纵横向变化不太剧烈的地区,并要预先给出风化层的速度;但是不能反映层内速度变化,不适合风化层速度变化剧烈或存在速度反转的地区.后者在算法上比前者复杂许多,适合地形起伏较大、近地表结构复杂、速度变化剧烈的地区,能够解决因速度横向剧烈变化和速度倒转引起的静校正问题;但是计算量大,计算时间长,存在多解性.日常工作中往往会认为层析静校正一定比折射波静校正效果好,其实不然.本文在归纳了这两种方法适用条件和优缺点的基础上,结合低山区和黄土塬区两种典型地表结构的实际资料进行了静校正处理.从模型和剖面对比上可以证明,并不是算法越复杂处理效果就越好,只有在满足算法适用条件的前提下才能充分体现技术优势.     

Seismic refraction traveltime inversion for static corrections in a glaciated shield rock environment: a case study

A traveltime inversion technique is applied to model the upper ∼40 m of the subsurface of a glaciated shield rock area in order to calculate static corrections for a multi-azimuth multi-depth walk-away vertical seismic profile and a surface seismic reflection profile. First break information from a seismic refraction survey is used in conjunction with a ray-tracing program and an iterative damped least-squares inversion algorithm to create a two-dimensional model of the subsurface. The layout of the seismic survey required crooked seismic lines and substantial gaps in the source and receiver coverage to be accounted for. Additionally, there is substantial topographical variation and a complex geology consisting of glaciofluvial sediment and glacial till overlying a crystalline bedrock. The resolution and reliability of the models is measured through a parameter perturbation technique, normalized χ² values, root means square traveltime residuals and comparison to known geology.     

3D acoustic waveform inversion in the Laplace domain using an iterative solver

In this paper we propose a 3D acoustic full waveform inversion algorithm in the Laplace domain. The partial differential equation for the 3D acoustic wave equation in the Laplace domain is reformulated as a linear system of algebraic equations using the finite element method and the resulting linear system is solved by a preconditioned conjugate gradient method. The numerical solutions obtained by our modelling algorithm are verified through a comparison with the corresponding analytical solutions and the appropriate dispersion analysis. In the Laplace‐domain waveform inversion, the logarithm of the Laplace transformed wavefields mainly contains long‐wavelength information about the underlying velocity model. As a result, the algorithm smoothes a small‐scale structure but roughly identifies large‐scale features within a certain depth determined by the range of offsets and Laplace damping constants employed. Our algorithm thus provides a useful complementary process to time‐ or frequency‐domain waveform inversion, which cannot recover a large‐scale structure when low‐frequency signals are weak or absent. The algorithm is demonstrated on a synthetic example: the SEG/EAGE 3D salt‐dome model. The numerical test is limited to a Laplace‐domain synthetic data set for the inversion. In order to verify the usefulness of the inverted velocity model, we perform the 3D reverse time migration. The migration results show that our inversion results can be used as an initial model for the subsequent high‐resolution waveform inversion. Further studies are needed to perform the inversion using time‐domain synthetic data with noise or real data, thereby investigating robustness to noise.     

Palaeomagnetic distortion modelling and possible recovery by inversion

A robust finite-element technique is presented for computation of both the internal demagnetization effects and magnetic terrain effects in bodies with arbitrary shape and arbitrary susceptibility distribution. This method facilitates a flexible analysis of the palaeomagnetic deflection problem. Tests on geologically realistic settings of highly magnetic rocks demonstrate that deflections of several degrees may occur even for relatively simple two-dimensional models. Similarly, the magnetic intensity may well be biased by 5-15% by demagnetization effects. The present paper focuses on deflections and intensity variations inside the magnetized body, where we find a systematic shallowing of inclination for bodies with a horizontal elongation. Because the bodies sampled at a typical palaeomagnetic site will have a dominant direction of elongation, the magnetic deflection effect will tend to impose a systematic bias which doesn’t average out. An inversion-based procedure for elimination of the deflection effect is presented. It requires that the magnetic body is quite homogeneous and that its surface geometry is known, as may be the case for historical lava flows. Tests demonstrate that in order to recover both ambient palaeofield direction and the effective susceptibility at blocking temperature it is necessary to sample near strong topographic elements in the magnetic body. Since the surface geometry rarely is known it is proposed as an alternative to inversion that an effective susceptibility is assessed and a horizontal slab correction is applied for samples taken far from topographical features. When shape geometry is unknown and no correction applied, palaeomagnetic conclusions must take into account the possible bias from internal demagnetization and magnetic terrain effects.     

Use of traveltime skips in refraction analysis to delineate velocity inversion1

First arrival refraction data does not normally provide any indication of the velocity inversion problem. However, under certain favourable circumstances, when the low-velocity layer (LVL) is considerably thicker than the overlying higher-velocity layer (HVL), the velocity inversion can be seen in the form of a traveltime skip. Model Studies show that in such cases the length of the HVL traveltime branch can be used to determine the thickness of the HVL and the magnitude of the traveltime skip in order to determine the thickness of the LVL. This is also applicable in the case of field data.     

Recursive seismic ray modelling: applications in inversion and VSP

The recursive nature of rays in blocky models can be exploited to solve some difficult problems in seismic modelling. Each segment of a ray travels from an initial point up to a reflecting interface, where it is split into reflected and transmitted ray segments, which each continue in a similar way. The tree structure that thus emanates is conveniently handled by a recursive scheme. Recursion allows an automatic generation of all phases on a seismogram, together with all information necessary to analyse or select them. By operating recursively with a ray cell, bounded by a pair of vicinal rays in 2D, or a triplet of vicinal rays in 3D, and two successive isochrons, the two-point ray-tracing problem is reduced to a simple interpolation. Also, the cellular approach allows for a stable and robust evaluation of dynamic ray quantities without any paraxial tracing, which is cumbersome in blocky models of realistic complexity. Geometric shadows are filled by recursively generated diffractions. The recursive ray tracer has found applications in the fast computation of Green's functions in target-oriented inversion and in phase identification in VSP.     

Elongated horizontal and vertical receivers in three-dimensional electromagnetic modelling and inversion

Electromagnetic geophysical methods often rely on measurements of naturally occurring or artificially impressed electric fields. It is technically impossible, however, to measure the electric field directly. Instead, the electric field is approximated by recording the voltage difference between two electrodes and dividing the obtained voltage by the distance between the electrodes. Typically, modelling and inversion algorithms assume that the electric fields are obtained over infinitely short point-dipoles and thus measured fields are assigned to a single point between the electrodes. Such procedures imply several assumptions: (1) The electric field between the two electrodes is regarded as constant or being a potential field and (2) the receiver dimensions are negligible compared to the dimensions of the underlying modelling grid. While these conditions are often fulfilled for horizontal electric fields, borehole sensors for recordings of the vertical electric field have dimensions in the order of ≈100 m and span several modelling grid cells. Observations from such elongated borehole sensors can therefore only be interpreted properly if true receiver dimensions and variations of electrical conductivity along the receiver are considered. Here, we introduce a numerical solution to include the true receiver geometry into electromagnetic modelling schemes, which does not rely on such simplifying assumptions. The algorithm is flexible, independent of the chosen numerical method to solve Maxwell's equations and can easily be implemented in other electromagnetic modelling and inversion codes. We present conceptual modelling results for land-based controlled source electromagnetic scenarios and discuss consideration of true receiver geometries for a series of examples of horizontal and vertical electric field measurements. Comparison with Ez data measured in an observation borehole in a producing oil field shows the importance of both considering the true length of the receiver and also its orientation. We show that misalignment from the vertical axis as small as 0.1° may seriously distort the measured signal, as horizontal electric field components are mapped into the desired vertical component. Adequate inclusion of elongated receivers in modelling and inversion can also help reducing effects of static shift when interpreting (natural source) magnetotelluric data.     

Non‐uniqueness with refraction inversion – the Mt Bulga shear zone

The tau‐p inversion algorithm is widely employed to generate starting models with many computer programs that implement refraction tomography. However, this algorithm can frequently fail to detect even major lateral variations in seismic velocities, such as a 50 m wide shear zone, which is the subject of this study. By contrast, the shear zone is successfully defined with the inversion algorithms of the generalized reciprocal method. The shear zone is confirmed with a 2D analysis of the head wave amplitudes, a spectral analysis of the refraction convolution section and with numerous closely spaced orthogonal seismic profiles recorded for a later 3D refraction investigation. Further improvements in resolution, which facilitate the recognition of additional zones with moderate reductions in seismic velocity, are achieved with a novel application of the Hilbert transform to the refractor velocity analysis algorithm. However, the improved resolution also requires the use of a lower average vertical seismic velocity, which accommodates a velocity reversal in the weathering. The lower seismic velocity is derived with the generalized reciprocal method, whereas most refraction tomography programs assume vertical velocity gradients as the default. Although all of the tomograms are consistent with the traveltime data, the resolution of each tomogram is comparable only with that of the starting model. Therefore, it is essential to employ inversion algorithms that can generate detailed starting models, where detailed lateral resolution is the objective. Non‐uniqueness can often be readily resolved with head wave amplitudes, attribute processing of the refraction convolution section and additional seismic traverses, prior to the acquisition of any borehole data. It is concluded that, unless specific measures are taken to address non‐uniqueness, the production of a single refraction tomogram that fits the traveltime data to sufficient accuracy does not necessarily demonstrate that the result is either correct, or even the most probable.     

Advances in the combined interpretation of seismics with magnetotellurics1

Five examples, obtained during exploration for hydrocarbons in the Pannonian Basin of Hungary, are used to show how the interpretation of seismic sections can be usefully complemented by results from MT surveys. Selection of the most appropriate MT quantities, considered to be proper ‘MT attributes’ for the purpose of visualization as well as recognition of the subsurface structures and the different inversions of MT data is essential for practical integration of seismic and MT surveys. A new technique providing a semiquantitive MT-attribute pseudosection for the purpose of visualization of the subsurface structures is proposed. The procedure utilizes derivative functions of the phase of MT impedance for visualization and derives estimated depths from the Bostick transformation of Cagniard apparent resistivities. On the basis of the MT-attribute pseudosections, constructed from the phase derivatives and transformed resistivity data, depths are estimated for interfaces between geological formations with significant resistivity contrast. In particular examples, the interface between the Tertiary sediments and the older basement rocks as well as tectonic fracture zones with decreased resistivity can be resolved.     

基于模型正演的折射静校正方法研究与应用

随着地震勘探区域地质条件复杂化,静校正在资料处理中的作用日显重要.本文模拟研究了折射渡在不同地质条件中的特性,完成了正演模型下的折射波静校正量计算,并将静校正后记录与理论模型记录对比,分析该方法的准确性,总结折射波静校正方法的技术难点,提出一种更准确优越的解决方法.最后,运用该方法选取最合理的参数,并配合先进的处理软件...     

Rock physics modelling and inversion for saturation-pressure changes in time-lapse seismic studies

Time-lapse seismic data are generally used to monitor the changes in dynamic reservoir properties such as fluid saturation and pore or effective pressure. Changes in saturation and pressure due to hydrocarbon production usually cause changes in the seismic velocities and as a consequence changes in seismic amplitudes and travel times. This work proposes a new rock physics model to describe the relation between saturation-pressure changes and seismic changes and a probabilistic workflow to quantify the changes in saturation and pressure from time-lapse seismic changes. In the first part of this work, we propose a new quadratic approximation of the rock physics model. The novelty of the proposed formulation is that the coefficients of the model parameters (i.e. the saturation-pressure changes) are functions of the porosity, initial saturation and initial pressure. The improvements in the results of the forward model are shown through some illustrative examples. In the second part of the work, we present a Bayesian inversion approach for saturation-pressure 4D inversion in which we adopt the new formulation of the rock physics approximation. The inversion results are validated using synthetic pseudo-logs and a 3D reservoir model for CO₂ sequestration.     

2.5D modelling, inversion and angle migration in anisotropic elastic media

2.5D modelling approximates 3D wave propagation in the dip‐direction of a 2D geological model. Attention is restricted to raypaths for waves propagating in a plane. In this way, fast inversion or migration can be performed. For velocity analysis, this reduction of the problem is particularly useful. We review 2.5D modelling for Born volume scattering and Born–Helmholtz surface scattering. The amplitudes are corrected for 3D wave propagation, taking into account both in‐plane and out‐of‐plane geometrical spreading. We also derive some new inversion/migration results. An AVA‐compensated migration routine is presented that is simplified compared with earlier results. This formula can be used to create common‐image gathers for use in velocity analysis by studying the residual moveout. We also give a migration formula for the energy‐flux‐normalized plane‐wave reflection coefficient that models large contrast in the medium parameters not treated by the Born and the Born–Helmholtz equation results. All results are derived using the generalized Radon transform (GRT) directly in the natural coordinate system characterized by scattering angle and migration dip. Consequently, no Jacobians are needed in their calculation. Inversion and migration in an orthorhombic medium or a transversely isotropic (TI) medium with tilted symmetry axis are the lowest symmetries for practical purposes (symmetry axis is in the plane). We give an analysis, using derived methods, of the parameters for these two types of media used in velocity analysis, inversion and migration. The kinematics of the two media involve the same parameters, hence there is no distinction when carrying out velocity analysis. The in‐plane scattering coefficient, used in the inversion and migration, also depends on the same parameters for both media. The out‐of‐plane geometrical spreading, necessary for amplitude‐preserving computations, for the TI medium is dependent on the same parameters that govern in‐plane kinematics. For orthorhombic media, information on additional parameters is required that is not needed for in‐plane kinematics and the scattering coefficients. Resolution analysis of the scattering coefficient suggests that direct inversion by GRT yields unreliable parameter estimates. A more practical approach to inversion is amplitude‐preserving migration followed by AVA analysis. SYMBOLS AND NOTATION A list of symbols and notation is given in Appendix D .     

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.     

Joint inversion of travel times and waveforms by means of annealing algorithms. Application to a seismic refraction experiment

An algorithm of annealing is applied to a joint inversion of travel times and waveforms belonging to a synthetically generated seismic refraction experiment. The medium (crust and upper mantle) is modelled by a set of plane stratified layers and a halfspace. The obtained structure (elastic parameters and depth of layers) shows that, in spite of contaminating seismic noise and poor knowledge of the seismic source, annealing methods are a good tool in these kinds of inversion problems. We think that many characteristics of the annealing process described here could be used with real data and more sophisticated media for the crust and upper mantle of the earth than the present example.