Theory

The theory of least-squares reverse time migration is well established (Symes and Carazzone, 1991; Mulder and Plessix, 2004; Dai et al., 2012). In this section, I will first review the theory of LSRTM assuming the constant density acoustic wave equation,

$$\frac{1}{c(\textbf{x})^2}\frac{\partial^2 p(\textbf{x},t;\textbf{... ...partial t^2}-\bigtriangledown^2 p(\textbf{x},t;\textbf{x}_s)=s(t;\textbf{x}_s),$$ (14)

where $c(\textbf{x})$ is the velocity distribution, and $p(\textbf{x},t;\textbf{x}_s)$ is the pressure field associated with the source term $s(t;\textbf{x}_s)$ . A perturbation in the velocity model $c(\textbf{x})\rightarrow c(\textbf{x})+\delta c(\textbf{x})$ will generate a wavefield $p(\textbf{x},t;\textbf{x}_s)\rightarrow p(\textbf{x},t;\textbf{x}_s)+\delta p(\textbf{x},t;\textbf{x}_s)$ , which obeys the equation

$$\frac{1}{(c(\textbf{x})+\delta c(\textbf{x}))^2}\frac{\partial^2 ... ...xtbf{x},t;\textbf{x}_s)+\delta p(\textbf{x},t;\textbf{x}_s))=s(t;\textbf{x}_s).$$ (15)

Expanding the velocity term according to

$$\frac{1}{(c(\textbf{x})+\delta c(\textbf{x}))^2} \approx \frac{1}{c(\textbf{x})^2}-\frac{2\delta c(\textbf{x})}{c(\textbf{x})^3},$$ (16)

and subtracting equation from equation yields the wave equation for the wavefield perturbation $\delta p(\textbf{x},t;\textbf{x}_s)$

$$\frac{1}{c(\textbf{x})^2}\frac{\partial^2 \delta p(\textbf{x},t;\... ...al t^2}\frac{2\delta c(\textbf{x})}{c(\textbf{x})^3}+O(\delta c(\textbf{x})^2).$$ (17)

Neglecting the higher order terms and defining the reflectivity model as $m(\textbf{x})=\frac{2\delta c(\textbf{x})}{c(\textbf{x})}$ , the above equation becomes

$$\frac{1}{c(\textbf{x})^2}\frac{\partial^2 \delta p(\textbf{x},t;\... ...},t;\textbf{x}_s)=m(\textbf{x})\bigtriangledown^2 p(\textbf{x},t;\textbf{x}_s).$$ (18)

Equations and will be used to derive the Born modeling operator. Numerically, the calculation of the reflection data $\delta p(\textbf{x},t;\textbf{x}_s)$ requires two finite-difference simulations: one to solve equation to obtain the wavefield $p(\textbf{x},t;\textbf{x}_s)$ , and one to solve equation for the reflection data $\delta p(\textbf{x},t;\textbf{x}_s)$ . The wavefield $\delta p(\textbf{x},t;\textbf{x}_s)$ will be recorded at the receiver position $\textbf{x}_g$ to give the shot gather $d(\textbf{x}_g,t;\textbf{x}_s)$ . By the adjoint state method (Plessix, 2006), the migration operation of a shot gather $d(\textbf{x}_g,t;\textbf{x}_s)$ requires two finite-difference simulations, one for the source-side wavefield and one for the receiver-side wavefield:

$$\frac{1}{c(\textbf{x})^2}\frac{\partial^2 p(\textbf{x},t;\textbf{... ...partial t^2}-\bigtriangledown^2 p(\textbf{x},t;\textbf{x}_s)=s(t;\textbf{x}_s),$$ (19)

$$\frac{1}{c(\textbf{x})^2}\frac{\partial^2 q(\textbf{x},t;\textbf{... ...\bigtriangledown^2 q(\textbf{x},t;\textbf{x}_s)=d(\textbf{x}_g,t;\textbf{x}_s),$$ (20)

where $q(\textbf{x},t;\textbf{x}_s)$ is the receiver-side wavefield. Note that the source-side wavefield $p(\textbf{x},t;\textbf{x}_s)$ propagates forward in time but the receiver-side wavefield $q(\textbf{x},t;\textbf{x}_s)$ propagates backward in time. The migration image associated with the shot at $\textbf{x}_s$ is produced by applying the imaging condition

$$m(\textbf{x};\textbf{x}_s)=\sum_t \bigtriangledown^2 p(\textbf{x},t;\textbf{x}_s) \cdot q(\textbf{x},t;\textbf{x}_s).$$ (21)

To simplify the formulas, matrix-vector notation will be used to represent the Born modeling operator

$$\textbf{d}_i=\textbf{L}_i\textbf{m},$$ (22)

where $\textbf{d}_i$ is the reflection data vector for the $ith$ shot, $\textbf{m}$ is a reflectivity model, and $\textbf{L}_i$ represents the Born modeling operator associated with the $ith$ shot. Similarly, the reverse time migration operator can be expressed as

$$\textbf{m}_{mig,i}=\textbf{L}^{T}_{i}\textbf{d}_{i}.$$ (23)

with $\textbf{m}_{mig,i}$ indicating the migration image for the $ith$ shot and $\textbf{L}^{T}_{i}$ representing the migration operator associated with the $ith$ shot.