Reverse Time Migration (RTM)

• RTM is a computational imaging approach that uses forward and backpropagated wavefields via the two-way wave equation to recover subsurface reflectivity with high resolution.
• It employs advanced numerical methods such as high-order finite-difference, discontinuous Galerkin, and Fourier-integral techniques, enhancing scalability and artifact suppression.
• Modern RTM integrates uncertainty quantification and machine learning techniques to enable efficient, target-oriented inversion and true-amplitude recovery.

Reverse Time Migration (RTM) is a high-fidelity computational imaging method for reconstructing perturbations within a background wave-propagating medium, foundational in seismic, ultrasonic, electromagnetic, and acoustic inverse methodologies. RTM is mathematically grounded in the two-way wave equation and forms the analytic and algorithmic core for linearized inverse scattering and least-squares migration, supporting quantitative amplitude analysis and microlocal resolution in complex media. Modern developments have extended RTM to target-oriented, uncertainty-aware, and hybrid inversion contexts, leveraging efficient numerical implementations and increasing integration with machine learning.

1. Mathematical Formulation and Operator Structure

RTM seeks to recover the spatial distribution of reflectivity $r(x)$ (or more general perturbation models such as $\delta m(x)$) in a domain where a background wave equation with speed $c(x)$ is known or estimated. The core forward model is the two-way wave equation, typically for acoustic pressure $p(x,t)$: $$\frac{\partial^2 p(x,t)}{\partial t^2} = c^2(x)\nabla^2 p(x,t) + s(x,t)$$ where $s(x, t)$ is a source term, often a point impulsive or band-limited excitation. RTM is built by propagating wavefields both forward (source field, $u_s$) and backward (receiver field, $u_r$) in time, and correlating them at every spatial location: $i_C(x) = \int_{t} u_s(x, t) u_r(x, t) dt$ This "zero-lag cross-correlation" is the classical imaging condition and is equivalent to the application of the adjoint of the linearized (Born) scattering operator. However, direct inversion via the normal operator $N=F^* F$ leads to artifacts because $N$ is highly ill-conditioned at low wavenumber.

Op 't Root, Stolk, and de Hoop introduce a factorized operator approach to RTM-based inverse scattering, establishing the operator chain: $i(x) = (G V_1 P F r)(x), \qquad R = G V_1 P F$ Here,

• $F$ is the forward (scattering) operator mapping reflectivity to data,
• $P$ is a revert operator encoding pseudodifferential (one-way) time reversal,
• $G$ is a pseudodifferential imaging operator,
• $V_1$ projects onto the incoming-wave component. Crucially, this framework avoids explicit formation of the normal operator $F^* F$, eliminating low-frequency artifacts and yielding a partial inverse (a parametrix) which is microlocally invertible on the visible set (Root et al., 2010).

2. Microlocal Analysis, Imaging Condition, and Artifact Suppression

Key microlocal and support assumptions for RTM’s invertibility are:

• The background velocity $c(x)$ is smooth and density is constant.
• Multipathing of the source wave is excluded (Source-Multipath Exclusion, SME).
• Direct waves from source to receiver are excluded (Direct-Source Exclusion, DSE).
• No grazing rays (cutoff near tangential incidence).
• The reflectivity $r(x)$ is compactly supported below the acquisition surface.

The operator factorization in RTM introduces a one-way decoupling and imaging condition: $$G = K H\,:\, w(x,t) \mapsto A_s(x)^{-1}\, \partial_t^{-(n+1)/2} [\partial_t + c\,\nabla T_s \cdot \nabla_x] w(x, t)\,\Big|_{t=T_s(x)}$$ which acts on reverse-propagated fields and normalizes by the source amplitude $A_s(x)$. This structure mitigates the frequency blowup in the conventional normal operator and removes the signature "low-frequency artifacts" present in classical RTM (Root et al., 2010).

The imaging condition is further generalized to robustly handle velocity discontinuities, multipath scenarios, and varying acquisition geometries by invoking the cross-correlation of appropriately processed wavefields. In electromagnetic and acoustic scattering RTM, similar positivity and stability properties of the imaging functional are obtained by using the imaginary part of the cross-correlation, as justified via Helmholtz–Kirchhoff identities and point-spread function analysis (Chen et al., 2014, Chen et al., 2014).

3. Computational Approaches and Scalability

RTM is computationally dominated by repeated solutions of the wave equation over large spatial domains and time intervals. Implementations typically employ high-order finite-difference (FD), discontinuous Galerkin (DG), or Fourier-integral-operator (FIO) methods:

• High-order FD: RTM discretizes the Laplacian with 4th–8th order finite-difference stencils, using regular grids and explicit time-stepping for efficiency (Paul et al., 2016, Assis et al., 2019).
• Nodal DG: On unstructured meshes and for heterogeneous media, nodal DG with local time-stepping and thin communication halos provides high accuracy, GPU efficiency, and excellent scalability (Modave et al., 2015).
• FIO/wave-packet RTM: The dyadic parabolic (wave-packet) decomposition accelerates RTM by treating the reverse-time continuation as a sum of localized FIOs, preserving caustics and multiscale features without full time-domain wavefield storage (Andersson et al., 2013).

