Seismic exploration data regularization processing method and device

Abstract

The invention discloses a seismic exploration data regularization processing method and device. The method comprises the following steps: converting to obtain a wave number spectrum data volume based on time-space domain seismic data; carrying out the following regularization processing on the wave number spectrum of each frequency slice in the wave number spectrum data volume: determining the maximum amplitude in the amplitude spectrum of the wave number spectrum of the current frequency slice according to the inclination angle weight function, carrying out the least square optimization on the wave number component corresponding to the maximum amplitude, and updating the wave number spectrum of the current frequency slice in combination with a Grubrum matrix; returning to continue to execute the step of the maximum amplitude until the preset number of iterations is reached, and obtaining the wave number spectrum after the regularization processing of the current frequency slice; and obtaining a regularized wave number spectrum data body based on the regularized wave number spectrums of all the frequency slices, and converting to obtain time-space domain seismic data. According to the method, the reconstruction result is accurate, the pseudo Gibbs phenomenon is improved, the spatial sampling attribute of the data is improved, and the seismic exploration data imaging quality is improved.

Description

Technical Field

[0001] The present invention relates to the technical field of seismic exploration data processing, and particularly to a method and device for regularizing seismic exploration data. Background Art

[0002] In current seismic exploration, the collected data often fails to meet the requirements of subsequent processing and imaging for the spatial regularity and spatial sampling density of seismic data. Sparse or irregular spatial sampling attributes will affect the inversion process and migration imaging effect of seismic data, while seismic data regularization technology can, to a certain extent, improve the spatial sampling attributes of the seismic data observation system, thereby improving the inversion and imaging effects of seismic data.

[0003] Currently, the commonly used seismic data regularization methods include algorithms such as the anti-leakage Fourier transform (ALFT). The ALFT algorithm is a high-dimensional seismic data regularization algorithm, which has better amplitude fidelity and can, to a certain extent, complement data gaps. Due to its good application effect and easy high-dimensional implementation, it has become the mainstream regularization method in the industry. Summary of the Invention

[0004] The inventors of the present application found that although the traditional ALFT algorithms can complement data gaps to a certain extent, they are generally only applicable to the regularization of uniformly and randomly spatially sampled seismic data. When there are regularly missing data and large gaps in seismic data, their processing effects are not ideal. On the one hand, due to the influence of periodic aliasing, reconstruction artifacts may occur and the data reconstruction results are inaccurate; on the other hand, pseudo-Gibbs phenomena will occur at the discontinuity points of the event axis. As a result, the processing effect of seismic data is poor, affecting the subsequent imaging quality and effect, and there are distortion phenomena in the regularized seismic data and imaging.

[0005] In view of the above problems, the present invention is proposed to provide a method and device for regularizing seismic exploration data that overcome the above problems or at least partially solve the above problems.

[0006] An embodiment of the present invention provides a method for regularizing seismic exploration data, including:

[0007] Based on the obtained time-space domain seismic data, convert to obtain the wavenumber spectrum data volume of the seismic data;

[0008] Perform the following regularization process on the wavenumber spectrum of each frequency slice in the wavenumber spectrum data volume:

[0009] According to the dip weight function, determine the maximum amplitude in the amplitude spectrum of the wavenumber spectrum of the current frequency slice; the dip weight function is determined according to the amplitude energy of the wavenumber spectrum data volume in different ray directions;

[0010] Perform least - squares optimization on the wavenumber component corresponding to the maximum amplitude according to the Gram matrix, and update the wavenumber spectrum of the current frequency slice according to the least - squares optimization data and the Gram matrix; the Gram matrix is pre - established based on the wavenumber spectrum.

[0011] Return to continue performing the step of determining the maximum amplitude until the preset number of iterations is reached, and obtain the wavenumber spectrum after regularization processing of the current frequency slice.

[0012] Based on the wavenumber spectra after regularization processing of all frequency slices, obtain a wavenumber spectrum data volume after regularization processing.

[0013] Perform conversion on the wavenumber spectrum data volume after regularization processing to obtain seismic data in the time - space domain after regularization processing.

[0014] In some alternative embodiments, the conversion of the obtained seismic data in the time - space domain to obtain the wavenumber spectrum data of the seismic data includes:

[0015] Perform Fourier transform on the obtained seismic data in the time - space domain along the time direction to obtain a frequency - space domain data volume.

[0016] Perform non - uniform Fourier transform on the frequency slices in the frequency - space domain data volume along the space direction to obtain the wavenumber spectrum of the frequency slice, and generate a wavenumber spectrum data volume based on the wavenumber spectra of all frequency slices.

[0017] In some alternative embodiments, the performing non - uniform Fourier transform on the frequency slices in the frequency - space domain data volume along the space direction to obtain the wavenumber spectrum of the frequency slice includes:

[0018] Perform non - uniform Fourier transform on the frequency slice f(x m ) in the frequency - space domain data volume along the space direction using the following formula to obtain the wavenumber spectrum F(k n ):

[0019]

[0020] where i is the imaginary unit, M is the number of sample points of the spatial point x, and N is the number of sample points of the wavenumber k;

[0021] The superscripts 1, 2, 3, 4 in the upper right corner respectively represent 4 different spatial dimensions, and the subscript m represents the sample point number of the spatial point x;

[0022] The superscripts 1, 2, 3, 4 in the upper right corner respectively represent 4 different wavenumber dimensions, and the subscript n represents the sample point number of the wavenumber k.

[0023] In some alternative embodiments, determining the maximum amplitude in the amplitude spectrum of the wavenumber spectrum of the current frequency slice according to the dip weight function includes:

[0024] Obtaining the amplitude spectrum of the wavenumber spectrum of the current frequency slice, searching for the maximum amplitude position based on the amplitude spectrum and the dip weight function, and obtaining the wavenumber position corresponding to the maximum amplitude in the amplitude spectrum.

[0025] In some alternative embodiments, searching for the maximum amplitude position based on the amplitude spectrum and the dip weight function includes:

[0026] The following formula is used to search for the maximum amplitude position:

[0027] |w(f,p j )C j (p j )|≥|w(f,n)C j (n)|, 1 ≤ n ≤ N

[0028] where the dip weight function is w(f,k) = (S -1 TS|D(f,k)|) q , and k in the dip weight functions on both sides of the above formula is p j and n respectively; p j and n represent wavenumber positions; C j (p j ) represents the wavenumber component at the wavenumber position p j ; C j (n) represents the wavenumber component at the wavenumber position n;

[0029] In the dip weight function, D(f,k) is the wavenumber spectrum data volume, T is the smoothing operator, S is the amplitude energy operator, representing the amplitude energy sum along the specified ray direction in the wavenumber spectrum data volume, and the exponent q represents the sharpness of the dip weight function: S -1 is the inverse operator of S; where the expression of the amplitude energy operator S is as follows:

[0030] S(θ 1 ,θ 2 ,θ 3 ,θ 4 ) = ∫ r |D(θ 1 ,θ 2 ,θ 3 ,θ 4 ,r)|dr

[0031]

[0032] (θ 1 ,θ2 , θ 3 , θ 4 , r) is the Cartesian coordinate in the frequency - wavenumber domain (f, k 1 , k 2 , k 3 , k 4 ), and the corresponding polar coordinate, r is the radial length, θ j represents the angle of the dimension corresponding to the j - th wavenumber of the wavenumber k; k j represents the dimension corresponding to the j - th wavenumber of the wavenumber k.

[0033] In some alternative embodiments, according to the Gram matrix, perform least - squares optimization on the wavenumber component corresponding to the maximum amplitude, and update the wavenumber spectrum of the current frequency slice according to the least - squares optimization data and the Gram matrix, including:

[0034] Perform least - squares optimization on the wavenumber component corresponding to the maximum amplitude according to the extracted Gram sub - square matrix and the wavenumber spectrum sub - vector extracted from the wavenumber spectrum to obtain least - squares optimization data; the Gram sub - square matrix is obtained by extracting specified elements from the Gram matrix according to the iteration number of the current frequency slice;

[0035] Update the wavenumber spectrum of the current frequency slice according to the least - squares optimization data and the extracted Gram sub - matrix; the Gram sub - matrix is obtained by extracting corresponding columns from the Gram matrix according to the iteration number.

[0036] In some alternative embodiments, perform least - squares optimization on the wavenumber component corresponding to the maximum amplitude according to the extracted Gram sub - square matrix and the wavenumber spectrum sub - vector extracted from the wavenumber spectrum to obtain least - squares optimization data, including:

[0037] Update the preset wavenumber component matrix V0 based on the wavenumber component data at the wavenumber position corresponding to the maximum amplitude;

[0038] Process the wavenumber component matrix V0 using the least - squares method according to the extracted Gram sub - square matrix A and the wavenumber spectrum sub - vector B to obtain the least - squares optimization matrix V; correspondingly, update the wavenumber spectrum of the current frequency slice according to the least - squares optimization data and the extracted Gram sub - matrix, including: calculating the product of the extracted Gram sub - matrix E and the least - squares optimization matrix V, and updating the wavenumber spectrum of the current slice according to the difference between the wavenumber spectrum of the current frequency slice and this product.

[0039] In some alternative embodiments, process the wavenumber component matrix V0 using the least - squares method according to the extracted Gram sub - square matrix A and the wavenumber spectrum sub - vector B to obtain the least - squares optimization matrix V, including:

[0040] The least - squares optimization matrix V is obtained by using the following least - squares equation:

[0041] (A + 0.01I)V=(B + 0.01V0)

[0042] where I is the identity matrix;

[0043] A is a p - order square matrix extracted from the Gram matrix at the position (p l ,p r ), l = 1,…,j, r = 1,…,j, j

[0044] B is a p - dimensional vector extracted from the wavenumber spectrum F(k) at the position p l , l = 1,…,j; j

[0045] j represents the number of iterations;

[0046] Correspondingly, calculating the product of the extracted Gram sub - matrix E and the least - squares optimization matrix V includes:

[0047] Obtain an N - row and p - column Gram sub - matrix E, where the Gram sub - matrix E is a matrix formed by extracting p columns from the Gram matrix at the position p j , l = 1,…,j, and j represents the number of iterations; l j

[0048] Calculate the product of the Gram sub - matrix E and the least - squares optimization matrix V.

[0049] In some alternative embodiments, converting the wavenumber - spectrum data volume after the regularization processing of seismic data to obtain the seismic data in the regularized time - space domain includes:

[0050] Performing non - uniform inverse Fourier transform on the wavenumber spectrum of each frequency slice in the wavenumber - spectrum data volume after the regularization processing of seismic data to obtain the regularized frequency - space domain data volume;

[0051] Performing inverse Fourier transform on the regularized frequency - space domain data volume to obtain the seismic data in the regularized time - space domain.

[0052] An embodiment of the present invention further provides a device for regularizing seismic exploration data, including:

[0053] A first data acquisition module, configured to convert the obtained seismic data in the time - space domain to obtain the wavenumber - spectrum data volume of the seismic data;

[0054] ​​​​A regularization processing module for performing the following regularization processing procedures on the wavenumber spectra of each frequency slice in the wavenumber spectrum data volume respectively:

[0055] Determine the maximum amplitude in the amplitude spectrum of the wavenumber spectrum of the current frequency slice according to the dip weight function; the dip weight function is determined according to the amplitude energy of the wavenumber spectrum data volume in different ray directions;

[0056] Perform least squares optimization on the wavenumber components corresponding to the maximum amplitude according to the Gram matrix, and update the wavenumber spectrum of the current frequency slice according to the least squares optimization data and the Gram matrix; the Gram matrix is established in advance based on the wavenumber spectrum;

[0057] Return to continue executing the step of the maximum amplitude until the preset number of iterations is reached, and obtain the wavenumber spectrum after regularization processing of the current frequency slice;

[0058] Based on the wavenumber spectra after regularization processing of all frequency slices, obtain a wavenumber spectrum data volume after regularization processing;

