Misfit Function Research Articles

On the adjoint state method for the gradient computation in full waveform inversion: a complete mathematical derivation for the (visco-)elastodynamics approximation

Summary High resolution seismic imaging at all scales using full waveform inversion is now routinely used in the industry and in the academy. One key element for the success of this approach is a numerical method, named adjoint state method, originally designed for optimization problems constrained by partial differential equations, a category to which full waveform inversion belongs. This method provides an efficient way to compute the gradient of the full waveform inversion misfit function, which is the most computationally demanding task in the implementation of full waveform inversion. While well known, the complete and rigorous mathematical derivation of the adjoint state method for full waveform inversion remains missing in the scientific bibliography. The aim of this study is to remedy this lack. The derivation is performed in general settings, that is in the elastodynamics approximation, with and without considering viscosity. Through the calculus, the mechanism of the adjoint state strategy makes clear the connection between the incident and adjoint fields, especially regarding their initial and boundary conditions. The impact of introducing the viscosity is carefully analyzed. The resulting gradient formulas are analyzed and shown to be consistent with already published ones. The generic approach which is adopted also makes it possible to derive misfit function gradients with respect to other quantities than the subsurface mechanical parameters, for instance with respect to the initial or the boundary conditions, which could be of interest for specific applications where the reconstructed parameters are not only volumetric mechanical parameters but boundary parameters or initial field values.

Optimal transport assisted full waveform inversion for multiparameter imaging of soft tissues in ultrasound computed tomography

Ultrasound computed tomography (USCT) has emerged as a promising platform for imaging tissue properties, offering non-ionizing and operator-independent capabilities. In this work, we demonstrate the feasibility of obtaining quantitative images of multiple acoustic parameters (sound speed and impedance) for soft tissues using full waveform inversion (FWI), which are justified with both numerical and experimental cases. A 3D reconstruction based on a series of 2D slice images is presented for the experimental case of ex vivo soft tissues. To improve the robustness of the reconstruction process, a hierarchical FWI strategy is adopted, gradually iterating from low to high frequencies. In parallel, we employ a graph-space optimal transport misfit function, avoiding convergence into local minima and minimizing inversion artifacts caused by skin-related supercritical reflections. Our method first carries out sound speed inversion based on transmitted waves in the low and middle frequency bands, and then uses all types of waves in the high frequency band for simultaneous inversion of both sound speed and impedance. Compared to conventional strategies, the proposed approach can accurately reconstruct physical models consistent with the actual soft tissue sample. These high-resolution ultrasound images of acoustic parameters are promising to allow for quantitative differentiation among different types of tissues (e.g., muscles and fats). These results have significant implications for advancing our understanding of tissue properties and for potentially contributing to disease diagnosis through USCT, which is a flexible and cost-effective alternative to X-ray computed tomography or magnetic resonance imaging at no significant sacrifices for resolution.

Towards a new standard for seismic moment tensor inversion containing 3-D earth structure uncertainty

SUMMARY Moment tensor (MT) inversion is a classical geophysical inverse problem that infers a force-equivalent model of a seismic source from seismological observations. Like other inverse problems, the accuracy of the inversion depends on the reliability of the forward problem simulating waveforms from the source location through an Earth structural model. Apart from errors in data, the error in forward waveform simulation, also known as theory error, is a significant source of error contributing to the misfit function between the predicted and observed waveforms. Here, we set up numerical experiments to comprehensively probe the sensitivity of the linearized MT inversion to 3-D regional earth model errors, a known predominant factor of the theory error. Using the Monte Carlo method, we estimate the empirical structural covariance matrices to characterize the waveform mismatch due to the imperfect knowledge of Earth's structure. First, although the inversion accuracy deteriorates with increasing model errors, incorporating the structural covariance matrices into the misfit function improves the accuracy of inversion results for all theorized error distributions. Secondly, we propose a slightly modified form of the structural covariance matrix, which further enhances the inversion outcome. Lastly, as the true structural errors are likely spatially correlated, we highlight the importance of adequately treating the correlation into the MT inversion because of its significant impact on inversion. Overall, as a preliminary effort in quantifying 3-D structural errors on MT inversion, this study proves the computational feasibility by means of numerical experiments and will hopefully provide a way forward for future work on this topic.

Optimal Transport Map With Prescribed Direction Indicator for Seismic Full‐Waveform Inversion

