A method for reconstructing high-dimensional seismic data based on irregular grids based on fast low-rank prediction

Abstract

The present invention discloses a method for reconstructing high-dimensional seismic data of irregular grids based on fast low-rank prediction, and belongs to the technical field of seismic data processing. The problem of slow matrix decomposition calculation efficiency of the conventional seismic data reconstruction method based on low-rank prediction in the prior art is solved; the present invention obtains irregular grid seismic data by interpolating regular network seismic data, combines the mapping gradient descent algorithm for iterative solution, obtains updated seismic data, performs one-dimensional Fourier transform on it, then performs forward fast Hank matrix-vector product, performs low-rank prediction processing, and then performs transposed fast Hank matrix-vector product, uses random QR decomposition instead of singular value decomposition, uses convolution to accelerate anti-diagonal averaging, obtains the result after low-rank prediction, performs inverse Fourier transform on the conjugate matrix of frequency domain seismic data, and iteratively solves until convergence conditions are met to obtain reconstructed seismic data. The present invention effectively improves the efficiency and accuracy of seismic data reconstruction.

Description

Technical Field

[0001] The present invention relates to a method for reconstructing high-dimensional seismic data of irregular grids, and in particular to a method for reconstructing high-dimensional seismic data of irregular grids based on fast low-rank prediction, and belongs to the technical field of seismic data processing. Background Art

[0002] Reconstruction of seismic data from irregular grids to regular grids is a long-standing problem in seismic data processing. Due to environmental factors in the acquisition process or irregular acquisition design based on compressed sensing theory, the sources or detectors are usually irregularly distributed. The goal of seismic reconstruction is to regularize the acquired irregular grid data onto a regular grid and provide a higher density of coverage times for subsequent seismic data processing, thereby improving imaging accuracy.

[0003] In essence, the irregularities generated by field seismic data acquisition can be summarized into two categories. One is objective irregularities, which are caused by natural obstacles in the acquisition process and are inevitable. For example, there are obstacles such as buildings in the acquisition area, which cause large gaps in the observed data, resulting in discontinuity or even distortion of the seismic image, affecting subsequent imaging. In order to solve this problem, many scholars have proposed many seismic data reconstruction methods, such as minimum weighted norm interpolation (MWNI), anti-leakage Fourier transform (ALFT), prediction error filter, sparse Fourier inversion and improved minimum weighted norm interpolation (MWNI). The other is subjective irregularities. Regularity and subjective irregularity are irregular sampling observation systems predefined based on compressed sensing theory. Through artificial subjective design, the irregularity of pseudo-random distribution is used to reduce the number of field acquisition points, thereby shortening the field acquisition time. Through subsequent data reconstruction, seismic data with the same accuracy as conventional acquisition is obtained. By utilizing the spatial coherence of seismic signals in the transformed domain, seismic records compressed below the Nyquist frequency can be perfectly restored through sparse inversion. Related issues of subjective irregularity can usually be processed by the following algorithms: interpolation compressed sensing method, non-uniform optimal sampling method (NUOS), extended convex set projection method (EPOCS), etc. for data reconstruction.

[0004] In recent years, the method of using low-rank prediction to reconstruct irregular grid data has attracted more and more attention. The data reconstruction model based on low-rank prediction method has significant advantages in the accuracy of seismic data reconstruction from irregular grids to regular grids. However, in the existing seismic data reconstruction methods based on low-rank prediction, there are disadvantages such as slow matrix decomposition calculation efficiency and high memory load requirements. To solve the above problems, a method for reconstructing irregular grid high-dimensional seismic data based on fast low-rank prediction is needed. Summary of the invention

[0005] A brief overview of the present invention is provided below in order to provide a basic understanding of certain aspects of the present invention. It should be understood that this overview is not an exhaustive overview of the present invention. It is not intended to identify key or important parts of the present invention, nor is it intended to limit the scope of the present invention. Its purpose is merely to present certain concepts in a simplified form as a prelude to a more detailed description discussed later.

[0006] In view of this, in order to solve the problem of slow matrix decomposition calculation efficiency of traditional seismic data reconstruction methods based on low-rank prediction in the prior art, the present invention provides a non-regular grid high-dimensional seismic data reconstruction method based on fast low-rank prediction.

