Migration velocity analysis for anisotropic models

Migration Velocity Analysis for Anisotropic parameters

Anisotropic MVA is a non-linear inversion process that aims to find the background anisotropic model that minimizes the residual image $\Delta {\bf I}$ . The residual image is derived from the background image ${\bf I}$ , which is computed with the current background model. To form the image, both the source and receiver wavefields are downward continued using the one-way wave equations. Assuming that the shear velocity is much smaller than the P-wave velocity, one way of formulating up-going and down-going one-way acoustic wave equations for VTI is shown as follows (Shan, 2009):

$$\left( \frac{\partial}{\partial z} \mp i\Lambda \right)P = 0 ,$$ (1)

where $P=P(x,y,z,\omega)$ is the wavefield in the space-frequency domain and $\Lambda$ describes the dispersion relationship in terms of P-wave vertical slowness $s_0$ and Thomsen parameters $\epsilon$ and $\delta$ (Thomsen, 1986):

$$\Lambda = \omega s_0\sqrt{\frac{\omega^2 s_0^2 - (1+2\epsilon)\vert{\bf k}\vert ^2}{\omega ^2 s_0^2 - 2 (\epsilon - \delta) \vert{\bf k}\vert ^2}},$$ (2)

where ${\bf k}=(k_x,k_y)$ is the spatial wavenumber vector.

Many authors (Alkhalifah and Tsvankin, 1995; Tsvankin and Thomsen, 1994) have shown that P-wave traveltime can be characterized by the NMO slowness, $s_n$ , and the anellipticity parameter $\eta$ . Therefore, the one-way wave-equation in terms of $s_n$ , $\eta$ and $\delta$ is:

$$\left( \frac{1}{\sqrt{1+2\delta}}\frac{\partial}{\partial z} \mp i\Lambda' \right)P = 0$$ (3)

where

$$\Lambda' = \omega s_n\sqrt{1-\frac{\vert{\bf k}\vert ^2}{\omega ^2 s_n^2 - 2 \eta \vert{\bf k}\vert ^2}}.$$ (4)

Notice that in the dispersion relationship in Equation 3, $\delta$ and the derivative in depth $\frac{\partial}{\partial z}$ , are coupled with each other. This is a theoretical proof of the well-accepted observation that $\delta$ cannot be determined by the surface seismic data. To constrain this parameter, we need well information (e.g. checkshots) to add the depth dimension into the inversion. Now, if we apply the change of variables

$$d \bar{z} = \sqrt{1+2\delta}dz$$ (5)

and neglect the derivatives of $\delta$ , Equation 3 becomes

$$\left( \frac{\partial}{\partial \bar{z}} \mp i\Lambda' \right)P = 0.$$ (6)

We can therefore formulate the image-space migration velocity analysis problem with NMO slowness $s_n$ and anisotropic parameters $\eta$ and $\delta$ , but we invert only for $s_n$ and $\eta$ assuming $\delta$ model is known from other source of information.

Notice that when $\eta=0$ , the dispersion relationship ( equation 4) is the same as the isotropic dispersion relationship, and the corresponding one-way wave equation ( equation 6) is almost the same as for the isotropic case, except for a depth stretch caused by $\delta$ . In other words, an elliptic anisotropic wavefield inversion is almost equivalent to an isotropic wavefield inversion. Plessix and Rynja (2010) reached the same conclusions for full-waveform inversion (FWI). Figure 1 compares the original NMO velocity to the stretched NMO velocity. Notice that the geological features are stretched downward for positive $\delta$ . Because we ignore $\delta$ in the inversion, we expect the inverted NMO velocity to have more similarity to the stretched NMO velocity than to the original one.

nmo
Figure 1. (a) Original NMO velocity for the anisotropic Hess model; (b) Stretched NMO velocity according to $\delta$ . The overlaid box denotes the part of model that we work on in the numerical test.

In general, the residual image is defined as (Biondi, 2008)

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

where ${\bf F}$ is a focusing operator. In the least-square sense, the tomographic 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.$$ (8)

To perform MVA for anisotropic parameters, we first need to extend the tomographic operator from the isotropic medium (Shen, 2004; Sava, 2004; Guerra et al., 2009) to the anisotropic medium. We define the wave-equation tomographic operator T for anisotropic models as follows:

$${\bf T} = \left.\frac{\partial{\bf I}}{\partial{\bf m}}\right\vert _{{\bf m}=\widehat{\bf m}}$$

$$= \left.\frac{\partial{\bf I}}{\partial{\bf s_n}}\right\vert _{{\bf... ...{\partial{\bf I}}{\partial{\bf \eta}}\right\vert _{{\bf\eta}=\widehat{\bf\eta}}$$ (9)

where m is the anisotropy model, which in this case includes NMO slowness ${\bf s_n}$ and anellipticity parameter ${\bf\eta}$ ; $\widehat{\bf m}$ is the background anisotropy model, consisting of the background NMO slowness $\widehat{\bf s_n}$ and background anellipticity $\widehat{\bf\eta}$ ; and I is the image. This wave equation tomographic operator T is a linear operator that relates the model perturbation ${\Delta {\bf m}}$ to the image perturbation ${\Delta {\bf I}}$ as follows:

$$\Delta {\bf I} = {\bf T} \Delta {\bf m}.$$ (10)

In the shot-profile domain, both source and receiver wavefields are downward continued using the one-way wave equation (6):

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

and

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

where $D({\bf x},{\bf x}_s)$ is the source wavefield at the 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)$ is the receiver wavefield at the image point ${\bf x}$ with the source located at ${\bf x}_s$ ; $f_s$ is the source signature, and $f_s\delta({\bf x}-{\bf x}_s)$ defines the point-source function at ${\bf x_s}$ , which serves as the boundary condition of equation 11; and $Q(x,y,z=0,{\bf x}_s)$ is the recorded shot gather for the shot located at ${\bf x}_s$ , which serves as the boundary condition of equation 12.

The dispersion relationship in equation (4) can be approximated with a rational function by Taylor series and Padé expansion analysis (Shan, 2009):

$$\Lambda' = \omega s_n\left(1 - \frac{a \vert{\bf k}\vert ^2}{\omega ^2 s_n^2 - b \vert{\bf k}\vert ^2}\right),$$ (13)

where, to the second order of the expansion, $a = 0.5, b = 0.25+2\eta$ . Equation (13) using binomial expansion can be further expanded to polynomials:

$$\Lambda' = \omega s_n - \frac{a}{\omega^2 s_n^2} \vert{\bf k}\vert ^2 - \frac{ab}{\omega^4 s_n^4} \vert{\bf k}\vert^4.$$ (14)

Now it is straightforward to take the derivative of $\Lambda'$ with respect to $s_n$ and $\eta$ .

The background 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)} U({\bf x}+{\bf h},{\bf x}_s),$$ (15)

where the overline stands for complex conjugate, and ${\bf h}=(h_x,h_y,h_z)$ is the subsurface half-offset. Perturbing the wavefields in equation (15) and ignoring the higher-order term, we can get the perturbed image as follows:

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

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

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

To evaluate the adjoint tomographic operator ${\bf T}^{\ast}$ , which maps from the image perturbation to the model perturbation, we first compute the wavefield perturbation from the image perturbation using the adjoint imaging condition:

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

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

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

$$\left\{ \begin{array}{l} \left( \frac{\partial}{\partial z}+i\Lam... ...\bf m}^{\ast}({\bf x})\\ \Delta D(x,y,z=0,{\bf x}_s) = 0 \end{array} \right. ,$$ (18)

and

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

where ${\bf m}$ is the row vector $[{\bf s_n} ~ {\bf\eta}]$ .

During the inversion, the model perturbation is unknown, and in fact must be estimated. Therefore, we obtain the image perturbation by applying a focusing operator (equation 7) to the current background image. Then the perturbed image is convolved with the background wavefields to get the perturbed wavefields (equation 17). The scattered wavefields are obtained by applying the adjoint of the one-way wave-equations (18) and (19). Finally, the model-space gradient is obtained by cross-correlating the upward propagated scattered wavefields with the modified background wavefields [the terms in the parentheses on the right-hand sides of equations (18) and (19)].

Migration velocity analysis for anisotropic models