Seismic data reconstruction method and computer readable storage medium

Abstract

The invention provides a seismic data reconstruction method, which comprises the following steps of: executing Radon transformation on seismic data to convert missing seismic data from a time domain to a Radon transformation domain; combining the seismic data converted to the Radon transform domain with an orthogonal polynomial representing seismic AVO characteristics to form high-order Radon transform; performing inverse transformation on the high-order Radon transformation to obtain missing part data of the seismic data; filling the seismic data with the missing part data to obtain reconstructed new seismic data; according to the method, high-order Radon transformation and repeated iteration are carried out on the seismic data, high-fidelity reproduction of the missing seismic data is realized, and reconstruction of the missing and irregular seismic data plays an important role in the seismic data processing process.

Description

Technical Field

[0001] The present invention belongs to the technical field of seismic data processing, and particularly relates to a seismic data reconstruction method and a computer-readable storage medium. Background Art

[0002] In the field of petroleum geophysical exploration, especially in seismic exploration technology, challenges in seismic data acquisition are often faced. When implementing a field seismic data acquisition project, problems such as poor working conditions of geophones and cables and poor coupling effect with the surface often lead to random missing of the acquired seismic data or the existence of strong-energy noise. In addition, obstacles or no-go areas such as rivers, lakes, roads, bridges, and important infrastructure frequently appearing in the field work area result in segmental missing of seismic data acquisition.

[0003] Existing mainstream seismic data processing technologies have relatively high requirements for the regularity and integrity of seismic records. Incomplete seismic data will seriously affect the processing effect. For example, in the application of stacking technology, the missing of data leads to insufficient stacking times for some common midpoint bins, affecting stacking energy and phase. In the wave equation migration technology, the missing and irregularity of data will cause spatial aliasing phenomena and reduce imaging accuracy.

[0004] Therefore, reconstructing missing and irregular seismic data plays an important role in the process of seismic data processing. The present invention proposes a seismic data reconstruction method, aiming to improve the integrity of seismic data and provide a better data basis for seismic data processing. This method is based on Radon transform and combines amplitude equalization processing, which can effectively reconstruct the missing part of seismic data and improve data quality.

[0005] The present invention performs hyperbolic Radon transform on seismic data with some missing seismic traces and performs amplitude equalization processing to obtain seismic data that conforms to its proper form in terms of morphology and amplitude. Through this method, while ensuring data quality, it can overcome the data acquisition problems brought by the complexity of the on-site environment and equipment limitations. Summary of the Invention

[0006] The purpose of the present application is to provide a seismic data reconstruction method and a computer-readable storage medium to solve the above problems.

[0007] The purpose of the present application is achieved by adopting the following technical solutions:

[0008] In the first aspect, the present application provides a seismic data reconstruction method, including the following steps: performing Radon transform on seismic data to transform the missing seismic data from the time domain to the Radon transform domain;

[0009] Combining the seismic data transformed into the Radon transform domain with orthogonal polynomials characterizing seismic AVO characteristics to form a high-order Radon transform;

[0010] Performing an inverse transform on the high-order Radon transform to obtain the missing data part of the seismic data;

[0011] Filling the missing data part into the seismic data to obtain a newly reconstructed seismic data.

[0012] The beneficial effects of this technical solution are as follows: By performing the Radon transform, this technical solution can effectively transform the missing seismic data from the time domain to the Radon transform domain, thereby making the data reconstruction more accurate and effective; Combining the seismic data transformed into the Radon transform domain with orthogonal polynomials characterizing seismic AVO characteristics to form a high-order Radon transform can better characterize the relationship between amplitude variation and offset variation in seismic data, not only improving the processing accuracy of seismic data but also enhancing the ability to identify the properties of underground rocks and fluids; By performing an inverse transform on the high-order Radon transform, the missing part in the seismic data can be accurately reconstructed, and a seismic waveform close to the real one can be restored, thus providing a more reliable data basis for subsequent interpretation and analysis; By filling the reconstructed missing data part back into the original data, a complete new seismic data set is generated.

[0013] Further, after the step of filling the missing data part into the seismic data to obtain a newly reconstructed seismic data, it further includes using the new seismic data as the missing seismic data for continuous iteration.