AbstractThe quadratic Wasserstein (W2) metric has been proposed as a promising misfit function to mitigate cycle‐skipping phenomena in full‐waveform inversion. Mathematically, we demonstrate that the smoothness of the W2‐based adjoint source is two orders of magnitude higher than that based on L2‐norm, which guarantees a larger convergence radius of related inverse problems. However, the oscillatory characteristics of seismic signals and subsequent operations of transforming them into probability densities would decrease the accuracy of the optimal transport map T(t) and exacerbate the nonconvexity of the misfit function. To tackle these challenges, we propose the concept of prescribed direction indicator, which indicates the properly matching direction from predictions to observations, in order to correct inaccurate T(t). 1D synthetic examples suggest that reasonable bijection can be constructed through the proposed method. Numerical experiments demonstrate that it works well during optimization procedures, including enlarging the convergence radius of the inverse problem, improving the computational efficiency and enhancing the reliability of inversion results.

Inversion strategies for Q estimation in viscoacoustic full-waveform inversion

Estimation of an attenuation parameter, represented by the quality factor Q, holds paramount importance in seismic exploration. One of the main challenges in Q estimation through viscoacoustic full-waveform inversion (FWI) is effectively decoupling Q from velocity. In this study, our objective is to enhance Q inversion by addressing critical aspects, such as gradient preconditioning, workflow, and misfit design. By developing a new preconditioner that approximates the diagonal of the Hessian, we facilitate automatic parameter tuning across different classes, ensuring comparable magnitudes of preconditioned gradients for velocity and Q. Moreover, our investigations confirm the efficacy of the two-stage hierarchical strategy in mitigating velocity- Q trade-offs, enabling more accurate Q estimation by first focusing on velocity reconstruction before jointly estimating velocity and Q. The analysis and numerical examples also highlight the importance of broadband data and long-offset acquisition for a reliable Q estimation. In addition, leveraging amplitude information can improve Q estimation to some extent, but careful consideration of frequency band and noise effects is necessary. We explore two misfit functions that capture amplitude variation with frequency in the time-frequency domain, noting their sensitivity to noise. To address this, we develop a differential strategy that can effectively mitigate the effects of low-frequency noise. This comprehensive study on enhancing Q estimation in viscoacoustic FWI offers valuable insights for multiparameter inversion in realistic scenarios.

Gradient-based joint inversion of point-source moment tensor and station-specific time-shifts

SUMMARY The misalignment of the observation and predicted waveforms in regional moment tensor inversion is mainly due to seismic models’ incomplete representation of the Earth's heterogeneities. Current moment tensor inversion techniques, allowing station-specific time-shifts to account for the model error, are computationally expensive. Here, we propose a gradient-based method to jointly invert moment-tensor parameters, centroid depth and unknown station-specific time-shifts utilizing the modern functionalities in deep learning frameworks. A $L_2^2$ misfit function between predicted synthetic and time-shifted observed seismograms is defined in the spectral domain, which is differentiable to all unknowns. The inverse problem is solved by minimizing the misfit function with a gradient descent algorithm. The method's feasibility, robustness and scalability are demonstrated using synthetic experiments and real earthquake data in the Long Valley Caldera, California. This work presents an example of fresh opportunities to apply advanced computational infrastructures developed in deep learning to geophysical problems.

Multimodal surface wave inversion with automatic differentiation

SUMMARY Investigating subsurface shear wave velocity (vs) structures using surface wave dispersion data involves minimizing a misfit function that is commonly solved through gradient-based optimization. Sensitivity kernels for model updates are commonly estimated using numerical differentiation, variational methods or implicit functions which however, may involve numerical instability and computational challenges when dealing with complex velocity models and large data sets. In this study, we propose a novel surface wave inversion framework in which error-free gradients are calculated by automatic differentiation (AD) and forward modelling is implemented by convenient computational graphs in the state-of-the-art deep learning framework. The AD-based inversion approach is first validated using two synthetic data sets. Then, the subsurface structures at three distinct locations, namely the Great Plains and the Long Beach in the US and Tong Zhou in China, are also derived using this method with seismic ambient noise data, which show nice consistency with those obtained using traditional methods. With the significantly improved computational efficiency, a great number of initial models can be inverted simultaneously to mitigate the impact of local minima and to estimate the uncertainty in the invert models. We have developed a new surface wave inversion package named ADsurf based on automatic differentiation and computational graphs in the deep learning framework, and its computational efficiency is also compared with the traditional finite-difference-based gradient estimation approach. While a great number of intriguing studies on the geophysical inverse problems have been conducted recently using deep learning for end-to-end mapping, the use of AD provided in the in the deep learning frameworks to assist and expedite the gradient computations are still underexploited in geophysics. Thus, it is expected that various geophysical inverse problems in many different areas beyond the surface wave inversion can also be tackled with this new paradigm in the future.

