Multivariate interpolation method for suppressing noise

Abstract

The application provides a multi-element interpolation method for suppressing noise, comprising the following steps: step 1, data preparation; step 2, spatial sampling design, a pre-interpolation spatial model is established according to a field observation mode; step 3, establishment of an estimated spectrum, seismic data is converted from a time-offset domain to a frequency-wavenumber domain, and then the estimated spectrum is established; step 4, multi-element interpolation in the frequency-wavenumber domain, interpolation is performed according to the pre-interpolation spatial model; step 5, removal of interference waves; step 6, removal of new interpolation data according to the pre-interpolation spatial model, only data after denoising of original seismic data is reserved; and step 7, discrete Fourier inverse transformation, seismic data after noise removal is inversely transformed to the time-offset domain. The multi-element interpolation method for suppressing noise simulates a wave field with sufficient sampling and complete wave field by using a multi-element interpolation method for data with insufficient sampling and incomplete wave field, so that the noise is further suppressed.

Description

Technical Field

[0001] This invention relates to the field of geophysical exploration technology, and in particular to a method for suppressing noise using multivariate interpolation. Background Technology

[0002] Desert seismic data commonly suffers from interference from first-arrival multiples and surface waves. These two types of interference waves have a significant impact on seismic data and are difficult to remove. Based on the characteristics of first-arrival multiples and surface waves, many suppression methods have been proposed, such as Radon transform, FK filtering, and regional filtering. These methods each have their own characteristics and have achieved significant results, but they all focus on only one aspect of the characteristics or rely on certain assumptions, and thus have their limitations. Radon transform, for example, suppresses noise by recognizing the difference in the apparent velocity sign and magnitude of effective waves and noise. While first-arrival multiples and surface waves are regular linear noise, their distribution in seismic records exhibits a severe broom-like pattern from shallow to deep, with variations in velocity and frequency. Transforming seismic data containing first-arrival multiples and surface waves into the Radon domain reveals that these are not isolated points, making complete removal difficult. FK filtering requires regular spatial sampling intervals and a sound wavefield; it is ineffective or even non-effective for data with irregular sampling or incomplete wavefields. Regional filtering, while avoiding the influence of low-cutoff filtering on signals outside the surface wave region, also filters effective low-frequency signals within the region, exhibiting boundary effects and significantly losing mid-to-deep low-frequency effective wave information. Existing denoising techniques for first-arrival multiples and surface wave interference all have certain limitations.

[0003] First-arrival multiples and surface waves are widely present in desert seismic data, both characterized by high energy and linear relationships. The presence of these two types of interference waves affects the signal-to-noise ratio and resolution of the data. If they are not effectively suppressed before stacking, the overall quality of data processing will inevitably be affected. How to avoid damaging the effective signal while suppressing these two types of interference waves has always been an important aspect of seismic data processing. Experiments show that to obtain a sound wavefield for first-arrival multiples and surface waves, the trace spacing needs to be within 10 meters. However, due to exploration cost constraints, the trace spacing of existing desert seismic data is mostly 50 meters or 25 meters, which leads to insufficient sampling and an incomplete wavefield.

[0004] Chinese patent application CN201611240359.1 discloses a method for suppressing noise in the curve wave domain using Radon transform in the Loess Plateau region. The method includes: performing a forward curve wave transform on the original seismic data to obtain matrices of curve wave coefficient domains with different directions and scales; performing a forward Radon transform on each matrix to transform the data in the curve wave coefficient domain to the Radon domain; performing threshold filtering on each matrix; since interference waves in the Radon domain appear as points or energy clusters with strong energy, and random noise appears as unfocused energy, while effective waves, although unfocused, still exhibit focused energy greater than random noise, threshold filtering can effectively eliminate interference waves; transforming the data back from the Radon domain to the curve wave coefficient domain using an inverse Radon transform; and performing an inverse curve wave transform on the data in the curve wave domain to obtain the final data.

