ABSTRACT

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. Path-summation 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

The path-integral formulation provides an efficient imaging method that does not require subjective and time-consuming velocity picking and produces an image of the subsurface without specifying a velocity model (Landa et al., 2006; Burnett et al., 2011). 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 that correspond to the true media features (Keydar, 2004; Landa, 2004). In other words, the 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 path-integral 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 imaging purposes in correspondence with path-integral formulation in quantum mechanics (Feynman and Hibbs, 1965). In this paper, we follow Burnett et al. (2011) and use the inherent property of the diffractions’ behavior under time-migration transformation — apex stationarity — for time-path-summation imaging implementation for diffraction imaging. We use the term “path summation” to highlight that no weighting is applied before stacking the images as opposed to the path-integral scheme, in which weights represent the 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 the lateral and vertical shifts of the event under velocity perturbation. The VC transformation is continuous and analytical: a velocity step between constant velocity images in theory is infinitely small (Claerbout, 1986; Fomel, 1994, 2003a, 2003b; Burnett and Fomel, 2011; Decker 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 (Burnett et al., 2011). If the constant-velocity images are weighted before stacking by their corresponding velocities, they produce a double path-summation image that can be used for velocity estimation (Schleicher and Costa, 2009; Santos et al., 2016).

The aforementioned 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, the produced images appear to be contaminated by artifacts compared with the images produced by picking velocities (Burnett et al., 2011). The artifacts come from incomplete cancellation of under- and overmigrated flanks of hyperbolas. We refer to these flanks as tails. Tails might overlap with the 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, in which the weighting describing the probability of each image and corresponding velocity is applied before summation.

In this paper, we propose a direct analytical approach in performing path-summation and double path-summation migration as opposed to the stacking of multiple constant-velocity images. Computational efficiency is improved significantly because we do not require multiple VC steps for the calculation of constant velocity images, but instead we compute the path integral in one step at the cost of only two fast Fourier transforms (FFTs). To improve the image quality, we 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 VC concept for migration (Fomel, 2003a). VC describes the 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, 1994, 2003a):  
2Ptv+vt22Px2=0.
(1)
After a substitution σ=t2 and a 2D Fourier transform, the analytical solution for equation 1 takes the following form (Fomel, 2003b):  
P˜(Ω,k,v)=P^0(Ω,k)eik2v216Ω,
(2)
where v is the VC velocity; Ω is the Fourier dual of σ; k is the wavenumber; P^0(Ω,k) is the input stacked and not migrated data, which correspond to v=0; and P˜(Ω,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 (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 the 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 stacking constant velocity images generated by VC (Burnett et al., 2011). Theoretically, the velocity step between images generated by VC can be infinitely small, and path-summation migration in this case will correspond to the following integral over the velocity dimension:  
IPI(t,x,va,vb)=vavbP(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 because it stacks all the constant velocity images generated by VC.
Note that the stacking of images can be performed in the (t,x) and (Ω,k) domains. For the latter case, it takes the following form:  
IPI^(Ω,k,va,vb)=vavbP^0(Ω,k)eik2v216Ωdv,
(4)
where P^0(Ω,k) corresponds to zero-offset or stacked data and does not depend on the migration velocity. Therefore, it can be taken outside of the integral:  
IPI^(Ω,k,va,vb)=P^0(Ω,k)vavbeik2v216Ωdv.
(5)
We can treat the remaining integral as a filter FPI(Ω,k,va,vb)=vavbe(ik2v2)/16Ωdv, which depends on the frequency Ω, wavenumber k, and velocity integration range limits va and vb. Moreover, the filter has an analytical expression: The integral of the complex exponential function is proportional to the imaginary error function (erfi):  
FPI(Ω,k,va,vb)=ei5π42Ωπkerfi(ei3π4kv4Ω)|vavb.
(6)
Thus, path-summation migration according to equation 5 amounts to filtering the data with the filter given in 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 the two FFTs needed to transform the data to the frequency-wavenumber domain and back to the time-space domain. The filtered image can be treated as if it was calculated by the 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 2a. The path-summation migration image contains artifacts because of incomplete cancelation of two hyperbolas, which are still prominent in the image: one corresponding to the lowest constant velocity va image (an undermigrated tail) and the other corresponding to the highest velocity vb image (overmigrated tail). Those two tails are not canceled out during stacking of constant-velocity images, whereas hyperbolas’ flanks within the velocity integration range interfere destructively and, therefore, suppress each other. When path-summation migration is applied to real diffraction images, the tails may interfere with the useful signal and create artifacts.

To extend analytical path-summation migration formulas to the 3D isotropic poststack case, one needs to substitute the absolute value |k| of the wavenumber vector k=(kx,ky) instead of a 1D wavenumber k in the expressions given above. Prestack 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. The corresponding analytical evaluation formulas are provided in Appendix B.

Attenuating tail artifacts

Artifacts in the path-summation images appear due to the absence of weights that describe the probability of the image migrated with a certain velocity and taper out “unlikely” solutions (Landa et al., 2006; Landa, 2009; Schleicher and Costa, 2009). We use analytical Gaussian weighting of images being stacked in path-summation 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).

An analytical weighting term can be added to the VC phase shift as follows:  
P˜(Ω,k,v,β,vbias)=P^0(Ω,k)eik2v216Ωβ(vbiasv)2,
(7)
and the corresponding integral expression describing path-summation migration takes the following form:  
IGPI^(Ω,k,va,vb,β,vbias)=vavbP^0(k,Ω)eik2v216Ωβ(vbiasv)2dv=P^0(Ω,k)vavbeik2v216Ωβ(vbiasv)2dv=P^0(Ω,k)·FGPI(Ω,k,va,vb,β,vbias),
(8)
where filter FGPI(Ω,k,va,vb,β,vbias) is evaluated analytically as follows:  
FGPI(Ω,k,va,vb,β,vbias)=πe(βvbias2+β2vbias2γ(Ω,k,β)2)2γ(Ω,k,β)erfi(γ(Ω,k,β)v+βvbiasγ(Ω,k,β))|vavb.
(9)
Here, γ(Ω,k,β)=iik2/16Ω+β, β is a parameter controlling the contribution of each velocity to the final image, and vbias is a velocity bias. By choosing a constant value of vbias for the whole section between va and vb, we can enhance the apexes of summed hyperbolas and decrease the tails’ contribution to the final image. Because in a section, velocity varies and vbias is constant, over- or undermigrated diffraction curves may get emphasized at locations where the true velocity is higher or lower than vbias, respectively. This is a possible drawback of the Gaussian weighting approach.

We analyze the response of a single diffraction hyperbola modeled with 1.5  km/s velocity (Figure 1a) to path-summation and Gaussian-weighted path-summation migrations. The synthetic model spectrum in the Ωk domain is shown in Figure 1b. The path-summation migration filter magnitude and its phase spectrum are shown in Figure 3a and 3c. This 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 phase-shift values. This area corresponds to the diffraction hyperbola apex remaining stationary under time-migration velocity perturbation. Nonzero wavenumber events have lower corresponding filter coefficients that taper out to the spectrum boundaries. Figure 2a shows the result of path-summation migration applied to diffraction hyperbola. Several lobes with nonzero wavenumber values appear to be noticeable in the Ωk domain (Figure 2c). Corresponding phase-shift values (Figure 3c) are not zero and lead to displacement of the corresponding events on the path-summation migration image (Figure 2a). These events are associated with tail artifacts.

The Gaussian weighting scheme spectrum magnitude and phase are shown in Figure 3b and 3d. The result of corresponding path-summation migration in Ωk domain is shown in Figure 2d. The Gaussian weighting scheme tapers nonstationary lobes and, therefore, it suppresses the tails. Figure 2b 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. The following sections illustrate the applicability and effectiveness of the proposed techniques.

FIELD DATA EXAMPLES

Reflection/diffraction separation

Diffraction amplitudes can be several orders of magnitude lower than those of reflections (Klem-Musatov, 1994). To highlight the 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), Berkovitch et al. (2009), Dell and Gajewski (2011), Tsingas et al. (2011), Rad et al. (2014), and Bauer et al. (2015) use optimal stacking of diffracted energy along diffraction traveltime curves. Harlan et al. (1984), Papziner and Nick (1998), Taner et al. (2006), Fomel et al. (2007), Reshef and Landa (2009), and Klokov and Fomel (2012) decompose full-waveform seismic records into diffracted and reflected components; Kozlov et al. (2004), Moser and Howard (2008), Koren and Ravve (2011), Sturzu et al. (2013), Klokov and Fomel (2013), and Popovici et al. (2015) modify a Kirchhoff migration kernel to eliminate specular energy coming from the first Fresnel zone, and they image only diffractions.

To perform reflection/diffraction separation, we used plane-wave destruction (PWD) 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 the 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 data set (Claerbout, 2005; Fomel et al., 2007). A stacked section of diffractions separated using a plane-wave 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. The integration limits corresponded to 1.42.6  km/s. Diffractions collapse to form an image of faults.

A path-summation migration image with Gaussian weighting evaluated with 1.42.6  km/s integration limits, vbias=2.0  km/s, and β=10 is shown in 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 the double path-summation integral by the path-summation integral, and the corresponding diffraction image generated with this velocity model are shown in Figure 4b and 4c. The 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 the double path-summation velocity has the highest resolution. A full image with the reflections and diffractions imaged using the same velocity field is shown in Figure 4d. Discontinuities in the 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 data set 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 geologic interpretation of diffraction images and their comparison to discontinuity-type attributes are provided in Tyiasning et al. (2016).

The input stack is shown in Figure 5a. The 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 data set. The extracted velocity distribution is shown in Figure 6. The velocity limits used for path-summation integration correspond to 2.93.1  km/s. They were chosen based on an a priori migration velocity distribution.

We used automatically extracted velocity to migrate initial (Figure 7) and diffraction data (Figure 8). Images shown in Figures 7 and 8 are extracted from the whole migrated volume along the target horizon/interface between a coal layer and tight gas sand. The diffraction image (Figure 8) appears to have 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 68.5  km/s and crosslines 0.5–2.5 km is prominent in both images, but it is emphasized more 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 with the double path-summation velocity image (Figure 9a). However, the latter one appears to be noisier.

Several 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’ parts 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 prestack time migration (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 — from 1.5 to 4.5  km/s. The image appears to be much less focused; however, the pattern of the elliptical structure can still be identified. The unweighted path-summation migration image (Figure 9d) appears to be significantly blurred due to the overlap between the tails and the useful signal. The Gaussian weighting path-summation scheme allows us to significantly suppress tails and highlight actual apexes of hyperbolas (Figure 9f). We used a 3.1  km/s velocity bias with β=10. In comparison with Figure 9c, the Gaussian weighting path-summation image has lower spatial resolution. This observation is connected with the wider integration range used for Gaussian weighting path-summation migration in comparison to the image in Figure 9c acquired with 2.93.1  km/s velocity limits. When a wider integration range is used, the 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 cancelation of hyperbolas’ flanks corresponding to the lowest and highest velocities. We propose to apply Gaussian weighting scheme (equation 8) 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 path-summation 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 a 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 (Appendix B) by a 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 the numerical summation. The proposed approach is based on VC, which describes continuous image transformation with perturbation of the migration velocity, and in theory it 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 removed by a Gaussian weighting scheme. In this case, path-summation migration can still be performed analytically. Analytical evaluation formulas for double-path summation allows for efficient migration velocity extraction based on diffractions.

ACKNOWLEDGMENTS

We thank the TCCS sponsors and the Geofrac Consortium for financial support and W. Burnett, D. Cooke, and E. Landa for inspiring discussions. We thank two anonymous reviewers for their constructive suggestions.

PRESTACK PATH-SUMMATION MIGRATION

Keydar et al. (2009) describe prestack data imaging using path-integral formulation in which the curvature of the diffraction hyperbola acts as a parameter describing different images being stacked. In this appendix, we derive a direct and analytical formula for unweighted prestack path-summation imaging based on the extension of the VC technique to the prestack domain. Fomel (2003a, 2003b) describes VC phase shift for a prestack case:  
P˜(Ω,k,v)=hP˜0(Ω,k,h,v0)eik2(v02v2)16Ω+4iΩh2(1v21v02),
(A-1)
where h is a constant half-offset value and v0 is a velocity for the first migration run applied before cascade of VC transformations. The integral describing prestack path-summation migration takes the following form:  
I(Ω,k,va,vb)=vavbhP˜0(Ω,k,h,va)eik2(va2v2)16Ω+4iΩh2(1v21va2)dv=hP˜0(Ω,k,h,va)vavbeik2(va2v2)16Ω+4iΩh2(1v21va2)dv=hP˜0(Ω,k,h,va)F(Ω,k,h,va,vb),
(A-2)
where F(Ω,k,h,va,vb) is a prestack path-summation migration filter that can be evaluated analytically as  
F(Ω,k,h,va,vb)=eiπ4Ωπkeiva2k216Ω4iΩh2va2[e2λ(Ω,k)μ(Ω,h)erf(λ(Ω,k)v+μ(Ω,h)v)+e2λ(Ω,k)μ(Ω,h)erf(λ(Ω,k)vμ(Ω,h)v)]|vavb,
(A-3)
where erf is the error function, λ(Ω,k)=eiπ4k/4Ω and μ(Ω,h)=2ieiπ4hΩ.

DOUBLE PATH-SUMMATION MIGRATION

Double path-integral formulation was proposed by Schleicher and Costa (2009) and Santos et al. (2016). In unweighted path-summation migration framework, it can be described by the following equation:  
IDPI(t,x,va,vb)=vavbvP(t,x,v)dv.
(B-1)
The velocity distribution can be acquired through a smooth division of the double path-summation integral by the path-summation integral:  
v(t,x)=IDPI(t,x,va,vb)/IPI(t,x,va,vb).
(B-2)
We propose to evaluate the double path-summation integral analytically as well:  
IDPI(Ω,k,va,vb)=vavbvP^0(k,Ω)eik2v216Ωdv=P^0(Ω,k)vavbveik2v216Ωdv=P^0(Ω,k)·FDPI(Ω,k,va,vb),
(B-3)
where the filter FDPI(Ω,k,va,vb) has the following analytical expression:  
FDPI(Ω,k,va,vb)=8iΩk2eik2v216Ω|vavb.
(B-4)

To extract the migration velocity as described by equation B-2, we need to evaluate the double path-summation integral and the path-summation integral. These evaluations can be done analytically as described by equations 5, 6, and B-3. A smooth division procedure (Fomel, 2007) is applied in equation B-2 to acquire a regularized time-migration velocity estimate, which can then be applied to image the subsurface using any available time imaging algorithm.

Freely available online through the SEG open-access option.