Iterative deblending using the POCS algorithm in the approximate flattened domain

We proposed an improved method to eliminate the interference generated by simultaneous-source acquisition, which can help shorten the acquisition period and improve the quality of seismic data. An iterative mathematical framework is devised, which uses the projection onto convex sets algorithm to estimate the blending noise subtracted from the pseudo-deblended data to separate the blended data in an iterative way. Differently to the conventional method using the coherent-promoting operator only based on the curvelet transform, we combine the curvelet transform and the approximate flattened operator (AFO) to improve the deblended result, which can flatten seismic events approximately to preserve the details of useful signals. This is the first time that the AFO and the curvelet transform are combined to enhance the effect of the coherent-promoting operator and improve the performance of deblending. To display the advantages of the improved method, we use both simulated synthetic data and field data examples to compare and analyse the deblended results using our method and the conventional method, and confirm that the improved method can perform better.


Introduction
Simultaneous-source techniques have attracted significant attention from researchers because of their economic benefits and technical challenges (for example, Beasley et al. 1998;Berkhout 2008;Mahdad 2012;Wu 2014). However, the benefits are compromised by the intense interference between different shots, which makes the acquired seismic data difficult to process in a conventional seismic processing workflow. Without a careful design of the processing workflow for the simultaneous-source data, many subsurface seismic tasks, e.g. surface-related multiple attenuation (SRME) (Kelamis & Verschuur 2000), reverse time migration (Dai et al. 2011), full waveform inversion (Wang & Rao 2009), and hydro-carbon reservoir characterization (Liu et al. 2017), cannot be performed appropriately. Generally, there are two ways to deal with this issue. One way is to use a two-step strategy, first separating (which is also referred to as deblending) and then conventional processing (for example, Beasley et al. 2012;Wapenaar et al. 2012;Li et al. 2013;Ibrahim & Sacchi 2014;Berkhout & Blacquire 2014). The other approach is to image directly (for example, Dai et al. 2011;Verschuur & Berkhout 2011;Guitton & Daz 2012;Jiang & Abma 2010;Soni & Verschuur 2015), which strongly depends on an fairly acceptable initial macro velocity model. Currently, deblending is the dominant way to remove the interference noise and the result is better. Therefore, we mainly focus on deblending of simultaneous sources in this paper.
Different filtering and inversion methods have been used to deblend blended records. Filtering approaches can be applied to separate the simultaneous sources according to the assumptions that useful signal is coherent and blending noise is incoherent in some domains other than the common shot domain. Researchers have proposed many different filtering methods to separate the interference noise (for example, Zhang et al. 2015), which can also be referred to as a noise attenuation problem. Different types of noise attenuation methods, e.g. the thresholding based methods , rank-reduction based methods (Zhou & Zhang 2017), stacking based methods (Neelamani et al. 2006), decomposition based methods (Han & van der Baan 2015), mathematical morphology based methods (Duan et al. 2010), and prediction based methods (Liu & Chen 2013), can be applied in a straightforward manner to suppress the blending noise. Kim et al. (2009) used a simple approach to model the interference from the other source in the common-offset domain. Mahdad et al. (2011) separated the simultaneous source in the common-receiver domain using the iterative F-K filter. Huo et al. (2012) attenuated the blending noise in the common midpoint 1105 Figure 2. Diagram of reconstruction error for the synthetic example in the curvelet domain. The red line corresponds to synthetic data with large curvature, and the blue line corresponds to synthetic data with small curvature. domain using a median filter. Yang & Ying (2014) developed a more powerful curvelet transform filtering method and utilized it to improve the deblending performance. Hampson et al. (2008) applied a simple dip filter to remove the blending noise according to the slope differences between simultaneous sources in the common shot domain.
The inversion methods treat the deblending processing as an estimate problem (for example, Doulgeris et al. 2012;Wason et al. 2014;Kumar et al. 2015;Cheng & Sacchi 2016), which aim at estimating the ideal unknown unblended data by attenuating the blending noise. Because of the ill-posed property of inversion problems, some constraints should be added into the inversion framework to regularize the inversion problem. Akerberg et al. (2008) used the sparse constraints in the Radon domain to regularize the inversion. Abma et al. (2010) proposed using F-K domain sparsity as a constraint. Mahdad et al. (2011) proposed a coherencebased inversion method to deblend the simultaneous-source data. Chen et al. (2017) proposed a general iterative deblending framework via curvelet-domain shaping regularization. Cheng & Sacchi (2013) developed a rank-reduction algorithm based on singular spectrum analysis (SSA) that is capable of suppressing the interference noise generated by simultaneous-source acquisition. There are also some other researchers improving the rank-reduction algorithm to obtain superior denoising performance. Li et al. (2014) used iterative deblending of the simultaneous-source method to perform a field trial. Zhou (2017) proposed an inversion with multiple convex constraints to improve the deblending performance based on the projection onto convex sets (POCS) iterative framework.  proposed an approach using a wavelet transform to deterministically separate the primary signal from the noise, including simultaneous-source interference. Chen et al. (2017) used a sparse constraint in the curvelet transform to remove the interference noise and estimate the unblended seismic data in an iterative way. Furthermore, the inversion result is better than that using an iterative seislet thresholding algorithm. We follow the sparse constraint in the curvelet transform and then propose a novel deblending approach in this paper, in which we modified the coherentpromoting operator by adding an operator called the approximate flattened operator (AFO). It can enhance the sparsity of the curvelet coefficients of useful signals, making the inversion result better by using AFO to approximately flatten seismic events based on the theory of normal moveout (NMO) correction. Subsequently, we use the novel deblending approach to separate the simultaneous source in this paper. Both synthetic data and numerically blended field data examples demonstrate that the proposed method performs better than the conventional method.

