## ABSTRACT

Diffraction imaging aims to emphasize small-scale subsurface heterogeneities, such as faults, pinch-outs, fracture swarms, channels, etc. and can help seismic reservoir characterization. The key step in diffraction imaging workflows is based on the separation procedure suppressing higher energy reflections and emphasizing diffractions, after which diffractions can be imaged independently. Separation results often contain crosstalk between reflections and diffractions and are prone to noise. We have developed an inversion scheme to reduce the crosstalk and denoise diffractions. The scheme decomposes an input full wavefield into three components: reflections, diffractions, and noise. We construct the inverted forward modeling operator as the chain of three operators: Kirchhoff modeling, plane-wave destruction, and path-summation integral filter. Reflections and diffractions have the same modeling operator. Separation of the components is done by shaping regularization. We impose sparsity constraints to extract diffractions, enforce smoothing along dominant local event slopes to restore reflections, and suppress the crosstalk between the components by local signal-and-noise orthogonalization. Synthetic- and field-data examples confirm the effectiveness of the proposed method.

## INTRODUCTION

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, the reflection and diffraction separation procedure is a key step in diffraction imaging workflows (Harlan et al., 1984; Khaidukov et al., 2004; Fomel et al., 2007).

Diffraction imaging methods can be classified based on the separation technique being used. Methods based on optimal stacking of diffracted energy and suppression of reflections are described by Kanasewich and Phadke (1988), Landa and Keydar (1998), Berkovitch et al. (2009), Dell and Gajewski (2011), Tsingas et al. (2011), and Rad et al. (2014). Wavefield-separation methods aim to decompose conventional full-wavefield seismic records into different components representing reflections and diffractions (Harlan et al., 1984; Papziner and Nick, 1998; Taner et al., 2006; Fomel et al., 2007; Reshef and Landa, 2009; Klokov and Fomel, 2012b; Tyiasning et al., 2016; Klokov et al., 2017; Merzlikin et al., 2017a, 2017b). Kozlov et al. (2004), Moser and Howard (2008), Koren and Ravve (2011), Klokov and Fomel (2013), and 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 the 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 significant computational effort. Huang and Schuster (2014) propose to use resonant multiples to improve the resolution of diffraction images. Their method requires nonlinear least-squares reverse time migration or full wave inversion 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 et al. (2010a, 2010b), Klokov and Fomel (2012a), and Burnett et al. (2015) enhance signal-to-noise ratio of the diffractions in the dip-angle gather (DAG) 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 a 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 can allow for diffraction separation from noise and reflections. Merzlikin and Fomel (2016) perform least-squares migration chained with plane-wave destruction (PWD) (Fomel, 2002; Fomel et al., 2007) and path-summation integral filtering while enforcing sparsity in a diffraction model. Yu et al. (2016) use common-offset Kirchhoff least-squares migration with a sparse model regularization to emphasize diffractions. Yu et al. (2017a) extract diffractions based on PWD 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 the DAG domain in least-squares migration.

In this paper, we propose a workflow to decompose the input seismic data into reflections, diffractions and noise. To restore 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, PWD, 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. 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 us 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 of the developments can be extended to the depth and prestack 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.

## METHOD

We start with a brief review of the path-summation imaging framework. Then, we give a formulation of the least-squares inverse problem for diffraction imaging. We then proceed with a description of an optimization strategy and introduce a regularization framework.

### Analytical path-summation migration

### Objective function

Reflections restored by equation 10 ($mrw$) are subtracted from the full-wavefield data $d$ to allow for diffraction-only restoration $mdw$. No interference is expected between reflections $mrw$ and diffractions $mdw$ provided accurate starting models from the previous inversions (equations 6 and 8).

Instead of implementing $RS\sigma ,T\lambda (mr,md)$, $RS\sigma ,T\lambda ,Omr(mr,md)$, $RS\sigma (mrw)$, and $RT\lambda (mdw)$ as penalty terms for reflections and diffractions correspondingly, the problem is reformulated in the shaping framework (see the “Regularization” subsection for details).

