Seismic data denoising method utilizing cyclic matrix convolution acceleration

Abstract

The invention discloses a seismic data denoising method utilizing cyclic matrix convolution acceleration, and relates to the technical field of oil exploration, and the method comprises the following steps: s1, inputting noise-containing three-dimensional seismic data; s2, constructing a conical curved surface area division operator L based on conical propagation geometric characteristics of noise in a three-dimensional space, and mapping d to a conical curved surface area division domain m; s3, constructing an objective function, and solving by adopting a rapid optimization algorithm to obtain a three-dimensional data volume seismic channel with noise removed; s4, judging whether all the data volumes are denoised or not, and if not, returning to the step s2; and if so, outputting the three-dimensional seismic data volume with the noise removed, and ending the noise removal. According to the efficient seismic data denoising method, the cyclic matrix characteristic of a transformation operator is adopted, and a transformation operator matrix originally formed by integral elements is converted into a product form; compared with a conventional method, the calculation efficiency is remarkably improved.

Description

Technical Field

[0001] This invention relates to the field of petroleum exploration technology, and in particular to a method for denoising seismic data using cyclic matrix convolution acceleration. Background Technology

[0002] Surface waves, as typical high-energy, low-frequency, and long-duration coherent noise in seismic exploration, exhibit conical dispersion in 3D records, significantly reducing data resolution and signal-to-noise ratio, and causing serious problems for subsequent interpretation. Therefore, effective suppression of surface waves is a crucial step in data processing. Traditional τ-p-type Radon methods use plane wave assumptions as basis functions, which are mismatched with the geometric shape of conical surface waves, making it difficult to focus energy and resulting in significant denoising residuals. In recent years, researchers have proposed a generalized Radon transform with a conical surface as its core, which can map conical noise to extremely narrow energy clusters, achieving high-precision separation. However, this transform operator is highly nonlinear, and the operator-data coupling dimension is huge. When directly solving on industrial-grade 3D volumes, the storage and computational costs expand superlinearly, and real-time processing still faces severe challenges.

[0003] Therefore, there is an urgent need to study an efficient method for denoising seismic data. Summary of the Invention

[0004] To address the aforementioned problems, this invention discloses a seismic data denoising method that utilizes cyclic matrix convolution for acceleration. Compared to traditional methods, this invention's method improves computational efficiency by more than 1000 times and increases resolution.

[0005] To achieve the above objectives, the present invention adopts the following technical solution:

[0006] A seismic data denoising method using cyclic matrix convolution for acceleration includes the following steps:

[0007] s1. Input noisy 3D seismic data;

[0008] s2. Based on the conical propagation geometric characteristics of noise in three-dimensional space, a conical surface integral operator L is constructed to map the time-domain three-dimensional seismic data d to the conical surface integral domain m;

[0009] s3. Construct the objective function and use a fast optimization algorithm to obtain the noise-removed 3D data volume seismic trace;

[0010] s4. Determine if all data volumes have been denoised. If not, return to step s2. If they have been denoised, output the denoised 3D seismic data volume and the denoising process ends.

[0011] Optionally, in step s2, the expression for the conical surface integral is:

[0012] (1);

[0013] In the formula, This is 3D seismic data, where t is the recording time and (x,y) are the spatial coordinates. These are the positive transform coefficients. For intercept time, Let x be the slowness in the x-direction of the survey line. Let be the slowness in the y-direction of the measurement line. Represents the Dirac function, These are basis functions.

[0014] Optionally, step s3 specifically includes:

[0015] s31. Perform a one-dimensional Fourier transform on both sides of equation (1) with respect to t:

[0016] (2);

[0017] In the formula, d represents the three-dimensional seismic data in the frequency domain. Let ω be the frequency domain basis function, i be the angular frequency, and i be the imaginary unit.

[0018] s32. Introduce the proportional relationship between ray parameters p and q, that is, let:

[0019] p = k·q (3);

[0020] In the formula, k∈R + It is a positive integer;

[0021] Under this condition, the original matrix is ​​transformed into the following form:

[0022] (4);