Deblending via then POCS algorithm
To be convenient, we assume that one simultaneous source of data is blended using two independent sources, which correspond to two shooting sources in ocean bottom cable (OBC) acquisition. Therefore, the blending forward processing is given by in the common-receiver domain, where d is the blended data, d 1 and d 2 are the first and second independent gathers, respectively, and Γ 1 and Γ 2 denote the dithering operators corresponding to d 1 and d 2 , respectively. In the process of deblending, different phase-encoding blending operators have been introduced to eliminate the crosstalk noise, such as a linear phase term and a quadratic phase term. Here, the linear phase encoding that expresses the time delays in the time domain was applied (e.g. Mueller et al. 2015). The blending operator can be written in the following form, where Γ i is the dithering operator of the ith common-receiver record in the time domain, P i represents the phase matrix of the ith common-receiver record, t i, m is the delay time for the mth trace of the ith common-receiver record, w is the angular frequency, and j represents the imaginary unit. Note that we use random time delay in this paper.
To retrieve individual deblended shot records from blended data, we take the inverse of Γ i (i = 1, 2) to equation (1) and obtain the pseudo-deblending formulas, which indicate that the record of the main source (such as d 1 or d 2 ) is coherent, and the other term in equations (4) or (5) is incoherent in the pseudo-deblended data. Then we exchange the position of the unblended data (such as d 1 and d 2 ) and the pseudo-deblended data (such as Γ −1 1 d and Γ −1 2 d) in equations (4) and (5), and obtain two new equations as follows, where the pseudo-deblended data are not equal to the unblended data because of the extra interference noise (such as . Therefore, if we can predict the interference noise, the unblended data can be obtained by subtracting the interference noise from the pseudo-deblended data. Even though the unblended data cannot be known, we can estimate them in an iterative way. Accordingly, we provide the following general mathematical framework for blended multi-common-receiver records, where D ' n+1 is the estimated unblended data at the nth iteration, D is the blended data, Γ -1 is the time dithering operator, Γ -1 D represents the pseudo-deblended data, and D ' represents the predicted interference noise at the nth iteration. For equation (8), we can use a POCS algorithm (for example, Bauschke & Borwein 1996) to produce significantly improved results, which can provide a coherent-promoting filter based on the assumptions that signal is coherent and the interference noise is incoherent. Then equation (8) can be reformulated as follows, where S is a coherent-promoting operator, which can be chosen as multi-dimensional Fourier transform (F-K transform) Mahdad et al. 2011), Curvelet transform (Kontakis & Verschuur 2017), shearlet transform or learning-based sparse dictionary (Chang et al. 2016) and so on. The F-K transform is simple and fast, but its result is not very good. The seislet transform needs to estimate the local slopes of the blended data. It is generally difficult to obtain accurate local slopes for complicated data, and the computation is usually large. The learning-based sparse dictionary has a relatively higher computational cost. Even though we can increase the computational efficiency, the performance of denoising will be degraded. The curvelet transform directly applies the characteristics of multi-scale and sparsity to provide a better coherent-promoting operator, regardless of the local slopes, and the calculation amount is moderate. Consequently, we choose the curvelet transform to obtain 1108 Figure 5. Source and receiver geometry for the simulated synthetic data. The green and red dots denote the two sources, the green and red arrows indicate the directions of their travel, and the blue triangle denotes the receiver. Two corresponding shots in the two sources shoot nearly simultaneously, and the time schedule ranges from -400 to 400 ms, which can generate strong interference in the recorded common-receiver gathers.
better deblended results in this paper. However, according to many numerical tests, we find that the operator consisting of the curvelet transform cannot preserve the details of the inversion results well, in particular for seismic data with large curvature. In order to solve this issue, we analyse the sparsity of the curvelet coefficients of events with large and small curvature.
We first generate two simple synthetic datasets with large curvature and small curvature (figures 1a and 1b). The data with small curvature can be obtained using the NMO method to correct the data with large curvature. Figure 2 shows the reconstruction error generated in the process of using different percentages of the coefficients to reconstruct the synthetic data, which can illustrate the sparsity of the synthetic data in the curvelet domain. It is obvious that the data with small curvature have the smaller reconstruction error than those with large curvature. Therefore, the data with small curvature are sparser in the curvelet domain, which only need fewer rectangular windows to approximate the data if considering one event, because the base function of the discrete curvelet transform is an anisotropic rectangle. Furthermore, it can also be understood that the curvelet coefficients of the data with small curvature just exist in fewer directions, as shown in figure 1d, than those of the data with large curvature, as shown in figure 1c, in the curvelet domain, and at the extreme if the data are horizontal, the curvelet coefficients are distributed only in one direction, which can make the distributions of curvelet coefficients between the useful signal and noise more obvious if considering the noise.
In order to illustrate this characteristic, we synthesize the blended data, analyse the distributions of curvelet coefficients and obtain the inversion results. Figures 3a and 3b show the blended data corresponding to figures 1a and 1b, respectively. Figures 3c and 3d show the curvelet coefficients of blended data. It can be found that the distinctions between the curvelet coefficients of useful signal and blending noise are more apparent in figure 3d than those in figure 3c. Then we use the conventional inversion method to gain the inversion results as shown in figures 4a and 4b. Comparing figures 4a and 4b, the result is better for blended data with small curvature, in which the blending noise can be attenuated and simultaneously the details of signal can be preserved well in the process of inversion. Figures 4c and 4d show the corresponding noise sections, achieved by computing the differences between the deblended sections and the blended sections. There are obvious useful signals leaking out in figure 4c. However, there are not useful signals left in figure 4d. In addition, we also obtain the estimation error sections by computing the differences between the deblended sections and the unblended sections, shown in figures 4e and 4f. It can be seen that the estimation error for figure 4a is comparatively large, because of the existence of large curvature events.
From the above analysis, we can find that the event with small curvature can be deblended better than that with large curvature. Accordingly, in order to reach better inversion results, especially for the blended data with large curvature, we use an operator to transform the seismic data into a domain to approximately flatten the seismic data and reduce the curvature of the events. Then the coherent-promoting operator S can be written as, where A and A -1 indicate the forward and inverse AFOs, C and C -1 indicate the forward and inverse curvelet transforms, and the T represents a thresholding operator with an input parameter in the transform domain. Because of the differences between the curvelet coefficients of useful signals and interference noise in the curvelet domain, we can use a thresholding operator to achieve better separation between the coherent events and less coherent events. The thresholding operator T can be either a soft-thresholding operator or a hard-thresholding operator. In this paper, we choose the hard-thresholding operator (e.g. Davies 2008, 2009;Gao et al. 2013), which can achieve a higher signal-to-noise ratio (SNR) than the soft-thresholding function according to many numerical tests.

