Adjoint Wavefield Research Articles

Elastic time-reverse imaging of backscattered surface-wave data

Imaging shallow near-vertical or short-wavelength structure using seismic data is an important goal in near-surface geophysical applications. Numerous seismic analysis approaches use backscattered analysis of surface waves for these purposes; however, there is a need for new associated frameworks that can take advantage of these events to generate interpretable images in the depth-migrated domain. By drawing a connection with elastic time-reverse imaging (E-TRI) approaches developed for microseismic event analysis, an efficient elastic migration approach is developed that uses two adjoint wavefields defined by partitioned outward-propagating transmission and inward-propagating scattered surface-wave data. Subsurface images are generated by applying the energy-norm imaging condition that correlates temporally and spatially colocated energy throughout the two time-reversed adjoint wavefield domains. Synthetic tests highlight the ability of the E-TRI framework to image shallow (quasi)vertical subsurface discontinuities. Results from a field test using near-surface seismic data acquired in the Majes region of Arequipa, Peru, demonstrate the ability of the E-TRI framework to highlight short-wavelength structure identified in a parallel elastic full-waveform inversion analysis. The results suggest that the E-TRI technique should be applicable in a range of near-surface geotechnical and environmental applications.

Compensating attenuation effects in full-waveform inversion with dissipation-dispersion decoupling

Seismic attenuation poses challenges to velocity model building from acoustic/elastic full-waveform inversion (FWI). In particular, when constructing the FWI gradient for velocity inversion, recent studies have indicated that the high-attenuation structure could distort the gradient by damping the amplitudes and shifting the kinematic phases, resulting in an imbalanced update and thus unreliable velocity model. These Q effects are particularly significant in the reflection acquisition geometry due to the “double-damping” issue. Here, we develop a Q-compensated FWI algorithm for constructing a Q-free FWI gradient. By using a recently developed viscoacoustic wave propagator, this compensation can be done conveniently by keeping the dispersion term and flipping the dissipation term in the wave equation when we simulate the forward and time-reversed adjoint wavefields. The resultant gradient obtained by interacting these two wavefields has correct kinematics and Q-free amplitudes. This Q-compensated FWI can balance the update between Q- and no- Q-areas, which we determine using synthetic examples. In addition, we illustrate how to take advantage of the dissipation-dispersion decoupling to determine the anomalous Q value in the Q model building workflow via Q-compensated reverse time migration.

Modeling and inversion in acoustic-elastic coupled media using energy-stable summation-by-parts operators

Seismic data acquired at the seafloor are valuable in characterizing the subsurface and monitoring producing hydrocarbon fields. To fully use such data, a stable, accurate, and efficient numerical scheme is needed that accounts for acoustic and elastic wave propagation, their interaction at the seafloor interface, and for sources and receivers placed on either side of that interface. Existing methods either make incorrect assumptions or have high implementation and computational costs. We have developed a high-order finite-difference summation-by-parts framework for the acoustic-elastic wave equations in second-order form (in terms of displacements in the solid and an extension of velocity potential in the fluid). Our modified discretization of the elastic operator overcomes the dispersion errors known to plague displacement-based schemes in the high VP/ VS limit. We weakly impose boundary and interface conditions using simultaneous approximation terms, leading to an energy-stable numerical scheme that rigorously handles point injection and extraction. The fully discrete system is self-adjoint after a time reversal and a sign flip and is furthermore a high-order accurate discretization of the continuous problem. The self-adjointness ensures that forward and adjoint wavefields are computed with similar accuracy and simplifies gradient computation for inversion purposes. We find that our numerical scheme achieves accuracy comparable to a more computationally expensive spectral-element method and demonstrate its application to full-waveform inversion using the Marmousi2 model.

A new scheme of wavefield decomposed elastic least-squares reverse time migration

Elastic least-squares reverse time migration (ELSRTM) describes the reflectivity of the underground media more accurately than acoustic LSRTM in theory while suffering from the P- and S-waves crosstalk artifacts. We propose a new wavefield decomposed ELSRTM scheme to alleviate these crosstalk artifacts, which is different from conventional methods. In our new scheme, we implement the wavenumber domain elastic wavefield vector decomposition equivalently in the time-space domain to decompose source wavefield without Fourier transform, but with high precision. Then we decompose adjoint wavefield by constructing the shear component in a decoupled adjoint wave equation. Finally, based on elastic impedance parameterization, we derive the gradients with respect to elastic reflectivity in the wavefield-decomposed ELSRTM. Numerical examples show that our method is feasible even when applied to models with complex and uncorrelated P- and S-wave velocity structures.