[0014] The beneficial effects of this technical solution are as follows: By continuously iterating and processing the newly reconstructed seismic data, the seismic data can be further optimized and refined. Each iteration may improve the data quality and make it closer to the real geological situation; The iteration process enables the data reconstructed each time to be improved based on the previous one, thereby gradually reducing the reconstruction error. For data collected in a complex geological environment, a single reconstruction may not be able to fully recover the missing information, and the iterative processing can better adapt to this complexity and gradually improve the accuracy and reliability of data reconstruction; The iteration process provides a dynamic adjustment and feedback mechanism, allowing the processing strategy to be adjusted in a timely manner according to the quality of the reconstruction results.

[0015] Further, the expression of the Radon transform is:

[0016]

[0017]

[0018] Its discrete form is:

[0019]

[0020]

[0021] Among them, τ is the intercept time; p is the reciprocal of the velocity, i.e., slowness; x is the shot-receiver distance or offset; m(τ, p) is the seismic data after being transformed into the Radon domain; d(t, x) is the two-dimensional common-shot gather seismic data, and i and k are ordinal numbers.

[0022] The beneficial effects of this technical solution are as follows: The discrete form of the Radon transform allows the effective conversion of continuous seismic data into a discrete representation in the Radon domain, making the data easier to operate and understand during processing and interpretation; through the precise Radon transform, it is possible to better identify and process noise and anomalies in the seismic data, thereby improving the quality and usability of the finally reconstructed data.

[0023] Furthermore, the step of combining the seismic data transformed into the Radon transform domain with the orthogonal polynomials characterizing the seismic AVO characteristics to form a high-order Radon transform includes calculating the conjugate solution of the high-order Radon transform and calculating the high-order Radon transform using the matching pursuit algorithm;

[0024] Among them, the specific steps for calculating the conjugate solution of the high-order Radon transform are as follows:

[0025] {q j (x), j = 0, 1, …, N} is a set of orthonormal polynomials determined by N offset coordinates x

[0026] q j (x) is an algebraic polynomial of order j;

[0027] The amplitude change of the seismic data d(t, x) at a certain moment can be fitted with orthonormal polynomials, that is:

[0028]

[0029] Among them, c j (t) is the decomposition coefficient of the orthonormal polynomial of order j for the amplitude change with the offset at time t, which is called the orthonormal polynomial coefficient spectrum. Combining Equation (2) and Equation (5) can obtain:

[0030]

[0031] Since {q j (x), j = 0, 1, …, N} is a set of orthonormal polynomials, multiplying both sides of Equation (6) by q j (x) can obtain the conjugate solution of the high-order Radon transform:

[0032]

[0033] The beneficial effects of this technical solution are as follows: The orthogonal polynomial combined with the seismic AVO characteristics can more accurately reflect the reflection characteristics of seismic waves at different rock interfaces. The use of this high-order Radon transform helps to improve the interpretation ability of seismic data, especially when identifying geological structures such as oil and gas reservoirs; by combining the seismic AVO characteristics, this method optimizes the data processing flow and improves the processing efficiency, which is particularly important for processing large-scale seismic data sets, saving time and reducing the consumption of computing resources; the high-order Radon transform can better adapt to complex geological environments, such as heterogeneous rock formations and complex stratigraphic structures, making the method effective in diverse exploration environments.

[0034] Furthermore, the specific steps for calculating the high-order Radon transform using the matching pursuit algorithm are as follows:

[0035] Represent the matrix of the Radon transform as:

[0036] M = DL (8)

[0037] where L is the integral path function;

[0038] S1. Let the seismic data D be the initial residual D resi , and the Radon parameter M be 0;

[0039] S2. Calculate the conjugate solution M of the Radon transform of the residual data using Equation (7) i ;

[0040] S3. Calculate the Radon energy distribution using Equation (9), and determine the subspace S through the local maximum of the energy i ;

[0041]

[0042] S4. The energy represents the energy distribution of the in-phase axis along different curvature directions at different times;

[0043] S5. Minimize the function J(m) = ‖D resi - LM i ‖ 2 using the conjugate gradient method to obtain the least squares solution M of the Radon transform in the above subspace S i ; i ;

[0044] S6. Update the parameter M = M + M i , and the residual data D resi = D resi - LM i ;

[0045] S7. Determine whether the predefined number of iterations is reached. If the predefined number of iterations is reached, output all subspaces S i Construct a high-order Radon transform. If the predefined number of iterations is not reached, return to S1 for continued iteration.

