Method and device for multiple wave suppression of three-dimensional seismic data based on a rotating ellipse model

Abstract

The application relates to a three-dimensional seismic data multiple wave suppression method based on a rotating ellipse model, rotating ellipse model parameters with a characteristic of changing with a direction are added, and a spatially-varying ellipse model is used to correctly characterize a seismic wave field; first wave cutting is carried out in a transform domain, a multiple wave model is obtained through inverse transformation, and then the multiple wave model is subtracted to obtain a result after multiple wave suppression. In order to improve the resolution of the transform domain and better distinguish the first wave and the multiple wave in the transform domain, a low-frequency constraint method is used to improve the resolution, and the calculation result of the low-frequency data is used as a constraint matrix of the next high-frequency. Since the low-frequency data has no aliasing, the calculation result has higher space-time domain resolution, and is more conducive to the distinction between the first wave and the multiple wave and the suppression of the multiple wave. The method has better effect in the case that the underground medium is relatively complex and the wave field has anisotropy.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of geophysical exploration, in particular to a method and device for multiple wave suppression of three-dimensional seismic data based on a rotating elliptical model. BACKGROUND

[0002] With the vigorous development of oil and gas industry, three-dimensional seismic has become the main means of exploration. Compared with two-dimensional, three-dimensional exploration can more accurately and effectively suppress various interference waves, enhance effective waves, and more truly and delicately reflect the geological conditions underground, providing high-precision and high-quality processing results for structural lithology interpretation, reservoir inversion and oilfield development monitoring. Among them, wide-azimuth seismic exploration, due to its wide azimuth, can obtain more complete wave field information, which is beneficial to improve the imaging effect and is also more favored by people.

[0003] Compared with narrow-azimuth data, wide-azimuth seismic data has the characteristic of changing with the azimuth. Seismic traces with the same or similar offset but different azimuths may have large time difference due to the influence of different direction strata dip or strata anisotropy. When the dynamic corrected seismic data is arranged in order of increasing scalar offset to form a three-dimensional common midpoint gather similar to the two-dimensional case, the data of different azimuths will have a large time difference, and the time difference of adjacent traces will change dramatically, showing the characteristics of phase axis sample jumping. The conventional two-dimensional Radon algorithm cannot characterize this data, and cannot effectively separate the effective wave and the interference wave, so the suppression effect is poor. The three-dimensional common midpoint gather can be arranged according to the offset distance in x and y directions, i.e. the vector offset distance, to form a common midpoint gather composed of small gathers of different azimuths. The three-dimensional common midpoint gather composed of small gathers describes the propagation of three-dimensional wave field more intuitively. However, when the underground strata are inclined, the shape of the gather of different azimuths will have large differences, and the time difference of multiple waves of different azimuths will also have certain differences due to the influence of streamer drift and the difference of strata characteristics of different azimuths. Therefore, the two-dimensional transformation with the same parameters cannot achieve good results in multiple wave suppression of small gathers. Therefore, it is necessary to consider the propagation characteristics of three-dimensional seismic wave field, and use three-dimensional Radon theory to accurately characterize the data of different azimuths for better processing.

[0004] However, when the three-dimensional Radon transform is used to characterize the multi-azimuth gather seismic data, it is assumed that the characteristics of the medium in different directions are the same, the propagation characteristics of the three-dimensional wave field in different directions are also the same, and the time difference between the multiple waves and the effective waves in different azimuth gathers is also the same. The different azimuth gathers can be characterized by using the standard elliptical model because the time difference between the different azimuth gathers is the same. In the case that the dip angle of the underground structure is small and the fluctuation is gentle, the existing method can achieve good application effect; however, when the underground structure is complex, and under the influence of the drift of the acquisition streamer, the anisotropy of the stratum and the like, the different azimuth gathers have great difference, the existing method cannot accurately describe the three-dimensional propagation of the wave field, the characterization effect of the multiple waves in the transform domain is poor, the multiple waves transformed back have great difference with the original multiple waves in the data, the effective suppression of the multiple waves cannot be realized, and serious interference is caused to the later migration imaging and interpretation. SUMMARY

