Wave-equation tomography using image-space phase-encoded data

image-space wave-equation tomography

Image-space wave-equation tomography is a non-linear inverse problem that tries to find an optimal background slowness that minimizes the residual field, $\Delta {\bf I}$, defined in the image space. The residual field is derived from the background image, $\bf I$, which is computed with a background slowness. The general form of the residual field is (Biondi, 2008)

$${\Delta {\bf I}} = {\bf I} - {\bf F}({\bf I}),$$ (1)

where $\bf F$ is a focusing operator, which measures the focusing of the migrated image. In particular, in the Differential Semblance Optimization (DSO) method (Shen, 2004), the focusing operator takes the form:

$${\bf F}({\bf I}) = \left( {\bf 1} - {\bf O} \right) {\bf I},$$ (2)

where ${\bf 1}$ is the identity operator and ${\bf O}$ is the DSO operator either in the subsurface offset domain or in the angle domain (Shen, 2004).

Under $\ell_2$ norm, the tomography objective function can be written as follows:

$$J = \frac{1}{2}\vert\vert{\Delta {\bf I}}\vert\vert _2 = \frac{1}{2}\vert\vert{\bf I}-{\bf F}({\bf I})\vert\vert^2.$$ (3)

The gradient of $J$ with respect to the slowness ${\bf s}$ is

$${\bf\nabla} J = \left(\frac{\partial {\bf I}}{\partial {\bf s}} -... ...f I})}{\partial {\bf s}} \right)^{\ast} \left( {\bf I}-{\bf F}({\bf I})\right),$$ (4)

where $^{\ast}$ denotes the adjoint.

The linear operator $\left.\frac{\partial{\bf I}}{\partial{\bf s}}\right\vert _{{\bf s}=\widehat{\bf s}}$, which defines a linear mapping from the slowness perturbation ${\Delta {\bf s}}$ to the image perturbation ${\Delta {\bf I}}$, can be computed by expanding the image ${\bf I}$ around the background slowness ${\widehat {\bf s}}$. Keeping only the zeroth and first order terms, we get the linear operator $\left.\frac{\partial{\bf I}}{\partial{\bf s}}\right\vert _{{\bf s}=\widehat{\bf s}}$ as follows:

$${\Delta {\bf I}} = \left. \frac{\partial{\bf I}}{\partial{\bf s}}\right\vert _{{\bf s}=\widehat{\bf s}} {\Delta {\bf s}} = {\bf T}{\Delta {\bf s}},$$ (5)

where ${\Delta {\bf I}} = {\bf I} - {\widehat {\bf I}}$, ${\widehat {\bf I}}$ is the background image computed with the background slowness ${\widehat {\bf s}}$ and ${\Delta {\bf s}} = {\bf s} - {\widehat {\bf s}}$. ${\bf T} = \left. \frac{\partial{\bf I}}{\partial{\bf s}}\right\vert _{{\bf s}=\widehat{\bf s}}$ is the wave-equation tomographic operator. The tomographic operator can be evaluated either in the source and receiver domain (Sava, 2004) or in the shot-profile domain (Shen, 2004).

In the shot-profile domain, both source and receiver wavefields are downward continued with the one-way wave equations (Claerbout, 1971)