### Optimization strategy

It is important to highlight that the same operator is used to update the images of reflections $mr$ and diffractions $md$. The implication of this strategy is that without regularization, the same updates are attributed to both models. For simplicity, consider conjugate direction $sn$ at the zero iteration $n=0$. It has the form of $sn=0=[LTDdataTPTr,LTDdataTPTr]T$ (PWD ($Ddata$), and the path-summation integral $P$ is disabled for objective functions in equations 8 and 9 correspondingly), where $r$ corresponds to the residual between the observed data $d$ and the data modeled from the initial guess model $[mrinit,mdinit]T$ (initialized by zeros for the first inversion [equation 6] and by the output of the first inversion [equation 6] for optimization of the objective function in equation 8). The conjugate direction $sn=0$ is equal to the negative gradient of the objective function $\u2212\u2207J(mrinit,mdinit)$. It is obvious that the residual $r$ in the conjugate direction $sn=0$ is mapped to the reflection and diffraction image using the same operator — the adjoint of the “chain.” The same is true for all of the iterations. In this case, previous directions $sn\u22121=[srn\u22121,sdn\u22121]T$, which are also the same for both models, participate in the update. To perform reflection and diffraction separation, regularization is required.

### Regularization

As soon as PWD is removed from the misfit term in the second inversion cascade (equation 8), the residual gets dominated by reflections. Thus, before disabling PWD, we need to make sure that most of the diffracted energy has been extracted by the first inversion (equation 6). Because of the high energy of the reflections, the sparsity constraints alone can be insufficient to remove the reflected component leaked to the diffraction image after disabling PWD. Reflection removal will require higher threshold levels leading to diffraction suppression.

To address this issue, a signal and noise orthogonalization workflow is used (Chen and Fomel, 2015). We follow the approach and measure the similarity between the update to diffractivity $sdn$ and reflectivity obtained on the previous external iteration $mrk\u22121$. Events in the diffraction image update $sdn$ caused by reflection energy leakage have a high similarity to reflectivity. Based on the local similarity estimate (Fomel, 2007a), orthogonalization operator is formed $Omrk\u22121$. It is then applied to the update of the diffractions $Omrk\u22121(sdn)$ removing reflected energy with a high similarity to reflectivity as noise. This procedure can be treated as additional regularization allowing for component separation in the second cascade of the inversions with PWD disabled (equation 8). Orthogonalization can also be applied to the diffractivity model itself rather than to its update. Although both ways are mathematically equivalent to each other in the examples shown below, update orthogonalization demonstrates better reflection crosstalk suppression. After update orthogonalization, sparsity constraints prove to be efficient to denoise diffractions and suppress the remaining crosstalk, major part of which has been eliminated in the update.

### Workflow

The workflow we propose to decompose the input data into reflections, diffractions, and noise is outlined in the flowchart (Figure ^{1}) and is as follows:

- •
Run the first inversion (equation 6).

- •
Run the second inversion (equation 8) and use the output from the first inversion cascade as the starting model.

- •
Restore high wavenumbers in reflections (equation 10); use the output of the second inversion as the starting model.

- •
Restore high wavenumbers in diffractions (equation 11), use the output of the second inversion as the starting model for diffractivity, and use reflectivity from the third step to subtract it from the data.

For shaping regularization of reflections $S\sigma $, dominant local slopes should be estimated in the image domain (conventional Kirchhoff migration image should suffice). For PWD the filter ($Ddata$), dominant local slopes should be estimated in the data domain (zero-offset or stacked section). Both slopes should be precomputed before the inversion.

## SYNTHETIC DATA EXAMPLE

