*2.1. Theory and Algorithm*

In an isotropic acoustic medium, the wave equation is as follows [8,10]:

$$\begin{cases} \left(\frac{1}{v\_0^2} \frac{\partial^2}{\partial t^2} - \rho\_0 \nabla \frac{1}{\rho\_0} \cdot \nabla \right) p\_0(\mathbf{x}; t; \mathbf{x}\_\mathbf{s}) = \delta(\mathbf{x} - \mathbf{x}\_\mathbf{s}) \delta(t) \\\ d(\mathbf{x}\_\mathbf{r}; t; \mathbf{x}\_\mathbf{s}) = p\_0(\mathbf{x} = \mathbf{x}\_\mathbf{r}; t; \mathbf{x}\_\mathbf{s}) \end{cases} \tag{1}$$

where *v*<sup>0</sup> and *ρ*<sup>0</sup> are velocity and density, respectively; *p*0(**x**; *t*; **x***s*) is the pressure wavefield at any location of **x** due to a source located at **x***s*; and *d*(**x***r*; **x***s*; *t*) is the recorded data, i.e., the measured value of the pressure wavefield at a receiver position **x** = **x***r*.

For a second medium with small perturbations in velocity (*δv*) and density (*δρ*) compared to the previous medium, the velocity and density are *v*<sup>0</sup> + *δv* and *ρ*<sup>0</sup> + *δρ*, respectively. The pressure wavefield in the perturbed medium *p*<sup>0</sup> + *δp* satisfied the same acoustic wave equation as follows:

$$\begin{cases} \left(\frac{1}{\left(\mathbf{r}\_0 + \delta\mathbf{v}\right)^2 \delta t^2} - \left(\rho\_0 + \delta\rho\right) \nabla \frac{1}{\rho\_0 + \delta\rho} \cdot \nabla\right) \left(p\_0 + \delta p\right) = \delta(\mathbf{x} - \mathbf{x}\_\delta)\delta(t) \\ \qquad \delta d\left(\mathbf{x}\_\mathbf{r}; t; \mathbf{x}\_\delta\right) = \delta p\left(\mathbf{x} = \mathbf{x}\_\mathbf{r}; t; \mathbf{x}\_\delta\right) \end{cases} \tag{2}$$

where *δp* is the wavefield perturbation, and the observed data perturbation is represented by *δd*.

Based on the Born approximation, the following equation for *δp* can be derived by subtracting Equation (1) from Equation (2):

$$
\delta \left( \frac{1}{v\_0^2} \frac{\partial^2}{\partial t^2} - \rho\_0 \nabla \frac{1}{\rho\_0} \cdot \nabla \right) \delta p(\mathbf{x}; t; \mathbf{x}\_s) \approx \left( \frac{2 \delta v}{v\_0^3} \frac{\partial^2}{\partial t^2} - \left( \nabla \frac{\delta \rho}{\rho\_0} \right) \cdot \nabla \right) p\_0(\mathbf{x}; t; \mathbf{x}\_s) \tag{3}
$$

By following Zhang et al. [10] and using the asymptotic approximation, we obtained the ray-based relationship between *δd*, *δv* and *δρ*:

$$\delta d(\mathbf{x}\_{\mathbf{r}};\omega;\mathbf{x}\_{\mathbf{s}}) = -\int \frac{2\omega^2}{\mathbf{v}\_0(\mathbf{x})} \left(\frac{\delta v(\mathbf{x})}{v\_0(\mathbf{x})} + \cos^2\theta \frac{\delta\rho(\mathbf{x})}{\rho\_0(\mathbf{x})}\right) \times A(\mathbf{x}\_{\mathbf{s}};\mathbf{x}\_{\mathbf{r}};\mathbf{x}) e^{i\omega \mathcal{T}(\mathbf{x}\_{\mathbf{s}};\mathbf{x}\_{\mathbf{r}})} d\mathbf{x} \tag{4}$$

