Analytical path-summation imaging of seismic diffractions

Diffraction imaging aims to emphasize small subsurface objects, such as faults, fracture swarms, and channels. Similar to classical reflection imaging, velocity analysis is crucially important for accurate diffraction imaging. Pathsummation migration provides an imaging method that produces an image of the subsurface without picking a velocity model. Previous methods of path-summation imaging involve a discrete summation of the images corresponding to all possible migration velocity distributions within a predefined integration range and thus involve a significant computational cost. We have developed a direct analytical formula for path-summation imaging based on the continuous integration of the images along the velocity dimension, which reduces the cost to that of only two fast Fourier transforms. The analytic approach also enabled automatic migration velocity extraction from diffractions using a double path-summation migration framework. Synthetic and field data examples confirm the efficiency of the proposed techniques.


INTRODUCTION
Path-integral formulation provides an efficient imaging method, which does not require subjective and time-consuming velocity picking and produces an image of the subsurface without specifying a velocity model . The concept of path-integral imaging is based on the fact that an image can be obtained from summing (stacking) of wavefields along wavefronts: the image accumulates contributions from only those events that are nearly in phase and which correspond to the true media features (Keydar, 2004;Landa, 2004). In other words, path-integral method constructs the final image by summing a set of images obtained from the range of different velocity distributions without searching for the "optimal" one. Landa et al. (2006) and Landa (2009) propose to measure the flatness of common image gathers and use it to weight images before stacking them to perform pathintegral imaging. These weights are used to enhance coherent summation of images migrated with velocities in the vicinity of the true velocity and to cancel contributions from incorrect velocities. This approach brings path-integral application for seismic Analytical path-summation imaging imaging purposes in correspondence with path-integral formulation in quantum mechanics (Feynman and Hibbs, 1965). In this paper, we follow  and employ the inherent property of diffractions' behavior under time-migration transformation -apex stationarity -for time path-summation imaging implementation for diffraction imaging. We use a term "path summation" in order to highlight that no weighting is applied before stacking the images as opposed to path-integral scheme, where weights represent assumed probability of each individual image.
A path-summation image can be calculated through a summation of a set of constant velocity images, generated by velocity continuation (VC), which accounts for both lateral and vertical shifts of the event under velocity perturbation. The VC transformation is continuous and analytical: a velocity step between constant velocity images can be in theory infinitely small (Claerbout, 1986;Fomel, 1994Fomel, , 2003aDecker and Fomel, 2014). Diffractions are sensitive to velocity perturbation (Novais et al., 2006): hyperbola flanks change their shape under VC transformation whereas their apexes remain stationary. Therefore, stacking constant velocity images can superimpose diffraction apexes constructively, cancel out hyperbola flanks and generate an image . If the constant velocity images are weighted before stacking by their corresponding velocities, they produce a double path-summation image, which can be used for velocity estimation (Schleicher and Costa, 2009;Santos et al., 2016).
The forementioned approach for path-summation migration requires a sequence of VC steps to generate a set of constant velocity images and therefore may consume significant computational resources. Moreover, produced images appear to be contaminated by artifacts compared to the images produced by picking velocities . The artifacts come from incomplete cancellation of underand overmigrated flanks of hyperbolas. We refer to these flanks as tails. Tails might overlap with useful signal and generate spurious features in the images. Note that tails are an intrinsic feature of unweighted path-summation images, as opposed to path-integral images, where the weighting describing the probability of each image and corresponding velocity is applied before summation.
In this paper, we propose a direct analytical approach to performing path-summation and double path-summation migration as opposed to stacking of multiple constantvelocity images. Computational efficiency is improved significantly because we do not require multiple velocity continuation steps for calculation of constant velocity images but instead compute the path integral in one step at the cost of only two Fast Fourier Transforms (FFTs). To improve the image quality, we propose to apply a Gaussian weighting scheme to eliminate tails in path-summation migration images. We test the effectiveness and the robustness of the proposed workflow on synthetic and field data examples.

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, 1994Fomel, , 2003a: After a substitution σ = t 2 and a 2D Fourier transform, the analytical solution for equation (1) takes the following form (Fomel, 2003b): where v is VC velocity, Ω is the Fourier dual of σ, k is the wavenumber,P 0 (Ω, k) is input stacked and not migrated data, which corresponds to v = 0, andP (Ω, 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.  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 (Larner and Beasley, 1987;Mikulich and Hale, 1992;Yilmaz et al., 2001). 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 . Path-summation migration can be performed by simply stacking constant velocity images generated by velocity continuation . 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: 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 (Ω, k) domains. For the latter case it takes the form:

TCCS-12
whereP 0 (Ω, k) corresponds to zero-offset or stacked data and does not depend on migration velocity. Therefore, it can be taken outside of the integral: We can treat the remaining integral as a filter dv, which depends on frequency Ω, 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): 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 poststack case one needs to substitute the absolute value |k| of the wavenumber vector 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 (Schleicher and Costa, 2009;Santos et al., 2016) allows for automatic migration velocity extraction. Corresponding analytical evaluation formulae are provided in Appendix B.