Approximate flattened operator (AFO)
Some researchers have studied the methods to flatten seismic data, such as the prediction painting method. The prediction painting method uses the local slope estimated by the plane-wave destruction (PWD) method to obtain a plane-wave prediction operator as a flattened operator. However, it usually generates very large errors in the process of forward and inverse flattened transform, because of the existence of interference noise affecting the accuracy 1109 Figure 6. Simulated synthetic data test. (a) The blended data, (b) the transformed blended data obtained by transforming the blended data into the approximate flattened domain, (c) the curvelet coefficients corresponding to the unblended data, which synthesize the blended data corresponding to figures 6a, 6d the curvelet coefficients corresponding to the transformed data generated by transforming the unblended data to the approximate flattened domain, (e) the curvelet coefficients corresponding to figures 6a, and 6f the curvelet coefficients corresponding to figure 6b.
1110 Downloaded from https://academic.oup.com/jge/article-abstract/15/4/1104/5710220 by guest on 21 January 2020 of the local slope. In addition, our method does not have to flatten seismic events. Therefore, we provide an AFO based on the theory of NMO regardless of the local slope, which can be applied in the iterative framework to preserve the details of inversion results, and the operator is easily implemented. The NMO can be used to correct the reflection time of different offsets to zero-offset time, which corresponds to each two-way zero-offset traveltime to flatten seismic events. We can model the traveltime t p at a source-receiver distance x as the following equation (11), assuming that the near-offset seismic events follow a hyperbolic shape, where t 0 is the traveltime at zero offset, and v p is the velocity. Then the moveout can be calculated as follows, where v p = v nmo , and Δ t represents the moveout. Then bringing equation (12) into equation (11), we can acquire the result after applying NMO correction, From equation (13), it can be found that the common midpoint seismic records can be flattened after precise NMO correction with a precise NMO velocity, obtained using the scanning method, which is time-consuming and easily generates an NMO stretch phenomenon for large-offset gathers. Consequently, we employ the theory of NMO correction to approximately flatten seismic records, avoiding the need for a precise velocity. We just use several constant velocities to determine the approximate flattened velocity with the criterion of maximum auto-correlation in the phenomenon of no stretching. Then the AFO can be obtained, to be applied in operator S (equation (10)) to improve the inversion result. Note that the AFO can be used not only in the common midpoint gather, but also in other gathers, except for the common shot gather.