Seismic Imaging of Complex Velocity Structures by 2D Pseudo-Viscoelastic Time-Domain Full-Waveform Inversion

In the presented study, multi-parameter inversion in the presence of attenuation is used for the reconstruction of the P- and the S- wave velocities and the density models of a synthetic shallow subsurface structure that contains a dipping high-velocity layer near the surface with varying thicknesses. The problem of high-velocity layers also complicates selection of an appropriate initial velocity model. The forward problem is solved with the finite difference, and the inverse problem is solved with the preconditioned conjugate gradient. We used also the adjoint wavefield approach for computing the gradient of the misfit function without explicitly build the sensitivity matrix. The proposed method is capable of either minimizing the least-squares norm of the data misfit or use the Born approximation for estimating partial derivative wavefields. It depends on which characteristics of the recorded data—such as amplitude, phase, logarithm of the complex-valued data, envelope in the misfit, or the linearization procedure of the inverse problem—are used. It showed that by a pseudo-viscoelastic time-domain full-waveform inversion, structures below the high-velocity layer can be imaged. However, by inverting attenuation of P- and S- waves simultaneously with the velocities and mass density, better results would be obtained.

Parsimonious truncated Newton method for time-domain full-waveform inversion based on the Fourier-domain full-scattered-field approximation

To exploit Hessian information in full-waveform inversion (FWI), the matrix-free truncated Newton method can be used. In such a method, Hessian-vector product computation is one of the major concerns due to the huge memory requirements and demanding computational cost. Using the adjoint-state method, the Hessian-vector product can be estimated by zero-lag crosscorrelation of the first-/second-order incident wavefields and the second-/first-order adjoint wavefields. Different from the implementation in frequency-domain FWI, Hessian-vector product construction in the time domain becomes much more challenging because it is not affordable to store all of the time-dependent wavefields. The widely used wavefield recomputation strategy leads to computationally intensive tasks. We have developed an efficient alternative approach to computing the Hessian-vector product for time-domain FWI. In our method, discrete Fourier transform is applied to extract frequency-domain components of involved wavefields, which are used to compute wavefield crosscorrelation in the frequency domain. This makes it possible to avoid reconstructing the first- and second-order incident wavefields. In addition, a full-scattered-field approximation is proposed to efficiently simplify the second-order incident and adjoint wavefield computation, which enables us to refrain from repeatedly solving the first-order incident and adjoint equations for the second-order incident and adjoint wavefields (re)computation. With our method, the computational time can be reduced by 70% and 80% in viscous media for Gauss-Newton and full-Newton Hessian-vector product construction, respectively. The effectiveness of our method is also verified in the frame of a 2D multiparameter inversion, in which our method almost reaches the same iterative convergence of the conventional time-domain implementation.

Attenuation compensation for wavefield‐separation‐based least‐squares reverse time migration in viscoelastic media

ABSTRACTElastic least‐squares reverse time migration can image the multicomponent seismic recordings. However, on the one hand, without considering the intrinsic attenuation of the subsurface, it may produce blurred reflectivity images with incorrect positions of reflectors for seismic recordings with strong attenuation. On the other hand, the crosstalk artefacts created by the different wave modes may severely degrade the imaging quality. To alleviate these crosstalk artefacts and compensate for the attenuation effects, we have developed a time‐domain wavefield‐separation‐based least‐squares reverse time migration approach in viscoelastic media, which utilizes the viscoelastic wave equation on the basis of the standard linear solid model to simulate intrinsic subsurface attenuation. The key point of the proposed approach is that we utilize the separated gradient contribution of PP‐ and PS‐wave modes based on the wavefield separation technique in viscoelastic media to construct the P‐wave and S‐wave velocity images, respectively. Unlike the conventional wavefield‐separation‐based least‐squares reverse time migration approach, which generally uses the new stress–velocity equations to formulate, our proposed wavefield separation scheme fully depends on the conventional viscoelastic wave equation and its adjoint wave equation. By introducing the pure P‐wave stress in the forward and adjoint wavefields, the coupled wavefields can be well decomposed, which requires much less computational costs than the wavenumber‐domain wavefield separation scheme. Numerical examples using the layered and Marmousi‐2 models have shown that the proposed approach can improve the migration image quality under geological conditions with strong attenuation where elastic least‐squares reverse time migration may produce blurred and unfocused events. Meanwhile, the proposed least‐squares reverse time migration approach with the wavefield separation scheme has a better convergence rate and produces fewer crosstalk artefacts than that without the wavefield separation scheme.