[0005] Chinese patent application CN201510379743.9 discloses a method for suppressing multiples in deep-water OBC data based on model statistical picking, comprising the following steps: 1) preprocessing of raw seismic data; 2) sorting of preprocessed seismic data gathers; 3) determining the multiple phase axes to be suppressed; 4) reading a gather and picking out the multiple phase axes to be suppressed from that gather; 5) smoothing the picked multiple phase axes; 6) using the smoothed multiple phase axes as the center... 7) Static time shift is performed on each trace in the common offset gather to flatten the smoothed multiple phase axis. Then, the amplitude characteristics of multiple waves in each trace within the time window are statistically analyzed along the multiple phase axis to establish a multiple model; 8) The multiple model is subtracted from the flattened seismic data; 9) The primary wave data obtained after subtracting the multiple model is statically time shifted in the opposite direction along the centerline of the multiple phase axis in step 6); 10) Steps 4) to 8) are repeated on the next gather using the smoothed multiple phase axis obtained in step 5) as a reference.

[0006] Chinese patent application CN201510883681.5 discloses a method for suppressing sliding scan harmonic interference based on frequency division processing. After acquiring seismic data, depending on the amount of data, the data for the entire area is grouped, either by one beamline or multiple beamlines, re-sorting the seismic data from the shot domain to the common center point domain. Because sliding scan harmonic interference has strong energy characteristics and is concentrated near the offset on a single shot, its identification and suppression in the shot domain are difficult. However, after the data is sorted to other data domains, the distribution of harmonic interference is relatively dispersed. Based on its strong energy characteristics, corresponding anomalous energy identification and suppression methods can be designed to suppress it.

[0007] The existing technologies described above are significantly different from the present invention and have failed to solve the technical problem we want to address. Therefore, we have invented a new multivariate interpolation method for suppressing noise. Summary of the Invention

[0008] The purpose of this invention is to provide a multivariate interpolation method for noise suppression that simulates a fully sampled and complete wave field from insufficiently sampled and incomplete data using multivariate interpolation, thereby further suppressing noise.

[0009] The objective of this invention can be achieved through the following technical measures: a multivariate interpolation method for suppressing noise, wherein the first step of the multivariate interpolation method for suppressing noise is: data preparation;

[0010] Step 2: Perform spatial sampling design and establish a pre-interpolated spatial model based on field observation methods;

[0011] Step 3: Establish the estimation spectrum by converting the seismic data from the time-offset domain to the frequency-wavenumber domain, and then establishing the estimation spectrum.

[0012] Step 4: Perform frequency-wavenumber domain multivariate interpolation according to the pre-interpolation space model;

[0013] Step 5: Perform interference removal, and perform linear interference removal on the well-sampled and sound wavefield seismic data;

[0014] Step 6: Remove the new interpolated data according to the pre-interpolated spatial model, and only retain the data after denoising the original seismic data;

[0015] Step 7: Perform a discrete inverse Fourier transform to transform the noise-removed seismic data back to the time-offset domain.

[0016] The objective of this invention can also be achieved through the following technical measures:

[0017] In step 1, data preparation is performed by loading the seismic data collected in the field into the observation system and performing static correction.

[0018] In step 1, the data preparation includes loading the seismic data acquired in the field into the observation system to obtain information such as the coordinates, elevation, and trace distance of the shot point and receiver point, while calculating information such as azimuth and offset, and adding a three-dimensional grid to the seismic data.

[0019] In step 1, if the seismic data has undergone tomographic static correction or refraction static correction, the tomographic or refraction static correction values ​​should be applied; if no tomographic or refraction static correction has been performed, the elevation static correction values ​​should be calculated and applied.

[0020] In step 2, spatial sampling design is carried out. Based on the positions of the shot points and receiver points in the field observation system, new shot points and receiver points for pre-interpolation are designed. The field observation system and the pre-interpolated shot points and receiver points are merged to form a pre-interpolated spatial model.

[0021] In step 3, an estimation spectrum is established by converting the seismic data from the time-offset domain to the frequency-wavenumber domain. For a single-frequency seismic data, the corresponding wave number is first obtained from the maximum Fourier spectrum. Then, the wave number is removed from the single-frequency seismic data to obtain a new single-frequency seismic data. The same method is then used to obtain the wave number corresponding to the maximum Fourier spectrum for this new single-frequency seismic data. This process is repeated until the maximum Fourier spectrum is less than 0.05. All the obtained wave numbers are then combined to obtain the entire estimation spectrum.