[0007] The technical solution is as follows: A method for reconstructing irregular grid high-dimensional seismic data based on fast low-rank prediction, comprising the following steps:

[0008] S1. Regularizing and preprocessing the irregular grid seismic data by using an interpolation operator to obtain regular grid seismic data;

[0009] S2. According to the regular grid seismic data, the mapping gradient descent algorithm is combined to perform iterative solution, and the exponential decay strategy is adopted for the parameter iteration step to ensure the convergence of the mapping gradient descent algorithm to obtain the updated seismic data;

[0010] S3. Perform one-dimensional Fourier transform on the updated seismic data to obtain frequency domain seismic data and its dimension size;

[0011] S4. Perform a forward fast Hank matrix-vector product according to the frequency domain seismic data and its dimension to obtain a matrix Y;

[0012] Specific:

[0013] S41. Obtain a randomly generated vector Ω according to the dimension of the frequency domain seismic data and a predetermined rank p;

[0014] S42. Perform one-dimensional Fourier transform on the input secondary Hank matrices c and v to obtain updated secondary Hank matrices C and V;

[0015] S43. Perform inverse Fourier transform on the updated secondary Hank matrices C and V to obtain a matrix r;

[0016] S44. truncate and reshape the matrix r to obtain the matrix Y;

[0017] S5. Perform low-rank prediction processing to obtain a low-rank approximation of the matrix Y

[0018] Specific:

[0019] S51. Perform random projection on the matrix Y to obtain the matrix M;

[0020] S52. Perform QR decomposition on the matrix M to obtain a standard orthogonal matrix Q and an upper triangular matrix R;

[0021] S53. According to the standard orthogonal matrix Q and the conjugate transposed matrix Q of the standard orthogonal matrix Q H , and obtain a low-rank approximation of the matrix Y

[0022] S6. Perform transposed fast Hank matrix-vector product on the matrix q obtained by reshaping the standard orthogonal matrix Q and the frequency domain seismic data to obtain a matrix z;

[0023] S7. According to the matrix q and the matrix z, convolution is used to accelerate the anti-diagonal average to obtain the result after low-rank prediction;

[0024] S8. Based on the results of low-rank prediction, the conjugate matrix of frequency domain seismic data is Perform inverse Fourier transform to obtain preliminary reconstructed seismic data, and define steps S3 to S8 as a mapping operation P;

[0025] S9. Repeat steps S1 to S8 until the convergence condition is reached to obtain the reconstructed seismic data.

[0026] Furthermore, in S1, the interpolation formula is expressed as:

[0027] u(t,x,y)=wd(t,x,y)

[0028] Among them, u(t,x,y) is the observed seismic data of the initial irregular grid, t is the sampling time of the seismic data of the initial irregular grid, (x,y) is the coordinate of the seismic data of the initial irregular grid, is the interpolation operator, and d(t,x,y) is the regular network seismic data.

[0029] Furthermore, in S2, the updated seismic data m(t,x,y) is expressed as:

[0030]

[0031] Among them, d k (t,x,y) is the regular network seismic data at the kth moment, is the adjoint operator of the interpolation operator w, λ is the iteration step size, and an exponential decay strategy is used. The initial step size is given by Decision, e max for The maximum eigenvalue of .

[0032] Furthermore, in S3, the frequency domain seismic data D(ω,x,y) is expressed as:

[0033]

[0034] in, is the one-dimensional Fourier transform;

[0035] The dimension size (N) of frequency domain seismic data D(ω,x,y) t ,N1,N2) is expressed as:

[0036] (N t ,N1,N2)=size(D(ω,x,y))

[0037] Among them, size(·) is the function to solve the matrix size.

[0038] Further, in S41, the vector Ω is expressed as:

[0039] Ω=rand(M1,M2,p)

[0040] M1=[N1 / 2]+1

[0041] M2=[N2 / 2]+1

[0042] Among them, M1 and M2 are parameters used to specify the size of vector Ω;

[0043] In S42, the input secondary Hank matrices are c and v, and the dimension size of the matrix c is obtained;

[0044] The dimension size of matrix c (L x ,L y ) is expressed as:

[0045] (L x ,L y )=size(c)

[0046] The dimensions of the secondary Hank matrices c and v are expressed as:

[0047]

[0048]