Time-dispersion correction for arbitrary even-order Lax-Wendroff methods and the application on full-waveform inversion

A second-order accurate finite-difference (FD) approximation is commonly used to approximate the second-order time derivative of a wave equation. The second-order accurate FD scheme may introduce time dispersion in wavefield extrapolation. Lax-Wendroff methods can suppress such dispersion by replacing the high-order time FD error terms with space FD error-correcting terms. However, the time dispersion cannot be completely eliminated and the computation cost dramatically increases with increasing order of (temporal) accuracy. To mitigate the problem, we have extended the existing time-dispersion correction scheme for the second- or fourth-order Lax-Wendroff method to a scheme for arbitrary even-order methods, which uses the forward and inverse time-dispersion transform (FTDT and ITDT) to add and remove the time dispersion from synthetic data. We test the correction scheme using a homogeneous model and the Sigsbee2A model. The modeling examples suggest that the use of derived FTDT and ITDT pairs on high-order Lax-Wendroff methods can effectively remove time-dispersion errors from high-frequency waves while using longer time steps than allowed in low-order Lax-Wendroff methods. We investigate the influence of the time dispersion on full-waveform inversion (FWI) and show an antidispersion workflow. We apply the FTDT to source terms and recorded traces before inversion, resulting in the source and adjoint wavefields containing equal time dispersion from source-side wave propagation and the modeled and observed traces accumulating equal time dispersion from source- and receiver-side wave propagation. The inversion results reveal that the antidispersion workflow is capable of increasing the accuracy of FWI for arbitrary even-order Lax-Wendroff methods. In addition, the high-order method can obtain better inversion results compared to the second-order method with the same antidispersion workflow.

Efficient snapshot-free reverse time migration and computation of multiparameter gradients in full-waveform inversion

Computing images in reverse time migration and model parameter gradients from adjoint wavefields in full-waveform inversion (FWI) requires the correlation of a forward-propagated wavefield with another reverse-propagated wavefield. Although in theory only two wavefield propagations are required, one forward propagation and one reverse propagation, it requires storing the forward-propagated wavefield as a function of time to carry out the correlations, which is associated with significant input/output (I/O) cost. Alternatively, three wavefield propagations can be carried out to reverse propagate the forward-propagated wavefield in tandem with the reverse-propagated wavefield. We have determined how highly accurate reverse time migrated images and FWI model parameter gradients for anisotropic elastic FWI can be efficiently computed without significant disk I/O using two wavefield propagations by means of the principle of superposition.

Estimation of micro-earthquake source locations based on full adjoint P and S wavefield imaging

SUMMARY Locating micro-earthquakes with high resolution and accuracy is a challenge for traveltime inversion, which has uncertainty on the order of a Fresnel zone (many wavelengths). We develop a wave-equation imaging method to increase resolution and reduce location errors to less than a wavelength, but requires very densely deployed receiver arrays with wide aperture and considerable computational cost. Instead of using acoustic data or direct P wave arrivals only, we use elastic multicomponent data and present a new method that uses the full P and S adjoint wavefields to image the microseismic source locations. We separate the P and S waves from the data, and extrapolate the P and S wavefields of each receiver subarray by solving the P and S adjoint wave equations in parallel. We formulate three source imaging conditions by multiplying over subarrays the adjoint P wavefield (IP), S wavefield (IS) and cross-correlated P and S wavefields (IPS). We perform numerical experiments on the highly realistic SEG SEAM4D reservoir model using surface acquisition array geometries. Results for 2-D and 3-D microseismic source estimations show clean images without noisy artefacts at shallow depths. In particular, IPS provides the highest resolution source location image, while IP is limited by the P wavelength and IS is influenced by small coda artefacts. The major-axis alignment and resolution of the source location image are determined by the hypocentral location with respect to the receiver array and illumination-angle coverage, respectively. We discuss the impacts of S-wave attenuation and frequency bandwidth on the source location images. Noise tests indicate that the imaging results are relatively insensitive to ambient noise, as is observed for the surface monitoring data. Using smoothed velocity models, the imaging results are similar to the results using the true realistically heterogeneous velocity model. The 90 per cent confidence ellipse of the source location due to Gaussian-distributed velocity errors shows a larger depth error as the source becomes deeper, while the horizontal error does not change as much.