[0022] In step 4, frequency-wavenumber domain multivariate interpolation is performed. Based on multivariate information such as the horizontal and vertical line numbers, time, offset, and azimuth of the seismic data, independent wavenumber scanning is conducted to determine the optimal interpolation wavenumber, and interpolation is performed according to the pre-interpolation spatial model.

[0023] In step 4, when determining the optimal interpolation wavenumber, for interpolation at any location of the pre-interpolation spatial model, the maximum Nyquist wavenumber is first determined based on the spatial sampling interval of the seismic data. The spatial wavenumber is then continuously reduced at certain intervals, and the minimum root mean square error corresponding to these wavenumbers is calculated respectively. The maximum wavenumber that meets the set minimum root mean square error threshold value is selected, and then the interpolation weight coefficient corresponding to the wavenumber is calculated.

[0024] In step 5, when removing interference waves, the seismic data interpolated in the previous step is in the frequency-wavenumber domain and is a well-sampled seismic data with a sound wavefield. Therefore, the linear interference relationship between the first arrival multiple waves and surface waves is clear. Based on this linear relationship, apparent velocity filtering can be performed to remove this interference wave.

[0025] In step 7, a discrete inverse Fourier transform is performed to transform the noise-removed seismic data to the time-offset domain, resulting in the final noise-suppressed data.

[0026] The multivariate interpolation method for noise suppression in this invention falls under the technical category of removing linear interference waves and improving the signal-to-noise ratio in seismic data processing. It achieves the goal of maximizing the protection of effective information when suppressing first-to-multiple waves and surface waves. The beneficial effects of this invention are as follows:

[0027] First, the multivariate interpolation method effectively solves the noise reduction problem caused by insufficient seismic data sampling and incomplete wavefield due to exploration costs.

[0028] Second, by establishing an estimated spectrum after converting the seismic data from the time-offset domain to the frequency-wavenumber domain, instead of directly using the original Fourier spectrum, the spurious frequency phenomenon caused by uneven spatial sampling is effectively avoided.

[0029] Third, using the optimal interpolation wavenumber in the multivariate interpolation process can effectively solve the complex spectrum problem, thereby ensuring the accuracy of interpolation.

[0030] Fourth, the multivariate interpolation method for noise suppression can effectively remove first arrival multiples and surface waves, improve the signal-to-noise ratio of seismic data, and provide more detailed basic data for subsequent well location verification. Attached Figure Description

[0031] Figure 1 This is a flowchart of a specific embodiment of the noise suppression method using multivariate interpolation of the present invention;

[0032] Figure 2 This is a raw single-shot data diagram containing first to multiple waves and surface waves in a specific embodiment of the present invention;

[0033] Figure 3 This is a spectrum analysis diagram of the surface wave of the original single-shot data in a specific embodiment of the present invention;

[0034] Figure 4 This is a frequency-wavenumber spectrum of the original single-shot data in a specific embodiment of the present invention;

[0035] Figure 5 This is a data graph generated by multivariate interpolation in a specific embodiment of the present invention;

[0036] Figure 6 This is a frequency-wavenumber spectrum of the interpolated data in a specific embodiment of the present invention;

[0037] Figure 7 This is a single-shot data diagram of noise suppression after multivariate interpolation according to a specific embodiment of the present invention;

[0038] Figure 8 This is a specific embodiment of the present invention, showing the final single-shot data image after processing with the multivariate interpolation method for noise suppression.

[0039] Figure 9 This is a single-shot data image obtained by conventional denoising of the original single-shot data in a specific embodiment of the present invention. Detailed Implementation

[0040] It should be noted that the following detailed descriptions are exemplary and intended to provide further illustration of the invention. Unless otherwise specified, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains.

