P‐wave stacking‐velocity tomography for VTI media

A major complication caused by anisotropy in velocity analysis and imaging is the uncertainty in estimating the vertical velocity and depth scale of the model from surface data. For laterally homogeneous VTI (transversely isotropic with a vertical symmetry axis) media above the target reflector, P‐wave moveout has to be combined with other information (e.g. borehole data or converted waves) to build velocity models for depth imaging. The presence of lateral heterogeneity in the overburden creates the dependence of P‐wave reflection data on all three relevant parameters (the vertical velocity V_P0 and the Thomsen coefficients ε and δ) and, therefore, may help to determine the depth scale of the velocity field. Here, we propose a tomographic algorithm designed to invert NMO ellipses (obtained from azimuthally varying stacking velocities) and zero‐offset traveltimes of P‐waves for the parameters of homogeneous VTI layers separated by either plane dipping or curved interfaces. For plane non‐intersecting layer boundaries, the interval parameters cannot be recovered from P‐wave moveout in a unique way. Nonetheless, if the reflectors have sufficiently different azimuths, a priori knowledge of any single interval parameter makes it possible to reconstruct the whole model in depth. For example, the parameter estimation becomes unique if the subsurface layer is known to be isotropic. In the case of 2D inversion on the dip line of co‐orientated reflectors, it is necessary to specify one parameter (e.g. the vertical velocity) per layer. Despite the higher complexity of models with curved interfaces, the increased angle coverage of reflected rays helps to resolve the trade‐offs between the medium parameters. Singular value decomposition (SVD) shows that in the presence of sufficient interface curvature all parameters needed for anisotropic depth processing can be obtained solely from conventional‐spread P‐wave moveout. By performing tests on noise‐contaminated data we demonstrate that the tomographic inversion procedure reconstructs both the interfaces and the VTI parameters with high accuracy. Both SVD analysis and moveout inversion are implemented using an efficient modelling technique based on the theory of NMO‐velocity surfaces generalized for wave propagation through curved interfaces.     

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.     

Seismic traveltime inversion of 3D velocity model with triangulated interfaces

Seismic traveltime tomographic inversion has played an important role in detecting the internal structure of the solid earth. We use a set of blocks to approximate geologically complex media that cannot be well described by layered models or cells. The geological body is described as an aggregate of arbitrarily shaped blocks, which are separated by triangulated interfaces. We can describe the media as homogenous or heterogeneous in each block. We define the velocities at the given rectangle grid points for each block, and the heterogeneous velocities in each block can be calculated by a linear interpolation algorithm. The parameters of the velocity grid positions are independent of the model parameterization, which is advantageous in the joint inversion of the velocities and the node depths of an interface. We implement a segmentally iterative ray tracer to calculate traveltimes in the 3D heterogeneous block models. The damped least squares method is employed in seismic traveltime inversion, which includes the partial derivatives of traveltime with respect to the depths of nodes in the triangulated interfaces and velocities defined in rectangular grids. The numerical tests indicate that the node depths of a triangulated interface and homogeneous velocity distributions can be well inverted in a stratified model.     

METHOD FOR DETERMINATION OF VELOCITY AND DEPTH FROM SEISMIC REFLECTION DATA1

The estimation of velocity and depth is an important stage in seismic data processing and interpretation. We present a method for velocity-depth model estimation from unstacked data. This method is formulated as an iterative algorithm producing a model which maximizes some measure of coherency computed along traveltimes generated by tracing rays through the model. In the model the interfaces are represented as cubic splines and it is assumed that the velocity in each layer is constant. The inversion includes the determination of the velocities in all the layers and the location of the spline knots. The process input consists of unstacked seismic data and an initial velocity-depth model. This model is often based on nearby well information and an interpretation of the stacked section. Inversion is performed iteratively layer after layer; during each iteration synthetic travel-time curves are calculated for the interface under consideration. A functional characterizing the main correlation properties of the wavefield is then formed along the synthetic arrival times. It is assumed that the functional reaches a maximum value when the synthetic arrival time curves match the arrival times of the events on the field gathers. The maximum value of the functional is obtained by an effective algorithm of non-linear programming. The present inversion algorithm has the advantages that event picking on the unstacked data is not required and is not based on curve fitting of hyperbolic approximations of the arrival times. The method has been successfully applied to both synthetic and field data.     

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.     