[0059] A second data acquisition module for converting the wavenumber spectrum data volume after regularization processing to obtain seismic data in the regularized time - space domain.

[0060] An embodiment of the present invention also provides a computer storage medium, in which computer - executable instructions are stored, and when the computer - executable instructions are executed by a processor, the above - mentioned seismic exploration data regularization processing method is implemented.

[0061] An embodiment of the present invention also provides a terminal device, including: a memory, a processor, and a computer program stored on the memory and executable on the processor, and when the processor executes the program, the above - mentioned seismic exploration data regularization processing method is implemented.

[0062] The beneficial effects of the above - mentioned technical solutions provided by the embodiments of the present invention at least include:

[0063] The frequency - wavenumber domain seismic data regularization method provided by the embodiments of the present invention performs regularization processing on the wavenumber spectra of each frequency slice in the wavenumber spectrum data volume of the seismic data transformed into the frequency - wavenumber domain. During the processing, the regularization processing of seismic data with regular missing and large gaps is improved. An inclination weight function is generated by inclination scanning in the frequency - wavenumber domain, and the inclination weight function is used to constrain the energy picking process during each iteration, enhancing the anti - spatial aliasing ability of the regularization algorithm, thereby effectively avoiding the possible generation of reconstruction artifacts and inaccurate data reconstruction results. After selecting the wavenumber component corresponding to the maximum energy each time, all the selected wavenumber components are updated through a least - squares optimization, improving the pseudo - Gibbs phenomenon at the discontinuity points of the event axis. Through data regularization processing, the spatial sampling attributes of seismic exploration data are improved, thereby improving the imaging quality of seismic exploration data.

[0064] Other features and advantages of the present invention will be described in the following specification, and, in part, will be obvious from the specification, or will be understood by implementing the present invention. The objectives and other advantages of the present invention can be achieved and obtained by the structures specifically pointed out in the written specification, claims, and drawings.

[0065] The technical solutions of the present invention will be further described in detail below through the drawings and embodiments. Description of the Drawings

[0066] The drawings are used to provide a further understanding of the present invention, and constitute a part of the specification. They are used to explain the present invention together with the embodiments of the present invention, and do not constitute a limitation to the present invention. In the drawings:

[0067] Figure 1 is a flowchart of the seismic exploration data regularization processing method in the embodiments of the present invention;

[0068] Figure 2 is a flowchart of the regularization processing of a frequency slice in the embodiments of the present invention;

[0069] Figure 3 is an example diagram of seismic data in the embodiments of the present invention;

[0070] Figure 4 is an example diagram of the data regularization result without inclination weight constraint and without least - squares optimization in the embodiments of the present invention;

[0071] Figure 5 is an example diagram of the data regularization result with inclination weight constraint and without least - squares optimization in the embodiments of the present invention;

[0072] Figure 6 is an example diagram of the data regularization result with inclination weight constraint and introduction of least - squares optimization in the embodiments of the present invention;

[0073] Figure 7 This is a schematic structural diagram of the seismic exploration data regularization processing device in the embodiments of the present invention. Detailed implementation manners

[0074] The exemplary embodiments of the present disclosure will be described in more detail below with reference to the accompanying drawings. Although the exemplary embodiments of the present disclosure are shown in the drawings, it should be understood that the present disclosure can be implemented in various forms and should not be limited by the embodiments set forth herein. On the contrary, these embodiments are provided so that the present disclosure can be more thoroughly understood and the scope of the present disclosure can be fully conveyed to those skilled in the art.

[0075] In order to solve the problem that in the prior art, when traditional ALFT algorithms process seismic data with regular missing and large gaps, due to the influence of periodic aliasing, reconstruction artifacts may occur, and pseudo-Gibbs phenomena may occur at the discontinuity points of the isochrones, resulting in poor processing effects of seismic data and affecting the subsequent imaging quality and effects. The embodiments of the present invention provide a frequency-wavenumber domain seismic data regularization method. An inclination weight function is generated by inclination scanning in the frequency-wavenumber domain, and the inclination weight function is used to constrain the energy picking process during each iteration, enhancing the anti-spatial aliasing ability of the regularization algorithm; after each wavenumber component corresponding to the maximum energy is selected, all the selected wavenumber components are updated through a least squares optimization, improving the pseudo-Gibbs phenomenon at the discontinuity points of the isochrones, thereby improving the spatial sampling attributes of seismic data and enhancing the imaging quality of seismic data.

[0076] The embodiments of the present invention provide a frequency-wavenumber domain seismic data regularization method, and its process is as Figure 1 shown, including the following steps:

[0077] Step S11: Based on the acquired time-space domain seismic data, convert to obtain the wavenumber spectrum data volume of the seismic data.

[0078] Step S12: Perform the following regularization processing process on the wavenumber spectrum of each frequency slice in the wavenumber spectrum data volume to obtain the wavenumber spectrum after frequency slice regularization processing.

[0079] Regularize the wavenumber spectra for each frequency slice in the wavenumber spectrum data volume. For each frequency slice, determine the maximum amplitude in the amplitude spectrum of the wavenumber spectrum of the current frequency slice according to the dip weight function. By setting the number of iterations, perform least-squares optimization on the wavenumber component corresponding to the maximum amplitude according to the Gram matrix. Update the wavenumber spectrum of the current frequency slice according to the least-squares optimized data and the Gram matrix. Through multiple iterations with a preset number of iterations, obtain the regularized wavenumber spectrum of the current frequency slice. After processing all the slices in the wavenumber spectrum data volume, obtain the wavenumber spectra of all frequency slices.

[0080] Step S13: Based on the regularized wavenumber spectra of all frequency slices, obtain a regularized wavenumber spectrum data volume.

[0081] Step S14: Transform the regularized wavenumber spectrum data volume to obtain regularized time-space domain seismic data.

[0082] In some alternative embodiments, the above Step S11 transforms the obtained time-space domain seismic data to obtain the wavenumber spectrum data of the seismic data, including:

