Interpretation of seismic reflection records: Direct calculation of interval velocities and layer thicknesses from travel times

In reflection surveys and velocity analysis, calculations of interval velocities and layer-thicknesses of a multilayered horizontal structure are often based on Dix's equation which requires the travel times at zero offsets and a prior estimate of the root mean squared velocities.In this paper a method is presented which requires only the reflection travel-time data. A set of equations are derived which relate the interval velocity and thickness of a layer to the reflection travel time from the top and the bottom of that layer, the offset distances and the ray parameter. It is shown that the difference of the offset distances and the difference of the picked travel times of any reflected rays with the same value of ray parameter from the top and the bottom of a horizontal layer can be used to calculate the interval velocity and thickness of that layer.     

Seismic velocity estimation from post-critical wide-angle reflections in layered structures

The travel time inversion of wide-angle seismic data is a technique commonly used in the deep seismic sounding. We propose an application of this technique to a smaller scale of a sedimentary layer, where the characteristics of seismic observations changes significantly. Field observations confirmed by synthetic analysis recognize the dominant amplitudes of wide-angle post-critical reflections. A case study is presented in this paper, of a joint interpretation of conventional reflection seismic with reflection imaging, combined with the wide-angle travel time inversion of additional full-spread observations. A joint interpretation results in a precise recognition of the seismic velocity distribution, that is further used for the seismic depth conversion with the uncertainty analysis of the depth of the reflecting horizons. Despite the salt layer in the studied structure this method is able to precisely recognize the seismic velocities of the sub-salt structures.     

Constraining the anisotropy structure of the crust by joint inversion of seismic reflection travel times and wave polarizations

In the context of wide-angle seismic profiling, the determination of the physical properties of the Earth crust, such as the elastic layer depth and seismic velocity, is often performed by inversion of P- and/or S-phases propagation data supplying the geometry of the medium (reflector depths) or any other structural parameter (P- or S-wave velocity, density...). Moreover, the inversion for velocity structure and interfaces is commonly performed using only seismic reflection travel times and/or crustal phase amplitudes in isotropic media. But it is very important to utilize more available information to constrain the non-uniqueness of the solution. In this paper, we present a simultaneous inversion method of seismic reflection travel times and polarizations data of transient elastic waves in stratified media to reconstruct not only layer depth and vertical P-wave velocity but also the anisotropy feature of the crust based on the estimation of the Thomsen’s parameters. We carry out a checking with synthetic data, comparing the inversion results obtained by anisotropic travel-time inversion to the results derived by joint inversion of seismic reflection travel times and polarizations data. The comparison proves that the first procedure leads to biased anisotropic models, while the second one fits nearly the real model. This makes the joint inversion method feasible. Finally, we investigate the geometry, P-wave velocity structure and anisotropy of the crust beneath Southeastern China by applying the proposed inversion method to previously acquired wide-angle seismic data. In this case, the anisotropy signature provides clear evidence that the Jiangshan-Shaoxing fault is the natural boundary between the Yangtze and Cathaysia blocks.     

Scanning anisotropy parameters in horizontal transversely isotropic media

The horizontal transversely isotropic model, with arbitrary symmetry axis orientation, is the simplest effective representative that explains the azimuthal behaviour of seismic data. Estimating the anisotropy parameters of this model is important in reservoir characterisation, specifically in terms of fracture delineation. We propose a travel‐time‐based approach to estimate the anellipticity parameter η and the symmetry axis azimuth ? of a horizontal transversely isotropic medium, given an inhomogeneous elliptic background model (which might be obtained from velocity analysis and well velocities). This is accomplished through a Taylor's series expansion of the travel‐time solution (of the eikonal equation) as a function of parameter η and azimuth angle ?. The accuracy of the travel time expansion is enhanced by the use of Shanks transform. This results in an accurate approximation of the solution of the non‐linear eikonal equation and provides a mechanism to scan simultaneously for the best fitting effective parameters η and ?, without the need for repetitive modelling of travel times. The analysis of the travel time sensitivity to parameters η and ? reveals that travel times are more sensitive to η than to the symmetry axis azimuth ?. Thus, η is better constrained from travel times than the azimuth. Moreover, the two‐parameter scan in the homogeneous case shows that errors in the background model affect the estimation of η and ? differently. While a gradual increase in errors in the background model leads to increasing errors in η, inaccuracies in ?, on the other hand, depend on the background model errors. We also propose a layer‐stripping method valid for a stack of arbitrary oriented symmetry axis horizontal transversely isotropic layers to convert the effective parameters to the interval layer values.     