DG-based joint transmission-reflection traveltime tomography and its application of borehole seismic

Abstract The limitations of the coverage range and density of transmission wave often result in less-than-ideal results in traveltime tomography. In contrast, joint transmission-reflection traveltime tomography can not only recover deep structures that transmission tomography cannot detect but also optimize its inversion results. In this article, we perform joint tomography on borehole seismic (VSP, RVSP and crosswell seismic) data to obtain near-wellbore structures. In the forward part, we solve the factored equation by the discontinuous Galerkin (DG) method to calculate the transmission/reflection traveltime. Due to the large wavefront curvature near the source point, the traveltime errors generated by the numerical simulation will propagate from the source to all the calculation domains. According to the factorization principle, the equation solution is decomposed into two parts to solve the point-source singularity. To further improve the accuracy of solving traveltime, we use the DG method to solve the factored eikonal equation with additive factors (the factored DG method), obtaining second-order accuracy solution. The adjoint-state method is employed in the inversion section to calculate the gradient of the misfit function. And we use the traveltime difference observed inside the model to define the misfit function, which is more suitable for borehole seismic and avoids the influence of surface normal vectors on gradients. Numerical tests applied on models indicate that the joint tomography method has the potential to accurately inverse the seismic structure information near the well and recover the deep underground structure.

Novel Phase-Sensitive Full-Waveform Tomography for Seismic Imaging

Abstract Full-waveform tomography (FWT) is increasingly recognized as a pivotal technique for delineating high-resolution subsurface properties. Despite its significant potential, practical applications of FWT encounter persistent challenges, particularly in dealing with local minima and cycle-skipping problems. These difficulties often arise and are intensified by the least-squares (L2) norm’s intrinsic insensitivity to phase mismatches. To address these challenges, we have redefined the traditional L2 norm misfit function by incorporating a time shift within the synthetic waveform. This shift is determined by the temporal discrepancies between the observed and synthetic waveforms, identified through a cross-correlation technique. This approach, termed phase-sensitive FWT, integrates phase differences into the new misfit function, thus significantly mitigating the cycle-skipping problem. Numerical experiments demonstrate that PSFWT reduces dependence on the initial model and achieves more accurate inversion results compared with the traditional L2 norm method, highlighting its potential for enhancing the precision and reliability of seismic imaging.

Deep-unrolling architecture for image-domain least-squares migration

Deep-image prior (DIP) is a novel approach to solving ill-posed inverse problems whose solution is parameterized with an untrained deep neural network and cascaded with the forward modeling operator. A key component to the success of such a method is represented by the choice of the network architecture, which must act as a natural prior to the inverse problem at hand and provide a strong inductive bias toward the desired solution. Inspired by the close link between neural networks and iterative algorithms in classical optimization, we apply an unrolled version of the gradient descent (GD) algorithm as our DIP network architecture, denoted as the deep-unrolling (DU) architecture. Each layer of the unrolled network is comprised of two parts: the first part corresponds to the GD step of the data-fidelity term, whereas the second part, formed by a six-layer convolutional neural network (CNN), plays the role of a regularizer in the associated objective function. The developed DU architecture is applied to the problem of image-domain least-squares migration (IDLSM) to invert migrated seismic images for their underlying reflectivity and is denoted as DU-IDLSM. As such, the DU architecture parameterizes the reflectivity, and the input of each layer of the unrolled network is the reflectivity at the previous layer. Similar to the classical DIP approach, the parameters of the DU architecture are optimized in an unsupervised fashion by minimizing the data misfit function itself. Through experiments with a part of the Sigsbee2A model and a marine field data set, we test the effectiveness of the DU-IDLSM approach and highlight two key benefits. First, the DU architecture can effectively regularize the inversion process, resulting in reflectivity estimates with fewer artifacts and higher image resolution than those produced by conventional IDLSM approaches. Second, we indicate that by including dropout layers in the CNN architecture, DU-IDLSM can produce a qualitative measure of the uncertainty associated with the least-squares migration process.

Ocean bottom dual-sensorQ-compensated elastic least-squares reverse time migration based on acoustic and separated-viscoelastic coupling equations