[0023] In the formula, ω is the angular frequency, j is the imaginary unit which remains constant in each row of the matrix, and the offsets x and y extend from 0 to X. max and Y max At this point, matrix A has a Hermitian Toeplitz structure, where the matrix elements are:

[0024] (5);

[0025] In the formula, It is the element of matrix A located in the p-th row and q-th column. Let A be the element located in the (p+1)th row and (q+1)th column of matrix A;

[0026] Define matrix element a n for:

[0027] (6);

[0028] In the formula, n=pq;

[0029] matrix It can be written in the form of a cyclic matrix as follows:

[0030] (7);

[0031] In the formula, circ(·) represents a cyclic matrix with dimensions N*N, and the matrix elements are constants along each diagonal, depending only on the index difference n=pq; since this matrix is ​​a cyclic matrix, its product with any vector or matrix V is equivalent to the first column of the cyclic matrix. Circular convolution with V:

[0032] (8);

[0033] According to the convolution theorem, circular convolution in the time domain corresponds to multiplication in the frequency domain:

[0034] (9);

[0035] In the formula, Represents the Discrete Fourier Transform. This indicates element-wise multiplication;

[0036] therefore:

[0037] (10);

[0038] This property can be used to significantly accelerate matrix-vector multiplication operations. Based on this, the iterative steps in the original reconstruction algorithm are accelerated; before applying the fast algorithm, the iterative update formula was:

[0039] (11);

[0040] After applying the fast algorithm, the iterative update formula becomes:

[0041] (12);

[0042] In the formula, Let be the conjugate transpose of the conical surface integral operator, m be the conical surface integral domain, and D be the time-domain 3D seismic data. Represents the soft threshold function. Here, T is the threshold parameter, T is the iteration step size, and B represents the number of iterations. To update the step size, y B x B As an intermediate variable;

[0043] (13);

[0044] The above fast algorithm is used to solve formula (13) to obtain m.

[0045] s33. Combine m with Multiplying these results in a noise-removed 3D data volume seismic trace.

[0046] The beneficial effects of this invention are that the efficient seismic data denoising method of this invention uses the cyclic matrix characteristics of the transformation operator to transform the transformation operator matrix, which is originally composed of integral elements, into a product form; compared with conventional methods, it significantly improves computational efficiency. Attached Figure Description

[0047] Figure 1 This is a flowchart illustrating a seismic data denoising method accelerated by cyclic matrix convolution according to the present invention.

[0048] Figure 2 It is a volume of noisy 3D seismic data;

[0049] Figure 3 This is a zero-intercept time slice obtained using the three-dimensional seismic data noise suppression method of the present invention;

[0050] Figure 4 This is a zero-intercept time slice obtained using traditional methods for noise suppression of 3D seismic data;

[0051] Figure 5 This is a three-dimensional data volume with suppressed noise obtained using the three-dimensional seismic data noise suppression method of the present invention;

[0052] Figure 6 The noise is separated using the three-dimensional seismic data noise suppression method of this invention;

[0053] Figure 7 This is a noise-suppressed 3D data volume obtained using traditional 3D seismic data noise suppression methods;

[0054] Figure 8 This refers to the noise separation obtained using traditional methods for noise suppression in 3D seismic data. Detailed Implementation

[0055] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, not all of them. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention. Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to represent selected embodiments of the invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0056] A seismic data denoising method that utilizes cyclic matrix convolution for acceleration, such as... Figure 1 As shown, it includes the following steps:

[0057] s1. Input noisy 3D seismic data.

[0058] s2. Based on the conical propagation geometry of noise in three-dimensional space, a conical surface integral operator L is constructed to map the time-domain three-dimensional seismic data d to the conical surface integral domain m; the expression of the conical surface integral is:

[0059] (1);

[0060] In the formula, This is 3D seismic data, where t is the recording time and (x,y) are the spatial coordinates. These are the positive transform coefficients. For intercept time, Let x be the slowness in the x-direction of the survey line. Let be the slowness in the y-direction of the measurement line. Represents the Dirac function, These are basis functions.

