THE CAGNIARD METHOD IN COMPLEX TIME REVISITED1

The Cagniard-de Hoop method is ideally suited to the analysis of wave propagation problems in stratified media. The method applies to the integral transform representation of the solution in the transform variables (s, p) dual of the time and transverse distance. The objective of the method is to make the p-integral take the form of a forward Laplace transform, so that the cascade of the two integrals can be identified as a forward and inverse transform, thereby making the actual integration unnecessary. That is, the exponent (–sw(p)) is set equal to –sτ, with τ varying from some (real) finite time to infinity. As usually presented, the p-integral is deformed onto a contour on which the exponent is real and decreases to –∞ as p tends to infinity. We have found that it is often easier to introduce a complex variable τ for the exponent and carry out the deformation of contour in the complex τ-domain. In the τ-domain the deformation amounts to ‘closing down’ the contour of integration around the real axis while taking due account of singularities off this axis. Typically, the method is applied to an integral that represents one body wave plus other types of waves. In this approach, the saddle point of w(p) that produces the body wave plays a crucial role because it is always a branch point of the integrand in the τ-domain integral. Furthermore, the paths of steepest ascent from the saddle point are always the tails of the Cagniard path along which w(p) →∞. That is, the image of the pair of steepest ascent paths in the p-domain is a double covering of a segment of the Re τ-axis in the τ-domain. The deformed contour in the p-domain will be the only pair of steepest ascent paths unless the original integrand had other singularities in the p-domain between the imaginary axis and this pair of contours. This motivates the definition of a primary p-domain, i.e. the domain between the imaginary axis and the steepest descent paths, and its image in the τ-domain, the primary τ-domain. In terms of these regions, singularities in the primary p-domain have images in the primary τ-domain and the deformation of contour on to the real axis in the τ-domain must include contributions from these singularities. This approach to the Cagniard-de Hoop method represents a return from de Hoop's modification to Cagniard's original method, but with simplifications that make the original method more tractable and straightforward. This approach is also reminiscent of van der Waerden's approach to the method of steepest descents, which starts exactly the same way. Indeed, after the deformation of contour in the τ-domain, the user has the choice of applying asymptotic analysis to the resulting ‘loop’ integrals (Watson's lemma) or continuing to obtain the exact, although usually implicit, time-domain solution by completing the Cagniard-de Hoop analysis. In developing the method we examine the transformation from a frequency-domain representation of the solution (ω) to a Laplace representation (s). Many users start from the frequency-domain representation of solutions of wave propagation problems. In this case issues arising from the movement of singularities under the transformation from ω to s must be considered. We discuss this extension in the context of the Sommerfeld half-plane problem.     

FREQUENCY WAVENUMBER APPROACH OF THE τ-p TRANSFORM: SOME APPLICATIONS IN SEISMIC DATA PROCESSING*

The τ-p transform is an invertible transformation of seismic shot records expressed as a function of time and offset into the τ (intercept time) and p (ray parameter) domain. The τ-p transform is derived from the solution of the wave equation for a point source in a three-dimensional, vertically non-homogeneous medium and therefore is a true amplitude process for the assumed model. The main advantage of this transformation is to present a point source shot record as a series of plane wave experiments. The asymptotic expansion of this transformation is found to be useful in reflection seismic data processing. The τ-p and frequency-wavenumber (or f-k) processes are closely related. Indeed, the τ-p process embodies the frequency-wavenumber transformation, so the use of this technique suffers the same limitations as the f-k technique. In particular, the wavefield must be sampled with sufficient spatial density to avoid wavenumber aliasing. The computation of this transform and its inverse transform consists of a two-dimensional Fast Fourier Transform followed by an interpolation, then by an inverse-time Fast Fourier Transform. This technique is extended from a vertically inhomogeneous three-dimensional medium to a vertically and laterally inhomogeneous three-dimensional medium. The τ-p transform may create artifacts (truncation and aliasing effects) which can be reduced by a finer spatial density of geophone groups by a balancing of the seismic data and by a tapering of the extremities of the seismic data. The τ-p domain is used as a temporary domain where the attack of coherent noise is well addressed; this technique can be viewed as ‘time-variant f-k filtering’. In addition, the process of deconvolution and multiple suppression in the τ-p domain is at least as well addressed as in the time-offset domain.     

