Wave-equation migration velocity analysis for VTI media using optimized implicit finite difference

Impulse response of the image-space wave-equation tomography operator

The core of anisotropic image-space wave-equation migration velocity analysis is the tomography operator ${\bf T}$ , which relates the perturbation in the anisotropic models ( $\Delta {\bf m}$ ) to to the perturbation in the image ($\Delta I$ ) and vice versa. Namely,

$$\Delta I = {\bf T} \Delta {\bf m},$$ (7)

$$\Delta {\bf m} = {\bf T'} \Delta I,$$ (8)

where ${\bf m} = [v_v~\eta]$ .

I refer the readers to Li and Biondi (2010) for a detailed derivation for the tomographic operator.

A different approximation to the exact dispersion relation leads to a different perturbed wave fields due to a perturbation in the model parameters. When the only available data come from surface seismic surveys, parameter $\delta$ is the least constrained (Plessix and Rynja, 2010; Li and Biondi, 2011b). Therefore, I assume the $\delta$ model is perfectly obtained from other sources of data and keep it fixed throughout the inversion. I will invert for $v_v$ and $\eta$ in this study.

In the downward extrapolation, the wavefield at the next depth ($P_{z+1}$ ) can be computed from the wavefield at the current depth ($P_z$ ) according to the following equation:

$$P_{z+1} = P_z e^{i k_z dz},$$ (9)

where $i = \sqrt{-1}$ , $dz$ is the extrapolation distance in depth and $k_z$ can be obtained from the first-order approximation of the dispersion relation 5:

$$k_z = \frac{w}{v_v} \left ( 1 - \frac{\alpha \frac{k_r^2}{(w/v_v)^2}}{1 - \beta \frac{k_r^2}{(w/v_v)^2} } \right ).$$ (10)

Dispersion relation 10 can be further simplified to polynomials using Taylor expansion:

$$k_z = \frac{w}{v_v} \left ( 1 - \alpha \frac{k_r^2}{(w/v_v)^2} \left (1 + \beta \frac{k_r^2}{(w/v_v)^2} \right ) \right )$$ (11)

$$= \frac{w}{v_v} \left ( 1 - \alpha \frac{k_r^2}{(w/v_v)^2} - \alpha \beta \frac{k_r^4}{(w/v_v)^4} \right ).$$

Therefore, the perturbed wavefield is

$$\Delta P_{z+1} = e^{ik_zdz} i dz P_z \Delta k_z,$$ (12)

with

$$\Delta k_z = \frac{\partial k_z}{\partial v_v} \Delta v_v + \frac{\partial k_z}{\partial \eta} \Delta \eta,$$ (13)

$$\frac{\partial k_z}{\partial v_v} = - \frac{w}{v_v^2} \left( 1 + \alpha \frac{k_r^2}{(w/v_v)^2} + 3 \alpha \beta \frac{k_r^4}{(w/v_v)^4} \right),$$ (14)

and

$$\frac{\partial k_z}{\partial \eta} = - \frac{w}{v_v} \left( \frac... ... \frac{\partial \beta}{\partial \eta} \right ) \frac{k_r^4}{(w/v_v)^4} \right).$$ (15)

Since the finite difference parameters $\alpha$ and $\beta$ are obtained by optimization, the derivatives in Equation 15 are obtained numerically by taking derivatives along the $\eta$ axis in Figure 2. The tables of the derivatives of the coefficients with respect to $\eta$ are shown in Figure 3.

dr-coef
Figure 3. (a) Table for $\frac {\partial \alpha }{\partial \eta }$ and (b) table for $\frac {\partial \beta }{\partial \eta }$ at background $\eta$ and $\delta$ locations.

I test the implementation of the adjoint tomographic operator using this optimized implicit finite difference scheme in a homogeneous background VTI medium with $v_v = 2000$ km, $\eta=0.09$ and $\delta=0.05$ . The synthetic data is produced by Born modeling with a horizontal reflector at the depth of 1500 km. The input of the adjoint tomographic operator is a spike in the image space $\Delta I = \delta(x,y,z=1500)$ . The dominant frequency of the source wavelet is $20$ Hz, and the samplings in all directions are $10$ m.

2dkernel
Figure 4. 2D impulse responses for vertical velocity (left column) and $\eta$ (right column). Top row: zero offset impulse responses; middle row: impulse responses when source-receiver offset is $4$ km; bottom row: summation of the two rows above.

I first test the adjoint operator in 2D. A source and receiver pair is collocated at $x=y=z=0$ . The top row in Figure 4 shows the back-projected vertical velocity $v_v$ gradient and $\eta$ gradient when source-receiver offset is zero. These back projections are often referred as banana-donut kernels in the literature when transmission waves are under study (eg. Marquering et al. (1998); Rickett (2000); Marquering et al. (1999)). Similar reflection tomography sensitivity kernel analysis for isotropic WEMVA operator can be found in Sava (2004) and Xie and Yang (2009).

Compared with the $\eta$ gradient, the $v_v$ gradient has a nearly uniform strength with depth, while the $\eta$ gradient fades away as the wavepath moves away from the source and the receiver location. Also, the dominant energy of the $\eta$ gradient points to the opposite direction of the $v_v$ gradient points. In fact, the $\eta$ gradient is not reliable and should be ignored because the wave that travels in the vertical direction is not sensitive to $\eta$ .

When the source-receiver offset is $4$ km, the gradients are shown in the middle row in Figure 4. Clearly, the back projections are spread along the wavepaths from the source to the perturbed image point and from the perturbed image point to the receiver. In this case, the gradients in both $v_v$ and $\eta$ point in the same direction. Comparing the gradients in the cases of zero and nonzero offset, one can see that the vertical waves are more sensitive to $v_v$ , and the waves traveling at a large angle ($36^\circ$ to the vertical in this case) are more sensitive to $\eta$ . The bottom row in Figure 4 shows the summation of the gradients in these two cases, and confirms these observations.

The 3D extension of this method is straightforward. The sensitivity kernels for $v_v$ and $\eta$ in 3D are shown in Figures 5 and 6. A source and receiver pair with $4$ km offset are located at $y=0$ . The 3D sensitivity kernels carry the same characteristics as the 2D kernels, only expanding to the crossline direction.

3dkernel-vel-new
Figure 5. 3D $v_v$ kernel.

3dkernel-eta-new
Figure 6. 3D $\eta$ kernel.

Wave-equation migration velocity analysis for VTI media using optimized implicit finite difference