[0083] 1) Perform a fast Fourier transform on the obtained time-space domain seismic data along the time direction to obtain a frequency-space domain data volume.

[0084] Optionally, perform a fast Fourier transform on the obtained time-space domain seismic data along the time direction to obtain a frequency-space domain data volume. Optionally, other transformation methods can also be used as long as the time-space domain seismic data can be transformed to the frequency-space domain.

[0085] 2) Perform a non-uniform Fourier transform on each frequency slice f(x) in the frequency-space domain data volume along the space direction to obtain the wavenumber spectrum F(k) of the frequency slice, and generate a wavenumber spectrum data volume D(f, k) based on the wavenumber spectra of all frequency slices.

[0086] For the frequency slice f(x m ) in the frequency-space domain data volume, perform a non-uniform Fourier transform along the space direction using the following formula to obtain the wavenumber spectrum F(k n ) of the above frequency slice:

[0087]

[0088] where i is the imaginary unit, M is the number of samples of the spatial point x, and N is the number of samples of the wavenumber k;

[0089] The upper superscripts 1, 2, 3, and 4 respectively represent four different spatial dimensions, and the lower subscript m represents the sample point serial number of the spatial point x; that is, the four dimensions corresponding to the four spatial directions are distinguished by the upper superscripts.

[0090] The upper superscripts 1, 2, 3, and 4 respectively represent four different wavenumber dimensions, and the lower subscript n represents the sample point serial number of the wavenumber k.

[0091] The seismic data can use synthetic 3D seismic records. Taking this as an example below, the four spatial directions adopted are the x and y coordinates of the common midpoint, the offset, and the offset azimuth. For the convenience of explanation, in this embodiment, the number of spatial sample points before and after regularizing the offset and the offset azimuth directions is both 1, the sample point serial number in the x - coordinate direction is cmp, and the sample point serial number in the y - coordinate direction is line. Figure 3 For the display of seismic traces of the data volume, where the missing regular 3D seismic traces are displayed in the cmp - line order, and there are 9 empty traces out of every 10 traces in the line direction of this data volume.

[0092] In the above step S12, the wavenumber spectra of each frequency slice in the wavenumber spectrum data volume are respectively regularized. The process of regularizing a frequency slice is shown in Figure 2 as follows:

[0093] Step S121: According to the dip - weight function, determine the maximum amplitude in the amplitude spectrum of the wavenumber spectrum of the current frequency slice; among them, the dip - weight function is determined according to the amplitude energy of the wavenumber spectrum data volume in different ray directions.

[0094] In this step, obtain the amplitude spectrum of the wavenumber spectrum of the current frequency slice, search for the position of the maximum amplitude based on the amplitude spectrum and the dip - weight function, and obtain the wavenumber position corresponding to the maximum amplitude in the amplitude spectrum. This step uses the dip - weight function to constrain the energy picking process during each iteration to select the wavenumber component corresponding to the maximum energy, enhancing the anti - spatial aliasing ability of the regularization algorithm.

[0095] Step S122: According to the Gram matrix, perform least - squares optimization on the wavenumber component corresponding to the maximum amplitude, and update the wavenumber spectrum of the current frequency slice according to the least - squares optimization data and the Gram matrix; among them, the Gram matrix is established in advance based on the wavenumber spectrum.

[0096] In this step, according to the extracted Gram submatrix and the wavenumber spectrum subvector extracted from the wavenumber spectrum, least-squares optimization is performed on the wavenumber component corresponding to the maximum amplitude to obtain least-squares optimization data; the Gram submatrix is obtained by extracting specified elements from the Gram matrix according to the iteration number of the current frequency slice; according to the least-squares optimization data and the extracted Gram submatrix, the wavenumber spectrum of the current frequency slice is updated; the Gram submatrix is obtained by extracting the corresponding columns from the Gram matrix according to the iteration number. In this step, by performing least-squares optimization to update all the selected wavenumber components after each selection of the wavenumber component corresponding to the maximum energy, the pseudo Gibbs phenomenon at the discontinuity points of the isophase axis is improved.

[0097] Step S123: Determine whether the preset number of iterations is reached. If so, execute step S124; otherwise, return to step S121 to continue the step of determining the maximum amplitude.

[0098] Step S124: Obtain the regularized wavenumber spectrum of the current frequency slice.

[0099] In some alternative embodiments, before regularizing the frequency slice, a Gram matrix and a dip weight function can be established in advance.

[0100] The Gram matrix is established based on the wavenumber spectrum. The Gram matrix can be an N-order square matrix G. The N-order square matrix G = (g1,..., g N ), where the nth (n ≤ N) column vector g n of G is as follows:

[0101]

[0102] The dip weight function is determined according to the amplitude energy of the wavenumber spectrum data volume in different ray directions. Specifically, for the wavenumber spectrum data volume D(f, k), the amplitude energy sum is obtained along different ray directions passing through the origin, and then the dip weight function w(f, k) is obtained according to the amplitude energy sum:

[0103] w(f, k) = (S -1 TS|D(f, k)|) q

[0104] where D(f, k) is the wavenumber spectrum data volume, T is a smoothing operator, S is an amplitude energy operator, representing the amplitude energy sum along the specified ray direction in the wavenumber spectrum data volume, and the exponent q represents the sharpness of the dip weight function: S -1 is the inverse operator of S; among them, the expression of the amplitude energy operator S is as follows:

[0105] S(θ 1 , θ 2 , θ 3 , θ4 ) = ∫ r |D(θ 1 , θ 2 , θ 3 , θ 4 , r)| dr

[0106]

[0107] where (θ 1 , θ 2 , θ 3 , θ 4 , r) is the polar coordinate corresponding to the Cartesian coordinate (f, k 1 , k 2 , k 3 , k 4 ) in the frequency - wavenumber domain, r is the radial length, and θ j represents the angle of the j - th wavenumber corresponding to the wavenumber k in the corresponding dimension; k j represents the j - th wavenumber corresponding to the wavenumber k in the corresponding dimension. The operator S represents obtaining the sum of ray amplitude energies along the angular ray direction in the frequency - wavenumber domain; the operator T represents processes such as smoothing, and the exponent q is used to control the sharpness of the dip weight function w. S -1 represents the inverse process of the transformation S, which means repeatedly placing the energy values in the ray direction corresponding to (θ 1 , θ 2 , θ 3 , θ 4 ) into the corresponding (k 1 , k 2 , k 3 , k 4 ) positions at different frequencies f.