$$\left\{ \begin{array}{l} \left( \frac{\partial}{\partial z}+i\Lam... ...{\bf x}_s,\omega) = {f_s(\omega)\delta({\bf x}-{\bf x}_s)} \end{array} \right.,$$ (6)

and

$$\left\{ \begin{array}{l} \left( \frac{\partial}{\partial z}+i\Lam... ... U(x,y,z=0,{\bf x}_s,\omega) = Q(x,y,z=0,{\bf x}_s,\omega) \end{array} \right.,$$ (7)

where $D({\bf x},{\bf x}_s,\omega)$ is the source wavefield for a single frequency $\omega$ at image point ${\bf x}=(x,y,z)$ with the source located at ${\bf x}_s=(x_s,y_s,0)$; $U({\bf x},{\bf x}_s,\omega)$ is the receiver wavefield for a single frequency $\omega$ at image point ${\bf x}$ for the source located at ${\bf x}_s$; $f_s(\omega)$ is the frequency dependent source signature, and ${f_s(\omega)\delta({\bf x}-{\bf x}_s)}$ defines the point source function at ${\bf x}_s$, which serves as the boundary condition of Equation 6. $Q(x,y,z=0,{\bf x}_s,\omega)$ is the recorded shot gather for the shot located at ${\bf x}_s$, which serves as the boundary condition of Equation 7. $\Lambda$ is the square-root operator

$$\begin{array}{l} \Lambda = \sqrt{\omega ^2 s^2({\bf x})-\vert{\bf k}\vert ^2}, \end{array}$$ (8)

where $s({\bf x})$ is the slowness at ${\bf x}$; ${\bf k}=(k_x,k_y)$ is the spatial wavenumber vector. The image is computed by applying the cross-correlation imaging condition:

$$I({\bf x},{\bf h}) = \sum_{{\bf x}_s}\sum_{\omega} \overline{D({\bf x}-{\bf h},{\bf x}_s,\omega)} U({\bf x}+{\bf h},{\bf x}_s,\omega),$$ (9)

where the overline stands for complex conjugate; $D({\bf x},{\bf x}_s,\omega)$ is the source wavefield for a single frequency $\omega$ at image point ${\bf x}=(x,y,z)$ with the source located at ${\bf x}_s=(x_s,y_s,0)$; $U({\bf x},{\bf x}_s,\omega)$ is the receiver wavefield and ${\bf h}=(h_x,h_y,h_z)$ is the subsurface half-offset.

The perturbed image can be derived by the application of the chain rule to Equation 9:

$$\Delta I({\bf x},{\bf h}) = \sum_{{\bf x}_s}\sum_{\omega} \left( \overline{\Delta D({\bf x}-{\bf h},{\bf x}_s,\omega)} {\widehat U}({\bf x}+{\bf h},{\bf x}_s,\omega) + \right.$$

$$\left. \overline{{\widehat D}({\bf x}-{\bf h},{\bf x}_s,\omega)} \Delta U ({\bf x}+{\bf h},{\bf x}_s,\omega) \right),$$ (10)

where ${\widehat D}({\bf x}-{\bf h},{\bf x}_s,\omega)$ and ${\widehat U}({\bf x}+{\bf h},{\bf x}_s,\omega)$ are the background source and receiver wavefields computed with the background slowness ${\widehat s}({\bf x})$; $\Delta D({\bf x}-{\bf h},{\bf x}_s,\omega)$ and $\Delta U({\bf x}+{\bf h},{\bf x}_s,\omega)$ are the perturbed source wavefield and perturbed receiver wavefield, which are the results of the slowness perturbation $\Delta s({\bf x})$.

To evaluate the adjoint of the tomographic operator, ${\bf T}^{\ast }$, we first apply the adjoint of the imaging condition to get the perturbed source and receiver wavefields, $\Delta D ({\bf x}+{\bf h},{\bf x}_s,\omega)$ and $\Delta U({\bf x}+{\bf h},{\bf x}_s,\omega)$, as follows

$$\Delta D({\bf x},{\bf x}_s,\omega) = \sum_{\bf h} \Delta I({\bf x},{\bf h}) {\widehat U}({\bf x}+{\bf h},{\bf x}_s,\omega)$$

$$\Delta U({\bf x},{\bf x}_s,\omega) = \sum_{\bf h} \Delta I({\bf x},{\bf h}) {\widehat D}({\bf x}-{\bf h},{\bf x}_s,\omega).$$ (11)

The perturbed source and receiver wavefields satisfy the following one-way wave equations, linearized with respect to slowness:

$$\left\{ \begin{array}{l} \left( \frac{\partial}{\partial z}+i\Lam... ...elta s({\bf x})\\ \Delta D(x,y,z=0,{\bf x}_s,\omega) = 0 \end{array} \right. ,$$ (12)

and

$$\left\{ \begin{array}{l} \left( \frac{\partial}{\partial z}+i\Lam... ...Delta s({\bf x})\\ \Delta U(x,y,z=0,{\bf x}_s,\omega) = 0 \end{array} \right..$$ (13)

When solving the optimization problem, the gradient of the objective function is obtained by computing the perturbed wavefields using the adjoint of the imaging operator (equation 11), where the image perturbation results from the application of a focusing operator (equation 1) on the background image; then the scattered wavefields are obtained by applying the adjoint of the one-way wave equations 12 and 13; and, finally, the adjoint scattering operator cross-correlates the upward propagated the scattered wavefields with the modified background wavefields (term in the parenthesis on the right-hand side of equations 12 and 13). Figure 1 displays the image-space wave-equation tomography flowchart. The gray box on the left represents the process of obtaining the image perturbation, while the gray box on the right corresponds to the application of the adjoint of the wave-equation tomography operator. WE stands for wavefield extrapolation. The light gray boxes contain the wavefields, images and slowness perturbation. The processes and operators are represented as white boxes. More detailed information on how to evaluate the forward and adjoint operators can be found in Tang et al. (2008).

ISWET
Figure 1. Image-space wave-equation tomography flowchart. The gray box on the left represents the process of obtaining the image perturbation, while the gray box on the right corresponds to the application of the adjoint of the wave-equation tomography operator.[NR]

Wave-equation tomography using image-space phase-encoded data