利用地震回折波资料反演界面位置与速度分布

论述了利用地震回折波资料反演界面位置与速度分布的方法，推导了地震波走时对于界面位置偏导数的计算公式。数值模拟和实测资料的计算结果表明了该方法的有效性和编制的计算程序的实用性。该方法最突出的特点是充分地利用了透射波资料中所含的界面位置的信息。界面位置的分辨率与界面两边的速度反差有关，速度差别越大，则分辨率越高。     

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.     

STOCHASTIC INVERSION BY RAY CONTINUATION: APPLICATION TO SEISMIC TOMOGRAPHY1

The conventional tomographic inversion consists in minimizing residuals between measured and modelled traveltimes. The process tends to be unstable and some additional constraints are required to stabilize it. The stochastic formulation generalizes the technique and sets it on firmer theoretical bases. The Stochastic Inversion by Ray Continuation (Sirc ) is a probabilistic approach, which takes a priori geological information into account and uses probability distributions to characterize data correlations and errors. It makes it possible to tie uncertainties to the results. The estimated parameters are interval velocities and B -spline coefficients used to represent smoothed interfaces. Ray tracing is done by a continuation technique between source and receivers. The ray coordinates are computed from one path to the next by solving a linear system derived from Fermat's principle. The main advantages are fast computations, accurate traveltimes and derivatives. The seismic traces are gathered in CMPs. For a particular CMP, several reflecting elements are characterized by their time gradient measured on the stacked section, and related to a mean emergence direction. The program capabilities are tested on a synthetic example as well as on a field example. The strategy consists in inverting the parameters for one layer, then for the next one down. An inversion step is divided in two parts. First the parameters for the layer concerned are inverted, while the parameters for the upper layers remain fixed. Then all the parameters are reinverted. The velocity-depth section computed by the program together with the corresponding errors can be used directly for the interpretation, as an initial model for depth migration or for the complete inversion program under development.     

COMPUTATION OF INTERVAL VELOCITIES FROM COMMON REFLECTION POINT MOVEOUT TIMES FOR n LAYERS WITH ARBITRARY DIPS AND CURVATURES IN THREE DIMENSIONS WHEN ASSUMING SMALL SHOT-GEOPHONE DISTANCES*

It is well known that interval velocities can be determined from common-reflection-point moveout times. However, the mathematics becomes complicated in the general case of n homogeneous layers with curved interfaces dipping in three dimensions. In this paper the problem is solved by mathematical induction using the second power terms only of the Taylor series which represents the moveout time as a function of the coordinate differences between shot and geophone points. Moreover, the zero-offset reflection times of the nth interface in a certain area surrounding the point of interest have to be known. The n—I upper interfaces and interval velocities are known too on account of the mathematical induction method applied. Thus, the zero-offset reflection raypath of the nth interface can be supposed to be known down to the intersection with the (n—1)th interface. The method applied consists mainly in transforming the second power terms of the moveout time from one interface to the next one. This is accomplished by matrix algebra. Some special cases are discussed as e.g. uniform strike and small curvatures.     

OBTAINING INTERVAL VELOCITIES FROM STACKING VELOCITIES WHEN DIPPING HORIZONS ARE INCLUDED *

Common-depth-point stacking velocities may differ from root-mean-square velocities because of large offset and because of dipping reflectors. This paper shows that the two effects may be treated separately, and proceeds to examine the effect of dip. If stacking velocities are assumed equal to rms velocities for the purpose of time to depth conversion, then errors are introduced comparable to the difference between migrated and unmigrated depths. Consequently, if the effect of dip on stacking velocity is ignored, there is no point in migrating the resulting depth data. For a multi-layered model having parallel dip, a formula is developed to compute interval velocities and depths from the stacking velocities, time picks, and time slope of the seismic section. It is shown that cross-dip need not be considered, if all the reflectors have the same dip azimuth. The problem becomes intractable if the dips are not parallel. But the inverse problem is soluble: to obtain, stacking velocities; time picks, and time slopes from a given depth and interval velocity model. Finally, the inverse solution is combined with an approximate forward solution. This provides an iterative method to obtain depths and interval velocities from stacking velocities, time picks and time slopes. It is assumed that the dip azimuth is the same for all reflectors, but not necessarily in the plane of the section, and that the curvature of the reflecting horizons is negligible. The effect of onset delay is examined. It is shown that onset corrections may be unnecessary when converting from time to depth.     