The ocean bottom dual-sensor (OBD) observation system places water detectors (pressure sensors) and land detectors on the irregular seabed interface to receive acoustic pressure fields and elastic wavefields, respectively. These data are usually applied to effectively remove ghost waves and ringing interference in the seawater layer, rather than being jointly used for migration. We propose an ocean bottom dual-sensor viscoacoustic and separated-viscoelastic coupled Q-compensated least-squares reverse time migration (OBD-VSV- QLSRTM) in curvilinear coordinates. In this method, to solve the Q-compensated elastic inverse problem, we jointly use the recorded acoustic pressure and elastic multicomponent seismic data from the OBD observation system. A new misfit function in the sense of least-squares inversion is constructed based on viscoacoustic and separated viscoelastic wavefields. Acoustic and separated viscoelastic equations are used to describe the acoustic wavefield in the acoustic medium and the Q-compensated P and S wavefields in the underlying viscoelastic medium, respectively. On the acoustic-viscoelastic interface, we implement a stable governing equation to ensure continuous and steady transmission of Q-compensated stresses in the subseabed and pressure in the seawater. To address the irregular seabed problem, we mesh the acoustic-viscoelastic models onto curvilinear grids and apply coordinate transformations to map the models and equations to a new curvilinear coordinate system. We derive viscoacoustic and separated-viscoelastic coupled adjoint and demigration operators in the curvilinear coordinates, enabling linear acoustic-viscoelastic modeling and stable receiver-side back-propagated wavefield continuation. Numerical examples conducted on two typical acoustic-viscoelastic models demonstrate that our proposed OBD-VSV- QLSRTM in curvilinear coordinates can produce highly accurate P- and S-wave images efficiently.

Robust and efficient waveform-based velocity model building by optimal transport in the pseudotime domain: An ocean-bottom cable case study in the North Sea

Full-waveform inversion (FWI) in the North Sea has demonstrated its imaging power starting from low-resolution models obtained by traveltime tomography, enriching them with geologically interpretable fine-scale details. However, building a traveltime-based kinematically accurate starting model for FWI is a time-consuming and rather subjective process requiring phase identification and selection. The two main problems faced by FWI starting from noninformative initial models are the susceptibility to cycle skipping and a lack of sensitivity to low wavenumbers in the deep subsurface not sampled by turning waves. On a North Sea ocean-bottom cable 3D data set, a novel [Formula: see text] building methodology is applied that addresses those issues by jointly inverting reflections and refractions (joint full-waveform inversion [JFWI]) using a robust misfit function in the vertical traveltime domain (pseudotime). Pseudotime addresses reflectivity-velocity coupling and attenuates phase ambiguities at short offsets, whereas a graph-space optimal transport (GSOT) objective function with dedicated data windowing averts cycle skipping at intermediate-to-long offsets. A fast and balanced reflectivity reconstrution is obtained prior to JFWI thanks to an asymptotic-preconditioned impedance waveform inversion ([Formula: see text]WI). Starting from a linearly increasing one-dimensional model, GSOT-pseudotime JFWI is effective at obtaining a meaningful P-wave velocity macromodel down to depths sampled by reflections only, without phase identification and picking. P-wave FWI, starting from the JFWI-based model, injects the high wavenumbers missing in the JFWI solution, attaining apparent improvements in shallow and deep model reconstruction and imaging compared with the previous studies in the literature, and a satisfactory prediction of the ground-truth logs.

CO2 Storage Monitoring via Time-Lapse Full Waveform Inversion with Automatic Differentiation.

In the field of CO2 capture utilization and storage (CCUS), recent advancements in active-source monitoring have significantly enhanced the capabilities of time-lapse acoustical imaging, facilitating continuous capture of detailed physical parameter images from acoustic signals. Central to these advancements is time-lapse full waveform inversion (TLFWI), which is increasingly recognized for its ability to extract high-resolution images from active-source datasets. However, conventional TLFWI methodologies, which are reliant on gradient optimization, face a significant challenge due to the need for complex, explicit formulation of the physical model gradient relative to the misfit function between observed and predicted data over time. Addressing this limitation, our study introduces automatic differentiation (AD) into the TLFWI process, utilizing deep learning frameworks such as PyTorch to automate gradient calculation using the chain rule. This novel approach, AD-TLFWI, not only streamlines the inversion of time-lapse images for CO2 monitoring but also tackles the issue of local minima commonly encountered in deep learning optimizers. The effectiveness of AD-TLFWI was validated using a realistic model from the Frio-II CO2 injection site, where it successfully produced high-resolution images that demonstrate significant changes in velocity due to CO2 injection. This advancement in TLFWI methodology, underpinned by the integration of AD, represents a pivotal development in active-source monitoring systems, enhancing information extraction capabilities and providing potential solutions to complex multiphysics monitoring challenges.