Numerically blended synthetic data
In this section, we first synthesize a complex numerical example to test the effectiveness of our proposed method, which is generated by convolving reflectivity with a Ricker wavelet and then blending via the time dithering ranging from -400 to 400 ms. The source and receiver geometry is shown in figure 5. There are two sources traveling in opposite directions, shooting nearly simultaneously, and recorded in the common receiver (the blue triangle). Because of the symmetric geometry, the two common-receiver gathers are the same. Therefore, we only need to display the inversion result for one common-receiver gather. Figure 6a shows the synthetic blended data, which can be transformed into the approximate flattened domain according to equation (13), as shown in figure 6b. The velocity of transforming is chosen as 2500 m/s. Figures 6e and 6f show the curvelet coefficients corresponding to figures 6a and 6b, respectively. Figures 6c and  6d show the curvelet coefficients corresponding to the unblended data (composing the blended data shown in figure  6a) and the transformed data obtained by transforming the unblended data into the approximate flattened domain, respectively. As shown in figure 6c, the distributions of curvelet coefficients of the useful signals are not concentrated, so that some of the curvelet coefficients in figure 6e cannot be judged as to whether they are useful, because the energy differences between the curvelet coefficients of the useful signals and blending noise are not obvious and are mixed together on the same scale. Different to figure 6c, the curvelet coefficients of useful signals shown in figure 6d are more concentrated in some directions, which can help suppress the curvelet coefficients of blending noise in figure 6f. Figure 7 shows the reconstruction errors for the unblended data and the transformed data. It shows that the transformed data are sparser than the unblended data in the curvelet domain, further explaining the reason why the distributions of the curvelet coefficients of useful signals in figure 6d are more concentrated than those in figure 6c. Then if we use the same thresholding to attenuate the curvelet coefficients of blending noise in the curvelet domain, compared with figure 6e, the result of figure 6f will be more favorable and protect the curvelet coefficients of useful signals better. In other words, the blending noise in the 1111 Figure 8. Simulated synthetic data test. (a) The inversion result using our proposed method, (b) the inversion result using the conventional method, (c) the blending noise using our proposed method, (d) the blending noise using the conventional method, (e) the estimation error using our proposed method, and (f) the estimation error using the conventional method. approximate flattened domain can easily be muted out than that before the transformation with the same thresholding.
We use our improved method and the conventional method to separate the blended data. Figures 8a and 8b show the deblended results using the improved method and the conventional method, respectively. It is obvious that the deblended result using our improved method (figure 8a) is cleaner without any evident blending noise and all useful events are recovered better, compared with using the conventional method (figure 8b). The corresponding noise sections are shown in figures 8c and 8d. There are hardly any leakages of useful energy (actually hidden in the noise), as demonstrated from the blending noise sections. Figures 8e  and 8f show the corresponding estimation error sections. In figure 8e, our improved method causes a smaller estimation error, meaning an acceptable inversion result. However, the estimation error for the conventional method is comparatively large in figure 8f, especially for the components with large dip, because the energy differences between the curvelet coefficients of useful signals and blending noise are so smaller that we cannot separate them well and damage the useful signals. Furthermore, the zoomed sections corresponding to figures 8a and 8b are shown in figures 9a and 9b, respectively, in which it is easy to observe the effect of deblending. There is no doubt that the deblending result using our improved method outperforms the conventional method, with a clearer event. We also demonstrate the two single-trace comparisons, shown in figures 10a and 10b,  where the red line corresponds to the improved method, the blue line corresponds to the conventional method, and the black dotted line corresponds to the unblended data. On the whole, these three lines have the same trend. However, the red line is mostly overlapped with the black dotted line. This indicates that the improved method can preserve the details better.
We use the signal-to-noise ratio (SNR) as the evaluation criterion to analyse the deblending performance quantitatively, 1113 Figure 11. Diagram of SNR for the synthetic example. The blue line corresponds to our improved method, and the red line corresponds to the conventional method.

