Lowrank finite-differences and lowrank Fourier finite-differences for seismic wave extrapolation in the acoustic approximation

Numerical Examples

Our first example is wave extrapolation in a 2-D, smoothly variable velocity model. The velocity ranges between 500 and 1300 m/s, and is formulated as

$$v(x,z) = 500+1.2\times10^{-4}(x-800)^2+10^{-4}(z-500)^2;$$ (20)

$0\le x\le2560,\;0\le z\le2560.$ A Ricker-wavelet source with a 20 Hz dominant frequency ($f_{d}$) is located at the center of the model. The maximum frequency ($f_{max}$) is around 60 Hz. The amplitude corresponding to $f_{max}$ is about $10^{-5}$ of that of $f_{d}$. For numerical simulations based on this model, we use the same grid size: $\Delta x = 5$ m and $\Delta t = 2$ ms. We use $\alpha=v_{max} \Delta t/\Delta x$ to specify the stability condition and $\beta=v_{min}/(f_{max}\Delta x)$ as the dispersion factor, where $v_{max}$ and $v_{min}$ are the maximum and minimum velocities of the model. The dispersion factor $\beta$ indicates the number of sampling points for the minimum wavelength. For simulations with the above parameters, $\alpha \approx 0.52$ and $\beta \approx 1.67$.

wavfd7sn,wavel
Figure 4. Wavefield snapshot in a varable velocity field by: (a) conventional 4th-order FD method (b) Lowrank method.

It is easy to observe obvious numerical dispersions on the snapshot computed by the 4th-order FD method (Figure 4(a)). The lowrank FD method with the same order exhibits higher accuracy and fewer dispersion artifacts (Figure 5(a)). The approximation rank decomposition in this case is $N=3 M=4$, with the expected error of less than $10^{-4}$. Figure 5(b) displays the snapshot by the 10th-order LFD method with a larger time step: $\Delta t = 2.5$ ms, $\alpha \approx 0.65$. Note that the result is still accurate. However, the regular FD method becomes unstable in this case. For comparison, Figure 4(b) displays the snapshot by the lowrank method with the same time step. The approximation rank decomposition in this case is $N=M=3$, with the expected error of less than $10^{-4}$. Thanks to the spectral nature of the algorithm, the result appears accurate and free of dispersion artifacts.

wavapp7sn,wavapp10sn
Figure 5. Wavefield snapshot in a varable velocity field by: (a) the 4th-order lowrank FD method (b) the 10th-order lowrank FD method. Note that the time step is 2.5 ms and the LFD result is still accurate. However, the FD method is unstable in this case.

Next, we test the lowrank FD method in a complex velocity model. Figure 6 shows a part of the BP velocity model (Billette and Brandsberg-Dahl, 2004), which is a complicated model containing a salt body and sharp velocity contrasts on the flanks of the salt body. We use a Ricker-wavelet at a point source. The dominant frequency is 17 Hz ( $f_{max} \approx 54$). The horizontal grid size $\Delta x$ is 12.5 m, the vertical grid size $\Delta z$ is 12.5 m, and the time step is 1.5 ms. Thus $\alpha \approx 0.57$ and $\beta \approx 2.2$. The approximation rank decomposition in this case is $N=4 M=5$, with the expected error of less than $10^{-4}$. In this case, we adopt a disk-shaped compact scheme (8th-order) for LFD with a 4-point radius ( $\vert\mathbf{\xi}\vert \le 4$, $L=25$). Figure 7 displays a wavefield snapshot in the above velocity model. The snapshot is almost free of dispersions. This experiment confirms that the lowrank FD method is able to handle sharp velocity variations.

sub1
Figure 6. Portion of BP 2004 synthetic velocity model.

wavsnapabc
Figure 7. Wavefield snapshot by the 8th-order lowrank FD (compact scheme) in the BP Model shown in Figure 6.

