Iterative Deblending of off-the-Grid Simultaneous Source Data

Simultaneous source acquisition can enhance the seismic data quality or improve the field acquisition efficiency. However, one of the disadvantages is that the simultaneous source data are often obtained on a non-uniform sampled grid in realistic acquisition. Except for the blending noise introduced by the temporal overlap, the non-uniform samplings usually cause serious artifacts which also greatly reduce the data quality and imaging resolution, such that data regularization must be implemented during the deblending. At present, most deblending algorithms are suitable for the uniformly sampled data. Thus, we propose the nonequispaced fast discrete curvelet transform-based deblending approach to deal with the non-uniformly sampled simultaneous source data. In order to enhance the deblending efficiency, we introduce the iterative shrinkage thresholding algorithm to solve the undetermined problem via a soft thresholding operator. Once the coefficients in the curvelet domain are acquired, uniformly sampled deblended data can be obtained via inverse curvelet transform. During the iterations, the threshold is decreased from a large value to a small value according to the new exponential function. Compared with the existing deblending methods that are applied on uniform grids, the key contribution is that the proposed method can directly work on the non-uniformly sampled simultaneous source data. The numerical examples are used to demonstrate the successful performance of the proposed method in attenuating blending noise and correcting the distorted events.


I. INTRODUCTION
Simultaneous source acquisition, or blended acquisition, which can save the acquisition cost and increase the data quality, has attracted significant attention from both academia and industry [1]- [3]. However, the blending noise caused by simultaneous source acquisition will reduce the benefits [4]- [7]. To eliminate the adverse impact of the strong blending noise on those subsequent processing tasks, a main approach for processing simultaneous source data is to separate the blended records first rather than direct imaging [8]- [10], although some direct migration or migration velocity analysis algorithms for the blended data are proposed [11]- [13]. Then separated seismic gathers can be used in a conventional seismic processing workflow. At present, there are many reported effective deblending methods in simultaneous sources separation (or deblending) The associate editor coordinating the review of this manuscript and approving it for publication was Abdel-Hamid Soliman . either on theoretical or field deblending tests. All of these methods can be almost classified into two categories of deblending approaches that are carried out in some sorted gathers as follows: 1) treating it as a noise suppression problem [14]- [16]. According to the assumptions that useful signal is coherent and blending noise is incoherent in some sorted gathers, the blending noise can be removed by the noise filtering method, such as the median filtering and others [17], [18]. Filtering-based deblending method is efficient and easy to implement. However, it cannot achieve a good performance when the blended data is very complex. 2) treating it as an inversion problem [19]- [23]. This method aims at estimating the unknown unblended data from the blended data. Because of the ill-posed property of the deblending problem, inversion-based methods require some constraints, such as sparsity constraint and compressible constraint [24]- [27], and a regularization term often need to be introduced to ensure the uniqueness and stability of solution [28]. These two categories of deblending methods play VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ important roles in different stages. Usually, inversion-based deblending methods can lead to better separation results compared with the filtering-based deblending methods in suppressing blending noise by using some prior knowledge reasonably [29], [30], especially in the complex subsurface structure. In this paper, we also solely focus on inversionbased approaches.
In simultaneous source acquisition, the sampling grid can be classified into uniform grid and off-the-grid along spatial coordinate in practice. Naturally, deblending methods can be also divided into two categories according to what grid it is designed on. Here, we use 'off-the-grid' in other places in order to distinguish from 'irregular or uniform sampling on regular or uniform grid' [31]. For uniform sampling grid, the traces are placed equidistantly, and the dataset is simple to store and compute in computers. Recently, there are many effective deblending methods for uniformly sampled data, including irregularly sampled muti-source data deblending and reconstruction [32]. All the aforementioned deblending methods belong to this category. Nevertheless, detectors and sources are seldomly distributed at uniform spatial sampling grid as originally designed due to the presence of obstacles, feathering, dead traces and cost limitations [3], [33]. Namely, the field simultaneous source datasets are mostly recorded at off-the-grid positions, especially in complex environments. In this situation, both dataset and position information are needed to be stored in computers. Although, simultaneous source acquisition can enhance the acquisition efficiency by firing more than one source with short time delays, the deviations of detectors and sources between the actual locations and planned locations often lead to off-thegrid blended data, which can produce severe artifacts that aforementioned deblending methods cannot deal with [34]. Nevertheless, subsequent processing steps including velocity analysis, three-dimensional (3D) surface-related multiple elimination, wave equation migration, full-waveform inversion, etc., need a uniform sampling grid data [35]- [37]. In addition, off-the-grid blended data, which unevenly illuminates the subsurface, can induce a distorted acquisition footprint. If not properly deal with it, it will lead to inaccuracies in the subsequent subsurface image, and eventually deteriorates the resolution of the seismic profile.
Therefore, to deal with the issue of off-the-grid data, and meet the requirements of some certain approaches, it is preferable to regularize data beforehand. At present, a common and widely used practice in the field of geophysics is data binning. In this method, off-the-grid traces are simply projected to the nearest bin centres with no amplitude and phase corrections. However, data binning ignores the accurate recording positions of the detectors and sources, which will eventually reduce the subsequent data processing quality [38]. Xu et al. [39] also proposed a binning method with parial normal moveout (NMO), and can provide better regularization results. However, this method still suffers from the simple assumption of NMO. Another approach is the wavefield-based inversion method based on an assumption of a known velocity, which is recognized as a better regularization method because it can combine the exact positions of off-the-grid traces with the wave propagation operator to create new uniform grid traces [40], [41]. However, the calculation speed of this method is relatively slow when dealing with a large off-the-grid blended data set. Obviously, these techniques are not optimally designed to deal with the offthe-grid data.
In recent years, to accurately recover the signal from offthe-grid to the regular grid, many researchers also developed the new regularization procedure or reconstruction method to solve this off-the-grid sampled problem. The band-limited seismic data reconstruction is proposed using nonequispaced discretized Fourier transform (NDFT) by Duijndam et al. [42]. Hindriks and Duijndam [43] extended it to process a 3D seismic data along two spatial directions, showing significant improvement over the conventional binning and stacking methods. In order to handle aliased sampled seismic data, a Fourier reconstruction-based sparse inversion method is proposed using NDFT by Zwartjes and Sacchi [44]. To further improve the reconstruction efficiency, Meng et al. [45] proposed irregularly sampled seismic data reconstruction based on nonequispaced fast Fourier transform (NFFT), and Jin [46] developed a novel damped least-norm inversion-based 5D regularization method by introducing the NFFT. In order to avoid the assumption of linear or quasi-linear events for the Fourier transform-based method, Hennenfent et al. [47] proposed a 2D nonequispaced fast discrete curvelet transform (NFDCT)based reconstruction method with the help of NFFT. Zhang et al. [34] extended it to 3D data reconstruction for off-the-grid missing data along the two spatial directions, and improve the reconstruction performance. In the low rank domain, matrix completion and rank reduction methods are widely used in missing traces reconstruction [48], [49]. Based on a linear event assumption, Lopezet et al. [50] extended matrix completion to seismic data regularization on unstructured grids with discretized Fourier transform (DFT) and NDFT, then it can project the data on the regular grid to the irregular grid. For 3D off-the-grid data regularization, Da Silva and Herrmann [51] proposed a tensor completion, and it shows a good performance in stationary field data.
However, all the aforementioned methods can only recover the traces from off-the-grid to the regular grid. They can not separate the off-the-grid blended data. As a matter of fact, for the separation of off-the-grid simultaneous source data, the conventional approach is to first regularize traces from off-the-grid to the regular grid, and then performs the deblending method. It will reduce the deblending effect and increases the processing time because the blending noise will affect the regularization process. Thus, we need to develop a simultaneous regularization and deblending method for the off-the-grid blended data. However, we can see that few deblending algorithms are proposed to separate off-the-grid simultaneous source data. Nevertheless, simultaneous source data are typically sampled on irregular, non-Carteisian grid 4924 VOLUME 9, 2021 along spatial coordinates due to logistic, topographical, and economic constraints in field [31], [34]. Thus, it is very important to find an efficient deblending method for the off-the-grid simultaneous source data.
The fast discrete curvelet transform (FDCT) has provided remarkable results for deblending of simultaneous-source seismic data on the regular grid [52], [53]. Comparison to other transforms, the localization and sparsity of curvelets can help obtain separated gather accurately for conventional seismic processing workflow. Aiming at the defect of FDCT that is not suitable for the off-the-grid data, Hennenfent et al. [47] proposed the 2D non-equispaced fast discrete curvelet transform (NFDCT) for seismic data reconstruction, and it is effective and robust. However, in their papers, they solved the inverse problem with a spectral projected-gradient method (SPGL1). SPGL1 is difficult and complex to implement on industry scale problems since it is computationally expensive [54]. In addition, the NFDCT is only tested on 2D data reconstruction, not on off-the-grid blended data.
In this paper, we propose a novel method for off-thegrid simultaneous source data using NFDCT via iterative shrinkage thresholding algorithm (ISTA), which can remove blending noise and preserve useful signal effectively while projecting data from off-the-grid to the uniform grid. In order to improve the efficiency, we propose a new exponential thresholding formula decreased from a large value to a small value instead of the conventional thresholding formula, and then applying a soft thresholding operator in the curvelet domain to provide a coherency-based constraint for the iterated model. At each iteration, the effect of constraint is weakened by decreasing the threshold in the curvelet domain, then the most useful signals are well recovered through ISTA. The paper is organized as follows: the simultaneous source acquisition, especially for off-the-grid data acquisition, and its research progress are first introduced. Theory of deblending is then reviewed, and off-the-grid deblending algorithm using the NFDCT is proposed. Next, synthetic and numerically blended data are used to demonstrate the validity of the proposed deblending method. Finally, the discussions and conclusions are summarized.