where *T*(**x***s*; **x***r*; **x**) = *τ*(**x**; **x***s*) + *τ*(**x***r*; **x**) and *A*(**x***s*; **x***r*; **x**) = *A*(**x**; **x***s*)*A*(**x***r*; **x**) represent the travertine summation and the amplitude product of the Green's function from the source location **x***<sup>s</sup>* to the image point **x** and reflected back to the receiver location **x***r*, respectively, and *θ* is the subsurface reflection angle. Equation (4) represents the forward modeling

formulation under a high-frequency assumption. The detailed derivations of Equations (3) and (4) can be found in [10].

Following a similar method to [5] and [8], we can invert Equation (4) for the composite model parameter perturbation, *<sup>δ</sup>v*(**x**) *<sup>v</sup>*0(**x**) <sup>+</sup> cos2*<sup>θ</sup> δρ*(**x**) *<sup>ρ</sup>*0(**x**). In 2D, the composite model parameter perturbation can be described in terms of the perturbed wavefield as follows (detailed derivation is given in Appendix A):

$$\begin{array}{c} \sin^2\theta \frac{\delta\boldsymbol{v}}{\boldsymbol{v}\_0} + \cos^2\theta \frac{\delta(\boldsymbol{\rho}\boldsymbol{v})}{\rho\_0 \boldsymbol{v}\_0} = -\iint \iint \frac{32\cos^2\theta'}{|\boldsymbol{\omega}|} \frac{\cos\theta\_\mathbf{r}}{\boldsymbol{v}(\mathbf{x}\_\theta)} \frac{\cos\theta\_\mathbf{s}}{\boldsymbol{v}(\mathbf{x}\_\theta)} A(\mathbf{x};\mathbf{x}\_\theta) A(\mathbf{x}\_\mathbf{r};\mathbf{x}) \\ \times (\theta'-\theta)\delta d(\mathbf{x}\_\mathbf{r};\omega;\mathbf{x}\_\mathbf{s}) e^{-i\omega\boldsymbol{v}\boldsymbol{T}(\mathbf{x}\_\theta;\mathbf{x}\_\mathbf{r})} d\mathbf{x}\_\mathbf{r} d\mathbf{x}\_\mathbf{s} d\omega d\theta' \end{array} \tag{5}$$

where β<sup>s</sup> and β<sup>r</sup> represent the takeoff angles at the source and receivers, respectively.

Within the framework of amplitude-preserving RTM, the asymptotic forms of *pF*(**x**; *ω*; **x***s*) and *pB*(**x**; **x***s*; *ω*) for 2D acoustic case are given as follows [10,12]:

$$p\_F^\*(\mathbf{x}; \omega; \mathbf{x}\_s) = -2\frac{\cos \beta\_s}{v(\mathbf{x}\_s)} \sqrt{|\omega|} A(\mathbf{x}; \mathbf{x}\_s) e^{-i\omega \tau(\mathbf{x}; \mathbf{x}\_s) - i(\frac{\pi}{4}) \text{sgn}(\omega)} \tag{6}$$

and

$$p\_{\mathcal{B}}(\mathbf{x};\omega;\mathbf{x}\_{\mathfrak{k}}) = 2 \int \frac{\cos \beta\_{\mathfrak{r}}}{v(\mathbf{x}\_{\mathfrak{r}})} \sqrt{|\omega|} A(\mathbf{x}\_{\mathfrak{r}};\mathbf{x}) e^{-i\omega \tau(\mathbf{x}\_{\mathfrak{r}};\mathbf{x}) - i\left(\frac{\mathfrak{a}}{\mathfrak{r}}\right) \text{sgn}(\omega)} \delta d(\mathbf{x}\_{\mathfrak{r}};\omega;\mathbf{x}\_{\mathfrak{k}}) d\mathbf{x}\_{\mathfrak{r}} \tag{7}$$

Substituting Equations (6) and (7) for the terms on the right-hand side of Equation (5), we can obtain:

$$\sin^2\theta \frac{\delta\upsilon}{\upsilon\_0} + \cos^2\theta \frac{\delta(\rho\upsilon)}{\rho\_0\upsilon\_0} = \iiint \frac{8\cos^2\theta'}{\omega^2} \delta\left(\theta'-\theta\right) p\_F^\*(\mathbf{x};\omega;\mathbf{x}\_s) p\_B(\mathbf{x};\omega;\mathbf{x}\_s) \mathrm{d}\omega \mathrm{d}\theta' \mathrm{d}\mathbf{x}\_s \tag{8}$$