Viscoacoustic least squares reverse-time migration using L1-2 norm sparsity constraint

Abstract Least-squares reverse-time migration (LSRTM) has become an advanced technique for complex structures imaging of the subsurface, as it can provide a higher resolution and more balanced amplitude migrated image than conventional reverse-time migration (RTM). However, the intrinsic attenuation of subsurface introduces amplitude attenuation and phase dispersion of seismic wavefield, which leads to the inverted image kinematically and dynamically inexactitude. Moreover, the imperfect geometry, limited bandwidth of seismic data, and inappropriate modeling kernel etc., would inevitably introduce two side-effects in migrated image, resulting in degradation of LSRTM imaging potential. To alleviate above issues, we present a data-domain sparsity constraint viscoacoustic least-squares reverse-time migration algorithm in this paper. In particular, we utilize the decoupled constant Q fractional Laplacians (DFLs) viscoacoustic wave equation as the modeling kernel to describe the attenuation effects of the subsurface, while a model constraint constructed in the misfit function via L1-2 norm is carried out to clear the migrated artefacts and boost the imaging resolution. Thanks to the excellent performance in sparsity, the drawbacks of unconstraint LSRTM can be effectively mitigated by the L1-2 norm-based regularization. In this paper, we adopt the alternating direction of multipliers method (ADMM) to iteratively address the constrained L1-2 minimization problem by implementing a proximal operator, and three synthetic examples are hired to evaluate the effectiveness and practicability of the proposed strategy. Migration results prove that the proposed scheme can effectively compensate the attenuation effects, improve the resolution, and suppress the migration artifacts of inverted images even in the complex imaging situations.

Constrained nonlinear amplitude-variation-with-offset inversion for reservoir dynamic changes estimation from time-lapse seismic data

We develop a novel elastic amplitude-variation-with-offset (AVO) inversion process to estimate the fluid saturation and effective pressure or variations in these properties from time-lapse seismic data sets. These changes occur in oil and gas reservoirs caused by fluid injection or hydrocarbon production that leads to changes in the elastic wave properties, reflectivity, and seismic response. Our method is based on a seismic forward model that consists of a linearized AVO equation and a rock-physics model. The AVO equation links the elastic wave properties to seismic reflection amplitudes, whereas the rock-physics model maps the saturation and pressure into seismic velocities and density. The inversion approach relies on the gradient descent technique to estimate the unknown variables by searching for the minimum of the least-squares misfit between the observed and modeled data. The first-order gradient equations of the least-squares data misfit function with respect to effective pressure and water saturation are derived by using the adjoint-state method and the chain rule. The optimization method used to minimize the misfit function and obtain the best optimal solution is a limited-memory quasi-Newton algorithm. This inversion process allows us to incorporate prior constraints by using the logistic function to map the model variables to a bounded range. To achieve a stable solution to ill-posed inversion problems, optimal regularization weights are applied. The application of the developed workflow on 1D synthetic and real well-log data from the Edvard Grieg oil field simulating various saturation-pressure conditions during production demonstrates the validity of the approach with different noise levels. The inversion is then applied to a 2D synthetic data set modeled from the reservoir model of the Smeaheia field, a potential site for a large-scale offshore CO2 storage field located in the North Sea. Our results illustrate that our inversion method efficiently and accurately estimates reservoir saturation and pressure variations.

Rayleigh Wave Dispersion and Inversion for Shallow Surface with a High-velocity Rigid Pavement

We investigate the rigid pavement influence on phase velocity dispersion of the Rayleigh wave by using numerical modeling. An overestimation of Rayleigh wave phase velocity is observed more than 10% deviation caused by high-velocity rigid pavement. Based on the generalized reflection and transmission (R/T) coefficient method, the linear eigenvalue problem is solved for vertically heterogeneous media with a rigid pavement structure. The secular function demonstrates the multi-valuedness which is different phase velocity at the same frequency for the single Rayleigh mode. To address mode misidentification problem, we adopt a new misfit function for multimodal dispersion inversion. A Markov chain Monte Carlo (MCMC) sampling algorithm is used to solve dispersion inversion. The synthetic test shows that Rayleigh wave dispersion inversion incorporating a pavement layer can improve estimation accuracy of shear wave velocity.