A. INVERSION-BASED DEBLENDING
In the blending case of two independent sources, we assume primary source remains the same and secondary is blended using a random dithering operator, which means that each shot in primary source has a random time dithering compared with the corresponding shot in secondary source [55], [56].
In the common-receiver domain, the blending problem can be formulated as follows: where d denotes the blended data; d 1 and d 2 denote the seismic records to be separated; denotes the dithering operator of the second source, and its expression can be expressed as where F and F −1 indicate forward and inverse Fourier transforms, respectively, and P is an diagonal operator given by where t k represents the random dithering time for the kth shot in the secondary source. Note that, the signal is coherent and the blending noise is incoherent for each blended record in the common-receiver domain, then the signal can be recovered using the iterative method based on the sparsity constraint in some sparse transform domains [57], [58]. We apply the inverse dithering operator −1 for equation (1), and can obtain another equation: (4) −1 d indicates the pseudodeblended data in which the signal becomes incoherent blending noise and the blending noise becomes coherent signal. By combining equations (1) and (4), we can get an augmented estimation problem: For equation (5), assuming sparse coefficient x = x 1 x 2 is a sparse representation ofd in the curvelet domain C. Then, equation (5) can be reformulated as where the superscript H represents the conjugate transpose.
Then the sparse solution of x can be obtained by solving the following L1-minimum optimization problem [59].
In these expressions, the tilde represents estimated quantities.