[0005] In view of the above problems, the purpose of the present application is to provide a three-dimensional seismic data multiple wave suppression method and device based on a rotating elliptical model, which can solve the problem that the traditional algorithm cannot cope with the azimuth anisotropy of the seismic data.

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

[0007] The three-dimensional seismic data multiple wave suppression method based on the rotating elliptical model comprises the following steps:

[0008] The time-space domain seismic data d t , x , y ) is transformed into the frequency domain by using the Fourier transform to obtain the frequency-space domain seismic data D ω ( ω , x , y );

[0009] Based on the frequency-space domain seismic data D ω ( ω , x , y ), the frequency-space domain seismic data vector D vec ( ω , G ) is calculated, and the operator L is calculated;

[0010] The first vector D vec (ω1, G) of the frequency-space domain seismic data vector is calculated by using the least square method to obtain the first vector M vec (ω1, K) of the frequency Radon domain seismic data, and the calculation result is saved in the weighted matrix W;

[0011] The other vectors M of the frequency Radon domain seismic data are calculated by using the low frequency constraint method on the other frequency spatial domain seismic data vectors vec ( ω i , K ), i = 2 ~ N, and the calculation result is saved in the weighted matrix W, so as to complete the calculation of all frequencies and obtain all vectors M of the frequency Radon domain seismic data vec ( ω , K );

[0012] Wherein, t is time; x is the in-line offset; y is the cross-line offset; ω is the angular frequency; G is the in-line and cross-line comprehensive shot-receiver factor after the matrix D ω is converted into a vector; K is the in-line and cross-line direction comprehensive curvature factor; ω1 is the angular frequency of the first vector; ω i is the angular frequency of the first vector; i is the angular frequency of the first vector; i is the vector number.

[0013] The three-dimensional seismic data multiple wave suppression method, preferably, the time-space domain seismic data d( t , x , y ) is transformed into the frequency domain by using Fourier transform to obtain the frequency spatial domain seismic data D ω ( ω , x , y ) The specific steps are as follows:

[0014] Based on the three-parameter three-dimensional Radon inverse transform of the rotating ellipse model, the time-space domain seismic data d( t , x , y ) can be written as:

[0015] (1)

[0016] The formula (1) corresponds to the discrete form of the frequency domain:

[0017] (2)

[0018] Wherein, q x is the in-line direction curvature parameter; qy The curvature parameter in the transverse survey line direction; q xy For the parameters of the rotated ellipse; t and All times are time; m For three-parameter Radon domain data in the spatiotemporal domain; M ω This is the three-parameter Radon domain data in the frequency spatial domain.