[0049] Where N is the row size of the secondary Hank matrix, M is the column size of the secondary Hank matrix, and when flag = forward, the size of the secondary Hank matrix is ​​expressed as N c 、M c When flag = adjoint, the size of the secondary Hank matrix is ​​expressed as N r 、M r ;

[0050] The updated secondary Hank matrices C and V are expressed as:

[0051]

[0052]

[0053] Among them, padding(·) is the function of filling zero elements, reverse(·) is the inversion function, To reverse the elements of each dimension of the matrix v in sequence and then fill them with zeros, its size is [N-1,M-1];

[0054] In S43, the matrix r is expressed as:

[0055]

[0056] in, is element-wise multiplication, is the inverse Fourier transform, conj(·) is the conjugate, δ(C) is the indicator function, and its expression is:

[0057]

[0058] In S44, the matrix Y is expressed as:

[0059] Y=reshape(r[N:L x ,M:L y ],[],1)

[0060] Among them, reshape is the reshaping function, and [] indicates the adaptive size of the matrix.

[0061] Further, in S51, according to the size N Y ×N Y The matrix Y of size N consists of p independent vectors Y ×p random matrix T to obtain the matrix M, where p<<N Y ;

[0062] The matrix M is represented as:

[0063] M=YT

[0064] In S52, the standard orthogonal matrix Q and the upper triangular matrix R are expressed as:

[0065] [Q,R]=qr(M)

[0066] Where qr(·) is the QR decomposition;

[0067] In S53, the low-rank approximation of the matrix Y It is expressed as:

[0068]

[0069] Furthermore, in S6, the standard orthogonal matrix Q is reshaped to obtain a matrix q;

[0070] The matrix q is expressed as:

[0071] q=reshape(Q,L1,L2)

[0072] L1=N1+1-M1

[0073] L2=N2+1-M2

[0074] Among them, L1 and L2 are parameters used to reshape the size of the standard orthogonal array Q;

[0075] Perform transposed fast Hank matrix-vector product on matrix q and frequency domain seismic data D(ω,x,y) to obtain matrix z;

[0076] The matrix z is represented as:

[0077]

[0078] Furthermore, in S7, the result after low-rank prediction is expressed as:

[0079]

[0080] in, is the result of the k+1th low-rank prediction, is the result of the k-th low-rank prediction.

[0081] Furthermore, in S8, the conjugate matrix of the frequency domain seismic data D(ω,x,y) The inverse Fourier transform of is expressed as:

[0082]

[0083] Among them, d k+1 is the initial reconstructed seismic data, d k+1 (t,x,y)=P(m(t,x,y)).

[0084] Furthermore, in S9, the convergence condition is that the relative error res is less than a fixed value, that is, res<1e -6 ;

[0085] The relative error res is expressed as:

[0086]

[0087] Among them, ||·|| F is the Frobenius norm.

[0088] The beneficial effects of the present invention are as follows: the present invention applies the I-FMSSA algorithm to seismic data reconstruction, effectively utilizes the information and structural characteristics of the original data, and thus improves the reconstruction efficiency and accuracy of the seismic data; the present invention uses an interpolation operator instead of a sampling operator to map the observed irregular grid to the desired regular grid, avoiding the introduction of phase and amplitude reconstruction errors; the present invention uses fast Fourier transform to calculate the Hank matrix-vector product, uses random QR decomposition instead of singular value decomposition, and uses convolution to accelerate anti-diagonal averaging processing, which improves the calculation efficiency and reduces the memory load compared to traditional methods. BRIEF DESCRIPTION OF THE DRAWINGS

[0089] The drawings described herein are used to provide a further understanding of the present invention and constitute a part of the present invention. The exemplary embodiments of the present invention and their descriptions are used to explain the present invention and do not constitute an improper limitation of the present invention. In the drawings:

[0090] Figure 1 A flowchart of an algorithm for reconstructing irregular grid high-dimensional seismic data based on fast low-rank prediction;

[0091] Figure 2 The simulated seismic data diagram for test 1, where (a) is simulated seismic data with irregular distribution, and (b) is simulated seismic data with regular grid;

