A seismic data denoising method accelerated by cyclic matrix convolution

Abstract

The application discloses a kind of seismic data denoising method using cyclic matrix convolution acceleration, it is related to oil exploration technical field, including the following steps: s1. input noisy three-dimensional seismic data;S2. based on the conical propagation geometric characteristics of noise in three-dimensional space, construct conical curved area integral operator L, d is mapped to conical curved area integral domain m;S3. build objective function, solve and obtain three-dimensional data body seismic trace of removing noise using fast optimization algorithm;S4. whether all data bodies are denoised, if not completed, return to step s2;If completed, output three-dimensional seismic data body of removing noise, remove noise end.The application efficient seismic data denoising method, using the cyclic matrix characteristics of transformation operator, the transformation operator matrix originally formed by integral element is converted into product form;Compared with conventional method, the calculation efficiency is significantly improved.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of oil exploration, and particularly relates to a seismic data denoising method accelerated by cyclic matrix convolution. BACKGROUND

[0002] As a typical strong energy, low frequency and long-lasting coherent noise in seismic exploration, the surface wave is dispersed in a conical shape in three-dimensional recording, significantly reducing the data resolution and signal-to-noise ratio, and causing serious problems to subsequent interpretation, so its effective suppression becomes a key link in data processing. The traditional Tau-p type Radon method takes the plane wave assumption as the basis function, which is mismatched with the geometric shape of the conical surface wave, so that the energy is difficult to focus, and the denoising residual is significant. In recent years, researchers have proposed a generalized Radon transform with a conical curved surface as the kernel, which can map the conical noise to a very narrow energy group, realizing high-precision separation. However, the transform operator is highly nonlinear, and the operator-data coupling dimension is huge, so when directly solving on an industrial three-dimensional body, the storage and calculation amount is super-linearly expanded, and real-time processing still faces severe challenges.

[0003] Therefore, it is urgent to study an efficient seismic data denoising method. SUMMARY

[0004] To solve the above problems, the present application discloses a seismic data denoising method accelerated by cyclic matrix convolution, which has a calculation efficiency improved by more than 1000 times compared with the traditional method and increased resolution.

[0005] To achieve the above purpose, the present application adopts the following technical scheme:

[0006] A seismic data denoising method accelerated by cyclic matrix convolution, comprising the following steps:

[0007] s1. inputting three-dimensional seismic data containing noise;

[0008] s2. based on the conical propagation geometric characteristics of the noise in the three-dimensional space, constructing a conical curved surface integral operator L, and mapping the three-dimensional seismic data d in the time domain to the conical curved surface integral domain m;

[0009] s3. constructing an objective function, and solving to obtain a three-dimensional data body seismic trace without noise by using a fast optimization algorithm;

[0010] s4. judging whether all data bodies are denoised, if not, returning to step s2; if yes, outputting the three-dimensional seismic data body without noise, and ending the denoising.

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

[0012] (1);

[0013] where, is the three-dimensional seismic data, t is the recording time, (x, y) is the spatial coordinate value, is the positive transform coefficient, is the intercept time, is the slowness in the x-line direction, is the slowness in the y-line direction, denotes the Dirac function, is the base function.

[0014] Optionally, the step s3 specifically comprises:

[0015] s31. One-dimensional Fourier transform is performed on both sides of equation (1) with respect to t:

[0016] (2);

[0017] where, d is the three-dimensional seismic data in the frequency domain, is the base function in the frequency domain, ω is the angular frequency, and i is the imaginary unit;

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

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

[0020] where, k∈R + is a positive integer;

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

[0022] (4);

[0023] where, ω is the angular frequency, j is the imaginary unit, which remains constant in each row of the matrix, the range of offset x and y respectively extends from 0 to X max and Y max , at this time, the matrix A has a Hermitian Toeplitz structure, and the matrix element is:

[0024] (5);

[0025] where, is the element of matrix A located in the pth row and the qth column, is the element of matrix A located in the (p+1)th row and the (q+1)th column;

[0026] The matrix element a n is defined as:

[0027] (6);

[0028] where n = p - q;

[0029] matrix is written as a circulant matrix as follows:

[0030] (7);

[0031] where circ(·) denotes a circulant matrix with dimension N*N, and the matrix elements are constant along each diagonal, depending only on the index difference n = p - q; since the matrix is a circulant matrix, its product with any vector or matrix V is equivalent to the first column of the circulant matrix and the circulant convolution with V:

[0032] (8);

[0033] According to the convolution theorem, the circulant convolution in the time domain corresponds to the product in the frequency domain:

[0034] (9);

[0035] where denotes the discrete Fourier transform, denotes element-wise multiplication;

[0036] Therefore:

[0037] (10);

[0038] Using this property, the multiplication operation of a matrix and a vector can be significantly accelerated. On this basis, the iteration step in the original reconstruction algorithm is accelerated; before applying the fast algorithm, the iteration update formula is:

[0039] (11);

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

[0041] (12);