The left-hand side of Equation (8) is the composite form of relative velocity perturbation and impedance perturbation. The right-hand side has the same form as RTM. This shows that the composite parameter on the left-hand side can be estimated using the RTM framework. In order to invert each individual parameter, Zhang et al. [10] pointed out that the near-angle stacked image can output impedance perturbation, <sup>δ</sup>(*ρv*) *<sup>ρ</sup>*0*v*<sup>0</sup> , while the far-angle stacked image can be used to estimate the velocity perturbation, <sup>δ</sup>*<sup>v</sup> <sup>v</sup>*<sup>0</sup> . They separated the effects of impedance and velocity by first generating RTM angle-domain common-image gathers, and then using stacked images within different angle sections for velocity and impedance estimations. However, it could be computationally expensive to compute RTM angle gathers, especially in 3D.

Note that the right-hand side of Equation (8) shares a similar form as the true amplitude imaging principle proposed by Kiyashchenko et al. [13], where the cos2θ term is approximated by the ray theoretical slowness vectors in the Fourier domain and computed by applying the time derivatives and spatial gradients of wavefields. The time derivatives and spatial gradients of wavefields were also utilized in the inverse scattering imaging condition proposed by Whitmore and Crawley [11]. Next, we investigated the inverse scattering imaging condition, and further proposed a modified inverse scattering imaging condition to output the near-angle stacked images without computing the angles. Our method can reduce the computational costs compared with the method using angle gathers, and still utilize the imaging capabilities of RTM for the relative impedance estimation.

The inverse scattering imaging condition proposed by Whitmore and Crawley [11]:

$$I(\mathbf{x}; \mathbf{x}\_{\sf s}) = I\_{\sf \nabla}(\mathbf{x}; \mathbf{x}\_{\sf s}) + B(\mathbf{x}; \mathbf{x}\_{\sf s}) I\_{\sf d}(\mathbf{x}; \mathbf{x}\_{\sf s}) \tag{9a}$$

where

$$I\_{\nabla}(\mathbf{x}; \mathbf{x}\_{\sf s}) = \int \nabla p\_{\sf F}(\mathbf{x}; t; \mathbf{x}\_{\sf s}) \cdot \nabla p\_{\sf B}(\mathbf{x}; t; \mathbf{x}\_{\sf s}) dt \tag{9b}$$

and

$$I\_{dt}(\mathbf{x}; \mathbf{x}\_s) = -\int \frac{1}{\upsilon^2(\mathbf{x})} \frac{\partial p\_F(\mathbf{x}; t; \mathbf{x}\_s)}{\partial t} \frac{\partial p\_B(\mathbf{x}; t; \mathbf{x}\_s)}{\partial t} dt \tag{9c}$$

and where *B*(**x**; **x***s*) is a weighting function to attenuate backscattered energy. Note that here we used the backpropagated wavefields *pB*(**x**; *t*; **x***s*) in terms of wavefield propagation time *t* (instead of *T* − *t*, where *T* is the maximum extrapolation time), and used the relationship *∂pB*(**x**;*T*−*t*;**x***s*) *<sup>∂</sup><sup>t</sup>* <sup>=</sup> <sup>−</sup>*∂pB*(**x**;*t*;**x***s*) *<sup>∂</sup><sup>t</sup>* .

The relationship between the product of the time derivatives of *pF* and *pB*, and the products of their spatial gradients is given by the following equation [11,14,15]:

$$-\frac{\partial p\_F}{\partial t}\frac{\partial p\_B}{\partial t}\frac{\cos(2\theta)A(\mathbf{x})}{v^2(\mathbf{x})} = \nabla p\_F \cdot \nabla p\_B(\mathbf{x}; t; \mathbf{x}\_\mathbf{s})\tag{10}$$

where *A*(**x**) is the parameter to compensate for far field approximation. We can ignore it by assuming far field approximation and obtain the following equations:

$$\nabla p\_F(\mathbf{x};t;\mathbf{x}\_s) \cdot \nabla p\_B(\mathbf{x};t;\mathbf{x}\_s) \approx -\frac{\cos(2\theta)}{v^2(\mathbf{x})} \frac{\partial p\_F(\mathbf{x};t;\mathbf{x}\_s)}{\partial t} \frac{\partial p\_B(\mathbf{x};t;\mathbf{x}\_s)}{\partial t} \tag{11a}$$

$$-\frac{1}{v^{2}(\mathbf{x})}\frac{\partial p\_{\mathbf{F}}(\mathbf{x};t;\mathbf{x}\_{\star})}{\partial t}\frac{\partial p\_{\mathbf{B}}(\mathbf{x};t;\mathbf{x}\_{\star})}{\partial t} + \nabla p\_{\mathbf{F}}(\mathbf{x};t;\mathbf{x}\_{\star}) \cdot \nabla p\_{\mathbf{B}}(\mathbf{x};t;\mathbf{x}\_{\star}) \approx -\frac{2\cos^{2}\theta}{v^{2}(\mathbf{x})}\frac{\partial p\_{\mathbf{F}}(\mathbf{x};t;\mathbf{x}\_{\star})}{\partial t}\frac{\partial p\_{\mathbf{B}}(\mathbf{x};t;\mathbf{x}\_{\star})}{\partial t} \tag{11b}$$

$$-\frac{1}{\upsilon^{2}(\mathbf{x})}\frac{\partial p\_{\mathbf{F}}(\mathbf{x};t;\mathbf{x}\_{s})}{\partial t}\frac{\partial p\_{\mathbf{B}}(\mathbf{x};t;\mathbf{x}\_{s})}{\partial t} - \nabla p\_{\mathbf{F}}(\mathbf{x};t;\mathbf{x}\_{s}) \cdot \nabla p\_{\mathbf{B}}(\mathbf{x};t;\mathbf{x}\_{s}) \approx -\frac{2\sin^{2}\theta}{\upsilon^{2}(\mathbf{x})}\frac{\partial p\_{\mathbf{F}}(\mathbf{x};t;\mathbf{x}\_{s})}{\partial t}\frac{\partial p\_{\mathbf{B}}(\mathbf{x};t;\mathbf{x}\_{s})}{\partial t} \tag{11c}$$

Choosing B(**x**) = 1, Equation (9a) can be expressed in the frequency domain as:

$$\begin{split} \mathcal{I}\_{\nabla}(\mathbf{x};\mathbf{x}\_{s}) &+ I\_{dt}(\mathbf{x};\mathbf{x}\_{s}) \\ = \int \left( -\frac{1}{\upsilon^{2}(\mathbf{x})} \frac{\partial p\_{F}(\mathbf{x};\mathbf{x}\_{s})}{\partial t} \frac{\partial p\_{B}(\mathbf{x};\mathbf{x}\_{s})}{\partial t} + \nabla p\_{F}(\mathbf{x};t;\mathbf{x}\_{s}) \cdot \nabla p\_{B}(\mathbf{x};t;\mathbf{x}\_{s}) \right) dt \\ \approx & - \int \frac{2\cos^{2}\theta}{\upsilon^{2}(\mathbf{x})} \frac{\partial p\_{F}(\mathbf{x};\mathbf{x}\_{s})}{\partial t} \frac{\partial p\_{B}(\mathbf{x};t;\mathbf{x}\_{s})}{\partial t} dt \\ = & - \int \frac{2\omega^{2}\cos^{2}\theta}{\upsilon^{2}(\mathbf{x})} p\_{F}^{\*}(\mathbf{x};\omega;\mathbf{x}\_{s}) p\_{B}(\mathbf{x};\omega;\mathbf{x}\_{s}) d\omega \end{split} \tag{12}$$

