Seismic data have two dimensions of time and space. The seismic traces are usually uniformly sampled in the temporal domain, and the Nyquist–Shannon sampling theorem is satisfied. However, due to existing obstacles and economic restrictions, there is always missing traces in seismic data, resulting in non-uniformly sampled, even aliased in the spatial domain. Moreover, seismic data are inevitably polluted by random noise. Complete and clean seismic data need to be reconstructed from recorded noisy and incomplete seismic data.

The relationship between the recorded seismic data and the complete data can be formalized as: y = Lx , (1) where y and x are the recorded (incomplete) seismic data and complete data, respectively. Both y and x are constructed by arranging all the traces of the recorded (Y) and complete (X) data matrix into vectors. L represents the measurement matrix, which is a block diagonal matrix and is defined as L = [ I I 0 I ⋱ I ] (each 0 represents a missing trace, and each I represents a recorded seismic trace).

It should note that both y and x are infected by random noise. The purpose of the proposed method is to reconstruct the expected clean and complete data x ∗ from y. Solving Eq. 1 can only give rise to complete noisy data. However, worse still, seismic wavefield reconstruction aims to recover a high-dimensional signal cube from a smaller number of measurements. formula (1) is underdetermined and has infinitely many solutions. The proposed countermeasure to address these two bottlenecks will be described in detail in the following subsections.

It is possible to reconstruct sparse signals from highly incomplete sets, which can be done by sparse representation. The seismic data are compressible and sparse in the curvelet and Fourier domains. The sparsity adds effective constraints to the underdetermined inverse problem, thereby reducing the distress of the multiplicity. So, an effective strategy for the underdetermined problem (1) is that we can transform x into the Fourier or curvelet domain (as shown in Eq. 2) and then formulate the reconstruction problem as an L₁-norm regularization (as shown in Eq. 3) (Mansour et al., 2013). y = Am , A = L F H , (2) where F is a sparsifying operator; that is, Fourier or curvelet transform, m is the Fourier or curvelet coefficient, and superscript H represents the conjugate transpose. arg min m ‖ m ‖ 1   s . t .   ‖ y − Am ‖ 2 2 ≤ ε , (3) where ‖ ⋅ ‖ 1 and ‖ ⋅ ‖ 2 denote L₁- and L₂-norms, respectively; ε is a fitting error. We have learned that Eq. 3 is a basis pursuit denoise (BPDN) problem; it can be optimized by the spectral projection gradient method. However, L₁-norm is non-differentiable and the augmented model for Eq. 2 is adopted (Yin 2010): m ∗ = arg min m ‖ m ‖ 1 + 1 2 α ‖ m ‖ 2   s . t .   ‖ y − Am ‖ 2 2 ≤ ε , (4) where typically α = 1 .

Solving Eq. 4 can only get complete noisy data, that is, m ∗ is corrupted by noise. Consequently, the next step is to denoise the reconstruction m ∗ . We first rearrange m * into a 2D matrix M.

Due to the repetitive texture structure, the prestack seismic gathers are low-rank matrixes (Ma 2013). In a transformed domain, that is, Fourier and curvelet domains, the coefficient matrix M is also characterized by the low-rank feature (Nazari Siahsar et al., 2016). Suppose M is sullied by Gaussian noise, we have the noisy coefficient matrix M about a clean low-rank matrix M den of interest. A natural approach to this problem is singular value thresholding (SVT), which formalizes noise suppression as: M den = SV T λ ( M ) = arg min M den λ rank ( M den ) + 1 2 ‖ M − M den ‖ F 2 , (5) where ‖ ⋅ ‖ F is the Frobenius norm and M = UΣ V H ( U = [ u 1 , u 2 , ⋯ u n ] and V = [ v 1 , v 2 , ⋯ v n ] are left and right singular vectors, respectively; Σ = d i a g ( σ 1 , σ 2 , ⋯ , σ n ) is a diagonal matrix with singular values of M on its diagonal) is the singular value decomposition (SVD) of M. The solution of (5)is given by retaining only the part of the expansion with singular values exceeding λ (Candes et al., 2013): M den = ∑ ( σ i − λ ) + u i v i H , (6) where x + = max ( x , 0 ) and λ is the threshold value which can be estimated by Stein’s unbiased risk estimate (SURE).

The clean and complete seismic section is reconstructed from the denoised coefficient matrix M den just via matrix-vector multiplication, which can be formulated as y den = F m den , (7) where m den is the vectorization of M den .