B. INVERSION-BASED OFF-THE-GRID DEBLENDING
Generally speaking, the operator C H indicates the conventional FDCT inverse operator, denoted as C H U . When the blended data is uniform grid data, we can get the uniform curvelet coefficients of seismic records to be separated by sloving equation (8) using sparsity-promoting inversion method. At present, most of the existing deblending methods using the FDCT belong to this category. Note that the current implementation of the FDCT still assumes a uniformly sampled grid along all axes. For a non-uniformly sampled grid, FDCT can no longer detect wavefronts because there is no spatial continuity. When blended data is off-the-grid sampled data, it is obvious that deblending via the conventional FDCT inverse operator would not recast the data from off-the-grid to VOLUME 9, 2021 the uniform grid, which will cause severe distortions. To deal with off-the-grid blended data, it is necessary to construct NFDCT inverse operator, denoted as C N . That is A = FC N , then we can not only separate off-the-grid blended seismic, but also project the blended data from off-the-grid to the uniform grid by solving an L1-optimization problem.
C. NONEQUISPACED CURVELET TRANSFORM Candès et al. (2006) proposed 2D FDCT via wrapping specially selected Fourier Transform [52]. The main steps of its forward operator are as follows: (1) obtaining the Fourier coefficients by applying the forward 2D FFT, forming the angular wedges and wrapping each wedge around the origin in the Fourier domain; (2) obtaining the curvelet coefficients by applying the inverse 2D FFT to each wedge. In our paper, the forward operator of FDCT can be formulated as follows: where F is the forward 2D FFT operator defined in step 1.
T is the curvelet tiling operator that links the f -k domain to the curvelet coefficients defined in step 2. Due to the tight-frame property of the FDCT (C U C H U = I), the inverse operator C H U can be simply expressed as As the curvelet construction outlined above is defined in the equispaced Fourier domain. It is obvious that the FDCT cannot be readily used for the off-the-grid data. As the the inverse 2D FFT operator is the Kronecker product, denoted by the symbol ⊗, we can use F H t as the inverse 1D FFT operator along the temporal axis and N x as the inverse 1D NFFT operator along the spatial axes and form a new inverse operator as which projects blended traces from off-the-grid to the uniformly sampled curvelet coefficients. Substituting the operator into equation (8) and replacing the operator C H with operator C N , we can get the value