Equation (12) is similar to the energy norm imaging condition proposed by Rocha et al. [16], which was used to attenuate reflections with opening angles close to 180◦. Comparing Equations (8) and (12), we can see that the composite parameter of the relative impedance and velocity perturbation (sin2*θ <sup>δ</sup><sup>v</sup> <sup>v</sup>*<sup>0</sup> <sup>+</sup> cos2*<sup>θ</sup> <sup>δ</sup>*(*ρv*) *<sup>ρ</sup>*0*v*<sup>0</sup> ) can be computed using the inverse scattering imaging condition (Equation (12)) after the proper preprocessing of source and receiver wavefields. To estimate the impedance perturbation, we must further separate the effects of impedance and velocity. Inspired by Rocha et al. [16], who proposed an exponential weighting function to select far-angle reflections for the purpose of tomographic inversion, we applied the following weighting function for Equation (12) to select near-angle reflections [15]:

$$I\_{\text{small angle}} = \int w \left( -\frac{1}{v^2} \frac{\partial p\_F}{\partial t} \frac{\partial p\_B}{\partial t} + \nabla p\_F \nabla p\_B \right) dt,\\ \text{with } w = e^{-\frac{\kappa \left( -\frac{1}{p^2} \frac{\partial p\_F}{\partial t} \frac{\partial p\_B}{\partial t} - \nabla p\_F \nabla p\_B \right)}{-\frac{2}{v^2} \frac{\partial p\_F}{\partial t} \frac{\partial p\_B}{\partial t}} \,, \ \alpha > 1 \tag{13}$$

where the weighting term *w* is approximate to *e*−*αsin*2*<sup>θ</sup>* following Equation (11c). This weighting function equals to 1 when *θ* = 0◦, and rapidly approaches 0 when *θ* increases, which is designed to select small reflection angles. More details about this weighting function and the choice of *α* are presented in the later section. Equation (8) shows that small-angle reflection produces an estimate of the impedance perturbation. Therefore, from Equations (8) and (13), we derived our proposed imaging condition for the relative impedance estimation from RTM:

$$\frac{\delta(\rho v)}{\rho\_0 v\_0} = -\iint 4v^2 w \left( -\frac{1}{v^2} \frac{\partial p\_F'}{\partial t} \frac{\partial p\_B'}{\partial t} + \nabla p\_F' \cdot \nabla p\_B' \right) dt d\mathbf{x}\_6, \text{ with } w = c \tag{14}$$

where *p <sup>F</sup>* and *<sup>p</sup> <sup>B</sup>* are forward and backward wavefields, respectively, scaled by 1/*ω*<sup>2</sup> in the frequency domain.

Figure 1 shows the workflow of our proposed method of estimating the relative impedance perturbation from RTM using the modified inverse scattering imaging condition. Compared with conventional RTM, the only difference is that we replaced the conventional cross-correlation imaging condition with our proposed imaging condition in Equation (14).

**Figure 1.** Workflow of the proposed method.

### *2.2. Comparison of Imaging Conditions*

To illustrate how the proposed imaging condition in Equation (13) attenuates large angle reflections, we performed simple experiments on the prestack impulse response. The synthetic data were generated by a two-reflector layered model. Figure 2a shows the impulse response using *Idt* for a single trace at an offset of 4000 m, from which we can see the image of two reflection events, as well as the backscattered events. The imaging condition *I*dt + *I*<sup>∇</sup> attenuated backscattered events, as shown in Figure 2b. By using the proposed imaging condition in Equation (13), Figure 2c only preserved the small-angle reflections. Figure 2d–f show the corresponding imaging comparison for a single trace at an offset of 0 m (one source and one receiver, both at 0 m on the surface). Using our proposed small-angle imaging condition, Figure 2f preserved all the zero-angle reflections available in Figure 2d, while attenuating the backscattered energy, which is mostly notable at location 0 m.

Figure 3 compares the imaging results of the single trace at zero offset using Equations (13) and (14). From the comparison, we can see that the factor 1/*ω*<sup>2</sup> in Equation (14) boosts the low-frequency components in Figure 3b.

**Figure 2.** Imaging results of a single trace (offset = 4000 m) showing (**a**) *I*dt; (**b**) *I*dt + *I*∇; and (**c**) *I*small angle; Imaging results of a single trace (offset = 0 m) showing (**d**) *I*dt; (**e**) *I*dt + *I*∇; and (**f**) *I*small angle.

**Figure 3.** Comparison of imaging results of a single trace (offset = 0 m) (**a**) using Equation (13) and (**b**) using Equation (14).