[0108] In the above step S12, for each frequency slice, the regularization process of steps S121 - 124 is performed. For a frequency slice, the number of iterations can be set, and the regularization process is achieved through multiple iterations. For example: Let the number of iterations j = 1, copy the wavenumber spectrum F(k) of the frequency slice in the data volume D(f, k) as C(k), and perform the regularization process on it as the current frequency slice.

[0109] Optionally, in the above step S121: When determining the maximum amplitude in the amplitude spectrum of the wavenumber spectrum of the current frequency slice according to the dip weight function, the amplitude spectrum of the wavenumber spectrum of the current frequency slice can be obtained, and the maximum amplitude position is searched based on the amplitude spectrum and the dip weight function to obtain the wavenumber position corresponding to the maximum amplitude in the amplitude spectrum. The following formula can be used for the maximum amplitude position search:

[0110] |w(f, p j )C j (p j)|≥|w(f,n)C j (n), 1 ≤ n ≤ N

[0111] where the dip angle weight function is w(f,k) = (S -1 TS|D(f,k)|) q , where k in the dip angle weight functions on both the left and right sides of the above formula is p j and n respectively; p j and n represent the wavenumber positions; C j (p j ) represents the wavenumber component at the wavenumber position p j ; C j (n) represents the wavenumber component at the wavenumber position n.

[0112] Continuing with the above example, for the amplitude spectrum of the current frequency slice C(k), based on the amplitude spectrum and the dip angle weight function, the following formula is used to search for the maximum amplitude position, obtaining the wavenumber position p j corresponding to the maximum amplitude |C j (p j ). The wavenumber component C j at the wavenumber position p j (p j ) is placed at the corresponding position in the vector V0.

[0113] Optionally, in the above step S122, according to the Gram matrix, least squares optimization is performed on the wavenumber component corresponding to the maximum amplitude. Based on the least squares optimization data and the Gram matrix, the wavenumber spectrum of the current frequency slice is updated, including the following implementation process:

[0114] 1) Extract the specified elements from the Gram matrix according to the iteration number of the current frequency slice to obtain the Gram sub-square matrix, and extract the wavenumber spectrum sub-vector from the wavenumber spectrum

[0115] According to the iteration number j of the current frequency slice, extract the specified elements from the Gram matrix G to obtain the sub-square matrix A and extract the sub-vector B from the wavenumber spectrum, where:

[0116] A is a p l -order square matrix extracted from the Gram matrix according to the (p r , p j ) position, l = 1,..., j, r = 1,..., j; j represents the iteration number;

[0117] B is a p l -dimensional vector extracted from the wavenumber spectrum F(k) according to the p j position, l = 1,..., j.

[0118] 2) According to the extracted Gram submatrix and the wavenumber spectrum subvector extracted from the wavenumber spectrum, perform least squares optimization on the wavenumber component corresponding to the maximum amplitude to obtain least squares optimization data.

[0119] Update the preset wavenumber component matrix V0 based on the wavenumber component data at the wavenumber position corresponding to the maximum amplitude;

[0120] According to the extracted Gram submatrix A and the wavenumber spectrum subvector B, process the wavenumber component matrix V0 using the least squares method to obtain the least squares optimization matrix V; the following least squares equation can be used to obtain the least squares optimization matrix V:

[0121] (A + 0.01I)V = (B + 0.01V0)

[0122] where I is the identity matrix.

[0123] After obtaining the least squares optimization matrix V, the matrix V can be copied as V0, that is, update V0 with the matrix V.

[0124] 3) Extract the corresponding columns from the Gram matrix according to the number of iterations to obtain the Gram submatrix.

[0125] Obtain an N-row p j column Gram submatrix E. The Gram submatrix E is a matrix formed by extracting p l columns from the Gram matrix according to the p j positions, l = 1,..., j, where j represents the number of iterations; specifically, according to the number of iterations j of the current frequency slice, extract an N-row p l column matrix from the Gram matrix to obtain the matrix E, and N is the number of wavenumber samples.

[0126] 4) Update the wavenumber spectrum of the current frequency slice according to the least squares optimization data and the extracted Gram submatrix.

[0127] Calculate the product of the extracted Gram submatrix E and the least squares optimization matrix V, and update the wavenumber spectrum of the current slice according to the difference between the wavenumber spectrum of the current frequency slice and this product. That is, calculate the product of the N-row p j column matrix E and the p j dimensional vector V, and update the wavenumber spectrum C(k) with the result of subtracting the product result from the wavenumber spectrum F(k) of the current frequency slice.

[0128] If it is determined in step S123 that the preset number of iterations has not been reached, let the number of iterations j = j + 1, and repeat steps S121 and S122 until the preset number of iterations is reached, that is, all wavenumber components of this frequency slice are obtained, and then the regularized wavenumber spectrum of the current frequency slice can be obtained. Then, store the obtained least squares optimization matrix V in the corresponding position of the output matrix H.

[0129] The output matrix H is a matrix storing the regularized wavenumber spectra of all frequency slices. After the wavenumber spectrum of each frequency slice is regularized, the matrix V is stored in the corresponding position. After all frequency slices are stored, the regularized wavenumber spectrum data volume can be obtained from this matrix H.

[0130] In some optional embodiments, the above step S14 transforms the regularized wavenumber spectrum data volume to obtain regularized time-space domain seismic data, including: performing a non-uniform inverse Fourier transform on the wavenumber spectrum of each frequency slice in the regularized wavenumber spectrum data volume of the seismic data to obtain a regularized frequency-space domain data volume; performing an inverse Fourier transform on the regularized frequency-space domain data volume to obtain regularized time-space domain seismic data. Continuing with the above example, after obtaining all wavenumber components of all frequency slices, the output matrix H is then inverse-transformed back to the time-space domain as the data regularization result.