D. ITERATIVELY SOLVING THE INVERSE PROBLEM
To solve equation (8), however, is very challenging, because equation (8) is very underdetermined. We not only need to separate the blended data, but also need to project data from the off-the-grid to the uniform grid. We plan to use simple and robust ISTA to deal with the undermined ill-posedness problem [60], [61]. The iteration format is as follows: where l n is the step length. nand N denote the nth iteration and the maximum iteration, respectively. T λ denotes the thresholding operator with λ as an input threshold parameter. General speaking, the thresholding operator can be divide into soft and hard thresholding operators which are corresponding to L1-norm and L0-norm sparse constraints [62]. The expression of the soft thresholding operator is given by where x is the coefficient of off-the-grid blended data in the NFDCT domain. The hard thresholding operator can be expressed as Through many numerical tests, we find the soft thresholding operator can achieve a higher signal-to-noise ratio (SNR) than the hard thresholding operator. In this paper, we choose the soft thresholding operator. After achieving the maximum number of iterations, the separated curvelet coefficients can be obtained, and the uniformly sampled deblended data can be obtained as follows:m where C H U is the conventional inverse curvelet transform. The algorithm workflow of the proposed deblending method is shown in Fig. 1. After separating the blended data from the off-the-grid to the uniform grid, it will be beneficial for the subsequent conventional seismic data processing, including velocity analysis, 3D surface-related multiple elimination, wave equation migration, full-waveform inversion, etc.

E. THRESHOLDING STRATEGY
In order to improve the accuracy and efficiency of the proposed deblending method, it is very important to select the threshold parameter because different threshold parameters can produce different deblending results and computation cost. During the iterative process of the ISTA, the threshold value is decreased from a large value to a small value. According to amplitude difference between the incoherent blending noise and coherent useful signals. Most of the curvelet coefficients are removed in the previous iterations, leaving only the strong energy events. When the iteration goes on, the majority of useful curvelet coefficients are retained, only the incoherent blending noise of small energy is removed. The most useful deblended signals are then well recovered. The threshold parameter λ meets the criterion as following: where ε is the minimum threshold. Usually, its value is determined by the users according to the noise standard deviation estimated in the NFDCT domain. Different datasets have different minimum threshold value. Max represents the maximum threshold value. Abma and Kabir [63] proposed a linear threshold formula in the projections onto convex sets (POCS) algorithm, and proves its good performance in reconstructing the missing traces after 400 iterations. This threshold λ n is defined as However, the convergence rate of linear threshold is very slow. Gao et al. [64] proposed the threshold parameter according to the exponential decay, and demonstrated that the exponential thresholding scheme can significantly improve the convergence rate compared with the linear threshold scheme. This threshold λ n is expressed as In order to further improve the speed of convergence, we propose a new exponential thresholding formula [65]. This threshold λ n is expressed as Next, we will prove the superiority of the proposed threshold scheme.