Interval velocity and thickness estimate from wide-angle reflection data

A method to estimate interval velocities and thickness in a horizontal isotropic layered medium from wide-angle reflection traveltime curves is presented. The method is based on a relationship between the squared reflection traveltime differences and the squared offset differences relative to two adjacent reflectors. The envelope of the squared-time versus offset-difference curves, for rays with the same ray parameter, is a straight line, whose slope is the inverse of the square of the interval velocity and whose intercept is the square of the interval time. The method yields velocity and thickness estimates without any knowledge of the overlying stratification. It can be applied to wide-angle reflection data when either information on the upper crust and/or refraction control on the velocity is not available. Application to synthetic and real data shows that the method, used together with other methods, allows us to define a reliable 1D starting model for estimating a depth profile using either ray tracing or another technique.     

我国西北地区地壳中的高速夹层

在我国西北地区的柴达木盆地东部和甘肃地区,在距离炮点40互100公里处,能够接收到不少能量较强的地壳深界面反射波。另外还发现一种与一般反射波性质不同的波,其视速度特大,视速度随距离的变化不大,而且有较明显的终点;其吋距曲线与一般深界面反射波的时距曲线相交。根据它的特征可以判断地壳中存在具有速度梯度的高速夹层.求得的夹层参数为: 甘肃地区柴达木盆地东部覆盖层厚度 18.8公里 30.5公里覆盖层平均速度 5.5公里/秒 5.3公里/秒夹层厚度 6.0公里 3.2公里夹层速度 7.5-8.5公里/秒 7.5-8.0公里/秒夹层的上下界面均为强反射面,可以产生多次反射波。分別利用相邻两个反射波可以求得各层参数,并能避免射线折射的影响。甘肃地区和柴达木盆地东部的地壳厚度分別为51和52公里。地壳中有高速夹层的存在,可以更好地说明P~*速度分散的原因,而且也能够解释Lg波的传播机制。     

Subcrustal Low Velocity Layers in Central India and their Implications

A 2-D subcrustal velocity model for the central Indian continental lithosphere has been derived by travel time and relative amplitude modeling of a digitally normalized analog seismic record section of the Hirapur-Mandla DSS profile, using a ray-tracing technique. Some prominent wave groups with apparent velocities slightly higher than the Moho reflection phase (P_MP) are identified on the normalized record sections assembled with a reduction velocity of 6 km s⁻¹. We interpret these phases as the wide-angle reflections from subcrustal lithospheric boundaries. Comparison of synthetic seismograms with the observed record section shows that the observed phases cannot be explained either by multiples or by the P-to-S converted phase (P_MS) from the Moho. Subcrustal velocity models either with a velocity increase or with a single low velocity layer (LVL) also do not provide a satisfactory fit. We infer that a subcrustal velocity model with two alternate LVLs (velocity 7.2 km s⁻¹), separated by a 6-km thick high velocity layer (velocity 8.1 km s⁻¹), can satisfy both the observed travel times and amplitudes. The prominent reflection phases are modeled at depths of 49, 51, 57 and 60 km. It is inferred that the subcrustal lithosphere in the central Indian region has a lamellar structure with varying structural and mechanical properties. The alternating LVLs, occurring at relatively shallow depths below Moho, may be associated with the zones of weakness and lower viscosity suggesting continued mobility, with a possible thermal source in the upper mantle. This explains the source of observed high heat flow values in the central Indian region.     

Interpretation of first arrival travel times in seismic refraction work

The necessary condition for the seismic refraction method to succeed is that the refracted first arrivals from each layer in a multilayered earth system should be detected on a seismogram as first arrivals, and this is possible only when velocities of all underlying layers are successively greater. The usual procedure to interpret the refraction travel times is to fit such a data set with several intersecting straight lines by employing a visual technique which may lead to errors of subjective judgment, as the velocity model depends on the selection of various line segments through the data. To remove the visual fit we propose here a layer stripping method based on minimum intercept time, apparent velocity, rms residual, and maximum data points by least-squares fitting to yield several intersecting straight lines. Once data are segmented out, the conventional equations can be used to determine the velocity structure.     

