Anisotropic and attenuating properties of subsurface media cause amplitude loss and waveform distortion in seismic wave propagation, resulting in negative influence on seismic imaging. To correct the anisotropy effect and compensate amplitude attenuation, a compensated-amplitude vertical transverse isotropic (VTI) least-squares reverse time migration (LSRTM) method is adopted. In this method, the attenuation term of an attenuated acoustic wave equation is extended to a VTI quasi-differential wave equation, which takes care of effects from anisotropy and attenuation. The finite-difference method is used to solve the equation, in which attenuation terms are solved in the wavenumber domain, and other terms are solved in the space or wavenumber domain. Stable regularization operators are derived and introduced to the equations to eliminate severe numerical noise in high-frequency components during backward propagation. Meanwhile, a demigration operator, migration operator, and gradient formula for attenuated VTI media are derived to implement the amplitude-compensated VTI LSRTM. Test of a homogeneous model proves the accuracy of the attenuated VTI quasi-differential equations and the effectiveness of the regularization operators. A numerical example for a modified Marmousi model verifies the accuracy and superiority to the amplitude-compensated VTI LSRTM. Our results show that the sensitivity to anisotropic parameters is much higher than that to the $Q$ parameters.

Elastic full-waveform inversion (FWI) provides a better description of the subsurface information than those given by the acoustic assumption. However, it suffers from a more serious cycle-skipping problem compared with the latter. Reflection waveform inversion (RWI) is able to build a good background model, which can serve as an initial model for elastic FWI. Because, in RWI, we use the model perturbation to explicitly fit reflections, such perturbations should include density, which mainly affects the dynamics. We applied Born modeling to generate synthetic reflection data using optimized perturbations of the P- and S-wave velocities and density. The inversion for the perturbations of the P- and S-wave velocities and density is similar to elastic least-squares reverse time migration. An incorrect background model will lead to misfits mainly at the far offsets, which can be used to update the background P- and S-wave velocities along the reflection wavepath. We optimize the perturbations and background models in an alternate way. We use two synthetic examples and a field-data case to demonstrate our proposed elastic RWI algorithm. The results indicate that our elastic RWI with variable density is able to build reasonably good background models for elastic FWI with the absence of low frequencies, and it can deal with the variable density, which is required in real cases.

This study derives a preconditioned stochastic conjugate gradient (CG) method that combines stochastic optimization with singular spectrum analysis (SSA) denoising to improve the efficiency and image quality of plane-wave least-squares reverse time migration (PLSRTM). This method reduces the computational costs of PLSRTM by applying a controlled group-sampling method to a sufficiently large number of plane-wave sections and accelerates the convergence using a hybrid of stochastic descent (SD) iteration and CG iteration. However, the group sampling also produces aliasing artifacts in the migration results. We use SSA denoising as a preconditioner to remove the artifacts. Moreover, we implement the preconditioning on the take-off angle-domain common-image gathers (CIGs) for better results. We conduct numerical tests using the Marmousi model and Sigsbee2A salt model and compare the results of this method with those of the SD method and the CG method. The results demonstrate that our method efficiently eliminates the artifacts and produces high-quality images and CIGs.

The traditional high-order finite-difference (FD) methods approximate the spatial derivatives to arbitrary even-order accuracy, whereas the time discretization is still of second-order accuracy. Temporal high-order FD methods can improve the accuracy in time greatly. However, the present methods are designed mainly based on the acoustic wave equation instead of elastic approximation. We have developed two temporal high-order staggered-grid FD (SFD) schemes for modeling elastic wave propagation. A new stencil containing the points on the axis and a few off-axial points is introduced to approximate the spatial derivatives. We derive the dispersion relations of the elastic wave equation based on the new stencil, and we estimate FD coefficients by the Taylor series expansion (TE). The TE-based scheme can achieve ([ ])th-order spatial and ([ ])th-order temporal accuracy ([ ]). We further optimize the coefficients of FD operators using a combination of TE and least squares (LS). The FD coefficients at the off-axial and axial points are computed by TE and LS, respectively. To obtain accurate P-, S-, and converted waves, we extend the wavefield decomposition into the temporal high-order SFD schemes. In our modeling, P- and S-wave separation is implemented and P- and S-wavefields are propagated by P- and S-wave dispersion-relation-based FD operators, respectively. We compare our schemes with the conventional SFD method. Numerical examples demonstrate that our TE-based and TE + LS-based schemes have greater accuracy in time and better stability than the conventional method. Moreover, the TE + LS-based scheme is superior to the TE-based scheme in suppressing the spatial dispersion. Owing to the high accuracy in the time and space domains, our new SFD schemes allow for larger time steps and shorter operator lengths, which can improve the computational efficiency.

Finite difference (FD) is widely used for modeling seismic‐wave propagation in theoretical and applied seismology. Traditional FD methods with the cross stencil exhibit spatial general‐order and temporal second‐order accuracy, causing strong temporal dispersion and poor stability. We propose two novel time–space domain staggered‐grid FD (SFD) schemes, in which the spatial derivatives are approximated through a weighted summation of function values from the axial and some off‐axial nodes. We derive the dispersion relation of acoustic wave equation and estimate FD coefficients by the Taylor series expansion (TE). The TE‐based SFD scheme with the new stencil can obtain ($2N$)th‐order accuracy in both time and space domains along all the propagation directions for $N≤4$. In contrast, it reaches ($2N$)th‐order accuracy along four propagation directions and eighth‐order accuracy along other propagation directions for $N>4$. Furthermore, we investigate an SFD scheme using a combination of TE and least squares (LS). FD coefficients at the off‐axial nodes are computed by TE, and then FD coefficients at the axial nodes are optimized by LS. We compared our schemes with the conventional, the temporal fourth‐ and sixth‐order, and several optimal SFD schemes. Dispersion analysis, stability analysis, and 2D and 3D numerical examples show that our TE‐based SFD scheme has smaller temporal dispersion and better stability than other schemes. Besides, our TE+LS‐based SFD scheme can further suppress spatial dispersion with slightly degraded stability. We also analyze the efficiency of different schemes. Our new SFD schemes allow for shorter operator lengths and larger time steps owing to higher accuracy and less stringent stability conditions. Numerical results reveal that the computational cost of our TE+LS‐based scheme is 49.1% and 85.7% of that of the temporal fourth‐order and optimal SFD schemes, respectively, under the same accuracy demand.

Gaussian-beam migration (GBM) is flexible and adaptable for imaging geologically complicated areas. It avoids some limitations, such as amplitude singularity in caustic zones and an inability to image multiple arrivals, of traditional ray-based migration approaches. Previous studies on GBM mainly focused on acoustic media. We have developed a 2D isotropic elastic GBM scheme for common-shot multicomponent records. Our method extrapolates P- and S-mode wavefields simultaneously using the Kirchhoff-Helmholtz integral solution of the isotropic elastodynamic equation. We separate the extrapolated wavefields into compressional and shear modes using the Helmholtz decomposition, and we then implement a modified dot-product imaging condition. This approach enables us to produce clear PP-images and avoid polarity reversal issues for PS-images. In addition, based on the theory of wavefield approximation in the effective vicinity of central rays, we derive a formula to compute the propagation angles of paraxial rays, which can be used to extract angle-domain common-image gathers.