Mesh-free radial-basis function generated finite differences and its application in full-waveform inversion

As the key point of full-waveform inversion (FWI), numerical solution of wave equation is required to simulate the forward propagated and adjoint wavefield. When traditional regular grids are used for forward modelling, scattering artifacts may occur due to the stepped approximation of layer interfaces and rugged topography. Unstructured mesh or irregular grids method can achieve certain geometric flexibility, however, its algorithm is quite complex. Mesh-free FWI can effectively reduce the scattered artifacts under regular grids and avoid the extra computation in the process of irregular grids generation. For the implementation of mesh-free FWI method, an algorithm with fast generation of node distributions is used to discretize the velocity model, radial-basis function generated finite difference is used to realise seismic wave propagation numerical simulation, and Limited-memory BFGS (L-BFGS) algorithm is used for iteration. The mesh-free FWI method we proposed achieves flexibility of simulation region and abundant wavefield information. It reduces the storage required for FWI and illustrates the applicability of high-precision reconstruction of underground velocity in the case of rugged topography, which can provide more accurate velocity information for oil and gas exploration under complex geological conditions.

Correlative least-squares reverse time migration in viscoelastic media

Standard elastic least-squares reverse time migration (LSRTM) generally ignores the subsurface attenuation effects which can lead to blurred and unfocused migration images. To correct these attenuation effects, we develop a correlative Q-LSRTM approach in viscoelastic media. We use the viscoelastic wave equation based on the standard linear solid model to simulate the subsurface attenuation. The Born modeling and adjoint state equations in viscoelastic media can be derived by the Born approximation and adjoint state method, respectively. Instead of the conventional amplitude-matching misfit function, the proposed Q-LSRTM approach utilizes the normalized zero-lag crosscorrelation misfit function to emphasize the phase measure between the predicted and observed data. The gradients of the crosscorrelation misfit function with regard to reflectivity images of P- and S-wave impedance are estimated by a zero-lag crosscorrelation between the forward background wavefields and the reversal adjoint wavefields. In view of computational efficiency, we also derive an analytical step-length formula for the correlative Q-LSRTM approach, which only requires to solve the Born modeling operator once. Numerical examples using several 2D models have determined that the proposed Q-LSRTM can improve better image quality, compared with the standard elastic LSRTM approach. Furthermore, the crosscorrelation misfit function is less sensitive to amplitude errors.

Rytov-approximation-based wave-equation traveltime tomography

Most finite-frequency traveltime tomography methods are based on the Born approximation, which requires that the scale of the velocity heterogeneity and the magnitude of the velocity perturbation should be small enough to satisfy the Born approximation. On the contrary, the Rytov approximation works well for large-scale velocity heterogeneity. Typically, the Rytov-approximation-based finite-frequency traveltime sensitivity kernel (Rytov-FFTSK) can be obtained by integrating the phase-delay sensitivity kernels with a normalized weighting function, in which the calculation of sensitivity kernels requires the numerical solution of Green’s function. However, solving the Green’s function explicitly is quite computationally demanding, especially for 3D problems. To avoid explicit calculation of the Green’s function, we show that the Rytov-FFTSK can be obtained by crosscorrelating a forward-propagated incident wavefield and reverse-propagated adjoint wavefield in the time domain. In addition, we find that the action of the Rytov-FFTSK on a model-space vector, e.g., the product of the sensitivity kernel and a vector, can be computed by calculating the inner product of two time-domain fields. Consequently, the Hessian-vector product can be computed in a matrix-free fashion (i.e., first calculate the product of the sensitivity kernel and a model-space vector and then calculate the product of the transposed sensitivity kernel and a data-space vector), without forming the Hessian matrix or the sensitivity kernels explicitly. We solve the traveltime inverse problem with the Gauss-Newton method, in which the Gauss-Newton equation is approximately solved by the conjugate gradient using our matrix-free Hessian-vector product method. An example with a perfect acquisition geometry found that our Rytov-approximation-based traveltime inversion method can produce a high-quality inversion result with a very fast convergence rate. An overthrust synthetic data test demonstrates that large- to intermediate-scale model perturbations can be recovered by diving waves if long-offset acquisition is available.