[0046] The beneficial effects of this technical solution are as follows: By calculating the conjugate solution of the Radon transform of the residual data and updating the parameters, this method can effectively eliminate noise and enhance the data signal. The iterative process allows continuous optimization of data processing, with each step adjusted based on the results of the previous step, contributing to gradually approaching the optimal data reconstruction result; By calculating the Radon energy distribution and determining the local maximum value, this method can better adapt to complex geological environments, especially for seismic data processing of heterogeneous and multi-layer complex strata; Using the conjugate gradient method to minimize the function and the selection of subspaces can improve the calculation efficiency and reduce the required computing resources.

[0047] Furthermore, the specific calculation method for the value range of p is as follows:

[0048] Assume the sampling rate of p is Δp, then

[0049] Δp = p k - p k-1 (10)

[0050] where p in the formula k and p k-1 are any two adjacent slowness values;

[0051] v max is the maximum effective reflection velocity in the seismic data, then Δp is expressed as:

[0052]

[0053] where Δv is the value of v k - v k-1 and

[0054]

[0055] After determining the value of Δp, according to p = 1 / v, give the initial value p0 and the maximum value p of p max The range between the initial value and the maximum value is the value range of p.

[0056] The beneficial effects of this technical solution are as follows: By calculating the value range of the slowness p, the parameter space for Radon transform processing can be accurately determined; This technical solution allows adjustment of the slowness range according to different seismic data characteristics to make it more suitable for specific geological conditions and exploration targets.

[0057] Further, it also includes restoring the amplitude through an amplitude equalization coefficient during the iteration process. The specific steps are as follows:

[0058] S8. Calculate the average amplitude of the non-missing traces near the missing seismic trace;

[0059] S9. Calculate the average amplitude of the reconstructed seismic trace;

[0060] S10. Calculate the amplitude equalization coefficient based on the average amplitude of the non-missing traces and the average amplitude of the reconstructed seismic trace;

[0061] S11. Correct the newly reconstructed seismic data based on the amplitude equalization coefficient.

[0062] 8. A method for reconstructing seismic data according to claim 7, wherein

[0063] The calculation method of the amplitude equalization coefficient w i is specifically as follows:

[0064]

[0065] where A n is the average amplitude of the non-missing traces, and A r is the average amplitude of the seismic trace obtained by reconstruction;

[0066]

[0067]

[0068] where N is the number of sampling points of this trace; f i,j is the seismic record data; i is the trace number; j is the time number.

[0069] The beneficial effects of this technical solution are as follows: By calculating and applying the amplitude equalization coefficient, the amplitude of the reconstructed seismic data can be effectively adjusted to make it closer to the true situation of the original seismic data; Amplitude equalization ensures that the reconstructed data is consistent with the non-missing data in terms of amplitude, maintaining the consistency of the data set; Through this systematic amplitude restoration step, the process of reconstructing seismic data can be simplified, making it more efficient and easier to operate.

[0070] Further, in step S11, specifically, multiply the uncorrected newly reconstructed seismic data by the amplitude equalization coefficient to obtain the corrected newly reconstructed seismic data.

[0071] In the second aspect, the present invention also provides a computer-readable storage medium, which stores a computer program. When the computer program is executed by a processor, it implements a method for reconstructing seismic data according to any one of the above. BRIEF DESCRIPTION OF THE DRAWINGS

[0072] The drawings herein are incorporated into and constitute a part of this specification, showing embodiments consistent with the present disclosure, and are used together with the specification to explain the principles of the present disclosure.

[0073] To more clearly illustrate the technical solutions in the embodiments of the present disclosure or the prior art, the following will briefly introduce the drawings required for use in the description of the embodiments or related technologies. Obviously, for those of ordinary skill in the art, without creative efforts, other drawings can also be obtained based on these drawings.

[0074] Figure 1 Schematically shows a flowchart of a method according to an embodiment of the present disclosure;

[0075] Figure 2 Schematically shows a seismic data map according to an embodiment of the present disclosure;

[0076] Figure 3 Schematically shows a missing seismic data map according to an embodiment of the present disclosure;

[0077] Figure 4 Schematically shows a Radon forward transform result map of seismic data according to an embodiment of the present disclosure;