[0131] The frequency-wavenumber domain seismic data regularization method provided by the embodiments of the present invention separately regularizes the wavenumber spectra of each frequency slice in the wavenumber spectrum data volume of the seismic data transformed to the frequency-wavenumber domain. During the processing, the regularization of seismic data with regular missing and large gaps is improved. An inclination weight function is generated by inclination scanning in the frequency-wavenumber domain, and the inclination weight function is used to constrain the energy picking process during each iteration, enhancing the anti-spatial aliasing ability of the regularization algorithm; after each wavenumber component corresponding to the maximum energy is selected, all selected wavenumber components are updated through a least squares optimization, improving the pseudo-Gibbs phenomenon at the discontinuity points of the isophase axis.

[0132] The above method adds an isophase axis inclination constraint to enhance the anti-spatial aliasing ability of the regularization algorithm, and adds a least squares optimization process to improve the pseudo-Gibbs phenomenon at the discontinuity points of the isophase axis, thereby improving the spatial sampling attributes of the seismic data and enhancing the imaging quality of the seismic data.

[0133] Using the above seismic exploration data regularization method to process synthetic 3D seismic record data can obtain good results. The following analyzes and explains according to the graphical results in different cases, including the comparison of the data regularization results in different cases such as with or without inclination constraint and with or without least squares optimization.

[0134] Figure 4It is the regularization result without dip weight constraint and without least - squares optimization. Without dip weight constraint means that each element value of the dip weight function w in step S121 takes 1, that is, no dip weight constraint is performed; without least - squares optimization means that the least - squares optimization in step S122 is not executed, and V0 is directly copied to V.

[0135] Figure 5 It is the regularization result with dip constraint and without least - squares optimization. That is, the least - squares optimization in step S122 is not executed, and V0 is directly copied to V.

[0136] Figure 6 It is the regularization result with dip constraint and after introducing least - squares optimization, that is, the result after performing the regularization process using the processing procedures in steps S121 - 124.

[0137] For Figure 4 and Figure 5 It can be seen that for seismic data with missing rules, due to the interference of periodic aliasing energy, without dip constraint, inaccurate reconstruction results will occur during the regularization process due to incorrect energy search. After adding the dip constraint of the event axis, the regularization effect is significantly improved, but pseudo - Gibbs phenomena will occur at the discontinuity points of the event axis.

[0138] Comparing Figure 5 and Figure 6 It can be seen that after introducing least - squares optimization, the pseudo - Gibbs phenomena at the discontinuity points of the event axis are suppressed, and the reconstruction effect is further improved.

[0139] The above practice proves that a frequency - wavenumber domain seismic data regularization method provided by an embodiment of the present invention can improve the spatial sampling attributes of seismic data and the imaging quality of seismic data after processing the data with missing rules. The regularization effect is significantly improved, and the reconstruction effect is further improved.

[0140] Based on the same inventive concept, an embodiment of the present invention further provides a seismic exploration data regularization processing device. This device can be set in a computer device with data - processing capabilities. The structure of this device is as Figure 7 shown, including: a first data acquisition module 101, a regularization processing module 102, and a second data acquisition module 103.

[0141] The first data acquisition module 101 is used to convert the acquired time - space domain seismic data into a wavenumber spectrum data volume of the seismic data;

[0142] The regularization processing module 102 is used to perform the following regularization processing procedures on the wavenumber spectrum of each frequency slice in the wavenumber spectrum data volume respectively:

[0143] Determine the maximum amplitude in the amplitude spectrum of the wavenumber spectrum of the current frequency slice according to the dip angle weight function; the dip angle weight function is determined according to the amplitude energy of the wavenumber spectrum data volume in different ray directions;

[0144] Perform least squares optimization on the wavenumber component corresponding to the maximum amplitude according to the Gram matrix, and update the wavenumber spectrum of the current frequency slice according to the least squares optimization data and the Gram matrix; the Gram matrix is pre-established based on the wavenumber spectrum;

[0145] Return to continue executing the step of the maximum amplitude until the preset number of iterations is reached, and obtain the regularized wavenumber spectrum of the current frequency slice;

[0146] Based on the regularized wavenumber spectra of all frequency slices, obtain the regularized wavenumber spectrum data volume;

[0147] The second data acquisition module 103 is used to convert the regularized wavenumber spectrum data volume to obtain the regularized time-space domain seismic data.

[0148] An embodiment of the present invention further provides a computer storage medium, in which computer-executable instructions are stored, and when the computer-executable instructions are executed by a processor, the above-mentioned seismic exploration data regularization processing method is implemented.

[0149] An embodiment of the present invention further provides a terminal device, including: a memory, a processor, and a computer program stored on the memory and executable on the processor, and when the processor executes the program, the above-mentioned seismic exploration data regularization processing method is implemented.

[0150] Regarding the seismic exploration data regularization processing device in the above embodiments, the specific manners in which each module performs operations have been described in detail in the embodiments related to the method, and will not be elaborated here.

[0151] Unless otherwise specifically stated, terms such as processing, computing, operation, determination, display, etc. may refer to the actions and / or processes of one or more processing or computing systems, or similar devices, and the actions and / or processes will represent data operations and conversions of physical (such as electronic) quantities in the registers or memories of the processing system into other data of physical quantities similarly represented in the memories, registers, or other such information storage, transmission, or display devices of the processing system. Information and signals can be represented using any of a variety of different technologies and methods. For example, the data, instructions, commands, information, signals, bits, symbols, and chips mentioned throughout the above description can be represented by voltages, currents, electromagnetic waves, magnetic fields or particles, optical fields or particles, or any combination thereof.

[0152] It should be understood that the specific order or hierarchy of steps in the disclosed process is an example of an exemplary method. Based on design preferences, it should be understood that the specific order or hierarchy of steps in the process can be rearranged without departing from the scope of the present disclosure. The appended method claims present the elements of the various steps in an exemplary order and are not intended to be limited to the specific order or hierarchy recited.