We use the following synthetic to test the proposed chain-inversion approach. Figure ^{2} shows a zero-offset modeling result in a medium with a constant vertical gradient of velocity. The starting velocity at the surface corresponds to $3\u2009\u2009km/s$. The section contains reflections with various dips, several layers of sparse diffractions, and an area in which the density of the scatterers is increased (toward the right edge of the synthetic). Random noise with a 30% magnitude of that of the diffractions and filtered to the signal frequency band is added to the section.

First, we illustrate the advantage of misfit minimization with a diffraction probability weight using the same synthetic (Figure ^{2}), but with no reflections (Figure ^{3}) to compare its performance with least-squares Kirchhoff migration. Reflection/diffraction separation is considered further. Figure ^{4} shows the result of least-squares Kirchhoff migration (equation 5) with no regularization. Noise leakage into the diffraction image domain can be easily noticed. The result of path-summation migration applied to the zero-offset synthetic data (Figure ^{3}) is shown in Figure ^{3}. The path-summation migration image can provide us with a probability of a scatterer at a certain sample location. We use this information by weighting the misfit term in the equation 5 by the path-summation migration operator ($\Vert P(d\u2212Lm)\Vert 22$), which forces the misfit minimization only at probable diffraction locations. No PWD ($Ddata$) is required to emphasize diffractions because there are no reflections. Because the input contains diffractions only, we need only one model for diffractions $m$. No reflection model is considered. Figure ^{4} shows the result of least-squares Kirchhoff migration with a misfit weighted by path-summation integral filter. No regularization has been applied. It can be seen that weighted minimization reduces high-frequency artifacts. The path-summation integral operator emphasizes diffraction apexes remaining stationary under migration velocity perturbation. It acts as a low-frequency filter guiding the inversion toward probable scatterer locations. Indeed, lower frequency content is noticeable and diffractions are more visible on the result of weighted minimization (Figure ^{4}) than on the “conventional” least-squares Kirchhoff migration image (Figure ^{4}). Incorporation of sparsity constraints will remove high-frequency artifacts from the least-squares Kirchhoff migration image and will reduce the area corresponding to a diffraction location in a weighted minimization result. Here, regularization is disabled to show the effectiveness of path-summation operator misfit weighting and avert the subjective choice of a regularization parameter, which will influence the results and make their unbiased comparison harder.

We then apply the proposed workflow to the synthetic containing all the components with regularization enabled. First, we estimate slopes of the dominant events, which are associated with reflections, in the data and image domains. We then generate the data for the first inversion cascade (Figure ^{5}): The zero-offset section (Figure ^{2}) is weighted by PWD ($Ddata$) and path-summation integral filter $P$. We run 100 outer iterations and five inner iterations. The convergence curves are shown in Figure ^{6}. The results of the first inversion are shown in Figure ^{5}. Although subtle reflection remainders are present (e.g., 0.9–1.0 s two-way traveltime [TWT] and 6.0–6.5 km distance), most of the diffracted energy has been extracted with high accuracy.

To reduce these remainders and restore reflections trapped in the null space of the chain of operators, we run the second inversion cascade with no PWD. First, we generate the data to be fit (Figure ^{7}): PWD is disabled, and only the path-summation integral operator is used to weight the misfit. The result of the first inversion cascade (Figure ^{5}) is used as the starting model for the second inversion. We run 10 outer and five inner iterations to restore the reflections and reduce the reflection-diffraction crosstalk. Orthogonalization is required. Figure ^{8} shows the updates to reflections and diffractions after the first external iteration of the second inversion (equation 8). Orthogonalization is applied to remove reflections from the diffraction image update: We can see that reflections in the update are significantly stronger than diffractions and are hard to suppress with thresholding only. Orthogonalization significantly reduces the energy of reflections and diffraction-related update becomes more visible (Figure ^{8}). The orthogonalized update is then added to the diffraction model followed by shaping regularization. Reflection and diffraction images shown in Figure ^{7} are generated. Reflections are recovered from the null space introduced by PWD. Crosstalk is reduced (Figure ^{7}, e.g., 0.9–1.0 s TWT and 6.0–6.5 km distance). Different choices of regularization parameters, including the balance between the numbers of internal and external iterations in the first and second inversions, can be tweaked to improve the results even more.

