Reconstruction and denoising of high-dimensional seismic data via Frobenius-nuclear mixed norm constraints

Abstract

1. Introduction

Field seismic data collection faces various constraints, including complex near surface conditions, rugged terrain, and the existence of industrial and agricultural facilities. These factors can lead to irregular spatial distributions or missing data in the collected seismic data. Moreover, the collection process is susceptible to environmental factors, resulting in diverse interference noises being introduced into the seismic data. These noises negatively affect the signal-to-noise ratio (SNR), sometimes obscuring seismic signals. Consequently, this can detrimentally affect subsequent processing outcomes such as seismic wavelet estimation, surface-related multiple suppression, first arrival picking, and parameter inversion. Hence, noise suppression and data reconstruction are imperatively important preprocessing steps to enhance seismic data quality.

In recent years, numerous studies have investigated various methods on improving the SNR and increasing the precision of reconstructing seismic data. These approaches are primarily categorized into filtering, sparse transformation, low-rank matrix/tensor decomposition, and deep learning techniques. Filtering methods predominantly exploit the predictability of signals across different domains such as |$t - x$|⁠, |$f - x$| and |$f - k$| to mitigate random noise or interpolate seismic data (Spitz 1991, Gülünay 2000, 2017, Naghizadeh and Sacchi 2009, Curry 2010, Zhang et al. 2019). In addition to the aforementioned methods, statistical-based filtering approaches are utilized for denoising by assuming that data within a local window follow a specific probability distribution. Examples of such methodologies comprise Gaussian filtering (Fu et al. 2018), median filtering (Liu et al. 2022), and mean filtering (Bonar and Sacchi 2012). Sparse transformation methods exploit the sparsity of data following transformation through the use of basis functions to reduce noise or assist in reconstruction. Basis functions in the transformation domain are typically divided into two categories: fixed basis functions such as Fourier (Xu et al. 2005, Dobróka et al. 2017), radon (Trad et al. 2002, Zhang and Lu 2014, Xue et al. 2016, Wang et al. 2017, Li and Sacchi 2021, Geng et al. 2022), curvelet and seislet bases (Starck et al. 2002, Fomel and Liu 2010, Liu et al. 2016, Liu et al. 2018, Naghizadeh and Sacchi 2018, Wang et al. 2018, Zhao et al. 2023), and data-driven basis functions, with dictionary bases being the most commonly used (Beckouche and Ma 2014, Zhu et al. 2017, Zhou et al. 2024).

In addition to recognizing the sparsity of seismic data, it is assumed that complete data without noise exhibits low-rank characteristics, while random missing data and noise significantly elevate the rank of the data. Therefore, low-rank tensor completion methods are utilized for the reconstruction and noise suppression in 2D or higher-dimensional seismic data. Ulrych et al. (1999) introduced a feature image approach that uses singular value decomposition (SVD) after NMO. Multichannel singular spectrum analysis is a data-driven technique rooted in the examination and investigation of multichannel time series. The multichannel singular spectrum analysis method (Sacchi 2009) typically entails reorganizing the data in the frequency-space domain into a Hankel matrix, followed by reducing the rank of the Hankel matrix using the SVD technique. Subsequently, the reconstructed Hankel matrix is summed along the anti-diagonal direction and averaged. Kreimer and Sacchi (2012) employed higher-order SVD (HOSVD) to reconstruct 5D pre-stack seismic data after NMO. Gao et al. (2017) proposed a rapid and efficient low-rank tensor completion algorithm—parallel submatrix factorization (PSMF)—which unfolds the fourth-order tensor in the time-space domain along the mode to approximate a square matrix, thereby enabling more accurate data tensor recovery. When seismic data is contaminated by abnormal noise, Carozzi and Sacchi (2019) developed a robust parallel matrix factorization tensor completion method, which uses the |${{L}_1}/{{L}_2}$| norm to measure the data error term and applies it to the reconstruction of 5D seismic data. Overall, by selecting different low-rank approximation methods (Battaglino et al. 2018, Cavalcante and Porsani 2021, Liu et al. 2022, Wang et al. 2023, Zhang et al. 2023) or input data types, approaches based on low-rank tensor completion have been successfully applied in seismic data reconstruction.