Accelerating numerical wave propagation by wavefield adapted meshes. Part II: full-waveform inversion

SUMMARYWe present a novel full-waveform inversion (FWI) approach which can reduce the computational cost by up to an order of magnitude compared to conventional approaches, provided that variations in medium properties are sufficiently smooth. Our method is based on the usage of wavefield adapted meshes which accelerate the forward and adjoint wavefield simulations. By adapting the mesh to the expected complexity and smoothness of the wavefield, the number of elements needed to discretize the wave equation can be greatly reduced. This leads to spectral-element meshes which are optimally tailored to source locations and medium complexity. We demonstrate a workflow which opens up the possibility to use these meshes in FWI and show the computational advantages of the approach. We provide examples in 2-D and 3-D to illustrate the concept, describe how the new workflow deviates from the standard FWI workflow, and explain the additional steps in detail.

Elastic Full Waveform Inversion With Source-Independent Crosstalk-Free Source-Encoding Algorithm

Elastic full waveform inversion (FWI) is more suitable to process multicomponent seismic data and can provide more subsurface medium information than acoustic FWI often with lower efficiency. Except for the parallel algorithms, source-encoding methods are usually adopted to improve the efficiency of FWI, but it often includes crosstalk noise. Besides, the additional source estimation process, critical for a successful FWI, would counteract the high-efficiency advantage of the source-encoding algorithm. We propose an elastic FWI with source-independent crosstalk-free encoding algorithm to solve the above problems. Arbitrary-phase harmonic sine functions are used as new source wavelets to perform the time-domain wavefield simulation regardless of the true wavelet. Treating the harmonic wavelet as the encoding operator and based on the orthogonality of trigonometric functions within integer periods, the amplitude and phase of each source are recovered from the blended source and adjoint wavefields so that the influence of crosstalk noise is avoided. With the deblended data, the proposed algorithm can be naturally applied to unfixed-spread acquisition systems. Moreover, we can conveniently perform the multiscale inversion by controlling the frequencies of simultaneous-source signals as conventional frequency-domain FWI does. Synthetic examples show that the proposed algorithm has high efficiency and accuracy with a strong robustness to the incorrect wavelets.

Compressed Imaging to Reduce Storage in Adjoint-State Calculations

To generate seismic images of the subsurface with adjoint-state methods such as reverse-time migration (RTM) and full-waveform inversion (FWI), the gradient of a misfit function is computed efficiently by applying what is referred to as an imaging condition to the forward propagated source wavefield and the backward propagated adjoint wavefield. In order to reduce the storage in adjoint-state calculations, we evaluate the imaging condition only at a randomly selected subset of spatial grid points (compressed imaging) and then efficiently reconstruct the full image from the imaged points via compressed sensing theory, which combines the compressibility characteristics of seismic images and convex optimization tools for the reconstruction. We use the second-order total variation regularization for the reconstruction and, using different numerical tests from RTM and FWI, we show that the new method allows a significant reduction in wavefield storage, while still recovering the full image accurately. Furthermore, regularization, applied on the gradient during the reconstruction stage improves the convergence of the FWI algorithm.

Lamb Waves Topological Imaging of Multiple Blind Defects in an Isotropic Plate

The study investigates the feasibility of the Lamb wave topological imaging method for detecting multiple blindholes in an isotropic plate. The topological imaging method is performed based on the computations of two wave fields, a forward and an adjoint, in the defect-free reference medium using different emitting sources. The image is computed by multiplying the forward and adjoint wave fields together and integrating them over time or frequency. The interferences of multimode aliasing and the scattering effect can thus be eliminated at the defectfree positions with an improved image resolution. To investigate the physical mechanism, the refocusing process of the multimode Lamb waves at the defect positions is presented by a face-to-face comparison between the snapshots of the forward and adjoint wave fields using the finite element simulation. The Lamb wave topological imaging method is numerically and experimentally verified to identify multiple blind-holes in an isotropic aluminium plate. The results demonstrate that the topological imaging method enables the suppression of the sartefacts resulting from the mode conversion and achieve high-resolution imaging of the blind defects