THREE-TERM TAYLOR SERIES FOR t2 - x2-CURVES OF P- AND S-WAVES OVER LAYERED TRANSVERSELY ISOTROPIC GROUND*

The arrival-time curve of a reflection from a horizontal interface, beneath a homogeneous isotropic layer, is a hyperbola in the x - t-domain. If the subsurface is one-dimensionally inhomogeneous (horizontally layered), or if some or all of the layers are transversely isotropic with vertical axis of symmetry, the statement is no longer strictly true, though the arrival-time curves are still hyperbola-like. In the case of transverse isotropy, however, classical interpretation of these curves fails. Interval velocities calculated from t² - x²-curves do not always approximate vertical velocities and therefore cannot be used to calculate depths of reflectors. To study the relationship between velocities calculated from t² - x²-curves and the true velocities of a transversely isotropic layer, we approximate t² - x²-curves over a vertically inhomogeneous transversely isotropic medium by a three-term Taylor series and calculate expressions for these terms as a function of the elastic parameters. It is shown that both inhomogeneity and transverse isotropy affect slope and curvature of t² - x²-curves. For P-waves the effect of transverse isotropy is that the t² - x²-curves are convex upwards; for SV-waves the curves are convex downwards. For SH-waves transverse isotropy has no effect on curvature.     

INVERSION METHODS FOR τ-p MAPS OF NEAR OFFSET DATA—LINEAR INVERSION*

Conventional velocity analysis, based on the ideas of rms velocity and hyperbolic reflection events in the x-t domain, is restricted in validity to near vertical incidence. Thus analysis of near-offset datasets usually requires the muting of wide-angle reflections from shallow interfaces before the rms velocities are determined. The ray-theoretical integral for the delay time τ, which depends on the slowness p and the velocity function, is valid for all angles. The wide-angle reflections can be used to improve the accuracy of the derived velocity function in the near surface region, if the recorded x-t data are mapped into the τ-p domain. By representing the velocity function between reflectors as a series of gradient zones, i.e. regions with a uniform increase in velocity with depth, the recovery of the velocities may be posed as a matrix linear inverse problem for the slopes of the gradient zones. In order to convert the problem to a linear one, the velocity discontinuities at the reflecting interfaces must be fixed in advance. Their positions are based on the behaviour of the τ-p map of the data. Finding a stable velocity model may require several iterations with the reflecting interfaces at different positions. An understanding of the workings of the inversion algorithm allied with an analysis of the causes of instability aids the search for a stable model.     

Anisotropic velocity analysis1

Results from walkaway VSP and shale laboratory experiments show that shale anisotropy can be significantly anelliptic. Heterogeneity and anellipticity both lead to non-hyperbolic moveout curves and the resulting ambiguity in velocity analysis is investigated for the case of a factorizable anisotropic medium with a linear dependence of velocity on depth. More information can be obtained if there are several reflectors. The method of Dellinger et al. for anisotropic velocity analysis in layered transversely isotropic media is examined and is shown to be restricted to media having relatively small anellipticity. A new scheme, based on an expansion of the inverse-squared group velocity in spherical harmonics, is presented. This scheme can be used for larger anellipticity, and is applicable for horizontal layers having monoclinic symmetry with the symmetry plane parallel to the layers. The method is applied to invert the results of anisotropic ray tracing on a model Sand/shale sequence. For transversely isotropic media with small anisotropy, the scheme reduces to the method of Byun et al. and Byun and Corrigan. The expansion in spherical harmonics allows the P-phase slowness surface of each layer to be determined in analytic form from the layer parameters obtained by inversion without the need to assume that the anisotropy is weak.     

GUIDED LOW-FREQUENCY NOISE FROM AIRGUN SOURCES*

