Gradient of image-space wave-equation tomography by the adjoint-state method

Gradient of ISWET by the adjoint-state method

Image-space wave-equation tomography aims to iteratively solve for the slowness model, $s=s({\bf x})$, that minimizes the nonlinear objective function:

$$J(s)=\frac{1}{2} \left \vert \left \vert \Delta r(s) \right \vert... ...1}{2} \left \vert \left \vert r(s) - {\bf M} r(s) \right \vert \right \vert ^2,$$ (7)

where $\Delta r = \Delta r({\bf x},{\bf h})$ is the image perturbation that measures the accuracy of the slowness model. To compute $\Delta r$, a differential residual-focusing operator $\bf M$ is applied to the image $r = r({\bf x},{\bf h})$ obtained with the current slowness (Biondi, 2008), using either differential residual prestack migration (Sava and Biondi, 2004a,b) or differential-semblance optimization (DSO) operators (Shen and Symes, 2008). In this paper, operators are represented by bold capital letters.

If the differential residual-focusing operator $\bf M$ is independent of the slowness, the gradient of this objective function evaluated at the current slowness ${\widehat s}={\widehat s}(\bf x)$ is

$$\nabla J(s)= \left. \left ( \frac{\partial r}{\partial s} \right )'\right \vert _{s=\widehat s} \left ({\bf I - M}' \right ) {\Delta \widehat r},$$ (8)

where the prime denotes the adjoint, $\bf I$ is the identity operator, and $\Delta \widehat r=\Delta \widehat r({\bf x},{\bf h})$ is the perturbed image obtained with the current slowness model. The linear operator $\frac{\partial r}{\partial s}$ defines the mapping $\frac{\partial r}{\partial s} \Delta s = \Delta r$ between the slowness perturbation $\Delta s$ and the image perturbation $\Delta r$; it is called the image-space wave-equation tomographic operator.

Because the image-space wave-equation tomographic operator is composed of different operators, it is difficult to envision from equation 8 which operations are performed to compute the gradient. Therefore, for a clear explanation of the operators involved, I use the adjoint-state method to derive the gradient of the objective function (equation 7).

In migration with generalized sources or shot-profile migration, the source and receiver wavefields are propagated independently, and the image $r_z = r_z({\bf x},{\bf h})$ at a depth level $z$, is computed by the crosscorrelation

$$r_z({\bf x},{\bf h}) = \sum_{\omega} d_z^*({\bf x}-{\bf h},\omega) u_z({\bf x}+{\bf h},\omega),$$ (9)

where $d_z({\bf x},\omega)$ is the source wavefield for a single frequency $\omega$ at horizontal coordinates ${\bf x}=(x,y)$ ; $u_z({\bf x},\omega)$ is the receiver wavefield, ${\bf h}=(h_x,h_y)$ is the subsurface half-offset, and ``$^*$'' stands for the complex-conjugate. An additional summation over shots is required when migrating more than one shot. Hereafter, letters $d$ and $u$ stand for source and receiver wavefields, respectively, irrespective of the migration scheme.

In a more compact notation, not explicitly writing the dependencies on $\bf x$ and $\bf h$, equation 9 can be re-written as follows:

$$r_z = {\bf {S}}{\bf {D}}_z'(\omega) u_z(\omega) = {\bf {S}}{\bf {U}}_z(\omega) d_z^*(\omega),$$ (10)

where ${\bf {D}}$ and ${\bf {U}}$ are convolutional operators composed of $(h_x,h_y)$-shifted versions of $d_z({\bf x},\omega)$ and $u_z({\bf x},\omega)$, respectively. Operator ${\bf {S}}$ corresponds to the summation over frequency.

For subsequent depth levels, $d({\bf x},\omega)$ is computed by means of the recursive downward propagation