Weighted Adjoint-state Crosshole Traveltime Tomography

The adjoint-state method of first arrival traveltime tomography, featured with implicit ray tracing calculation, was used to compute the gradient of the misfit function without introducing Fréchet derivatives. However, the majority of velocity gradients in the adjoint-state method focus on the high velocity area due to high raypath density. Thus, the current adjoint-state method faces challenges of imaging low velocity anomalies in the near-surface by crosshole seismic acquisition. In this paper, we proposed a weighted adjoint-state Eikonal equation-based crosshole traveltime tomography method. A preconditioning operator is provided to fix the ray path density influence. The weighted method was tested on a homogeneous model, a checkerboard model, and a reduced-scale Marmousi model. Both misfit cures and velocity logs are used to make a comparison between the weighted method and the traditional method. The results show the weighted method has better resolution, faster convergence speed, and lower final residual, which demonstrate the method can mitigate the ray path effect and reduce the inverted velocity error of low velocity anomalies.

Statistical inversion of normal‐mode interference spectral parameter in ocean waveguide

AbstractA statistical inversion method of estimating normal mode interference spectral parameters (including waveguide invariance and interference frequency) from one‐dimensional broadband sound intensity spectrum (1D BSIS) is presented. First, the interference phase of 1D BSIS is compensated by a non‐linear factor for each pair of unknown parameters in a prior parameter space. Second, the singular value decomposition is performed for a data matrix, which is constructed by the phase‐compensated BSIS. Then, a misfit function is defined based on both the maximum singular value and the corresponding singular vector of the data matrix. Finally, the posterior probability distributions of unknown parameters are analysed based on statistical inversion theory. This method is suitable for the situation using a single receiver and without any environmental information. Numerical simulations illustrate the validity of the presented method.

Resolution and trade-offs in global anelastic full-waveform inversion

Summary Improving the resolution of seismic anelastic models is critical for a better understanding of the Earth’s subsurface structure and dynamics. Seismic attenuation plays a crucial role in estimating water content, partial melting, and temperature variations in the Earth’s crust and mantle. However, compared to seismic wavespeed models, seismic attenuation tomography models tend to be less resolved. This is due to the complexity of amplitude measurements and the challenge of isolating the effect of attenuation in the data from other parameters. Physical dispersion caused by attenuation also affects seismic wavespeeds, and neglecting scattering/defocusing effects in classical anelastic models can lead to biased results. To overcome these challenges, it is essential to account for the full 3D complexity of seismic wave propagation. Although various synthetic tests have been conducted to validate anelastic Full-Waveform Inversion (FWI), there is still a lack of understanding regarding the trade-off between elastic and anelastic parameters, as well as the variable influence of different parameter classes on the data. In this context, we present a synthetic study to explore different strategies for global anelastic inversions. To assess the resolution and sensitivity for different misfit functions, we first perform mono-parameter inversions by inverting only for attenuation. Then, to study trade-offs between parameters and resolution, we test two different inversion strategies (simultaneous and sequential) to jointly constrain the elastic and anelastic parameters. We found that a sequential inversion strategy performs better for imaging attenuation than a simultaneous inversion. We also demonstrate the dominance of seismic wavespeeds over attenuation, underscoring the importance of determining a good approximation of the Hessian matrix and suitable damping factors for each parameter class.

Staining algorithm for least-squares reverse time migration

The staining algorithm-based reverse time migration (RTM) technique has shown promise in improving imaging results for certain complex structures, such as subsalt and steeply dipping layers, by incorporating the complex domain. However, RTM suffers from low-resolution images as it utilizes the adjoint of the forward modeling operator rather than its inverse operator. To address this limitation, we propose a new approach called stained least-squares RTM (ST-LSRTM) that enhances the gradient through the integration of the staining algorithm. In ST-LSRTM, we first develop the wave equation and Born approximation, transitioning from the real to the complex domains. To simplify the linear perturbation operator, we constrain the magnitude of the imaginary velocity. Next, we compute the gradient in the complex domain. The enhanced gradient is obtained by replacing the traditional gradient with the normalized stained gradient within the stained region. Finally, we construct a misfit function based on the L2 norm and calculate the step size and gradient direction using the conjugate gradient method. Model tests demonstrate that ST-LSRTM achieves high-resolution imaging results even when the stained area is imprecise. Comparisons of waveform curves and normalized residual curves confirm that ST-LSRTM approaches true reflection coefficient model more closely compared to traditional LSRTM.