Original Wavefield Research Articles

Elastic-wave reverse-time migration based on decoupled elastic-wave equations and inner-product imaging condition

Polarity reversal in converted wave images of elastic reverse time migration destructs the reflection events after stacking multi-shot migration profile. We derive a new imaging method for elastic reverse-time migration to automatically circumvent polarity reversal. Instead of obtaining scalar P-wave and vector S-wave potentials from the wavefield by using Helmholtz decomposition as in conventional methods, we obtain vector P- and S-wave displacement wavefields based on a decoupled elastic wave equation, which denotes the displacement component along the propagation direction and perpendicular to the propagation direction, respectively. Based on this decomposition method, the vector P- and S-wavefields preserve the amplitude and phase attributes of the original wavefield. As for the vector wavefields (vector P- and S-wave displacement wavefields), the inner-product imaging condition is proposed to extract reflectivity of specified wave modes at interfaces. The analysis of the imaging kernel demonstrates this imaging condition is valid not only for pure-mode imaging (PP and SS), but also for converted wave imaging (PS and SP) of ground-based seismic exploration. With this new method, we do not have to correct the polarity reversal in converted wave images, which is an essential step in the conventional method with expensive computation costs. Numerical examples with synthetic data have shown that the inner product imaging method works and the quality of the imaging events is effectively improved.

Efficient angle-domain common-image gathers using Cauchy-condition-based polarization vectors

Because angle-domain common-image gathers (ADCIGs) from reverse time migration (RTM) are capable of obtaining the correct illumination of a subsurface geologic structure, they provide more reliable information for velocity model building, amplitude-variation versus angle analysis, and attribute interpretation. The approaches for generating ADCIGs mainly consist of two types: (1) indirect approaches that convert extended image gathers into ADCIGs and (2) direct approaches that first obtain propagating angles of wavefronts and then map the imaging result to the angle domain. In practice, however, generation of ADCIGs usually incurs high computational cost, poor resolution, and other drawbacks. To generate efficient ADCIGs using RTM methods, we have introduced a novel approach to obtain polarization vectors — directions of particle motion — from the Cauchy wavefield (CWF) and an efficient localized plane-wave decomposition algorithm to implement the angle-domain imaging condition. The CWF is a wavefield constructed from the Cauchy condition of the wave equation at any given time, and it only contains negative frequencies of the original wavefield so that the polarization vector is obtained from the local CWF in the wavenumber domain. With polarization vectors at our disposal, we have further developed an efficient localized plane-wave decomposition algorithm to implement the angle-domain imaging condition. Numerical examples have indicated that the new approach is able to handle complex wave phenomenon and has advantages in illuminating subsurface structure.

Prestack exploding reflector modelling and migration for anisotropic media

ABSTRACTThe double‐square‐root equation is commonly used to image data by downward continuation using one‐way depth extrapolation methods. A two‐way time extrapolation of the double‐square‐root‐derived phase operator allows for up and downgoing wavefields but suffers from an essential singularity for horizontally travelling waves. This singularity is also associated with an anisotropic version of the double‐square‐root extrapolator. Perturbation theory allows us to separate the isotropic contribution, as well as the singularity, from the anisotropic contribution to the operator. As a result, the anisotropic residual operator is free from such singularities and can be applied as a stand alone operator to correct for anisotropy. We can apply the residual anisotropy operator even if the original prestack wavefield was obtained using, for example, reverse‐time migration. The residual correction is also useful for anisotropic parameter estimation. Applications to synthetic data demonstrate the accuracy of the new prestack modelling and migration approach. It also proves useful in approximately imaging the Vertical Transverse Isotropic Marmousi model.

Non‐convex compressed sensing with frequency mask for seismic data reconstruction and denoising