III. SYNTHETIC DATA EXAMPLES A. ANALYSIS AND SETTING OF PARAMETERS
Here, we plan to use two artificially blended synthetic data to demonstrate the validity of the proposed novel deblending method. The used random dithering time series is shown in Fig. 2, which can construct three different blending operators to be used later. To quantitatively measure the deblending performance, the SNR= 20 log 10 f 0 2 / f − f 0 2 (dB) is used as a measure. f 0 denotes the clean unblended data, and fdenotes the deblended data after n iterations. After data deblending, if the deblended data are consistent with the reference ones, and has a high SNR, it can demonstrate the effectiveness of the deblending method. During iterative deblending of off-the-grid simultaneous source data, in order to obtain optimal performance, the scale of NFDCT and FDCT is selected as five and the angle of coarse scale is selected as sixteen after many tests.
In the case of ensuring the deblending accuracy, appropriate threshold parameter can improve the computational efficiency and save calculation time. From three threshold formulas, the new exponential threshold formula has the fastest decay rate, and the linear threshold has the slowest decay rate under the same number of iterations. In order to compare the deblending efficiency and accuracy of three thresholds parameters, we perform the deblending comparison using three different threshold parameters for the offthe-grid blended synthetic data (Fig.7a). Firstly, we select the maximum number of iterations to be 50, and get the relationship curve between every iteration and the SNR after deblending, as shown in Fig. 3. Here, the SNR of the separated primary sources is used for drawing due to the SNRs of the separated primary and secondary sources are almost same. It can be seen that the SNR after deblending is highest at each iteration using the proposed new exponential threshold, followed by conventional exponential threshold, and the last is linear threshold.
Second, we select different maximum number of iterations increasing at a rate of 5, then SNR curve after deblending is generated with maximum number of iterations increasing from 5 to 60 (Fig. 4). It is obviously that the deblending performance of our proposed new exponential threshold is the best among three different threshold parameters. In order to achieve a certain deblended result with the SNR equal to 15 dB, the linear threshold requires 56 iterations, the conventional exponential threshold requires 13 iterations, whereas the proposed threshold only requires about 10 iterations. Therefore, the proposed threshold can reduce  the number of iterations while ensuring the deblending accuracy. In general, the deblending SNR of three threshold parameters increases with iterations number, However, too many iterations require too much computation time, and we find the SNR increases little with the number of iterations when the number of iterations is greater than 30, especially the conventional exponential threshold and the proposed exponential threshold. Therefore, in order to balance time and accuracy, 30 iterations are used in the following synthetic and filed examples.

B. FIRST SYNTHETIC EXAMPLE
To test the iterative performance of the proposed deblending method, we first synthesize a seismic data with four hyperbolic events using the ricker wavelet of 25Hz. Figs. 5a and 5b are two original unblended common receiver gathers uniformly sampled along one spatial coordinate, which contain 256 traces with 2048 samples per trace. The trace space and time sampling rate are 10 m and 1 ms, respectively. In order to get artificially off-the-grid data, we use a forward FFT and an inverse NFFT with a random deviation in the spatial axis of 5 m, as shown in Figs. 5c and 5d. For inspecting the details, Figs. 6a-6d shows zoomed areas which are highlighted by the green dashed boxes in Figs. 5a-5d, respectively. Because the uniformly sampled grid of traces  are seriously discarded, the continuity along wavefronts in off-the-grid data is severely distorted compared with Figs. 5a and 5b, leading to the measured SNRs of only 14.92 dB and 14.45 dB. In this situation, although the nominal sampling grid is 10 m in the spatial direction, the traces are actually located from 0 m to 20 m. In other words, the minimum distance is close 0 m, and the maximum distance is close 20 m between two adjacent traces.
In order to obtain the synthetic blended data non-uniformly sampled in the common receiver domain, a random time dithering codes varying within [−250 250] ms, as shown in Fig. 2a, are applied to construct the blending operator. Then we can obtain the off-the-grid blended seismic data, as shown in Fig. 7a. Fig. 7b is its corresponding F-K spectrum. As shown in Figs. 7a and 7b, we can see that the signal for the primary source is coherent. However, it is contaminated by strong spike-like blending noise caused by the secondary source. The blending noise is relatively easy to be removed for most of deblending methods. Nevertheless, the continuity along wavefronts in off-the-grid blended data is further destroyed, and off-the-grid blended data produces some artifacts, which most of deblending methods can not deal with.  two zoomed areas. Figs. 9e and 9f display corresponding residual sections. It is obvious that the deblended data is very consistent with the unblended data (Figs. 5a and 5b), and the distorted events are effectively corrected, especially for weak events. We cannot see any obvious blending noise, which indicates that all the coherent useful signals are well recovered and the blending noise is effectively suppressed. Furthermore, from the deblending residual sections, we can see that there is almost no loss of the effective energy for the proposed method, which further demonstrates that the deblended result of our proposed method outperforms the one of conventional method. This example shows that if we want to obtain a good perfermance in the separation of off-the-gird blended data, we must use the off-the-grid deblending method, and treat the off-the-gird data correctly in the curvelet domain. Otherwise, it could cause seriously VOLUME 9, 2021 deteriorated deblended results if simply assuming uniformly sampled blended data.
In order to further demonstrate the effectiveness of our proposed method, we randomly decimated data with 50% missing traces in two original unblended data, as shown in Figs. 10a and 10b. Figs. 10c and 10d show two zoomed areas which are highlighted by two green dashed boxes in Figs. 10a and 10b, respectively. Their corresponding F-K spectra are shown in Figs. 10e and 10f. Since there are only 90 traces left, the nominal spatial sampling remains 20 m. The minimum distance is close 0 m, and the maximum distance is close 50 m between two consecutive traces. From their corresponding F-K spectra, we can see that the unblended data is spatially aliased above a temporal frequency of about 35 Hz, which is likely caused by the use of sparse sampling interval in the spatial direction. Moreover, off-the-grid traces produce some artifacts and further enhance the aliasing. Then we obtain the synthetic blended data non-uniformly sampled in the common receiver domain, as shown in Fig. 11a. Fig. 11b is its corresponding F-K spectrum. It can be shown that the blending noise contaminates the useful signals seriously, and off-the-grid sampled data will further affect the deblending process. To remove incoherent blending noise and reduce aliasing in deblending, we test a denser sampling gird. Figs. 12a and 12b are the estimated primary and secondary sources with a spatial uniform sampling internal of 10 m, and their recovered SNRs are 14.88 dB and 14.53 dB, respectively. Figs. 12c and 12d show two zoomed areas in detail. Figs. 12e and 12f show their corresponding residual sections. It is obvious that the separated wavefronts become continuous, and the distorted events are well corrected. We cannot see any evident blending noise, which demonstrates the validity 4930 VOLUME 9, 2021 FIGURE 11. Artificially blended data and its F-K spectrum in the common receiver domain (sampling grid: 20 m). (a) off-the-grid blended data (b) F-K spectrum of Fig. 11a. of the proposed deblending algorithm in separating off-thegrid blended data and improving the resolution with a denser sampling grid. However, due to 50% missing traces, the residuals are still relatively large. Thus, there is a little damage to the useful signal during deblending from off-the-grid to the uniform grid.