[0078] Figure 5 Schematically shows a reconstructed new seismic data map according to an embodiment of the present disclosure. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0079] To make the objectives, technical solutions, and advantages of the embodiments of the present disclosure clearer, the following will clearly and completely describe the technical solutions in the embodiments of the present disclosure with reference to the drawings in the embodiments of the present disclosure. Obviously, the described embodiments are part of the embodiments of the present disclosure, rather than all of them. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present disclosure without creative efforts belong to the scope of protection of the present disclosure.

[0080] It should be noted that in this text, relational terms such as "first" and "second" are only used to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any actual relationship or order between these entities or operations. Moreover, the terms "include", "comprise" or any other variants thereof are intended to cover non-exclusive inclusion, so that a process, method, article or device comprising a series of elements not only includes those elements, but also includes other elements not expressly listed, or also includes elements inherent to such process, method, article or device. Without further limitation, an element defined by the statement "comprising a..." does not exclude the existence of additional identical elements in the process, method, article or device comprising the said element.

[0081] Next, a simple introduction to one of the application fields of the embodiments of the present application (i.e., seismic data reconstruction) will be given first.

[0082] The Radon transform is an integral transform widely used in fields such as image processing, computer vision, seismology, and medical imaging (such as CT scans). This transform can convert a function (such as an image or seismic data) from its original space to a parameter space, providing an effective means for analyzing and processing such data; the Radon transform works by mapping a two-dimensional function (e.g., an image or a seismic record) onto a set of line integrals. In seismology, this typically involves converting seismic data from the time-space domain (e.g., time and offset) to the parameter domain (such as intercept time and slowness). The Radon transform can effectively separate and extract specific signals from complex data. For example, it can separate direct waves and reflected waves from seismic data. There are various variants of the Radon transform, such as the filtered backprojection algorithm (commonly used in CT scans) and the higher-order Radon transform (used in seismic data processing to handle more complex wavefield characteristics).

[0083] The Radon inverse transform is the process of converting the data obtained through the Radon transform back to its original space (such as the time-space domain of an image or seismic data). The basic objective of the Radon inverse transform is to reconstruct the original function (such as an image or a seismic record) from the data in the Radon transform domain (parameter space). In the Radon transform, the original data is mapped onto a series of line integrals; while in the inverse transform process, the original data is reconstructed from these line integrals. The Radon inverse transform can be implemented through various mathematical methods. One common method is the filtered backprojection algorithm, whose basic idea is to perform appropriate filtering on the data in the Radon transform domain and then backproject (i.e., reverse map) it back to the original space. In seismology, the Radon inverse transform is used to reconstruct the original seismic record from the data in the transform domain, especially in the migration and velocity analysis of seismic data.

[0084] Seismic AVO (Amplitude Versus Offset) characteristics is a method in seismology for analyzing and interpreting how the reflection amplitude of seismic waves changes with the offset or incident angle; when seismic waves pass through different geological interfaces, their reflection amplitudes change with the incident angle of the waves, and these changes can reveal important information about the properties of underground rocks and fluids; in seismic exploration, the offset refers to the horizontal distance between the seismic source (such as the explosion point) and the geophone, and as the offset increases, the incident angle of the seismic waves also changes, thus affecting the amplitude of the reflected wave.

[0085] The Matching Pursuit (MP) algorithm is a signal processing technique used to decompose a signal into a linear combination of a series of predefined atoms or basis functions; the core idea of the matching pursuit algorithm is to select a series of atoms (basis functions) from a large dictionary that can best approximate the original signal, and the algorithm iteratively selects a single atom, with each selection based on the current residual (i.e., the difference between the original signal and the current approximation). In seismology, matching pursuit is used to extract specific waveforms or events from seismic data. However, for a large dictionary, the algorithm may take a long time to find the best-matching atom, and the algorithm performance highly depends on the choice of the dictionary. An inappropriate dictionary may lead to poor performance, that is, in the present invention, an accurate value of p is required to support the efficient calculation of the matching pursuit algorithm.

[0086] As Figure 1 shown, a method for reconstructing seismic data includes the following steps,

[0087] Perform a Radon transform on the seismic data (as Figure 2 shown) to transform the missing seismic data (as Figure 3 shown) from the time domain to the Radon transform domain (as Figure 4 shown); perform an inverse transform on the high-order Radon transform to obtain the missing part of the seismic data; fill the missing part of the data into the seismic data to obtain the newly reconstructed seismic data (as Figure 5 shown).

