Separating and imaging diffractions of seismic waves in the full-azimuth dip-angle domain


 Seismic diffractions are ideal carriers of information on small-scale, discontinuous objects and can therefore be used to detect these geologic objects. However, recognizing diffractions is difficult because specular reflection with strong energy masks the weak diffraction. In this study, we propose a diffraction separation and imaging method based on a Mahalanobis-based and phase-based attenuation function used to modify the Kirchhoff migration formula in the full-azimuth dip-angle domain. In this domain, reflections are restricted to within the first Fresnel zone and are distributed in the vicinity of the stationary point, while diffractions are located across a wide range of azimuth and dip angles. Synthetic and field data applications suggest that this new method can effectively separate and image diffractions. The results also demonstrate the efficiency of the new method in clarifying subsurface small-scale objects, which can provide finer information about these structures for seismic interpretation.


Introduction
Diffractions are generated from small-size geologic structures, such as faults, fractures and collapsed columns, when seismic waves encounter them. These subsurface geologic discontinuities are related to oil and gas exploration (Laubach et al. 2000), and may also lead to water inrush, gas outburst and other geologic hazards related to coal exploration (Peng & Zhang 2007). As an ideal carrier, diffractions are more effective than reflections in revealing these geologic discontinuities because of the superior illuminations. However, standard seismic processing is based on enhancing reflections and suppressing diffractions to some degree (Khaidukov et al. 2004). Therefore, an effective method specific to imaging diffractions is necessary to detect these discontinuities (Rieber 1936;Krey 1952;Hagedoorn 1954). Although the energy of diffractions is weaker than that of specular reflections, the difference is about one or two orders of magnitude (Klem-Musatov 1994;Klokov et al. 2010). Therefore, the strong reflections have a great masking effect on the diffraction imaging results and need to be removed to enable extraction of the details of weak diffractions.
Diffraction imaging methods have been extensively reported. The first kind of diffraction imaging method suppresses reflections and enhances diffractions simultaneously in the imaging process. Considering the phase correlation of diffractions and the corresponding statistical criteria, Landa et al. (1987) proposed a diffraction detection method in the common-offset domain. By stacking the energy emitted from diffraction edges, Kanasewich & Phadke (1988) identified the edge diffractor using moveout and amplitude corrections. For detecting local heterogeneities, Landa & Keydar (1998) used the common-diffraction-point section (D-section) to enhance the amplitudes of diffractions.  Using the multifocusing method, Berkovitch et al. (2009) developed a local-time-correction formula to focus diffraction energy for imaging local heterogeneities. Based on the energy-angle distribution of diffractions and reflections, the Kirchhoff migration method was modified by Zhu & Wu (2010) to form a local-angle domain diffraction imaging method. Similarly, Koren & Ravve (2011) presented a new imaging technique in the angle domain by modifying the classic Kirchhoff migration formula to obtain high-resolution information. Using the dynamic patterns of reflections and diffractions, a Mahalanobis-based diffraction imaging method was developed by Zhao et al. (2016), which can effectively remove reflections and automatically image diffractions in the shot domain.
Some scholars have suggested another type of diffraction imaging method: separating diffractions in advance and then imaging them. Harlan et al. (1984) removed specular reflections and separated diffractions using the local slant stack method in the post-stack seismic data. Based on the different geometrical focusing characteristics, a focusingmuting-defocusing strategy was proposed by Khaidukov et al. (2004) to suppress reflections and image diffractions. In addition, the plane-wave destruction (PWD) method is recognized as a popular way of enhancing diffractions    (Claerbout 1992). Based on the PWD method, Lin et al. (2018a) extended the PWD method with a regularization constraint to enhance the stability while separating diffractions.
The third kind of diffraction imaging method is based on the different behaviors of reflections and diffractions in the dip-angle domain. When the migration velocity is accurate, reflections show a concave shape with a stationary point at the real dip angle and diffractions appear flat when observed at their real diffractor locations (Audebert et al. 2002;Qin et al. 2005;Reshef & Landa 2009;Lin et al. 2018b). According to this difference in characteristics between diffractions and reflections, Zhang & Zhang (2014) developed a migration scheme by muting the events in the Fresnel zones in the shot and opening-angle gathers. In addition, other effective methods for extracting angle domain common imaging gathers (CIGs) can improve the diffraction images (Cheng et al. 2012;Liu et al. 2018aLiu et al. , 2018b. Information on the Fresnel zone and the stationary point is indispensable in many diffraction imaging methods, and the method for obtaining this information is complicated (Sun et al. 2019). In this study, the proposed diffraction imaging method can automatically suppress reflections without information about the Fresnel zone and stationary point. Considering the different characteristics of diffractions and reflections in the full-azimuth dip-angle domain, we present a Mahalanobis-based and phase-based attenuation function to destroy reflections and enhance diffractions. Then, we modify the Kirchhoff migration algorithm with this function to image diffractions. This study is an extension of a Mahalanobis-based diffraction imaging method (Zhao et al. 2015). It extends the previous method to the fullazimuth dip-angle domain by analyzing the behaviors of reflections and diffractions in the angle domain. That allows the proposed method to be used to image diffractions in the 3D case. Furthermore, diffraction events with the same instantaneous phase are distributed over a wide range of azimuth-dip angles. Therefore, a higher value of the stack results of the instantaneous phase is generated for the diffraction. This behavior is used to modify the attenuation function to enhance diffractions. In this paper, first, we review and describe the kinematic and dynamic characteristics of reflections and diffractions in the full-azimuth dip-angle domain. Then, the modified Kirchhoff migration algorithm with the attenuation function is provided. It is based on the Mahalanobis distance and the Hilbert transform to suppress reflections and enhance diffractions. Finally, one synthetic and one field example show the effectiveness of the new method in clarifying faults and other small-scale scattering objects.

Basic geometrical theory of reflections and diffractions
The kinematic characteristics of reflections and diffractions are introduced using three sketches (figures 1-3), and the zero-offset survey is adopted for simplicity. In these sketches, Letters I and O are the image point and the corresponding vertical projection, respectively. Let M and the ON axes denote the zero-offset point and the north direction, respectively. The observed dip angle is the angle between IO and IM, and the angle between OM and the north direction ON is defined as the observed azimuth angle . Based on Snell's law, we observe the reflection and locate the stationary point when IM is perpendicular to the reflecting interface (figure 1). In other words, the stationary point is obtained when the observed dip and azimuth angle are equal to the real dip and azimuth angle of the reflector. In accordance with the geometrical theory of diffraction (Keller 1962), the edge diffraction fronts are rounded when a ray is normally incident on the reflector edge, as shown in figure 2, and appear conical when the incident ray hits the edge obliquely. Therefore, the edge diffraction can be observed when the observation direction OM is perpendicular to the boundary. The response of the edge diffraction is therefore limited by the azimuth angle. In accordance with the Huygens principle, the point diffraction can be generated at different dip and azimuth angles, as shown in figure 3; the point diffraction can be observed from a wide range of dip and azimuth angles.

Kinematic and dynamic characteristics of reflections and diffractions in the full-azimuth dip-angle domain
We construct a 3D geologic model (shown in figure 4) to further to illustrate the kinematic and dynamic characteristics of reflections and diffractions. The geologic model consists of a horizontal reflector and a point diffractor at a depth of 600 m with a constant velocity of 2000 m s −1 . We assume that the reflectivity is generated by density variations and print reflection coefficients in the xz section. Through the Kirchhoff modeling method, the zero-offset data are generated using 401 × 201 receivers and shots placed at intervals of 10 m in the x and y directions. Then, we extract the fullazimuth dip-angle gather at the position of x = 0.5 km and y = 0.5 km to illustrate the amplitude characteristics of the reflector responses in figure 5a. In the corresponding time slice (figure 5b), a reflection pattern in the full-azimuth dipangle gathers is constructed and limited by the dip angle near the stationary point. This area depends on the Fresnel zone in the angle domain. In the azimuth angle slice (figure 5c), the reflection event has the shape of a smile with a stationary point.
The full-azimuth dip-angle gather at x = 1 km and y = 0.5 km above the edge diffractor is shown in figure 6a. In the time slice (figure 6b), the weak edge diffractions are obviously masked by the strong reflections and have the opposite polarity (red arrow) outside the energy band. The edge diffraction is distributed near the azimuth angle, which is equivalent to the reflector azimuth (Keller 1962). The diffraction is divided into two parts in the azimuth angle slice (figure 6c); the shadow part under the reflection (red arrow) has an opposite polarity to of the reflection, while the light part (yellow arrow) has the same polarity as that of the reflection. Figure 7a shows the full-azimuth dip-angle gather at x = 1.5 km and y = 0.5 km, the same position as that of the point diffractor. In the time slice (figure 7b), the amplitude is zero (red arrows) due to the limited range of observation. The point diffraction energy is attenuated with increase in the dip angle because of the geometric spreading, and unaltered when the azimuth angle changes. Note that the point diffraction event does not have an energy band and behaves flat in the azimuth angle slice (figure 7c). These slices indicate that we can observe the point diffraction from a wide range of dip and azimuth angles.
In the full-azimuth dip-angle domain, reflections are restricted to within the first Fresnel zone and are distributed in the vicinity of the stationary point, but diffractions are located in a wide range of azimuth-dip angles. Furthermore, the reflection energy is focused on the first Fresnel zone, which is much stronger than the energy outside of the Fresnel zone. In contrast, the diffraction energy is distributed over a wide range of azimuth-dip angles. These differences between reflections and diffractions are revealed in full-azimuth dip-angle gathers.

Mahalanobis-based attenuation function
The conventional Euclidean distance (ED) is the most commonly used distance measures. However, the ED cannot rea-sonably measure the difference between sample and sample/whole when the properties of the samples have different orders of magnitudes. Fortunately, by eliminating the effect of the orders of magnitude, the Mahalanobis distance (MD) can measure the difference reasonably when the number of samples is larger than the dimensions of the samples. Therefore, the MD is a feasible way to measure the distance between a sample point and a whole dataset and can be used to estimate whether a sampling point belongs to a dataset or not (Maesschalck et al. 2000). The distance can be formulated as follows: where ⃗ x is a simple point that can be defined as ⃗ x = [x 1 , x 2 , ..., x n ], ⃗ u = [u 1 , u 2 , ..., u n ] is the mean vector and S is the covariance matrix. In this study, the MD can be used to measure the distance between a sample point x (i,j) in the full-azimuth dip-angle gather and the corresponding time slice dataset of this gather. For simplicity, the MD is redefined as: where u and are the mean and the standardized deviation of the time slice dataset, respectively. In the time slice dataset, the reflection energy is concentrated in the Fresnel zone. Thus, the average energy of the corresponding time slice would be much lower than the reflection energy. A high value would be obtained when we measure the MD between the reflection and the corresponding time slice dataset. On the contrary, the diffraction energy is distributed over wide azimuth-dip angles in the time slice dataset. The average energy of the corresponding time slice is about same as the diffraction energy and, thus, a small MD would be attained. Based on the difference between reflections and diffractions, we construct an attenuation function to suppress reflections. Considering the dynamic characteristics of diffractions (Zhao et al. 2015), the attenuation function is written as: where dm(x)∕dm(x) is a decreasing attenuation coefficient when the MD is larger than the mean of the MD (dm(x)). To further suppress reflections, we use the instantaneous phase calculated by the Hilbert transform to develop the attenuation function (Taner et al. 1979), which can be formulated as  (Barnes 1996): ) .
The A(i, j, t) represents the full-azimuth dip-angle gather, and the A * (i, j, t) is the corresponding results obtained using the Hilbert transform. For diffraction, a higher value of the stack results of the instantaneous phase would be generated as the diffraction events are distributed over a wide range of azimuth-dip angles. Therefore, equation (3) can be modified as: where (t) represents the stack results from the instantaneous phase in time t of the full-azimuth dip-angle domain. The 3D example is applied to show the feasibility of the proposed method in suppressing reflections and strength-ening diffractions. The attenuation result of the reflection is displayed in figure 8a. In the corresponding time slice (figure 8b), the reflection in the Fresnel zone is destroyed and its concave shape disappears in the azimuth angle slice (figure 8c). Figure 9a shows the attenuation results of the edge diffraction. In the time slice (figure 9b), the reflection and the shadow part (yellow arrow) of the edge diffraction are smeared, and the light part (red arrow) of the edge diffraction is strengthened. However, the polarities of the two parts of the edge diffractions are reversed, so the summation of diffraction would be enhanced after stacking. The disappearance of the reflection in the azimuth angle slice (figure 9c) also shows the good efficiency of the proposed method in suppressing reflections. The point diffractions in figure 10a-c are better strengthened compared to those in figure 7a-c. The example demonstrates that the attenuation function is feasible for removing reflections and enhancing diffractions.