[0042] where is the conjugate transpose matrix of the conical curved surface integral operator, m is the conical curved surface integral domain, D is the three-dimensional seismic data in the time domain, denotes the soft threshold function, is the threshold parameter, T is the iteration step, B represents the number of iterations, is the update step, y B , x B are intermediate variables;

[0043] (13);

[0044] Solve formula (13) by using the above fast algorithm to obtain m.

[0045] s33. Multiply m with to obtain the three-dimensional data body seismic trace after removing noise.

[0046] The beneficial effects of the present application are that the efficient seismic data denoising method of the present application uses the cyclic matrix characteristics of the transform operator to convert the transform operator matrix originally composed of integral elements into a product form; compared with conventional methods, the calculation efficiency is significantly improved. BRIEF DESCRIPTION OF DRAWINGS

[0047] Figure 1 is a flowchart of a seismic data denoising method using cyclic matrix convolution acceleration according to the present application;

[0048] Figure 2 is a three-dimensional seismic data body containing noise;

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

[0050] Figure 4 is a zero-intercept time slice obtained by using the conventional three-dimensional seismic data noise suppression method;

[0051] Figure 5 is a three-dimensional data body after noise suppression obtained by using the three-dimensional seismic data noise suppression method of the present application;

[0052] Figure 6 is separated noise obtained by using the three-dimensional seismic data noise suppression method of the present application;

[0053] Figure 7 is a three-dimensional data body after noise suppression obtained by using the conventional three-dimensional seismic data noise suppression method;

[0054] Figure 8 is separated noise obtained by using the conventional three-dimensional seismic data noise suppression method. DETAILED DESCRIPTION

[0055] In order to make the purposes, technical solutions and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the drawings in the embodiments of the present application. Obviously, the described embodiments are only part of the embodiments of the present application, rather than all the embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by a person of ordinary skill in the art without creative labor fall within the protection scope of the present application. Therefore, the following detailed description of the embodiments of the present application provided in the drawings is not intended to limit the scope of the claimed present application, but only represents selected embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by a person of ordinary skill in the art without creative labor fall within the protection scope of the present application.

[0056] A seismic data denoising method accelerated by cyclic matrix convolution, as shown in Figure 1 , comprises the following steps:

[0057] s1. Input the three-dimensional seismic data containing noise.

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

[0059] (1);

[0060] In the formula, is three-dimensional seismic data, t is recording time, (x, y) is spatial coordinate value, is a positive transform coefficient, is an intercept time, is the slowness in the x-line direction, is the slowness in the y-line direction, represents the Dirac function, is a basis function.

[0061] s3. Constructing a target function, using a fast optimization algorithm to obtain the three-dimensional data body seismic trace after removing noise; specifically including:

[0062] s31. One-dimensional Fourier transform is performed on both sides of equation (1) with respect to t:

[0063] (2);

[0064] In the formula, d is three-dimensional seismic data in the frequency domain, is a basis function in the frequency domain, ω is an angular frequency, and i is an imaginary unit;

[0065] s32. Introducing a proportional relationship between the ray parameters p and q, i.e. let:

[0066] p = k q (3);

[0067] where k e R + is a positive integer;

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

[0069] (4);

[0070] where w is the angular frequency and j is the imaginary unit, which remains constant in each row of the matrix, the ranges of the offsets x and y extend from 0 to X max and Y max respectively, at this time, the matrix A has a Hermitian Toeplitz structure, where the matrix elements are:

[0071] (5);

[0072] where is the element of the matrix A located in the pth row and the qth column, is the element of the matrix A located in the (p+1)th row and the (q+1)th column;

[0073] The matrix element a n is defined as:

[0074] (6);

[0075] where n = p - q;

[0076] The matrix is written in the following circulant matrix form:

[0077] (7);

[0078] where circ(·) represents a circulant matrix with a dimension of N*N, and the matrix elements along each diagonal line are constant, only depending on the index difference n = p - q; Since this matrix is a circulant matrix, its product with any vector or matrix V is equivalent to the first column of the circulant matrix and the circulant convolution of V:

[0079] (8);

[0080] According to the convolution theorem, the circulant convolution in the time domain corresponds to the product in the frequency domain:

[0081] (9);

[0082] where denotes a discrete Fourier transform, denotes an element-wise multiplication;

[0083] Therefore:

[0084] (10);

[0085] With this property, the multiplication operation of matrix and vector can be significantly accelerated. On this basis, the iteration step in the original reconstruction algorithm is accelerated; the iteration update formula before applying the fast algorithm is:

[0086] (11);

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

[0088] (12);

[0089] In the formula, is the conjugate transpose matrix of the conical curved surface integral operator, m is the conical curved surface integral domain, D is the three-dimensional seismic data in the time domain, denotes a soft threshold function, is a threshold parameter, T is an iteration step, B represents the number of iterations, is an update step, y B , x B are intermediate variables;

[0090] (13);

[0091] The formula (13) is solved by using the above fast algorithm to obtain m.

[0092] s33. Multiply m with to obtain a three-dimensional data body seismic trace with noise removed.