Our next example is wave propagation in a TTI model with a tilt of $45 ^{\circ }$ and smooth velocity variation ($v_x$: 800-1225.41 m/s, $v_z$: 700-883.6 m/s). Figure 8(a) shows wavefield snapshots at different time steps by a 16th-order LFD operator in the TTI model. The space grid size is 5 m and the time step size is 2 ms. So $\alpha \approx 0.49$ and $\beta \approx 2.3$. The approximation rank decomposition in this case is $N=6 M=6$, with the expected error of less than $10^{-4}$. For TTI model, we adopt a high-order (16th order) LFD operator in order to reduce dispersions. The scheme is compact and shaped as a disk with a radius of 8 points ($L=99$).

snapshotlfd,snapshotlffd
Figure 8. Wavefield snapshots in a TTI medium with a tilt of $45 ^{\circ }$ by: (a) Lowrank FD method; (b) Lowrank FFD method. $v_x(x,z)=800+10^{-4}(x-1000)^2+10^{-4}(z-1200)^2$; $v_z(x,z)=700+10^{-4}(z-1200)^2$; $\eta =0.3$; $\theta =45 ^{\circ }$.

Song and Fomel (2011) showed an application of FFD method for TTI media. However, the example wavefield snapshot by FFD method still had some dispersion caused by the fact that the FD scheme in the FFD operator is derived from Taylor's expansion around zero wavenumber. It was apparent that 4th-order FD scheme is not accurate enough for TTI case and requires denser sampling per wavelength ( $\beta \approx 4.6$). We first apply lowrank approximation to the mixed-domain velocity correction term in FFD. The rank is $N=9 M=9$, with the expected error of less than $10^{-6}$. Then we propose to replace that 4th-order FD operator with an 8th-order LFD compact scheme. The scheme has the shape of a disk with a radius of 4 points ($L=25$), the same as the one for LFD in the above BP model. Figure 8(b) shows wavefield snapshots by the proposed LFFD operator. The time step size is 1.5 ms ( $\alpha \approx 0.37$). Note that the wavefront is clean and almost free of dispersion with $\beta \approx 2.3$. Because we use the exact dispersion relation, Equation 17 for TTI computation, there is no coupling of q-SV wave and q-P wave (Zhang et al., 2009; Grechka et al., 2004; Duveneck et al., 2008) in our snapshots by either LFD or LFFD methods.

vp0,vx0,yita0,theta0
Figure 9. Partial region of the 2D BP TTI model. a: $v_z$. b: $v_x$. c: $\eta$. d:$\theta$.

Next we test the LFD and LFFD methods in a complex TTI model. Figure 9(a)-9(d) shows parameters for part of the BP 2D TTI model (Shah, 2007). The dominant frequency is 15 Hz ( $f_{max} \approx 50$). The space grid size is 12.5 m and the time step is 1 ms. Thus $\alpha \approx 0.42$ and $\beta \approx 2.4$. The approximation rank decomposition for LFD method is $N=22 M=22$, with the expected error of less than $10^{-6}$. For FFD, $N=24 M=30$, with the expected error of less than $10^{-6}$. Both methods are able to simulate an accurate qP-wave field in this model as shown in Figure 10(a) and 10(b).

snapshotslfd,snapshotslffd
Figure 10. Scalar wavefield snapshots by LFD and LFFD methods in the 2D BP TTI model.

It is difficult to provide analytical stability analysis for LFD and LFFD operators. In our experience, the values of $\alpha$ are around 0.5 for 2D LFD and LFD methods appear to allow for a larger time step size than that of the LFFD method. In TTI case, the conventional FD method for acoustic TTI has known issues of instability caused by shear-wave numerical artifacts or sharp changes in the symmetry-axis tilting (Zhang et al., 2009; Grechka et al., 2004; Duveneck et al., 2008). Conventional methods may place limits on anisotropic parameters, smooth parameter models or include a finite shear-wave velocity to alleviate the instability problem(Fletcher et al., 2009; Yoon et al., 2010; Zhang and Zhang, 2008). Both LFD and LFFD methods are free of shear-wave artifacts. They require no particular bounds for anisotropic parameters and can also handle sharp tilt changes.

Lowrank finite-differences and lowrank Fourier finite-differences for seismic wave extrapolation in the acoustic approximation