The modified Kirchhoff migration method
The classic Kirchhoff migration method can be expressed as (Hubral et al. 1996): For imaging diffractions, we incorporate the proposed damping function in the Kirchhoff migration method to suppress reflections. The modified migrations method can be expressed as: where u(⃗ r, t) represents the seismic records. The attenuation factor a(x) is calculated using equation (5) and the sample x is formed by In this study, the one-way traveltimes t s and t r are estimated using the Cheng's method (Cheng et al. 2012) and the factor w(⃖⃖ ⃗ m,⃗ r) for geometric spreading referees to Lee et al. (2004). The illumination-azimuth is estimated using (Cheng et al. 2011) where (x I , y I ), (x s , y s ) and (x r , y r ) are the horizontal locations of the image point, the source and the receiver, respectively.  The illumination-dip is obtained by (Cheng et al. 2012) where ⃖⃖ ⃗ p s and ⃖⃖ ⃗ p r are the incident and scattering slowness vectors of the image point and ⃗ z is a vertical downward unit vector. The azimuth is defined as the angle between the north direction and the horizontal component of ⃖⃖ ⃗ p s + ⃖⃖ ⃗ p r . The proposed attenuation function can destroy strong reflections without the known apertures and stationary point position, which are important information in the conventional imaging method. In this study, the information on wide azimuth-dip angles is used to image diffractions. Therefore, the proposed imaging method with the attenuation function can suppress reflections and image diffractions effectively.