ABSTRACTCompressed Sensing has recently proved itself as a successful tool to help address the challenges of acquisition and processing seismic data sets. Compressed sensing shows that the information contained in sparse signals can be recovered accurately from a small number of linear measurements using a sparsity‐promoting regularization. This paper investigates two aspects of compressed sensing in seismic exploration: (i) using a general non‐convex regularizer instead of the conventional one‐norm minimization for sparsity promotion and (ii) using a frequency mask to additionally subsample the acquired traces in the frequency‐space () domain. The proposed non‐convex regularizer has better sparse recovery performance compared with one‐norm minimization and the additional frequency mask allows us to incorporate a priori information about the events contained in the wavefields into the reconstruction. For example, (i) seismic data are band‐limited; therefore one can use only a partial set of frequency coefficients in the range of reflections band, where the signal‐to‐noise ratio is high and spatial aliasing is low, to reconstruct the original wavefield, and (ii) low‐frequency characteristics of the coherent ground rolls allow direct elimination of them during reconstruction by disregarding the corresponding frequency coefficients (usually bellow 10 Hz) via a frequency mask. The results of this paper show that some challenges of reconstruction and denoising in seismic exploration can be addressed under a unified formulation. It is illustrated numerically that the compressed sensing performance for seismic data interpolation is improved significantly when an additional coherent subsampling is performed in the domain compared with the domain case. Numerical experiments from both simulated and real field data are included to illustrate the effectiveness of the presented method.

On the numerical implementation of time-reversal mirrors for tomographic imaging

A general approach for constructing numerical equivalents of time-reversal mirrors is introduced. These numerical mirrors can be used to regenerate an original wavefield locally within a confined volume of arbitrary shape. Though time-reversal mirrors were originally designed to reproduce a time-reversed version of an original wavefield, the proposed method is independent of the time direction and can be used to regenerate a wavefield going either forward in time or backward in time. Applications to computational seismology and tomographic imaging of such local wavefield reconstructions are discussed. The key idea of the method is to directly express the source terms constituting the time-reversal mirror by introducing a spatial window function into the wave equation. The method is usable with any numerical method based on the discrete form of the wave equation, for example, with finite difference (FD) methods and with finite/spectral elements methods. The obtained mirrors are perfect in the sense that no additional error is introduced into the reconstructed wavefields apart from rounding errors that are inherent in floating-point computations. They are fully transparent as they do not interact with waves that are not part of the original wavefield and are permeable to these. We establish a link between some hybrid methods introduced in seismology, such as wave-injection, and the proposed time-reversal mirrors. Numerical examples based on FD and spectral elements methods in the acoustic, the elastic and the visco-elastic cases are presented. They demonstrate the accuracy of the method and illustrate some possible applications. An alternative implementation of the time-reversal mirrors based on the discretization of the surface integrals in the representation theorem is also introduced. Though it is out of the scope of the paper, the proposed method also apply to numerical schemes for modelling of other types of waves such as electromagnetic waves.

Error estimates for Gaussian beam superpositions

Gaussian beams are asymptotically valid high frequency solutions to hyperbolic partial differential equations, concentrated on a single curve through the physical domain. They can also be extended to some dispersive wave equations, such as the Schrodinger equation. Superpositions of Gaussian beams provide a powerful tool to generate more general high frequency solutions that are not necessarily concentrated on a single curve. This work is concerned with the accuracy of Gaussian beam superpositions in terms of the wavelength epsilon. We present a systematic construction of Gaussian beam superpositions for all strictly hyperbolic and Schrodinger equations subject to highly oscillatory initial data of the form Ae(i Phi/) (epsilon). Through a careful estimate of an oscillatory integral operator, we prove that the k-th order Gaussian beam superposition converges to the original wave field at a rate proportional to epsilon(k/2) in the appropriate norm dictated by the well-posedness estimate. In particular, we prove that the Gaussian beam superposition converges at this rate for the acoustic wave equation in the standard, epsilon-scaled, energy norm and for the Schrodinger equation in the L-2 norm. The obtained results are valid for any number of spatial dimensions and are unaffected by the presence of caustics. We present a numerical study of convergence for the constant coefficient acoustic wave equation in R-2 to analyze the sharpness of the theoretical results.

Downward and upward continuation of 2-D seismic data to eliminate ocean bottom topography’s effect