[0061] s3. Construct the objective function and use a fast optimization algorithm to obtain the noise-removed 3D data volume seismic trace; specifically including:

[0062] s31. Perform a one-dimensional Fourier transform on both sides of equation (1) with respect to t:

[0063] (2);

[0064] In the formula, d represents the three-dimensional seismic data in the frequency domain. Let ω be the frequency domain basis function, i be the angular frequency, and i be the imaginary unit.

[0065] s32. Introduce the proportional relationship between ray parameters p and q, that is, let:

[0066] p = k·q (3);

[0067] In the formula, k∈R + It is a positive integer;

[0068] Under this condition, the original matrix is ​​transformed into the following form:

[0069] (4);

[0070] In the formula, ω is the angular frequency, j is the imaginary unit which remains constant in each row of the matrix, and the offsets x and y extend from 0 to X. max and Y max At this point, matrix A has a Hermitian Toeplitz structure, where the matrix elements are:

[0071] (5);

[0072] In the formula, It is the element of matrix A located in the p-th row and q-th column. Let A be the element located in the (p+1)th row and (q+1)th column of matrix A;

[0073] Define matrix element a n for:

[0074] (6);

[0075] In the formula, n=pq;

[0076] matrix It can be written in the form of a cyclic matrix as follows:

[0077] (7);

[0078] In the formula, circ(·) represents a cyclic matrix with dimensions N*N, and the matrix elements are constants along each diagonal, depending only on the index difference n=pq; since this matrix is ​​a cyclic matrix, its product with any vector or matrix V is equivalent to the first column of the cyclic matrix. Circular convolution with V:

[0079] (8);

[0080] According to the convolution theorem, circular convolution in the time domain corresponds to multiplication in the frequency domain:

[0081] (9);

[0082] In the formula, Represents the Discrete Fourier Transform. This indicates element-wise multiplication;

[0083] therefore:

[0084] (10);

[0085] This property can be used to significantly accelerate matrix-vector multiplication operations. Based on this, the iterative steps in the original reconstruction algorithm are accelerated; the iterative update formula before applying the fast algorithm is:

[0086] (11);

[0087] After applying the fast algorithm, the iterative update formula becomes:

[0088] (12);

[0089] In the formula, Let be the conjugate transpose of the conical surface integral operator, m be the conical surface integral domain, and D be the time-domain 3D seismic data. Represents the soft threshold function. Here, T is the threshold parameter, T is the iteration step size, and B represents the number of iterations. To update the step size, y B x B As an intermediate variable;

[0090] (13);

[0091] The above fast algorithm is used to solve formula (13) to obtain m.

[0092] s33. Combine m with Multiplying these results in a noise-removed 3D data volume seismic trace.

[0093] s4. Determine if all data volumes have been denoised. If not, return to step s2. If they have been denoised, output the denoised 3D seismic data volume and the denoising process ends.

[0094] This invention discloses a seismic data denoising method that utilizes cyclic matrix convolution for acceleration. When applied to 3D seismic data volumes, it achieves ideal computational results.

[0095] Figure 2 For a noisy 3D seismic data volume, from Figure 2 As can be seen, in addition to the reflected waves, there is also a large amount of noise. Figure 3 This is a zero-intercept time slice obtained using the three-dimensional seismic data noise suppression method of this invention; as a comparison, Figure 4The diagram shows a zero-intercept time slice obtained by a conventional method. A comparison reveals that the method of this invention yields a higher resolution and a more clearly visible effective signal. Therefore, the efficient seismic data denoising method utilizing cyclic matrices of this invention has higher resolution and faster computation speed in the slow-motion direction compared to conventional methods.

[0096] Figure 5 This is a three-dimensional data volume with suppressed noise obtained using the method of the present invention. Figure 6 The separated noise obtained using the method of the present invention; as a comparison, Figure 7 This shows a three-dimensional data volume with suppressed noise obtained using conventional methods. Figure 8 The separated noise obtained using conventional methods is shown. By comparison, it is demonstrated that the method of the present invention achieves better results in eliminating three-dimensional noise and reducing reflected wave loss, and has a faster calculation speed.