Validity of the long-wave approximation in periodically layered media

In seismic modelling, a stack of thin layers is often replaced by an effective equivalent anisotropic homogeneous slab. For waves with finite wavelength, this is an approximation, and the error thus introduced can be quantified by considering the relative error in the phase velocity between the layer stack and the effective medium. For periodic layering, the relative phase-velocity error can be expressed in closed form as a function of wavelength, reflection coefficients and layer thicknesses. By comparing the relative phase-velocity error with laboratory measurements and numerical simulations, we find that the difference in seismic response between a periodic layer stack and an equivalent effective medium depends not only on wavelength, but it also depends significantly on reflection coefficients and the ratio between layer thicknesses. For a 1% relative error in the phase velocity, and if all layers have the same thickness measured in vertical traveltime, we find that the wavelength must be larger than approximately three times the layer period for a reflection coefficient of 0.1, but this increases to 13 times the layer period for a reflection coefficient of 0.9, which is highly unrealistic in a geological setting.     

Piecewise‐continuous velocity inversion

We suggest a new method to determine the piecewise‐continuous vertical distribution of instantaneous velocities within sediment layers, using different order time‐domain effective velocities on their top and bottom points. We demonstrate our method using a synthetic model that consists of different compacted sediment layers characterized by monotonously increasing velocity, combined with hard rock layers, such as salt or basalt, characterized by constant fast velocities, and low velocity layers, such as gas pockets. We first show that, by using only the root‐mean‐square velocities and the corresponding vertical travel times (computed from the original instantaneous velocity in depth) as input for a Dix‐type inversion, many different vertical distributions of the instantaneous velocities can be obtained (inverted). Some geological constraints, such as limiting the values of the inverted vertical velocity gradients, should be applied in order to obtain more geologically plausible velocity profiles. In order to limit the non‐uniqueness of the inverted velocities, additional information should be added. We have derived three different inversion solutions that yield the correct instantaneous velocity, avoiding any a priori geological constraints. The additional data at the interface points contain either the average velocities (or depths) or the fourth‐order average velocities, or both. Practically, average velocities can be obtained from nearby wells, whereas the fourth‐order average velocity can be estimated from the quartic moveout term during velocity analysis. Along with the three different types of input, we consider two types of vertical velocity models within each interval: distribution with a constant velocity gradient and an exponential asymptotically bounded velocity model, which is in particular important for modelling thick layers. It has been shown that, in the case of thin intervals, both models lead to similar results. The method allows us to establish the instantaneous velocities at the top and bottom interfaces, where the velocity profile inside the intervals is given by either the linear or the exponential asymptotically bounded velocity models. Since the velocity parameters of each interval are independently inverted, discontinuities of the instantaneous velocity at the interfaces occur naturally. The improved accuracy of the inverted instantaneous velocities is particularly important for accurate time‐to‐depth conversion.     

CALCULATION OF VELOCITY FROM SEISMIC REFLECTION AMPLITUDE *

It is not possible to determine accurate geological velocities from seismic velocity analysis for thin layers or complex structural features, especially under an unconformity. Instead, we can approach the problem of interval velocity with seismic amplitudes analysis and compute the reflection coefficient along the unconformity surface. An error estimation has been made on a model to test the possibility of such a method and to choose the best parameters to be used. The method has been applied on an actual case: the computed interval velocities show good correlation with the values obtained by a sonic log.     

PARAMETER OPTIMIZATION OF VELOCITY DEPTH FUNCTIONS OF GIVEN FORM BY USE OF ROOT-MEAN-SQUARE VELOCITIES*

From seismic surveys zero offset reflection times and root-mean-square velocities are obtained. By use of Dix-Krey's formula, the interval velocities can be calculated. If no well velocity survey exists, the interval velocities and T(o) times are the only available information. The suggested way to get a regionally valid velocity distribution is to select N“leading horizons”, where a major change in the velocity parameters occurs and to compute the parameters of the selected velocity depth function (in most cases linear increase with depth) by a special approximation for the interval between two adjacent “leading horizons”. Herewith all reflection horizons within the interval are taken into account.     

Tutorial: unified 1D inversion of the acoustic reflection response

