In elastic reverse time migration (RTM), wavefield separation is an important step to remove crosstalk artifacts and improve imaging quality. State-of-the-art techniques for wavefield separation in isotropic elastic media include using the Helmholtz decomposition and introducing an auxiliary wave equation. Although these two approaches produce pure-mode vector wavefields with correct amplitudes, phases, and physical units, their computational costs are still high under current computational capability, especially for 3D large-scale problems. Based on the P- and S-wave dispersion relations, we have developed an efficient wavefield separation strategy for elastic RTM. Instead of solving a vector Poisson’s equation in the Helmholtz decomposition, we modify the phases of source wavelet as well as multicomponent records and scale the amplitudes of extrapolated wavefields with the squares of P- and S-wave velocities. This operation allows us to produce vector P- and S-wavefields with the same phases and amplitudes as the input coupled wavefields while significantly reducing computational costs. With the separated vector wavefields, we implemented a modified dot-product imaging condition for elastic RTM. In comparison with the previously proposed dot-product imaging condition, this modified imaging condition enables us to eliminate the effects of multiplication with a cosine function and hence produces migrated images with accurate amplitudes. Several 2D and 3D numerical examples are used to demonstrate the feasibility and robustness of our method for imaging complex subsurface structures.

Using adjoint-based elastic reverse time migration, it is difficult to produce high-quality reflectivity images due to the limited acquisition apertures, band-limited source time function, and irregular subsurface illumination. Through iteratively computing the Hessian inverse, least-squares migration enables us to reduce the point-spread-function effects and improve the image resolution and amplitude fidelity. By incorporating anisotropy in the 2D elastic wave equation, we have developed an elastic least-squares reverse time migration (LSRTM) method for multicomponent data from the vertically transversely isotropic (VTI) media. Using the perturbed stiffness parameters C33 and C55 as PP and PS reflectivities, we linearize the elastic VTI wave equation and obtain a Born modeling (demigration) operator. Then, we use the Lagrange multiplier method to derive the corresponding adjoint wave equation and reflectivity kernels. With linearized forward modeling and adjoint migration operators, we solve a linear inverse problem to estimate the subsurface reflectivity models for ΔC33 and ΔC55. To reduce the artifacts caused by data over-fitting, we introduce total-variation regularization into the reflectivity inversion, which promotes a sparse solution in terms of the model derivatives. To accelerate the convergence of LSRTM, we use source illumination to approximate the diagonal Hessian and use it as a preconditioner for the misfit gradient. Numerical examples help us determine that our elastic VTI LSRTM method can improve the spatial resolution and amplitude fidelity in comparison to adjoint migration.

We have developed a new approach to perform Bayesian linearized amplitude-variation-with-offset (AVO) inversion in the depth domain using nonstationary wavelets. Bayesian linearized AVO inversion, a hybrid approach combining physics-based models with statistical learning, has gained immense popularity because of its superior computational speed and ability to estimate uncertainties in inverted model parameters. Bayesian linearized AVO inversion is performed on time-domain seismic data; therefore, depth-imaged seismic cannot be inverted directly using this method and would require depth-to-time conversion before AVO inversion can be done. Subsequently, time-to-depth conversion of the inverted volumes would be required for reservoir modeling and well placement. Domain conversions introduce additional uncertainty in geophysical workflows. In conventional AVO inversion, the seismic wavelet is assumed to be stationary, and this assumption leads to a restriction in the length of the time window over which the inversion can be performed. Therefore, AVO inversion is usually restricted to a narrow time window around the target of interest, and if multiple targets are present at different depths, multiple inversions must be run on the same volume. Depth-domain amplitude inversion is a recent development and has been previously presented in an iterative formulation. Implementing linearized Bayesian inversion directly in the depth domain using nonstationary wavelets is a convenient new approach that takes advantage of superior computational speed and uncertainty quantification without compromising the accurate spatial location that depth imaging provides. Combining these two ideas creates a novel, unique, and powerful seismic inversion technique that can be useful for quantitative interpretation and reservoir characterization.

Finite-frequency traveltime inversion offers higher accuracy for velocity model building than ray-based traveltime inversion. The adjoint force is the key for the computation of inversion kernels. Starting at the definition of inversion kernels for the acoustic wave equation, we have derived the explicit formula for the spectral distribution density function used in the adjoint force computation. Two formulations are developed for the computation of adjoint forces for receiver-side extrapolation, frequency-domain representation, and time-domain representation. The accuracy of finite-frequency traveltime inversion kernels is benchmarked with the analytical solutions for homogeneous isotropic media. We use wavefront construction to compute the first Fresnel zones for kernel conditioning. Based on dynamic ray tracing, we design a processing procedure guided by synthetic data tests to extract the desired events for wavefield backward extrapolation from the data. Unlike ray-based velocity tomography, finite-frequency inversion can resolve the velocity structures comparable with the size of Fresnel zones as we demonstrate on a marine salt model using ocean-bottom node acquisition geometry. Despite the fact that the inversion kernels are based on Born approximation, velocities with errors up to 20% can be well-resolved. For practical purposes, a simple formulation is given for the determination of the shot spacing. Our workflow for finite-frequency inversion is efficient and converges in only very few iterations.