In order to eliminate the effect of ocean bottom topography on seismic wave field, we transformed curved (x, z) coordinate system grids into rectangular (ξ, η) coordinate system grids and derived a 2-D scalar acoustic wave equation in the ξ, η domain. The seismic wave field collected at the sea surface was downward continued to the ocean bottom by the inverse finite difference method with the water velocity and then was reversely continued to the ocean surface by the finite difference method using the layer velocity from just below the ocean bottom in the (ξ, η) domain. Simulation calculations and practical application show that this method can not only remove the reflection travel time distortion but also correct the dynamic parameter changes caused by the ocean bottom topography. The inverted velocity after wave field continuation is much more accurate than before continuation and the image section was greatly improved compared to the original wave field.

2D and 3D elastic wavefield vector decomposition in the wavenumber domain for VTI media

A pragmatic decomposition of a vector wavefield into P- and S-waves is based on the Helmholtz theory and the Christoffel equation. It is applicable to VTI media when the plane-wave polarization is continuous in the vicinity of a given wavenumber and is uniquely defined by that wavenumber, except for the kiss singularities on the VTI symmetry axis. Unlike divergence and curl, which separate the wavefield into a scalar and a vector field, the decomposed P- and S-wavefields are both vector fields, with correct amplitude, phase, and physical units. If the vector components of decomposed wavefields are added, they reconstruct those of the original input wavefield. Wavefield propagation in any portions of a VTI medium that have the same polarization distribution (i.e., the same eigenvector) in the wavenumber domain have the same decomposition operators and can be recon-structed with a single 3D Fourier transform for each operator (e.g., one for P-waves and one for S-waves).This applies to isotropic wavefields and to VTI anisotropic wavefields, if the polarization distribution is constant, regardless of changes in the velocity. Because the anisotropic phase polarization is local, not global, the wavefield decomposition for inhomogeneous anisotropic media needs to be done separately for each region that has a different polarization distribution. The complete decomposed vector wavefield is constructed by combining the P-, SV-, and SH-wavefields in each region into the corresponding composite P-, SV-, and SH-wavefields that span the model. Potential practical applications include extraction of separate images for different wave types in prestack reverse time migration, inversion, or migration velocity analysis, and calculation of wave-propagation directions for common-angle gathers.

Optimized two-step phase-shifting algorithm applied to image encryption

The two-step phase-shifting algorithm proposed previously is optimized, the original object wave field can be reconstructed by only one phase shift value in (0, π) and two interferograms, with the removal (or suppression) of background intensity (or dc term), and the additional measurements such as the object wave intensity, reference wave intensity, etc., are no longer required. Together with double random phase encoding technique in the Fresnel domain, the optimized two-step phase-shifting algorithm is then applied to image encryption system. The feasibility of the proposed scheme is verified by computer simulation. Furthermore, the sensitivity of geometrical keys has also been tested and analyzed.

Recovery of High Frequency Wave Fields for the Acoustic Wave Equation

Computation of high frequency solutions to wave equations is important in many applications and notoriously difficult in resolving wave oscillations. Gaussian beams are asymptotically valid high frequency solutions concentrated on a single curve through the physical domain, and superposition of Gaussian beams provides a powerful tool to generate more general high frequency solutions to PDEs. An alternative way to compute Gaussian beam components such as phase, amplitude, and Hessian of the phase is to capture them in phase space by solving Liouville-type equations on uniform grids. Following [H. Liu and J. Ralston, Multiscale Model. Simul., to appear] we present a systematic construction of asymptotic high frequency wave fields from computations in phase space for acoustic wave equations; the superposition of phase space based Gaussian beams over two moving domains is shown to be necessary. Moreover, we prove that the kth order Gaussian beam superposition converges to the original wave field in the energy norm at the rate of $\epsilon^{\frac{k}{2}+\frac{1-n}{4}}$ in dimension n.

Recovery of High Frequency Wave Fields from Phase Space–Based Measurements

Computation of high frequency solutions to wave equations is important in many applications, and notoriously difficult in resolving wave oscillations. Gaussian beams are asymptotically valid high frequency solutions concentrated on a single curve through the physical domain, and superposition of Gaussian beams provides a powerful tool for generating more general high frequency solutions to PDEs. An alternative way to compute Gaussian beam components such as phase, amplitude, and Hessian of the phase is to capture them in phase space by solving Liouville-type equations on uniform grids. In this work we review and extend recent constructions of asymptotic high frequency wave fields from computations in phase space. We give a new level set method of computing the Hessian and higher derivatives of the phase. Moreover, we prove that the kth order phase space–based Gaussian beam superposition converges to the original wave field in $L^2$ at the rate of $\epsilon^{\frac{k}{2}-\frac{n}{4}}$ in dimension n.

