Modelling and migration with orthogonal isochron rays

For increasing time values, isochrons can be regarded as expanding wavefronts and their perpendicular lines as the associated orthogonal isochron rays. The speed of the isochron movement depends on the medium velocity and the source-receiver position. We introduce the term equivalent-velocity to refer to the speed of isochron movement. In the particular case of zero-offset data, the equivalent velocity is half of the medium velocity. We use the concepts of orthogonal isochron-rays and equivalent velocity to extend the application of the exploding reflector model to non-zero offset imaging problems. In particular, we employ these concepts to extend the use of zero-offset wave-equation algorithms for modelling and imaging common-offset sections. In our imaging approach, the common-offset migration is implemented as a trace-by-trace algorithm in three steps: equivalent velocity computation, data conditioning for zero-offset migration and zero-offset wave-equation migration. We apply this methodology for modelling and imaging synthetic common-offset sections using two kinds of algorithms: finite-difference and split-step wavefield extrapolation. We also illustrate the isochron-ray imaging methodology with a field-data example and compare the results with conventional common-offset Kirchhoff migration. This methodology is attractive because it permits depth migration of common-offset sections or just pieces of that by using wave-equation algorithms, it extends the use of robust zero-offset algorithms, it presents favourable features for parallel processing, it permits the creation of hybrid migration algorithms and it is appropriate for migration velocity analysis.     

Numerical tests of 3D true-amplitude zero-offset migration1

An amplitude-preserving migration aims at imaging compressional primary (zero-or) non-zero-offset reflections into 3D time or depth-migrated reflections so that the migrated wavefield amplitudes are a measure of angle-dependent reflection coeffcients. The principal objective is the removal of the geometrical-spreading factor of the primary reflections. Various migration/inversion algorithms involving weighted diffraction stacks proposed recently are based on Born or Kirchhoff approximations. Here, a 3D Kirchhoff-type zero-offset migration approach, also known as a diffraction-stack migration, is implemented in the form of a time migration. The primary reflections of the wavefield to be imaged are described a priori by the zero-order ray approximation. The aim of removing the geometrical- spreading loss can, in the zero-offset case, be achieved by not applying weights to the data before stacking them. This case alone has been implemented in this work. Application of the method to 3D synthetic zero-offset data proves that an amplitude-preserving migration can be performed in this way. Various numerical aspects of the true-amplitude zero-offset migration are discussed.     

COMMON OFFSET PLANE MIGRATION (COPMIG)*

The quality of results of migration before stack is sensitive to inaccuracies in the velocity field applied. This does not hold if only traces of similar sources-receiver distances (common offset traces) enter the migration process. In this case, velocity deviations generate minor shifts in travel times of migrated interfaces but no deterioration in quality. These time shifts are proportional to both the velocity error and the square of the source-receiver distance. The above observations suggest the following migration scheme: migrate separately the traces of the various common offset planes or groups of neighbouring common offset planes; for every common midpoint plane and as a function of travel-time perform a residual NMO search to find trajectories t) =t)_o+px)² of maximum coherency along which migrated events are aligned; correct for residual NMO and stack the migration results obtained in the various common offset planes to obtain the final migration result. This process not only takes care of inaccurate migration velocities but also corrects partly for effects of refraction. It is shown by means of an example that good migration results are generated even with a considerably deviating velocity field.     

INVERSION METHOD IN THE SLANT STACK DOMAIN USING AMPLITUDES OF REFLECTION ARRIVALS*

Almost all ray-tracing methods ignore the analysis of the amplitudes of seismic arrivals and therefore utilize only half of the available information. We propose a method which is a combination of ray-tracing imaging and transformation of the amplitudes of wide-aperture data. Seismic data in the conventional X-T domain are first transformed to the domain of intercept time τ and ray parameter p to recover the plane wave response. The next step is the derivation of a series of plane wave reflection coefficients, which are mapped as a function of τ and p. The reflection coefficients R(τ, p) for two arbitrarily chosen traces can then be used in our inversion method to derive a slowness-depth and a density-depth profile. It is shown that the inclusion of amplitudes of seismic arrivals (in this method, we consider the acoustic case) makes the inverse method highly stable and accurate. In a horizontally stratified medium one can recover separate profiles of velocity and density. Since this method utilizes large-offset data, it can be used for separate recovery of velocity and density to a greater depth.     

Wavefield downward extrapolation for migration velocity analysis on Marmousi data-set