The presence of the water layer in marine seismic prospecting provides an effective waveguide for acoustic energy trapped between the sea-bed and the sea-surface. This energy persists to large ranges and can be the dominant early feature on far-offset traces. On airgun records, there is commonly a lower frequency set of arrivals following the water-trapped waves. These arrivals are not as obvious with higher frequency watergun sources. By using a combination of intercept-time/slowness (τ—p) mapping on observational data and theoretical modelling, we are able to identify the origin of the events. If a very rapid increase in a seismic wavespeed occurs beneath the sea-bed sediments, a new waveguide is formed bounded by the sea surface and this transition zone. The low frequency waves are principally guided within this thicker waveguide. Numerical filtering in the τ—p domain followed by trace reconstruction is very effective in removing the low frequency noise.     

On Christoffel roots for nondetached slowness surfaces

The only restriction on the values of the elasticity parameters is the stability condition. Within this condition, we examine the Christoffel equation for nondetached qP slowness surfaces in transversely isotropic media. If the qP slowness surface is detached, each root of the solubility condition corresponds to a distinct smooth wavefront. If the qP slowness surface is nondetached, the roots are elliptical but do not correspond to distinct wavefronts; also, the qP and qSV slowness surfaces are not smooth.     

Study on the simulation of acoustic logging measurements in horizontal and deviated wells

The conventional acoustic logging interpretation method, which is based on vertical wells that penetrate isotropic formations, is not suitable for horizontal and deviated wells penetrating anisotropic formations. This unsuitability is because during horizontal and deviated well drilling, cuttings will splash on the well wall or fall into the borehole bottom and form a thin bed of cuttings. In addition, the high velocity layers at different depths and intrinsic anisotropy may affect acoustic logging measurements. In this study, we examine how these factors affect the acoustic wave slowness measured in horizontal and deviated wells that are surrounded by an anisotropic medium using numerical simulation. We use the staggered-grid finite difference method in time domain (FDTD) combined with hybrid-PML. First, we acquire the acoustic slowness using a simulated array logging system, and then, we analyze how various factors affect acoustic slowness measurements and the differences between the effects of these factors. The factors considered are high-velocity layers, thin beds of cuttings, dipping angle, formation thickness, and anisotropy. The simulation results show that these factors affect acoustic wave slowness measurements differently. We observe that when the wavelength is much smaller than the distance between the borehole wall and high velocity layer, the true slowness of the formation could be acquired. When the wavelengths are of the same order (i.e., in the near-field scenarios), the geometrical acoustics theory is no longer applicable. Furthermore, when a thin bed of cuttings exists at the bottom of the borehole, Fermat's principle is still applicable, and true slowness can be acquired. In anisotropic formations, the measured slowness changes with increments in the dipping angle. Finally, for a measurement system with specific spacing, the slowness of a thin target layer can be acquired when the distance covered by the logging tool is sufficiently long. Based on systematical simulations with different dipping angles and anisotropy in homogenous TI media, slowness estimation charts are established to quantitatively determine the slowness at any dipping angle and for any value of the anisotropic ratio. Synthetic examples with different acoustic logging tools and different elastic parameters demonstrate that the acoustic slowness estimation method can be conveniently applied to horizontal and deviated wells in TI formations with high accuracy.     

Receiver function method in reflection seismology

The receiver function method was originally developed to analyse earthquake data recorded by multicomponent (3C) sensors and consists in deconvolving the horizontal component by the vertical component. The deconvolution process removes travel path effects from the source to the base of the target as well as the earthquake source signature. In addition, it provides the possibility of separating the emergent P and PS waves based on adaptive subtraction between recorded components if plane waves of constant ray parameters are considered. The resulting receiver function signal is the local PS wave's impulse response generated at impedance contrasts below the 3C receiver.We propose to adapt this technique to the wide‐angle multi‐component reflection acquisition geometry. We focus on the simplest case of land data reflection acquisition. Our adapted version of the receiver function approach consists in a multi‐step procedure that first removes the P wavefield recorded on the horizontal component and next removes the source signature. The separation step is performed in the τ?p domain while the source designature can be achieved in either the τ?p or the t?x domain. Our technique does not require any a priori knowledge of the subsurface. The resulting receiver function is a pure PS‐wave reflectivity response, which can be used for amplitude versus slowness or offset analysis. Stack of the receiver function leads to a high‐quality S wave image.     