We then run the inversions (equations 10 and 11) to restore high wavenumbers with the PWD and path-summation integral operators disabled. For diffraction-only recovery (equation 11), besides using the starting model from the second inversion as an additional regularization, the model is masked by the locations of the diffractors. This mask can easily be created from the diffraction image generated as the output of the second inversion cascade (Figure ^{7}). Figure ^{9} shows the reflection and diffraction components restored by the inversions, which were generated by applying Kirchhoff forward modeling operator to the reflectivity and diffractivity output by wavenumber-restoration inversions (equations 10 and 11). We then subtract the recovered reflections and diffractions from the input data to generate the noise estimate. The noise and diffractions are plotted with the same gain. Although the noise section contains some diffracted energy, the approach is highly effective. We then compare the results with the true components used to generate the synthetic. The high correlation between the restored and the “true” components of the synthetic also supports the high accuracy of the workflow. Signal and noise orthogonalization helps to account for amplitude differences between the true and the restored diffractions (Figure ^{9} and ^{9}): Although a high restoration quality is achieved, subtle amplitude differences can be noticed. Because diffraction shapes are restored, this energy can be brought from the noise section back to the diffraction image domain by measuring the similarity between the corresponding components.

The simultaneous inversion for the reflections and diffractions allows to account for the crosstalk between the components and thus improves the confidence of the scatterer position identification. We proceed with the workflow application to the field-data example.

## FIELD-DATA EXAMPLE

We apply the proposed approach to the Viking Graben data set (Keys and Foster, 1998). Klokov and Fomel (2012a) perform diffraction imaging of the Viking Graben data set using the hybrid Radon transform in the dip-angle common image gathers’ domain. Gonzalez et al. (2016) perform diffraction-based depth velocity analysis on the data set.

Figure ^{10} illustrates the zero-offset section. It shows that gradually dipping continuous reflections are interrupted by faults encoded in the diffracted wavefield. Additional diffractions are associated with rough reflector surfaces. Diffractions exhibit significantly lower amplitudes than reflections, and a separation procedure aimed at emphasizing diffractions must be applied before proceeding along the workflow. We separated diffractions (Figure ^{10}) from reflections using the PWD filter (Fomel et al., 2007), and we first imaged them using conventional poststack Kirchhoff time migration (Figure ^{11}). A shallow high-diffractivity layer can be seen in the upper part of the image. Fault surfaces and diffractivity associated with reflector roughness stand out in comparison with the conventional migration image shown in Figure ^{11}. However, because of the overlap of diffraction responses and the presence of noise, the parts of the image between the faults are hard to interpret.

We apply the chain of operators inversion approach to improve the diffraction signal-to-noise ratio and increase the diffraction image reliability. First, we prepare the data to be fit by the inversion (equation 6). The zero-offset section after application of PWD in the data domain ($Ddata$) followed by the path-summation integral filter $P$ is shown in Figure ^{12}. These are the data to be fit, which can be treated as a diffraction probability. Slope values required for PWD in the data domain ($Ddata$) and for shaping regularization in the model space ($S\sigma $) are estimated using the zero-offset section (Figure ^{10}) and conventional Kirchhoff migration image (Figure ^{11}) correspondingly. We run 10 outer and five inner iterations to restore the reflectivity and diffractivity shown in Figure ^{12}. It can be seen that diffractions are significantly more highlighted than in the Kirchhoff migration image of the section after PWD (Figure ^{11}) and appear to be denoised. At the same time, the crosstalk between reflections and diffractions is strong and a lot of strong reflection energy is in the null space of the operator. For example, compare the reflectivity restored (Figure ^{12}) with the conventional Kirchhoff image (Figure ^{11}) at small TWT, e.g., 0.5–0.6 s. For reflection-diffraction crosstalk, for example, refer to the fragment delimited by the 1.75–2.0 s TWT and 19.5–20.0 km distance. Continuous layering associated with the reflections leaks into the diffraction image domain and overlaps with discontinuities along the faults. Events with similar behavior correspond to unconformities at approximately $2s$ TWT between the 18 and 19 km distance.