The adjoint method in seismology—: II. Applications: traveltimes and sensitivity functionals

Sensitivity functionals which allow us to express the total derivative of a physical observable with respect to the model parameters, are defined on the basis of the adjoint method. The definition relies on the existence of Green’s functions for both the original and the adjoint problem. Using the acoustic wave equation in a homogeneous and unbounded medium, it is shown that the first derivative of the wavefield u with respect to the model parameter p = c − 2 ( c is the wave speed) does not contain traveltime information. This property also influences objective functions defined on u and in particular the least squares objective function. Therefore, a waveform inversion should either be complemented by a traveltime tomography or work with initially very long wavelengths that decrease in the course of the iteration. The definition of the sensitivity functionals naturally introduces waveform sensitivity kernels. Analytic examples are shown for the case of an isotropic, elastic and unbounded medium. In the case of a double couple source there are three classes of sensitivity kernels, one for each of the parameters λ , μ (Lamé parameters) and ρ (density). They decompose into P → P, P → S, S → P and S → S kernels and can be described by third-order tensors incorporating the radiation patterns of the original wavefield and the adjoint wavefield. An analysis of the sensitivity kernels for density suggests that a waveform inversion procedure should exclude either the direct waves or the source and receiver regions. S-wave sensitivity kernels, corresponding to conversions from or to S-waves, are larger than the kernels corresponding to P-waves only. This implies that S-wave residuals caused by parameter differences between the true Earth model and the numerical model will dominate a waveform inversion that does not account for that effect.

Statistics of breaking waves and its applications to estimation of air-sea fluxes (I)

Wave breaking statistics, such as the whitecap coverage and average volume of broken seawater, are evaluated in terms of wave parameters by use of wave breaking model (Yuan et al., 1988) taking the fifth order Stokes's wave as the analog of the original wave field. Based on the observed fact that breaking waves play an important role in the exchange of mass, momentum and energy between the atmosphere and the ocean, the influence of wave breaking on air-sea fluxes of heat and moisture is investigated. Theoretical expressions of bubble-volume flux and sea spray spectrum at the sea surface and models for bubble-induced and spray droplet-induced heat and moisture fluxes are established. This work can be taken as the basis for further understanding the mechanism of air-sea coupling and parameterization models.

Statistics of breaking waves and its applications to estimation of air-sea fluxes (I) ??Theoretical models

Wave breaking statistics, such as the whitecap coverage and average volume of broken seawater, are evaluated in terms of wave parameters by use of wave breaking model (Yuan et al., 1988) taking the fifth order Stokes’s wave as the analog of the original wave field. Based on the observed fact that breaking waves play an important role in the exchange of mass, momentum and energy between the atmosphere and the ocean, the influence of wave breaking on air-sea fluxes of heat and moisture is investigated. Theoretical expressions of bubble-volume flux and sea spray spectrum at the sea surface and models for bubble-induced and spray droplet-induced heat and moisture fluxes are established. This work can be taken as the basis for further understanding the mechanism of air-sea coupling and parameterization models.

Patterning of particulate films using Faraday waves

Faraday waves were used to create a patterned particulate film on the surface of a glass substrate. The process involves deposition of a thin liquid film on the substrate, where the liquid is comprised of water and a concentrated suspension of particles. Faraday waves were created on the surface of the liquid film by subjecting the substrate and film to a vertical oscillation. Because the films are thin, the standing waves create velocity fluctuations that are not negligible near the substrate surface. As the particles in the liquid settled to the bottom, the flow field caused by the waves caused preferential particle deposition organized in a pattern determined by the Faraday wave field. Subsequent evaporation of the liquid film left a particulate film on the solid substrate that retained the pattern of the original wave field. We refer to these particulate films as Faraday films. The length scale of the patterns in these films was on the order of millimeters in the experiments described here. Use of higher frequency oscillations could permit the formation of Faraday films having smaller spatial scales.