Acoustic inversion in one-dimension gives impedance as a function of travel time. Inverting the reflection response is a linear problem. Recursive methods, from top to bottom or vice versa, are known and use a fundamental wave field that is computed from the reflection response. An integral over the solution to the Marchenko equation, on the other hand, retrieves the impedance at any vertical travel time instant. It is a non-recursive method, but requires the zero-frequency value of the reflection response. These methods use the same fundamental wave field in different ways. Combining the two methods leads to a non-recursive scheme that works with finite-frequency bandwidth. This can be used for target-oriented inversion. When a reflection response is available along a line over a horizontally layered medium, the thickness and wave velocity of any layer can be obtained together with the velocity of an adjacent layer and the density ratio of the two layers. Statistical analysis over 1000 noise realizations shows that the forward recursive method and the Marchenko-type method perform well on computed noisy data.     

Normal moveout in focus1

Improving the accuracy of NMO corrections and of the corresponding interval velocities entails implementing a better approximation than the formula used since the beginning of seismic processing. The exact equations are not practical as they include many unknowns. The approximate expression has only two unknowns, the reflection time and the rms velocity, but becomes inaccurate for large apertures of the recording system and heterogeneous vertical velocities. Several methods of improving the accuracy have been considered, but the gains do not compensate for the dramatic increase in computing time. Two alternative equations are proposed: the first containing two parameters, the reflection time and the focusing time, is not valid for apertures much greater than is the standard formula, but has a much faster computing time and does not stretch the far traces; the other, containing three parameters, the reflection time, like focusing time and the tuning velocity, retains high frequencies for apertures about twice those allowed by the standard equation. Its computing time can be kept within the same limits. NMO equations, old and new, are designed strictly for horizontal layering, but remain reliable as long as the rays travel through the same layers in both the down and up directions. An equation, similar to Dix's formula, is given to compute the interval velocities. The entire scheme can be automated to produce interval-velocity sections without manual picking.     

The amplitude ratioPcP/P and the core-mantle boundary

Summary Using the Haskell matrix formulation, theoretical reflection coefficient curves have been calculated for a multi-layered core-mantle boundary for comparison with observational data. Two cases are considered, first when the shear velocity in the core is equal to zero and second when the core has a finite rigidity. If the velocity contrast is large between the imbedded layer and the mantle, the reflection coefficient curves for the multi-layered medium are irregular in shape as compared to those for two half-spaces, representing the core and the mantle, respectively. The reflection coefficient curves show an oscillatory character if the imbedded layer is thick and has a high velocity contrast.The observational data consist of short-period vertical-component seismograph records ofP andPcP from nuclear explosions in the Aleutian chain, Nevada, Novaya Zemlya, Kazakh and Sahara. Attenuation and geometrical spreading are taken into consideration. Four different models for the quality factorQ are applied to the observational data. The data are found to be much affected by theQ-model used for the corrections.Based on proposedQ-values, a model for the core-mantle boundary is found, characterized by two low-velocity layers at the bottom of the mantle. The thicknesses are 16.10 km (outer layer) and 19.96 km (inner layer), the compressional wave velocities 12.17 km/sec and 10.94 km/sec and the shear wave velocities are 6.29 km/sec and 5.33 km/sec, respectively. A better fit to this model is found when in addition the shear velocity in the outer core is 2.20 km/sec and the density ratio at the core-mantle boundary is 1.07. In other words, the observations favour a layer of finite rigidity in the outer core rather than a fluid one.     

NON-LINEAR AVO INVERSION FOR A STACK OF ANELASTIC LAYERS1

Parameters in a stack of homogeneous anelastic layers are estimated from seismic data, using the amplitude versus offset (AVO) variations and the travel-times. The unknown parameters in each layer are the layer thickness, the P-wave velocity, the S-wave velocity, the density and the quality factor. Dynamic ray tracing is used to solve the forward problem. Multiple reflections are included, but wave-mode conversions are not considered. The S-wave velocities are estimated from the PP reflection and transmission coefficients. The inverse problem is solved using a stabilized least-squares procedure. The Gauss-Newton approximation to the Hessian matrix is used, and the derivatives of the dynamic ray-tracing equation are calculated analytically for each iteration. A conventional velocity analysis, the common mid-point (CMP) stack and a set of CMP gathers are used to identify the number of layers and to establish initial estimates for the P-wave velocities and the layer thicknesses. The inversion is carried out globally for all parameters simultaneously or by a stepwise approach where a smaller number of parameters is considered in each step. We discuss several practical problems related to inversion of real data. The performance of the algorithm is tested on one synthetic and two real data sets. For the real data inversion, we explained up to 90% of the energy in the data. However, the reliability of the parameter estimates must at this stage be considered as uncertain.     

