Lowrank one-step wave extrapolation for reverse-time migration

Wave propagation in TTI media

vp,vx,eta,theta
Figure 11. BP-2007 anisotropic benchmark model. (a) P-wave phase velocity along the axis of symmetry $V_z$ . (b) P-wave phase in the symmetry plane $V_x$ . (c) Anellipticity parameter $\eta$ . (d) Tilt angle of the symmetry axis $\theta$ .

wavesnaps-4,wavesnaps-10,wavesnaps-20,wavesnaps-30
Figure 12. Wavefield snapshots taken at different times showing P-wave propagation through the 2007 BP TTI model. Propagation-direction-dependent absorbing boundary conditions allow little energy to be reflected, even at wide-incident angles. The four wavefield snapshots calculated using different time step sizes are almost identical: (a) $\Delta t = 4\;ms$ ; (b) $\Delta t=10\;ms$ ; (c) $\Delta t=20\;ms$ ; (d) $\Delta t=30\;ms$ .

swavesnaps-2,swavesnaps-4,swavesnaps-10,swavesnaps-20
Figure 13. Wavefield snapshots taken at different times showing S-wave propagation through the 2007 BP TTI model. The four wavefield snapshots calculated using different time step sizes are almost identical: (a) $\Delta t=2\;ms$ ; (b) $\Delta t = 4\;ms$ ; (c) $\Delta t=10\;ms$ ; (d) $\Delta t=20\;ms$ .

To show that the proposed method handles anisotropic media accurately, we propagate qP- and qSV-waves in a portion of the 2007 anisotropic benchmark model from BP. The model has a horizontal spacing of $50\;m$ and a vertical spacing of $25\;m$ . As shown in Figure 11, the model exhibits a strong tilted transverse isotropy (TTI). We restrict the modeled area between $34\;km$ and $60\;km$ . Since the original model does not include S-wave velocity, we compute a vertical S-wave velocity profile simply as $0.3$ times the original vertical P-wave velocity.

To calculate the qP- and qSV-wave phase velocities, we employ the exact formulas defined by Gassmann (1964):

$$V^2_{P}(\mathbf{n},\mathbf{x})\vert\mathbf{k}\vert^2 = \frac{1}{2} \left[(c_{11}+c_{55}) \hat{k}_x + (c_{33}+c_{55}) \hat{k}_z\right] +$$

$$\frac{1}{2} \sqrt{\left[(c_{11}-c_{55}) \hat{k}_x - (c_{33}-c_{55}) \hat{k}_z\right]^2 + 4 (c_{13}+c_{55})^2 \hat{k}_x \hat{k}_z}\;;$$ (31)

$$V^2_{SV}(\mathbf{n},\mathbf{x})\vert\mathbf{k}\vert^2 = \frac{1}{2} \left[(c_{11}+c_{55}) \hat{k}_x + (c_{33}+c_{55}) \hat{k}_z\right] -$$

$$\frac{1}{2} \sqrt{\left[(c_{11}-c_{55}) \hat{k}_x - (c_{33}-c_{55}) \hat{k}_z\right]^2 + 4 (c_{13}+c_{55})^2 \hat{k}_x \hat{k}_z}\;,$$ (32)

where $c_{11},c_{13},c_{33}$ and $c_{55}$ are the density-normalized components of the elastic tensor, $\mathbf{k}$ is the wavenumber vector, $\hat{k}_x$ and $\hat{k}_z$ are the wavenumbers evaluated in the rotated coordinate aligned with the symmetry axis:

$$\begin{array}{*{20}c} \hat{k}_x=k_x\cos{\theta}+k_z\sin{\theta}\;, \hat{k}_z=k_z\cos{\theta}-k_x\sin{\theta}\;, \end{array}$$ (33)

where $\theta$ is the tilt angle measured with respect to the horizontal direction. Alternatively, the qP-wave phase velocity can be calculated using accurate approximate formulas (Alkhalifah, 1998,2000; Sripanich and Fomel, 2015; Fomel, 2004).

To mimic an unbounded medium, we apply the absorbing boundary condition in equation 26. First, a time step size of $4\;ms$ for qP-wave and $2\;ms$ for qSV-wave is used. A Ricker wavelet source with a peak frequency of $16\;Hz$ is located at $X=47.3\;km$ and $Z=0.175\;km$ . Figures 12a and 13a illustrates wavefield snapshots taken at different times overlaid on top of each other. Because of the accumulative effect of the large time step, model complexity and extra anisotropy introduced by the absorbing boundary conditions, the lowrank approximation took $N=19$ for qP-wave and $N=15$ for qSV-wave for an accuracy threshold of $\epsilon=10^{-4}$ .

To test the effect of using large $\Delta t$ for wave propagation, we choose a series of increasing $\Delta t$ while enforcing the same accuracy requirement ( $\epsilon=10^{-4}$ ). The qP-wave wavefields (Figure 12) generated using $\Delta t=4\;ms, 10\;ms, 20\;ms$ and $30\;ms$ are almost identical. Similarly, the qSV-wave wavefields (Figure 13) generated using $\Delta t=2\;ms, 4\;ms, 10\;ms$ and $20\;ms$ have very little difference. This shows that the stability and accuracy of wave extrapolation is not compromised by using very large time steps, even beyond the Nyquist limit of $10\;ms$ .

To investigate the effects of time step size ($\Delta t$ ) and accuracy level ($\epsilon$ ) of lowrank approximation on the computational cost, which is directly controlled by the rank of the approximation, we conduct a series of experiments using different values of $\Delta t$ and $\epsilon$ for qP-wave propagation. No absorbing boundary condition is applied in this test. Table 1 lists rank $N$ required by different time step sizes to achieve different accuracy levels. We can observe that the lowrank algorithm requires a higher rank to maintain the same level of accuracy when the time step size increases. On the other hand, with the same time step size, a higher rank may be needed to achieve a higher accuracy. In practice, we find that $\epsilon=10^{-4}$ is small enough to guarantee accurate wave kinematics. Since making a larger step in time requires a smaller number of steps, the optimal time step size $\Delta t$ can be selected to be the one that minimizes the total number of Fast Fourier Transforms (FFTs) required to accomplish the modeling or imaging task at hand. In practice, we found that this minimization problem usually has a straightforward solution: the larger $\Delta t$ is, the less computation is required. Therefore, the optimal solution will be the largest $\Delta t$ that satisfies other constraints, such as the imaging condition when performing RTM.

$\epsilon$ $\Delta t$ (ms)	0.5	1	2	4	8	16	32	64
$10^{-3}$	4	7	7	10	14	19	33	53
$10^{-4}$	7	9	13	15	22	30	46	76
$10^{-5}$	11	14	18	23	30	43	59	84

Table 1. Rank $N$ calculated from the lowrank approximation of the propagation matrix for qP-wave propagation in the BP 2007 TTI model using different time step sizes $\Delta t$ and accuracy levels $\epsilon$ .

Lowrank one-step wave extrapolation for reverse-time migration