Estimation of Seismic Velocities from Deep Reflections in the τ-p Domain

—In deep reflection seismics the estimation of seismic velocities is hampered in most cases due to the low signal level with respect to noise. In the τ-p domain, it is possible to perform the velocity analysis even under such unfavorable signal conditions. This is achieved by making use of special properties of the transform, which enhance the signal-to-noise ratio. Further noise suppression is realized by incorporating filter procedures into the transform algorithm. The velocity analysis itself is also done in the τ-p domain by calculating and evaluating constant velocity gathers. The results can be directly used in the time domain. A mute algorithm, implemented into the τ-p velocity analysis procedure, further reduces noise. This velocity estimation method is discussed with synthetic data and applied to DEKORP data.     

Anomalous Behaviour of Seismic Waves at Interfaces between Anisotropic Media

An account of possible anomalous effects in reflection and refraction of elastic waves at an interface between anisotropic media is presented. These effects are due to anisotropy and they cannot occur at an interface between isotropic media. The shape of the slowness surface (its local deviations from spherical symmetry) is the decisive factor for appearance of these effects. A numerical example of such anomalous behaviour of elastic waves at a free boundary of the crystal of spinel is presented.     

Anomalous behaviour of seismic waves at an interface between anisotropic media

Summary An account of possible anomalous effects in reflection and refraction of elastic waves at an interface between anisotropic media is presented. These effects are due to anisotropy and they cannot occur at an interface between isotropic media. The shape of the slowness surface (its local deviations from spherical symmetry) is the decisive factor for appearance of these effects. A numerical example of such anomalous behaviour of elastic waves at a free boundary of the crystal of spinel is presented.     

Reflection/transmission laws for slowness vectors in viscoelastic anisotropic media

The reflection/transmission laws (R/T laws) of plane waves at a plane interface between two homogeneous anisotropic viscoelastic (dissipative) halfspaces are discussed. Algorithms for determining the slowness vectors of reflected/transmitted plane waves from the known slowness vector of the incident wave are proposed. In viscoelastic media, the slowness vectors of plane waves are complex-valued, p = P + iA, where P is the propagation vector, and A the attenuation vector. The proposed algorithms may be applied to bulk plane waves (A = 0), homogeneous plane waves (A ≠ 0, P and A parallel), and inhomogeneous plane waves (A ≠ 0, P and A non-parallel). The manner, in which the slowness vector is specified, plays an important role in the algorithms. For unrestricted anisotropy and viscoelasticity, the algorithms require an algebraic equation of the sixth degree to be solved in each halfspace. The degree of the algebraic equation decreases to four or two for simpler cases (isotropic media, plane waves in symmetry planes of anisotropic media). The physical consequences of the proposed algorithms are discussed in detail. vcerveny@seis.karlov.mff.cuni.cz     

Inhomogeneous Harmonic Plane Waves in Viscoelastic Anisotropic Media