With the advancement in computing capabilities and the evolution of big data analytics, deep learning has experienced widespread application in geophysics. As a novel data-driven technology, deep learning obviates the need for data to adhere to specific prior assumptions (such as sparsity, low rank, etc.) and can model more intricate data, thus making headway in seismic data reconstruction and noise suppression. Yu et al. (2019) introduced convolutional neural networks for seismic data noise attenuation. Wang et al. (2019) devised an eight-layer residual network (ResNet) for interpolating regularly missing data, demonstrating high accuracy in noise-free conditions. In the U-Net network (Ronneberger et al. 2015), there is a skip connection cascaded with the corresponding upsampling for each downsampling, which can fuse features of different scales and improve image segmentation performance. Therefore, the U-Net network can be introduced into the noise suppression and reconstruction of seismic data (Fang et al. 2021, Nakayama and Blacquière 2021, Park et al. 2021) to improve the noise suppression and reconstruction results. Tang et al. (2023) introduced a novel approach named rank-reduction U-Net (RRU-net), which extends the U-Net framework. Through the replacement of t-x domain seismic data with a Hankel matrix, the RRU-net exhibits superior performance compared to conventional U-Net methods in reconstructing and denoising DAS-VSP seismic data. The aforementioned deep learning methods belong to supervised learning, which requires a large number of labeled samples for network training with disadvantage of time-consuming. Liu et al. (2023a) proposed an unsupervised deep learning method that eliminates the requirement for training labels meanwhile achieving ground roll and scattering noise attenuation. In short, more and more deep learning-based methods are being used in seismic data reconstruction and denoising (Liu et al. 2022, 2023b, Liao et al. 2023, Yang et al. 2023, Dong et al. 2024), but traditional methods have certain advantages in label generation, accelerating network learning, and production applications, and can complement deep learning methods.

In situations where the underground geologic structure exhibits complexity or lateral velocity variation, the conventional reflection events resulted from dynamic correction of the collected common midpoint (CMP) data may fail to meet linear assumptions while still retaining low-rank characteristics. Therefore, this paper proposes a novel approach for high-dimensional seismic data reconstruction and denoising based on Frobenius/nuclear mixed norm constraints, with low-rank tensor completion as its core principle. By incorporating non-convex Frobenius/nuclear mixed norms onto the mode |$(m,n)$| unfolding matrix of the 4D data tensor within the frequency-space domain, the proposed model can more accurately approximate the rank function of the matrix. This technique is employed in the noise suppression and regularization reconstruction process of high-dimensional seismic data, aimed at enhancing the quality of subsequent seismic wave migration images and inversion imaging. Consequently, it could contribute toward improving reservoir characterization for drilling success during the costly drilling operation.

2. Methodology

2.1. Low-rank tensor completion

Low-rank modeling of noisy or incomplete data is often expressed as a low-rank matrix approximation or completion problem, with widespread applications in fields such as image restoration, system recommendation, and seismic data reconstruction. In the context, |$x \in$|£, |${{\bf x}} \in$|£ |${^I}$|⁠, |${{\bf X}} \in$|£|${^{{{I}_1} \times {{I}_2}}}$|⁠, |$\mathcal{X} \in$|£|${^{{{I}_1} \times {{I}_2} \times {{I}_3} \times \cdots {{I}_N}}}$| refers to different types of tensors: x represents a scalar (zero-order tensor), |${{\bf x}}$| represents a vector (first-order tensor) with size I, |${{\bf X}}$| represents a matrix (second-order tensor) with size |${{I}_1} \times {{I}_2}$|⁠, and |$\mathcal{X}$| represents an N-order tensor (N > 2) with size |${{I}_1} \times {{I}_2} \times {{I}_3} \times \cdots {{I}_N}$| in the complex Euclidean space, respectively. Taking matrices as an example, the low-rank structure of a matrix is defined by the singular values' sparsity. The problem of low-rank matrix approximation can typically be formulated as

where |$f({{\bf X}})$| represents the data approximation term; |${{\| {{\bf X}} \|}_*}{\rm{ = }}\sum\nolimits_i {{{\sigma }_i}( {{\bf X}} )} $| is the nuclear norm of matrix |${{\bf X}}$|⁠, which is the low-rank regularization term; and |$\lambda $| is the regularization parameter. Low-rank regularization terms are generally divided into two types: (i) applying different types of constraint to the singular values obtained from the SVD decomposition of the matrix; and (ii) applying constraints to the factor matrices obtained through bilinear factorization of the matrix; as this method is more efficient in solving problems, this paper focuses on the second type of method. Based on bilinear factorization (as shown in Fig. 1a), Equation (1) can be rewritten as