The authors present a method for estimation of interval velocities using the downward continuation of the wavefield to perform layer-stripping migration velocity analysis. The generalized, phase-shift migration MG(F-K) in wavenumber-frequency domain was used for fulltime downward extrapolation of the wavefield. Such downward depth extrapolation accounts for strong changes of velocity in lateral and vertical directions and helps in correct positioning of the wavefield image in complex structures. Determination of velocity is the recursive process which means that the wavefield on depth level z _n−1 (n = 0, 1, ...) is an input data-set for determination of velocity on level z _n . The velocity ν [x, z _n − z _n−1] can be thus treated as interval velocity in Δz _n = z _n − z _n−1 step. This method was tested on synthetic Marmousi data-set and showed satisfactory results for complex, inhomogeneous media.     

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.     

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

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

波场延拓深度滤波方法

消除面波是地震数据处理中的一个重要内容．本文提出了基于15°波动方程的深度滤波方法．由于面波与有效反射波具有不同的传播深度,可利用波场向下延拓方法将二者进行波场分离．把向下延拓后的波场中集中在地表附近的面波能量切除后,再将波场重新延拓回原始的观测面,达到去除干扰的目的．实际资料处理显示：方法计算稳定,消除面波能力强,能更好地保持波场的有效成份和幅值,符合波场的实际传播状态,表明该方法正确可行．     

Combined structure and velocity stacks via the tau–p transform

Reflection and refraction data are normally processed with tools designed to deal specifically with either near- or far-offset data. Furthermore, the refraction data normally require the picking of traveltimes prior to analysis. Here, an automatic processing algorithm has been developed to analyse wide-angle multichannel streamer data without resorting to manual picking or traveltime tomography. Time–offset gathers are transformed to the tau–p domain and the resulting wavefield is downward continued to the depth–p domain from which a velocity model and stacked section are obtained. The algorithm inputs common-depth-point (CDP) gathers and produces a depth-converted stacked section that includes velocity information. The inclusion of long-offset multichannel streamer data within the tau–p transformation enhances the signal from high-velocity refracted basalt arrivals. Downward continuation of the tau–p transformed wavefield to the depth–p domain allows the reflection and refraction components of the wavefield to be treated simultaneously. The high-slowness depth–p wavefield provides the velocity model and the low-slowness depth–p wavefield may be stacked to give structural information. The method is applied to data from the Faeroe Basin from which sub-basalt velocity images are obtained that correlate with an independently derived P-wave model from the line.     

FIELD RESULTS WITH DOWNWARD CONTINUATION TECHNIQUE IN INDUCED POLARIZATION PROFILING USING POINT ELECTRODES1

The method of downward continuation is well known to those working in gravity, magnetic, SP and low-frequency electromagnetic exploration. It is demonstrated that the method of continuation can also be usefully employed in the interpretation of induced polarization gradient profiling using point electrodes to determine target depth. The apparent resistance R_a and chargeability M_a measurements obtained with point electrode excitation of the ground have been used to compute the values of (R_a)_l and (M_a)_l that would be obtained with a linear array. Continuation of the apparent polarizability profile thus obtained with the linear array gives a value for the depth of the target which agrees closely with that obtained by the continuation of the SP profile. On the other hand, continuation of the profile of the secondary transient signal (V_S)_L alone, yields a depth of the target which is in agreement with the borehole information. However, it is seen that the secondary transient voltage profiling response splits into two anomalies which fall on either side of the SP and/or (M_a)_l anomaly centre, and does not coincide with that of the latter.     

Common-reflection-point time migration

Since the early days of seismic processing, time migration has proven to be a valuable tool for a number of imaging purposes. Main motivations for its widespread use include robustness with respect to velocity errors, as well as fast turnaround and low computation costs. In areas of complex geology, in which it has well-known limitations, time migration can still be of value by providing first images and also attributes, which can be of much help in further, more comprehensive depth migration. Time migration is a very close process to common-midpoint (CMP) stacking and, more recently, to zero-offset commonreflection- surface (CRS) stacking. In fact, Kirchhoff time migration operators can be readily formulated in terms of CRS parameters. In the nineties, several studies have shown advantages in the use of common-reflection-point (CRP) traveltimes to replace conventional CMP traveltimes for a number of stacking and migration purposes. In this paper, we follow that trend and introduce a Kirchhoff-type prestack time migration and velocity analysis algorithm, referred to as CRP time migration. The algorithm is based on a CRP operator together with optimal apertures, both computed with the help of CRS parameters. A field-data example indicates the potential of the proposed technique.     

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.     