Synthetic data set
The performance of the proposed method is tested with a 3D geologic model, shown in figure 11. This model includes six faults whose displacements range from about 80 to 200 m and six point diffractors all with a radius of 5 m. There are 601 × 201 receivers/shots placed at intervals of 10 m in both the x and y directions. The time length is 1.5 s, and the sample interval is 1 ms. The common-offset data are obtained using the Kirchhoff modeling method with a Ricker wavelet of the dominant frequency of 30 Hz. Figures 12 and 13 are the two profiles of the synthetic data. The reflection is generated by the smooth part of the reflector, and the edge diffraction and the point diffractions are generated from the reflector boundary and the point diffractor, respectively. The boundary of the geologic model also caused several other waves (shown with red arrows). show that the reflection events are suppressed and the diffraction events are enhanced.
The modified imaging method is tested using the common-offset data generated from the 3D geologic model. Two profiles of the imaging result are provided in figures 16 and 17, respectively. In the two profiles, the reflection has been suppressed as expected, and the residual reflection energy is relatively weak. The edge of the reflector (fault breakpoint) can be clearly identified as indicated (yellow  arrows). Six point diffractors (red arrows) are highlighted and revealed accurately. The disturbing waves as mentioned above do not affect the imaging result, indicating that the diffraction imaging method can enhance and image diffractions even in the background of strong reflections. The synthetic data example verifies the feasibility and effectiveness of the new imaging method in revealing faults and small-size scattering objects.
However, the diffraction energy is so weak that imaging diffraction is difficult, especially in the case of a low signalto-noise ratio (SNR). We add the Gaussian noise in the synthetic data to further test the proposed method in the low SNR case. Figure 18 parts a and b display the time slice and azimuth angle slice with the same positions, as shown in figure 14 parts b and c, respectively. It was found by comparing figures 14 and 18 that the noise has little effect on the dominant reflection events but seriously disrupts diffraction events. The results of the proposed method are shown in figure 19a and b. The diffraction was damaged to some degree, but it can still be distinguished using the proposed method. This shows the proposed method can enhance diffraction in the case of low SNR data.