Now we recall the simultaneous reconstruction and denoising methods as follows:

1. Noise level estimation: estimating the parameter ε for (4) (discussed in Subsection 2.5).

The original Bregman method was proposed to restore noisy and blurry images via total variation (TV) regularization (Osher et al., 2005), and then be formulized to solving the basis pursuit (BP) problem (Yin et al., 2008), where the author had proved that Bregman iterative method was equivalent to the augmented Lagrangian method.

The Bregman distance with respect to a convex function J between two points u and v is defined as D J p ( u , v ) : = J ( u ) − J ( v ) − ≥ 〈 p , u − v 〉 , (8) where p ∈ ∂ J ( v ) ( ∂ J ( v ) denotes the subdifferential of J at v), and symbol 〈 ⋅ , ⋅ 〉 denotes the inner product. The original Bregman method is an iterative regularization procedure which focuses on the problem arg min x J ( x ) s . t .   Ax = b , that solves a sequence of convex problems, which can be formulated as x k + 1 = min x D J p k ( x , x k ) + 1 2 ‖ Ax − b ‖ 2 2 . (9)

As a variant of the original Bregman method, the linearized Bregman algorithm imposes two modifications to the original method: (1) substitute the last quadratic penalty term ‖ Ax − b ‖ 2 2 with its linearization 〈 A T ( A x k − b ) , x 〉 (superscript T denotes the transposition), (2) add a new quadratic penalty term 1 2 α ‖ x − x k ‖ 2 2 . Then, we get the following updates for the linearized Bregman algorithm (Yin, 2010): { x k + 1 = min x D J p k ( x , x k ) + 〈 A T ( A x k − b ) , x 〉 + 1 2 α ‖ x − x k ‖ 2 2 p k + 1 = p k − A T ( A x k − b ) − 1 α ( x k + 1 − x k ) , (10)

Based on the equivalence between the linearized Bregman algorithm and the gradient descent, Yin (2010) generalized the linearized Bregman by integrating the gradient-based optimization techniques to make the method more accurate and much faster.

To address Eq. 4, we take J ( x ) = ‖ x ‖ 1 , and the problem can be solved with the linearized Bregman algorithm by the following iterative scheme: { v k + 1 = v k + δ k ( b − α A S h r i n k ( A T x k ) ) x k + 1 = α S h r i n k ( A T v k + 1 ) , (11) where δ k is the step size at iteration k.

In this subsection, we focus on the vital issue of threshold parameter selection because the unreasonable threshold inevitably limits the denoising ability even to generate artifacts. Specifically, excessive shrinkage results in a large bias, while excessively low threshold value results in a high variance. Common methods of parameter selection include discrepancy principle, cross-validation (Li et al., 2019), L-curve (Hansen 1992), and estimation of the mean-squared error (MSE) (Batu and Cetin 2011), and so on. Ramani et al. (2012) summarized the advantages and disadvantages of these methods and recommended using MSE to optimize the regularization parameter quantitatively. A classical unbiased estimator of the MSE is SURE (Stein 1981), which was reported to be used in quantitative parameter optimization (Ramani et al., 2012; Candes et al., 2013; Rasti et al., 2015).

MSE-estimation-based methods tune the regularization parameter λ to minimize the MSE or risk: R λ = E ‖ M den − M 0 ‖ F 2 , (12) where M 0 is the ground truth.

When the seismic data are corrupted by additive random noise, it is possible to obtain an unbiased estimate of the risk via SURE (Candes et al., 2013): R λ = − a b τ + ∑ i = 1 min ( a , b ) min ( λ 2 , σ i 2 ) + 2 τ 2 div ( SVT λ ( M ) ) , (13) where τ is the standard variance of the additive random noise; a and b are the numbers of rows and columns, respectively; symbol “div” denotes the divergence operator. From Candes et al. (2013), it is clear that div ( SVT λ ( M ∗ ) ) = ( 2 | a − b | + 1 ) ∑ i = 1 min ( a , b ) ( 1 − λ σ i ) + ∑ i = 1 min ( a , b ) I ( σ i > λ ) + ∑ i ≠ j , i , j = 1 min ( a , b ) σ i ( σ i − λ ) + σ i 2 − σ j 2 , (14) where I is the indicator function.

Noise level is an extremely important parameter to achieve a guaranteed performance of both sparse representation and SVT. Specifically, in Eqs 3, 4, parameter ε prescribes the desired fit in Am ≈ y , which can be substituted by noise variance τ 2 in the random noise scenario. Meanwhile, a reliable estimate of the noise level is needed for SURE to get the proper threshold value. In this work, we choose the multiple regression theory-based approach (Roger and Arnold 1996; Chein-I and Qian 2004; Bioucas-Dias and Nascimento 2008) for noise level estimation. The seismic gather, especially common-offset gather, exhibits high spatial coherence between neighboring traces in both the time and frequency domains, which guarantees the successful application of the multiple regression theory-based approach.

Let Y = [ y 1 , y 2 , ⋯ , y M ] denote a M × N matrix holding the seismic gather, with each column y i accommodating a single trace. Let Y ∂ i = [ y 1 , y 2 , ⋯ , y i − 1 , y i + 1 , ⋯ , y M ] be the explanatory data matrix of y i . Let us assume that y i can be linearly represented by other seismic traces in the gather. Therefore, we can express y i as y i = Y ∂ i β i + ξ i , (15) where β i denotes the regression vector of length M − 1 , ξ i denotes the modeling error of length N. The noise for each trace can be derived from the following formula: ξ ^ i = y i − Y ∂ i β i . (16)

Then, it is easy to calculate the noise variance using the standard formula.

In this section, we validate the effectiveness of the proposed method by simulated data. We first consider a classic problem of sparse solution recovery to emphasize the importance of ε in (3) and 4. We then provide a layered model and perform wavefield modeling, followed by randomly decimated data with 25% missing traces. A set of experiments, that is, noise level estimation and simultaneous reconstruction and denoising, are conducted. All experiments deal with common-shot gathers and common-offset gathers separately, and the purpose is to compare different method combinations.

In (3) and 4, parameter ε prescribes the fitting error. However, the importance of parameter ε is often ignored in literature or software packages. ε = 1 e − 3 is usually recommended. The parameter ε is essentially an estimate of the noise level of the seismic signal and can be quantified by an upper bound on its norm. In this subsection, we make a digression to emphasize the importance of ε, with a compressive sensing problem based on random Gaussian matrices. The size of the matrices is 2000 × 4000 . We compare SNRs of reconstructed solutions on two sets of values for the parameter ε: 1 e − 3 and L₂-norm of deterministic noise. We first fix the sparse signal with 5% non-zero entries. The results of the model (4) with noisy observations of different SNRs are plotted in Figure 1A. On the contrary, SNRs of the observations are fixed and equal to 10 dB, we change the sparsity of the sparse solutions, and plot the calculation results of Eq. 4 in Figure 1B. Both Figure 1A,B indicate that the proper parameter value of ε is a sufficient condition for recovering the sparse solution to high accuracy. Taking L₂-norm of deterministic noise can be potentially much better than the default. We fix sparsity  =   2  % and SNR of the observation is 10 dB and plot the computed coefficients (red) on top of the ground truth (blue) in Figure 2. Model (4) finds better solutions when ε takes the L₂-norm of noise.

We first use a five-layer model (Figure 3A) to perform wavefield modeling. The synthetic line has 48 shot gathers, and 48 traces per each shot gather and CMP gather. Source-station and receiver-station spacings are both 50 m. The data are recorded for 1.000 s, with a sample interval of 0.002 s. Figure 3B,C show the reflectivity and seismogram gathers of the first CMP gather. Figure 3D presents the 24th common-offset gather (offset is 50 m). CMP gather is characterized by hyperbolic events, while the common-offset gather consists of flat events, which share approximately equal two-way times. Compared with CMP gather, each trace in common-offset gather shows better similarity and simpler texture structure. This makes it easier to deal with common-offset gather in both the time-space domain and Fourier domain.

The field-recorded seismic data are inevitably perturbed by random noise. Therefore, we add random noise at various levels to the synthetic data to make the SNR of each noisy trace vary within a range of 1–10 dB. The first noisy CMP and the 24th noisy common-offset gathers are shown in Figure 4A,B, respectively. To simulate the irregularly sampled seismic data, we randomly zero out 25% of the receivers and artificially create some gaps to increase the difficulty of the experiments. The incomplete noisy CMP and common-offset gathers are shown in Figure 4C,D, respectively.

To inspect the amplitudes of the seismic data at different frequencies and wave-numbers, we perform 2-D Fourier transform on both complete noisy CMP and common-offset gathers to transform them into frequency-wavenumber (f-k) domain, as indicated in Figure 5A,B, respectively. Note that obvious spatial aliasing exists in the f-k spectra for the CMP gather that indicates insufficient sampling of the data along the space axis. However, COG is free from aliasing.

Most relevant literature (Zhang et al., 2019; Cao et al., 2020) characterizes the sparse representation of seismic data based on curvelet transform. We use the approved algorithm (curvelet domain approach) to deal with the common-offset and the CMP gathers, separately. The experiment served two purposes: (1) Spatial aliasing has severe effects on the performance of curvelet domain approach in the time-offset domain (common-shot and CMP gathers); (2) Noise estimation and wavefield reconstruction can be done better in the time-midpoint domain (common-offset gather).

Determining an optimal threshold value is an essential first step for sparse representation and the implicit reconstruction and denoising. We first apply the multiple regression theory-based approach to evaluate the noise level both in time-midpoint and time-offset domain, respectively. In implementing this step, we get poor result with CMP gather, as demonstrated in Figure 6A. The estimated L₂-norm of the noise aliasing in each trace is much lower than the actual noise level. This is basically because the multiple regression theory-based method assumes that each trace can be linearly represented by other seismic traces in the gather. CMP gather is characterized by hyperbolic events that have disparate TWT. They share a low correlation between neighboring traces that dissatisfy the multiple regression theory-based method. The procedure is also applied to common-offset gather, and the obtained result is plotted in Figure 6B that indicates that the estimated value is approximately correct.

To fill in the gaps and denoise the incomplete noisy data using sparse representation, we use the curvelet transform as “sparse basis” to create matrices A in Eq. 2. The linear relation between the y and m inspires us to consider estimating complete sparse and low-rank coefficient vector m ∗ by solving Eq. 4 via the linearized Bregman method. We take ε as the Frobenius norm of noise to balance the fidelity and sparsity. Figure 7A,B are based upon the processing of the CMP gather, and they illustrate reconstructed data (first CMP gather) using the curvelet dictionary and corresponding residual section. It does not seem plausible that the curvelet dictionary is to blame for false artifacts, because the curvelet transform has only a weak ability to anti-aliasing, even though it has proven to have an outstanding ability that enables a nearly optimal representation of seismic data. However, such spatial aliasing can be maximumly remitted in the time-midpoint domain because the common-offset gather is characterized by flat, even approximate horizontal events. The reconstructed results of the first CMP and the 24th common-offset gathers using the time-midpoint domain approach are shown in Figure 7C,E, respectively. Figure 7D,F show the corresponding residuals of the CMP and common-offset gathers, respectively. As expected, the curvelet reconstruction performs much better in the time-midpoint domain than in the time-offset domain. All missing seismic traces are perfectly reconstructed. The defect is that random noise is not well suppressed.

We next examine whether Fourier transform operator does it as good as the curvelet transform operator. The results from the Fourier transform operator associated with both the time-offset and the time-midpoint domains approaches are plotted in Figure 8A,C, respectively. The corresponding residual images are shown in Figure 8B,D, respectively. There are no apparent differences between curvelet and Fourier transform operators. However, to our surprise, the Fourier transform operator does slightly better than the curvelet transform operator both in the time-offset domain and the time-midpoint domain. The possible reason is that the simulation data used in this experiment are spatially under-sampled. The curvelet transform operator tended to express the local tectonics of the seismic section, and its anti-aliasing ability is weak. The opposite effect is observed for the Fourier transform operator, which points to the ability of the Fourier transform operator to address issues that arose in the field of aliased data. There is no doubt that the computational complexity of the Fourier transform is much lower than the curvelet transform. Therefore, the Fourier transform operator is the best option.

We have experimentally verified the superiority of the time-midpoint domain approach and Fourier transform operator in the previous subsection. Naturally, a combination of the two is also appreciated. The assay is performed as described in Subsection 2.1. Figure 9A shows the recovered common-offset gather in the time-midpoint domain, and Figure 9B is the corresponding residual image. Both two objectives of reconstruction and denoising have been achieved. Next, we sort the data into common-shot gathers to obtain the final result shown in Figure 9C, which is very close to the clean synthetic seismogram (Figure 3C). Figure 9D is the corresponding residual image. The proposed approach improves the SNR of seismic data as compared with the previous two experiments.

To quantitatively compare the reconstruction performance in detail, we plot the noise-free single trace which is a missing trace deleted in the preceding step and the reconstructions with different methods in Figure 10. Figure 10A,E plot the clean and noisy signal (SNR = 4.421 dB), respectively. Results of three independent experiments are all obtained by time-midpoint domain approach, as plotted in Figure 10B–D. Full recovery of the missing trace is perfectly achieved in each of the three experiments. Fourier and curvelet transform operators suppress the random noise and improve the SNR of data to some extent, as plotted in Figure 10B,C. The SNRs of reconstructed traces of the two methods are 6.994 and 7.812 dB, respectively. Comparing the two methods, we can see that the superiority of the proposed method is clear (Figure 10D). SNR of the proposed method is equal to 18.243 dB.

In this section, we apply the proposed simultaneous reconstruction and denoising method to a prestack data set. The original raw data consist of 1,001 shot records with 120 channels per shot record. The minimum and maximum offsets are 262 and 3237 m, respectively. Both the source station interval and receiver station interval are 25 m. This study extracts 48 shot gathers from the raw data to ensure that the source station interval is 50 m, exactly twice the original source station interval. Then, each shot gather is down-sampled to reduce the number of channels to 60; that is, the receiver station interval is increased to 50 m. Figure 11A,B show the observed shot gather and common-offset gather, respectively. Intuitively, the shot gather consisted of hyperbolic reflected events and gradient direct and seabed reflection waves. However, the common-offset gather is characterized by much flatter events, which is very well adapted to reconstruction and denoising. Such down-sampling may give rise to spatial aliasing. Figure 12A shows the f-k spectra for the shot record with strong spatial aliasing, as indicated by the solid arrows. As mentioned previously, processing seismic data in the time-shotpoint domain is a simple and efficient strategy. Figure 12B demonstrates the f-k spectra for the common-offset gather shown in Figure 11B, and spatial aliasing is not visible as of the common-shot gather. To simulate the observed incomplete field data, we just take out the columns with missing receivers. 25% of the receivers are deleted randomly in this study. Figure 11C,D show the incomplete shot and common-offset gathers, respectively.

We perform three shortened versions of the simultaneous reconstruction and denoising experiments on the truncated seismic data, that is, curvelet transform operator, Fourier transform operator, and the proposed sparse and low-rank regularization approach. The complete workflow of each experiment consisted of three flows: (1) transforms the seismic data to the time-shotpoint (common-offset) domain, (2) noise level estimation, and (3) reconstructs the clean and complete seismic section by different methods. The reconstructed common-offset gathers of three different methods are contrasted in Figure 13. The reconstructed result using the curvelet transform operator and the associated difference section are shown in Figure 13A,B, respectively. The narrow gaps are well recovered, while the reconstructed results of the big gaps are not as good as it would be expected (indicated by the arrows in Figure 13A,B). Recovered common-offset gather using Fourier transform operator and the associated difference section are shown in Figure 13C,D, respectively. The missing traces were well recovered, in particular at the big gaps. However, reconstructed data seem to have significant residual noise in the big gap at 0–0.5s. Figure 13E,F demonstrate the interpolated result and the corresponding residuals of the proposed strategy. As can be seen, the recovery is comparable to the Fourier transform operator. A further advantage is that we can obtain much cleaner data that does not contain the artifacts before the first arrivals.

The reconstructed and denoised common-offset gathers are transformed into the common-shot gathers. Figure 14 compares the recovered shot gathers of all the three methods. Figure 14A,C,E indicate the recovered shot gathers of the three different approaches, respectively. Figure 14B,D,F show the corresponding residuals. Similar conclusions can also be obtained based on the comparison of reconstructed results and corresponding residuals of the three methods. The proposed method is more effective in seismic wavefield reconstruction and denoising than the curvelet and Fourier transforms operators. To better compare the details of reconstructed results by the above three methods, we compare the recoveries of a single trace, as plotted in Figure 15. Figure 15A plots the real noisy trace (Trace No. 39 from the noisy shot shown in Figure 11A). Figure 15B–D plot the recovered traces of the three different approaches, respectively. Figure 15E–G correspond to the removed noise associate to each method. These results can be explained as follows. The result of the curvelet transform operator is not very good. The result of the Fourier transform operator is close to the exact solution. However, we should note that the presence of noise before the first arrival (at 0–1.5 s, as indicated by the arrows in Figure 15C,F) is caused by not very good localization of Fourier transform. The proposed method compensates for the defect exactly, and the recovered signal does not contain the evident artifacts and noise. Therefore, the proposed sparse and low-rank regularization approach should be the best choice owing to its effectiveness and fidelity.