Elastic wave-equation migration velocity analysis preconditioned through mode decoupling

Multicomponent seismic data acquisition can reveal more information about geologic structures and rock properties than single component acquisition. Full elastic wave seismic imaging, which uses multicomponent seismic to its full potential, is promising because it provides more opportunities to understand the material properties of the earth by the joint use of P- and S-waves. A prerequisite of seismic imaging is the availability of a reliable macrovelocity model. Migration velocity analysis for P-waves, which can fill that requirement for the P-wave velocity, has been well-studied, especially under the acoustic approximation. However, a reliable estimation of the S-wave velocities remains troublesome. Elastic wave-equation migration velocity analysis has the potential to build P- and S-wave velocity models together, but it inevitably suffers from the effects of mode coupling and conversion in the forward and adjoint wavefield reconstructions. We have developed a differential semblance optimization approach to sequentially invert the background P- and S-wave velocity models from extended PP- and PS-images in the subsurface offset domain. Preconditioning of the gradients with respect to the S-wave velocity through mode decoupling can improve the reliability of the optimization. Numerical investigations with synthetic examples demonstrate the effectiveness of gradient preconditioning and the feasibility of our migration velocity analysis approach for elastic wave imaging.

Elastic model low- to intermediate-wavenumber inversion using reflection traveltime and waveform of multicomponent seismic data

Low-, intermediate-, and high-wavenumber components of P- and S-wave velocities jointly influence the elastic wave propagation and scattering in an isotropic medium. By taking advantage of all information in the data, elastic full-waveform inversion (E-FWI) has the potential to recover these model components. However, if the transmitted wave data are insufficient to illuminate the deeper part of the subsurface, we should rely on the solutions using reflection data. To reduce the nonlinearity of waveform inversion, we choose to decouple the effects of the model background and perturbation on the reflected waves within a linearized inversion framework. This resorts to three stages aiming to gradually fit the traveltimes and waveforms of the reflected PP and PS waves based on data or gradient preconditioning through P/S mode decomposition. For the first two stages, once the multicomponent seismograms have been separated into PP and PS reflection recordings, reflection traveltime inversion using an acoustic wave propagator (A-RTI) can successively recover the low-wavenumber components of P- and S-wave velocities. In the last stage, starting from the models having reliable low-wavenumber components, elastic reflection waveform inversion (E-RWI) can easily get out of the local minima and continue to retrieve the increasing wavenumber features sensitive to the waveform and amplitude variations. This is supported by gradient preconditioning through P/S mode decomposition of the extrapolated normal and adjoint wavefields, and alternately updating model background and high-wavenumber components in terms of linearized least-squares inversion. Numerical examples have demonstrated the performance of our E-RWI approach and the validity of the three-stage inversion workflow.

Frozen Gaussian approximation for 3D seismic tomography

Three-dimensional (3D) wave-equation-based seismic tomography is computationally challenging in large scales and high-frequency regime. In this paper, we apply the frozen Gaussian approximation (FGA) method to compute 3D sensitivity kernels and seismic tomography of high-frequency. Rather than standard ray theory used in seismic inversion (e.g. Kirchhoff migration and Gaussian beam migration), FGA is used to compute the 3D high-frequency sensitivity kernels for travel-time or full waveform inversions. Specifically, we reformulate the equations of the forward and adjoint wavefields for the purpose of convenience to apply FGA, and with this reformulation, one can efficiently compute the Green’s functions whose convolutions with source time function produce wavefields needed for the construction of 3D kernels. Moreover, a fast summation method is proposed based on local fast Fourier transform which greatly improves the speed of reconstruction as the last step of FGA algorithm. We apply FGA to both the travel-time adjoint tomography and full waveform inversion (FWI) on synthetic crosswell seismic data with dominant frequencies as high as those of real crosswell data, and confirm again that FWI requires a more sophisticated initial velocity model for the convergence than travel-time adjoint tomography. We also numerically test the accuracy of applying FGA to local earthquake tomography. This study paves the way to directly apply wave-equation-based seismic tomography methods into real data around their dominant frequencies.