where |${{\bf U}} \in$|£|${^{m \times r}},{{\bf V}} \in$|£|${^{n \times r}}$| is the factor matrix and |$\| {*} \|_F^{}\,{\rm{ = }}\,\sqrt {\sum\nolimits_{i = 1}^m {\sum\nolimits_{j = 1}^n {{{{| {{{a}_{ij}}} |}}^2}} } } $| denotes the Frobenius norm of the matrix. Since the nuclear norm can be equivalently written as |${{\| {{\bf X}} \|}_*} = \mathop {\min }\limits_{{{\bf X}} = {{\bf U}}{{{{\bf V}}}^T}} \frac{1}{2}( {\| {{\bf U}} \|_F^2 + \| {{\bf V}} \|_F^2} )$|⁠, the problems expressed by Equations (1) and (2) are equivalent. When data |${{\bf Y}} \in$|£|${^{m \times n}}$| is missing (assuming it is randomly missing) and contaminated with noise, the missing data can be effectively recovered through the low-rank matrix completion method as expressed by Equation (2).

Figure 1.

(a) Bilinear factorization of a matrix. (b) Matrix completion for randomly missing data.

Open in new tabDownload slide

2.2. Low-rank tensor completion using Frobenius/nuclear mixed norm

In practical data processing, low-rank matrix/tensor completion methods are commonly employed to reconstruct data following NMO. This choice arises from the transformation of field data, especially high-dimensional seismic data, which initially exhibit hyperbolic moveout relationship concerning reflection events along the offset distance directions |${{H}_x}$| and |${{H}_y}$|⁠. These relationships evolve into linear or parabolic events post-transformation. Consequently, the previously coupled dimensions in the frequency-space domain become separable, thereby reducing the rank of seismic data and facilitating a more straightforward selection of tensor rank. As detailed in Section 2.1, for 5D seismic data, regularization constraints can be imposed on the mode-|$(m,n)$| unfolding matrix of the data tensor (Fig. 2a) to reconstruct data |$\mathcal{X}$| with a low-rank structure from missing and noisy data |$\mathcal{Y}$|⁠:

Figure 2.

(a) Schematic diagram illustrating the mode (m, n) unfolding matrix of a fourth-order tensor. (b) Schematic diagram depicting the bilinear factor decomposition of the mode (m, n) unfolding matrix of a fourth-order tensor.

Open in new tabDownload slide

In the equation, P represents the sampling operator, and |${{{{\bf X}}}_{[n]}} = {{{{\bf U}}}_n}{{\bf V}}_n^T$| denotes the bilinear factor matrix decomposition of the mode-|$(m,n)$| unfolding matrix of the fourth-order tensor |$\mathcal{X}$| (Fig. 2b). Since tensor unfolding into an approximate square matrix is feasible only when the data dimension |$N \ge 4$|⁠, Equation (3) can be utilized to model 5D seismic data (where data reconstruction operates on 4D frequency slice data), thus facilitating the reconstruction and denoising of high-dimensional seismic data.

When the subsurface medium is complex, the CMP data after NMO correction may not adhere to the linear or quasi-linear assumptions within local data windows, thereby affecting the efficacy of the nuclear norm approximation to the rank function. To mitigate this challenge, a non-convex Frobenius/nuclear mixed norm minimization model has been proposed (Shang et al. 2018). This model is presumed to offer a superior approximation to the rank function and produce improved outcomes when completing or restoring low-rank matrices. For any matrix |${{\bf X}} \in {{\mathbb{R}}^{m \times n}}$| of rank R, the Frobenius/nuclear mixed norm is defined as

Consequently, the |${{\| {{\bf X}} \|}_{F - N}}$| norm can be utilized to replace the nuclear norm approximation function in Equation (3). Thus, the low-rank tensor completion (Frobenius-nuclear mixed norm constraint, FN-TC) problem under the constraint of FN mixed norms for high-dimensional seismic data reconstruction and denoising can be formulated as follows:

To solve the constrained optimization problem described in Equation (5), auxiliary variables |${{{{\hat{\bf V}}}}_n}$| and Lagrange multiplier variables |${{{{\bf D}}}_n}$| and |${{{{\bf F}}}_n}$| are introduced |$1 \le n \le 3$|⁠. By converting the equality-constrained optimization problem into an unconstrained optimization problem, we derive the augmented Lagrangian function for Equation (5) as follows:

The alternating direction multiplier method (ADMM) algorithm can be employed to solve the unconstrained optimization problem presented in Equation (6). This iterative approach enables the determination of variables |${{{{\bf U}}}_n},{{{{\bf V}}}_n},{{{{\hat{\bf V}}}}_n},\mathcal{X},{{{{\bf D}}}_n},{{{{\bf E}}}_n}$|⁠. Next, a comprehensive explanation of the solution procedure for each variable in Equation (6) will be provided.

• The subproblem in Equation (6) concerning the factor matrix |${{{{\bf U}}}_n}$| is given by

  $$\begin{eqnarray} \mathop {\min }\limits_{{{{{\bf U}}}_n}} \frac{1}{3}\left\| {{{{{\bf U}}}_n}} \right\|_F^2 + \frac{\mu }{2}\left\| {{{{{\bf X}}}_{[n]}} - {{{{\bf U}}}_n}{{\bf V}}_n^T + {{{{{{\bf D}}}_n}} \mathord{\left/ {\vphantom {{{{{{\bf D}}}_n}} \mu }} \right. } \mu }} \right\|_F^2.\end{eqnarray}$$

  (7)

The least squares solution for the matrix |${{{{\bf U}}}_n}$| is given by:

• The formulation of the subproblem associated with the factor matrix |${{{{\bf V}}}_n}$| within Equation (6) is as follows:

  $$\begin{eqnarray} \mathop {\min }\limits_{{{{{\bf V}}}_n}} \frac{1}{2}\left\| {{{{{\bf X}}}_{[n]}} - {{{{\bf U}}}_n}{{\bf V}}_n^T + {{{{{{\bf D}}}_n}} \mathord{\left/ {\vphantom {{{{{{\bf D}}}_n}} \mu }} \right. } \mu }} \right\|_F^2 + \frac{1}{2}\left\| {{{{{\bf V}}}_n} - {{{{{\hat{\bf V}}}}}_n} + {{{{{{\bf F}}}_n}} \mathord{\left/ {\vphantom {{{{{{\bf F}}}_n}} \mu }} \right. } \mu }} \right\|_F^2.\end{eqnarray}$$

  (9)

The solution for the matrix |${{{{\bf V}}}_n}$| via least squares is provided as follows:

• With respect to the subproblem in Equation (6) concerning the matrix |${{{{\hat{\bf V}}}}_n}$|⁠, it can be articulated as

  $$\begin{eqnarray} \mathop {\min }\limits_{{{{{{\hat{\bf V}}}}}_n}} \frac{2}{3}{{\left\| {{{{{{\hat{\bf V}}}}}_n}} \right\|}_ * } + \frac{\mu }{2}\left\| {{{{{\bf V}}}_n} - {{{{{\hat{\bf V}}}}}_n} + {{{{{{\bf E}}}_n}} \mathord{\left/ {\vphantom {{{{{{\bf E}}}_n}} \mu }} \right. } \mu }} \right\|_F^2.\end{eqnarray}$$

  (11)

The solution for the matrix |${{{{\hat{\bf V}}}}_n}$| is given by

where the symbol |${{\mathcal{T}}_{{2 \mathord{/ {\vphantom {2 {3\mu }}} } {3\mu }}}}( x )$| denotes the soft thresholding operator, defined as follows:

• The subproblem for the tensor |$\mathcal{X}$| in Equation (6) is formulated as

  $$\begin{eqnarray} \mathop {\min }\limits_{\mathcal{X}} \frac{\lambda }{2}\left\| {\mathcal{P} * \mathcal{X} - \mathcal{Y}} \right\|_F^2 + \frac{\mu }{2}\sum\limits_{n = 1}^N {\left\| {{{{{\bf X}}}_{[n]}} - {{{{\bf U}}}_n}{{\bf V}}_n^T + {{{{{{\bf D}}}_n}} \mathord{\left/ {\vphantom {{{{{{\bf D}}}_n}} \mu }} \right. } \mu }} \right\|_F^2} .\end{eqnarray}$$

  (14)

The least squares solution for the tensor |$\mathcal{X}$| is determined by

In the given context, the tensor |${{\mathcal{Z}}_n}$| is referred to as |$mat2ten( {\mu {{{{\bf U}}}_n}{{\bf V}}_n^T - {{{{\bf D}}}_n},{{{{\bf j}}}_n}} )$|⁠. The operation |$mat2ten( {{{\bf X}},{{{{\bf j}}}_n}} )$| signifies the process of folding matrix |${{\bf X}}$| along mode-|${{{{\bf j}}}_n}$|⁠, which is the converse of the mode-|${{{{\bf j}}}_n}$| matrix unfolding for tensor |$\mathcal{X}$|⁠.

