Separation of diffracted waves via SVD filter

Diffracted seismic waves may be used to help identify and track geologically heterogeneous bodies or zones. However, the energy of diffracted waves is weaker than that of reflections. Therefore, the extraction of diffracted waves is the basis for the effective utilization of diffracted waves. Based on the difference in travel times between diffracted and reflected waves, we developed a method for separating the diffracted waves via singular value decomposition filters and presented an effective processing flowchart for diffracted wave separation and imaging. The research results show that the horizontally coherent difference between the reflected and diffracted waves can be further improved using normal move-out (NMO) correction. Then, a band-rank or high-rank approximation is used to suppress the reflected waves with better transverse coherence. Following, separation of reflected and diffracted waves is achieved after the filtered data are transformed into the original data domain by inverse NMO. Synthetic and field examples show that our proposed method has the advantages of fewer constraints, fast processing speed and complete extraction of diffracted waves. And the diffracted wave imaging results can effectively improve the identification accuracy of geological heterogeneous bodies or zones.


Introduction
Accurate identification of faults, thin-outs, karsts, lens bodies, collapse columns, fracture zones and other heterogeneous regions is one of the important, though difficult, goals of seismic data processing and interpretation (Yilmaz 2001;Khaidukov et al. 2004). Reflected waves are a reflection of the interface morphology of the subterranean strata, which is mainly characterized by lateral continuity. However, the interface of geological heterogeneous bodies is rough or there is no horizontal continuous interface at all. When the downing wave fields meet the geological heterogeneity, it will cause the seismic energy to diverge to the surroundings. Therefore, the seismic response characteristics of geological heterogeneous bodies (such as breakpoints, pinch-outs, karsts and collapse columns) often appear as diffracted waves (Bansal and Imhof 2005;Fomel et al. 2007;Decker et al. 2015). The diffracted waves are seismic responses caused by uneven geological bodies in the strata, which carry high-resolution, potentially even super-highresolution, geological information (Neidell 1997;Khaidukov et al. 2004;Rad et al. 2005;Sturzu et al. 2015). Therefore, diffracted waves are one type of wave field used to effectively identify and track geologically heterogeneous bodies or zones. Since the early 1950s, researchers recognized the potential of using diffracted waves to identify geological heterogeneous zones and began using diffracted waves to detect small faults (e.g., Krey 1952;Angona 1960;Harper 1965;Kovalevsky 1971;Landa and Maximov 1980;Landa et al. 1987;Kanasewich and Phadks 1988), characterize karst edges (e.g., Decker et al. 2015) and identify fractures (e.g., Popovici et al. 2014;Sturzu et al. 2014Sturzu et al. , 2015 and uneven geological bodies (e.g., Landa and Keydar 1998). The research results of Sturzu et al. (2014Sturzu et al. ( , 2015 show that the single diffracted wave imaging results can effectively depict fractures in carbonate reservoirs. Decker et al. (2015) believed that diffracted wave imaging results can effectively Edited by Jie Hao * Hong-Yan Shen shenhongyan@xsyu.edu.cn 1 improve the horizontal resolution of a single karst cave and effectively identify heterogeneous regions below the reflection resolution. However, diffracted waves have weak signals, which are usually concealed under the background of reflected waves and other noises (Klem-Musatov et al. 1994). Therefore, to effectively use diffraction wave information to characterize the subsurface, one must first separate it from reflection signals. In addition to surface waves and random noise, reflected waves have become the main interference noise of diffracted waves. How to completely separate diffracted waves from reflected waves is one of the key issues that determine whether diffracted waves can be fully and effectively utilized. Researchers have been working on this for a long time. Khaidukov et al. (2004) separated diffracted waves by focusing and removing reflected waves and then anti-focusing the remaining seismic signals. Taner et al. (2006) and Kong et al. (2017) extracted diffracted waves by suppressing smooth and continuous quasilinear reflected waves based on the difference in the time distance curves between diffracted and reflected waves in plane wave records. Moser and Howard (2008) built the anti-stability phase function to suppress the reflected wave field while enhancing diffracted waves. Landa et al. (2008) and Zhu et al. (2013) used local a dip filter to achieve diffraction and reflection separation in the migration dip domain. Based on the theory of diffraction coherence summation, Berkovitch et al. (2009) proposed using the local time correction formula to parameterize the diffraction travel time curve and then stacked diffraction events while suppressing reflected waves. Klokov and Fomel (2012) used the Radon transform to achieve the separation of diffracted and reflected waves with common imaging point gathers. Zhang and Zhang (2014) and Li et al. (2018) proposed removing the Fresnel zone in common imaging point gathers to enhance diffracted waves. Although a variety of diffracted wave separation methods have been developed to date, each method has its own applicable conditions, and the wave field separation quality still needs to be further improved. Therefore, developing new and more efficient methods and techniques for diffracted wave separation is still one of the challenging problems in the use of diffracted waves.
As a method of matrix factorization, singular value decomposition (SVD) was introduced into the field of seismic signal processing in the 1980s. Freire and Ulrych (1988) first tried to separate the upgoing and downgoing waves of VSP data via SVD. Jackson et al. (1991) extensively analyzed the principle of SVD processing seismic data. Bekara and Baan (2007) used SVD filtering to improve the signalto-noise ratio (SNR) of seismic data. Porsani et al. (2009) used SVD to filter ground roll waves. Gao et al. (2013) developed a two-step SVD transform method for separating upgoing and downgoing waves of zero-offset VSP data. Zhu and Wu (2010) used SVD and correlation to realize diffracted wave extraction in the local imaging matrix. Shen et al. (2009Shen et al. ( , 2016 proposed a seismic wave field separation and denoising method via SVD in a linear domain to separate P-P waves and P-S converted waves. Shen et al. (2019) also developed a method to enhance GPR diffracted waves using SVD filtering. Recently, Lin et al. (2020) developed a method based on the multichannel singular-spectrum analysis algorithm to suppress time-linear signals (reflections) and separate weaker time-nonlinear signals (diffractions) in the common-offset or poststack domain. In fact, the SVD filter is a method based on the difference in the lateral coherence between different signals to achieve seismic wave field separation and denoising. Because of the difference in travel time between diffracted and reflected waves, SVD filters should also be applied to the separation of reflected and diffracted waves. Based on this understanding, we developed a method to separate diffracted waves via SVD filters. The core of the technique is to first flatten the reflected waves through NMO to make the difference in transverse coherence between the reflected and diffracted waves clear. Then, SVD filtering is implemented to suppress the reflected waves with strong transverse coherence while extracting the diffracted waves. The extracted diffracted waves can be used for the identification and tracking of geological heterogeneous bodies (zones) after migration imaging.

SVD filter
A singular value decomposition (SVD) is a linear algebra tool that has been applied in seismic signal processing. SVD filter uses singular values as an orthogonal matrix and then orthogonally decomposes and low-rank approximation in signal space to separate signals with different lateral coherence (Freire and Ulrych 1988;Jackson et al. 1991;Vrabie et al. 2006;Shen et al. 2016Shen et al. , 2019. Let the 2D seismic data with any gather type (such as a shot gather) be X, the total number of traces be m and the number of sampling points per trace be n. Then, the SVD of the m × n matrix X can be transformed decomposed via the m × m orthonormal matrix U, the m × n orthonormal matrix A and the n × n orthonormal matrix V as follows: where U is composed of the eigenvectors of XX T , V is composed of the eigenvectors of X T X, and A is composed of singular values (Eq. 2).
where σ is a singular value labeled by i (where i = 1, 2, 3, …, r). Singular values occur on the main matrix diagonal from the largest to smallest values (σ 1 ≥ σ 2 ≥ … ≥ σ i ≥ … ≥ σ r > 0). The number of nonzero singular values is equal to the matrix rank r (r ≤ min{m, n}).
Since seismic traces can be normalized to ensure a null mean, the covariance matrix XX T of the 2D seismic data X can be represented by the correlation matrix R x (Jones and Levy 1987) as follows: Equation 3 can be further developed, The diagonal of SVD in this case is simply the variance of the signal (if it is zero mean). This is equivalent to autocorrelation function of X(t) for a delay △t = 0. If the signals between seismic traces are completely consistent in the lateral direction, then the rank (r) of the matrix X is equal to 1. And the lower the coherence of the signals between the traces in the lateral direction is, the greater the rank of the matrix X will be. X has eigenvalues only if n = m. If the square root of the eigenvalue of the matrix X is defined as the singular value of the matrix, this indicates that the strength of the seismic energy is related to the amplitude of the singular value σ i , or the sum of the eigenvalues reflects the sum of the energy E x of the seismic signals (Eq. 5).
It can be seen from Eq. 5 that the larger the component of the singular value is, the greater the power contribution in the seismic signals, and the most principal component will be the most coherent one. From the framework of lowrank approximation, seismic signals with different lateral coherence can be separated by different rank approximation. This is the basic principle of the SVD filter (Freire and Ulrych 1988;Jackson et al. 1991;Shen et al. 2016Shen et al. , 2019. Based on the above principles, assume r is greater than or equal to 3, the labeled i of the singular value (4) series σ i can be divided into three intervals ( Fig. 1), and four dividing schemes can be formed: Correspondingly, four SVD filters can be designed: low-rank approximation filter (Eq. 6), band-rank approximation filter (Eq. 7), high-rank approximation filter (Eq. 8) and band-stop-rank approximation filter (Eq. 9). The low-rank approximation filter can extract the signals with better lateral coherence. The band-rank approximation filter can separate the signals with certain coherence. The high-rank approximation filter can separate the signals with poor lateral coherence or no coherence. The band-stop-rank approximation filter can be used to filter the signal carrying a certain horizontal coherence (Shen et al. 2016(Shen et al. , 2019.
where the subscript LR indicates low-rank approximation, BR indicates band-rank approximation, HR indicates highrank approximation, BSR indicates band-stop-rank approximation, p ∈ Z, q ∈ Z and 1 ≤ p ≤ q ≤ r.

Data processing flowchart
Since reflected and diffracted waves have different moveout curves, it is possible to separate the two types of seismic signals by the SVD filters. However, the direct use of the SVD filters to separate the diffracted waves may not achieve a good filtering effect. The reason the SVD filters are good at identifying the transverse coherent signals is that admit low-rank approximations. Therefore, it is necessary to process seismic data to achieve the best horizontal coherence difference before using SVD filtering. One possible approach is to perform NMO processing to flatten the reflected waves. Reflected waves have the best transverse coherence after NMO, while diffracted waves still have poor transverse coherence due to the residual time difference. Then, band-rank or high-rank approximation filtering can effectively separate the diffracted waves after singular value spectrum analysis to determine the filter factors. The separation of reflected and diffracted waves is realized so that the filtered data are transformed into the original data domain after inverse NMO. To achieve thorough processing of the diffracted waves, we presented a complete set of diffracted wave separation and imaging processing flowcharts (Fig. 2). For comparison purposes, the simplified processing flows of conventional prestack and poststack migration imaging for full-wave field are also given in the processing flowchart. Specifically, after inputting the seismic records, the original seismic data must be preprocessed, including static correction, energy compensation, gathering and interference noise suppression. Then, velocity analysis in the common middle point (CMP) domain and NMO in the common shot point (CSP) domain is performed. Next, singular value spectrum analysis is performed on the shot records after NMO to obtain the SVD filtering factors, and then, a band-rank or high-rank approximation filtering is implemented to extract the diffraction waves. The extracted diffracted waves may contain waveform distortion at large Fig. 2 Flowchart of diffraction wave separation and imaging. For comparison purposes, the simplified processing flows of conventional prestack and poststack migration imaging for full-wave field are also given in the processing flowchart offsets caused by NMO, so it is necessary to remove these distorted noises. The diffraction seismic records can be obtained after inverse NMO. The seismic imaging sections for geological interpretation can be obtained after diffraction prestack migration imaging, full-wave field poststack imaging and prestack migration imaging.
It should be emphasized that separation of diffracted waves via SVD filtering can be carried out in the common shot point (CSP) domain, common receive point (CRP) domain or common midpoint (CMP) domain. The processing flow for the separation of diffracted waves in different processing domains is basically the same. The only difference is that the gathering is required before SVD filtering. The separation of diffracted waves in any processing domain must satisfy the condition that a survey line with many gathers can be processed continuously by one gather after another, that is, the number of channels in each gather must be the same. We selected the common shot point (CSP) domain to separate diffracted waves in this paper. In addition, a seismic line with multiple gathers can be processed continuously with the same filter parameters when the lateral fluctuations of the stratum interface are small. Otherwise, it is necessary to take measures of segmentation processing, that is, different filtering parameters are used to process the gathers of different formation interface characteristic sections.

Geological model and seismic wave forward modeling
To assist the test of this method, a geological model with one thin-out, two caves and four reverse faults was generated. The identification of the thin-out is A, the caves are P1 and P2, and the faults are F1, F2, F3 and F4. The geological model parameters are shown in Fig. 3, the thin-out and cave elements are shown in Table 1, and the fault elements are shown in Table 2. The finite difference decomposition wave equation technique was used to simulate seismic waves. The Ricker wavelet was selected as the source wavelet with a main frequency of 60 Hz. An observation system with a fixed geophone and a moved shot point was adopted. The total number of geophones was 241, the trace spacing was 5 m, the shot spacing was 10 m, the sampling rate was 0.5 ms, the recording length was 0.7 s, and the total number of shots was 121. Figure 4 shows the 1st (Fig. 4a), 60th (Fig. 4b) and 121st shot (Fig. 4c) seismic records. It can be seen from the original seismic records that the diffracted waves caused by the pinch-out, caves and faults are clearly visible, they develop between 200 and 700 ms, and the diffracted wave energy is also weakened as the size of the faults or tip points decreases.

Diffracted wave separation
To test the robustness of wave field separation, Gaussian random noise was added to the seismic records (Fig. 5a), and the SNR was approximately 8.0. SVD filtering was directly   applied to the seismic data containing noise. It can be found from the obtained singular value spectrum (Fig. 5b) that the energy of the singular value decays rapidly, there are three σ decay characteristic curves, and two inflection points can be defined that are i = 24 and i = 60. We selected one of the inflection points (i = 24) as the demarcation point of the filtering factor according to the SVD filtering principle. It can be seen from the filtering results that the low-rank approximation filtering (Fig. 6a) and high-rank approximation filtering (Fig. 6b) not only fail to clearly separate the diffracted waves or reflected waves, but also seriously damage the seismic signals. The failure of filtering can be attributed to the fact that the transverse coherence difference between the different signals in the original seismic records is not obvious. The seismic data were processed again according to the diffracted wave processing flowchart (Fig. 2). After velocity analysis, NMO was carried out for the noise-containing seismic records. It can be seen from the NMO results (Fig. 7a) that the reflected waves have been flattened, the diffraction  waves still have a residual time difference, and the waveform distortion at the large offset is caused by NMO stretch. Compared with the singular value spectrum before NMO (Fig. 5b), the spectrum decays more rapidly after NMO correction indicating the system admits a low-rank or bandrank approximation filtering (Fig. 7b). There are also three σ decay characteristic curves, and the inflection points are defined at i = 12 and i = 24. According to the distribution law of different seismic signals in SVD domain, it can be judged that the signals before inflection point 1 are represented reflected waves, the signals between inflection point 1 and inflection point 2 are expressed as diffracted waves, and the signals after inflection point 2 are mainly denoted as irregular noise and partial diffracted waves. In general, inflection point method is the basic method to determine the filtering parameters. However, when the seismic data is complex, such as low SNR or weak energy of diffracted waves, it is not necessarily able to accurately determine where the inflection point is. In this case, we choose different filter parameters to determine the final filter parameters through multiple filter tests.
Based on the above key information, we tested the filter parameters. First, i = 12 was used as the cut-off point (Fig. 8a) to perform high-rank approximation filtering in Fig. 7a. The filtering results in Fig. 8b show that the diffracted waves were separated while the random noise has been also separated. Then, i = 18 was used as the cut-off point (Fig. 8c) to perform high-rank approximation filtering in Fig. 7a. The filtering results of Fig. 8d show that part of the diffracted waves has been lost. Finally, i = 200 was used as the cut-off point (Fig. 8e) to perform high-rank approximation filtering in Fig. 7a. The filtering results of Fig. 8f show that only the noise has been separated. Based on the filtering test results (Fig. 8), we finally determined that the band-rank approximation filtering parameters were 12 ≤ i ≤ 200 (Fig. 9a). It can be seen from the final filtering results (Fig. 9b) that the diffracted waves were completely extracted with a band-rank approximation filtering, while the strong reflected wave interference was effectively suppressed. After removing the waveform distortion caused by NMO and inverse NMO, the diffracted wave seismic  Fig. 5a. a Low-rank approximation (i < 24) and b high-rank approximation (i ≥ 24) results. The filtering results are unable to separate the diffracted waves  Fig. 7a. a High-rank approximation filtering spectrum with i ≥ 12 and corresponding filtering results b, c high-rank approximation filtering spectrum with i ≥ 18 and corresponding filtering results d, e high-rank approximation filtering spectrum with i ≥ 200 and corresponding filtering results f 1 3 records shown in Fig. 9c were obtained. Figure 9d shows the obtained seismic records of the reflected waves after a low-rank approximation filtering (i < 12) and inverse NMO (waveform distortion removal cased by NMO). It can be seen by comparing Fig. 9c with Fig. 9d that the diffracted waves were fully extracted, while the reflected waves were effectively suppressed.

Imaging results comparison and analysis
The prestack migration imaging of the separated diffracted waves (Fig. 9c) was performed using the Kirchhoff migration method and compared with the results of the full-wave field imaging results (including stacked and prestack migration). The imaging results show that the seismic response characteristics of the thin-out, caves and faults can be identified in the stacked section (Fig. 10a). However, the diffracted waves are not focused, the vertical and horizontal resolutions are low, and it is difficult to accurately determine the specific location and scale of the thin-out, caves and faults. Compared with the stacked section, the imaging quality was significantly improved on the prestack migration imaging (Fig. 10b). The diffracted waves caused by the thin-out (A), caves (P1 and P2) and faults (F1, F2 and F3) were effectively focused, and the vertical and horizontal resolutions were also significantly improved. However, only the P1 and P2 caves, the F1, F2 and F3 faults can be identified, while the thin-out (A) and the small-scale F4 fault cannot be identified. The reason is mainly attributed to the weak seismic response caused by the small size of thin-out (A) and fault F4. Moreover, the seismic response is seriously affected by reflected waves with strong energy. The migration imaging results of the extracted diffracted waves (Fig. 10c) are significantly better than the full-wave field imaging sections. The seismic responses of the thin-out, caves and faults have been effectively imaged. Especially the thin-out (A) and the F4 fault can be clearly identified, and the imaging positions accurately correspond to the thin-out, cave and fault positions in the geological model. The processed results of the model seismic dataset show that the proposed diffracted wave separation method is effective. And the results of diffracted wave imaging have an absolute advantage in the identification and tracking of small-scale faults, pinch-out, karsts and other heterogeneous bodies (zones).

Diffracted wave separation
We further applied the method to process a field seismic dataset that was acquired using multicoverage rolling acquisition observation system. The number of receiving channels was 96, the channel spacing was 5 m, the shot spacing was 5 m, the sampling rate was 1.0 ms, the recording length was 1.0 s, and the total number of shots was 121. Figure 11a shows typical original single-shot seismic data (Shot No. 86). The surface waves, which are the main source of interference noise, are more developed with stronger energy. The reflected and diffracted waves hidden under the surface waves are highlighted after denoising (Fig. 11b). Figure 11c shows the singular value spectrum of a shot of seismic records (Shot No. 86). Singular value spectral analysis showed that the band-rank approximation filtering factors were 10 ≤ i ≤ 85. The diffracted waves shown in Fig. 11d were obtained by using the data processing flowchart in Fig. 2, in which the diffracted waves were  Figure 12 shows the imaging results of the field seismic dataset. It can be seen from the stacked section of the fullwave field (Fig. 12a) that the structure of the section is relatively complex. There are two thin-outs between CMP200 and CMP500 with several small interstratal faults. There is an unconformity plane on the right of CMP400. Due to the difference in erosion levels, several diffraction events with weak energy are developed on the seismic section that are almost submerged in the reflected waves with strong energy. On the right side of the section, there is a group of diffracted waves with a wide spatial distribution range (from CMP600 to CMP1000) and a long depth span (from 100 ms to 800 ms). It is difficult to determine what causes the generation of diffracted waves only from the stacked section. The diffracted waves have been effectively converged in the poststack migration section of the full-wave field (Fig. 12b). The diffracted wave group on the right of the section almost converged into a seismic wave anomaly area. Combined with the events of the fault on both sides of the abnormal area, it can be preliminarily judged that a large reverse fault developed in this area, which caused the formation to break and form a heterogeneous zone. On the prestack migration section of the full-wave field (Fig. 12c), both the vertical and horizontal resolutions are significantly higher than those of the poststack migration results. However, due to the influence of Fig. 11 Diffracted wave separation results of field seismic dataset. a The original single-shot seismic records (Shot No. 86), b the seismic records after denoising, c the band-rank approximation filtering spectrum and d the extracted diffracted waves. The diffracted waves have been enhanced, while the reflected waves have been effectively suppressed strongly reflected energy, the thin-outs, breakpoints and convex points of unconformity are almost completely covered under the reflected waves. On the prestack migration section of the diffraction wave (Fig. 12d), the thin-outs, breakpoints and convex points are clearly imaged. Based on the disordered seismic response characteristics, three geological heterogeneous zones can be circled and labeled as I, II and III. Zone I clearly shows the information of the thin-outs and convex points caused by the weathering crust. Zone II shows the fault and the fracture heterogeneity caused by the fault. Zone III is a heterogeneous area that may be caused by stratigraphic heterogeneity.

Conclusions
The complete separation of diffracted and reflected waves is an important basis for the effective utilization of diffracted waves. Based on the travel time difference between diffracted and reflected waves and the characteristics of SVD filtering, we developed a diffracted wave separation method. The conclusions and understandings are as follows: 1. The amplitude of the singular component determines its power contribution in 2D seismic data. From the framework of low-rank approximation, seismic signals with different transverse coherence characteristics can be separated by different rank approximation.
2. The time distance characteristics of diffracted and reflected waves are inconsistent. After NMO, the reflected waves can be flattened, but the transverse coherence of the diffraction waves is still poor due to the existence of the residual time difference. Then, the diffracted waves can be effectively extracted by a bandrank or high-rank approximation filtering. 3. The extracted diffracted waves can be used for prestack migration imaging after inverse NMO. If the imaging results are interpreted in combination with the reflected wave imaging results, the accuracy of identification and tracking of geological heterogeneous bodies (zones) can be further improved.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.

Fig. 12
Comparison of field seismic data imaging results. a The stacked section, b poststack migration section and c prestack migration section of the full-wave field. d The prestack migration section of the extracted diffracted waves. The thin-outs, breakpoints and convex points have been clearly imaged, and three geological heterogeneous zones labeled I, II and III are clearly presented