[0092] Figure 3 Figure 1 shows the reconstruction of noise-free synthetic data for test 1, where (a) is the clean regular grid seismic data, (b) is the initial irregular grid seismic data, (c) is the irregular grid seismic data after the loss, (d) and (e) are the irregular reconstruction results of noise-free synthetic data obtained by the FISTA algorithm, and (f) and (g) are the irregular reconstruction results of noise-free synthetic data obtained by the I-FMSSA algorithm;

[0093] Figure 4 Figure 2 is the reconstruction diagram of noisy data for test 2, where (a) is the clean regular grid data, (b) is the initial noisy irregular grid data, (c) is the noisy irregular grid data after deletion, (d) is the irregular reconstruction result of noisy synthetic data obtained by the FISTA algorithm, (e) is the residual comparison between the reconstruction result obtained by the FISTA algorithm and the clean regular grid data, (f) is the irregular reconstruction result of noisy synthetic data obtained by the I-FMSSA algorithm, and (g) is the residual comparison between the reconstruction result obtained by the I-FMSSA algorithm and the clean regular grid data;

[0094] Figure 5 Figure 3 shows the complex synthetic data in test 3, where (a) is the observed irregular grid coordinates of the shot points and (b) is the expected regular grid output;

[0095] Figure 6 Reconstructed slice diagram of complex synthetic data in test three, where (a) is the original regular grid data, (b) is the observed irregular grid data, (c) is the irregular grid data after missing, (d) is the reconstruction result of FISTA algorithm, (e) is the difference between the original regular grid data and the reconstruction result of FISTA algorithm, (f) is the reconstruction result of I-FMSSA algorithm, and (g) is the difference between (a) and (f);

[0096] Figure 7 The shot point coordinate diagram for the 3D seismic data acquisition in test 4, where (a) is the observed shot point coordinates of the irregular grid and (b) is the expected output regular grid;

[0097] Figure 8 Reconstructed slices of 3D real data in test 4, where (a) is the observed irregular grid data, (b) is the irregular grid data after loss, (c) is the reconstruction result of FISTA algorithm, and (d) is the reconstruction result of I-FMSSA algorithm;

[0098] Fig. 9 This is the coordinate diagram of the shot points for the five-dimensional seismic data acquisition in test 5, where (a) is the coordinate distribution of the shot points and receiver points for the actual five-dimensional seismic data acquisition, and (b) is the expected regular grid shot points and receiver point distribution;

[0099] Fig.10 Slices of the shot point Y direction data volume for the fixed shot point and the detection point X direction in test 5, where (a) is the incomplete data containing considerable noise, (b) is the reconstruction result of the FISTA algorithm, and (c) is the reconstruction result of the I-FMSSA algorithm. DETAILED DESCRIPTION

[0100] In order to make the technical solutions and advantages of the embodiments of the present invention more clearly understood, the exemplary embodiments of the present invention are further described in detail below in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than an exhaustive list of all the embodiments. It should be noted that the embodiments of the present invention and the features in the embodiments can be combined with each other without conflict.

[0101] refer to Figure 1-10 The present embodiment is described in detail, and the irregular grid high-dimensional seismic data reconstruction method based on fast low-rank prediction specifically includes the following steps:

[0102] S1. Regularizing and preprocessing irregular network seismic data by using interpolation operators to obtain regular network seismic data;

[0103] S2. According to the regular grid seismic data, the mapping gradient descent algorithm is combined to perform iterative solution, and the exponential decay strategy is adopted for the parameter iteration step to ensure the convergence of the mapping gradient descent algorithm to obtain the updated seismic data;

[0104] S3. Perform one-dimensional Fourier transform on the updated seismic data to obtain frequency domain seismic data and its dimension size;

[0105] S4. Perform a forward fast Hank matrix-vector product according to the frequency domain seismic data and its dimension to obtain a matrix Y;

[0106] Specific:

[0107] S41. Obtain a randomly generated vector Ω according to the dimension of the frequency domain seismic data and a predetermined rank p;

[0108] S42. Perform one-dimensional Fourier transform on the input secondary Hank matrices c and v to obtain updated secondary Hank matrices C and V;

[0109] S43. Perform inverse Fourier transform on the updated secondary Hank matrices C and V to obtain a matrix r;

[0110] S44. truncate and reshape the matrix r to obtain the matrix Y;

[0111] S5. Perform low-rank prediction processing to obtain a low-rank approximation of the matrix Y