A STABLE AND FLEXIBLE PROCEDURE FOR THE INVERSE MODELLING OF SEISMIC FIRST ARRIVALS1

Multiple coverage reflection seismic data provide an important source of information concerning the subsurface. However, due to the stacking and migration techniques used in the processing, the first arrivals are muted and details about the upper part of the sections are generally lost. This paper describes a computerized method for the inverse modelling of laterally varying velocities and shallow depths which are not sufficiently resolved in the reflection seismic processing. The method minimizes, in a least-squares manner, the difference between the observed first arrivals, picked from the reflection traces, and a set of synthetic traveltimes, calculated by ray tracing in a cell model. An initial model, e.g. from a priori knowledge or the application of a conventional interpretation method, is refined iteratively until no further essential improvement can be achieved. Traditional first-arrival inversion methods cannot, in general, provide such flexible modelling. The technique is successfully tested on synthetic data as well as on first arrivals picked automatically from the records of a reflection seismic survey in North Jutland, Denmark.     

TIME MIGRATION—SOME RAY THEORETICAL ASPECTS*

Using an elementary theory of migration one can consider a reflecting horizon as a continuum of scattering centres for seismic waves. Reflections arising at interfaces can thus be looked upon as the sum of energy scattered by interface points. The energy from one point is distributed among signals upon its reflection time surface. This surface is usually well approximated by a hyperboloid in the vicinity of its apex. Migration aims at focusing the scattered energy of each depth point into an image point upon the reflection time surface. To ensure a complete migration the image must be vertical above the depth point. This is difficult to achieve for subsurface interfaces which fall below laterally in-homogeneous velocity media. Migration is hence frequently performed for these interfaces as well by the Kirchhoff summation method which systematically sums signals into the apex of the approximation hyperboloid even though the Kirchhoff integral is in this case not strictly valid. For a multilayered subsurface isovelocity layer model with interfaces of a generally curved nature this can only provide a complete migration for the uppermost interface. Still there are various advantages gained by having a process which sums signals consistently into the minimum of the reflection time surface. The position of the time surface minimum is the place where a ray from the depth point emerges vertically to the surface. The Kirchhoff migration, if applied to media with laterally inhomogeneous velocity, must necessarily be followed by a further time-to-depth migration if the true depth structure is to be recovered. Primary normal reflections and their respective migrated reflections have a complementary relationship to each other. Normal reflections relate to rays normal to the reflector and migrated reflections relate to rays normal to the free surface. Ray modeling is performed to indicate a new approach for simulating seismic reflections. Commonly occuring situations are investigated from which lessons can be learned which are of immediate value for those concerned with interpreting time migrated reflections. The concept of the ‘image ray’ is introduced.     

Ray path of head waves with irregular interfaces

Head waves are usually considered to be the refracted waves propagating along flat interfaces with an underlying higher velocity. However, the path that the rays travel along in media with irregular interfaces is not clear. Here we study the problem by simulation using a new approach of the spectral-element method with some overlapped elements (SEMO) that can accurately evaluate waves traveling along an irregular interface. Consequently, the head waves are separated from interface waves by a time window. Thus, their energy and arrival time changes can be analyzed independently. These analyses demonstrate that, contrary to the case for head waves propagating along a flat interface, there are two mechanisms for head waves traveling along an irregular interface: a refraction mechanism and transmission mechanism. That is, the head waves may be refracted waves propagating along the interface or transmitted waves induced by the waves propagating in the higher-velocity media. Such knowledge will be helpful in constructing a more accurate inversion method, such as head wave travel-time tomography, and in obtaining a more accurate model of subsurface structure which is very important for understanding the formation mechanism of some special areas, such as the Tibetan Plateau.     

