A five-dimensional seismic data reconstruction method based on preconditioned riemannian gradient descent

Abstract

This invention discloses a five-dimensional seismic data reconstruction method based on preconditioned Riemann gradient descent, belonging to the field of seismic data processing technology. The method involves performing a Hankel transform on fixed-frequency four-dimensional seismic data to construct a fourth-order block Hankel matrix; performing hard thresholding on the fourth-order block Hankel matrix to obtain a low-rank approximation matrix, and calculating the Euclidean gradient at the low-rank approximation matrix in the current iteration; calculating a preconditioner based on the diagonal part of the Euclidean gradient outer product; calculating the projection of the Euclidean gradient onto the tangent space of the Riemann manifold based on the Euclidean gradient and the preconditioner, and updating the low-rank approximation matrix based on the projected gradient to obtain the low-rank approximation matrix for the next iteration; iterative updates continue until a termination condition is met; the low-rank approximation matrix from the last iteration is inversely transformed into a fourth-order tensor, and the reconstructed five-dimensional seismic data is obtained based on the fourth-order tensor and the frequency components in the five-dimensional seismic data. This method can improve the reconstruction efficiency and quality of five-dimensional seismic data.

Description

Technical Field

[0001] This application relates to the field of seismic data processing technology, and in particular to a five-dimensional seismic data reconstruction method based on preconditional Riemann gradient descent. Background Technology

[0002] The quality of seismic exploration data directly determines the imaging accuracy of subsequent subsurface media and the reliability of geological interpretation. Seismic data acquisition is often limited by factors such as surface conditions, equipment deployment, and acquisition costs, leading to problems such as data loss and noise interference. These factors result in insufficient spatial sampling and low signal-to-noise ratios in the acquired seismic data, thus affecting imaging accuracy and the reliability of geological interpretation. Therefore, how to effectively recover seismic data under limited acquisition conditions and develop a seismic data reconstruction technology that combines high robustness and high computational efficiency has become an urgent problem to be solved in the field of seismic data processing.

[0003] As exploration targets extend deeper and into complex tectonic zones, traditional low-dimensional data is no longer sufficient to meet the demands of high-resolution imaging. High-dimensional seismic data possesses stronger spatial correlations; therefore, five-dimensional seismic data extracted from three-dimensional seismic data can provide richer wavefield information across multiple dimensions, including time, space, offset, and azimuth. Reconstruction based on five-dimensional seismic data can achieve better reconstruction quality.

[0004] Complete, noise-free seismic data is mathematically low-ranked. Noise and missing data increase the rank of the seismic data matrix or tensor. Therefore, existing technologies mainly employ low-rank methods for seismic data denoising and reconstruction. However, while traditional reconstruction methods based on the low-rank assumption are effective, they often rely on computationally expensive large-scale singular value decomposition (SVD), which limits their processing efficiency. Summary of the Invention

[0005] Therefore, it is necessary to provide a five-dimensional seismic data reconstruction method based on preconditional Riemann gradient descent to address the aforementioned technical problems.

[0006] The following technical solution is adopted in this specification: This specification provides a five-dimensional seismic data reconstruction method based on preconditional Riemann gradient descent, including: Five-dimensional seismic data in frequency domain representation is obtained, and fixed-frequency slices are made from the five-dimensional seismic data to obtain four-dimensional seismic data. The four-dimensional seismic data is subjected to Hankel transformation to obtain a fourth-order block Hankel matrix; Hard thresholding is performed on the fourth-order block Hankel matrix to obtain a low-rank approximate matrix; The Euclidean gradient at the low-rank approximation matrix at the current iteration point is calculated based on the fourth-order block Hankel matrix, the low-rank approximation matrix, and the sampling operator. The Euclidean gradient represents the residual between the low-rank approximation matrix and the fourth-order block Hankel matrix. The sampling operator represents the linear operator that selects known observation data from the fourth-order block Hankel matrix. The first preconditioner is determined by the diagonal matrix formed by the outer product of the row vectors of the Euclidean gradient, and the second preconditioner is determined by the diagonal matrix formed by the outer product of the column vectors of the Euclidean gradient. Based on the first and second preconditioners, the Euclidean gradient is calculated to a value of rank r and magnitude r. × The projection of the matrix into the tangent space of the Riemann manifold is obtained, and the low-rank approximation matrix is ​​updated according to the gradient after projection. The updated result is then mapped back to the Riemann manifold to obtain the low-rank approximation matrix for the next iteration. The low-rank approximation matrix is ​​updated again for the next iteration until the square of the Frobenius norm of the difference matrix between the low-rank approximation matrices of two adjacent iterations is less than or equal to the iteration stopping error or the maximum number of iterations is reached. The low-rank approximation matrix from the last iteration is transformed into a fourth-order tensor using the inverse Hankel transformation. Based on the fourth-order tensor and the frequency components in the five-dimensional seismic data, the reconstructed five-dimensional seismic data is determined.

[0007] Optionally, the Euclidean gradient at the low-rank approximation matrix in the current iteration is calculated based on the fourth-order block Hankel matrix, the low-rank approximation matrix, and the sampling operator, including: Calculate the product of the low-rank approximation matrix of the current iteration and the sampling operator, and determine the difference between this product and the fourth-order block Hankel matrix as the Euclidean gradient at the low-rank approximation matrix of the current iteration.

[0008] Optionally, a first preconditioner is determined based on the diagonal matrix formed by the outer products of the row vectors of the Euclidean gradient, and a second preconditioner is determined based on the diagonal matrix formed by the outer products of the column vectors of the Euclidean gradient, including: Calculate the product of the Euclidean gradient and its transpose to obtain the first matrix, and extract the diagonal elements of the first matrix to form the first diagonal matrix. The first preconditioner is obtained based on the first diagonal matrix and the identity matrix weighted by the regularization parameter; Calculate the product of the transpose of the Euclidean gradient and itself to obtain the second matrix, and extract the diagonal elements of the second matrix to form a second diagonal matrix; The second preconditioner is obtained from the second diagonal matrix and the identity matrix weighted by the regularization parameter.

[0009] Optionally, based on the first and second preconditioners, the projection of the Euclidean gradient onto the tangent space of the Riemannian manifold composed of low-rank matrices of rank r is calculated, including: The weighted inner product of the Euclidean gradient is calculated using the first and second preconditioners. Calculate the rank of the Euclidean gradient under the weighted inner product as r, with a magnitude of r. × The projection of the tangent space of the Riemannian manifold composed of matrices.

[0010] Optionally, the low-rank approximation matrix is ​​updated based on the projected gradient, including: Given a fixed constant step size and the rank of the Euclidean gradient under the weighted inner product as r, and a magnitude of... n 1× n The tangent space projection of the Riemannian manifold composed of matrices of size 2 is used to update the low-rank approximation matrix, yielding the updated result; the update formula for the low-rank approximation matrix is: in, Indicates the first The update result of the low-rank approximation matrix in the next iteration. Indicates the first The low-rank approximation matrix of the next iteration. Indicates a fixed constant step size. This represents the projection of the Euclidean gradient onto the tangent space under the weighted inner product. This represents the tangent space projection operator under the weighted inner product.

[0011] Optionally, the update result is mapped back to the Riemannian manifold to obtain the low-rank approximation matrix for the next iteration, including: Perform a hard threshold operation on the update result of the current iteration to obtain the low-rank approximation matrix for the next iteration.

[0012] Optionally, five-dimensional seismic data in frequency domain representation are acquired, including: Acquire noisy and missing five-dimensional seismic data in the time domain; The five-dimensional seismic data in the time domain, which contains noise and has missing data, is subjected to Fourier transform to obtain the five-dimensional seismic data in the frequency domain.

[0013] Optionally, the four-dimensional seismic data is subjected to a Hankel transform to obtain a fourth-order block Hankel matrix, including: Four-dimensional seismic data is embedded into a first-order block Hankel matrix using all components of the first dimension of the four-dimensional seismic data. The first-order Hankel matrix is ​​embedded into the second-order block Hankel matrix using all components of the second dimension in the four-dimensional seismic data. The second-order Hankel matrix is ​​embedded into a third-order block Hankel matrix using all components of the third dimension in the four-dimensional seismic data. The third-order Hankel matrix is ​​embedded into a fourth-order block Hankel matrix using all components of the fourth dimension in the four-dimensional seismic data.