[0112] Specific:

[0113] S51. Perform random projection on the matrix Y to obtain the matrix M;

[0114] S52. Perform QR decomposition on the matrix M to obtain a standard orthogonal matrix Q and an upper triangular matrix R;

[0115] S53. According to the standard orthogonal matrix Q and the conjugate transposed matrix Q of the standard orthogonal matrix Q H , and obtain a low-rank approximation of the matrix Y

[0116] S6. Perform transposed fast Hank matrix-vector product on the matrix q obtained by reshaping the standard orthogonal matrix Q and the frequency domain seismic data to obtain a matrix z;

[0117] S7. According to the matrix q and the matrix z, convolution is used to accelerate the anti-diagonal average to obtain the result after low-rank prediction;

[0118] S8. Based on the results of low-rank prediction, the conjugate matrix of frequency domain seismic data is Perform inverse Fourier transform to obtain preliminary reconstructed seismic data, and define steps S3 to S8 as a mapping operation P;

[0119] S9. Repeat steps S1 to S8 until the convergence condition is reached to obtain the reconstructed seismic data.

[0120] Furthermore, in S1, the interpolation formula is expressed as:

[0121]

[0122] Among them, u(t,x,y) is the observed seismic data of the initial irregular grid, t is the sampling time of the seismic data of the initial irregular grid, (x,y) is the coordinate of the seismic data of the initial irregular grid, is the interpolation operator, and d(t,x,y) is the regular network seismic data to be reconstructed.

[0123] Furthermore, in S2, the updated seismic data m(t,x,y) is expressed as:

[0124]

[0125] Among them, d k (t,x,y) is the regular network seismic data at the kth moment, is the interpolation operator The adjoint operator, λ is the iteration step size, and the exponential decay strategy is adopted. The initial step size is Decision, e max for The maximum eigenvalue of .

[0126] Furthermore, in S3, the frequency domain seismic data D(ω,x,y) is expressed as:

[0127]

[0128] in, is the one-dimensional Fourier transform;

[0129] The dimension size (N) of frequency domain seismic data D(ω,x,y) t ,N1,N2) is expressed as:

[0130] (N t ,N1,N2)=size(D(ω,x,y))

[0131] Among them, size(·) is the function to solve the matrix size.

[0132] Further, in S41, the vector Ω is expressed as:

[0133] Ω=rand(M1,M2,p)

[0134] M1=[N1 / 2]+1

[0135] M2=[N2 / 2]+1

[0136] Among them, M1 and M2 are parameters used to specify the size of vector Ω;

[0137] In S42, the input secondary Hank matrices are c and v, and the dimension size of the matrix c is obtained;

[0138] The dimension size of matrix c (L x ,L y ) is expressed as:

[0139] (L x ,L y )=size(c)

[0140] The dimensions of the secondary Hank matrices c and v are expressed as:

[0141]

[0142]

[0143] Where N is the row size of the secondary Hank matrix, M is the column size of the secondary Hank matrix, and when flag = forward, the size of the secondary Hank matrix is ​​expressed as N c 、M c When flag = adjoint, the size of the secondary Hank matrix is ​​expressed as N r 、M r ;

[0144] The updated secondary Hank matrices C and V are expressed as:

[0145]

[0146]

[0147] Among them, padding(·) is the function of filling zero elements, reverse(·) is the inversion function, To reverse the elements of each dimension of the matrix v in sequence and then fill them with zeros, its size is [N-1,M-1];

[0148] In S43, the matrix r is expressed as:

[0149]

[0150] in, is element-wise multiplication, is the inverse Fourier transform, conj(·) is the conjugate, δ(C) is the indicator function, and its expression is:

[0151]

[0152] In S44, the matrix Y is expressed as:

[0153] Y=reshape(r[N:L x ,M:L y ],[],1)

[0154] Among them, reshape is the reshaping function, [] indicates the adaptive size of the matrix;

[0155] Specifically, steps S42 to S44 are: The specific implementation form of It is called the fast two-level Hank matrix-vector product function.

[0156] Further, in S51, according to the size N Y ×N Y The matrix Y of size N consists of p independent vectors Y ×p random matrix T to obtain the matrix M, where p<<N Y ;

[0157] The matrix M is represented as:

[0158] M=YT

[0159] In S52, the standard orthogonal matrix Q and the upper triangular matrix R are expressed as:

[0160] [Q,R]=qr(M)

[0161] Where qr(·) is the QR decomposition;

[0162] In S53, the low-rank approximation of the matrix Y It is expressed as:

[0163]

[0164] Furthermore, in S6, the standard orthogonal matrix Q is reshaped to obtain a matrix q;

[0165] The matrix q is expressed as:

[0166] q=reshape(Q,L1,L2)

[0167] L1=N1+1-M1

[0168] L2=N2+1-M2

[0169] Among them, L1 and L2 are parameters used to reshape the size of the standard orthogonal array Q;

[0170] Perform transposed fast Hank matrix-vector product on matrix q and frequency domain seismic data D(ω,x,y) to obtain matrix z;

[0171] The matrix z is represented as:

[0172]

[0173] Furthermore, in S7, the result after low-rank prediction is expressed as:

[0174]

[0175] in, is the result of the k+1th low-rank prediction, is the result of the k-th low-rank prediction.

[0176] Furthermore, in S8, the conjugate matrix of the frequency domain seismic data D(ω,x,y) The inverse Fourier transform of is expressed as:

[0177]

[0178] Among them, d k+1 is the initial reconstructed seismic data, d k+1 (t,x,y)=P(m(t,x,y)).

[0179] Furthermore, in S9, the convergence condition is that the relative error res is less than a fixed value, that is, res<1e -6 ;

[0180] The relative error res is expressed as:

[0181]

[0182] Among them, ||·|| F is the Frobenius norm.

[0183] Specifically, in this embodiment, five groups of tests are performed to analyze the reconstruction effect of the present invention by comparing the fast iterative shrinkage threshold algorithm (FISTA algorithm) and the interpolation fast low-rank prediction algorithm (I-FMSSA algorithm) proposed in the present invention;

[0184] Test 1: The FISTA algorithm and the I-FMSSA algorithm were used to operate the noise-free simulated seismic data of 50% irregular grid containing three seismic traces to obtain the reconstructed seismic data respectively;

[0185] The restoration quality is quantified according to the signal-to-noise ratio, and the signal-to-noise ratios of the FISTA algorithm and the I-FMSSA algorithm are obtained under the condition of clean synthetic data;

[0186] The signal-to-noise ratio formula is expressed as:

[0187]

[0188] Among them, ||·|| F represents the Frobenius norm, d true is the complete earthquake data, d r To reconstruct the seismic data, the running time and signal-to-noise ratio after reconstruction are shown in Table 1;

[0189] Table 1

[0190]

[0191] From Table 1, we can observe that the I-FMSSA algorithm has shorter calculation time and higher signal-to-noise ratio than the FISTA algorithm, that is, better data recovery quality;

[0192] refer to Figure 2 , a simulated seismic data with a data size of 2500 traces and irregularly distributed data with 50% missing, with a trace spacing of 20m, which is synthesized using a Ricker wavelet with a main frequency of 30Hz;

[0193] refer to Figure 3 ,It can be observed that the reconstruction effect of I-FMSSA algorithm is better;

[0194] Test 2: The FISTA algorithm and the I-FMSSA algorithm were used to operate the simulated seismic data with 50% irregular sampling containing three seismic phase axes and Gaussian noise with a signal-to-noise ratio of S / N=1, and the reconstructed seismic data were obtained respectively;

[0195] According to the signal-to-noise ratio formula, the reconstructed running time and signal-to-noise ratio obtained in Test 2 are integrated into Table 2;

[0196] Table 2

[0197]

[0198] refer to Figure 3 and Figure 4 ,And it can be observed from Table 2 that the I-FMSSA algorithm requires less running time, has clearer ,reconstruction results, higher signal-to-noise ratio, and less signal leakage,compared with the FISTA algorithm;

[0199] Test 3: Use the FISTA algorithm and the I-FMSSA algorithm to reconstruct seismic data containing complex seismic trace data, and analyze the effect of reconstructing seismic data;

[0200] According to the signal-to-noise ratio formula, the reconstructed running time and signal-to-noise ratio obtained in test three are integrated into Table 3;

[0201] Table 3

[0202]

