First-order Born scattering operators

Imaging by wavefield extrapolation (WE) is based on recursive continuation of wavefields $\UU$ from a given depth level to the next by means of an extrapolation operator $\Eop$.At every extrapolation step, we can write that  

$$\UU\atzz= \Eop_z\lb \UU\atzo \rb,$$ (92)

where $\UU\atzo$ is the wavefield at the top of the slab, and $\UU\atzz$ is the wavefield at the bottom of the slab. The operator $\Eop$ involves a spatially-dependent phase shift described by:

$$\Eop_z\lb \rb = e^{i \kzz \dz},$$ (93)

where $\kzz$ represents the depth wavenumber, and $\dz$ the wavefield extrapolation depth step. The relation () corresponds to the analytical solution of the differential equation  

$$\UU'(z) = i\kzz \UU(z)$$ (94)

which describes depth extrapolation of monochromatic plane waves (24). The ' sign represents a derivative with respect to the depth z. The depth wavenumber $\kzz$ is given by the one-way wave equation, also known as the single square root (SSR) equation  

$$\kzz = \sqrt{\omega^2 {s}^2 - \vert\bf k\vert^2},$$ (95)

where $\ww$ is the temporal frequency, s is the laterally variable slowness of the medium, and $\bf k$ is the horizontal wavenumber. I use the laterally variable s and the horizontal wavenumber $\bf k$ in SSR just for conciseness, although such a notation not mathematically correct in laterally varying media.

Since downward continuation by Fourier-domain phase shift can be applied for slowness models that only vary with depth, we need to split the operator $\Eop$ into two parts: a constant slowness continuation operator applied in the $\ww-{\bf k}$ domain, which accounts for the propagation in depth, and a screen operator applied in the $\ww-{\bf x}$ domain, which accounts for the wavefield perturbations due to the lateral slowness variations. In essence, we approximate the vertical wavenumber $\kzz$ with its constant slowness counterpart ${\kzz}_0$, corrected by a term describing the spatial variability of the slowness function (81).

Furthermore, we can separate the depth wavenumber $\kzz$ into two components, one which corresponds to the background medium $\widetilde{\kzz}$ and one which corresponds to a perturbation of the medium:

$$\kzz = \widetilde{\kzz}+ \DEL \kzz\;.$$ (96)

In a first-order approximation, we can relate those two depth wavenumbers by a Taylor series expansion:

&& + . d d s |_s=ss - s
&& + s^2 s^2 - |k|^2 s - s,

where $s \lp\zz,\xx\rp$ is the slowness corresponding to the perturbed medium, and $\tilde{s}\lp\zz,\xx\rp$ is the background slowness.

Within any depth slab, we can extrapolate the wavefield from the top either in the perturbed or in the background medium. The wavefields at the bottom of the slab, ${\widetilde{\UU}}\atzz=\UU\atzo e^{i \widetilde{\kzz}\dz}$ and $\UU\atzz=\UU\atzo e^{i \kzz \dz}$ are related by the relation  

$$\UU\atzz\approx {\widetilde{\UU}}\atzz e^{i\DEL \kzz\dz} \;.$$ (97)

rytov.w is a direct statement of the Rytov approximation (56), since the wavefields at the bottom of the slab correspond to different phase shifts related by a linear equation.

The wavefield perturbation $\DEL \VV$ at the bottom of the slab is obtained by subtracting the background wavefield $\widetilde{\UU}$ from the perturbed wavefield $\UU$:

&& -
&& e^i-1
&& e^i e^i . d d s |_s=ss -1 ,

where $\Delta s=s-\tilde{s}$ is the perturbation between the correct and the background slownesses at depth z.

In operator form we can write

$${\DEL \VV}\atzz= \Eop_z\lb \Nop_z\lp{\widetilde{\UU}}\atzo\rp \lb{\Delta s}\atzo\rb \rb,$$ (98)

where $\Eop_z$ represents the downward continuation operator at depth z, and $\Nop_z$ represents the Rytov scattering operator which is dependent on the background wavefield ${\widetilde{\UU}}\atzo$ and the slowness perturbation ${\Delta s}\atzo$ at that depth level:  

$$\Nop_z\lp{\widetilde{\UU}}\atzo\rp \lb{\Delta s}\atzo\rb = \... ...=\tilde{s}}{\Delta s}\atzo\dz} -1 \rp {\widetilde{\UU}}\atzo\;.$$ (99)

In this approximation, we assume that the scattered wavefield is generated only by the background wavefield and we ignore all multi-scattering effects. For the Born approximation (56), we further assume that the wavefield differences are small, such we can linearize the exponential according to the relation $e^{i\Delta \phi} \approx 1+i\Delta \phi$. With this new approximation, the expression for the downward-continued scattered wavefield becomes:  

$${\DEL \VV}\atzz\approx e^{i \widetilde{\kzz}\dz} \lp i \left... ...rt _{s=\tilde{s}}{\Delta s}\atzo\dz \rp {\widetilde{\UU}}\atzo.$$ (100)

In operator form, we can write the scattered wavefield at z as

$${\DEL \VV}\atzz= \Eop_z\lb \Sop_z\lp{\widetilde{\UU}}\atzo\rp \lb{\Delta s}\atzo\rb \rb$$ (101)

where $\Eop_z$ represents the downward continuation operator at depth z, and $\Sop_z$ represents the Born scattering operator which is dependent on the background wavefield and operates on the slowness perturbation at that depth level.

The linear scattering operator $\Sop$ is a mixed-domain operator similar to the extrapolation operator $\Eop$. This operator depends on the background wavefield and background slowness by the expression:  

$$\Sop_z\lp{\widetilde{\UU}}\atzo\rp \lb{\Delta s}\atzo\rb \ap... ...vert _{s=\tilde{s}}\dz {\Delta s}\atzo{\widetilde{\UU}}\atzo\;.$$ (102)

In practice, we can implement the scattering operator described by born.op in different ways.

One option is to implement the Born operator born.op in the space domain using an expansion (50) like  

$$\left. \frac{d \kzz} {d s } \right\vert _{s=\tilde{s}}\approx \ww \lp 1+ \SQR4exp{\S}+\dots \rp \;.$$ (103)

In practice, the summation of the terms in born.exp involves forward and inverse Fast Fourier Transforms (FFT and IFT) and multiplication in the space domain with the spatially variable $\tilde{s}$:

$${\DEL \VV}\atzo= i \ww \dz {\Delta s}\atzo\lb 1+\sum_{j=1,..... ...rt^{2j} \texttt{FFT} \lb {\widetilde{\UU}}\atzo \rb \rb \rb\;,$$ (104)

where $c_j=\frac{1}{2},\frac{3}{8},\dots$Another option is to implement the Born operator born.op in the Fourier domain relative to the constant reference slowness in any individual slab. In this case, we can write  

$$\left. \frac{d \kzz} {d s } \right\vert _{s=s_o}\approx \ww... ... s_o}{\sqrt{\omega^2 {s_o}^2-(1-i\eta)^2\vert\bf k\vert^2}} \;,$$ (105)

where $\eta$ as a damping parameter which avoids division by zero (48). In practice, the implementation of born.sqr involves forward and inverse Fast Fourier Transforms (FFT and IFT):

$${\DEL \VV}\atzo= i \dz \; \texttt{IFT} \lb \left. \frac{d \k... ...xttt{FFT} \lb {\widetilde{\UU}}\atzo{\Delta s}\atzo \rb \rb \;.$$ (106)