Fast automatic wave-equation migration velocity analysis using encoded simultaneous sources

theory

I pose the velocity estimation problem as an optimization problem that tries to maximize the image stack power across the reflection angle, taking advantage of the fact that seismic events should be aligned and hence most constructively stacked in the angle domain, if migrated using an accurate velocity model (Soubaras and Gratacos, 2007). Instead of solving it as a maximization problem, I actually solve it as a minimization problem that minimizes the negative image stack power. Because the reflection-angle stacked section is equivalent to the zero-subsurface offset image, the objective function that I use to minimize is therefore defined as follows:

$$J = -\frac{1}{2}\sum_{{\bf x}}m_{\rm mig}^2({\bf x}),$$ (1)

where $m_{\rm mig}({\bf x})$ is the zero-subsurface-offset image at image point ${\bf x}=(x,y,z)$ , obtained by crosscorrelating the forward propagated source wavefield with the backward propagated receiver wavefield as follows:

$$m_{\rm mig}({\bf x}) = \sum_{{\bf x}_s}\sum_{\omega}S({\bf x},{\bf x}_s,\omega)R({\bf x},{\bf x}_s,\omega),$$ (2)

where $S({\bf x},{\bf x}_s,\omega)$ and $R({\bf x},{\bf x}_s,\omega)$ are the source and receiver wavefield at image point ${\bf x}$ , respectively, for a source located at ${\bf x}_s=(x_s,y_s,0)$ and at angular frequency $\omega$ . If a one-way extrapolator is used, $S$ and $R$ satisfy the following one-way wave equations:

$$\left\{ \begin{array}{l} \left( \frac{\partial}{\partial z}+i\sqr... ...0,{\bf x}_s,\omega) = \delta(x-x_s,y-y_s){f_s^{*}(\omega)} \end{array} \right.,$$ (3)

and

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

where $^{*}$ denotes taking the adjoint; $v({\bf x})$ is the velocity at image point ${\bf x}$ ; $f_s(\omega)$ is the source signature; $\delta(\cdot)$ is the Dirac delta function; $\nabla^2=\frac{\partial^2}{\partial x^2}+\frac{\partial ^2}{\partial y^2}$ is the Laplacian operator. $Q$ is the observed data mapped onto the computation grid, which is defined as follows:

$$Q(x,y,{\bf x}_s,\omega) = \sum_{{\bf x}_r} W({\bf x}_r,{\bf x}_s)\delta(x-x_r,y-y_r)d_{\rm obs}({\bf x}_r,{\bf x}_s,\omega),$$ (5)

where $d_{\rm obs}$ is the observed data recorded at ${\bf x}_r=(x_r,y_r,0)$ due to a source at ${\bf x}_s$ ; $W({\bf x}_r,{\bf x}_s)$ is the acquisition mask operator, which contains ones where we record data and zeros where we do not.

Since flat angle gathers generate the most coherent stack, the negative image-stack-power minimization objective function defined by equation 1 is intuitive to understand. Objective function 1, however, has an alternative interesting interpretation as shown in Appendix A, which proves that under the least-squares assumption, minimization of objective function 1 is equivalent to the data-domain Born wavefield inversion, which minimizes the differences between the modeled and observed primary reflections.

Objective function 1 is usually minimized using local optimization techniques, which require explicit calculation of the gradient. The gradient is obtained by taking the derivative of $J$ with respect to velocity $v({\bf y})$ (${\bf y}$ is the velocity coordinates) as follows:

$$g({\bf y}) = \frac{\partial J}{\partial v({\bf y})} = - \sum_{{\b... ...}\frac{\partial m_{\rm mig}({\bf x})}{\partial v({\bf y})}m_{\rm mig}({\bf x}),$$ (6)

where the sensitivity kernel, or tomographic operator, $\frac{\partial m_{\rm mig}({\bf x})}{\partial v({\bf y})}$ , can be easily obtained as follows:

$$\frac{\partial m_{\rm mig}({\bf x})}{\partial v({\bf y})}= \sum_{... ...mega) \frac{\partial R({\bf x,{\bf x}_s,\omega})}{\partial v({\bf y})} \right).$$ (7)

Note the summations over ${\bf x}_s$ in equations 2 and 7. This means that the computation for the image $m_{\rm mig}$ and the gradient $g$ needs to be carried out for each source independently, resulting in a cost proportional to the number of sources. The gradient $g$ is usually calculated using the adjoint-state technique without explicitly constructing the sensitivity kernel (Sava and Vlad, 2008; Shen, 2004; Tang et al., 2008).

For encoded simultaneous-source WEMVA, the objective function to be minimized is defined as follows:

$${\widetilde J} = -\frac{1}{2}\sum_{{\bf x}}{\widetilde m}_{\rm mig}^2({\bf x}),$$ (8)

where the zero-subsurface-offset image ${\widetilde m}_{\rm mig}$ is obtained by crosscorrelating the encoded source wavefield, ${\widetilde S}$ , and the encoded receiver wavefield, ${\widetilde R}$ , as follows:

$${\widetilde m}_{\rm mig}({\bf x}) = \sum_{\omega}{\widetilde S}({\bf x},\omega){\widetilde R}({\bf x},\omega).$$ (9)

The encoded source and receiver wavefields satisfy the following one-way wave equations:

$$\left\{ \begin{array}{l} \left( \frac{\partial}{\partial z}+i\sqr... ...ta(x-x_s,y-y_s)f_s^{*}(\omega)\alpha^{*}({\bf x}_s,\omega) \end{array} \right.,$$ (10)

and

$$\left\{ \begin{array}{l} \left( \frac{\partial}{\partial z}+i\sqr... ..._{{\bf x}_s}Q(x,y,{\bf x}_s,\omega)\alpha({\bf x}_s,\omega)\end{array} \right.,$$ (11)

where $\alpha({\bf x}_s,\omega)$ is the phase encoding function. In this paper, I mainly focus on random phase encoding, therefore $\alpha$ is defined as follows:

$$\alpha({\bf x}_s,\omega) = e^{i\gamma({\bf x}_s,\omega)},$$ (12)

where $i = \sqrt{-1}$ , and $\gamma({\bf x}_s,\omega)$ is a uniformly distributed random sequence from 0 to $2\pi$ . Tang (2011) shows that with this choice of random phase function, $\alpha$ has a zero expectation. Note that the source encoding can be applied to data recorded from arbitrary types of acquisition geometries. The simultaneous-source migrated image ( ${\widetilde m}_{\rm mig}$ ) will always converge to the separate-source migrated image ( $m_{\rm mig}$ ) as long as the encoding function satisfies $\alpha^{*}({\bf x}_s,\omega)\alpha({\bf x}_s',\omega)\approx\delta({\bf x}_s-{\bf x}_s')$ .

The gradient of objective function 8 is

$${\widetilde g}({\bf y}) = \frac{\partial {\widetilde J}}{\partial... ...e m}_{\rm mig}({\bf x})}{\partial v({\bf y})}{\widetilde m}_{\rm mig}({\bf x}),$$ (13)

where the tomographic operator, $\frac{\partial {\widetilde m}_{\rm mig}({\bf x})}{\partial v({\bf y})}$ , in the encoded-source domain is defined as follows:

$$\frac{\partial {\widetilde m}_{\rm mig}({\bf x})}{\partial v({\bf... ...a) \frac{\partial {\widetilde R}({\bf x,\omega})}{\partial v({\bf y})} \right).$$ (14)

Note that equations 9 and 14 do not have a summation over the sources. Therefore, the cost of computing the image ${\widetilde m}_{\rm mig}$ and the gradient ${\widetilde g}$ is independent of the number of sources, as opposed to the separate-source case. The gradient ${\widetilde J}$ is also calculated using the adjoint-state technique using encoded simultaneous sources (Tang et al., 2008).

Although the computational cost of WEMVA is significantly reduced, encoded simultaneous sources add crosstalk artifacts into the gradient. This becomes clear if we express the encoded source and receiver wavefield as follows using the fact that wavefield propagation is linear with respect to sources:

$${\widetilde S}({\bf x},\omega) = \sum_{{\bf x}_s}S({\bf x},{\bf x}_s,\omega)\alpha^{*}({\bf x}_s,\omega),$$ (15)

and

$${\widetilde R}({\bf x},\omega) = \sum_{{\bf x}_s}R({\bf x},{\bf x}_s,\omega)\alpha({\bf x}_s,\omega).$$ (16)

Substituting equations 15 and 16 into 14 and using the fact that $\alpha^{*}({\bf x}_s,\omega)\alpha({\bf x}_s',\omega)=1$ if ${\bf x}_s={\bf x}_s'$ yield

$${\widetilde g}({\bf y}) = g({\bf y}) + g_c({\bf y}),$$ (17)

where $g_c$ is the crosstalk:

$$g_c({\bf y}) = \sum_{\omega} \sum_{{\bf x}_s}\sum_{{\bf x}_s'\ne{\bf x}_s} \left(\frac{\partial S({\bf x},{\bf x}_s,\omega)}{\partial v({\bf... ...\omega)\frac{\partial R({\bf x},{\bf x}_s',\omega)}{\partial v({\bf y})}\right)$$

$$\times \alpha^{*}({\bf x}_s,\omega) \alpha({\bf x}_s',\omega).$$ (18)

A way to mitigate the influence of crosstalk is to change the random encoding function at each iteration (Krebs et al., 2009), so that the crosstalk will be destructively stacked over WEMVA iterations and consequently converge to zero because it has a zero expectation. It is important to point out that regeneration of the random code will result in the objective function (equation 8) changing at each iteration. Therefore, the objective function may not be monotonically decreasing over iterations, as opposed to the case in conventional separate-source WEMVA. The optimization algorithm using encoded simultaneous sources is summarized in Algorithm 1.

Fast automatic wave-equation migration velocity analysis using encoded simultaneous sources