• The update formulas for matrices |${{{{\bf D}}}_n}$| and |${{{{\bf F}}}_n}$|⁠, along with the variable |$\mu $|⁠, are provided next:

The parameter |$\rho $| is selected from the interval |$(1.0,1.8]$|⁠, and the initial value of |$\mu $| is typically set to 0.01, with |${{\mu }_{\max }}$| being set to 10⁸. Table 1 provides a detailed description of the update process for each variable in the FN-TC method.

It is noteworthy that the norm of each bilinear factor matrix in Equation (5) is convex, simplifying handling and extension in contrast to the Schatten p norm (0 < p < 1). Additionally, to ensure simplicity and clarity in the formula derivation, explicit assumptions about noise distribution are omitted. In practical scenarios, seismic data are typically represented as |${\mathcal{Y}} = {\mathcal{X}} + {\mathcal{E}}$|⁠, where |$\boldsymbol{\mathcal{E}}$| signifies noise. Different regularization constraints can be applied based on the probabilistic distribution characteristics of noise |$\boldsymbol{\mathcal{E}}$|⁠. For example, anomalous noise |${{\| \boldsymbol{\mathcal{E}} \|}_{{{L}_1}}}$| following a Laplace distribution and random noise |${{\| \boldsymbol{\mathcal{E}} \|}_{{{L}_2}}}$| following a Gaussian distribution can both be solved using ADMM algorithm. This methodology enhances the algorithm's robustness.

3. Examples

This section confirms the effectiveness of the method introduced in this article by using synthetic and field data for reconstruction and noise suppression. The synthetic data is generated by convolving a Ricker wavelet, which has a dominant frequency of 30 Hz, and reflection coefficients. Subsequently, the outcomes of the reconstruction and noise suppression are assessed using

Here, |${{{\mathcal{X}}}_{\textit{true}}}$| represents the data that is free of noise, and |${{{\mathcal{X}}}_{est}}$| stands for the reconstructed data (noise-reduced data).

The PSMF method conducts bilinear factor decomposition on the mode |$(m,n)$| unfolding matrix of the 4D tensor, followed by solving the factor matrix utilizing the alternating least squares algorithm. Conversely, the FN-TC method, introduced in this paper, imposes a non-convex low-rank constraint on the unfolded matrix and subsequently solves each variable through the ADMM algorithm. On the other hand, the HOSVD method derives the factor matrix by performing SVD on each mode n unfolding matrix of the 4D tensor. Therefore, the PSMF and HOSVD methods, being representative, are chosen as comparative approaches to evaluate the efficacy of the FN-TC method proposed in this paper.

3.1. Theoretical numerical examples

In this study, we analyze 5D synthetic seismic data, which comprises a |$301 \times 16 \times 16 \times 16 \times 16$| grid. The time sampling interval of the data is 2 ms, with a sampling interval of 12.5 m in every direction. We proceed with the assumption that the data are affected by Gaussian white noise. The synthetic data is theoretically composed of four events: two parabolic and two linear. These events are characterized by corresponding amplitude values of 1.8, 1.2, 1.0, and 0.6, respectively. Specifically, the dataset exhibits parabolic events in the direction of offset distance and linear events in the direction of CMP. Figure 3 depicts the distribution of normalized singular values of the matrices resulting from unfolding along various modes for the 30 Hz frequency slice data. Remarkably, the matrices unfolded along mode |$(m,n)$| of the 4D tensor consistently exhibit low-rank characteristics, with all unfolded matrices possessing a rank of 4. Initially, we evaluated the reconstruction efficacy of the HOSVD, PSMF, and FN-TC methods without noise, while varying the proportion of missing data (ranging from 20% to 90%). The rank and iteration number (Niter) were set to 4 and 30, respectively, for all three methods. In the PSMF method, the parameter |$\alpha $| was fixed at 1, whereas in the FN-TC method, the regularization parameter |$\lambda $| was denoted as |${{10}^3}$|⁠. Figure 4 illustrates the data reconstruction outcomes under noise-free conditions across different missing ratios. It is evident that as the missing ratio increases to 90%, the reconstruction quality notably deteriorates for both the HOSVD and PSMF methods, while the FN-TC method exhibits a slower decline. This suggests that all three methods necessitate additional iterations to converge.

Figure 3.