Figure 12.
Source and receiver geometry for the simulated field data test. The green and red dots denote the two sources, the green and red arrows indicate the directions of their travel, and the blue triangle denotes the receiver. Two corresponding shots in the two sources shoot nearly simultaneously, and the time schedule ranges from -400 to 400 ms, which can generate strong interference in the recorded common-receiver gathers. The differences between figure 5 and figure 12 are the locations of the selected receiver and the travel directions of the two sources.
where d n denotes the deblended data after n iterations, d denotes the unblended data, and ∥⋅∥ 2 2 denotes the squared L 2 norm of a function. Figure 11 shows the diagram of SNR corresponding to the improved method and the conventional method. It demonstrates a superior performance of our improved method in comparison with the conventional method.

Numerically blended field data
In this section, we use field data to test the proposed method. The acquisition geometry is shown in figure 12. Different to figure 5, the receiver is put on one side of the survey area in the middle of the two source lines, and the travel directions of the two sources are the same. We can acquire data with the same shapes from this acquisition geometry. Therefore, we can only display the deblended result for one commonreceiver gather, same as the test of blended synthetic data. Figure 13a shows the blended field data, whose time dithering code ranges from -400 to 400 ms. Figure 13b shows the transformed blended data generated by transforming the blended data (figure 13a) into the approximate flattened domain according to equation (13). The velocity of transforming is chosen as 2500 m/s. Figures 13e and 13f show the curvelet coefficients corresponding to figures 13a and 13b. Figure 13c shows the curvelet coefficients of the unblended data, which simulates the blended data shown in figure 13a. Then we transform the unblended data into the approximate flattened domain, called the transformed data, whose curvelet coefficients are shown in figure 13d. Even though the field data are relatively more complicated, we can still find that the energy distributions of curvelet coefficients in figure 13d are more concentrated, compared to those in figure 13c. Consequently, the curvelet coefficients (figure 13f) of blending noise can be attenuated better than those in figure 13e with the same thresholding, because of the obvious energy distribution differences between the curvelet coefficients of the useful signals and blending noise in figure 13f. Here, we do not display the diagram of the reconstruction errors of the unblended data and the transformed data, because the field data contains some noise in addition to the useful signals, due to which we cannot use reconstructive errors to represent the sparsity of useful signals.
The deblended results using our improved method and the conventional method are shown in figures 14a and 14b, respectively, from which we can find that our improved method performs better than the conventional method. There is a lot of residual blending noise in the deblended data as shown in figure 14b. However, the deblended result in figure 14a is relatively cleaner with less residual noise. Then we compute the corresponding noise sections as shown in figures 14c and 14d to describe the effect of deblending. Comparing figures 14c and 14d, we can find that the differences between them are not distinct, except for the boundary of the data. In order to explain the effect of our improved method more intuitively, we determine the corresponding estimation error sections shown in figures 14e and 14f. Even though the field data contain some noise, causing the estimation error sections to not be equal to zero even with perfect deblending results, we can evaluate the effect of deblending through residual information. It is obvious that the residual information (figure 14e) is less than that in figure 14f. In addition, we show the zoomed sections  Figure 14. Simulated field data test. (a) The inversion result using our proposed method, (b) the inversion result using the conventional method, (c) the blending noise using our proposed method, (d) the blending noise using the conventional method, (e) the estimation error using our proposed method, and (f) the estimation error using the conventional method.
in figures 15a and 15b corresponding to the two deblended datasets (figures 14b and 14b). Here it is more obvious that the result using our improved method is superior to that using the conventional method, with less residual noise. We also display the comparison of two single-traces in figures 16a and 16b. Because of the complexity of real single-source data containing some noise, the red line cannot be exactly equal to the black dotted line. However, the red line is closer to the black dotted line than the blue line, which may represent the advantages of the improved method. Moreover, the diagram of the SNR of our improved method also demonstrates superior behavior, as shown in figure 17. The above results reveal that our improved method excels in the conventional method.