RADON DEPTH MIGRATION1

A depth migration method is presented that uses Radon-transformed common-source seismograms as input. It is shown that the Radon depth migration method can be extended to spatially varying velocity depth models by using asymptotic ray theory (ART) to construct wavefield continuation operators. These operators downward continue an incident receiver-array plane wave and an assumed point-source wavefield into the subsurface. The migration velocity model is constrained to have longer characteristic wavelengths than the dominant source wavelength such that the ART approximations for the continuation operators are valid. This method is used successfully to migrate two synthetic data examples:

• 1 a point diffractor, and
• 2 a dipping layer and syncline interface model.

It is shown that the Radon migration method has a computational advantage over the standard Kirchhoff migration method in that fewer rays are computed in a main memory implementation.     

A new direction for shallow refraction seismology: integrating amplitudes and traveltimes with the refraction convolution section

The refraction convolution section (RCS) is a new method for imaging shallow seismic refraction data. It is a simple and efficient approach to full‐trace processing which generates a time cross‐section similar to the familiar reflection cross‐section. The RCS advances the interpretation of shallow seismic refraction data through the inclusion of time structure and amplitudes within a single presentation. The RCS is generated by the convolution of forward and reverse shot records. The convolution operation effectively adds the first‐arrival traveltimes of each pair of forward and reverse traces and produces a measure of the depth to the refracting interface in units of time which is equivalent to the time‐depth function of the generalized reciprocal method (GRM). Convolution also multiplies the amplitudes of first‐arrival signals. To a good approximation, this operation compensates for the large effects of geometrical spreading, with the result that the convolved amplitude is essentially proportional to the square of the head coefficient. The signal‐to‐noise (S/N) ratios of the RCS show much less variation than those on the original shot records. The head coefficient is approximately proportional to the ratio of the specific acoustic impedances in the upper layer and in the refractor. The convolved amplitudes or the equivalent shot amplitude products can be useful in resolving ambiguities in the determination of wave speeds. The RCS can also include a separation between each pair of forward and reverse traces in order to accommodate the offset distance in a manner similar to the XY spacing of the GRM. The use of finite XY values improves the resolution of lateral variations in both amplitudes and time‐depths. The use of amplitudes with 3D data effectively improves the spatial resolution of wave speeds by almost an order of magnitude. Amplitudes provide a measure of refractor wave speeds at each detector, whereas the analysis of traveltimes provides a measure over several detectors, commonly a minimum of six. The ratio of amplitudes obtained with different shot azimuths provides a detailed qualitative measure of azimuthal anisotropy and, in turn, of rock fabric. The RCS facilitates the stacking of refraction data in a manner similar to the common‐midpoint methods of reflection seismology. It can significantly improve S/N ratios.Most of the data processing with the RCS, as with the GRM, is carried out in the time domain, rather than in the depth domain. This is a significant advantage because the realities of undetected layers, incomplete sampling of the detected layers and inappropriate sampling in the horizontal rather than the vertical direction result in traveltime data that are neither a complete, an accurate nor a representative portrayal of the wave‐speed stratification. The RCS facilitates the advancement of shallow refraction seismology through the application of current seismic reflection acquisition, processing and interpretation technology.     

Generalized migration in frequency-wavenumber domain (MGF-K) in anisotropic media

In this paper, the background of MGF-K migration in dual domain (wavenumber-frequency K-F and space-time) in anisotropic media is presented. Algorithms for poststack (zero-offset) and prestack migration are based on downward extrapolation of acoustic wavefield by shift-phase with correction filter for lateral variability of medium’s parameters. In anisotropic media, the vertical wavenumber was determined from full elastic wavefield equations for two dimensional (2D) tilted transverse isotropy (TTI) model. The method was tested on a synthetic wavefield for TTI anticlinal model (zero-offset section) and on strongly inhomogeneous vertical transverse isotropy (VTI) Marmousi model. In both cases, the proper imaging of assumed media was obtained.     

A FILTER,DELAY AND SPREAD TECHNIQUE FOR 3D DMO1