[0203] refer to Figure 5 and Figure 6 ,It is observed that the I-FMSSA method is superior to the FISTA method, and the reconstruction results have less noise and less signal leakage, and the superiority of the I-FMSSA method in terms of accuracy and efficiency is quantitatively confirmed from Table 3;

[0204] Test 4: Use the FISTA algorithm and the I-FMSSA algorithm to reconstruct the coordinates of the real 3D seismic data acquisition shot points, and analyze the effect of reconstructing the seismic data;

[0205] refer to Figure 7 ,In the real 3D seismic data acquisition shot point coordinates, the observed irregular grid shot points include 2721 shot points, and the expected output regular grid contains 6480 shot points;

[0206] refer to Figure 8 ,For 3D seismic data, although the FISTA algorithm ,simultaneously denoises and reconstructs seismic data, the denoising effect is poor and it is ,prone to signal leakage, while the I-FMSSA algorithm reconstructs the results with clearer ,precision, and less signal leakage;

[0207] Test 5: Use the FISTA algorithm and the I-FMSSA algorithm to reconstruct the actual five-dimensional seismic data and analyze the effect of reconstructing the seismic data;

[0208] refer to Fig. 9 and Fig.10 , it is observed that the reconstruction results of the I-FMSSA method are clearer and the event axes are more continuous.

[0209] Although the present invention has been described according to a limited number of embodiments, it will be apparent to those skilled in the art, with the benefit of the above description, that other embodiments may be envisioned within the scope of the invention thus described. In addition, it should be noted that the language used in this specification is selected primarily for readability and teaching purposes, rather than for explaining or defining the subject matter of the present invention. Therefore, many modifications and variations will be apparent to those skilled in the art without departing from the scope and spirit of the appended claims. The disclosure of the present invention is illustrative, not restrictive, with respect to the scope of the present invention, which is defined by the appended claims.

Claims

1. A method for reconstructing high-dimensional seismic data from irregular grids based on fast low-rank prediction, characterized in that: The following steps are involved: S1. Regularizing and preprocessing the irregular grid seismic data by using an interpolation operator to obtain regular grid seismic data; S2. According to the regular grid seismic data, the mapping gradient descent algorithm is combined to perform iterative solution, and the exponential decay strategy is adopted for the parameter iteration step to ensure the convergence of the mapping gradient descent algorithm to obtain the updated seismic data; S3. Perform one-dimensional Fourier transform on the updated seismic data to obtain frequency domain seismic data and its dimension size; S4. Perform a forward fast Hank matrix-vector product according to the frequency domain seismic data and its dimension to obtain a matrix Y; Specific: S41. Obtain a randomly generated vector Ω according to the dimension of the frequency domain seismic data and a predetermined rank p; S42. Perform one-dimensional Fourier transform on the input secondary Hank matrices c and v to obtain updated secondary Hank matrices C and V; S43. Perform inverse Fourier transform on the updated secondary Hank matrices C and V to obtain a matrix r; S44. truncate and reshape the matrix r to obtain the matrix Y; S5. Perform low-rank prediction processing to obtain a low-rank approximation of the matrix Y Specific: S51. Perform random projection on the matrix Y to obtain the matrix M; S52. Perform QR decomposition on the matrix M to obtain a standard orthogonal matrix Q and an upper triangular matrix R; S53. According to the standard orthogonal matrix Q and the conjugate transposed matrix Q of the standard orthogonal matrix Q H , and obtain a low-rank approximation of the matrix Y S6. Perform transposed fast Hank matrix-vector product on the matrix q obtained by reshaping the standard orthogonal matrix Q and the frequency domain seismic data to obtain a matrix z; S7. According to the matrix q and the matrix z, convolution is used to accelerate the anti-diagonal average to obtain the result after low-rank prediction; S8. Based on the results of low-rank prediction, the conjugate matrix of frequency domain seismic data is Perform inverse Fourier transform to obtain preliminary reconstructed seismic data, and define steps S3 to S8 as a mapping operation P; S9. Repeat steps S1 to S8 until convergence conditions are reached to obtain reconstructed seismic data; In S1, the interpolation formula is expressed as: Among them, u(t,x,y) is the observed seismic data of the initial irregular grid, t is the sampling time of the seismic data of the initial irregular grid, (x,y) is the coordinate of the seismic data of the initial irregular grid, is the interpolation operator, d(t,x,y) is the regular network seismic data; In S2, the updated seismic data m(t,x,y) is expressed as: Among them, d k (t,x,y) is the regular network seismic data at the kth moment, is the interpolation operator The adjoint operator, λ is the iteration step size, and the exponential decay strategy is adopted. The initial step size is Decision, e max for The maximum eigenvalue of In S41, the vector Ω is expressed as: Ω=rand(M1,M2,p) M1=[N1 / 2]+1 M2=[N2 / 2]+1 Among them, M1 and M2 are parameters used to specify the size of the vector Ω. The dimension size of the frequency domain seismic data D(ω,x,y) is (N t ,N1,N2); In S42, the input secondary Hank matrices are c and v, and the dimension size of the matrix c is obtained; The dimension size of matrix c (L x ,L y ) is expressed as: (: x ,L y )=size(c) The dimensions of the secondary Hank matrices c and v are expressed as: Where N is the row size of the secondary Hank matrix, M is the column size of the secondary Hank matrix, and when flag = forward, the size of the secondary Hank matrix is ​​expressed as N c 、M c When flag = adjoint, the size of the secondary Hank matrix is ​​expressed as N r 、M r ; The updated secondary Hank matrices C and V are expressed as: Among them, padding(·) is the function of filling zero elements, reverse(·) is the inversion function, To reverse the elements of each dimension of the matrix v in sequence and fill them with zeros, its size is [N-1, M-1], is the one-dimensional Fourier transform; In S43, the matrix r is expressed as: in, is element-wise multiplication, is the inverse Fourier transform, conj(·) is the conjugate, δ(C) is the indicator function, and its expression is: In S44, the matrix Y is expressed as: Y=reshape(r[N:L x ,M:L y ],[],1) Among them, reshape is the reshaping function, [] indicates the adaptive size of the matrix; In S8, the conjugate matrix of the frequency domain seismic data D(ω,x,y) The inverse Fourier transform of is expressed as: Among them, d k+1 is the initial reconstructed seismic data, d k+1 (t,x,y)=P(m(t,x,y)).

