RTM using effective boundary saving: A staggered grid GPU implementation

Overview of RTM and its computation

In the case of constant density, the acoustic wave equation is written as

$$\frac{1}{v^2(\textbf{x})}\frac{\partial^2 p(\textbf{x},t;\textbf{... ...t^2}=\nabla^2 p(\textbf{x},t;\textbf{x}_s)+f(t)\delta(\textbf{x}-\textbf{x}_s),$$ (1)

where $p(\textbf{x},t;\textbf{x}_s)$ is the wavefield excited by the source at the position $\textbf{x}=\textbf{x}_s$ , $v(\textbf{x})$ stands for the velocity in the media, $\nabla^2=\nabla\cdot\nabla=\partial_{xx}+\partial_{zz}$ , $f(t)$ denotes the source signature. For the convenience, we eliminate the source term hereafter and use the notation $\partial_u=\frac{\partial}{\partial u}$ and $\partial_{uu}=\frac{\partial}{\partial u^2}$ , $u=x,z$ . The forward marching step can be specified after discretization as

$$p^{k+1}=2p^{k}-p^{k-1}+v^2\Delta t^2 \nabla^2 p^{k}.$$ (2)

Based on the wave equation, the principle of RTM imaging can be interpreted as the cross-correlation of two wavefields at the same time level, one computed by forward time recursion, the other computed by backward time stepping (Symes, 2007). Mathematically, the cross-correlation imaging condition can be expressed as

$$I(\textbf{x})=\sum_{s=1}^{ns}\int_{0}^{t_{\max}}\mathrm{d}t \sum_{g=1}^{ng} p_s(\textbf{x},t;\textbf{x}_s)p_g(\textbf{x},t;\textbf{x}_g),$$ (3)

where $I(\textbf{x})$ is the migrated image at point ${\bf x}$ ; and $p_s(\cdot)$ and $p_g(\cdot)$ are the source wavefield and receiver (or geophone) wavefield. The normalized cross-correlation imaging condition is designed by incorporating illumination compensation:

$$I(\textbf{x})=\sum_{s=1}^{ns}\frac{\int_{0}^{t_{\max}}\mathrm{d}t... ...max}}\mathrm{d}t p_s(\textbf{x},t;\textbf{x}_s)p_s(\textbf{x},t;\textbf{x}_s)}.$$ (4)

There are some possible ways to do RTM computation. The simplest one may be just storing the forward modeled wavefields on the disk, and reading them for imaging condition in the backward propagation steps. This approach requires frequent disk I/O and has been replaced by wavefield reconstruction method. The so-called wavefield reconstruction method is a way to recover the wavefield via backward reconstructing or forward remodeling, using the saved wavefield snaps and boundaries. It is of special value for GPU computing because saving the data in device variables eliminates data transfer between CPU and GPU. By saving the last two wavefield snaps and the boundaries, one can reconstruct the wavefield of every time step, in time-reversal order. The checkpointing technique becomes very useful to further reduce the storage (Dussaud et al., 2008; Symes, 2007). Of course, it is also possible to avert the issue of boundary saving by applying the random boundary condition, which may bring some noises in the migrated image (Liu et al., 2013b; Clapp, 2009; Liu et al., 2013a; Clapp et al., 2010).

RTM using effective boundary saving: A staggered grid GPU implementation