[0093] s4. Determine whether all data bodies are denoised, if not completed, return to step s2; if completed, output the three-dimensional seismic data body with noise removed, and the noise removal is completed.

[0094] The application discloses a seismic data denoising method using cyclic matrix convolution acceleration, which is applied to three-dimensional seismic data body, and ideal calculation effect is achieved.

[0095] Figure 2 is a three-dimensional seismic data body with noise, from Figure 2 It can be seen from the formula that, in addition to the reflected wave, there is a large amount of noise. Figure 3 is a zero-intercept time slice obtained by using the three-dimensional seismic data noise suppression method; as a comparison, Figure 4The zero-intercept time slice obtained by the traditional method is shown, and it can be seen by comparison that the resolution of the result obtained by the method of the application is higher, and the effective signal is more obvious. Therefore, the method for efficiently denoising seismic data by using the cyclic matrix of the application has higher resolution in the slowness direction than the traditional method, and the calculation speed is faster.

[0096] Figure 5 The three-dimensional data body obtained by suppressing noise by the method of the application, Figure 6 The separated noise obtained by the method of the application; as a comparison, Figure 7 The three-dimensional data body obtained by suppressing noise by the traditional method is shown, Figure 8 The separated noise obtained by the traditional method is shown, and it is proved by comparison that the method of the application has better effect in eliminating three-dimensional noise and reducing reflection loss, and the calculation speed is faster.

[0097] Table 1 is a comparison of the calculation time of the application and the traditional method in the GPU environment, and it can be seen that the calculation efficiency of the method of the application is significantly improved when processing the same set of data.

[0098] Table 1 is a comparison of the calculation time of the application and the traditional method in the GPU environment,

[0099]

[0100] The method of the application adopts the cyclic matrix characteristics of the cone transform operator, and converts the cone transform operator matrix originally composed of integral elements into a product form; compared with the conventional method, the calculation time is reduced by 3864.4 times, and the calculation efficiency is significantly improved.

[0101] Of course, the above description is not a limitation of the application, and the application is not limited to the above examples, and the changes, modifications, additions or replacements made by the person skilled in the art within the essential scope of the application should also belong to the protection scope of the application.

Claims

1. A method for seismic data denoising accelerated by circulant matrix convolution, characterized in that, The method comprises the following steps: s1. inputting three-dimensional seismic data containing noise; s2. constructing a conical curved surface integral operator L based on the conical propagation geometry of noise in three-dimensional space, and mapping the three-dimensional seismic data D in time domain to a conical curved surface integral domain m; s3. constructing an objective function, and solving the three-dimensional data body seismic trace after noise removal by using a fast optimization algorithm; s4. judging whether all data bodies are denoised, if not, returning to step s2; if yes, outputting the three-dimensional seismic data body after noise removal, and ending the noise removal; In step s2, the conical curved surface integral expression is as follows: （1）； wherein is the time domain three-dimensional seismic data, t is the recording time, (x, y) are spatial coordinate values, is the positive transform coefficient, is the intercept time, is the slowness in the x-line direction, is the slowness in the y-line direction, denotes the Dirac function, is the basis function; Step s3 specifically comprises: s31. performing one-dimensional Fourier transform on t on both sides of equation (1): （2）； where d is the three-dimensional seismic data in the frequency domain, is the base function in the frequency domain, ω is the angular frequency, and i is the imaginary unit. s32. introducing the proportional relationship between the ray parameters p and q, that is, let: p=k·q (3); where k e R + is a positive integer; The original matrix is transformed into the following form: （4）； where j is the imaginary unit, which remains constant in each row of the matrix, the offset distances x and y range from 0 to X, respectively max and Y max At this time, the matrix A has a Hermitian Toeplitz structure, where the matrix elements are: （5）； wherein is the element of matrix A located in the pth row and qth column, is the element of matrix A located in the p+1th row and q+1th column; Definition of matrix element a n is: （6）； In the formula, n=p-q; matrix written as a circulant matrix of the form: （7）； where circ(•) denotes a circulant matrix of dimension N*N with constant matrix elements along each diagonal, depending only on the index difference n = p - q; since this matrix is a circulant matrix, its product with any vector or matrix V is equivalent to the circular convolution of the first column of the circulant matrix with V Circulant convolution with V: （8）； According to the convolution theorem, the circular convolution in the time domain corresponds to the product in the frequency domain: （9）； wherein denotes a discrete Fourier transform, denotes element-wise multiplication; Therefore: （10）； Before applying the fast algorithm, the iterative update formula is: （11）； After applying the fast algorithm, the iterative update formula becomes: （12）； In the formula, is the conjugate transpose matrix of the conical curved surface integral operator, m is the conical curved surface integral domain, D is the three-dimensional seismic data in the time domain, represents a soft threshold function, is a threshold parameter, T is an iteration step, and B represents the number of iterations, is an update step, y B , x B is an intermediate variable; （13）； The formula (13) is solved by using the above fast algorithm to obtain m; s33. multiplying m with to obtain a noise-removed 3D data volume seismic trace.

Images