[0088] First, perform a Radon transform on the seismic data with missing data to transform it from the time domain to the Radon transform domain; combine the seismic data transformed to the Radon transform domain with the orthogonal polynomials characterizing seismic AVO characteristics to form a high-order Radon transform; this step stacks the seismic data through an integration path, thereby mapping the data in the time domain to the Radon transform domain, and is implemented by integrating and stacking along the integration path of the hyperbolic Radon transform; it can more clearly express and extract important features in the seismic data, laying a foundation for subsequent processing.

[0089] Combine the seismic data in the Radon transform domain with orthogonal polynomials characterizing the seismic AVO characteristics to form a high-order Radon transform. By adding orthogonal polynomials describing the variation of amplitude with offset to the Radon transform to improve data representation, this step provides a method that can more finely characterize the characteristics of seismic data, especially when analyzing the amplitude variation in seismic data.

[0090] Perform an inverse transform on the obtained high-order Radon transform to recover the missing part of the seismic data. The inverse transform is the process of converting the data in the Radon transform domain back to the original time-space domain. Through this step, the missing part in the seismic data can be effectively reconstructed and the integrity of the data can be restored.

[0091] Fill the reconstructed missing part data back into the original seismic data to obtain the newly reconstructed seismic data, and continue to iterate this process until a satisfactory reconstruction result is obtained. The reconstructed data is used to replace the missing part in the original dataset, and then continue to perform the Radon transform and inverse transform processes. This iterative process helps to gradually improve the quality and accuracy of the reconstructed data, especially when dealing with seismic data containing complex geological features.

[0092] After the step of filling the missing part data into the seismic data to obtain the newly reconstructed seismic data, it further includes continuing to iterate with the newly reconstructed seismic data as the missing seismic data.

[0093] Among them, transform it from the time domain to the Radon transform domain. The Radon transform is defined as:

[0094]

[0095]

[0096] Its discrete form is:

[0097]

[0098]

[0099] In the formula: τ is the intercept time; p is the slowness (the reciprocal of the velocity); x is the offset or the migration offset; is the data after transformation to the Radon domain; d(t, x) is the two-dimensional common shot gather seismic data.

[0100] The key to the forward Radon transform is to find basis functions that better match the actual data from the space of orthogonal polynomial surfaces with different curvatures, so as to obtain high-resolution results. The present invention uses the matching pursuit method to implement the high-resolution Radon transform. The matching pursuit algorithm is an iterative algorithm. Its basic principle is to find the basis function that best matches the original data at each iteration and perform the transform in this subspace; repeat the above matching steps for the residual data until the matching error is minimized.

[0101] The variation of seismic amplitude with offset can be represented by polynomials. The orthogonal polynomial transform characterizes the seismic AVO characteristics with fewer polynomial coefficients. Therefore, a transform space describing the amplitude variation is introduced, which together with the integral path space forms an overcomplete space. Orthogonal polynomials representing the variation of amplitude with offset are added to the overcomplete dictionary, and the Radon transform is combined with the orthogonal polynomial transform to form the high-order Radon transform.

[0102] The step of combining the seismic data transformed into the Radon transform domain with the orthogonal polynomials characterizing the seismic AVO characteristics to form the high-order Radon transform includes calculating the conjugate solution of the high-order Radon transform and calculating the high-order Radon transform using the matching pursuit algorithm;

[0103] Among them, the specific steps for calculating the conjugate solution of the high-order Radon transform are as follows:

[0104] {q j (x), j = 0, 1, …, N} is a set of orthonormal polynomials determined by N offset coordinates x

[0105] q j (x) is an algebraic polynomial of order j;

[0106] The amplitude variation of the seismic data d(t, x) at a certain moment can be fitted by orthogonal polynomials, that is:

[0107]

[0108] Among them, c j (t) is the decomposition coefficient of the orthogonal polynomial of order j for the amplitude variation with offset at time t, which is called the orthogonal polynomial coefficient spectrum. Combining Equation (2) and Equation (5) can obtain:

[0109]

[0110] Since {q j (x), j = 0, 1, …, N} is a set of orthogonal polynomials, multiplying both sides of Equation (6) by q j (x) can obtain the conjugate solution of the high-order Radon transform:

[0111]

[0112] The specific steps for calculating the high-order Radon transform using the matching pursuit algorithm are as follows:

[0113] Represent the matrix of the said Radon transform as:

[0114] M = DL (8)

[0115] where L is the integral path function;

[0116] S1. Let the seismic data D be the initial residual D resi , and the Radon parameter M be 0;