In viscoelastic media, the slowness vector p of plane waves is complex-valued, p = P + iA. The real-valued vectors P and A are usually called the propagation and the attenuation vector, repectively. For P and A nonparallel, the plane wave is called inhomogeneousThree basic approaches to the determination of the slowness vector of an inhomogeneous plane wave propagating in a homogeneous viscoelastic anisotropic medium are discussed. They differ in the specification of the mathematical form of the slowness vector p. We speak of directional specification, componental specification and mixed specification of the slowness vector. Individual specifications lead to the eigenvalue problems for 3 × 3 or 6 × 6 complex-valued matrices.In the directional specification of the slowness vector, the real-valued unit vectors N and M in the direction of P and A are assumed to be known. This has been the most common specification of the slowness vector used in the seismological literature. In the componental specification, the real-valued unit vectors N and M are not known in advance. Instead, the complex-valued vactorial component p of slowness vector p into an arbitrary plane with unit normal n is assumed to be known. Finally, the mixed specification is a special case of the componental specification with p purely imaginary. In the mixed specification, plane represents the plane of constant phase, so that N = ±n. Consequently, unit vector N is known, similarly as in the directional specification. Instead of unit vector M, however, the vectorial component d of the attenuation vector in the plane of constant phase is known.The simplest, most straightforward and transparent algorithms to determine the phase velocities and slowness vectors of inhomogeneous plane waves propagating in viscoelastic anisotropic media are obtained, if the mixed specification of the slowness vector is used. These algorithms are based on the solution of a conventional eigenvalue problem for 6 × 6 complex-valued matrices. The derived equations are quite general and universal. They can be used both for homogeneous and inhomogeneous plane waves, propagating in elastic or viscoelastic, isotropic or anisotropic media. Contrary to the mixed specififcation, the directional specification can hardly be used to determine the slowness vector of inhomogeneous plane waves propagating in viscoelastic anisotropic media. Although the procedure is based on 3 × 3 complex-valued matrices, it yields a cumbersome system of two coupled equations.     

QUASI-ANALYTIC CONVOLUTION SOLUTION OF THE ELECTROMAGNETIC FIELD*

The objective of this study is to generate the separation-distance-domain (r-domain) transformation of the theoretically calculated wave number domain (m-domain) electromagnetic induction field component B_z(m, ω) of a stratified medium and to search for interpretive information which has been absent in the previously achieved numerical solutions of the problem. The r-domain kernel R?(r, ω) function defining the induction field appears to adequately reflect the layering and electrical properties of the medium if it is expressed as a function of the frequency if the source-receiver separation r is small with respect to the thickness of the first layer. However, exact values of the conductivity cannot be distinguished from those of the neighboring values unless a resistive basement layer is present. This feature is the result of the truncation in series representation of the kernel function R?(m, ω). However, this truncation is regarded as significant in the case of a conductive first layer. In m-domain static-zone studies, a conductive first layer slightly influences its r-domain correspondent. Although the computational cost of obtaining the kernel B(r, ω) by evaluation of the convolution in a cylindrical coordinate system is high, this semi-analytic solution is still superior to those based on the asymptotic assumptions.     

Asymptotic Elastodynamic Green Function in the Kiss Singularity in Homogeneous Anisotropic Solids

The far-field asymptotic formula is derived for the elastodynamic Green function in the kiss singularity in homogeneous anisotropic solids. In contrast to standard asymptotics in regular directions the derived formula is more complex and expressed in the form of a 1-D integral. This integral is specified for the kiss singularity along the symmetry axis in transverse isotropy and along the fourfold symmetry axes in tetragonal and cubic symmetries. The shape of the slowness surface in the singularity is regular in transverse isotropy and the amplitude of the Green function is expressed by means of the Gaussian curvature of this surface in the singularity. However, the shape of the slowness surface is irregular and the Gaussian curvature is not defined in the singularity in tetragonal or cubic symmetries. In this case, the amplitude of the Green function is expressed by means of the generalized Gaussian curvature.     

τ-p SEISMIC DATA FOR VISCOELASTIC MEDIA – PART 1: MODELLING1