[0041] It should be noted that the terminology used herein is for the purpose of describing particular embodiments only and is not intended to limit the exemplary embodiments of the present invention. As used herein, the singular form is intended to include the plural form as well, unless the context clearly indicates otherwise. Furthermore, it should be understood that when the terms "comprising" and / or "including" are used in this specification, they indicate the presence of features, steps, operations, and / or combinations thereof.

[0042] The noise suppression method using multivariate interpolation of the present invention includes the following steps:

[0043] Step 1: Data preparation. Load the seismic data collected in the field into the observation system and perform static correction.

[0044] The data preparation includes loading the seismic data collected in the field into the observation system to obtain information such as the coordinates, elevation, and trace distance of the shot point and receiver point, while calculating information such as azimuth and offset, and adding a three-dimensional grid to the seismic data.

[0045] If the seismic data has undergone tomographic static correction or refraction static correction, the tomographic or refraction static correction values ​​should be applied; if no tomographic or refraction static correction has been performed, the elevation static correction values ​​should be calculated and applied.

[0046] Step 2: Spatial sampling design, establishing a pre-interpolated spatial model based on field observation methods;

[0047] Spatial sampling design is carried out. Based on the positions of the shot points and receiver points in the field observation system, new shot points and receiver points for pre-interpolation are designed. The field observation system and the pre-interpolated shot points and receiver points are merged to form a pre-interpolated spatial model.

[0048] Step 3: Establish the estimation spectrum by converting the seismic data from the time-offset domain to the frequency-wavenumber domain, and then establishing the estimation spectrum.

[0049] To establish an estimation spectrum, the seismic data is converted from the time-offset domain to the frequency-wavenumber domain. For a single-frequency seismic data set, the corresponding wave number is first obtained from the maximum Fourier spectrum. Then, this wave number is removed from the single-frequency seismic data set to obtain a new single-frequency seismic data set. The same method is then used to obtain the wave number corresponding to the maximum Fourier spectrum for this new single-frequency seismic data set. This process is repeated until the maximum Fourier spectrum is less than 0.05. All the obtained wave numbers are then combined to obtain the entire estimation spectrum.

[0050] Step 4: Frequency-wavenumber domain multivariate interpolation, interpolating according to the pre-interpolation space model;

[0051] Frequency-wavenumber domain multivariate interpolation involves independently scanning wavenumbers based on multivariate information from seismic data, including line numbers, time, offset, and azimuth, to determine the optimal interpolation wavenumber. Interpolation is then performed according to a pre-interpolation spatial model. Determining the optimal interpolation wavenumber involves first determining the maximum Nyquist wavenumber based on the spatial sampling interval of the seismic data when interpolating at any location in the pre-interpolation spatial model. The spatial wavenumber is then progressively reduced at certain intervals, and the minimum root mean square error (RMSE) corresponding to each wavenumber is calculated. The maximum wavenumber that meets the set RMSE threshold is then selected, and the interpolation weighting coefficients corresponding to this wavenumber are calculated.

[0052] Step 5: Interference removal. Perform linear interference removal on well-sampled seismic data with a sound wavefield.

[0053] Interference wave removal: The seismic data after the previous interpolation step is in the frequency-wavenumber domain and is a well-sampled seismic data with a sound wavefield. Therefore, the linear interference relationship between the first arrival multiple waves and surface waves is clear. Based on this linear relationship, apparent velocity filtering can be performed to remove this interference wave.

[0054] Step 6: Remove the new interpolated data according to the pre-interpolated spatial model, and only retain the data after denoising the original seismic data;

[0055] The new interpolated seismic data is removed according to the pre-interpolated spatial model, and only the denoised data of the original seismic data is retained, which can ensure the fidelity of the seismic data.

[0056] Step 7: Discrete Fourier transform, transforming the noise-removed seismic data inversely to the time-offset domain.

[0057] The Inverse Discrete Fourier Transform (IDFT) transforms the noise-reduced seismic data back to the time-offset domain, yielding the final noise-suppressed data. The IFT avoids transform aliasing caused by incomplete or irregular sampling of field data.

[0058] In a specific embodiment 1 of the present invention, the multivariate interpolation method for suppressing noise includes the following steps:

[0059] (1) Data preparation: This includes loading the seismic data acquired in the field into the observation system to obtain information such as the coordinates, elevation, and trace spacing of the shot points and receivers. At the same time, the azimuth and offset information are calculated using the coordinates of the shot points and receivers, and a three-dimensional grid is added to the seismic data. If the seismic data has undergone tomographic static correction or refraction static correction, the tomographic or refraction static correction values ​​should be applied; if no tomographic or refraction static correction has been performed, the elevation static correction values ​​should be calculated and applied.

[0060] (2) Perform spatial sampling design to prepare for the subsequent interpolation point locations. Based on the shot point and receiver point locations in the field observation system, design new shot points and receiver points for pre-interpolation. The field observation system and the pre-interpolated shot points and receiver points are merged to form a pre-interpolation spatial model.

[0061] (3) Establish the estimation spectrum. Convert the seismic data from the time-offset domain to the frequency-wavenumber domain. For a single-frequency seismic data, first obtain the corresponding wave number based on the maximum Fourier spectrum, then remove this wave number from the single-frequency seismic data to obtain a new single-frequency seismic data. Then, use the same method to obtain the wave number corresponding to the maximum Fourier spectrum for this new single-frequency seismic data. Repeat this process until the maximum Fourier spectrum is less than 0.05. Combine all the obtained wave numbers to obtain the estimation spectrum of a single-frequency seismic data. Then, combine the estimation spectra of all single-frequency seismic data to obtain the entire estimation spectrum.

[0062] (4) Frequency-Wavenumber Domain Multivariate Interpolation. Based on multivariate information such as the horizontal and vertical line numbers, time, offset, and azimuth of the seismic data, independent wavenumber scanning is performed to determine the optimal interpolation wavenumber. Interpolation is then performed according to the pre-interpolation spatial model. To determine the optimal interpolation wavenumber, when interpolating at any location in the pre-interpolation spatial model, the maximum Nyquist wavenumber is first determined based on the spatial sampling interval of the seismic data. The spatial wavenumber is then continuously reduced at certain intervals, and the minimum root mean square error corresponding to each wavenumber is calculated. The maximum wavenumber that meets the set minimum root mean square error threshold is selected, and then the interpolation weight coefficients corresponding to that wavenumber are calculated.

[0063] (5) Removal of interference waves such as first arrival multiples and surface waves. The seismic data after the previous interpolation is in the frequency-wavenumber domain and has been interpolated into a seismic data with sufficient sampling and a sound wavefield. Therefore, the linear interference relationship of first arrival multiples and surface waves is clear. Based on this linear relationship, apparent velocity filtering can be performed to remove such interference waves.

[0064] (6) Removal of new interpolated data. New interpolated seismic data are removed according to the pre-interpolated spatial model, and only the data after denoising the original seismic data is retained. This can ensure the fidelity of the seismic data.

[0065] (7) Discrete Fourier Transform. The noise-removed seismic data is inversely transformed back to the time-offset domain to obtain the finally suppressed noise data. The Discrete Fourier Transform can avoid transform aliasing caused by incomplete and irregular sampling of field data.

[0066] In a specific embodiment 2 of the present invention, the multivariate interpolation method for suppressing noise includes the following steps:

[0067] Load the observation system, calculate information such as azimuth and offset, and apply static correction values;

[0068] Spatial sampling design yields a pre-interpolated spatial model;

[0069] The estimated spectrum of single-frequency seismic data is calculated in the frequency-wavenumber domain, and the estimated spectra of single-frequency seismic data are combined to obtain the entire estimated spectrum.

[0070] By utilizing multivariate information, independent wavenumber scans are performed to determine the optimal interpolation wavenumber, and interpolation is performed according to the pre-interpolation space model.

[0071] Remove interference waves such as first-arrival multiples and surface waves;

[0072] Remove the new interpolated data and retain only the data in the original space;

[0073] The discrete inverse Fourier transform returns the data to the time-offset domain, yielding the denoised data.