C. SECOND SYNTHETIC EXAMPLE
Here, we use a more complex synthetic example generated using an acoustic finite-difference modeling method. Figs. 13a and 13b are the unblended data which represent two different sources uniformly sampled along one spatial coordinate. In each shot gather, there exist 256 receivers and 750 sampling per trace. The trace interval is 12 m and the temporal sampling interval is 4 ms. In order to get artificially off-the-grid data, 256 new non-uniformly sampled traces can be obtained by a forward FFT and an inverse NFFT, with a random deviation in the spatial axis of 6 m, as shown in Figs. 13c and 13d. The nominal spatial sampling still remains 12 m. Figs. 14a-14d show zoomed areas which are highlighted by the green dashed boxes in Figs. 13a-13d, respectively. Because the actual sampling locations of traces are seriously discarded, the continuity along wavefronts in off-the-grid data is severely destroyed compared with Figs. 13a and 13b, yielding a measured SNR of only 12.17 dB and 12.24 dB. In fact, although the nominal sampling interval is 12 m in the spatial direction, the traces are actually located every 12m ± 6 m. In other words, the minimum distance is close 0 m, and the maximum distance is close 24 m between two consecutive traces.
After applying a random time dithering codes varying within [−300 300] ms, as shown in Fig. 2b, we can obtain the synthetic blended data non-uniformly sampled in the common receiver domain, as shown in Fig. 15a. Fig. 15b is its corresponding F-K spectrum. From the Figs. 15a and 15b, we can see that the signal for the primary source is coherent. However, the signal is contaminated by strong spikelike blending interference caused by the secondary source, and the simultaneous sources separation problem is transformed into a spiky noise attenuation problem. Meanwhile, the continuity along wavefronts in off-the-grid blended data is further destroyed, and off-the-grid blended data produces some artifacts, which also affects the deblending process. Relatively speaking, it is very challenging to separate off-thegrid blended data.
Firstly, we also use the conventional FDCT-based deblending method to separate off-the-grid blended data. Figs. 16a and 16b show the deblending performances of the primary and secondary sources, with the corresponding recovered SNRs as 12.40 dB and 12.63 dB, respectively. Figs. 16c and 16d show two zoomed areas of the green dashed boxes in Figs. 16a and 16b in detail. Figs. 16e and 16f show the corresponding residual sections. Clearly, although the conventional FDCT-based method can separate the blended data, it cannot handle off-the-grid blended data  when ignoring the real positions of the sources, resulting in relatively low SNR. Then, the NFDCT-based deblending method is employed to separate off-the-grid blended data.  (Figs. 13a and 13b), and the distorted events are well corrected, especially for weak events. We cannot see any evident blending noise, which indicates that all the coherent useful signals are recovered and blending noise is effectively suppressed. Furthermore, from the deblending residual sections, we can see that there is almost no loss of the effective energy for the proposed method, which further demonstrates that the deblended result using our proposed method outperforms that using the conventional method. This example also shows the importance of treating off-the-gird blended data correctly in the curvelet domain. Otherwise, it could cause seriously deteriorated deblended results if simply assuming uniformly sampled blended data.
In order to further compare the performance between NFDCT-based deblending method and FDCT-based deblending method using the ISTA, the deblended SNRs are compared as shown in Fig. 18, which illustrates that the recovered SNRs increase with iterations for two kinds of deblending methods. However, it can be seen that the recovered SNRs of the NFDCT-based deblending method are always higher than that of the FDCT-based deblending method under the same number of iterations. Moreover, the proposed method can achieve ideal recovered results with only 30 iterations, and little improvements can be obtained beyond that number of iterations. On the other hand, the FDCT-based deblending method needs 20 iterations to get a relatively high SNR at least, and the highest SNR is about 12.10 dB. In order to get this highest SNR, the proposed method requires only 8 iterations, which shows that the proposed method has a fast convergence speed as well as high deblending accuracy.
Meanwhile, we also randomly decimated data with 40% missing traces in the original unblended data, as shown in Figs. 19a and 19b. Figs. 19c and 19d show two zoomed areas which are highlighted by two green dashed boxes in Figs. 19a and 19b, respectively. Here, the nominal spatial sampling remains approximately 20 m, the minimum distance is close 0 m, and the maximum distance is close 48 m between two consecutive traces. From their corresponding F-K spectra shown in Figs. 19e and 19f, we can see that off-the-grid traces produces some artifacts and enhance the aliasing. Then we obtain the synthetic blended data non-uniformly sampled in the common receiver domain using the random blending operator, as shown in Fig. 20a. Fig. 20b is its corresponding F-K spectrum. It can be shown that blending noise contaminates the useful signals seriously, and off-the-grid sampled data will further affect the deblending process. To remove incoherent blending noise and reduce artifacts in deblending, we plan to reduce sampling gird and get more traces during the iterative deblending from the off-ther-grid to the uniform grid. Figs. 21a and 21b are the estimated primary and secondary sources with a spatial sampling internal of approximately 8 m. In other words, there are 384 traces in deblended data. Figs. 21c and 21d show two zoomed areas in detail. Figs. 21e and 21f show corresponding residual sections. It is obvious that the separated wavefronts become continuous, and the distorted events are well corrected, especially for weak events. We hardly see obvious blending noise. This example also demonstrates the proposed deblending method can not only separate the off-the-grid blend data, but also improve the resolution with a denser sampling grid. To further prove its validity, one field data example is then tested and analyzed.