[0117] S2. Calculate the conjugate solution M of the Radon transform of the residual data using Equation (7) i ;

[0118] S3. Calculate the Radon energy distribution using Equation (9), and determine the subspace S through the local maximum of the energy i ;

[0119]

[0120] S4. The said energy represents the energy distribution of the in-phase axis along different curvature directions at different times;

[0121] S5. Minimize the function J(m) = ‖D resi - LM i ‖ 2 using the conjugate gradient method to obtain the least squares solution M of the Radon transform in the above subspace S i ; i ;

[0122] S6. Update the parameter M = M + M i , and the residual data D resi = D resi - LM i ;

[0123] S7. Determine whether the predefined number of iterations is reached. If the predefined number of iterations is reached, output all subspaces S i that constitute the high-order Radon transform. If the predefined number of iterations is not reached, return to S1 to continue the iteration.

[0124] To avoid aliasing during data reconstruction and obtain a high-resolution sampling result, it is necessary to reasonably select the maximum and minimum values of the slowness parameter, i.e., p, and the sampling interval. The selection criteria for the parameters are discussed below. The sampling rate of the slowness parameter is:

[0125] Δp = p k - p k-1 (10)

[0126] where p k and p k-1 are any two adjacent slowness values.

[0127] Given p = 1 / v, equation (3) can be transformed into:

[0128]

[0129] Let Δv = v k - v k-1 , we get:

[0130]

[0131] When v k and v k-1 are relatively close, equation (5) can become:

[0132]

[0133] Suppose the maximum effective reflection velocity in the original data is v max . To make equation (6) hold for all effective wave fields, then v k in the equation should be replaced by the maximum effective reflection velocity v max of the original data. At the same time, when performing the transformation on the CSP or CMP gather, the selected range of Δv is usually 50 - 100 m / s. Thus, the critical value relationship is as follows:

[0134]

[0135] After determining the value of Δp, according to p = 1 / v, the initial value p0 and the maximum value p max of p are given, and the reconstruction accuracy can be guaranteed without increasing the computational amount of data reconstruction.

[0136] In this process, even if there is a small amount of missing seismic data, the inverse Radon transform can restore the seismic data to a complete hyperbolic form. The inverse Radon transform formula is:

[0137]

[0138] The previously missing part of the seismic traces will generate seismic data. Use the reconstructed data of these missing positions to replace the missing positions in the original data, and then perform the transformation and iterate continuously until a better effect is obtained. However, in the process of the forward transform and the inverse transform, the reconstructed seismic data will lose amplitude information. To compensate for this loss, the present invention uses the amplitude equalization coefficient w i to solve this problem:

[0139]

[0140] where A n is the average amplitude of the non-missing traces; A r is the average amplitude of the seismic traces obtained through reconstruction.

[0141] The average amplitude of the non-missing traces near the missing seismic trace is as follows:

[0142]

[0143] N is the number of sampling points of this trace; f i,j is the seismic recording data; i is the trace number; j is the time number. Similarly, the average amplitude of the reconstructed seismic trace can be written as:

[0144]

[0145] Using the reconstructed trace processed by amplitude equalization for calculation during the iteration process can restore the proper amplitude size and can also improve the iteration efficiency.

[0146] The embodiment of the present application also provides a computer-readable storage medium. The computer-readable storage medium stores a computer program. When the computer program is executed by a processor, the steps of any of the above methods are implemented. Its specific implementation manner is the same as the implementation manner and the achieved technical effects recorded in the above method embodiment, and some contents will not be elaborated.

[0147] The above is only a preferred specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any person skilled in the art within the technical scope disclosed by the present invention, according to the technical solution of the present invention and its inventive concept, makes equivalent substitutions or changes, and should be covered by the protection scope of the present invention.

Claims

1. A method for reconstructing seismic data, characterized in that, It includes the following steps: Perform a Radon transform on the seismic data to transform the missing seismic data from the time domain to the Radon transform domain; Combine the seismic data transformed to the Radon transform domain with orthogonal polynomials characterizing the seismic AVO characteristics to form a high-order Radon transform; Perform an inverse transform on the high-order Radon transform to obtain the missing data of the seismic data; Fill the missing data into the seismic data to obtain a newly reconstructed seismic data.

2. A method for reconstructing seismic data according to claim 1, wherein: After the step of filling the missing data into the seismic data to obtain a newly reconstructed seismic data, it further includes iterating the newly reconstructed seismic data as the missing seismic data.