[0019] The aforementioned three-dimensional seismic data multiple suppression method preferably involves calculating the frequency spatial domain seismic data vector D. vec ( ω , G The operator L is specifically:

[0020] The frequency matrix D on both sides of formula (2) and M Written in vector form, denoted as D vec M vec ,Right now:

[0021] (3)

[0022] The exponent term in formula (2) Let L be the transformation operator, then equation (2) can be written in operator form:

[0023] (4)

[0024] Wherein, L is transformed by the longitudinal line transformation operator L x Transverse survey line transformation operator L y and rotation operator L qxy The Kronecker product of the three matrices is obtained as follows:

[0025] (5)

[0026] Among them, L x For longitudinal survey line transformation operator; L y For transverse survey line transformation operator; L qxy For rotation operators.

[0027] The three-dimensional seismic data multiple suppression method, preferably, involves suppressing the first vector D of the frequency spatial domain seismic data vector. vec ( ω 1, G), denoted as D vec_ω1 The first vector M of the frequency Radon domain seismic data is obtained by calculating using the least squares method. vec ( ω 1, K), which is denoted as M vec_ω1 Specifically:

[0028] (9)

[0029] and the calculation result is saved in the weighted matrix W through formula (10):

[0030] (10)

[0031] The three-dimensional seismic data multiple wave suppression method preferably uses a low-frequency constraint method to calculate other vectors M of the frequency Radon domain seismic data from other frequency spatial domain seismic data vectors vec ( ω i , K ), i = 2 ~ N, denoted as M vec_ωi , specifically:

[0032] (11)

[0033] and the calculation result M vec_ωi is saved in the weighted matrix W through formula (10).

[0034] The three-dimensional seismic data multiple wave suppression method preferably further comprises the following steps:

[0035] all vectors M vec ( ω , K ) of the frequency Radon domain seismic data are converted into a matrix M ω ( ω , q x , q y , q xy ) of the frequency Radon domain seismic data, and inverse Fourier transform is performed to obtain Radon data mr( t , q x , q y , q xy ) in the time-curvature domain;

[0036] a cut-off parameter qxcut is defined, and data of q x < qxcut is set to zero to obtain a multiple wave model mul_r( t , q x , q y ,q xy );

[0037] The time-curvature domain wave multiple model mul_r( t , q x , q y , q xy Perform a Fourier transform to the frequency domain to obtain the frequency-curvature domain multiple wave model matrix Mul_R. ω ( ω , q x , q y , q xy );

[0038] The multiple wave model matrix Mul_R in the frequency-curvature domain is calculated according to formula (3). ω ( ω , q x , q y , q xy The corresponding frequency-curvature domain multiple model vector Mul_R vec ( ω , K );

[0039] Based on operator L, the frequency-space domain multiple model vector mul_x is calculated using formulas (9) to (11). vec ( ω , G );

[0040] The frequency-space domain multiple model vector mul_x vec ( ω , G Transformed into a frequency-space domain multiple wave model matrix mul_x ω ( ω , x , y );

[0041] For the frequency space domain multiple wave model matrix mul_x ω ( ω , x , y Performing an inverse Fourier transform, we obtain the time-space domain multiple wave model mul_xt( t , x , y );

[0042] The time-space domain multiple model mul_xt( t , x , y ) is directly subtracted from the input data d( t , x , y ), to obtain the result after multiple suppression, i.e. primary wave p( t , x , y ).

[0043] The three-dimensional seismic data multiple suppression device based on the rotating elliptical model, comprises:

[0044] The first processing unit is used for transforming the time-space domain seismic data d( t , x , y ) into the frequency domain by using Fourier transform, to obtain the frequency-space domain seismic data D ω ( ω , x , y );

[0045] The second processing unit is based on the frequency-space domain seismic data D ω ( ω , x , y ), and is used for calculating the frequency-space domain seismic data vector D vec ( ω , G ) and calculating the operator L;

[0046] The third processing unit is used for calculating the first vector D vec (ω1, G) of the frequency-space domain seismic data vector, which is recorded as D vec_ω1 , by using the least square method, to obtain the first vector M vec (ω1, K) of the frequency Radon domain seismic data, which is recorded as M vec_ω1 , and the calculation result is saved in the weighted matrix W;

[0047] The fourth processing unit is used for calculating other vectors M vec ( ω i , K ) i = 2 ~ N, which is recorded as M vec_ωi , by using the low-frequency constraint method, and the calculation result is saved in the weighted matrix W, and the cycle is completed to calculate all frequencies, to obtain all vectors Mvec ω K

[0048] The computer storage medium of the application has a computer program stored thereon, and the computer program is executed by a processor to implement the method steps of the three-dimensional seismic data multiple wave suppression method based on the rotating elliptical model.

[0049] The computer equipment of the application comprises a memory, a processor and a computer program stored on the memory and executable on the processor, and the processor implements the method steps of the three-dimensional seismic data multiple wave suppression method based on the rotating elliptical model when executing the computer program.

[0050] The application has the following advantages due to the above technical scheme:

[0051] (1) better description of gathers with different azimuths and time difference changes;

[0052] (2) improved resolution of the transform domain by using a low-frequency constraint algorithm, which can improve the separation of three-dimensional Radon transform domain multiple waves and effective waves;

[0053] (3) deepening and expanding the three-dimensional Radon transform theory; at the same time, it can more thoroughly suppress multiple waves of complex regional seismic data, highlight effective signals, and improve the accuracy of seismic imaging, inversion and reservoir prediction. BRIEF DESCRIPTION OF DRAWINGS

[0054] Various other advantages and benefits will become apparent to those of ordinary skill in the art upon reading the following detailed description of the preferred embodiments. The accompanying drawings are included to provide a description of preferred embodiments, and are not meant to limit the present application. Throughout the drawings, the same reference designations are used to denote like elements. In the drawings:

[0055] Fig. 1 is the simulation data used for multiple wave suppression of the application; wherein (a) is a y-direction different facet gather; (b) is a scalar offset gather; (c) is a clearer display of the method effect;

[0056] Fig. 2 is the result display of the conventional three-dimensional least square algorithm for multiple wave suppression, wherein (a) is the estimated multiple wave facet gather of the method; (b) is the facet gather after multiple wave suppression; (c) is the scalar offset gather of the multiple wave suppression result; (d) is the estimated multiple wave among the three extracted small gathers; (e) is the facet gather after multiple wave suppression; (f) is the difference between the suppression result and the true primary wave;

[0057] ​​​Fig. 3 This invention provides the three-dimensional Radon transform result based on a rotated ellipse model.

[0058] Fig. 4 This is the result of the high-precision three-dimensional Radon transform based on the rotating ellipse model of the present invention. Detailed Implementation

[0059] Exemplary embodiments of the invention will now be described in more detail with reference to the accompanying drawings. While exemplary embodiments of the invention are shown in the drawings, it should be understood that the invention can be implemented in various forms and should not be limited to the embodiments set forth herein. Rather, these embodiments are provided to enable a more thorough understanding of the invention and to fully convey the scope of the invention to those skilled in the art.

[0060] This invention provides a method for suppressing multiple waves in 3D seismic data based on a rotating elliptical model. It adds parameters to the rotating elliptical model, which exhibits azimuth-dependent characteristics, and uses this spatially varying elliptical model to accurately represent the seismic wavefield. By performing primary wave shearing in the transform domain and inverse transforming to obtain the multiple wave model, the suppressed multiple wave model is then subtracted. To improve the resolution of the transform domain and better distinguish between primary and multiple waves, this invention employs a low-frequency constraint method. The low-frequency calculation results are used as the constraint matrix for the next high-frequency wave. Since low-frequency data is free of spurious frequencies, the resulting calculations have higher spatiotemporal resolution, which is more conducive to distinguishing between primary and multiple waves and suppressing multiple waves. This method performs better in cases with complex subsurface media and anisotropic wavefields.

[0061] The present invention provides a method for suppressing multiple waves of three-dimensional seismic data based on a rotated ellipse model, comprising the following steps:

[0062] S1. The spatiotemporal seismic data d( t , x , y The data is then transformed to the frequency domain using Fourier transform to obtain the frequency spatial domain seismic data D. ω ( ω , x , y );

[0063] S2. Based on frequency spatial domain seismic data D ω ( ω , x , y ), calculate the frequency spatial domain seismic data vector D vec ( ω , G ), and compute operator L;

[0064] S3. The first vector D of the frequency spatial domain seismic data vector. vec (ω1, G), denoted as D vec_ω1 The first vector M of the frequency Radon domain seismic data is obtained by calculating using the least squares method. vec (ω1, K), denoted as M vec_ω1 The calculation results are then stored in the weighting matrix W.

[0065] S4. Other frequency-domain seismic data vectors, such as M, are calculated using the low-frequency constraint method. vec ( ω i , K ), i = 2 ~ N, denoted as M vec_ωi The calculation results are then stored in the weighting matrix W. This process is repeated until all frequencies are calculated, yielding all vectors M of the frequency Radon domain seismic data. vec ( ω , K );

[0066] in, t For time; x y is the offset of the longitudinal survey line; y is the offset of the transverse survey line. ω Angular frequency; G To make matrix D ω After being converted into a vector, the combined shot-receiver distance factor for longitudinal and transverse survey lines; K ω1 is the combined curvature factor along the longitudinal and transverse survey lines; ω1 is the angular frequency of the first vector. ω i For the first i The angular frequency of a vector; i This is the vector index.

[0067] In the above embodiments, preferably, the step of transferring the spatiotemporal domain seismic data d( t , x , y The data is then transformed to the frequency domain using Fourier transform to obtain the frequency spatial domain seismic data D. ω ( ω , x , y The specific steps are as follows:

[0068] Based on the three-parameter three-dimensional inverse Radon transform of the rotating ellipse model, the spatiotemporal domain seismic data d( t , x , y It can be written as:

[0069] (1)

[0070] Formula (1) corresponds to the following frequency domain discrete form:

[0071] (2)

[0072] in, q x The curvature parameter in the longitudinal direction, q y The curvature parameter in the transverse direction. q xy For the parameters of the rotated ellipse; t and All values ​​are time; m represents the three-parameter Radon domain data in the spatiotemporal domain, M. ω This is the three-parameter Radon domain data in the frequency spatial domain.

[0073] In the above embodiments, preferably, the frequency spatial domain seismic data vector D is calculated. vec ( ω , G The operator L is specifically:

[0074] The frequency matrix D on both sides of formula (2) and M Written in vector form, denoted as D vec M vec ,Right now:

[0075] (3)

[0076] The exponent term in formula (2) Let L be the transformation operator. Then equation (2) can be written in operator form:

[0077] (4)

[0078] Wherein, L can be transformed by the longitudinal survey line transformation operator L x Transverse survey line transformation operator L y and rotation operator L qxy The Kronecker product of the three matrices is obtained as follows:

[0079] (5)

[0080] Among them, L x For longitudinal survey line transformation operator; L y For transverse survey line transformation operator; L qxy For rotation operators.

[0081] In the above embodiments, preferably, the first vector D of the frequency spatial domain seismic data vector...vec (ω1,G), denoted as D vec_ω1 , the first vector M of the frequency Radon domain seismic data is calculated by the least square method vec (ω1, K), denoted as M vec_ω1 , specifically:

[0082] (9)

[0083] and the calculation result is saved in the weighted matrix W through formula (10):

[0084] (10)

[0085] In the above embodiment, preferably, other frequency space domain seismic data vectors are calculated by using the low frequency constraint method to obtain other vectors M of the frequency Radon domain seismic data vec ( ω i , K ), i = 2 ~ N, denoted as M vec_ωi , specifically:

[0086] (11)

[0087] and the calculation result is saved in the weighted matrix W through formula (10):

[0088] In the above embodiment, preferably, the application further comprises the following steps:

[0089] S5. All vectors M of the frequency Radon domain seismic data vec ( ω , K ) are converted into a matrix M of the frequency Radon domain seismic data ω ( ω , q x , q y , q xy ), and inverse Fourier transform is performed to obtain Radon data mr( t , q x , q y , q xy ) in the time-curvature domain;

[0090] S6. A cut-off parameter qxcut is defined, and q ​x qxcut t q x q y q xy

[0091] S7. Fourier transform the time-curvature domain multiple model mul_r t q x q y q xy to the frequency domain to obtain the frequency-curvature domain multiple model matrix Mul_R ω ω q x q y q xy

[0092] S8. Calculate the frequency-curvature domain multiple model vector Mul_R ω ω q x q y q xy corresponding to the frequency-curvature domain multiple model matrix Mul_R vec ω K according to formula (3);

[0093] S9. Based on the operator L, calculate the frequency-space domain multiple model vector mul_x vec ω G using formula (9)~formula (11);

[0094] S10. Convert the frequency-space domain multiple model vector mul_x vec ω G to the frequency-space domain multiple model matrix mul_x ω ω x y ;

[0095] S11. Calculate the frequency-space domain multiple model matrix mul_x​​​​​​​​​​​​​​​​​​​​​​​​​​​​ω ω x y t x y

[0096] S12. The time-space domain multiple model mul_xt( t x y ) is directly subtracted from the input data d( t x y ), to obtain the result after multiple suppression, i.e., primary wave p( t x y ).

[0097] The application further provides a three-dimensional seismic data multiple suppression device based on a rotating elliptical model, comprising:

[0098] A first processing unit is configured to transform the time-space domain seismic data d( t x y ) into the frequency domain by using Fourier transform, to obtain frequency-space domain seismic data D ω ( ω x y ).

[0099] A second processing unit is configured to calculate the frequency-space domain seismic data vector D ω ( ω x y ) based on the frequency-space domain seismic data D vec ( ω G ), and calculate an operator L.

[0100] A third processing unit is configured to calculate the first vector D vec (ω1, G) of the frequency-space domain seismic data vector, which is recorded as D vec_ω1 , by using a least square method, to obtain the first vector M vec (ω1, K) of the frequency Radon domain seismic data, which is recorded as M vec_ω1 , and save the calculation result in a weighted matrix W.

[0101] A fourth processing unit is configured to calculate other vectors M​​​​​​​​​​​​​​​​​​​​vec ω i K i = 2 ~ N, denoted as M vec_ωi , and the calculation result is saved in the weighted matrix W, and the calculation is repeated to complete the calculation of all frequencies, to obtain all vectors M of the frequency Radon domain seismic data vec ω K

[0102] The application further provides a computer storage medium, which stores a computer program, and the computer program is executed by a processor to implement the method steps of the three-dimensional seismic data multiple wave suppression method based on a rotating ellipse model.

[0103] The application further provides a computer device, which comprises a memory, a processor, and a computer program stored in the memory and executable on the processor, and the processor implements the method steps of the three-dimensional seismic data multiple wave suppression method based on a rotating ellipse model when executing the computer program.

[0104] wherein, Fig. 1 is simulation data used for multiple wave suppression in the application; wherein (a) is a y-direction different facet gather; (b) is a scalar offset gather; (c) is three facet gathers extracted from (a) for clearer display of the method effect; it can be obviously seen from the simulation data that multiple waves of different directions have great differences.

[0105] Fig. 2 is a result display of multiple wave suppression by a conventional three-dimensional least square algorithm, wherein (a) is a multiple wave facet gather estimated by the method; (b) is a facet gather after multiple wave suppression, and it can be seen that, due to the fact that the method cannot cope with the azimuth anisotropy of seismic data, there is a great difference between the estimated multiple wave and the actual data, resulting in many residues of multiple wave suppression; (c) is a scalar offset gather of the multiple wave suppression result, and it can also be obviously seen that there are residues of multiple wave; (d) is the estimated multiple wave in the three extracted facet gathers, and it can be seen that there is a great difference between the multiple wave and the actual data, and there is also primary wave energy in it, which indicates that the conventional method cannot accurately estimate multiple wave and also causes damage to primary wave energy; (e) is a facet gather after multiple wave suppression, and it can be seen that there are a large number of residues of multiple wave after suppression; (f) is a difference between the suppression result and the true primary wave, and it can be obviously seen that there are multiple wave and primary wave energy, which indicates that the method is not complete and causes damage to primary wave.

[0106] Fig. 3 ​​​​​​The result is the result of the three-dimensional Radon transform based on the rotating ellipse model, i.e. the result obtained by using a least square algorithm similar to formula (9). It can be seen that compared with the traditional three-dimensional algorithm, the result is improved, but there are still some multiple wave residues and primary wave damage problems.

[0107] Fig. 4 The result of the high-precision three-dimensional Radon transform based on the rotating ellipse model of the application can obviously see that the method has very good effect, and the difference between the primary wave obtained after suppressing the multiple wave and the real primary wave, i.e. Fig. 4 f basically cannot see any energy, which fully illustrates the advantage of the method of the application.

[0108] Finally, it should be pointed out that: the above examples are only used to illustrate the technical solutions of the application, but not to limit them; although the application has been described in detail with reference to the foregoing examples, those skilled in the art should understand that: it can still modify the technical solutions recorded in the foregoing examples, or make equivalent replacement for part of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions of the embodiments of the application.

Claims

1. A method for suppressing multiple waves in three-dimensional seismic data based on a rotated ellipse model, characterized in that, Includes the following steps: S1: The spatiotemporal seismic data d( t , x , y The data is then transformed to the frequency domain using Fourier transform to obtain the frequency spatial domain seismic data D. ω ( ω , x , y ),in, t For time; x This is the offset distance of the longitudinal survey line; y This is the offset distance of the transverse survey line; ω To determine the angular frequency, the specific steps are as follows: Based on the three-parameter three-dimensional inverse Radon transform of the rotating ellipse model, the spatiotemporal domain seismic data d( t , x , y It can be written as: （1） Formula (1) corresponds to the following frequency domain discrete form: （2） in, q x The curvature parameter in the longitudinal direction; q y The curvature parameter in the transverse survey line direction; q xy For the parameters of the rotated ellipse; t and All times are time; m For three-parameter Radon domain data in the spatiotemporal domain; M ω The data are three-parameter Radon domain data in the frequency spatial domain; S2: Based on frequency spatial domain seismic data D ω ( ω , x , y ), calculate the frequency spatial domain seismic data vector D vec ( ω , G ), and compute operator L, where, G To make matrix D ω After being converted into a vector, the combined shot-receiver distance factor for longitudinal and transverse survey lines; Calculate the frequency spatial domain seismic data vector D vec ( ω , G The operator L is specifically: Rewrite the frequency matrices D and M on both sides of formula (2) in vector form, denoted as Di and Mm respectively. vec and M vec ,Right now: (3) The exponent term in formula (2) Let L be the transformation operator, then formula (2) can be written in operator form: (4) Wherein, L is transformed by the longitudinal line transformation operator L x Transverse survey line transformation operator L y and rotation operator L qxy The Kronecker product of the three matrices is obtained as follows: (5) S3: The first vector D of the frequency spatial domain seismic data vector vec (ω1, G) is calculated using the least squares method to obtain the first vector M of the frequency Radon domain seismic data. vec (ω1, K), and the calculation results are stored in the weighted matrix W, where, K ω1 is the combined curvature factor along the longitudinal and transverse survey lines; ω1 is the angular frequency of the first vector. ω i For the first i The angular frequency of a vector; i The vector index; S4: The low-frequency constraint method is used to calculate the other vectors M of the frequency Radon domain seismic data vectors. vec ( ω i , K ), i = 2 ~ N, and store the calculation results in the weighting matrix W. Repeat this process to complete the calculation for all frequencies, thus obtaining all vectors M of the frequency Radon domain seismic data. vec ( ω , K ); S5: M of all vectors from the frequency Radon domain seismic data vec ( ω , K The matrix M transformed into frequency Radon domain seismic data ω ( ω , q x , q y , q xy ), and perform an inverse Fourier transform to obtain Radon data mr in the time-curvature domain. t , q x , q y , q xy ); S6: Define resection parameters qxcut and will q x < qxcut By setting the data to zero, we obtain the time-curvature domain multiple wave model mul_r( t , q x , q y , q xy ); S7: The time-curvature domain multiple wave model mul_r( t , q x , q y , q xy Perform a Fourier transform to the frequency domain to obtain the frequency-curvature domain multiple wave model matrix Mul_R. ω ( ω , q x , q y , q xy ); S8: Calculate the multiple wave model matrix Mul_R in the frequency-curvature domain according to formula (3). ω ( ω , q x , q y , q xy The corresponding frequency-curvature domain multiple model vector Mul_R vec ( ω , K ); S9: Based on operator L, the frequency-space domain multiple model vector mul_x is calculated using formulas (9) to (11). vec ( ω , G ); S10: Convert the frequency-space domain multiple wave model vector mul_x vec ( ω , G Transformed into a frequency-space domain multiple wave model matrix mul_x ω ( ω , x , y ); S11: For the frequency space domain multiple wave model matrix mul_x ω ( ω , x , y Performing an inverse Fourier transform, we obtain the time-space domain multiple wave model mul_xt( t , x , y ); S12: The time-space domain multiple wave model mul_xt( t , x , y From input data d( t , x , y Subtracting directly yields the result after multiple wave suppression, i.e., the first wave p( t , x , y ).

2. A computer storage medium having a computer program stored thereon, characterized in that, When the computer program is executed by the processor, it implements the method steps of the three-dimensional seismic data multiple suppression method based on the rotating ellipse model as described in claim 1.

3. A computer device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the computer program, it implements the method steps of the three-dimensional seismic data multiple suppression method based on the rotating ellipse model as described in claim 1.

Images