IV. REAL DATA EXAMPLE
Two different sources of the field data contain both 180 traces with 25m as the spatial sampling interval, and 750 sampling per trace with 4 ms as the time sampling interval. Currently, we are unable to get the off-the grid field data; therefore, we randomly decimated uniformly sampled field data with 30% missing traces as original unblended field data (Figs. 22a and 22b). In this situation, the nominal spatial sampling is approximately 36 m. Figs. 22c and 22d show two zoomed areas which are highlighted by two green dashed boxes in Figs. 22a and 22b, respectively. Their corresponding F-K spectra are shown in Figs. 22e and 22f. We can see that the field reference data are spatially aliased above a temporal frequency of about 30 Hz, which is likely caused by the use of sparse sampling interval in the spatial direction. Moreover, off-the-grid traces produce some artifacts and further enhance the aliasing. To artificially simulate the blended acquisition in the common receiver domain, a random time dithering code varying within [−400 400] ms, as shown in Fig. 2c, is used to build the blending operator. After applying the blending operator to the unblended data of the second source, the blended seismic data is obtained in Fig. 23a, and its corresponding F-K spectrum is shown in Fig. 23b. It can be shown that blending noise contaminates the useful signals seriously, and off-the-grid sampled data will further affect the deblending process. To remove incoherent blending noise and reduce aliasing in deblending, we plan to test a denser sampling internal. Figs. 24a and 24b are the estimated primary and secondary sources with a spatial sampling internal of approximately 15 m. In other words, there are 300 traces in deblended data. To see the performance clearly, Figs. 24c and 24d show two zoomed areas in detail. It is obvious that the separated wavefronts become continuous, and the distorted events are well corrected except for some weak events. We cannot see any evident blending noise. Figs. 24e and 24f plot their corresponding the amplitude spectra. We can see that the aliasing above a temporal frequency of about 30 Hz almost disappear because of a denser spatial sampling interval used in the deblending stage, and the resolution of deblended seismic data is greatly improved. This example further demonstrates that our method can obtain improved deblending perfermances on continuities of the events while projecting separated data from off-the-grid to a denser regular grid.