Viscoelastic modelling reveals that the interaction of compressional-wave velocity C_p, compressional-wave quality factor Q_p, shear-wave velocity C_s, shear-wave quality factor Q_s and Poisson's ratio as a function of time intercept τ and ray parameter p, is complicated; however, distinct, potentially diagnostic behaviours are seen for different combinations of viscoelastic parameters. Synthetic seismograms for three viscoelastic reservoir models show that variations in the Poisson's ratio produce visible differences when compared to the corresponding elastic synthetic seismograms; these differences are attributable to interaction of the elastic parameters with Q_p and Q_s. When the P-wave acoustic impedance contrast is small, viscoelastic effects become more apparent and more useful for interpretation purposes. The corresponding amplitude and net phase spectra reveal significant differences between the elastic and the viscoelastic responses. When P-wave reflectivities are large, they tend to dominate the total response and to mask the Q reflectivity effects. The attenuation effects are manifested as an amplitude decay that increases with both time and ray parameter. The sensitivity of the computed seismic responses for various combinations of viscoelastic parameters suggests the opportunity for diagnostic interpretation of τ-p seismic data. The interpretation of the viscoelastic parameters can permit a better understanding of the rock types and pore fluid distribution existing in the subsurface.     

A LAYER-STRIPPING TECHNIQUE FOR THE SUPPRESSION OF WATER-BOTTOM MULTIPLE REFLECTIONS*

A new method to suppress water-bottom multiples (water-bottom reverberations) uses the fact that in the domain of intercept time and ray parameter (τ–p domain) the water-bottom reverberations are strictly periodical for a horizontal flat sea bottom. Using this property a comb filter can be designed. The window of the filter should be approximately equal to the duration of a source pulse. The algorithm finds the maximum of the periodical energy throughout the τ–p domain and then designs the comb filter which eliminates the water bottom reverberations from each trace in the τ– p domain. This process can be repeated for higher order reverberations. Finally the τ–p domain with attenuated multiples is transformed back to the conventional x -- t space. The method is illustrated on a variety of synthetic data and on a set of real marine CMP data acquired in the North Sea near the Norwegian shore.     

Slowness-domain kinematical characteristics for horizontally layered orthorhombic media. Part I: Critical slowness match

Kinematical characteristics of reflected waves in anisotropic elastic media play an important role in the seismic imaging workflow. Considering compressional and converted waves, we derive new, azimuthally dependent, slowness-domain approximations for the kinematical characteristics of reflected waves (radial and transverse offsets, intercept time and traveltime) for layered orthorhombic media with varying azimuth of the vertical symmetry planes. The proposed method can be considered an extension of the well-known ‘generalized moveout approximation’ in the slowness domain, from azimuthally isotropic to azimuthally anisotropic models. For each slowness azimuth, the approximations hold for a wide angle range, combining power series coefficients in the vicinity of both the normal-incidence ray and an additional wide-angle ray. We consider two cases for the wide-angle ray: a ‘critical slowness match’ and a ‘pre-critical slowness match’ studied in Parts I and II of this work, respectively. For the critical slowness match, the approximations are valid within the entire slowness range, up to the critical slowness. For the ‘pre-critical slowness match’, the approximations are valid only within the bounded slowness range; however, the accuracy within the defined range is higher. The critical slowness match is particularly effective when the subsurface model includes a dominant high-velocity layer where, for nearly critical slowness values, the propagation in this layer is almost horizontal. Comparing the approximated kinematical characteristics with those computed by numerical ray tracing, we demonstrate high accuracy.     

Approximate Relation Between the Ray Vector and the Wave Normal in Weakly Anisotropic Media

Determination of the ray vector (the unit vector specifying the direction of the group velocity vector) corresponding to a given wave normal (the unit vector parallel to the phase velocity vector or slowness vector) in an arbitrary anisotropic medium can be performed using the exact formula following from the ray tracing equations. The determination of the wave normal from the ray vector is, generally, a more complicated task, which is usually solved iteratively. We present a first-order perturbation formula for the approximate determination of the ray vector from a given wave normal and vice versa. The formula is applicable to qP as well as qS waves in directions, in which the waves can be dealt with separately (i.e. outside singular directions of qS waves). Performance of the approximate formulae is illustrated on models of transversely isotropic and orthorhombic symmetry. We show that the formula for the determination of the ray vector from the wave normal yields rather accurate results even for strong anisotropy. The formula for the determination of the wave normal from the ray vector works reasonably well in directions, in which the considered waves have convex slowness surfaces. Otherwise, it can yield, especially for stronger anisotropy, rather distorted results.