宽角反射地震波走时模拟的双重网格法

在研究地壳结构的人工源宽角反射地震资料解释中,常规宽角反射波走时和射线路径计算大都假定地壳模型为层状块状均匀介质.为了逼近实际地壳结构模型,要求模型尺度较大,为了提高地震资料解释的可靠性,须减小模型离散单元的尺寸,但同时计算量大大增加,使资料解释的效率较低.为此,本文尝试同时提高宽角反射地震资料解释效率和可靠性的方法,即使用双重网格计算宽角反射地震波走时和射线路径的最小走时树方法.双重网格法在均匀介质内部仅计算大网格节点,在速度变化点、震源点和检波点区域,同时计算小网格节点;在界面边界点使用比介质内部节点更大的子波传播区域.模型计算结果表明,对于大尺度的层状块状均匀介质模型,在保证精度的条件下,本文所提出的双重网格射线追踪方法的计算效率比单网格方法显著提高.     

THREE-DIMENSIONAL SEISMIC MODELING FOR ARBITRARY VELOCITY DISTRIBUTIONS *

The growing awareness of the significance of three-dimensional interpretation of seismic reflection data carries with it a need for better understanding of the role of the velocity configuration in the production of reflection patterns. Modeling the reflection responses of a simple dipping plane reflector through velocity models of various degrees of complexity demonstrates the importance of overlying velocities in determining the reflection pattern. Modeling is accomplished using a raytracing technique which determines total travel time of the normal incidence raypath through an arbitrary iso-velocity layer model.     

Moveout approximation for horizontal transversely isotropic and vertical transversely isotropic layered medium. Part II: effective model‡

We use residual moveouts measured along continuous full azimuth reflection angle gathers, in order to obtain effective horizontal transversely isotropic model parameters. The angle gathers are generated through a special angle domain imaging system, for a wide range of reflection angles and full range of phase velocity azimuths. The estimation of the effective model parameters is performed in two stages. First, the background horizontal transversely isotropic (HTI)/vertical transversely isotropic (VTI) layered model is used, along with the values of reflection angles, for converting the measured residual moveouts (or traveltime errors) into azimuthally dependent normal moveout (NMO) velocities. Then we apply a digital Fourier transform to convert the NMO velocities into azimuthal wavenumber domain, in order to obtain the effective HTI model parameters: vertical time, vertical compression velocity, Thomsen parameter delta and the azimuth of the medium axis of symmetry. The method also provides a reliability criterion of the HTI assumption. The criterion shows whether the medium possesses the HTI type of symmetry, or whether the azimuthal dependence of the residual traveltime indicates to a more complex azimuthal anisotropy. The effective model used in this approach is defined for a 1D structure with a set of HTI, VTI and isotropic layers (with at least one HTI layer). We describe and analyse the reduction of a multi‐layer structure into an equivalent effective HTI model. The equivalent model yields the same NMO velocity and the same offset azimuth on the Earth's surface as the original layered structure, for any azimuth of the phase velocity. The effective model approximates the kinematics of an HTI/VTI layered structure using only a few parameters. Under the hyperbolic approximation, the proposed effective model is exact.     

LIMITS OF STRUCTURAL INVERSION OF SEISMIC HORIZONS*

Time horizons can be depth-migrated when interval velocities are known; on the other hand, the velocity distribution can be found when traveltimes and NMO velocities at zero offset are known (wavefront curvatures; Shah 1973). Using these concepts, exact recursive inversion formulae for the calculation of interval velocities are given. The assumption of rectilinear raypath propagation within each layer is made; interval velocities and curvatures of the interfaces between layers can be found if traveltimes together with their gradients and curvatures and very precise V_NMO velocities at zero offset are known. However, the available stacking velocity is a numerical quantity which has no direct physical significance; its deviation from zero offset NMO velocity is examined in terms of horizon curvatures, cable length and lateral velocity inhomogeneities. A method has been derived to estimate the geological depth model by searching, iteratively, for the best solution that minimizes the difference between stacking velocities from the real data and from the structural model. Results show the limits and capabilities of the approach; perhaps, owing to the low resolution of conventional velocity analyses, a simplified version of the given formulae would be more robust.