V. DISCUSSION
The selection of threshold scheme is very important in separation of simultaneous source data. In this paper, the proposed new threshold formula is better compared with the other two threshold formulas in terms of deblending accuracy and computing efficiency. However, there must exist other better threshold formulas in the subsequent research on deblending, and different blended data have different optimal threshold formula, thus the new threshold formula should be developed. Meanwhile, the minimum threshold has a certain influence on the deblended result, and it can be determined by the 4934 VOLUME 9, 2021    manually or empirically to decide which minimum threshold is better according to the recovered SNR. On the other hand, through numerical simulation, we can find that the effect of soft thresholding operator is better than that of hard thresholding operator in separation of the off-the-grid simultaneous source data. However, the main purpose of this paper is not to propose the thresholding operator itself but focus on a novel deblending method of off-the-grid simultaneous source data, so we do not test it in detail.
Although the aforementioned NDFT-based or NFFT-based interpolation methods can recover the signal from off-thegrid to the regular grid with a high precision, they assume that the data consist of linear or quasi-linear events. These methods can not obtain an ideal regularization result when the data consist of complex curved seismic events. Furthermore, we do not find NDFT-based or NFFT-based deblending method to be reported in the literatures, or other similar published deblending methods for the off-the-grid blended data. In this paper, it is the first time to introduce the NFDCT to the deblending of simultaneous source data, and propose the novel iterative deblending of off-the-grid simultaneous source data. This proposed method in this paper can not only deal with curved event directly, but also separate the blended data from off-the-grid to a designated regular grid with a satisfactory deblending accuracy and computing efficiency. However, the proposed off-the-grid deblending method is not compared with other ones as no other off-the-grid deblending method is reported. In other words, our purpose is to merely present new avenues for the separation of off-the-grid simultaneous source data.
Also note that, although we performed the off-the-grid deblending in the common-receiver domain, the proposed method can also be applied in the domains other than common shot domain. However, our proposed method is only suitable for 2D seismic data, not for 3D seismic data. True 3D prestack seismic data are most naturally described in 5D domain (four spatial dimensions and one time dimension). In order to take all physical dimensions of seismic data into consideration for the deblending, we should use a 3D offthe-grid deblending method. However, it will lead to another issue of computational efficiency. Although we proposed a new exponential threshold formula, computational efficiency remains an important concern when dealing with a large seismic data set. Therefore, for 3D separation of off-thegrid blended data, it is required to present an approach that can significantly improve the computational efficiency in ensuring the separation accuracy. This is also a research topic that we need to solve in the future.

VI. CONCLUSION
In this paper, an effective iterative deblending algorithm projecting the blended data from off-the-grid to the uniform grid is proposed via NFDCT in the common receiver domain. By constructing NFDCT inverse operator, we can obtain accurate curvelet coefficients that link blended data from off-the-grid to the uniform grid. In the curvelet domain, the blending noise turns into low-amplitude coefficients when using the random dithering operator to construct the blending operator, which can be suppressed effectively using the ISTA via the soft thresholding operator. Once the uniform effective curvelet coefficients are obtained, separated gathers on uniform grids can be estimated through the conventional inverse curvelet transform. The successful performance of the proposed method is fully proved using both synthetic and artificially blended field data examples. Compared with the conventional FDCT-based deblending method, this proposed method can not only separate offthe-grid blended and aliased data, but also project separated data from off-the-grid to a desired uniform grid with a better deblending performance and fine resampling grid, which shows that it is a promising tool for off-the-grid data processing. 4936 VOLUME 9, 2021