Correlation-based wave-equation migration velocity analysis

Method

The first step in evaluating a tomographic operator is to linearize the image $\mathbf I$ around the background slowness $\mathbf s_0$ , as follows:

$$\mathbf I = \mathbf I_0 + \frac{\partial \mathbf I}{\partial \mathbf s}\vert_{s_0}(\mathbf s - \mathbf s_0) + ...,$$ (1)

where $\mathbf I_0$ is the background image and $\mathbf s$ is the slowness model. By neglecting the higher order terms in the image series, we can define the tomographic operator as follows:

$$\Delta \mathbf I = \frac{\partial \mathbf I}{\partial \mathbf s}\vert_{s_0} \Delta \mathbf s = \mathbf T \Delta \mathbf s,$$ (2)

where $\mathbf T$ is the tomographic operator. Now, we use the conventional imaging condition as follows:

$$I(\mathbf x, \mathbf h) = \sum_{\omega,\mathbf x_s,\mathbf x_r} G... ...^*(\mathbf x + \mathbf h,\mathbf x_r,\omega) d(\mathbf x_r,\mathbf x_s,\omega),$$ (3)

where $I$ is the image, $G$ is the Green's function, $d$ is the surface data, $\omega$ is frequency, $\mathbf x_s$ and $\mathbf x_s$ are the source and receiver coordinates, and $\mathbf h$ is the subsurface offset. To evaluate the tomographic operator $\mathbf T$ , I take the derivative of the imaging condition as follows:

$$\Delta I(\mathbf x, \mathbf h) = \sum_{\mathbf y} \frac{\partial I(\mathbf x, \mathbf h)}{\partial s(\mathbf y)} \Delta s(\mathbf y)$$

$$= \sum_{\omega,\mathbf x_s,\mathbf x_r,\mathbf y} \left\lbrace-2 ... ...hbf h,\mathbf x_r,\omega) d(\mathbf x_r,\mathbf x_s,\omega) \Delta s(\mathbf y)$$

$$+ \sum_{\omega,\mathbf x_s,\mathbf x_r,\mathbf y} \left\lbrace-2 ... ...bf y,\mathbf x_r,\omega) d(\mathbf x_r,\mathbf x_s,\omega) \Delta s(\mathbf y),$$ (4)

where $y$ is the slowness coordinate. The full derivation of the tomographic operator is presented in Almomin and Tang (2010).

After defining the tomographic operator, I use a cross-correlation function to estimate image perturbations:

$$f(\zeta, \gamma; \mathbf x) = \sum_{z} I_{\text {obs}}(z, \gamma; \mathbf x) I_{\text {cal}}(z + \zeta, \gamma; \mathbf x),$$ (5)

where $\zeta$ is the lag, $\gamma$ is the reflection angle, $\mathbf x$ is the surface coordinates, $z$ is depth, $I_{\text {obs}}$ is the angle-domain image using the observed data, and $I_{\text {cal}}$ is the angle-domain image using the calculated data, which is modeled with the background slowness. $I_{\text {cal}}$ is always going to be flat, since I create Born-modeled data using a reference model as the reflectivity and the background slowness and then migrate that data using the same background slowness. Next, I define $\xi$ to be the lag that maximizes the correlation function. Therefore, the derivative of the correlation function vanishes at that lag, as follows:

$$\mathbf g = \frac{\partial \mathbf f}{\partial \zeta}\vert_{\xi}=... ...t {obs}}(z, \gamma; \mathbf x) I_{\text {cal}}(z + \xi, \gamma; \mathbf x) = 0.$$ (6)

We can now use $\xi$ as our measure of the residual to minimize, casted as the following objective function:

$$J(\mathbf s) = \frac{1}{2} \Vert \xi(\gamma,\mathbf x) \Vert ^{2}.$$ (7)

Then, we take the derivative of the objective function with respect to slowness as follows:

$$\nabla J = \left( \frac{\partial \xi}{\partial \mathbf s} \right)^* \xi,$$ (8)

where $^*$ indicates an adjoint. By using the chain rule of differentiation, I relate the derivative of the maximum lag $\xi$ with respect to $\mathbf s$ to the derivative of the correlation function with respect to $\mathbf s$ as follows:

$$\left( \frac{\partial \xi}{\partial \mathbf s} \right)^*= \left( ... ...mathbf g}{\partial \mathbf s} \frac{\partial \xi}{\partial \mathbf g}\right)^*.$$ (9)

The second partial derivative in equation (9) is just a scalar that balances the energy between surface locations. The first partial derivative with respect to slowness can be calculated using equation (6) as follows:

$$\left( \frac{\partial \mathbf g}{\partial \mathbf s} \right)^* = \sum_{z} \left( \frac{\partial I_{\text {obs}}(z, \gamma; \math... ...ight)^* \frac{\partial}{\partial z} I_{\text {cal}}(z + \xi, \gamma; \mathbf x)$$

$$+ \sum_{z} \left( \frac{\partial I_{\text {cal}}(z + \xi, \gamma;... ...s} \right)^* \frac{\partial}{\partial z} I_{\text {obs}}(z, \gamma; \mathbf x).$$ (10)

The first tomographic operator in equation (10) can be computed as I described in equation (4). However, the second tomographic operator depends on how $I_{\text {cal}}$ is computed. If a fixed-reflectivity model is used, such as well data, and only the background slowness is updated, then this derivative will be very small and could be ignored, since changing the slowness updates does not change the reflectivity estimate. On the other hand, if we allow the modeling reflectivity to change location, i.e. we update the reflectivity model as we iterate, then this operator could have a significant component. However, evaluating this operator is much more expensive than the first tomographic operator, since it is a cascade of three operators. Therefore, I will assume that the first tomographic operator is sufficient and ignore the second term.

Correlation-based wave-equation migration velocity analysis