PARABOLIC AND HYPERBOLIC PARAXIAL TWO-POINT TRAVELTIMES IN 3D MEDIA1

The 4 × 4 T -propagator matrix of a 3D central ray determines, among other important seismic quantities, second-order (parabolic or hyperbolic) two-point traveltime approximations of certain paraxial rays in the vicinity of the known central ray through a 3D medium consisting of inhomogeneous isotropic velocity layers. These rays result from perturbing the start and endpoints of the central ray on smoothly curved anterior and posterior surfaces. The perturbation of each ray endpoint is described only by a two-component vector. Here, we provide parabolic and hyperbolic paraxial two-point traveltime approximations using the T -propagator to feature a number of useful 3D seismic models, putting particular emphasis on expressing the traveltimes for paraxial primary reflected rays in terms of hyperbolic approximations. These are of use in solving several forward and inverse seismic problems. Our results simplify those in which the perturbation of the ray endpoints upon a curved interface is described by a three-component vector. In order to emphasize the importance of the hyperbolic expression, we show that the hyperbolic paraxial-ray traveltime (in terms of four independent variables) is exact for the case of a primary ray reflected from a planar dipping interface below a homogeneous velocity medium.     

界面二次源波前扩展法全局最小走时射线追踪技术

以Moser方法为代表的最短路径射线追踪算法可以快速稳定地获得整个追踪区域的全局最小走时和路径,但它存在两个缺陷：一是射线大多由折线呈锯齿状相连,长度和位置偏离真实射线路径;二是在低变速区容易出现射线路径多值现象．本文提出的界面二次源波前扩展法全局最小走时射线追踪技术（以下简称界面源法）旨在解决上述两个问题．不同于Moser方法,界面源法只在物性分界面上设置子波源点,子波出射射线可以到达任何不穿越物性界面而直接到达的空间点和界面离散点,在均匀块体内或层内地震波以精确的射线路径传播．显然,界面源法的子波出射方向数远远大于传统方法,算法的追踪误差主要由界面离散引起的,因此,界面源法很好地解决了Moser法存在的问题,大大提高了追踪的精度．同时,由于界面源法的子波源点数远远小于Moser法,因而效率也很高．模型实算证实了该算法的高效性．     

Improving modelling and inversion in refraction seismics with a first-order Eikonal solver

A first-order Eikonal solver is applied to modelling and inversion in refraction seismics. The method calculates the traveltime of the fastest wave at any point of a regular grid, including head waves as used in refraction. The efficiency, robustness and flexibility of the method give a very powerful modelling tool to find both traveltimes and raypaths. Comparisons with finite-difference data show the validity of the results. Any arbitrarily complex model can be studied, including the exact topography of the surface, thus avoiding static corrections. Later arrivals are also obtained by applying high-slowness masks over the high-velocity zones. Such an efficient modelling tool may be used interactively to invert for the model, but a better method is to apply the refractor-imaging principle of Hagedoorn to obtain the refractors from the picked traveltime curves. The application of this principle has already been tried successfully by previous authors, but they used a less well-adapted Eikonal solver. Some of their traveltimes were not correct in the presence of strong velocity variations, and the refractor-imaging principle was restricted to receiver lines along a plane surface. With the first-order Eikonal solver chosen, any topography of the receiving surface can be considered and there is no restriction on the velocity contrast. Based on synthetic examples, the Hagedoorn principle appears to be robust even in the case of first arrivals associated with waves diving under the refractor. The velocities below the refractor can also be easily estimated, parallel to the imaging process. In this way, the model can be built up successively layer by layer, the refractor-imaging and velocity-mapping processes being performed for each identified refractor at a time. The inverted model could then be used in tomographic inversions because the calculated traveltimes are very close to the observed traveltimes and the raypaths are available.     

Review of the Anisotropic Interface Ray Propagator: Symplecticity,Eigenvalues, Invariants and Applications