2. The method for reconstructing irregular grid high-dimensional seismic data based on fast low-rank prediction according to claim 1, characterized in that: In S3, the frequency domain seismic data D(ω,x,y) is expressed as: in, is the one-dimensional Fourier transform; The dimension size (N) of frequency domain seismic data D(ω,x,y) t ,N1,N2) is expressed as: (N t ,N1,N2)=size(D(ω,x,y)) Among them, size(·) is the function to solve the matrix size.

3. The method for reconstructing irregular grid high-dimensional seismic data based on fast low-rank prediction according to claim 2, characterized in that: In S51, according to the size N Y ×N Y The matrix Y of size N consists of p independent vectors Y ×p random matrix T to obtain the matrix M, where p<<N Y ; The matrix M is represented as: M=YT In S52, the standard orthogonal matrix Q and the upper triangular matrix R are expressed as: [Q,R]=qr(M) Where qr(·) is the QR decomposition; In S53, the low-rank approximation of the matrix Y It is expressed as:

4. The method for reconstructing irregular grid high-dimensional seismic data based on fast low-rank prediction according to claim 3, characterized in that: In S6, the standard orthogonal matrix Q is reshaped to obtain a matrix q; The matrix q is expressed as: q=reshape(Q,L1,L2) L1=N1+1-M1 L2=N2+1-M2 Among them, L1 and L2 are parameters used to reshape the size of the standard orthogonal array Q; Perform transposed fast Hank matrix-vector product on matrix q and frequency domain seismic data D(ω,x,y) to obtain matrix z; The matrix z is represented as:

5. The method for reconstructing irregular grid high-dimensional seismic data based on fast low-rank prediction according to claim 4, characterized in that: In S7, the result after low-rank prediction is expressed as: in, is the result of the k+1th low-rank prediction, is the result of the k-th low-rank prediction.

6. The method for reconstructing irregular grid high-dimensional seismic data based on fast low-rank prediction according to claim 5, characterized in that: In S9, the convergence condition is that the relative error res is less than a fixed value, that is, res<1e -6 ; The relative error res is expressed as: Among them, ||·|| F is the Frobenius norm.

Images