Analytical path-summation imaging of seismic diffractions

Method

We deduce analytical formulas for path-summation migration using the velocity continuation concept for migration (Fomel, 2003a). Velocity continuation (VC) describes vertical and lateral shifts of time-migrated events under the change of migration velocity. It is a continuous process, which, in the zero-offset isotropic case, can be described by the following partial differential equation (Claerbout, 1986; Fomel, 2003a,1994):

$$\frac{\partial^2 P}{\partial t\,\partial v} + \frac{v\,t}{2}\,\frac{\partial^2 P}{\partial x^2}=0\;.$$ (1)

After a substitution $\sigma=t^2$ and a 2D Fourier transform, the analytical solution for equation (1) takes the following form (Fomel, 2003b):

$$\tilde{P}(\Omega,k,v) = \hat{P}_0 (\Omega,k)\,e^{-\frac{i k^2 v^2 }{16\Omega}}\;,$$ (2)

where $v$ is VC velocity, $\Omega$ is the Fourier dual of $\sigma$ , $k$ is the wavenumber, $\hat{P}_0 (\Omega,k)$ is input stacked and not migrated data, which corresponds to $v=0$ , and $\tilde{P}(\Omega,k,v)$ is a constant velocity image. Thus, an input image gets transformed to a specified velocity by a simple phase shift in the Fourier domain. Burnett and Fomel (2011) extend the formulation to the 3D azimuthally anisotropic case. A sequence of phase shifts corresponding to a predefined velocity range with a certain step can be calculated to generate a set of constant velocity images (Mikulich and Hale, 1992; Yilmaz et al., 2001; Larner and Beasley, 1987). The final image is built by picking the parts of constant velocity images with the highest focusing measure and sewing them together (Fomel et al., 2007). In other words, for each $(t,x)$ location we browse through constant velocity images and pick the velocity corresponding to the highest focus of diffraction events. However, generating focusing attributes can be computationally expensive, and picking is subjective.

Path-integral formulation aims to provide velocity-model-independent images of the subsurface (Landa et al., 2006). Path-summation migration can be performed by simply stacking constant velocity images generated by velocity continuation (Burnett et al., 2011). Theoretically, the velocity step between images generated by velocity continuation can be infinitely small, and path-summation migration in this case will correspond to the following integral over the velocity dimension:

$$I_{PI}(t,x,v_a,v_b) = \int_{v_a}^{v_b} P(t,x,v)\,dv\;,$$ (3)

where $P(t,x,v)$ is a constant velocity image. The resultant image is not associated with a particular velocity model since it stacks all the constant velocity images generated by VC.

Note that stacking of images can be performed both in $(t,x)$ and $(\Omega,k)$ domains. For the latter case it takes the form:

$$\hat{I_{PI}}(\Omega,k,v_a,v_b) = \int_{v_a}^{v_b} \hat{P}_0(\Omega,k)\ e^{-\frac{i k^2 v^2 }{16\Omega}}dv\;,$$ (4)

where $\hat{P}_0 (\Omega,k)$ corresponds to zero-offset or stacked data and does not depend on migration velocity. Therefore, it can be taken outside of the integral:

$$\hat{I_{PI}}(\Omega,k,v_a,v_b) = \hat{P}_0(\Omega,k)\ \int_{v_a}^{v_b} e^{-\frac{ i k^2 v^2 }{16\Omega}}dv \;.$$ (5)

We can treat the remaining integral as a filter $F_{PI}(\Omega,k,v_a,v_b) = \int_{v_a}^{v_b} e^{-\frac{ i k^2 v^2 }{16\Omega}}dv$ , which depends on frequency $\Omega$ , wavenumber $k$ and velocity integration range limits $v_a$ and $v_b$ . Moreover, the filter has an analytical expression: the integral of the complex exponential function is proportional to the imaginary error function (erfi):

$$F_{PI}(\Omega,k,v_a,v_b) = e^{i\frac{5\pi}{4}}\ \frac{2\sqrt{\Omega \pi}}{k}\ \big(e^{i\frac{3\pi}{4}}\ \frac{k\ v}{4\sqrt{\Omega}}\big) \bigg\vert _{v_a}^{v_b}.$$ (6)

Thus, path-summation migration according to the equation (5) amounts to filtering the data with the filter given by equation (6). Inverse Fourier transform and time-stretch removal produce the final image in the $(t,x)$ domain. The cost of the computation corresponds to the cost of two fast Fourier transforms needed to transform the data to the frequency-wavenumber domain and back to time-space domain. The filtered image can be treated as if it was calculated by stacking of constant velocity images with an infinitely small velocity discretization step.

As an illustration, a path-summation image corresponding to a point scatterer with 1.5 km/s velocity (Figure 1a), which has been calculated by application of equations (4) and (6), is shown in Figure 3c. The path-summation migration image contains artifacts because of incomplete cancellation of two hyperbolas, which are still prominent in the image: one corresponding to the lowest constant velocity ($v_a$ ) image (an undermigrated tail), and the other corresponding to the highest velocity ($v_b$ ) image (overmigrated tail). Those two tails are not cancelled out during stacking of constant velocity images whereas hyperbolas' flanks within the VC velocity range sum destructively. When path-summation migration is applied to real diffraction images the tails may interfere with useful signal and create artifacts.

To extend analytical path-summation migration formulae to 3D isotropic post-stack case one needs to substitute the absolute value $\left\vert \mathbf{k} \right\vert$ of the wavenumber vector $\mathbf{k} = (k_x,k_y)$ instead of a one-dimensional wavenumber $k$ in the expressions given above. Pre-stack path-summation migration is discussed in Appendix A. Double path-summation migration (Santos et al., 2016; Schleicher and Costa, 2009) allows for automatic migration velocity extraction. Corresponding analytical evaluation formulae are provided in Appendix B.

Analytical path-summation imaging of seismic diffractions