To address the null-space and crosstalk problems, we generate the data for the second inversion (equation 8), invert for reflections and diffractions (we run 50 outer and five inner iterations), and produce the models shown in Figure ^{13}. Orthogonalization of the diffraction update is shown in Figure ^{14}. The convergence curves are shown in Figure ^{15}. It can be seen that reflections are recovered from the null space (refer to the shallow events) and crosstalk between diffractions and reflections has also been reduced. In particular, the continuous layering events delimited by 1.75–2.0 s TWT and 19.5–20.0 km distance and approximately 2 s TWT between the 18 and 19 km distance are suppressed. There is still a remainder of events with elongated shapes in diffraction image (Figure ^{13}). These events can be connected with rough surface features, which are too irregular to be described by the chosen smoothing radius in shaping regularization penalizing reflections. Less aggressive smoothing and a more detailed dip field in the model space should allow for moving these events back to the reflection image domain. Stronger thresholding can also help.

Finally, we restore high wavenumbers and model reflections (Figure ^{16}) and diffractions (Figure ^{16}). As in the synthetic data example, diffractions in wavenumber-restoration inversion are additionally penalized by the locations’ mask determined from the output of the second inversion cascade (Figure ^{13}). We then subtract their sum from the zero-offset DMO section (Figure ^{10}) and arrive at the following noise estimate (Figure ^{16}). The noise section is plotted with the same gain as the “restored” diffractions. Signal and noise orthogonalization (Chen and Fomel, 2015) is applied to account for the aperture difference between observed (Figure ^{10}) and restored diffractions (Figure ^{16}) (for instance, often only a single diffraction flank can be observed in the data leading to a spuriously high difference or the “noise estimate”) and to bring back some of the energy accidentally leaked to the noise domain (Figure ^{16}). Some residuals that can be noticed in the lower part of the noise section, which can have a signature intermediate between the reflections and diffractions, were not attributed to either. The upper part of the noise section exhibits low-magnitude residuals showing the high accuracy of the method. Besides that, restored diffractions (Figure ^{16}) exhibit a high correlation with the PWD section (Figure ^{10}) also supporting the high efficiency of the method.

Overall impressions of the image shown in Figure ^{13} can be contradictory because of the continuity loss. However, diffractions have spiky and intermittent distribution, which correlates well with the faults and other discontinuities observed in the conventional image (Figure ^{11}). To highlight this statement, we overlay the conventional image (Figure ^{11}) with the estimated diffractivity (Figure ^{17}). A high correlation between continuous layering interruptions and diffractions can be observed (e.g., 1.75–2.0 s TWT and 19.5 and 20.0 km). Moreover, smaller scale faults, which are not evident from the conventional image, are highlighted and denoised with the help of the approach proposed. For instance, refer to the fragment delimited by 0.75–1 s TWT and 17–17.2 km distance. Some faults are not pronounced in the diffraction image. On the other hand, the approach allows for increasing the confidence of the diffractions present in the image and distinguishing them from noise. In other words, even though some events might be lost, we can be confident in the events that were successfully imaged.

## DISCUSSION

We have proposed a workflow to decompose full-waveform input data into reflections, diffractions, and noise by using three cascades of inversions.

We also highlight regularization through orthogonalization as a novel tool to aid the inversion, in which multiple models are estimated with a single modeling operator. Its ability to track down the leakage of the events from one model space to another appears to be crucial for the second cascade of the inversions. Orthogonalization, combined with shaping regularizations, allows for efficient reflection, diffraction, and noise separation.