[0074] For ease of explanation, actual seismic data from a certain exploration area were used, and noise suppression was performed using the multivariate interpolation method according to the operating steps of this invention to illustrate the application effect of this invention. The actual data was acquired using a 40-line, 14-shot observation system, with a seismic data duration of 3000 milliseconds, a time sampling interval of 1 millisecond, 310 traces per line, and a trace spacing of 50 meters.

[0075] In a specific embodiment 3 of the present invention, the invention is applied. Figure 1 The diagram shown is a flowchart of the noise suppression technology using multivariate interpolation according to the present invention.

[0076] Step 101, Data Preparation. This includes loading the seismic data acquired in the field onto the observation system to obtain information such as the coordinates, elevation, and trace spacing of the shot points and receivers. Simultaneously, using the coordinates of the shot points and receivers, azimuth and offset information are calculated, and a 3D grid is added to the seismic data. Tomographic static correction was performed on the seismic data for this area, and the tomographic static correction parameters are applied. Figure 2 For the raw single-shot data containing first-to-multiple waves and surface waves that have undergone data preparation, spectral analysis of the surface waves in the data is performed. Figure 3 The spectrum analysis shows that the surface wave frequency is very high, reaching 20 Hz. The frequency of the surface wave overlaps with the frequency of the effective wave by a large margin. Therefore, if conventional regional filtering methods are used, the effective wave will inevitably be damaged. Figure 4 As can be seen from the frequency-wavenumber spectrum, the original data, due to its relatively large track spacing of 50 meters, does not show any linear relationship between the first arrival multiple waves and surface waves.

[0077] Step 102, Spatial Sampling Design. Spatial sampling design is performed to prepare for subsequent interpolation point locations. Based on the shot and receiver locations in the field observation system of this area, the new receiver spacing for pre-interpolation is designed to be 10 meters. The field observation system and the new receivers for pre-interpolation are then merged to form a pre-interpolation spatial model.

[0078] Step 103: Establish the estimated spectrum. The original seismic data is in the time-offset domain. It is converted to the frequency-wavenumber domain, and then the estimated spectrum is calculated for the single-frequency seismic data. The corresponding wave number is obtained based on the maximum Fourier spectrum. This wave number is then removed from the single-frequency seismic data to obtain a new single-frequency seismic data. The same method is used to obtain the wave number corresponding to the maximum Fourier spectrum for this new single-frequency seismic data. This process is repeated until the maximum Fourier spectrum is less than 0.05. All the obtained wave numbers are combined to obtain the estimated spectrum of a single-frequency seismic data. Finally, the estimated spectra of all single-frequency seismic data are combined to obtain the entire estimated spectrum.

[0079] Step 104, Multivariate Interpolation. Based on multivariate information such as the horizontal and vertical line numbers, time, offset, and azimuth of the seismic data, an independent wavenumber scan is performed in the frequency-wavenumber domain to determine the optimal interpolation wavenumber. Interpolation is then performed according to the pre-interpolation spatial model. The determination of the optimal interpolation wavenumber involves first determining the maximum Nyquist wavenumber based on the spatial sampling interval of the seismic data when interpolating at any location in the pre-interpolation spatial model. The spatial wavenumber is then continuously reduced at certain intervals, and the minimum root mean square error (RMSE) corresponding to each wavenumber is calculated. The maximum wavenumber that meets the set minimum RMSE threshold is selected, and then the interpolation weight coefficients corresponding to this wavenumber are calculated. Figure 5 The interpolated data generated using the multivariate interpolation method of this invention clearly shows a linear relationship between the first-to-first-order multiple waves and the surface waves. Further frequency-wavenumber spectrum analysis of this interpolated data is performed, such as... Figure 6 As shown, the frequency-wavenumber spectrum linear relationship is very clear at this time, indicating that the interpolation effect of the method of the present invention is ideal.

[0080] Step 105, interference wave removal. The seismic data interpolated in the previous step is in the frequency-wavenumber domain and has been interpolated into a well-sampled seismic data with a sound wavefield. Therefore, the linear interference relationship between the first arrival multiples and surface waves is clear. Based on this linear relationship, apparent velocity filtering can be performed to remove this interference wave. Figure 7 This is a single-shot data graph of noise suppression after multivariate interpolation according to the present invention. As can be seen from the graph, linear interference waves such as first-to-last, multiple, and surface waves are completely removed.