A review of the 6 × 6 anisotropic interface ray propagator matrix in Cartesian coordinates and within the framework of the Hamiltonian formalism shows that there is one unique propagator satisfying the symplectic property. This is essential, since the symplecticity furnishes an exact inverse, while an eigenvalue analysis indicates that the propagator may be arbitrarily ill-conditioned. As such, the symplectic interface propagator naturally connects to symplectic ray integration algorithms for smooth media, designed to maintain accuracy. Moreover, several ray invariants for smooth media remain invariant across interfaces. It is straightforward to derive expressions for the interface propagator, both explicit and implicit. Symplecticity is equivalent to the condition that the propagator preserves the eikonal constraint across the interface. The symplectic interface propagator complies with phase matching of the incident and reflected/transmitted ray field, and is therefore in accordance with the earlier derived 4 × 4 matrix in ray-centred coordinates. The symplectic property is related to the symmetry of the second derivative matrix of the reflected/transmitted traveltime field. Thanks to the analytic expression of the symplectic interface propagator, relating interface curvature directly to second derivatives of traveltimes observed at a datum level, numerous applications are available in the area of processing and inversion.     

Inversion of travel times obtained during active seismic refraction experiments CELEBRATION 2000, ALP 2002 and SUDETES 2003

A series of kinematic inversions based on robust non-linear optimization approach were performed using travel time data from a series of seismic refraction experiments: CELEBRATION 2000, ALP 2002 and SUDETES 2003. These experiments were performed in Central Europe from 2000 to 2003. Data from 8 profiles (CEL09, CEL10, Alp01, S01, S02, S03, S04 and S05) were processed in this study. The goal of this work was to find seismic velocity models yielding travel times consistent with observed data. Optimum 2D inhomogeneous isotropic P-wave velocity models were computed. We have developed and used a specialized two-step inverse procedure. In the first “parametric” step, the velocity model contains interfaces whose shapes are defined by a number of parameters. The velocity along each interface is supposed to be constant but may be different along the upper and lower side of the interface. Linear vertical interpolation is used for points in between interfaces. All parameters are searched for using robust non-linear optimization (Differential Evolution algorithm). Rays are continuously traced by the bending technique. In the second “tomographic” step, small-scale velocity perturbations are introduced in a dense grid covering the currently obtained velocity model. Rays are fixed in this step. Final velocity models yield travel time residuals comparable to typical picking errors (RMS ∼ 0.1 s). As a result, depth-velocity cross-sections of P waves along all processed profiles are obtained. The depth range of the models is 35–50 km, the velocity varies in the range 3.5–8.2 km/s. Lowest velocities are detected in near-surface depth sections crossing sedimentary formations. The middle crust is generally more homogeneous and has typical P wave velocity around 6 km/s. Surprisingly the lower crust is less homogeneous and the computed velocity is in the range 6.5–7.5 km/s. The MOHO is detected in the depth ≈30–45 km.     

Multi-refractor imaging with stacked refraction convolution section

Multi‐refractor imaging is a technique for constructing a single two‐dimensional image of a number of refractors by stacking multiple convolved and cross‐correlated reversed shot records. The method is most effective with high‐fold data that have been obtained with roll‐along acquisition programs because the stacking process significantly improves the signal‐to‐noise ratios. The major advantage of the multi‐refractor imaging method is that all the data can be stacked to maximize the signal‐to‐noise ratios before the measurement of any traveltimes. However, the signal‐to‐noise ratios can be further increased if only those traces that have arrivals from the same refractor are used, and if the correct reciprocal times or traces are employed. A field case study shows that multi‐refractor imaging can produce a cross‐section similar to the familiar reflection cross‐section with substantially higher signal‐to‐noise ratios for the equivalent interfaces.     

COMPUTATION OF THE NORMAL MOVEOUT VELOCITY IN 3D LATERALLY INHOMOGENEOUS MEDIA WITH CURVED INTERFACES*

For a 3D velocity model of curved first order interfaces and layer velocities which are arbitrary smooth functions of the space coordinates, the normal moveout (NMO)-velocity can be computed by numerically integrating a system of first order ordinary differential equations for a hypothetical wavefront that originates at the normal incidence point of the normal ray and moves up along the ray to the common mid-point of the common datum point (CDP) profile.