Modeling

Given a background slowness model $\textbf{s}_{o}$ , the Green's function for the Helmholtz equation is the solution of

$$[\bigtriangledown^{2}+\omega^{2} \textbf{s}_{o}(\textbf{x})^{2}] G_{o}(\textbf{x}\vert\textbf{x}_{s})= -\delta(\textbf{x}-\textbf{x}_{s}),$$ (21)

where $G_{o}(\textbf{x}\vert\textbf{x}_{s})$ is the Green's function associated with the background slowness $\textbf{s}_{o}$ and an impulsive point source at $\textbf{x}_{s}$ , $\textbf{x}$ is the listening location, and $\omega$ is the angular frequency. For a point source at $\textbf{x}_{s}$ with spectrum $W(\omega)$ (corresponding to the source term $F(\textbf{x},\omega)=-\delta(\textbf{x}-\textbf{x}_{s})W(\omega)$ ), the solution is $\textbf{P}_{o}(\textbf{x}\vert\textbf{x}_{s})=W(\omega)G_{o}(\textbf{x}\vert\textbf{x}_{s})$ .

If the background slowness perturbation is represented by the amount of $\delta \textbf{s}(\textbf{x})$ , the true slowness model is described by $\textbf{s}(\textbf{x})=\textbf{s}_{o}+\delta \textbf{s}(\textbf{x})$ . The full wavefield is obtained by solving the Helmholtz equation with slowness model $\textbf{s}(\textbf{x})$ ,

$$[\bigtriangledown^{2}+\omega^{2} \textbf{s}(\textbf{x})^{2}]\textbf{P}=F,$$ (22)

with the same source term $F=-\delta(\textbf{x}-\textbf{x}_{s})W(\omega)$ as before. Our goal is to calculate the scattered wavefield $\textbf{P}_{1}=\textbf{P}-\textbf{P}_{o}$ induced by the slowness perturbation $\delta \textbf{s}(\textbf{x})$ with a linear modeling operator. Plugging $\textbf{s}(\textbf{x})=\textbf{s}_{o}+\delta \textbf{s}(\textbf{x})$ into equation (2.2), I get

$$[\bigtriangledown^{2}+\omega^{2} \textbf{s}_{o}(\textbf{x})^{2}+2\omega^{2} \textbf{s}_{o}(\textbf{x})\delta \textbf{s}(\textbf{x})] \textbf{P}= F,$$ (23)

where the high order term $O(\delta \textbf{s}^{2})$ is neglected. According to Green's theorem, moving the 3rd term in equation (2.3) to the right side, multiplying both sides with the Green's function $G_{o}(\textbf{x}\vert\textbf{x}')$ and integrating over the whole volume with index $\textbf{x}'$ , gives the Lippmann-Schwinger equation,

$$\textbf{P}(\textbf{x}) = \int G_{o}(\textbf{x}\vert\textbf{x}')F(\textbf{x}')d\textbf{x}' ... ...P}(\textbf{x}'\vert\textbf{x}_{s})G_{o}(\textbf{x}\vert\textbf{x}')d\textbf{x}'$$

$$= \textbf{P}_{o}(\textbf{x})+\omega^{2} \int \textbf{m}(\textbf{x}'... ...}(\textbf{x}'\vert\textbf{x}_{s})G_{o}(\textbf{x}\vert\textbf{x}')d\textbf{x}',$$ (24)

which is an integral equation with unknown $\textbf{P}(\textbf{x})$ on both sides and $\textbf{m}(\textbf{x}')=-2\textbf{s}(\textbf{x}')\delta \textbf{s}(\textbf{x}')$ representing the reflectivity model. Here, $\textbf{P}_{o}(\textbf{x})=\int \textbf{G}_{o}(\textbf{x}\vert\textbf{x}')F(\textbf{x}')d\textbf{x}'$ is the pressure field associated with the background velocity model. Defining $\textbf{P}(\textbf{x}'\vert\textbf{x}_{s})=W(\omega) G(\textbf{x}'\vert\textbf{x}_{s})$ and applying the Born approximation $G(\textbf{x}'\vert\textbf{x}_{s}) \approx G_{o}(\textbf{x}'\vert\textbf{x}_{s})$ to the right side and assuming that $\delta \textbf{s}(\textbf{x})$ is small gives the scattered field under Born approximation

$$\textbf{P}_{1} = \textbf{P}(\textbf{x})-\textbf{P}_{o}(\textbf{x})$$ (25)

$$\approx \omega^{2} \int W(\omega) \textbf{m}(\textbf{x}') G_{o}(\textbf{x}'\vert\textbf{x}_{s}) G_{o}(\textbf{x}\vert\textbf{x}')d\textbf{x}',$$ (26)

where formula (2.5) represents a non-linear equation for calculation of the scattered wavefield and equation (2.6) is a linear equation (Mulder and Plessix, 2004a). With $\textbf{P}_{o}(\textbf{x}')=W(\omega)G_{o}(\textbf{x}'\vert\textbf{x}_{s})$ , the linear modeling requires the solutions of

$$[\bigtriangledown^{2}+\omega^{2}\textbf{s}_{o}(\textbf{x})^{2}]\textbf{P}_{o}=F,$$

$$[\bigtriangledown^{2}+\omega^{2}\textbf{s}_{o}(\textbf{x})^{2}]\textbf{P}_{1}= \omega^2 \textbf{m}(\textbf{x}')\textbf{P}_{o}(\textbf{x}').$$ (27)

These fields can be computed by two finite-difference simulations: one with the original point source $F$ and background slowness model $\textbf{s}_{o}$ to generate $\textbf{P}_{o}$ ; The second finite-difference simulation also uses the background slowness model $\textbf{s}_{o}$ , but the source term is $\omega^{2}\textbf{m}(\textbf{x}')\textbf{P}_{o}(\textbf{x}')$ , where $\omega^{2}$ becomes the 2nd-order time derivative in the time domain. The adjoint of the linear operator is the reverse time migration (RTM) operator, so the RTM equation is

$$\textbf{m}_{mig}(\textbf{x})=\sum_{\textbf{x}_{s}} \int \omega^{2... ...extbf{x}\vert\textbf{x}_{s}) G^{*}_{o}(\textbf{x}\vert\textbf{x}')d\textbf{x}'.$$ (28)

In the following section, matrix-vector notation will be used to represent the operators, such that the non-linear modeling operator is defined as $\textbf{A}: \textbf{d}=\textbf{A}(\textbf{m})=\textbf{P}-\textbf{P}_{o}$ , the linear modeling operator is $\textbf{L}: \textbf{d}=\textbf{Lm}$ , and the reverse time migration operator is $\textbf{L}^{T}: \textbf{m}=\textbf{L}^{T}\textbf{d}$ . The modeling operator $\textbf{L}$ is called the Fréchet derivative of $\textbf{A}$ and $\textbf{L}^{T}$ is the adjoint of $\textbf{L}$ .