[0081] Step 106: Removal of new interpolated data. New interpolated seismic data are removed according to the pre-interpolated spatial model, retaining only the denoised original seismic data. The data then returns to the original 50-meter trace spacing, thus ensuring the fidelity of the seismic data.

[0082] Step 107, Discrete Fourier Transform. The noise-removed seismic data is inversely transformed back to the time-offset domain to obtain the final noise-suppressed data. The Discrete Fourier Transform can avoid transform aliasing caused by incomplete and irregular sampling of field data. Figure 8 This is the final single-shot data after processing using the multivariate interpolation noise suppression technique of this invention. Figure 9 The data is from a single gun using conventional noise reduction, and it can be seen that the present invention removes noise more effectively.

[0083] Finally, it should be noted that the above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

[0084] Except for the technical features described in the specification, all other technologies are known to those skilled in the art.

Claims

1. A noise suppression method using multivariate interpolation, characterized in that, The method includes: Step 101, data preparation, includes loading the seismic data acquired in the field into the observation system to obtain the coordinates, elevation, and trace distance information of the shot point and receiver point. At the same time, the azimuth and offset information are calculated using the coordinates of the shot point and receiver point. A three-dimensional grid is added to the seismic data. The seismic data is then subjected to tomographic static correction, and the tomographic static correction values ​​are applied. Step 102, Spatial sampling design: Based on the positions of shot points and receiver points in the field observation system of this area, design the new receiver point spacing for pre-interpolation to be 10 meters, and merge the field observation system with the new receiver points for pre-interpolation to form a pre-interpolation spatial model. Step 103: Establish the estimation spectrum. The original seismic data is in the time-offset domain. It is converted to the frequency-wavenumber domain. The estimation spectrum is then calculated for the single-frequency seismic data. The corresponding wave number is obtained based on the maximum Fourier spectrum. The wave number is then removed from the single-frequency seismic data to obtain a new single-frequency seismic data. The same method is used to obtain the wave number corresponding to the maximum Fourier spectrum for this new single-frequency seismic data. This process is repeated until the maximum Fourier spectrum is less than 0.

05. All the obtained wave numbers are combined to obtain the estimation spectrum of a single-frequency seismic data. The estimation spectra of all single-frequency seismic data are then combined to obtain the entire estimation spectrum. Step 104, Multivariate Interpolation: Based on the multivariate information of the seismic data, including the horizontal and vertical line numbers, time, offset, and azimuth, an independent wavenumber scan is performed in the frequency-wavenumber domain to determine the optimal interpolation wavenumber. Interpolation is then performed according to the pre-interpolation spatial model. When interpolating at any position in the pre-interpolation spatial model, the maximum Nyquist wavenumber is first determined based on the spatial sampling interval of the seismic data. The spatial wavenumber is then continuously reduced at certain intervals, and the minimum root mean square error corresponding to each of these wavenumbers is calculated. The maximum wavenumber that meets the set minimum root mean square error threshold is selected, and then the interpolation weight coefficients corresponding to that wavenumber are calculated. Step 105, interference wave removal: The seismic data interpolated in the previous step is in the frequency-wavenumber domain and has been interpolated into a well-sampled seismic data with a sound wavefield. Therefore, the linear interference relationship between the first arrival multiple waves and surface waves is clear. Based on this linear relationship, apparent velocity filtering is performed to remove this interference wave. Step 106: Removal of new interpolated data. The new interpolated seismic data is removed according to the pre-interpolated spatial model. Only the data after denoising the original seismic data is retained. The data returns to the original 50-meter trace spacing, thus ensuring the fidelity of the seismic data. Step 107, Discrete Fourier Transform, transforms the noise-removed seismic data back to the time-offset domain to obtain the final noise-suppressed data; Discrete Fourier Transform avoids transformation aliasing caused by incomplete and irregular sampling of field data.

Images