The normalized singular value distribution of matrices after unfolding the 30 Hz frequency slice data along different modes. (a) Singular value distribution of the mode n unfolding matrix; and (b) singular value distribution of the mode (m, n) unfolding matrix.

Open in new tabDownload slide

Figure 4.

The data reconstruction results using HOSVD, PSMF, and FN-TC methods under noise-free conditions and with varying percentages of missing data.

Open in new tabDownload slide

For the reconstruction of high-dimensional data with 60% trace missing and contaminated by −1dB random noise, the PSMF method's parameter |$\alpha $| and iteration number |${{N}_{\textit{iter}}}$| were set to 0.05 and 200, respectively. Meanwhile, the iteration number |${{N}_{\textit{iter}}}$| for both the HOSVD and FN-TC methods was set to 30. Additionally, the regularization parameter |$\lambda $| for the FN-TC method was set to 0.1. In Fig. 5, the reconstruction results on noisy missing data for the three methods are presented. Analysis of the reconstruction outcomes and data residual profiles reveals that the FN-TC method demonstrates superior noise suppression capabilities with minimal distortion to the effective signal, characterized by continuous events. Numerically, the FN-TC method achieves the highest data reconstruction SNR (24.5 dB) among the three methods, with a runtime of 34 s, which is 12 s longer than that of the HOSVD method and times longer than that of the PSMF method. Figure 6a showcases the reconstruction results of the HOSVD, PSMF, and FN-TC methods on different frequency slice data, while Fig. 6b displays the relative error convergence curves for the reconstruction of 30Hz frequency slice data by the three methods. It is observed that the FN-TC method surpasses the other two methods in reconstruction results across various frequency slice data. Furthermore, both the HOSVD and FN-TC methods achieve convergence after several iterations, whereas the PSMF method requires nearly 200 iterations to meet the convergence criteria, thereby explaining its longer runtime.

Figure 5.

In the comparison of HOSVD, PSMF, and FN-TC methods for reconstructing noisy and incomplete data composed of linear and parabolic common reflection axes, the following results were obtained: (a) a 3D subset of the original data after being preprocessed and fixed from the 5D dataset. (b) The input data with 60% missing traces and contaminated by −1 dB of random noise. (c) The reconstruction result using the HOSVD method, achieving a SNR of 20.9 dB with a processing time of 22 s. (d) The data residual between the result in (c) and the original data in (a). (e) The reconstruction result using the PSMF method, achieving an SNR of 22.9 dB with a processing time of 112 s. (f) The data residual between the result in (e) and the original data in (a). (g) The reconstruction result using the FN-TC method, achieving an SNR of 24.5 dB with a processing time of 34 s. (h) The data residual between the result in (g) and the original data in (a). Note that the color scale for data residual results (d), (f), and (h) ranges from −0.1 to 0.1, while the color scales for the other figures range from −1 to 1.

Open in new tabDownload slide

Figure 6.

(a) Reconstruction results of HOSVD, PSMF, and FN-TC methods at different frequencies. (b) Relative error convergence curves for the reconstruction of 30 Hz frequency slice data by the three methods, with the green and red dashed lines representing the relative errors of HOSVD and FN-TC methods at the 30th iteration, respectively.

Open in new tabDownload slide

The second numerical experiment aimed to assess the data reconstruction performance of the FN-TC method under conditions involving abnormal and random noise, along with data missingness. The process of introducing abnormal noise to individual trace data is detailed as follows: for each seismic trace, randomly select m sampling points, where m is a random number; m values are then generated from a uniform distribution function within the interval |$( - 1,1)$|⁠; these values are scaled by six times the maximum amplitude to yield reflection coefficients with abnormal values; subsequently, they are convolved with the Ricker wavelet to generate abnormal noise. This numerical experiment was conducted on the first synthetic dataset, consisting of four parabolic and linear events. To simulate data with abnormal noise, Gaussian white noise with a −1 dB intensity was initially added. Then, abnormal noise was introduced to 20% of the seismic traces, followed by random missingness affecting 30% of the seismic traces. During the reconstruction of the synthetic data with abnormal noise, the iteration number |${{N}_{\textit{iter}}}$| for the FN-TC method was set to 100, while the regularization parameters |$\beta $| and |$\lambda $| were both set to 1. Figure 7 showcases the data reconstruction outcomes in the CMP domain using the FN-TC method (in waveform display), achieving SNR = 23.4 dB after the reconstruction process. In Fig. 8a, the reconstruction results of the FN-TC method on various frequency slice data are presented, whereas Fig. 8b illustrates the relative error convergence curve for the reconstruction of the 30 Hz frequency slice data by the FN-TC method. It is noteworthy that the relative error curve does not exhibit monotonic decrease with increasing iterations, but eventually converges to a stable point.