Attenuating tail artifacts
Artifacts in the path-summation images appear due to the absence of weights, which describe probability of the image migrated with a certain velocity and taper out Analytical path-summation imaging "unlikely" solutions Landa, 2009;Schleicher and Costa, 2009). We propose to use analytical Gaussian weighting of images being stacked in pathsummation migration, in which velocities in the vicinity of the most likely one are assigned a higher weight. An alternative approach based on the adaptive subtraction method is described in our previous work (Merzlikin and Fomel, 2015).
The following analytical weighting term can be added to the velocity continuation phase-shift:P ( and corresponding integral expression describing path-summation migration takes the following form: where filter F GP I (Ω, k, v a , v b , β, v bias ) is evaluated analytically as follows: Here, γ(Ω, k, β) = i ik 2 /16Ω + β, β is a parameter controlling contribution of each velocity to the final image and v bias is a velocity bias. By choosing a constant value of v bias for the whole section between v a and v b we can enhance the apexes of summed hyperbolas and decrease tails contribution to the final image. Since in a section, velocity varies and v bias is constant, over-or undermigrated diffraction curves may get emphasized at locations where the true velocity is higher or lower than v bias respectively. This is a possible drawback of the Gaussian weighting approach.
We analyze a response of a single diffraction hyperbola modelled with 1.5 [km/s] velocity (Figure 1a) to both path-summation and Gaussian weighted path-summation migrations. Synthetic model spectrum in Ω − k domain is shown in Figure 1b. Pathsummation migration filter magnitude as well as its phase spectrum are shown in Figures 2a and 2c. It can be interpreted as a particular form of dip filtering. Wavenumbers in the vicinity of zero appear to be preserved and either remain stationary or exhibit small phaze-shift values. This area corresponds to the diffraction hyperbola apex remaining stationary under time-migration velocity perturbation. Non-zero wave-number events have lower corresponding filter coefficients, which taper out to the spectrum boundaries. Figure 3c shows the result of path-summation migration Analytical path-summation imaging applied to diffraction hyperbola. Several lobes with non-zero wavenumber values appear to be noticeable. Corresponding phase-shift values (Figure 2c) are not zero and lead to displacement of the corresponding events on the path-summation migration image (Figure 3c). These events are associated with tail artifacts.
Gaussian weighting scheme spectrum magnitude and phase are shown in Figures 2b and 2d. The result of corresponding path-summation migration in Ω − k domain is shown in Figure 3b. Gaussian weighting scheme tapers non-stationary lobes and, therefore, suppresses the tails. Figure 3d corresponds to the Gaussian weighting scheme for path-summation migration. The tail artifacts are successfully eliminated.
Proposed expressions for direct and analytical path-integral evaluation in one step allow us to drastically reduce computational costs. Following sections illustrate the applicability and the effectiveness of the proposed techniques.

Reflection/diffraction separation
Diffraction amplitudes can be several orders of magnitude lower than those of reflections (Klem-Musatov, 1994). In order to highlight subtle diffractivity features a reflection/diffraction separation procedure has to be applied to the data. There is a variety of methods based on different separation techniques. Kanasewich and Phadke (1988); Landa and Keydar (1998) To perform reflection/diffraction separation we used plane-wave destruction filter (Fomel, 2002;Fomel et al., 2007) for both field data examples. For both examples reflections appear to be continuous and laterally consistent. This behavior aids PWD reflection's prediction along estimated local slopes followed by their subtraction from the stack.

2D field data example
We tested the proposed workflow on a vintage Gulf of Mexico dataset (Claerbout, 2005;Fomel et al., 2007). A stacked section of diffractions separated using planewave destruction filter is shown in Figure 4a. The path-summation migration image, which was calculated directly according to equation (6), is shown in Figure 4e Figure 4f. The image in comparison with the regular path-summation migration image appears to have higher resolution. Tail artifacts appear to be suppressed as well. However, some velocity bias is noticeable in the upper part of the section.
The migration velocity field acquired by division of double path-summation integral by path-summation integral and the corresponding image of diffractions are shown in Figures 4c and 4b. Velocity values in the regions without diffraction energy are associated with noise. Strong velocity anomalies at the boundaries of the section are associated with edge effects. A possible extension of the workflow would include velocity model weighting by diffraction probability at selected locations.
The image derived from double-path-summation velocity has the highest resolution. A full image with both reflections and diffractions imaged using the same velocity field is shown in Figure 4d. Discontinuities in seismic reflections align with faults imaged by path-summation diffraction imaging.

3D FIELD DATA EXAMPLE
We applied the proposed path-summation diffraction imaging framework to a 3D land seismic dataset acquired to characterize a tight-gas sand reservoir in the Cooper Basin in Western Australia. The approach has been applied to a 3D stacked volume. Detailed geological interpretation of diffraction images as well as their comparison to discontinuity-type attributes are provided in Tyiasning et al. (2016).
The input stack is shown in Figure 5a. Diffraction-and-noise volume, which corresponds to the result of PWD applied to the stacked volume, is shown in Figure 5b. Diffraction energy associated with small scale subsurface heterogeneities is significantly emphasized in comparison with the initial stack (Figure 5a).
We have applied automatic migration velocity extraction by double path summation approach to the dataset. Extracted velocity distribution is shown in Figure 8. The velocity limits used for path-summation integration correspond to 2.9 [km/s] and 3.1 [km/s]. They were chosen based on an a priori migration velocity distribution.
We used automatically extracted velocity to migrate initial ( Figure 6) and diffraction data (Figure 7). Images shown in the Figures 6 and 7 are extracted from the whole migrated volume along the target horizon -interface between coal layer and tight gas sand. Diffraction image (Figure 7) appears to have a high spatial resolution. It should be noticed that some of the most visible features in the diffraction image can be identified in the conventional image as well. For instance, the elliptical structure limited by inlines 6 − 8.5 [km] and crosslines 0.5 − 2.5 [km] is prominent in both images, but is more emphasized in the diffraction image. This observation correlates with the fact that diffractions and reflections coexist on conventional images, but reflections exhibit higher energy and may mask smaller-scale scatterers.
We also performed unweighted path-summation migration. To illustrate the differences between the images we focus on the area of the image corresponding to 6 − 12 [km] inlines and 0.5 − 4.5 [km] crosslines. Unweighted path-summation image ( Figure 9c) has lower spatial resolution in comparison to double path-summation velocity image (Figure 9a). However, the latter one appears to be noisier.
A number of features such as small faults corresponding to 9.5 − 10.5 [km] inlines and 1 − 3 [km] crosslines and the elliptical structure in the lower left corner of the figures look more focused on the double-path-summation image (Figure 9a) than on the Figure 9e. This observation can be attributed to double-path-summation migration velocity extraction directly from diffractions as opposed to PSTM velocity, which has been estimated on conventional full-waveform data containing significantly stronger reflections masking diffractions. Figure 9b shows an image acquired with double-path-summation migration velocity and a wider integration range -1.5 − 4.5 [km/s]. The image appears to be much less focused, however the pattern of the elliptical structure can still be identified. Unweighted path-summation migration image (Figure 9d) appears to be significantly blurred due to overlap between tails and useful signal. Gaussian weighting path-summation scheme allows to significantly suppress tails and highlight actual apexes of hyperbolas (Figure 9f). We used 3.1 [km/s] velocity bias with β = 10. In comparison to Figure 9c Gaussian weighting path-summation image has lower spatial resolution. This observation is connected with wider integration range used for Gaussian weighting path-summation migration in comparison to the image 9c acquired with 2.9 − 3.1 [km/s] velocity limits. When a wider integration range is used, area around the apex, within which diffraction hyperbolas are stacked coherently, becomes larger leading to the resolution loss.

DISCUSSION
Path-summation images can be contaminated by tails, which come from an incomplete cancellation of hyperbolas' flanks corresponding to the lowest and highest velocities. We propose to apply Gaussian weighting scheme for tails' elimination. In that case, no additional computational cost is required besides the direct integral evaluation itself. However, the weighting approach may introduce some velocity bias into pathsummation images.
Resolution of path-summation images might appear inferior to that of conventional workflows based on picking "one optimal" velocity. The loss of resolution comes from the summation over a set of diffraction images: coherent summation is not limited to the very apex of diffraction hyperbola but rather to the flat area around it. "Blurred" path-summation images naturally incorporate the uncertainty of velocity estimation. Therefore, they can be treated as diffraction probability volumes. Analytical evaluation allows path-summation imaging framework to be used to constrain the model space for diffraction imaging and inversion problems.
On the other hand, images built with velocity models from division of double path-summation integral by not weighted path-summation integral have a relatively high resolution and are not blurred. In the examples provided in this paper those images appear to be optimally focused. In the analytical double-path-summation technique, no picking is required.

CONCLUSIONS
We have developed an approach for direct and analytical path-summation imaging, which improves the accuracy and significantly reduces the computation cost of numerical summation. The proposed approach is based on velocity continuation, which describes continuous image transformation with perturbation of migration velocity, and in theory can have an infinitely small velocity discretization step. As a result, path-summation migration amounts to filtering in the frequency-wavenumber domain and has a total cost equivalent to the cost of only two FFTs.
Tail artifacts inherent to path-summation migration images can be efficiently re-Analytical path-summation imaging    TCCS-12 Analytical path-summation imaging continuation phase-shift for a pre-stack case: where h is a constant half-offset value, and v 0 is a velocity for the first migration run applied before cascade of VC transformations. Integral describing pre-stack pathsummation migration will take the following form: is a pre-stack path-summation migration filter, which can be evaluated analytically as:

APPENDIX B: DOUBLE PATH-SUMMATION MIGRATION
Double path-integral formulation was proposed by Costa (2009), Santos et al. (2016). In unweighted path-summation migration framework it can be described by the following equation: Velocity distribution can be acquired through a smooth division of the double pathsummation integral by the path-summation integral: v(t, x) = I DP I (t, x, v a , v b )/I P I (t, x, v a , v b ) .
(A-5) Analytical path-summation imaging We propose to evaluate the double path-summation integral analytically as well: where the filter F DP I (Ω, k, v a , v b ) has the following analytical expression: To extract migration velocity as described by the equation (A-5) we need to evaluate both double path-summation integral and path-summation integral. These evaluations can be done analytically as described by equations 5 and A-6. A smooth division procedure (Fomel, 2007) is applied in equation A-5 to acquire a regularized time migration velocity estimate, which then can be applied to image the subsurface using any available time imaging algorithm.