Discussion
The mathematical framework used in our improved method is different from the general inversion framework, and is devised according to the idea of estimating the blending noise directly and then removing it from the pseudo-deblended data in an iterative way. This is in fact a kind of POCS method often used in interpolation problems. In our study, we do not test which mathematical framework is better. We just improve the coherent-promoting operator by adding the AFO to obtain better deblended results. The AFO is based on the theory of NMO with a constant rough velocity avoiding time-consumption, and just needs several tests, because our method does not require keeping the events absolutely  horizontal. Accordingly, the improved method does not cost more time than the conventional method and can be implemented in a straightforward manner.
In this paper, we first adopt a simple synthetic example to analyse the reason why the deblended result of the data with small curvature is better than that for the data with large cur-vature. We arrive at an observation that the curvelet coefficients of the data with small curvature are sparser than those for the data with large curvature in the curvelet domain. Note that this observation is just focused at before and after the approximate flattened transformation of data, not any specific forms of data. In addition, the theory of discrete curvelet 1117 Figure 17. Diagram of SNR for the synthetic field data example. The blue line corresponds to our improved method and the red line corresponds to the conventional method.
transform can also explain this. Because the basis function of the discrete curvelet transform is an anisotropic rectangle, it can therefore use fewer rectangle windows to describe the data with small curvature relative to those with large curvature, representing the characteristic of sparsity. Based on the above analysis, we proposed an improved method and then use simulated complex synthetic data and field data examples to test it. The deblended results, blending noise sections, estimation error sections, single-trace comparisons and the diagrams of SNR illustrate the superior performance of our method in contrast to the conventional method. Even though the simulated field data example is not a real test, the field data are relatively complex and have been processed by other researchers, enhancing the credibility of our improved method when applied to real-world blended acquisition.

Conclusions
We have proposed an improved inversion algorithm to suppress the blending noise from simultaneous sources in an iterative way. At each iteration, the blending noise can be estimated directly by applying the coherent-promoting operator to shape the estimated unblended data of the previous iteration in the equation of blending noise, and be removed from the pseudo-deblended data by simple subtraction. The coherent-promoting operator consists of the curvelet transform and the AFO based on the theory of NMO, which can preserve the details of useful signals and help to recover the unblended data better in the process of inversion. Because of the characteristic of the coherent-promoting operator, our improved method is superior to the conventional method, in which the coherent-promoting operator just consists of the curvelet transform. Finally, we have used both simulated synthetic data and field data tests to show that our method performs better.