Hybrid distributed implementations (MPI+OpenMP or MPI+X paradigms) decompose the spatial domain, overlap communications with computation, and auto-tune scheduler parameters for optimal cache use and load balance. Coupled simulated annealing for dynamic chunk-size adjustment yields 10–30% runtime improvements by minimizing L3 cache misses and scheduling overhead (Assis et al., 2019, Paul et al., 2016).

Storage of the full forward wavefield for the imaging condition is a major bottleneck in RTM. Recent advances implement random boundary conditions (RBC) and source-reconstruction equations, which enable recovery of the full wavefield during backpropagation without the I/O penalty of disk storage—yielding $>300 \times$ reduction in storage and halving computational time while maintaining imaging quality (Barbosa et al., 2022).

4. RTM-Based Inverse Scattering, LSRTM, and Target-Oriented Imaging

In linearized inverse problems, RTM forms the adjoint. For quantitative, true-amplitude recovery, Least-Squares RTM (LSRTM) solves the normal equations: $$\min_{\delta m} \frac12 \| L[m_0] \delta m - \delta d \|_2^2$$ where $L[m_0]$ is the linearized Born modeling operator. The gradient and Hessian structures are: $$\text{Gradient:}\quad \nabla J(\delta m) = L^T (L \delta m - \delta d) \qquad \text{Hessian:}\quad H = L^T L$$ As shown in comparative studies, data-domain LSRTM with full Hessian utilization achieves superior multi-parameter recovery (e.g., simultaneous inversion of velocity and impedance) relative to image-domain deblurring approaches, which are cost-effective for single-parameter imaging but prone to parameter crosstalk and limited amplitude fidelity (Yang et al., 14 Aug 2025).

Target-oriented LSRTM uses surface and subsurface boundary representations via reciprocity or Marchenko redatuming to confine inversion to a region of interest. Marchenko-based double-focusing rigorously accounts for all overburden multiples, enabling virtual data at the target’s boundary and substantial reduction in computational cost—on the order of $10\times$ for fields restricted to small spatial and temporal windows (Shoja et al., 2023, Shoja et al., 2023).

Modern LSRTM incorporates sparsity regularization, stochastic minibatch optimization, and (when the source-time function is unknown) variable-projection estimation, enabling true-amplitude images with low computational overhead and robust recovery under noise and uncertainty (Yang et al., 2020).

5. Generalizations: Physics, Targets, and Modalities

RTM methodology is extended to diverse physical settings:

• Elastic, electromagnetic, and mixed-physics scattering: RTM for Maxwell’s equations, including vector-valued fields and both penetrable and impenetrable obstacles, employs imaging functionals built from the imaginary part of cross-correlations with respect to the vector Helmholtz kernel, with rigorous analysis ensuring strict positivity and stability properties (Chen et al., 2014, Cai et al., 2023).
• Obstacles, periodic media, and rough interfaces: Single-frequency and periodic RTM is optimally formulated via Green's function approaches and Helmholtz–Kirchhoff identities, ensuring maximal localization at object boundaries, with explicit point-spread function characterization (Cai et al., 2023, Li et al., 2022, Li et al., 2022).
• Locally rough and embedded geometries: Modified RTM imaging functionals, based on background Green's functions tailored to layered or rough backgrounds, enable simultaneous imaging of interfaces and embedded inclusions in unbounded or semi-infinite domains (Li et al., 2022).

6. Integration with Learning, Uncertainty Quantification, and Signal Processing

Machine learning approaches augment or surrogate the RTM pipeline:

• Encoder–Decoder Surrogates: Trainable neural networks can replace computationally intensive RTM runs for uncertainty quantification, accurately mapping from ensembles of uncertain velocity fields to image ensembles, achieving $90$–$95\%$ reduction in computational expense in many-query settings (Freitas et al., 2020).
• Feature-Based Matching: Small convolutional neural networks within LSRTM compute residuals in a learned feature space (SiameseLSRTM), improving noise suppression and amplitude consistency with marginal extra cost (Mu et al., 13 May 2025).
• Physics-informed Inpainting: RTM is re-interpreted as a reconstruction-inpainting operator for high-dimensional sensor arrays in audio event classification, illustrating transferable utility beyond seismic inversion (Tonami et al., 20 Jan 2026).

7. Practical Performance, Limitations, and Guidance

RTM offers unexcelled resolution and structural fidelity, but at high computational and storage cost. The method benefits substantially from:

• Careful numerical choices—high-order schemes for reduced dispersion, checkpoint/reconstruction schemes for minimal I/O, and high arithmetic-intensity codes exploiting modern hardware.
• Parametrization in the log-domain when performing multi-parameter inversion to balance amplitudes and simplify sensitivity (Yang et al., 14 Aug 2025).
• Avoidance of normal operator formation and artifact-prone approximations. Source normalization, proper exclusion of direct/grazing events, and one-way operator design are essential for artifact suppression (Root et al., 2010).
• Joint use of physical reciprocity/Marchenko methods and machine learning for bounding computational effort and expanding interpretability.

Limiting factors remain in full-volume 3D inversions, incomplete illumination, parameter trade-offs, and the handling of complex rheologies (e.g., elasticity, anisotropy) and noise models. Continued innovation in hybrid physics-informed learning, efficient solvers, and tailored imaging functionals is moving RTM towards broader applicability and interpretative robustness across geophysical and signal processing domains.