Diffraction imaging aims to emphasize small-scale subsurface heterogeneities such as faults, pinch-outs, fracture swarms, channels and other small features and can play an important role in seismic reservoir characterization (Landa, 2012). Reflections and diffractions coexist on seismic records, however the latter typically have significantly lower energy than the former (Klem-Musatov, 1994). Moreover, most of the conventional seismic data processing is tuned to imaging and enhancing reflected waves, whereas diffractions on seismic images are being partially suppressed (Kozlov et al., 2004). Therefore, reflection and diffraction separation procedure is a key step in diffraction imaging workflows (Fomel et al., 2007; Khaidukov et al., 2004; Harlan et al., 1984).

Diffraction imaging methods can be classified based on the separation technique being employed. Methods based on optimal stacking of diffracted energy and suppression of reflections are described by Berkovitch et al. (2009); Tsingas et al. (2011); Rad et al. (2014); Dell and Gajewski (2011); Landa and Keydar (1998); Kanasewich and Phadke (1988). Wavefield separation methods aim to decompose conventional full-wavefield seismic records into different components representing reflections and diffractions (Papziner and Nick, 1998; Fomel et al., 2007; Klokov and Fomel, 2012b; Tyiasning et al., 2016; Reshef and Landa, 2009; Klokov et al., 2017; Taner et al., 2006; Merzlikin et al., 2017b; Harlan et al., 1984; Merzlikin et al., 2017a). Koren and Ravve (2011); Moser and Howard (2008); Kozlov et al. (2004); Klokov and Fomel (2013); Popovici et al. (2015) modify the Kirchhoff migration kernel to eliminate specular energy coming from the first Fresnel zone and image diffractions only.

Seismic waves diffract on heterogeneities with a size of a wavelength, whereas reflected energy is associated with larger objects in the subsurface corresponding to the size of the first Fresnel zone (Moser and Howard, 2008). Therefore, separating diffractions and imaging them independently can produce an image with a higher spatial resolution than that of the conventional migration. On the other hand, this resolution is limited by using conventional migration operators, which are tailored for reflection processing, in diffraction imaging workflows. To improve diffraction imaging resolution, Gelius et al. (2013) propose to use singular value decomposition of the data in windows along the Kirchhoff diffraction-stacking aperture. This approach significantly improves resolution of diffraction images but requires a significant computational effort. Huang and Schuster (2014) propose to use resonant multiples to improve resolution of diffraction images. Their method requires nonlinear least-squares reverse time migration (RTM) or full wave inversion (FWI) updates of the background reflectivity after each iteration to account for the contribution of multiples.

Robust diffraction separation from noise and reflections while accounting for the crosstalk between them is of utmost importance for a resolution increase, which cannot be accomplished without enhancing subtle features with magnitudes close to the noise level. The problem of diffraction and noise separation was first addressed by Harlan et al. (1984). Klokov and Fomel (2012a); Klokov et al. (2010b,a); Burnett et al. (2015) enhance diffractions' signal-to-noise ratio in the dip-angle gathers' domain.

Least-squares migration (Nemeth et al., 1999; Ronen and Liner, 2000) can produce higher resolution than conventional migration algorithms and suppress artifacts intrinsic to the latter. Effectively an operator for diffraction imaging can be constructed in an iterative fashion using least-squares migration framework. A priori knowledge about the signature of a scatterer's response to a seismic wavefield can be incorporated into the inversion in the form of regularization and allow for diffraction separation from noise and reflections. Merzlikin and Fomel (2016) perform least-squares migration chained with plane wave destruction (PWD) (Fomel et al., 2007; Fomel, 2002) and path-summation integral filtering while enforcing sparsity in a diffraction model. Yu et al. (2016) utilize common-offset Kirchhoff least-squares migration with a sparse model regularization to emphasize diffractions. Yu et al. (2017a) extract diffractions based on plane wave destruction and dictionary learning for sparse representation. Yu et al. (2017b) use two separate modeling operators for diffractions and reflections and impose a sparsity constraint on diffractions. Decker et al. (2017) denoise diffractions by weighting the model with semblance estimated in dip-angle gather domain (DAG) in least-squares migration.

In this paper, we propose a workflow to decompose the input seismic data into reflections, diffractions and noise. To restore both reflections and diffractions simultaneously we extend the least-squares migration workflow presented by Merzlikin and Fomel (2016), in which the inverted forward modeling operator is the chain of Kirchhoff modeling, plane wave destruction and path-summation integral filter operators. The proposed chain of operators allows for proper incorporation of diffraction component into the inversion as opposed to the conventional least-squares migration dominated by stronger reflections (Merzlikin and Fomel, 2016). Incorporation of reflection model space into the inversion reduces the crosstalk between reflections and diffractions. Both reflections and diffractions are modeled using the same forward modeling operator. Separation into the components is accomplished by shaping regularization (Fomel, 2007b,2008) and by local signal-and-noise orthogonalization (Chen and Fomel, 2015). Following Merzlikin and Fomel (2016) diffractions are penalized by sparsity constraints allowing to distinguish them from noise. Reflections are forced to be smooth along the dominant local slopes in the image domain. Local signal-and-noise orthogonalization is applied additionally to further minimize the leakage of reflections into the diffraction image domain: local similarity (Fomel, 2007a) is computed between two models and events following reflected energy pattern are removed from the diffraction image. Although we develop and illustrate the efficiency of the method for the zero-offset time-migration case, all the developments can be extended to depth and pre-stack domains.

We start with a theory section describing all the components of the workflow. We describe the workflow and then demonstrate its efficiency on synthetic and data examples and discuss its comparative advantages and disadvantages.