Figure 7.

Data reconstruction with abnormal and random noise and missing data. (a) Illustration of original data slices within the CMP domain (⁠|$CMPx = 8,\textit{CMPy} = 8,Hy = 2,5,8,10$|⁠); (b) input data contaminated with −1 dB random noise, 20% of seismic traces affected by abnormal noise, and 30% of seismic traces randomly missing; (c) reconstruction result using the FN-TC method, with SNR = 23.4 dB; and (d) error analysis between the reconstruction (c) and the original data (a).

Open in new tabDownload slide

Figure 8.

(a) Results of the FN-TC method's reconstruction at various frequencies. (b) Relative error convergence curve for the FN-TC method when reconstructing 30 Hz frequency slice data.

Open in new tabDownload slide

3.2. Field data examples

This subsection aims to evaluate the reconstruction performance of the FN-TC and PSMF methods using land-based 5D CMP data. In Fig. 9a, the spatial distribution of shot and receiver points for the 5D CMP data is presented, while Fig. 9b illustrates the frequency coverage at each CMP location. The field-collected CMP data covers a range of ∼−3000–2000 m in the x-direction of offset and ∼−2500–2500 m in the y-direction of offset. After NMO, irregular data is processed by binning onto a regular grid with grid spacings of 6.25 m in the inline and crossline directions, and 250 and 500 m in the x and y offsets, respectively. The binned seismic data forms a regular 5D tensor of size |$401 \times 201 \times 20 \times 24 \times 10$|⁠, with a time sampling interval of 4 ms, and 25% of the seismic traces containing missing data. For data reconstruction, the field data is generally divided into blocks for processing. Typically, under the condition of satisfying theoretical assumptions, larger block sizes yield better reconstruction quality, albeit at the cost of increased processing time. Meanwhile, for the entire dataset, the total time consumption equals the number of blocks multiplied by the time consumption of a single block. Therefore, in practical applications, a balance between effectiveness and efficiency is necessary before providing the criteria for selecting block size and stride. For the field data in this paper, a block size of |$401 \times 20 \times 20 \times 17 \times 10$| is chosen, with the data window moving with steps of 1, 15, 1, 12, and 1 along each dimension. This process results in the reconstruction of 28 data blocks, which are then averaged to obtain the final reconstruction result. In the FN-TC method, the rank R, regularization parameter |$\lambda $|⁠, and the number of iterations |${{N}_{\textit{iter}}}$| are set to 8, 0.1, and 30, respectively. In the PSMF method, the rank R and the number of iterations |${{N}_{\textit{iter}}}$| are set to 8 and 100, respectively.

Figure 9.

(a) The geographical arrangement of the shot points and receiver points in the field data is depicted with red dots indicating shot points and blue dots indicating receiver points. (b) The coverage counts for each CMP point.

Open in new tabDownload slide

Figures 10 and 11 depict the data reconstruction outcomes before and after applying the FN-TC and PSMF methods in the CMP domain. In Fig. 10a, a partial slice (⁠|$Hy = 4,5,6,7$|⁠) of a single CMP gather (⁠|$CMPx = 130,\textit{CMPy} = 10$|⁠) is displayed; Fig. 10 parts b and c showcase the reconstruction results using the FN-TC and PSMF methods, respectively. It is evident that both methods effectively reconstruct the field data, with strong energy events appearing more continuous and maintaining consistent energy levels pre- and post-reconstruction. Notably, for weak energy events, the FN-TC method better preserves continuity compared to the PSMF method. Moving to Fig. 11, the reconstruction results for two common offset profiles (⁠|$Hx = 10,Hy = 3,\textit{CMPy} = 8$| and |$Hx = 12,Hy = 5,\textit{CMPy} = 18$|⁠) are presented. The data reconstruction results of the FN-TC and PSMF methods are similar in most areas, but the PSMF method's data reconstruction profile exhibits more noise and discontinuities in the event sequences, as indicated by the red circles in the figure. Additionally, processing 28 data blocks with the PSMF method takes ∼56.4 minutes, while the FN-TC method's runtime is approximately only one-fifth of the PSMF method, taking ∼12.8 minutes. This analysis suggests that the FN-TC method outperforms the PSMF method in effectively reconstructing high-dimensional data.

Figure 10.