[0153] In the above detailed description, various features are combined in a single embodiment to simplify the present disclosure. This method of disclosure should not be interpreted as reflecting an intention that the embodiments of the claimed subject matter require more features than are expressly recited in each claim. Rather, as reflected in the appended claims, the invention lies in less than all of the features of a single disclosed embodiment. Accordingly, the appended claims are hereby expressly incorporated into the detailed description, with each claim standing on its own as a separate preferred embodiment of the invention.

[0154] Those skilled in the art should also understand that the various illustrative logical blocks, modules, circuits, and algorithm steps described in connection with the embodiments herein can be implemented as electronic hardware, computer software, or combinations thereof. To clearly illustrate the interchangeability of hardware and software, the various illustrative components, blocks, modules, circuits, and steps have been generally described in terms of their functionality. Whether such functionality is implemented as hardware or software depends upon the particular application and the design constraints imposed on the overall system. Skilled artisans may implement the described functionality in a flexible manner for each particular application, but such implementation decisions should not be interpreted as departing from the scope of the present disclosure.

[0155] The steps of a method or algorithm described in connection with the embodiments herein may be embodied directly in hardware, in a software module executed by a processor, or in a combination thereof. The software module may be located in RAM memory, flash memory, ROM memory, EPROM memory, EEPROM memory, registers, hard disk, a removable disk, a CD-ROM, or any other form of storage medium well known in the art. An exemplary storage medium is coupled to the processor such that the processor can read information from, and write information to, the storage medium. Of course, the storage medium may also be integral to the processor. The processor and the storage medium may be located in an ASIC. The ASIC may be located in a user terminal. Of course, the processor and the storage medium may also exist as discrete components in a user terminal.

[0156] For software implementation, the techniques described in this application can be implemented by modules (e.g., procedures, functions, etc.) that perform the functions described in this application. These software codes can be stored in a memory unit and executed by a processor. The memory unit can be implemented within the processor or outside the processor. In the latter case, it is communicatively coupled to the processor via various means, which are well known in the art.

[0157] The above description includes examples of one or more embodiments. Of course, it is not possible to describe all possible combinations of components or methods for the purpose of describing the above embodiments, but those of ordinary skill in the art should recognize that the various embodiments can be further combined and arranged. Therefore, the embodiments described herein are intended to cover all such changes, modifications, and variations that fall within the scope of the appended claims. In addition, with respect to the term "comprising" used in the specification or claims, the word is construed in a manner similar to the term "including," as "including" is interpreted when used as a transitional word in a claim. In addition, any use of the term "or" in the specification or claims of a patent is to mean "non-exclusive or."

Claims

1. A method for regularizing seismic exploration data, characterized in that Including: Based on the acquired time - space domain seismic data, converting to obtain a wavenumber spectrum data volume of the seismic data; Performing the following regularization processing procedure on the wavenumber spectrum of each frequency slice in the wavenumber spectrum data volume: Determining the maximum amplitude in the amplitude spectrum of the wavenumber spectrum of the current frequency slice according to the dip weight function; the dip weight function is determined according to the amplitude energy of the wavenumber spectrum data volume in different ray directions; Performing least - squares optimization on the wavenumber components corresponding to the maximum amplitude according to the Gram matrix, and updating the wavenumber spectrum of the current frequency slice according to the least - squares optimization data and the Gram matrix; The Gram matrix is pre - established based on the wavenumber spectrum; Returning to continue performing the step of determining the maximum amplitude until a preset number of iterations is reached, to obtain the regularized wavenumber spectrum of the current frequency slice; Based on the regularized wavenumber spectra of all frequency slices, obtaining a regularized wavenumber spectrum data volume; Performing conversion on the regularized wavenumber spectrum data volume to obtain regularized time - space domain seismic data.

2. The method according to claim 1, wherein The converting to obtain the wavenumber spectrum data of the seismic data based on the acquired time - space domain seismic data includes: Performing Fourier transform on the acquired time - space domain seismic data along the time direction to obtain a frequency - space domain data volume; Performing non - uniform Fourier transform on the frequency slices in the frequency - space domain data volume along the space direction to obtain the wavenumber spectra of the frequency slices, and generating a wavenumber spectrum data volume based on the wavenumber spectra of all frequency slices.

3. The method according to claim 2, wherein The performing non - uniform Fourier transform on the frequency slices in the frequency - space domain data volume along the space direction to obtain the wavenumber spectra of the frequency slices includes: For the frequency slice f(x m ) in the frequency - space domain data volume, perform non - uniform Fourier transform along the spatial direction using the following formula to obtain the wavenumber spectrum F(k n ) of the frequency slice: Where i is the imaginary unit, M is the number of samples of the spatial point x, and N is the number of samples of the wavenumber k; The superscripts 1, 2, 3, 4 in the upper right corner respectively represent four different spatial dimensions, and the subscript m represents the sample point number of the spatial point x; The superscripts 1, 2, 3, 4 in the upper right corner respectively represent four different wavenumber dimensions, and the subscript n represents the sample number of the wavenumber k.

4. The method according to claim 1, wherein Determining the maximum amplitude in the amplitude spectrum of the wavenumber spectrum of the current frequency slice according to the dip weight function includes: Obtaining the amplitude spectrum of the wavenumber spectrum of the current frequency slice, searching for the maximum amplitude position based on the amplitude spectrum and the dip weight function, and obtaining the wavenumber position corresponding to the maximum amplitude in the amplitude spectrum.