Shaping regularization allows for flexible model constraints because the models for reflections and diffractions are conditioned independently. For instance, if the regularization terms were incorporated as penalty terms in the objective functions (equations 6, 8, and 9), finding the trade-off between weights for all the terms would be challenging: The two terms would conflict with each other leading to stronger interference between the model spaces.

It should be mentioned that the workflow presented here targets the diffraction phenomenon in 2D rather than in 3D. In particular, in the case of three dimensions, sparsity constraints can be inadequate for edge-diffraction regularization. The latter features are continuous along the edges and act as reflections, but at the same time have a diffractive component in the direction perpendicular to the edge (Moser, 2011). To properly regularize edge-diffractions, sparsity constraints should be accompanied by an anisotropic smoothing operator enforcing continuity along the edges (Merzlikin et al., 2018). At the same time, PWD robustly serving the purpose of diffraction extraction in 2D, in 3D should be replaced by AzPWD — a workflow properly addressing the hybrid signature of edge diffractions (Merzlikin et al., 2016, 2017b). These are the modifications that have to be incorporated into the workflow presented in this paper for it to be extended to three dimensions. Meanwhile, all of the necessary modifications have already been derived and, thus, support the applicability of the reflection-diffraction crosstalk minimization approach presented in this paper in three dimensions.

Finally, the proposed approach (except for analytical time-domain path-summation migration describing diffraction probability) is compatible with any migration algorithms operating in the prestack, poststack, depth, and time domains. In the prestack domain, the PWD filter can be implemented using common-offset sorting of input data. Regularizations are performed in the image domain and, thus, are not dependent on a particular migration algorithm. Analytical expressions exist for path-summation imaging in the prestack time-migration domain (Merzlikin and Fomel, 2017a); however, path-summation integral migration in the depth domain remains challenging (Landa et al., 2006) and, thus, limits the applicability of the approach presented here. The possible workaround is to incorporate probability weights favoring diffractions (Merzlikin and Fomel, 2017b) into path-summation integral expression and use importance sampling strategies for computational load reduction associated with the absence of analytical expressions of path-integral evaluation in the depth domain.

## CONCLUSION

We have developed an approach to decomposing input data into reflection, diffraction, and noise components. The inverted forward modeling operator corresponds to the chain of Kirchhoff modeling, PWD, and path-summation integral migration operators. The chain of these operators accounts for the contribution of diffractions to the optimization, and it guides the inversion toward probable diffraction locations. The reflection model is estimated simultaneously. Separation into different components in an iterative inversion is achieved by regularization: Reflections are forced to be smooth along the dominant local slopes in the image domain, whereas diffractions are penalized by the sparsity constraints. Additional shaping regularization based on local signal-and-noise orthogonalization removes continuous events from diffractivity updates and reduces the crosstalk between the components during iterations. Incorporation of a reflection model into the inversion allows for reflected energy withdrawal from the diffraction image domain and increases the reliability of the diffraction images.

Numerical experiments with synthetic data emulating the high interference between reflections and diffractions demonstrate the high effectiveness of the approach. A field-data example confirms the robust separation of the components by the proposed approach.

The workflow has high computational efficiency. Shaping regularization allows for fast convergence, whereas the cost of the forward modeling operator at each iteration is predominated by Kirchhoff migration.

## ACKNOWLEDGMENTS

We thank E. Landa, L. Decker, R. Biswas, and J. Sun for stimulating discussions and the Texas Consortium for Computational Seismology (TCCS) sponsors for financial support. We also thank two anonymous reviewers for their valuable feedback.

## DATA AND MATERIALS AVAILABILITY

The data and the computational experiments’ scripts and codes are available through Madagascar open-source software package (http://www.ahay.org/wiki/Main_Page).