[0097] Table 1 compares the computation time of the present invention with that of the traditional method in a GPU environment. It can be seen that the method of the present invention significantly improves computational efficiency when processing the same set of data.

[0098] Table 1. Time comparison between the method of the present invention and the traditional method in a GPU environment.

[0099]

[0100] The method of this invention utilizes the cyclic matrix property of the cone transformation operator to transform the cone transformation operator matrix, which was originally composed of integral elements, into a product form; compared with the conventional method, the computation time is reduced by 3864.4 times, significantly improving computational efficiency.

[0101] Of course, the above description is not intended to limit the present invention, and the present invention is not limited to the examples given above. Any changes, modifications, additions or substitutions made by those skilled in the art within the scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. A seismic data denoising method using cyclic matrix convolution for acceleration, characterized in that, Includes the following steps: s1. Input noisy 3D seismic data; s2. Based on the conical propagation geometric characteristics of noise in three-dimensional space, a conical surface integral operator L is constructed to map the time-domain three-dimensional seismic data d to the conical surface integral domain m; s3. Construct the objective function and use a fast optimization algorithm to obtain the noise-removed 3D data volume seismic trace; s4. Determine if all data volumes have been denoised. If not, return to step s2. If they have been denoised, output the denoised 3D seismic data volume and the denoising process ends.

2. The seismic data denoising method using cyclic matrix convolution acceleration as described in claim 1, characterized in that, In step s2, the surface integral expression for the cone is: （1）； In the formula, This is 3D seismic data, where t is the recording time and (x,y) are the spatial coordinates. These are the positive transform coefficients. For intercept time, Let x be the slowness in the x-direction of the survey line. Let be the slowness in the y-direction of the measurement line. Represents the Dirac function, These are basis functions.

3. The seismic data denoising method using cyclic matrix convolution for acceleration as described in claim 2, characterized in that, Step s3 specifically includes: s31. Perform a one-dimensional Fourier transform on both sides of equation (1) with respect to t: （2）； In the formula, d represents the three-dimensional seismic data in the frequency domain. Let ω be the frequency domain basis function, i be the angular frequency, and i be the imaginary unit. s32. Introduce the proportional relationship between ray parameters p and q, that is, let: p = k·q (3); In the formula, k∈R + It is a positive integer; Under this condition, the original matrix is ​​transformed into the following form: （4）； In the formula, j is the imaginary unit, which remains constant in each row of the matrix, and the offsets x and y extend from 0 to X, respectively. max and Y max At this point, matrix A has a Hermitian Toeplitz structure, where the matrix elements are: （5）； In the formula, It is the element of matrix A located in the p-th row and q-th column. Let A be the element located in the (p+1)th row and (q+1)th column of matrix A; Define matrix element a n for: （6）； In the formula, n=pq; matrix It can be written in the form of a cyclic matrix as follows: （7）； In the formula, circ(·) represents a cyclic matrix with dimensions N*N, and the matrix elements are constants along each diagonal, depending only on the index difference n=pq; since this matrix is ​​a cyclic matrix, its product with any vector or matrix V is equivalent to the first column of the cyclic matrix. Circular convolution with V: （8）； According to the convolution theorem, circular convolution in the time domain corresponds to multiplication in the frequency domain: （9）； In the formula, Represents the Discrete Fourier Transform. This indicates element-wise multiplication; therefore: （10）； Before applying the fast algorithm, the iterative update formula was: （11）； After applying the fast algorithm, the iterative update formula becomes: （12）； In the formula, Let be the conjugate transpose of the conical surface integral operator, m be the conical surface integral domain, and D be the time-domain 3D seismic data. Represents the soft threshold function. Here, T is the threshold parameter, T is the iteration step size, and B represents the number of iterations. To update the step size, y B x B As an intermediate variable; （13）； The above fast algorithm is used to solve formula (13) to obtain m; s33. Combine m with Multiplying these results in a noise-removed 3D data volume seismic trace.

Images