[0014] Optionally, the reconstructed five-dimensional seismic data is determined based on the fourth-order tensor and the frequency components in the five-dimensional seismic data, including: Based on the frequency components in the fourth-order tensor and five-dimensional seismic data, the reconstructed frequency domain seismic data is determined; The reconstructed frequency domain seismic data is subjected to inverse Fourier transform to recover the complete time domain seismic data, i.e., the reconstructed five-dimensional seismic data.

[0015] This specification provides a five-dimensional seismic data reconstruction device based on preconditional Riemann gradient descent, including: The acquisition module is used to acquire five-dimensional seismic data in frequency domain representation and to slice the five-dimensional seismic data at fixed frequencies to obtain four-dimensional seismic data. The transformation module is used to perform Hankel transformation on four-dimensional seismic data to obtain a fourth-order block Hankel matrix. The computation module is used to perform hard thresholding on the fourth-order block Hankel matrix to obtain a low-rank approximate matrix. The iterative module is used in any iteration of the preconditional Riemann gradient descent method to calculate the Euclidean gradient matrix at the low-rank approximation matrix in the current iteration, based on the fourth-order block Hankel matrix, the low-rank approximation matrix, and the sampling operator. The Euclidean gradient matrix represents the residual between the low-rank approximation matrix and the fourth-order block Hankel matrix. The sampling operator represents a linear operator that selects known observation data from the fourth-order block Hankel matrix. The first preconditioner is determined based on the diagonal matrix formed by the outer products of the row vectors of the Euclidean gradient matrix, and the second preconditioner is determined based on the diagonal matrix formed by the outer products of the column vectors of the Euclidean gradient matrix. Based on the first and second preconditioners, the Euclidean gradient is calculated to a value of rank r and size r. × The projection of the matrix into the tangent space of the Riemann manifold is obtained, and the low-rank approximation matrix is ​​updated according to the gradient after projection. The updated result is then mapped back to the Riemann manifold to obtain the low-rank approximation matrix for the next iteration. The low-rank approximation matrix for the next iteration is updated again until the square of the Frobenius norm of the difference matrix between the low-rank approximation matrices of two adjacent iterations is less than or equal to the iteration stopping error or the maximum number of iterations is reached. The determination module is used to transform the low-rank approximation matrix of the last iteration into an inverse Hankel transformation into a fourth-order tensor, and to determine the reconstructed five-dimensional seismic data based on the fourth-order tensor and the frequency components in the five-dimensional seismic data.

[0016] This specification provides a computer-readable storage medium storing a computer program that, when executed by a processor, implements the above-described five-dimensional seismic data reconstruction method based on preconditional Riemann gradient descent.

[0017] This specification provides a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the program to implement the above-described five-dimensional seismic data reconstruction method based on preconditional Riemann gradient descent.

[0018] The above-mentioned technical solutions adopted in this specification can achieve the following beneficial effects: The five-dimensional seismic data reconstruction method based on preconditional Riemann gradient descent provided in this specification improves the convergence speed of the algorithm and reduces the number of iterations by constructing efficient preconditioners using the rows and columns of the gradient matrix. Furthermore, it avoids large-scale singular value decomposition by using subspace projection. This method can skip the subspace approximation process and directly calculate the subspace projection and the singular value decomposition of the small matrix, thereby significantly reducing the computational cost and improving the reconstruction efficiency of five-dimensional seismic data. Attached Figure Description

[0019] The accompanying drawings, which are included to provide a further understanding of this application and form part of this application, illustrate exemplary embodiments and are used to explain this application, but do not constitute an undue limitation of this application. In the drawings:

[0020] Figure 1 This document provides a flowchart illustrating a five-dimensional seismic data reconstruction method based on preconditional Riemann gradient descent. Figure 2 A comparison of computation time between PRGD and DRR methods for 5D data of different sizes; Figure 3 A comparison of reconstruction performance and reconstruction time between the PRGD and DRR methods under different data missing rates; Figure 4 For a single common center point (CMP) gather ( 5D composite data result comparison chart; Figure 5 For common offset gathers ( Comparison of 5D synthetic data results; Figure 6 Comparison chart of 5D composite data results for common center point gathers of data; Figure 7 Comparison chart of 5D composite data results from common offset gathers; Figure 8 A comparison of reconstruction performance and time for PRGD and DRR methods at different iteration numbers; Figure 9 For a common center point set 3D graphics comparison; Figure 10 A comparison plot of Fk spectra from field data; Figure 11 For fixed A comparison diagram of the two-dimensional unfolded set of common center point; Figure 12 For fixed Local comparison of common offset gathers; Figure 13 A schematic diagram illustrating the selection of step size under various missing rates and signal-to-noise ratio conditions; Figure 14 This specification provides a schematic diagram of a computer device for implementing a five-dimensional seismic data reconstruction method based on preconditional Riemann gradient descent. Detailed Implementation

[0021] To make the objectives, technical solutions, and advantages of this specification clearer, the technical solutions of this application will be clearly and completely described below in conjunction with specific embodiments and corresponding drawings. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of them. All other embodiments obtained by those skilled in the art based on the embodiments in this specification without creative effort are within the scope of protection of this application.

[0022] In recent years, researchers have proposed and developed various theories and methods to address this problem. The main reconstruction methods include seismic data reconstruction techniques based on sparse constraints, low-rank constraints, and deep learning.

[0023] Seismic data reconstruction methods based on sparsity constraints are mainly divided into two categories: reconstruction methods based on sparse transformations and reconstruction methods based on dictionary learning. The former achieves efficient reconstruction by mapping seismic data to a sparse domain, and its key lies in the selection of the sparse transformation. Commonly used transformation methods include Radon transform, Fourier transform, Curvelet transform, Seislet transform, and Dreamlet transform. Typical algorithms based on these transformations include Minimum Weighted Norm Interpolation (MWNI), Leakage-resistant Fourier Transform (ALFT), Projection OntoConvex Sets (POCS), and Matching Pursuit Regularization. Furthermore, researchers have proposed a five-dimensional shot-receiver distance vector patch (OVT) domain seismic data reconstruction method based on the Orthogonal Matching Pursuit (OMP) algorithm, achieving high-precision and high-fidelity data reconstruction. Researchers have also effectively improved the reconstruction performance of non-uniformly sampled seismic data based on the Iterative Soft Thresholding (IST) algorithm by introducing an approximate shrinkage operator based on Taylor series and supplementing it with a redundant Fourier transform inspired by compressed sensing. Sparse reconstruction methods based on dictionary learning can more flexibly capture seismic signal features by adaptively learning sparse representation bases from data. Researchers have proposed a fast dictionary learning method based on continuous generalized K-means (SGK), which significantly improves computational efficiency by replacing the time-consuming singular value decomposition in K-SVD with an arithmetic mean. - Orthogonal dictionary learning method with norm maximization ( This method (MODL) achieves sparse representation under orthogonal dictionary constraints by maximizing the MODL norm and uses the MSP algorithm for efficient dictionary updates, thereby improving the computational efficiency and quality of seismic data reconstruction.

[0024] Complete, noise-free seismic data is mathematically low-ranked; noise and missing data increase the rank of the seismic data matrix or tensor. Therefore, low-rank methods are effective for denoising and reconstructing seismic data. Based on this theory, seismic data reconstruction methods based on low-rank constraints have been developed. Low-rank methods can be divided into two categories. The first category employs dimensionality reduction techniques for multilinear arrays or tensors based on tensor completeness theory. For example, the High-Order Singular Value Decomposition (HOSVD) proposed by researchers. This method directly acts on the seismic data tensor to recover useful signals. Researchers have proposed a parallel matrix factorization (PMF) algorithm for seismic data reconstruction. This technique performs matrix factorization on different tensor expansions, thus avoiding the computation of Singular Value Decomposition (SVD) and significantly reducing computational complexity. For the reconstruction of 5D seismic data contaminated by anomalous noise, researchers have proposed a robust tensor completion algorithm. This algorithm introduces a robust loss function into the parallel matrix factorization framework to replace the traditional L2 norm, thereby significantly reducing the impact of outliers on the reconstruction results. Researchers introduced TT and TR decomposition into the PMF algorithm and accelerated computation using a randomized compressed sampling method, developing a series of PMF-based seismic data reconstruction algorithms. They also introduced Fully Connected Tensor Network (FCTN) decomposition based on frequency slice Hankel tensors for 3D seismic data reconstruction. FCTN decomposes the fourth-order tensor into four factor-shrinkable forms, overcoming the limitations of traditional tensor decomposition methods (such as CANDECOMP / PARAFAC(CP) and Tucker decomposition) which cannot establish connections between different factors and are ineffective in representing relationships. Fully connected tensor networks are also known as Complete Graph Tensor Networks (CGTN). Researchers designed a low-rank method (LRA-LTCGTN) for 5D seismic data reconstruction using CGTN and learnable transformations (LT), which exhibits stronger rank robustness than LRA-CGTN. Finally, researchers proposed a multi-component seismic data vector reconstruction method (QMF) based on quaternion matrix decomposition, which effectively preserves the vector structure characteristics of the seismic wavefield by jointly reconstructing three-component data. Another type of method is seismic data reconstruction based on matrix completeness theory. These methods typically ignore the inherent structure of tensors, reformulating the tensor completion problem as a matrix completion problem. Compared to tensor completion methods, these methods are more mature and easier to implement. A typical method is the Cadzow filter rank reduction technique, also known as singular spectral analysis (SSA). SSA has been extended to multichannel singular spectral analysis (MSSA) for the reconstruction and denoising of multidimensional seismic data. To improve the denoising and reconstruction results of seismic data, researchers have introduced damped rank reduction methods. This method can significantly improve the estimation quality of low-rank signal matrices even under conditions of strong noise and severe missing data. Its core principle is to use damped truncated singular value decomposition to more accurately recover the true signal.Researchers have proposed an Enhanced Low-Rank Matrix Estimation (ELRME) method to address the performance degradation of traditional damped rank reduction (DRR) methods under strong noise and high proportion of missing data. This algorithm constructs a novel proximity operator by combining a moving average filter with an arctangent penalty function, significantly optimizing the rank reduction process. Researchers have also proposed a two-step singular spectrum analysis (SSA) method for robust low-rank approximation and anomalous noise suppression of seismic data. This method significantly improves computational efficiency while enhancing the signal-to-noise ratio through a two-step "prediction-elimination" approach. Furthermore, researchers have proposed an efficient reconstruction method based on Geman function minimization. This method utilizes a non-convex Geman low-rank model (NCGL) to more accurately approximate the rank function and solves it using KKT conditions, achieving efficient seismic data recovery without introducing additional parameters. Finally, researchers have proposed a joint sparse and low-rank prior seismic data reconstruction method (JSLRP). This method models the data reconstruction problem of irregular missing traces as a joint sparse and low-rank matrix approximation problem and solves it efficiently using the alternating direction multiplier method. Researchers have proposed a seismic data reconstruction method based on nonlocal self-similarity (RRGR). This method effectively addresses the unreliable matching problem of similar blocks caused by missing traces through a two-stage iterative framework, combining pre-reconstruction guidance and low-rank matrix completion. Researchers have also proposed an In-the-Domain (I-FMSSA) reconstruction method, which directly handles irregular observation coordinates by introducing interpolation operators, significantly improving computational efficiency while maintaining spatial coordinate accuracy, thus avoiding explicit Hankel matrix construction. Although rank reduction methods are effective in reconstructing seismic data, their computational cost is typically high. To improve computational efficiency, researchers have developed various techniques, such as stochastic singular value decomposition, Lanczos bidiagonalization, and a Hankel matrix-free FMSSA method.

[0025] In recent years, deep learning has developed rapidly in the field of earthquake data reconstruction. Deep learning-based earthquake data reconstruction methods can be divided into two categories: supervised and unsupervised. Supervised deep learning methods require learning prior knowledge from the dataset to reconstruct missing data. Researchers have proposed a supervised learning method using support vector regression to reconstruct five-dimensional earthquake data. Various types of neural networks have also been used for earthquake data reconstruction, including autoencoders, generative adversarial networks, and model-driven networks. Researchers have proposed a five-dimensional convolutional neural network called CCNet-5D, which approximates five-dimensional convolution operations by concatenating two-dimensional and three-dimensional convolutions, thereby enabling the simultaneous extraction of multi-dimensional spatial features and achieving five-dimensional earthquake data reconstruction. Unlike supervised methods, unsupervised methods do not rely on labeled data, but their computational cost is usually high, making it difficult to meet real-time processing requirements. Researchers have proposed an unsupervised three-dimensional earthquake data reconstruction network based on deep priors (DPSI). This method inputs random noise into a multi-scale U-Net network and constrains the output through adaptive weighted Laplacian regularization to achieve unsupervised three-dimensional earthquake data reconstruction. Researchers have proposed a Noise2Void-inspired self-supervised reconstruction network (SSLI) that achieves high-precision reconstruction of 2D and 3D seismic data by extracting partial samples from missing data as labels and training a U-Net network. They have also proposed a seismic data reconstruction method combining unsupervised deep learning and Monte Carlo Dropout. Furthermore, they have proposed a frequency-space dependent unsupervised deep learning framework called FUDLInter, which achieves high-precision reconstruction of 3D and 5D seismic data in the frequency-space domain using a complex numerical convolutional neural network (CVU-Net). Finally, they have proposed a self-supervised 5D seismic data reconstruction method based on neural implicit representation (NeRSI), which achieves efficient mapping from coordinates to seismic profiles by inputting coordinates into a multilayer perceptron (MLP) and introducing a convolutional neural network as a decoder. Finally, they have proposed a 5D seismic data reconstruction method based on continuous implicit neural representation, called implicit seismic representation. This method maps the coordinates of the seismic wavefield to amplitude values ​​using a multilayer perceptron, achieving 5D seismic data reconstruction. Researchers have proposed an unsupervised 5D seismic data reconstruction method based on implicit neural representations (INR). This method utilizes an MLP with a sinusoidal activation function to achieve a direct mapping from coordinates to amplitude, compatible with both regular and irregular meshes, and achieving efficient reconstruction and noise suppression. The researchers also proposed an unsupervised deep learning framework for joint denoising and reconstruction of 5D seismic data. This method combines the Transformer architecture with the POCS algorithm, achieving high-dimensional data reconstruction through block processing and attention mechanisms.

[0026] To achieve efficient and high-precision five-dimensional seismic data reconstruction, this invention proposes a five-dimensional seismic data reconstruction method based on Preconditioned Riemannian Gradient Descent (PRGD) for joint denoising and reconstruction of five-dimensional seismic data. This method is essentially a non-convex low-rank recovery method based on Riemannian manifold optimization. The PRGD algorithm effectively improves upon traditional Riemannian Gradient Descent (RGD) by introducing a data-driven weighted metric (preconditioner) within the Riemannian manifold optimization framework. Its core idea is to construct a concise and efficient preconditioner using the norms of the gradient matrix's rows and columns, thereby dynamically adjusting the search direction during iteration. While maintaining the same computational complexity as RGD, it significantly improves convergence speed and numerical stability. Theoretically, this method avoids the over-parameterization problem caused by traditional factorization methods and has a linear convergence guarantee under the Restricted Isometry Property (RIP) condition; its convergence rate is independent of the condition number of the low-rank matrix.

[0027] Given the inherent low-rank and multidimensional structure of seismic exploration data, extending the PRGD algorithm to denoising and reconstruction of five-dimensional seismic data not only maintains stable reconstruction performance under conditions of strong noise and undersampling, but also significantly reduces computational burden while improving the quality of results. This method avoids large-scale singular value decomposition through subspace projection and can skip the subspace approximation process, directly calculating subspace projection and singular value decomposition of small matrices, thus significantly improving computational speed when processing large-scale seismic data. This invention compares the proposed method with the classical DRR method, using five-dimensional synthetic data and real five-dimensional seismic data with 85.23% missing data to verify the performance of the proposed method. Experiments with synthetic and field data show that, compared to the classical Damped Rank Reduction (DRR) method, the PRGD method, under conditions of severe data loss (e.g., 90%) and strong noise, not only achieves a higher reconstruction signal-to-noise ratio and better wavefield fidelity, but also significantly reduces computation time. This research provides a novel and reliable solution for efficient and high-precision processing of high-dimensional seismic data.

[0028] The technical solutions provided by the various embodiments of this application are described in detail below with reference to the accompanying drawings.

[0029] Figure 1 This is a schematic diagram of a five-dimensional seismic data reconstruction method based on preconditional Riemann gradient descent, as described in this specification. The method specifically includes the following steps: S101: Obtain five-dimensional seismic data in frequency domain representation, and slice the five-dimensional seismic data at fixed frequencies to obtain four-dimensional seismic data.

[0030] In one embodiment, obtaining five-dimensional seismic data in the frequency domain includes: obtaining noisy and missing five-dimensional seismic data in the time domain; and performing a Fourier transform on the noisy and missing five-dimensional seismic data in the time domain to obtain five-dimensional seismic data in the frequency domain.

[0031] Specifically, set This represents noisy and incomplete five-dimensional seismic data in the time domain. The five-dimensional seismic data in the time domain is then transformed into five-dimensional seismic data in the frequency domain through a single-channel Fourier transform. .in, Indicates time, and These represent the spatial coordinates of the common center point (CMP) along the main survey line (Inline) and the connecting survey line (Crossline), respectively. and These represent the offset components along the main survey line and the connecting survey line, respectively, with spatial variables being... , Represents frequency, with a fixed frequency, and samples are taken only from the spatial domain. reconstruction, , This refers to four-dimensional seismic data, with the four dimensions being... .

[0032] S102, perform Hankel transformation on the four-dimensional seismic data to obtain a fourth-order block Hankel matrix.

[0033] In one embodiment, performing a Hankel transform on four-dimensional seismic data to obtain a fourth-order block Hankel matrix includes the following steps: S201 embeds four-dimensional seismic data into a first-order block Hankel matrix using all components of the first dimension of the four-dimensional seismic data.

[0034] definition First use All components of the first dimension embed the seismic data into a first-order block Hankel matrix:

[0035] (1) The dimension of a first-order block Hankel matrix is .

[0036] S202 uses all components of the second dimension in the four-dimensional seismic data to embed a first-order Hankel matrix into a second-order block Hankel matrix.

[0037] The 2nd order block Hankel matrix is: (2) The dimension of the 2nd order block Hankel matrix is .

[0038] S203 uses all components of the third dimension in four-dimensional seismic data to embed a second-order Hankel matrix into a third-order block Hankel matrix.

[0039] The 3rd order block Hankel matrix is: (3) The dimension of a 3rd order block Hankel matrix is .

[0040] S204 uses all components of the fourth dimension in the four-dimensional seismic data to embed a 3rd-order Hankel matrix into a 4th-order block Hankel matrix.

[0041] The 4th-order block Hankel matrix is: (4) The dimension of a 4th-order block Hankel matrix is .

[0042] S103, perform hard thresholding on the fourth-order block Hankel matrix to obtain a low-rank approximation matrix.

[0043] Define the hard thresholding operation of a matrix as follows: , Indicates from 1 to r integers i Summation. Here... r It is the preset matrix rank, indicating that only the largest rank is retained. r Each singular value and its corresponding singular vector. For the first i A singular value, It is a left singular vector. This is the transpose of the right singular vector. Using the hard threshold operator... As a shrinkage operator, the fourth-order block Hankel matrix Mapped to rank r and size n 1× n The Riemannian manifold formed by matrices of size 2 yields a low-rank approximation matrix. The calculation formula can be expressed as:

[0044] (5) S104. In any iteration of the preconditional Riemann gradient descent method, the Euclidean gradient at the low-rank approximation matrix at the current iteration point is calculated based on the fourth-order block Hankel matrix, the low-rank approximation matrix, and the sampling operator. The first preconditioner is determined based on the diagonal matrix formed by the outer products of the row vectors of the Euclidean gradient, and the second preconditioner is determined based on the diagonal matrix formed by the outer products of the column vectors of the Euclidean gradient. The Euclidean gradient represents the residual between the low-rank approximation matrix and the fourth-order block Hankel matrix. The sampling operator represents the linear operator that selects known observation data from the fourth-order block Hankel matrix.

[0045] In one embodiment, calculating the Euclidean gradient at the low-rank approximation matrix in the current iteration, based on the fourth-order block Hankel matrix, the low-rank approximation matrix, and the sampling operator, includes: calculating the product of the low-rank approximation matrix in the current iteration and the sampling operator, and determining the difference between this product and the fourth-order block Hankel matrix as the Euclidean gradient at the low-rank approximation matrix in the current iteration. The formula for calculating the Euclidean gradient is:

[0046] (6) in, Indicates the first The Euclidean gradient of the next iteration. Represents the sampling operator. For element-wise product, Indicates the first The low-rank approximation matrix of the next iteration. This represents a fourth-order block Hankel matrix.

[0047] In one embodiment, determining a first preconditioner based on the diagonal matrix formed by the outer product of the row vectors of the Euclidean gradient, and determining a second preconditioner based on the diagonal matrix formed by the outer product of the column vectors of the Euclidean gradient, includes: calculating the product of the Euclidean gradient and its transpose to obtain a first matrix; extracting the diagonal elements of the first matrix to form a first diagonal matrix; obtaining the first preconditioner based on the first diagonal matrix and an identity matrix weighted by a regularization parameter; calculating the product of the transpose of the Euclidean gradient and itself to obtain a second matrix; extracting the diagonal elements of the second matrix to form a second diagonal matrix; and obtaining the second preconditioner based on the second diagonal matrix and an identity matrix weighted by a regularization parameter. The calculation formulas for the first and second preconditioners are as follows:

[0048] (7) in, Indicates the first The first preconditioner of the next iteration. Indicates the first The second preconditioner of the next iteration Represents the regularization parameter. Indicates the first The Euclidean gradient of the next iteration. Represents a diagonal matrix. express × The identity matrix, express × The identity matrix, express the number of rows, express The number of columns.

[0049] S105, based on the first and second preconditioners, calculate the Euclidean gradient to a value of rank r and magnitude r. × The projection of the matrix onto the tangent space of the Riemannian manifold is obtained, and the low-rank approximation matrix is ​​updated according to the gradient after projection. The updated result is then mapped back to the Riemannian manifold to obtain the low-rank approximation matrix for the next iteration.

[0050] In one embodiment, based on the first and second preconditioners, the Euclidean gradient is calculated to a value of rank r and magnitude r. n 1× n The projection onto the tangent space of a Riemannian manifold composed of matrices of size 2 includes: calculating the weighted inner product of the Euclidean gradient using the first and second preconditioners; and calculating the rank r and size of the Euclidean gradient under the weighted inner product. n 1× n The tangent space projection of a Riemannian manifold composed of matrices of size 2.

[0051] The low-rank approximation matrix is ​​updated based on the projected gradient, including: the rank of the weighted inner product is r, and the size is..., based on a fixed constant step size and the Euclidean gradient. n 1× n The tangent space projection of the Riemannian manifold composed of matrices of size 2 is used to update the low-rank approximation matrix, yielding the updated result.

[0052] Specifically, the update formula for the low-rank approximation matrix is: (8) in, Indicates the first The update result of the low-rank approximation matrix in the next iteration. Indicates the first The low-rank approximation matrix of the next iteration. Indicates a fixed constant step size. Denotes the tangent space of the Riemannian manifold composed of low-rank matrices of rank r and size n1×n2 under the weighted inner product of the Euclidean gradient. Projection on This represents the tangent space projection operator under the weighted inner product.

[0053] Mapping the update result back to the Riemannian manifold to obtain the low-rank approximation matrix for the next iteration includes: performing a hard thresholding operation on the update result of the current iteration to obtain the low-rank approximation matrix for the next iteration. For example, the first... The update result of the low-rank approximation matrix in the nth iteration is mapped back to the Riemannian manifold to obtain the nth... The low-rank approximation matrix of the nth iteration. for:

[0054] (9) S106, update the low-rank approximation matrix for the next iteration until the square of the Frobenius norm of the difference matrix between the low-rank approximation matrices of two adjacent iterations is less than or equal to the iteration stopping error or the maximum number of iterations is reached.

[0055] when Or, upon reaching the maximum number of iterations, output parameters: Or, the low-rank approximation matrix obtained by iteration at the maximum number of iterations. This represents the iteration stopping error.

[0056] S107 transforms the low-rank approximation matrix of the last iteration into an inverse Hankel transformation into a fourth-order tensor, and determines the reconstructed five-dimensional seismic data based on the fourth-order tensor and the frequency components in the five-dimensional seismic data.

[0057] In one embodiment, the reconstructed five-dimensional seismic data is determined based on the frequency components in the fourth-order tensor and the five-dimensional seismic data, including: determining the reconstructed frequency domain seismic data based on the frequency components in the fourth-order tensor and the five-dimensional seismic data; and performing an inverse Fourier transform on the reconstructed frequency domain seismic data to recover the complete time domain seismic data, i.e., the reconstructed five-dimensional seismic data.

[0058] The above process is used to reconstruct all frequency components, resulting in reconstructed five-dimensional seismic data.

[0059] In one embodiment, the principle of the preconditional Riemann gradient descent method is given: Converting single-frequency 4D hypercube data into a fourth-order block Hankel matrix can be represented using the following operators: (10) in, This represents the Hankelization operator. For convenience, it is omitted here. parameters Assuming that the complete, noise-free low-rank seismic data obtained under ideal conditions is Z, then Z and the observed data... The mathematical relationship between them can be expressed as follows:

[0060] (11) in, For element-wise multiplication, matrix The size is , , , This is the sampling operator.

[0061] Based on the above construction, this invention transforms the seismic data reconstruction problem into a low-rank matrix completion problem. Under Riemannian metric, this problem can be formulated as solving a least-squares problem on a manifold, and a non-convex model is established as follows:

[0062] (12) in, It is a smooth Riemannian manifold formed by all matrices of size n1×n2 and rank r. This invention can be solved using the Riemann gradient descent algorithm on the manifold, and its general framework is as follows:

[0063] (13) in, for In manifold gradient on, For a fixed constant step size, the shrinkage operator To ensure that the iteration results do not deviate from the manifold, For the first The solution for the next iteration.

[0064] To solve the above problem, we first need to calculate the objective function. In the Next iteration point The Euclidean gradient at point E. According to the matrix differentiation rule, this Euclidean gradient can be expressed as:

[0065] (14) For the sake of brevity, this Euclidean gradient is denoted as... Physically, This represents the residual between the current recovered data and the observed data at the sampling location.

[0066] To further accelerate Riemann gradient descent, the aforementioned Euclidean gradient is utilized. The data-driven adaptive preconditioner is designed as shown in equation (7). This preconditioner uses the second moment information of the gradient to reweight the update weights of rows and columns, which can effectively improve the condition number under ill-conditioned sampling. For any matrix Define the following weighted inner product:

[0067] (15) This invention restricts this inner product to a subspace. superior, Representing a smooth manifold exist The tangent space at a point. Matrices in this tangent space have a specific structure, and their rank does not exceed 2r. For simplified notation, use... To represent all tangent space At this point, the objective function is in the tangent space. The Riemann gradient on can be expressed as the projection of the Euclidean gradient onto the weighted inner product:

[0068] (16) in Let this be the tangent space projection operator under the new weighted inner product. Assume the current iterative solution... The singular value decomposition at point is Therefore, we construct weighted orthogonalized basis vectors. and as follows:

[0069] (17) Using the aforementioned basis vectors, the projection operator For any matrix The formula for its function is: (18) Define the hard thresholding operation of a matrix as follows: , Indicates from 1 to r integers i Summation. Here... r It is the preset matrix rank, indicating that only the largest rank is retained. r Each singular value and its corresponding singular vector. For the first i A singular value, It is a left singular vector. This is the transpose of the right singular vector. Using the hard threshold operator... As a shrinkage operator, the updated result is mapped back to the Riemannian manifold. The preconditional Riemann gradient descent algorithm (PRGD) is obtained as follows:

[0070] (19) In one embodiment, the five-dimensional seismic data reconstruction method of the present invention includes a projection operator. contraction operator ,gradient and preconditions and The calculation was performed using Riemannian manifolds. Structural characteristics, operators It can be computed efficiently.

[0071] In the above calculation process, According to the tangent space The structural form, Its rank is at most 2r. Therefore, in calculating In this case, there is no need to perform singular value decomposition on large-scale matrices; instead, it can utilize... The low-rank structure enables efficient computation. For any It can be calculated directly. :

[0072] (20) For matrix and Perform QR decomposition to obtain the following results: and Obviously, , . It can be rewritten as follows:

[0073] (twenty one) in, It is a 2r×2r matrix. Because and They are all orthogonal matrices. Singular value decomposition can be derived from The singular value decomposition is obtained, and the computational complexity of this decomposition is only... Therefore, in calculation At that time, due to the iterative algorithm The row space and column space can be used to approximate Since the tangent space is defined, the subspace approximation process in traditional subspace projection methods can be skipped, and the subspace projection can be calculated directly. The SVD of small matrices avoids large-scale matrix operations, significantly improving the computational efficiency of the algorithm.

[0074] In one embodiment, the process of the preconditioning Riemann gradient descent method provided by the present invention includes: In one embodiment, to verify the effectiveness of the method of the present invention, the PRGD algorithm proposed in this invention is compared and analyzed with the DRR method. A noise-free five-dimensional seismic dataset with four spatial dimensions is constructed, with a tensor size of L×L×L×L, where L increases from 10 to 15. Each seismic data track contains 100 time sampling points with a sampling interval of 4ms. Based on this dataset, the experimental analysis will focus on the following four aspects: first, the change in computational cost as the data size increases; second, the evaluation of reconstruction performance and computational efficiency under different sampling rates; third, the sensitivity analysis of the algorithm to noise; and fourth, the analysis of reconstruction effect and computation time under different iteration numbers.

[0075] First, a comparative experiment was conducted on computation time under different five-dimensional data volume sizes. Seismic traces with 90% random loss were used in the above dataset, and reconstruction was performed using both the DRR algorithm and the PRGD algorithm proposed in this invention. The general parameter settings for the two algorithms are as follows: number of iterations... ,rank Iteration stopping error limit The reconstructed frequency bandwidth is The damping factor of the DRR algorithm is Regularization parameters of the PRGD algorithm The step size is selected to obtain the optimal step size for each reconstruction result. The comparison results of computation time are as follows: Figure 2 As shown, Figure 2 A comparison of the computation time of PRGD and DRR methods for 5D data of different sizes (100×L×L×L×L, L=10, 11, ..., 15) shows that as the scale parameter L increases from 10 to 15, the computation time of the DRR algorithm exhibits a significant non-linear growth trend, increasing from approximately 80 seconds to approximately 1500 seconds. In contrast, the computation time of the PRGD algorithm increases more gradually, only increasing from approximately 10 seconds to approximately 200 seconds. Experimental results indicate that, under the same computational conditions, the PRGD algorithm has a shorter computation time and higher computational efficiency.

[0076] Secondly, a comparative experiment was conducted to examine the reconstruction results and computation time of the two methods at different sampling rates. Five-dimensional seismic data with L=10 was selected from the seismic dataset to test the reconstruction quality of the PRGD method; subsequent experiments will also be based on this data. For quantitative analysis, the signal-to-noise ratio (SNR) was used as the evaluation index (Zhang et al., 2017), defined as:

[0077] (twenty two) In the formula, and These are the original data after vectorization and the reconstructed data after vectorization, respectively. The experiment removed 50%-90% of the seismic traces from the original data. Both methods used the aforementioned parameters unchanged, and the experimental results are as follows: Figure 3 As shown, Figure 3 This is a comparison chart of the reconstruction performance and reconstruction time of the PRGD method and the DRR method under different data missing rates. Figure 3 Figure (a) shows the change in signal-to-noise ratio (SNR) with missing data: as the proportion of missing data increases, the SNR of both algorithms decreases, indicating that data integrity has a direct impact on the reconstruction effect. However, under high missing data conditions, the PRGD algorithm can still maintain a high SNR, demonstrating excellent robustness. Figure 3 Figure (b) in the diagram compares the computation time of the two algorithms: the computation time of the DRR algorithm increases significantly with the missing data rate, while the computational efficiency of PRGD remains stable. In summary, the proposed algorithm has significant advantages in terms of computational performance and stability when processing data with high missing data rates.

[0078] Select Figure 3 The proposed algorithm was used to reconstruct and compare the seismic data with the DRR method, which was 90% missing in the experimental conditions. To visually demonstrate the reconstruction details of the five-dimensional data, the data volume was displayed in slices. Figure 4 Demonstrates a single common centroid (CMP) gather ( ) Comparison of 5D synthetic data results. This gather is obtained through a fixed spatial location ( ) and longitudinal offset ( The data was extracted and mainly reflects the changes in data with the lateral offset distance ( ). The characteristics of change. Among them, Figure 4 Figure (a) shows the actual data for that slice. Figure 4 Figure (b) in the figure shows the corresponding 90% missing observation data. Figure 4 Figures (c) and (d) show the reconstruction results of the DRR method and the PRGD method, respectively. Figure 4 Figures (e) and (f) in the figure show the reconstruction errors of the DRR method and the PRGD method. Figure 5 Showing common offset gathers ( ) Comparison of 5D synthetic data results. This gather is obtained by fixing all offsets ( ) and vertical spatial coordinates ( The data was extracted and the seismic waveform was displayed as it changed with the lateral spatial location. ) continuous change. Figure 5 The order of the subgraphs and Figure 4Consistent. Reconstruction results show that the signal-to-noise ratios (SNRs) of the observed data, the DRR method, and the proposed algorithm are 0.46 dB, 5.27 dB, and 28.97 dB, respectively. The DRR method performs poorly under extreme missing data conditions and has a large residual. Figure 4 (e) diagram and Figure 5 (e) diagram); in contrast, the PRGD method ( Figure 4 Figure (b) and Figure 5 Figure (b) not only effectively recovered the missing wavefield information, but also better maintained the continuity of the phase axis and amplitude characteristics, demonstrating its advantages in high-dimensional data reconstruction.

[0079] To evaluate the proposed algorithm's sensitivity to noise, Gaussian random noise was added to the noise-free five-dimensional simulated seismic data. Considering that excessively high noise variance can easily lead to physical distortion of data in actual seismic data processing, the effective study range for noise variance was set between 0.1 and 0.3. Simultaneously, to verify the algorithm's recovery performance under extremely sparse conditions, tests were conducted using two high missing rate scenarios: 80% and 90%. For noisy seismic data, the damping factor in the DRR method... PRGD algorithm step size Other basic parameter settings remained consistent with the aforementioned experiments. The reconstruction results of the two algorithms under different noise levels and missing rates are shown in Table 1. Experimental data show that the reconstruction quality of the PRGD algorithm is significantly better than that of the DRR method under all test conditions. In a low-noise environment with a noise variance of 0.1, the PRGD algorithm exhibits a significant advantage: at 80% and 90% missing rates, its signal-to-noise ratio (SNR) reaches 17.14 dB and 13.83 dB, respectively, far exceeding the 9.90 dB and 4.68 dB of the DRR method. Even under strong interference conditions with a noise variance increasing to 0.3, the PRGD algorithm can still maintain an SNR of 4.80 dB at a high missing rate of 90%, while the DRR method drops to 3.17 dB at this point. In summary, the PRGD algorithm not only possesses high reconstruction accuracy in low-noise environments but also demonstrates stronger robustness and stability under complex noise and high missing rate conditions.

[0080] Table 1. Comparison of Reconstructed Signal-to-Noise Ratio (dB) under Different Noise Variances and Missing Rates To visually verify the conclusions in Table 1, experimental data with a noise variance of 0.1 and a missing rate of 80% were selected for visualization. Figure 6 and Figure 7 The results of the 5D composite data of the common center point gather and the common offset gather are presented respectively. Noisy data was generated by adding Gaussian noise with a variance of 0.1 to the real data. Figure 6 Figure (a) and Figure 7(See Figure (a)). Subsequently, 80% of the seismic traces were randomly selected from the noisy data to obtain noisy observed seismic data (80% of the seismic traces are missing), as shown in Figure (a). Figure 6 Figure (b) in the middle and Figure 7 The results of the DRR method reconstruction are shown in Figure (b). Figure 6 Figure (c) in the middle and Figure 7 As shown in Figure (c), the phase axis is not clear. Figure 6 (d) diagram and Figure 7 Figure (d) shows the PRGD reconstruction result. Figure 6 Figure (e) in the middle and Figure 7 Figure (e) shows the DRR reconstruction error. Figure 6 Figure (f) in the middle and Figure 7 Figure (f) shows the PRGD reconstruction error. The signal-to-noise ratios (S / N) of the noisy data, the observed data, the data recovered by the DRR method, and the data recovered by the proposed method are -6.68 dB, -1.36 dB, 9.90 dB, and 17.14 dB, respectively. Compared with the DRR method, the PRGD algorithm proposed in this invention exhibits superior reconstruction performance under conditions of strong noise and high sparse sampling. Although the DRR method can recover the main phase axes, the phase axes are blurred and the structural features are unclear, resulting in an S / N of only 9.90 dB. In contrast, the PRGD algorithm effectively suppresses the influence of noise, significantly improves the continuity and amplitude fidelity of the seismic phase axes, and achieves an S / N of 17.14 dB, which is about 5.65 dB higher than that of DRR. The results show that the PRGD algorithm has stronger robustness and reconstruction accuracy, and can achieve high-fidelity recovery of complex seismic data.

[0081] Finally, the convergence characteristics and computational efficiency of the proposed algorithm were further investigated. While keeping the aforementioned experimental conditions unchanged, the number of iterations was gradually increased from 2 to 20, and the changes in signal-to-noise ratio and computation time of the two algorithms were recorded. The results are as follows: Figure 8 As shown, Figure 8 A comparison of the reconstruction performance and time of PRGD and DRR methods under different iteration numbers. Figure 8 Figure (a) shows the signal-to-noise ratio (SNR) as a function of iterations. The PRGD algorithm (red line) exhibits extremely fast convergence, with a sharp increase in SNR during the first 5 iterations, reaching a peak around the 10th iteration (approximately 17.14 dB) and then stabilizing. In contrast, the DRR algorithm (blue line) shows a slow, linear increase in SNR, failing to reach the reconstruction level of PRGD in its 10th iteration by the 20th iteration, indicating its lower convergence efficiency. Figure 8Figure (b) illustrates the difference in time cost between the two algorithms. The computation time of the DRR algorithm increases significantly linearly with the number of iterations, exceeding 220 seconds for 20 iterations. In contrast, the PRGD algorithm's time curve is extremely flat, showing only a slight increase with the number of iterations, with the overall time remaining within 20 seconds. In summary, the PRGD algorithm achieves efficient convergence in approximately 10 iterations, and its computational cost per iteration is extremely low. This characteristic gives it significant computational advantages and practical value in processing large-scale multidimensional seismic data.

[0082] In one embodiment, the reconstruction performance of the PRGD method is verified by testing on five-dimensional real-world data. This data space has a dimension of 20×20×20×20, representing... The sampling interval was 2ms, and the number of sampling points was 540. 85.23% of the seismic traces were missing. To effectively recover this seismic data, the general parameter settings for both algorithms are as follows: number of iterations. ,rank Iteration stopping error limit The reconstructed frequency bandwidth is The damping factor of the DRR algorithm is Regularization parameters of the PRGD algorithm Step length . Figure 9 For a common center point set 3D graphics comparison Figure 10 A comparison plot of Fk spectra from field data; Figure 9 Figure (a) and Figure 10 Figure (a) shows the data before reconstruction and its spectrogram. This data was obtained by fixing... and The signal and spectrum reconstructed by the DRR algorithm are obtained. Figure 9 Figure (b) in the middle and Figure 10 Figure (b) in the middle. Figure 9 Figure (c) in the middle and Figure 10 Figure (c) shows the reconstruction result and spectrum obtained by the method proposed in this invention. Furthermore, this invention provides a two-dimensional slice representation of this common offset gather, as shown below. Figure 11 As shown, Figure 11 For fixed Comparison diagram of two-dimensional unfolded lookup of common center point set. Figure 11 (a) in the image represents the earthquake data (observation data) before reconstruction. Figure 11 Figures (b) and (c) show the reconstruction results of the DRR method and the method proposed in this invention, respectively. To present the details more clearly, Figure 12 The reconstruction results under different slices were further demonstrated. Figure 12 For fixed The local comparison images of the common offset gathers, and the reconstruction results of DRR and the algorithm proposed in this invention are shown below. Figure 12 As shown in Figures (a) and (d) in the document. Fixed. , get The changing DRR and PRGD reconstruction results, such as Figure 12 Figures (b) and (e) in Figure 12 are shown. Figures (c) and (f) in Figure 12 are magnified views of Figures (b) and (e) in Figure 11, respectively, within the black boxes. From the reconstruction results above, it can be seen that the proposed method effectively recovers the seismic phase axis, and the reconstructed signal exhibits good continuity and clear structural features.

[0083] In addition, step size This is a key hyperparameter affecting the convergence and reconstruction accuracy of the preconditioned Riemann gradient descent (PRGD) algorithm. To determine the optimal step size range and evaluate the algorithm's sensitivity to this parameter, this invention, based on the aforementioned synthetic data experiments, tested the signal-to-noise ratio (S / N) as a function of step size under different combinations of missing rates (70%-90%) and noise levels (0.1-0.3). The changing pattern, the results are as follows: Figure 13 As shown, Figure 13 This diagram illustrates the selection of step size under various missing rates and signal-to-noise ratio conditions. The trajectory in the diagram shows a significant correlation between reconstruction performance and step size. Within a relatively small step size range, the signal-to-noise ratio increases with... The increase in λ is approximately linear, indicating that an excessively small step size limits the efficiency of gradient descent, leading to insufficient algorithm convergence. As λ increases... Further increase ( The signal-to-noise ratio curves under all conditions gradually stabilized. Notably, despite changes in the missing rate and noise intensity within this range, the algorithm maintained relatively consistent optimal performance without exhibiting drastic oscillations or divergence.

[0084] Experimental results show that when the step size is within In this process, the PRGD algorithm can balance convergence speed and reconstruction accuracy. Furthermore, to ensure the numerical stability of the algorithm under finite sample and noisy environments, a regularization strategy is introduced to prevent the condition number of the matrix from becoming too large. In all experiments of this invention, the regularization parameter is fixed at [value missing]. ,Cooperate The step size setting ensures the robustness of the algorithm under different data conditions.

[0085] This method proposes a simultaneous denoising and reconstruction algorithm for five-dimensional seismic data by introducing a preconditional Riemann gradient descent algorithm. The PRGD algorithm achieves efficient recovery of low-rank matrices by combining manifold constraint optimization with data-driven preprocessing. The proposed method avoids large-scale matrix decomposition by using a fast low-rank approximation method based on tangent space projection. Experiments show that the PRGD method is significantly more computationally efficient than the traditional damped rank reduction (DRR) method, and its computation time increases more gradually with the data size. In terms of reconstruction quality, whether in synthetic or real data, PRGD can more effectively recover the continuity and structural characteristics of the seismic wavefield under conditions of high missing rates and strong noise, demonstrating stronger robustness.

[0086] The methods provided in this specification can be executed by a server, which can be a server set up on a business platform, or a device such as a desktop computer or laptop computer capable of executing the solutions in this specification. For ease of explanation, the following description will only use a server as the execution subject.

[0087] When applying the five-dimensional seismic data reconstruction method based on preconditional Riemann gradient descent provided in this manual, it is not necessary to consider... Figure 1 The steps shown are executed in sequence. The specific execution order of each step can be determined as needed, and this manual does not impose any restrictions on it.

[0088] The above describes a five-dimensional seismic data reconstruction method based on preconditional Riemann gradient descent, provided by one or more embodiments of this specification. Based on the same idea, this specification also provides a corresponding five-dimensional seismic data reconstruction device based on preconditional Riemann gradient descent, which includes: The acquisition module is used to acquire five-dimensional seismic data in frequency domain representation and to slice the five-dimensional seismic data at fixed frequencies to obtain four-dimensional seismic data. The transformation module is used to perform Hankel transformation on four-dimensional seismic data to obtain a fourth-order block Hankel matrix. The computation module is used to perform hard thresholding on the fourth-order block Hankel matrix to obtain a low-rank approximate matrix. The iterative module is used to calculate the Euclidean gradient matrix at the low-rank approximation matrix in any iteration of the preconditioned Riemann gradient descent method, based on the fourth-order block Hankel matrix, the low-rank approximation matrix, and the sampling operator. The Euclidean gradient matrix represents the residual between the low-rank approximation matrix and the fourth-order block Hankel matrix. The sampling operator represents the linear operator that selects known observation data from the fourth-order block Hankel matrix. The first preconditioner is determined based on the diagonal matrix formed by the outer products of the row vectors of the Euclidean gradient matrix, and the second preconditioner is determined based on the diagonal matrix formed by the outer products of the column vectors of the Euclidean gradient matrix. Based on the first and second preconditioners, the Euclidean gradient is calculated to a value of rank r and size r. × The projection of the matrix into the tangent space of the Riemann manifold is obtained, and the low-rank approximation matrix is ​​updated according to the gradient after projection. The updated result is then mapped back to the Riemann manifold to obtain the low-rank approximation matrix for the next iteration. The low-rank approximation matrix for the next iteration is updated again until the square of the Frobenius norm of the difference matrix between the low-rank approximation matrices of two adjacent iterations is less than or equal to the iteration stopping error or the maximum number of iterations is reached. The determination module is used to transform the low-rank approximation matrix of the last iteration into an inverse Hankel transformation into a fourth-order tensor, and to determine the reconstructed five-dimensional seismic data based on the fourth-order tensor and the frequency components in the five-dimensional seismic data.

[0089] Specific limitations regarding the five-dimensional seismic data reconstruction device based on preconditional Riemann gradient descent can be found in the limitations of the five-dimensional seismic data reconstruction method based on preconditional Riemann gradient descent described above, and will not be repeated here. Each module in the aforementioned five-dimensional seismic data reconstruction device based on preconditional Riemann gradient descent can be implemented entirely or partially through software, hardware, or a combination thereof. These modules can be embedded in or independent of the processor in a computer device, or stored in the memory of a computer device as software, so that the processor can call and execute the corresponding operations of each module.

[0090] This specification also provides a computer-readable storage medium storing a computer program that can be used to execute the above-described... Figure 1 The provided method is a five-dimensional seismic data reconstruction method based on preconditional Riemann gradient descent.

[0091] This instruction manual also provides Figure 14 The schematic diagram of the computer device shown is as follows: Figure 14 At the hardware level, the computer device includes a processor, internal bus, network interface, memory, and non-volatile memory, and may also include other hardware required for business operations. The processor reads the corresponding computer program from the non-volatile memory into memory and then runs it to achieve the above-mentioned functions. Figure 1 The provided method is a five-dimensional seismic data reconstruction method based on preconditional Riemann gradient descent.

[0092] Those skilled in the art will understand that all or part of the processes in the methods of the above embodiments can be implemented by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer-readable storage medium, and when executed, it can include the processes of the embodiments of the methods described above. Any references to memory, storage, databases, or other media used in the embodiments provided in this application can include at least one of non-volatile and volatile memory. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, or optical storage, etc. Volatile memory can include random access memory (RAM) or external cache memory. By way of illustration and not limitation, RAM can be in various forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM), etc.

[0093] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

Claims

1. A five-dimensional seismic data reconstruction method based on preconditional Riemann gradient descent, characterized in that, include: Five-dimensional seismic data in frequency domain representation is obtained, and fixed-frequency slices are made from the five-dimensional seismic data to obtain four-dimensional seismic data. The four-dimensional seismic data is subjected to Hankel transformation to obtain a fourth-order block Hankel matrix; Hard thresholding is performed on the fourth-order block Hankel matrix to obtain a low-rank approximate matrix; In any iteration of the preconditional Riemann gradient descent method, the Euclidean gradient at the low-rank approximation matrix in the current iteration is calculated based on the fourth-order block Hankel matrix, the low-rank approximation matrix, and the sampling operator. The Euclidean gradient represents the residual between the low-rank approximation matrix and the fourth-order block Hankel matrix. The sampling operator represents the linear operator that selects known observation data from the fourth-order block Hankel matrix. The first preconditioner is determined by the diagonal matrix formed by the outer product of the row vectors of the Euclidean gradient, and the second preconditioner is determined by the diagonal matrix formed by the outer product of the column vectors of the Euclidean gradient. Based on the first and second preconditioners, the Euclidean gradient is calculated to a value of rank r and magnitude r. × The projection of the matrix into the tangent space of the Riemann manifold is obtained, and the low-rank approximation matrix is ​​updated according to the gradient after projection. The updated result is then mapped back to the Riemann manifold to obtain the low-rank approximation matrix for the next iteration. The low-rank approximation matrix is ​​updated again for the next iteration until the square of the Frobenius norm of the difference matrix between the low-rank approximation matrices of two adjacent iterations is less than or equal to the iteration stopping error or the maximum number of iterations is reached. The low-rank approximation matrix from the last iteration is transformed into a fourth-order tensor using the inverse Hankel transformation. Based on the fourth-order tensor and the frequency components in the five-dimensional seismic data, the reconstructed five-dimensional seismic data is determined.

2. The method according to claim 1, characterized in that, Based on the fourth-order block Hankel matrix, the low-rank approximation matrix, and the sampling operator, calculate the Euclidean gradient at the low-rank approximation matrix in the current iteration, including: Calculate the product of the low-rank approximation matrix of the current iteration and the sampling operator, and determine the difference between this product and the fourth-order block Hankel matrix as the Euclidean gradient at the low-rank approximation matrix of the current iteration.

3. The method according to claim 1, characterized in that, The first preconditioner is determined based on the diagonal matrix formed by the outer product of the row vectors of the Euclidean gradient, and the second preconditioner is determined based on the diagonal matrix formed by the outer product of the column vectors of the Euclidean gradient, including: Calculate the product of the Euclidean gradient and its transpose to obtain the first matrix, and extract the diagonal elements of the first matrix to form the first diagonal matrix. The first preconditioner is obtained based on the first diagonal matrix and the identity matrix weighted by the regularization parameter; Calculate the product of the transpose of the Euclidean gradient and itself to obtain the second matrix, and extract the diagonal elements of the second matrix to form a second diagonal matrix; The second preconditioner is obtained from the second diagonal matrix and the identity matrix weighted by the regularization parameter.

4. The method according to claim 3, characterized in that, Based on the first and second preconditioners, the Euclidean gradient is calculated to a value of rank r and magnitude r. n 1× n The projections onto the tangent space of a Riemannian manifold composed of matrices of size 2 include: The weighted inner product of the Euclidean gradient is calculated using the first and second preconditioners. Calculate the rank of the Euclidean gradient under the weighted inner product as r, with a magnitude of r. n 1× n The tangent space projection of a Riemannian manifold composed of matrices of size 2.

5. The method according to claim 4, characterized in that, The low-rank approximation matrix is ​​updated based on the projected gradient, including: Given a fixed constant step size and the rank of the Euclidean gradient under the weighted inner product as r, and a magnitude of... n 1× n The tangent space projection of the Riemannian manifold composed of matrices of size 2 is used to update the low-rank approximation matrix, yielding the updated result; the update formula for the low-rank approximation matrix is: in, Indicates the first The update result of the low-rank approximation matrix in the next iteration. Indicates the first The low-rank approximation matrix of the next iteration. Indicates a fixed constant step size. Denotes the tangent space of the Riemannian manifold composed of low-rank matrices of rank r and size n1×n2 under the weighted inner product of the Euclidean gradient. Projection on This represents the tangent space projection operator under the weighted inner product.

6. The method according to claim 5, characterized in that, Mapping the updated result back to the Riemannian manifold yields the low-rank approximation matrix for the next iteration, including: Perform a hard threshold operation on the update result of the current iteration to obtain the low-rank approximation matrix for the next iteration.

7. The method according to claim 1, characterized in that, Acquire five-dimensional seismic data in frequency domain representation, including: Acquire noisy and missing five-dimensional seismic data in the time domain; The five-dimensional seismic data in the time domain, which contains noise and has missing data, is subjected to Fourier transform to obtain the five-dimensional seismic data in the frequency domain.

8. The method according to claim 1, characterized in that, Performing a Hankel transform on the four-dimensional seismic data yields a fourth-order block Hankel matrix, including: Four-dimensional seismic data is embedded into a first-order block Hankel matrix using all components of the first dimension of the four-dimensional seismic data. The first-order Hankel matrix is ​​embedded into the second-order block Hankel matrix using all components of the second dimension in the four-dimensional seismic data. The second-order Hankel matrix is ​​embedded into a third-order block Hankel matrix using all components of the third dimension in the four-dimensional seismic data. The third-order Hankel matrix is ​​embedded into a fourth-order block Hankel matrix using all components of the fourth dimension in the four-dimensional seismic data.

9. The method according to claim 1, characterized in that, Based on the fourth-order tensor and the frequency components in the five-dimensional seismic data, the reconstructed five-dimensional seismic data is determined, including: Based on the frequency components in the fourth-order tensor and five-dimensional seismic data, the reconstructed frequency domain seismic data is determined; The reconstructed frequency domain seismic data is subjected to inverse Fourier transform to recover the complete time domain seismic data, i.e., the reconstructed five-dimensional seismic data.

10. A five-dimensional seismic data reconstruction device based on preconditional Riemann gradient descent, characterized in that, include: The acquisition module is used to acquire five-dimensional seismic data in frequency domain representation and to slice the five-dimensional seismic data at fixed frequencies to obtain four-dimensional seismic data. The transformation module is used to perform Hankel transformation on four-dimensional seismic data to obtain a fourth-order block Hankel matrix. The computation module is used to perform hard thresholding on the fourth-order block Hankel matrix to obtain a low-rank approximate matrix. The iterative module is used to calculate the Euclidean gradient matrix at the low-rank approximation matrix in any iteration of the preconditioned Riemann gradient descent method, based on the fourth-order block Hankel matrix, the low-rank approximation matrix, and the sampling operator. The Euclidean gradient matrix represents the residual between the low-rank approximation matrix and the fourth-order block Hankel matrix. The sampling operator represents the linear operator that selects known observation data from the fourth-order block Hankel matrix. The first preconditioner is determined based on the diagonal matrix formed by the outer products of the row vectors of the Euclidean gradient matrix, and the second preconditioner is determined based on the diagonal matrix formed by the outer products of the column vectors of the Euclidean gradient matrix. Based on the first and second preconditioners, the Euclidean gradient is calculated to a value of rank r and size r. × The projection of the matrix into the tangent space of the Riemann manifold is obtained, and the low-rank approximation matrix is ​​updated according to the gradient after projection. The updated result is then mapped back to the Riemann manifold to obtain the low-rank approximation matrix for the next iteration. The low-rank approximation matrix for the next iteration is updated again until the square of the Frobenius norm of the difference matrix between the low-rank approximation matrices of two adjacent iterations is less than or equal to the iteration stopping error or the maximum number of iterations is reached. The determination module is used to transform the low-rank approximation matrix of the last iteration into an inverse Hankel transformation into a fourth-order tensor, and to determine the reconstructed five-dimensional seismic data based on the fourth-order tensor and the frequency components in the five-dimensional seismic data.

Images