Real data set
In coal mining, faults and collapsed columns may lead to water inrush, gas outburst and other geologic hazards (Peng & Zhang 2007). Predicting the location of these small-scale geologic structures is necessary to avoid potential accidents and reduce production cost. As good carriers of the detailed information, the diffraction can accurately locate the corresponding structures. However, weak diffractions are masked or disturbed by strong reflections or suppressed to some degree during conventional seismic processing. Therefore, we propose a diffraction imaging method with the Mahalanobisbased attenuation function for enhancing and imaging the weak diffractions.
The proposed imaging method was applied to 3D land data from West China. The inline and crossline intervals were 10 m, and the sample interval was set to 1 ms. A collapsed column was revealed by the interpreter, and was considered a highly credible result. The collapsed column penetrates the coal seam at a depth of about 600 m. When the coal seam is mined, the in situ stress is redistributed and further enhances the permeability of the collapsed column. The collapsed column may become a conduit for gas or confined water and lead to the water inrush or gas outbursts. The position of the collapsed column is 34 and 52 in the crossline and inline directions, as shown in figure 20. A crossline profile (crossline = 34) and an inline profile (inline = 52) with the collapsed column (red arrows) are displayed in figure 21a and b, respectively. The conventional Kirchhoff migration method was applied to the field data, as shown in figure 22a 14  and b. In this profile, the dominant reflection events are strong with good continuity, which can be used to reveal the reflector. However, information on small-scale objects is lacking and even masked by the strong reflections. This leads to the challenge of identifying the collapsed column (red arrows), so the proposed imaging method was adopted to obtain detailed information on the structure. The corresponding results are shown in figure 23a and b. The dominant reflection events appearing in conventional imaging methods are suppressed by the proposed method, and the collapsed column (red arrows) information is prominent in this profile.
Four time slices were extracted to further to verify the effectiveness of the proposed imaging method. Figure 24a and c show the results of the conventional imaging method for 318 and 352 ms, respectively. The two time slices lack finer information on the collapsed column, as the conventional method is based on enhancing reflections. It is difficult to obtain the boundary of the collapsed column. However, as shown in figure 24b and d, the proposed method can obtain the detailed information by strengthening diffractions. The result can be used to accurately identify and delimit the collapsed column during seismic interpretation. The field example demonstrates that the proposed diffraction imaging method is effective and reliable for detecting the subsurface discontinuities.

Conclusion
In the full-azimuth dip-angle domain, the reflection energy is concentrated on the Fresnel zone and the diffraction energy is distributed over wide azimuth-dip angles. This results in a higher MD value for reflections and a lower MD value for diffractions. Because of this difference, the attenuation of reflections can be performed in the full-azimuth dip-angle domain. Therefore, we propose the Mahalanobisbased and phase-based attenuation function. The function can suppress reflections and enhance diffractions effectively without calculating of the apertures and the stationary point position. For imaging diffractions, the modified Kirchhoff migration method is presented with the attenuation function.
The synthetic data and real data examples demonstrate the feasibility and reliability of the proposed method in highlighting diffractions. The proposed diffraction imaging method can accurately locate the subsurface discontinuities in the background of strong reflections. The method can be used to predict faults, crack zones, collapsed columns and other subsurface discontinuities. It offers a reliable way to avoid potential accidents during coal mining, and it can be further exploited in oil, gas and other resource exploration.