3. A method for reconstructing seismic data according to claim 2, wherein: The expression of the Radon transform is: Its discrete form is: where τ is the intercept time; p is the reciprocal of the velocity, i.e., slowness; x is the shot-receiver distance or offset; m(τ, p) is the seismic data transformed to the Radon domain; d(t, x) is the two-dimensional common-shot gather seismic data, and i and k are ordinals.

4. A method for reconstructing seismic data according to claim 2, wherein: The step of combining the seismic data transformed to the Radon transform domain with orthogonal polynomials characterizing the seismic AVO characteristics to form a high-order Radon transform includes calculating the conjugate solution of the high-order Radon transform and calculating the high-order Radon transform using the matching pursuit algorithm; Among them, the specific steps for calculating the conjugate solution of the high-order Radon transform are: {q j (x), j = 0, 1, …, N} is a set of orthonormal polynomials determined by N offset coordinates x q j (x) is an algebraic polynomial of order j; The amplitude change of the seismic data d(t, x) at a certain moment can be fitted with orthogonal polynomials, that is: where c j (t) is the decomposition coefficient of the j-th order orthogonal polynomial of the amplitude varying with the offset at time t, which is called the orthogonal polynomial coefficient spectrum. Combining Equation (2) and Equation (5) can obtain: Since {q j (x), j = 0, 1, …, N} is a set of orthogonal polynomials, multiplying both sides of Equation (6) by q j (x) can obtain the conjugate solution of the said high-order Radon transform:

5. A method for reconstructing seismic data according to claim 4, wherein: The specific steps for calculating the high-order Radon transform using the matching pursuit algorithm are: Represent the matrix of the Radon transform as: M = DL (8) where L is the integral path function; S1. Let the seismic data D be the initial residual D resi , and the Radon parameter M be 0; S2. Calculate the conjugate solution M of the Radon transform of the residual data using Equation (7) i ; S3. Calculate the Radon energy distribution using formula (9), and determine subspace S by the local maximum of the energy i ; S4. The energy represents the energy distribution of the in-phase axis along different curvature directions at different times; S5. Minimize the function J(m) = ‖D resi - LM i ‖ 2 , and obtain the least - squares solution M i of the Radon transform in the above subspace S i ; S6. Update the parameter M = M + M i , the residual data D resi = D resi - LM i ; S7. Determine whether the predefined number of iterations is reached. If the predefined number of iterations is reached, output all subspaces S i to form a high-order Radon transform. If the predefined number of iterations is not reached, return to S1 to continue the iteration.

6. A method for reconstructing seismic data according to claim 4, wherein: The specific method for calculating the value range of p is: Assume that the sampling rate of p is Δp, then Δp = p k - p k-1 (10) wherein, in the formula, p k and p k-1 are any two adjacent slowness values; v max is the maximum effective reflection velocity in the seismic data, then Δp is expressed as: where Δv is the value of v k -v k-1 and After determining the value of Δp, the initial value p0 and the maximum value p of p are given according to p = 1 / v max , and the range of values of p is between the initial value and the maximum value.

7. A method for reconstructing seismic data according to claim 4, wherein: It further includes restoring the amplitude through an amplitude equalization coefficient during the iteration process, and the specific steps are: S8. Calculate the average amplitude of the non-missing traces near the missing seismic trace; S9. Calculate the average amplitude of the reconstructed seismic trace; S10. Calculate the amplitude equalization coefficient based on the average amplitude of the non-missing traces and the average amplitude of the reconstructed seismic trace; S11. Correct the newly reconstructed seismic data based on the amplitude equalization coefficient.

8. A method for reconstructing seismic data according to claim 7, wherein: The calculation method of the amplitude equalization coefficient w i is specifically as follows: where A n is the average amplitude of the non-missing traces, and A r is the average amplitude of the reconstructed seismic traces; where N is the number of sampling points in this trace; f i,j is the seismic recording data; i is the trace number; j is the time number.

9. A method for reconstructing seismic data according to claim 8, wherein: In step S11, specifically, the uncorrected reconstructed new seismic data is multiplied by the amplitude equalization coefficient to obtain the corrected reconstructed new seismic data.

10. A computer-readable storage medium, characterized in that the computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, it implements the seismic data reconstruction method according to any one of claims 1-9.

Images