Application of time reverse acoustics to acoustic measurements in a noisy environment

For acoustic measurements made within noisy environments, time reversed acoustic methods can be used to improve the measurement. After reversing the array received time signal, rebroadcast of the reversed signals recreates the original source wave field at its original source location. In some measurement cases a background noise is created in the process of generating the source signal of interest; for example, wind and water tunnel acoustic experiments. By using time reverse acoustics, the signal-to-noise ratio can be improved. First, the source signal together with the background noise is recorded using an array of many reciprocal transducers. After the tunnel flow is stopped, the recorded signals are played in a time-reversed fashion through the reciprocal transducer elements of the array. If a receiver is located at the correct source location, only the source part of the rebroadcast signal is superimposed positively, with negligible influence from the background sources that originate elsewhere. Here, we describe noisy environment, free-field, time reversed acoustic simulations. The computer simulations show an improved signal-to-noise ratio. The advantages and disadvantages of using continuous or pulsed recorded signal are discussed. [Work supported by the Applied Research Laboratory E/F Program, and ONR, Code 333.]

Some remarks on surface multiple attenuation

The surface multiple attenuation algorithm discussed in this paper is a prestack inversion of a surface‐recorded, 2-D wavefield that aims to remove all orders of all surface multiples present within the wavefield. Although the algorithm requires no assumptions or modeling regarding the positions and reflection coefficients of the multiple‐causing reflectors, it does require complete internal physical consistency between primary and multiple events—something that exists only in ideal 2-D data sets. In field data sets the physical consistency between primaries and multiples is disturbed by phenomena such as variations in the acquisition wavelet, cable feathering, cross‐line dip, a finite near offset, and unequal or too coarse spatial sampling in source and receiver coordinates. Careful survey design can minimize the impact of those phenomena on surface multiple attenuation. If it is not too large, trace extrapolation can solve the finite near‐offset problem. Minor adjustments to the algorithm allow processing of data for which the source and receiver intervals differ by an integer multiple, although for those and other acquisition geometries, trace interpolation may be preferred. In the f-x domain, surface multiple attenuation can be formulated as an equation whose straightforward solution involves the inversion of a large matrix that is a function of the acquisition wavelet. Since that wavelet is generally unknown, solving this matrix equation becomes an optimization problem. Many matrix inversions are needed to estimate the acquisition wavelet that leads to the best multiple suppression, rendering the straightforward solution to the surface multiple attenuation equation quite costly. We offer two alternative approaches. In our first approach we compute an eigenvalue decomposition of the large matrix, allowing the equation to be recast so that the wavelet dependency appears in a diagonal matrix for which repetitive inversion is trivial. In our second approach we begin by using the surface multiple attenuation algorithm with a fixed, approximately correct wavelet to compute the surface multiple wavefield. We then filter the predicted multiples adaptively to match the actual multiples in the original wavefield and subtract these filtered multiples from the original wavefield. The second approach is relatively inexpensive and to some extent can cope with physical inconsistencies between primaries and multiples caused by field data set imperfections.

A method for direct computation of the differential seismogram with respect to the velocity change in a layered elastic solid

Abstract We derived a method to compute the differential seismogram directly in a layered elastic half-space. This method first computes a differential field caused by the velocity change in the layer and then multiplies it with the original elastic wave field. This product acts like a source, and the seismogram caused by this source term is the differential seismogram. The differential waves are propagated directly to the receiver in frequency/wavenumber domain using the generalized reflection and transmission coefficient method. In comparison with the typical finite-difference approach to compute a complete set of differential seismograms at one station in an N-layered earth model, our method not only saves about N/2 folds of computation time but also leads to a more accurate solution. Thus, it provides an efficient and direct procedure to compute the differential seismogram for a complete waveform fitting in a layered elastic earth model.

Phase determination by conjugate wave-front generation

A new method is proposed to solve the phase problem of the optical coherence theory. The method consists of generating a phase-conjugate wave using a suitable nonlinear process and superimposing it with the original wavefield. We discuss how the method could be applied to obtain the intensity distribution across an incoherent light source.