5. The method according to claim 4, characterized in that, Searching for the maximum amplitude position based on the amplitude spectrum and the dip weight function includes: Performing maximum amplitude position search using the following formula: |w(f,p j )C j (p j )|≥|w(f,n)C j (n)|, 1 ≤ n ≤ N Among them, the dip angle weight function is w(f,k) = (S -1 TS|D(f,k)|) q . In the above formula, k in the dip angle weight functions on both the left and right sides is p j and n respectively; p j and n represent the wavenumber positions; C j (p j ) represents the wavenumber component at the wavenumber position p j ; C j (n) represents the wavenumber component at the wavenumber position n In the dip angle weighting function, D(f,k) is the wavenumber spectrum data volume, T is the smoothing operator, S is the amplitude energy operator, representing the sum of the amplitude energy along the specified ray direction in the wavenumber spectrum data volume, and the exponent q represents the sharpness of the dip angle weighting function: S -1 is the inverse operator of S; among them, the expression of the amplitude energy operator S is as follows: S(θ 1 ,θ 2 ,θ 3 ,θ 4 ) = ∫ r |D(θ 1 ,θ 2 ,θ 3 ,θ 4 , r)| dr (θ 1 , θ 2 , θ 3 , θ 4 , r) is the polar coordinate corresponding to the Cartesian coordinate (f, k 1 , k 2 , k 3 , k 4 ) in the frequency-wavenumber domain, r is the radial length, and θ j represents the angle of the dimension corresponding to the j-th wavenumber of the wavenumber k; k j represents the dimension corresponding to the j-th wavenumber of the wavenumber k.

6. The method according to claim 1, characterized in that, Performing least - squares optimization on the wavenumber components corresponding to the maximum amplitude according to the Gram matrix, and updating the wavenumber spectrum of the current frequency slice according to the least - squares optimization data and the Gram matrix includes: Performing least - squares optimization on the wavenumber components corresponding to the maximum amplitude according to the extracted Gram sub - matrix and the wavenumber spectrum sub - vector extracted from the wavenumber spectrum to obtain least - squares optimization data; the Gram sub - matrix is obtained by extracting specified elements from the Gram matrix according to the number of iterations of the current frequency slice; Updating the wavenumber spectrum of the current frequency slice according to the least - squares optimization data and the extracted Gram sub - matrix; the Gram sub - matrix is obtained by extracting the corresponding columns from the Gram matrix according to the number of iterations.

7. The method according to claim 6, characterized in that, Performing least - squares optimization on the wavenumber components corresponding to the maximum amplitude according to the extracted Gram sub - matrix and the wavenumber spectrum sub - vector extracted from the wavenumber spectrum to obtain least - squares optimization data includes: Update a preset wavenumber component matrix V0 based on the wavenumber component data at the wavenumber position corresponding to the maximum amplitude; Based on the extracted Gram sub-square matrix A and the wavenumber spectrum sub-vector B, process the wavenumber component matrix V0 using the least squares method to obtain a least squares optimization matrix V; correspondingly, update the wavenumber spectrum of the current frequency slice according to the least squares optimization data and the extracted Gram sub-matrix, including: calculating the product of the extracted Gram sub-matrix E and the least squares optimization matrix V, and updating the wavenumber spectrum of the current slice according to the difference between the wavenumber spectrum of the current frequency slice and this product.

8. The method according to claim 7, wherein Based on the extracted Gram sub-square matrix A and the wavenumber spectrum sub-vector B, process the wavenumber component matrix V0 using the least squares method to obtain a least squares optimization matrix V, including: Use the following least squares equation to obtain the least squares optimization matrix V: (A + 0.01I)V = (B + 0.01V0) where I is the identity matrix; A is extracted from the Gram matrix according to the (p l , p r ) position as a p j -order square matrix, where l = 1, …, j, r = 1, …, j B is extracted from the wavenumber spectrum F(k) at position p l to obtain a p j -dimensional vector, where l = 1, …, j; j represents the number of iterations; Correspondingly, calculate the product of the extracted Gram sub-matrix E and the least squares optimization matrix V, including: Obtain N rows of p j Column Gram sub-matrix E, where the Gram sub-matrix E is l Extracted from the Gram matrix according to the p j Columns to form a matrix, l = 1, …, j, where j represents the number of iterations; Calculate the product of the Gram sub-matrix E and the least squares optimization matrix V.

9. The method according to any one of claims 1-8, characterized in that, Convert the wavenumber spectrum data volume after regularization processing of seismic data to obtain seismic data in the regularized time-space domain, including: Perform non-uniform inverse Fourier transform on the wavenumber spectrum of each frequency slice in the wavenumber spectrum data volume after regularization processing of seismic data to obtain the regularized frequency-space domain data volume; Perform inverse Fourier transform on the regularized frequency-space domain data volume to obtain seismic data in the regularized time-space domain.

10. An apparatus for regularizing seismic exploration data, characterized in that, Including: A first data acquisition module, configured to convert the acquired seismic data in the time-space domain to obtain the wavenumber spectrum data volume of the seismic data; A regularization processing module, configured to perform the following regularization processing process on the wavenumber spectrum of each frequency slice in the wavenumber spectrum data volume respectively: Determine the maximum amplitude in the amplitude spectrum of the wavenumber spectrum of the current frequency slice according to the dip weight function; the dip weight function is determined according to the amplitude energy of the wavenumber spectrum data volume in different ray directions; Perform least squares optimization on the wavenumber components corresponding to the maximum amplitude according to the Gram matrix, and update the wavenumber spectrum of the current frequency slice according to the least squares optimization data and the Gram matrix; The Gram matrix is pre-established based on the wavenumber spectrum; Return to continue executing the step of the maximum amplitude until the preset number of iterations is reached to obtain the wavenumber spectrum of the current frequency slice after regularization processing; Based on the wavenumber spectra of all frequency slices after regularization processing, obtain the wavenumber spectrum data volume after regularization processing; A second data acquisition module, configured to convert the wavenumber spectrum data volume after regularization processing to obtain seismic data in the regularized time-space domain.

11. A computer storage medium, characterized in that, The computer storage medium stores computer-executable instructions, and when the computer-executable instructions are executed by a processor, the seismic exploration data regularization processing method according to any one of claims 1-9 is implemented.

12. A terminal device, characterized in that, Including: A memory, a processor, and a computer program stored on the memory and executable on the processor, wherein when the processor executes the program, the method for regularizing seismic exploration data according to any one of claims 1-9 is implemented.

Images