The calculation of dip moveout involves spreading the amplitudes of each input trace along the source-receiver axis followed by stacking the results into a 3D zero-offset data cube. The offset-traveltime (x–t) domain integral implementation of the DMO operator is very efficient in terms of computation time but suffers from operator aliasing. The log-stretch approach, using a logarithmic transformation of the time axis to force the DMO operator to be time invariant, can avoid operator aliasing by direct implementation in the frequency-wavenumber (f–k) domain. An alternative technique for log-stretch DMO corrections using the anti-aliasing filters of the f–k approach in the x-log t domain will be presented. Conventionally, the 2D filter representing the DMO operator is designed and applied in the f–k domain. The new technique uses a 2D convolution filter acting in single input/multiple output trace mode. Each single input trace is passed through several 1D filters to create the overall DMO response of that trace. The resulting traces can be stacked directly in the 3D data cube. The single trace filters are the result of a filter design technique reducing the 2D problem to several ID problems. These filters can be decomposed into a pure time-delay and a low-pass filter, representing the kinematic and dynamic behaviour of the DMO operator. The low-pass filters avoid any incidental operator aliasing. Different types of low-pass filters can be used to achieve different amplitude-versus-offset characteristics of the DMO operator.     

Migration methods for imaging different-order multiples

Multiples contain valuable information about the subsurface, and if properly migrated can provide a wider illumination of the subsurface compared to imaging with VSP primary reflections. In this paper we review three different methods for migrating multiples. The first method is model-based, and it is more sensitive to velocity errors than primary migration; the second method uses a semi-natural Green's function for migrating multiples, where part of the traveltimes are computed from the velocity model, and part of the traveltimes (i.e., natural traveltimes) are picked from the data to construct the imaging condition for multiples; the third method uses cross-correlation of traces. The last two methods are preferred in the sense that they are significantly less sensitive to velocity errors and statics because they use “natural data” to construct part of the migration imaging conditions. Compared with the interferometric (i.e., crosscorrelation) imaging method the semi-natural Green's function method is more computationally efficient and is sometimes less prone to migration artifacts. Numerical tests with 2-D and 3-D VSP data show that a wider subsurface coverage, higher-fold and more balanced illumination of the subsurface can be achieved with multiple migration compared with migration of primary reflections only. However, there can be strong interference from multiples with different orders or primaries when multiples of high order are migrated. One possible solution is to filter primaries and different orders of multiples before migration, and another possible solution is least squares migration of all events. A limitation of multiple migration is encountered for subsalt imaging. Here, the multiples must pass through the salt body more than twice, which amplifies the distortion of the image.     

FULL WAVE EQUATION DOWNWARD CONTINUATION OF SEISMIC REFLECTION DATA*

A new method for suppressing multiple reflections in seismograms is developed. It is based on a downward continuation procedure which uses the full acoustic wave equation (hyperbolic form) as a downward continuation operator. We demonstrate that the downward continuation of the recorded wave field maps a reflectivity function without multiply reflected events. The method is applied successfully to individual traces of plane-wave decomposed (slant-stacked) synthetic and field data.     

VELOCITY-STACK PROCESSING1

A conventional velocity-stack gather consists of constant-velocity CMP-stacked traces. It emphasizes the energy associated with the events that follow hyperbolic traveltime trajectories in the CMP gather. Amplitudes along a hyperbola on a CMP gather ideally map onto a point on a velocity-stack gather. Because a CMP gather only includes a cable-length portion of a hyperbolic traveltime trajectory, this mapping is not exact. The finite cable length, discrete sampling along the offset axis and the closeness of hyperbolic summation paths at near-offsets cause smearing of the stacked amplitudes along the velocity axis. Unless this smearing is removed, inverse mapping from velocity space (the plane of stacking velocity versus two-way zero-offset time) back to offset space (the plane of offset versus two-way traveltime) does not reproduce the amplitudes in the original CMP gather. The gather resulting from the inverse mapping can be considered as the model CMP gather that contains only the hyperbolic events from the actual CMP gather. A least-squares minimization of the energy contained in the difference between the actual CMP gather and the model CMP gather removes smearing of amplitudes on the velocity-stack gather and increases velocity resolution. A practical application of this procedure is in separation of multiples from primaries. A method is described to obtain proper velocity-stack gathers with reduced amplitude smearing. The method involves a t²-stretching in the offset space. This stretching maps reflection amplitudes along hyperbolic moveout curves to those along parabolic moveout curves. The CMP gather is Fourier transformed along the stretched axis. Each Fourier component is then used in the least-squares minimization to compute the corresponding Fourier component of the proper velocity-stack gather. Finally, inverse transforming and undoing the stretching yield the proper velocity-stack gather, which can then be inverse mapped back to the offset space. During this inverse mapping, multiples, primaries or all of the hyperbolic events can be modelled. An application of velocity-stack processing to multiple suppression is demonstrated with a field data example.