$$\left\{ \begin{array}{l} d_{z+1}(\omega) = {\bf T}_z^{\downarrow}(\omega,s)d_z(\omega) \\ d_1(\omega) = q(\omega), \end{array}\right.$$ (11)

where ${\bf T}_z^{\downarrow}$ is the downward continuation operator, which is a function of the slowness $s$, and $q(\omega)$ is the source wavefield used as a boundary condition. In the case of conventional shot-profile migration, $q(\omega) = f_s(\omega)\delta({\bf x - x_s})$ is the source signature located at ${\bf x}_s=(x_s,y_s,0)$. If using the generalized sources in the image space, $q(\omega)$ represents the image-space phase-encoded source wavefield of equation 3.

The downward continuation of the receiver wavefield is performed by

$$\left\{ \begin{array}{l} u_{z+1}(\omega) = {\bf T}_z^{\downarrow}(\omega,s)u_z(\omega) \\ u_1(\omega) = w(\omega), \end{array}\right.$$ (12)

where $w(\omega)$ is the recorded data at the surface for shot-profile migration. If using generalized sources in the image space, $w(\omega)$ is the phase-encoded areal receiver wavefield of equation 4. In equations 11 and 12, I omitted the dependencies of the wavefield with respect to ${\bf x}$. The subscript ${\it 1}$ in equations 11 and 12 represents the surface for the shot-profile migration and the ``collection'' depth level, $z_c$, for the image-space phase-encoded wavefields.

In the image-space wave-equation tomography problem, the perturbed source and receiver wavefields, and the image perturbation are used to compute the slowness perturbation that updates the current slowness model. From the perturbation theory, we have $d = \widehat d + \Delta d$, $u = \widehat u + \Delta u$, and, consequently, $r = \widehat r + \Delta r$ are physical realizations with $s = \widehat s + \Delta s$, where the $hat$ refers to fields obtained with the background slowness. To the first order (Born approximation), these perturbed fields are given by

$$\Delta d_{z+1}(\omega) = {\bf T}_z^{\downarrow}(\omega,\widehat s... ... s_z %\Delta {\bf T}_z^{\downarrow}(\omega,\widehat s) {\widehat d}_z(\omega),$$ (13)

and

$$\Delta u_{z+1}(\omega) = {\bf T}_z^{\downarrow}(\omega,\widehat s... ... s_z. %\Delta {\bf T}_z^{\downarrow}(\omega,\widehat s) {\widehat u}_z(\omega)$$ (14)

The diagonal operators $\widetilde {\widetilde{\bf D}}_z$ and $\widetilde {\widetilde{\bf U}}_z$ have in the diagonal entries the scattered source and receiver wavefields, respectively. These wavefields are given by the action of the scattering operator $\Delta {\bf T}_z^{\downarrow}$ on the background wavefields:

$$\widetilde {\widetilde {\bf D}}_z(\omega) = \Delta {\bf T}_z^{\do... ...2 \widehat s^2 - \left \vert {\bf k} \right \vert ^2}} dz~ \widehat d_z(\omega)$$ (15)

and

$$\widetilde {\widetilde {\bf U}}_z(\omega) = \Delta {\bf T}_z^{\do... ... \widehat s^2 - \left \vert {\bf k} \right \vert ^2}} dz~ \widehat u_z(\omega).$$ (16)

The perturbed image is given by

$$\Delta r_{z} = {\bf S} \left ({\widehat {\bf U}}_z(\omega) \Delta d_z^*(\omega) + {\widehat {\bf D}}_z'(\omega) \Delta u_z(\omega) \right ).$$ (17)

The matrix representations of equations 13, 14, 17 are

$$\Delta {\underline {\bf d}} = {\bf T}^{\downarrow} \Delta {\underline {\bf d}} + {\widetilde {\widetilde {\bf P}}} {\bf S}' \Delta {\bf s},$$ (18)

$$\Delta {\underline {\bf u}} = {\bf T}^{\downarrow} \Delta {\underline {\bf u}} + {\widetilde {\widetilde{\bf U}}} {\bf S}' \Delta {\bf s},$$ (19)

and

$$\Delta {\bf {\underline r}} = {\bf S} \left ({\widehat {\bf U}} \... ...nderline {\bf d}}^* + {\widehat {\bf D}}' \Delta {\underline {\bf u}} \right ),$$ (20)

where ${\bf S}'$ is a spreading operator that replicates the slowness perturbation for every frequency.

Equations 18, 19, and 20 are the forward modeling equations of the image-space wave-equation tomography problem using the generalized sources or shot-profile schemes. They depend on the state variables $\Delta {\underline {\bf d}}$, $\Delta {\underline {\bf u}}$, and $\Delta {\underline {\bf r}}$. Plessix (2006) describes how to compute the adjoint states using the augmented functional methodology. By introducing the adjoint-state variables ${\underline \lambda}_d$, ${\underline \lambda}_u$, and ${\underline \lambda}_r$, the augmented Lagrangian reads

$$\begin{array}{l} \mathcal{L}(\Delta {\underline {\bf d}}, \De... ...a {\underline {\bf u}} \right ) \right \rangle. \\ \end{array}$$

The adjoint-state variables are computed by taking the derivative of $\mathcal{L}$ with respect to the state variables and equating to zero, which gives

Notice that

$$\left ( {\bf I} - {\bf T}^\downarrow \right )' = \left ( {\bf I} - {\bf T}^{\downarrow '} \right ) = \left ( {\bf I} - {\bf T}^{\uparrow} \right )$$ (22)

corresponds to the recursive upward propagation operator. Therefore, equations 21a and 21b can be written as

which correspond to the recursive upward propagation of the perturbed wavefields resulting from the convolution of the wavefields computed with the current slowness and the perturbed image.

Finally, the gradient of $J$ is

$$\nabla_s J({\bf s}) = {\bf S} \left ({\widetilde {\widetilde{\bf ... ... \lambda}_d + {\widetilde {\widetilde{\bf U}}}'{\underline \lambda}_u \right ).$$ (24)

To compute the gradient, the adjoint-state wavefields, ${\underline \lambda}_d$ and ${\underline \lambda}_u$, are upward propagated and cross-correlated in time with the scattered wavefields.

Gradient of image-space wave-equation tomography by the adjoint-state method