The results showing the reconstruction of field data on the CMP gather using both FN-TC and PSMF methods are as follows: (a) the slice results of the input data at a single CMP point (⁠|$CMPx = 130,\textit{CMPy} = 10,Hy = 4,5,6,7$|⁠); (b) the reconstruction result using the FN-TC method; and (c) the reconstruction result using the PSMF method.

Open in new tabDownload slide

Figure 11.

The results demonstrating the reconstruction of field data on the common offset gather using the FN-TC and PSMF methods before and after the process. (a) The single slice result (⁠|$CMPy = 8$|⁠) of the input data at the common offset point (⁠|$Hx = 10,Hy = 3$|⁠). (b) The reconstruction result after applying the FN-TC method. (c) The reconstruction result after applying the PSMF method.

Open in new tabDownload slide

4. Discussion

When reconstructing high-dimensional seismic data using the FN-TC method, four parameters need to be determined beforehand: the ranks |$({{R}_1},{{R}_2},{{R}_3})$| of each mode-|$(m,n)$| unfolding matrix and the regularization parameter |$\lambda $|⁠. This section investigates the impact of selecting different ranks |$({{R}_1},{{R}_2},{{R}_3})$| and parameters |$\lambda $| on the data reconstruction results. The numerical experiment is based on the first synthetic data (Fig. 5a), where the ranks |$({{R}_1},{{R}_2},{{R}_3})$| of each mode-|$(m,n)$| unfolding matrix were all set to 4. Thus, in the experiments, the ranks |${{R}_n}$| of each unfolding matrix were all set equal to R. First, Gaussian white noise with −1 dB was added to the data, followed by randomly removing 60% of the seismic traces. Figure 12a displays the data reconstruction results of the HOSVD, PSMF, and FN-TC methods with different ranks R and a fixed regularization parameter |$\lambda = 0.1$|⁠. From the figure, we see that the data reconstruction quality of both the HOSVD and PSMF methods decreases as the rank R increases, but the reconstruction quality of the HOSVD method deteriorates more rapidly (from 21 to 10 dB). However, the reconstruction quality of the FN-TC method does not exhibit significant changes, indicating its robustness to the selection of the rank R. Therefore, it is advisable to choose a larger rank R for data reconstruction in field data processing. Figure 12b illustrates the data reconstruction results of the FN-TC method with different regularization parameters |$\lambda $| and a fixed rank |$R = 4$|⁠. It can be seen from the figure that good reconstruction results can be obtained within a relatively large range of |$\lambda $| values (⁠|$0.1 \sim {{10}^3}$|⁠). From thisdiscussion, it can be concluded that the proposed FN-TC method exhibits robustness in parameter selection.

Figure 12.

For data with 60% missing seismic traces and contaminated by random noise with −1 dB, the reconstruction results of the HOSVD, PSMF, and FN-TC methods under different choices of rank R and regularization parameter |$\lambda $| are compared. (a) Reconstruction results with different ranks R. (b) Reconstruction results with different regularization parameters |$\lambda $|⁠.

Open in new tabDownload slide

5. Conclusions

In practical 3D seismic data acquisition, a 5D data volume is obtained, allowing for the selection of seismic data with varying dimensional sizes for modeling purposes. With increased data in dimension, the spatially structural information becomes richer, facilitating more precise signal modeling in high-dimensional space to accurately demonstrate the characteristics of high-dimensional seismic data. Hence, within the Bayesian framework and leveraging the low-rank nature of the data along with prior knowledge specific to seismic signals, this paper introduces a high-dimensional seismic data modeling approach termed the FN-TC. This method imposes non-convex FN mixed norm constraints on the mode (m, n) unfolding matrix of the 4D data tensor in the frequency-space domain. This constraint allows for a closer approximation to the rank function of the matrix, thereby enhancing the accuracy of seismic data reconstruction and denoising, consequently providing high-quality data support for fulfilling subsequent seismic processing tasks. In practice, the noise types included in the collected seismic data are highly complex. The method proposed in this paper assumes that the noise follows a Gaussian or Laplace distribution, with the L1 or L2 norm utilized for measurement. Therefore, an avenue for future research direction would involve analyzing the types of noise in seismic data and further modeling the noise through using a generalized Gaussian mixture model to handle more complicated scenarios in noise attenuation.

Conflict of interest statement. None declared.

Funding

This study is supported by the Joint Funds of the National Natural Science Foundation of China (grant no. U23B6010).

Data availability

There is no data availability information.

References