Adaptive stabilization

Mathematically speaking, the compensated viscoacoustic wave equations 2 and 3 are severely ill-posed due to the presence of the compensating term $+\tau \partial_t (-\nabla^2)^{\gamma+1/2} p(\mathbf{x},t)$. That is to say, amplitude compensation is a nonstationary process with energy exponentially amplified over travel time, which boosts high-frequency ambient noise and can even result in numerical instability. In the package, we apply an adaptive stabilization scheme for $Q$-RTM to suppress unwanted high-frequency artifacts, which is discussed in my previous work (Wang et al., 2018c). Here we briefly summarize the process.

I derive a k-space Green's function of equation 1 by enforcing a point source at time $t=t_0$ and $\mathbf{x}=\mathbf{x}_s$. The time-space harmonic Green's function $G(\mathbf{k},\omega)$ is the solution of the following Helmholtz equation

$$\left(\frac{\omega^2}{c^2}+\eta \vert\mathbf{k}\vert^{2\gamma+2}+... ...},\omega)=\frac{1}{(2\pi)^{d+1}}e^{-i \omega t_0} e^{i\mathbf{k} \mathbf{x}_s}.$$ (6)

Solving for Green's function $G(\mathbf{k},\omega)$, applying $d+1$ dimensional inverse Fourier transformation, and then integrating the kernel function with respect to $\omega$ based on Cauchy's residue theorem, we have the following time-domain attenuated Green's function

$$G_{att}(\mathbf{x},t)=\frac{c^2}{(2\pi)^{d}} \int_{\mathbb{C}^d} \frac{\mathrm{sin} (\xi_1 t)e^{-\xi_2 t}}{\xi_1} d\mathbf{k}.$$ (7)

The compensated Green's function can be obtained by reversing the absorption-related term $\tau$ in sign but leaving the other term $\eta$ unchanged:

$$G_{comp}(\mathbf{x},t)=\frac{c^2}{(2\pi)^{d}} \int_{\mathbb{C}^d} \frac{\mathrm{sin} (\xi_1 t)e^{\xi_2 t}}{\xi_1} d\mathbf{k},$$ (8)

where $G_{att}$ and $G_{comp}$ represent attenuated and compensated Green's functions, respectively. These two Green's functions lay the foundation for designing an adaptive stabilization operator. Inspired by stabilization in inverse $Q$ filtering (Irving and Knight, 2003; Wang, 2006,2002), we proposed a similar adaptive stabilization for $Q$-RTM, which can be defined as

$$\Lambda (\mathbf{k}, t)= \frac{\beta (\mathbf{k}, t)}{\beta^2 (\mathbf{k}, t)+\sigma^2} = \frac{e^{\xi_2 t}}{1+\sigma^2e^{2\xi_2 t}},$$ (9)

where the amplitude-attenuated operator $\beta(\mathbf{k}, t)=e^{-\xi_2 t}$. The final form of the proposed stabilization operator can be given by

$$s(\mathbf{k}, l \Delta t) = \left\{ \begin{array}{lll} \frac{... ...{2\xi_2 l \Delta t}}, & l=2,3, \dots, n. \\ \end{array}\right.$$ (10)