A multiple wave prediction method, suppression method, related device and equipment

Abstract

The present application relates to a kind of multiple wave prediction method, suppress method, related device and equipment.The prediction method includes: the prestack common receiver gathers are carried out Fourier positive transformation, obtain the transformed seismic trace set of three-dimensional frequency wavenumber domain;According to preset phase shift method, wave field continuation is carried out to transformed seismic trace set, obtain downgoing wave field;According to the migration velocity field obtained, the wave field weighting processing is carried out to each point in the continuation grid corresponding to the downgoing wave field of any depth in transverse direction, the interpolation downgoing wave field obtained is carried out Fourier inverse transformation, according to the seismic data of space-time domain obtained and the poststack depth domain migration profile obtained in advance, determine the seismic data excited by secondary source;The seismic data excited by secondary source is carried out Fourier transformation, and the three-dimensional frequency wavenumber domain wave field obtained is carried out wave field continuation according to preset phase shift method, obtain upgoing wave field;The upgoing wave field is carried out Fourier inverse transformation, and the multiple wave model data of space-time domain is obtained.

Description

Technical Field

[0001] This application belongs to the field of seismic data processing technology in oil and gas exploration, and in particular to a multiple wave prediction method, suppression method, related devices and equipment. Background Technology

[0002] With the increasing complexity of marine exploration, conventional narrow-azimuth towed cable seismic acquisition data can no longer meet the technical requirements for high-quality imaging. Ocean Bottom Nodes (OBNs), as an emerging seismic data acquisition method, represent an important development trend in marine seismic exploration due to their omnidirectional, multi-component, high-density, and high signal-to-noise ratio characteristics. However, the unique characteristics of OBN acquisition systems (i.e., the receiver is on the seabed, and the excitation point is on the sea surface) present numerous challenges for OBN seismic data processing. Multiples, being the most significant source of noise in OBN seismic data, must be eliminated. Improper processing of multiples in OBN seismic data can easily confuse and interfere with primary wave energy, severely impacting the final image quality. Previously, multiples in marine seismic data were generally suppressed using the Surface-Related Multiple Elimination (SRME) method. However, due to the unique characteristics of OBN acquisition systems—the receiver and excitation point are not on the same plane, and the receivers deployed on the seabed are relatively sparse—the assumptions of the SRME method—that the excitation and receiver points are on the same plane and in the same location—are no longer met. Currently, there is an urgent need for a new multiple wave suppression technology to solve the surface multiple wave problem in ocean OBN data. Summary of the Invention

[0003] In view of this, the present invention aims to provide a multiple prediction method, suppression method, related apparatus and equipment, which can solve the problem that the traditional SRME method cannot achieve multiple prediction and suppression in OBN seismic data, overcome the technical obstacles that the excitation point and receiver point of OBN seismic data are not in the same plane and the distribution of seabed nodes is sparse, and can predict surface multiples of all orders at once.

[0004] In a first aspect, this application provides a method for predicting multiple waves, including:

[0005] The pre-stack common receiver point seismic gathers were subjected to a positive Fourier transform to obtain the transformed seismic gathers in the three-dimensional frequency wavenumber domain.

[0006] The transformed seismic gathers are subjected to wavefield extension using a preset phase-shift method to obtain the downlink wavefield;

[0007] Based on the acquired offset velocity field, wave field weighting is performed on each point in the transverse grid corresponding to the downlink wave field at any depth to obtain the interpolated downlink wave field.

[0008] The interpolated downflow wave field is subjected to inverse Fourier transform to obtain seismic data in the spatiotemporal domain.

[0009] Based on the spatiotemporal seismic data and the pre-acquired post-stack depth domain migration profile, the seismic data excited by the secondary source are determined.

[0010] The seismic data excited by the secondary source are subjected to Fourier transform to obtain the wavefield in the three-dimensional frequency-wavenumber domain.

[0011] The three-dimensional frequency wavenumber domain wavefield is extended in the frequency wavenumber domain according to a preset phase shift method to obtain the upward wavefield.

[0012] The upward wave field is subjected to inverse Fourier transform to obtain multiple wave model data in the spatiotemporal domain.

[0013] In conjunction with the first aspect above, in one possible implementation, the step of performing wavefield weighting processing on each point in the transverse grid corresponding to the downlink wavefield at any depth based on the acquired offset velocity field to obtain the interpolated downlink wavefield includes:

[0014] Based on the acquired offset velocity field, the range of lateral velocity variation of the downflow wave field at any depth is determined, and the range of lateral velocity variation is divided into multiple velocity intervals.

[0015] Wavefield extension is performed using the endpoints of the multiple velocity intervals to obtain wavefield values ​​corresponding to the endpoints of the multiple velocity intervals, thus obtaining multiple wavefield intervals corresponding to the multiple velocity intervals;

[0016] The wave field interval corresponding to the velocity interval of each point in the horizontal direction of the extended grid is determined, and wave field weighting is performed based on the two endpoints of the wave field interval to obtain the interpolated downlink wave field.

[0017] In conjunction with the first aspect above, in one possible implementation, the step of determining the lateral velocity variation range of the downlink wavefield at any depth based on the acquired offset velocity field, and dividing the lateral velocity variation range into multiple velocity intervals, includes:

[0018] At any depth, the range of lateral velocity variation of the downlink wave field is obtained based on the lateral and longitudinal coordinates in the offset velocity field.

[0019] The velocity within the lateral velocity variation range is divided into multiple velocity end values ​​at equal intervals;

[0020] The range of speeds between any two adjacent speed values ​​is considered as a speed interval, resulting in multiple speed intervals.

[0021] In conjunction with the first aspect above, in one possible implementation, the step of performing wavefield extrapolation using the endpoints of the plurality of velocity intervals to obtain wavefield endpoints corresponding to the endpoints of the plurality of velocity intervals, and obtaining a plurality of wavefield intervals corresponding to the plurality of velocity intervals, includes:

[0022] Wavefield extension is performed using the multiple velocity endpoints to obtain wavefield endpoints corresponding to each velocity endpoint. The wavefield range between each two adjacent wavefield endpoints is taken as a wavefield interval, resulting in multiple wavefield intervals corresponding to the multiple velocity intervals.

[0023] In conjunction with the first aspect above, in one possible implementation, determining the wavefield interval corresponding to the velocity interval to which each point in the extended grid belongs, and performing wavefield weighting based on the two endpoints of the wavefield interval to obtain the interpolated downlink wavefield, includes:

[0024] Determine the velocity range in which the velocity value of each point in the horizontal direction of the extended grid lies;

[0025] Based on the correspondence between the multiple velocity intervals and the multiple wavefield intervals, the wavefield interval corresponding to the velocity value at each point in the lateral direction is determined. Wavefield weighting is then performed based on the following formula to determine the wavefield value at each point in the lateral direction, resulting in the interpolated downlink wavefield:

[0026] ;

[0027] In the formula, and Let be the endpoint value of the i-th wavefield interval. , These are the weighting coefficients.

[0028] In conjunction with the first aspect described above, in one possible implementation, the prediction method further includes obtaining the weighting coefficients in the following manner. and :

[0029] The weighting coefficients are determined based on the coordinate distance between each point in the horizontal direction of the extended grid and the endpoint of the wave field interval. and .

[0030] In conjunction with the first aspect above, in one possible implementation, the step of performing a positive Fourier transform on the acquired pre-stack common receiver point seismic gathers to obtain the transformed seismic gathers in the three-dimensional frequency-wavenumber domain includes:

[0031] Perform a one-dimensional fast Fourier transform along the time direction on each trace in the pre-stack common receiver point seismic trace set to obtain the frequency-spatial domain seismic trace set.

[0032] Perform a two-dimensional fast Fourier transform on each trace in the frequency-space domain seismic trace set along the lateral direction to obtain a two-dimensional frequency-wavenumber domain seismic trace set.

[0033] Perform a three-dimensional fast Fourier transform on each trace in the two-dimensional frequency-wavenumber domain seismic gather along the longitudinal direction to obtain the transformed seismic gather in the three-dimensional frequency-wavenumber domain.

[0034] In conjunction with the first aspect above, in one possible implementation, the step of performing wavefield extension on the transformed seismic gather according to a preset phase shift method to obtain the downlink wavefield includes:

[0035] Substituting the transformed seismic gathers into the following formula for downward wavefield extension, we obtain the downward wavefield:

[0036] ;

[0037] in, Indicates the downlink wave field. Represents the transformed seismic gathers; The angular frequency of the seismic data. for Wave number in direction, for Wave number in direction, These are the depth coordinates of the underground strata. This represents the step size of the seismic wavefield extension, where, , Indicates the speed of seawater. This indicates the maximum frequency of the pre-stack common receiver point seismic gather; Denotes the continuation operator for downward propagation, where This represents the wave number when the wave field extends downwards to a depth of z. In the formula, the symbol "+" indicates that the wave field propagates downwards. In the offset velocity field Velocity value at depth position.

[0038] In conjunction with the first aspect above, in one possible implementation, performing an inverse Fourier transform on the interpolated downlink wavefield to obtain seismic data in the spatiotemporal domain includes:

[0039] The interpolated downlink wave field Along Perform an inverse fast Fourier transform on the direction to obtain wavefield data in the frequency and wavenumber domains. ;

[0040] The wavefield data in the frequency wavenumber domain Along The direction is subjected to inverse fast Fourier transform to obtain seismic data in the frequency spatial domain. ;

[0041] The seismic data in the frequency spatial domain Performing a fast inverse Fourier transform along the frequency direction yields seismic data in the spatiotemporal domain. .

[0042] In conjunction with the first aspect above, in one possible implementation, determining the seismic data excited by the secondary source based on the spatiotemporal seismic data and the pre-acquired post-stack depth-domain migration profile includes:

[0043] Based on the spatiotemporal seismic data and the pre-acquired post-stack depth-domain migration profile, the seismic data excited by the secondary source are obtained according to the following formula:

[0044] ;

[0045] in, This represents the post-stack depth domain offset profile. This represents seismic data in the spatiotemporal domain.

[0046] In conjunction with the first aspect above, in one possible implementation, the step of performing wavefield extension on the three-dimensional frequency-wavenumber domain wavefield in the frequency-wavenumber domain according to a preset phase-shift method to obtain an upward wavefield includes:

[0047] Substituting the three-dimensional frequency-wavenumber domain wavefield into the following formula for upward wavefield extension, we obtain the upward wavefield:

[0048]

[0049] in, Indicates the upward wave field. Represents the three-dimensional frequency-wavenumber domain wavefield; The angular frequency of the seismic data. for Wave number in direction, for Wave number in direction, These are the depth coordinates of the underground strata. Indicates the step size of the seismic wavefield extension; This represents the extension operator that propagates upwards.

[0050] In conjunction with the first aspect above, in one possible implementation, performing an inverse Fourier transform on the ascending wave field to obtain multiple wave model data in the spatiotemporal domain includes:

[0051] Determine the total multiple wave field when the upward wave field ascends to a depth of 0;

[0052] Performing an inverse Fourier transform on the total multiple wave field yields the result at a depth of [missing information]. Multiple wave model data generated at the location .

[0053] In conjunction with the first aspect described above, in one possible implementation, the prediction method further includes: obtaining pre-stack common receiver point seismic gathers in the following manner:

[0054] If the seismic data acquired at sea is OBN acquired seismic data, the shot-receiver points of the OBN acquired seismic data are interchanged, and the pre-stack common receiver point gather is extracted from the seismic data obtained by the interchange.

[0055] If the seismic data acquired at sea is acquired by a towed cable, the pre-stack common receiver gather is extracted from the towed cable acquired seismic data.

[0056] Secondly, this application provides a method for suppressing multiple waves, including:

[0057] The pre-stack common receiver gathers are matched and subtracted from the multiple wave model data determined by the multiple wave prediction method to obtain the seismic data after suppressing multiple waves.

[0058] Thirdly, this application provides a multiple wave prediction device, comprising:

[0059] The first transformation module is used to perform a forward Fourier transform on the acquired pre-stack common receiver point seismic gathers to obtain the transformed seismic gathers in the three-dimensional frequency wavenumber domain.

[0060] The first extension module is used to extend the wavefield of the transformed seismic gather according to a preset phase shift method to obtain the downlink wavefield.

[0061] The interpolation module is used to perform wavefield weighting on each point in the transverse grid corresponding to the downlink wavefield at any depth based on the acquired offset velocity field, so as to obtain the interpolated downlink wavefield.

[0062] The second transformation module is used to perform an inverse Fourier transform on the interpolated downlink wave field to obtain seismic data in the spatiotemporal domain.

[0063] The determination module is used to determine the seismic data excited by the secondary source based on the seismic data in the spatiotemporal domain and the pre-acquired post-stack depth domain migration profile.

[0064] The third transformation module is used to perform Fourier transform on the seismic data excited by the secondary source to obtain the wavefield in the three-dimensional frequency wavenumber domain.

[0065] The second extension module is used to extend the three-dimensional frequency wavenumber domain wavefield in the frequency wavenumber domain according to a preset phase shift method to obtain an upward wavefield.

[0066] The fourth transformation module is used to perform an inverse Fourier transform on the ascending wave field to obtain multiple wave model data in the spatiotemporal domain.

[0067] Fourthly, this application provides a multiple wave suppression device, comprising:

[0068] The first transformation module is used to perform a forward Fourier transform on the acquired pre-stack common receiver point seismic gathers to obtain the transformed seismic gathers in the three-dimensional frequency wavenumber domain.

[0069] The first extension module is used to extend the wavefield of the transformed seismic gather according to a preset phase shift method to obtain the downlink wavefield.

[0070] The interpolation module is used to perform wavefield weighting on each point in the transverse grid corresponding to the downlink wavefield at any depth based on the acquired offset velocity field, so as to obtain the interpolated downlink wavefield.

[0071] The second transformation module is used to perform an inverse Fourier transform on the interpolated downlink wave field to obtain seismic data in the spatiotemporal domain.

[0072] The determination module is used to determine the seismic data excited by the secondary source based on the seismic data in the spatiotemporal domain and the pre-acquired post-stack depth domain migration profile.

[0073] The third transformation module is used to perform Fourier transform on the seismic data excited by the secondary source to obtain the wavefield in the three-dimensional frequency wavenumber domain.

[0074] The second extension module is used to extend the three-dimensional frequency wavenumber domain wavefield in the frequency wavenumber domain according to a preset phase shift method to obtain an upward wavefield.

[0075] The fourth transformation module is used to perform an inverse Fourier transform on the upward wave field to obtain multiple wave model data in the spatiotemporal domain.

[0076] The suppression module is used to match and subtract the pre-stack common receiver gather from the multiple wave model data to obtain the suppressed multiple wave seismic data.

[0077] Fifthly, this application provides a computer-readable storage medium storing instructions that, when executed on a terminal, cause the terminal to perform the multiple wave prediction method as described in the first aspect or the multiple wave suppression method as described in the second aspect.

[0078] Sixthly, this application provides a computer device, characterized in that it includes a processor, a communication interface, a memory, and a communication bus, wherein the processor, the communication interface, and the memory communicate with each other through the communication bus;

[0079] Memory, used to store computer programs;

[0080] When executing a program stored in memory, the processor implements the multiple wave prediction method as described in the first aspect or the multiple wave suppression method as described in the second aspect.

[0081] In a seventh aspect, this application provides a computer program product containing instructions that, when run on a computer device, cause the computer device to perform the multiple wave prediction method as described in the first aspect or the multiple wave suppression method as described in the second aspect.

[0082] Eighthly, this application provides a chip including a processor and a communication interface, the communication interface and the processor being coupled together, the processor being used to run computer programs or instructions to implement the multiple wave prediction method as described in the first aspect or the multiple wave suppression method as described in the second aspect.

[0083] Specifically, the chip provided in this application also includes a memory for storing computer programs or instructions.

[0084] It should be noted that the aforementioned computer instructions may be stored, in whole or in part, on a computer-readable storage medium. This computer-readable storage medium may be packaged together with the processor of the device, or it may be packaged separately from the processor of the device; this application does not impose any limitation on this.

[0085] Due to the adoption of the above technical solution, the beneficial effects of the present invention are:

[0086] This invention provides a multiple wave prediction method. Compared with the conventional SRME method, it is applicable to multiple wave prediction of marine seismic data acquired when the excitation point and receiver point are not in the same plane and the seabed nodes are sparsely distributed. Therefore, it can realize multiple wave prediction of both OBN seismic data and marine towed seismic data, thereby enabling the suppression of surface multiple waves in marine seismic data. It has a wide range of applications, strong versatility, and practical application value.

[0087] This invention provides a multiple wave prediction method. It employs a phase-shift wavefield extension approach to propagate seismic data downwards. Based on a pre-acquired post-stack depth-domain migration profile, it identifies secondary sources formed at strong-reflecting interfaces during downward propagation. When the secondary source wavefield propagates upwards, it captures the multiple waves of various orders formed by reflections between the sea surface and strong-reflecting interfaces, enabling the prediction of all multiple wave orders in a single operation. Furthermore, to overcome the limitation of phase-shift wavefield extension only propagating in horizontally layered media, a wavefield interpolation weighting method is introduced. Based on the acquired migration velocity field, wavefield weighting is applied to each point laterally in the extension grid corresponding to the downflow wavefield at any depth, resulting in the interpolated downflow wavefield. This effectively solves the problem of the immutable lateral velocity field and is applicable to multiple wave prediction of marine-acquired seismic data in complex underground strata, improving the accuracy of multiple wave prediction.

[0088] This invention provides a multiple wave prediction method that performs multiple wave prediction on acquired pre-stack common receiver point seismic gathers, realizing single-gather multiple wave prediction processing. Multiple wave prediction is achieved in the frequency wavenumber domain. Since frequency domain seismic data is decoupled and each frequency is independent of each other, multi-core parallel computing can be realized, resulting in fast computing speed and improved computing efficiency. Moreover, it does not require a large amount of memory space and is suitable for multiple wave prediction processing of large amounts of marine-acquired seismic data.

[0089] This invention provides a multiple wave prediction method. Compared with conventional multiple wave prediction methods, it does not require picking up the layer information of multiple wave generation during the multiple wave prediction process. It can predict the multiple wave information of all layers contained in a given profile, which is convenient for multiple wave prediction processing of large amounts of marine seismic data. Attached Figure Description

[0090] Figure 1 This is a flowchart illustrating the multiple wave prediction method provided in an embodiment of the present invention;

[0091] Figure 2 This is a simplified schematic diagram of the offset and reverse offset process provided in the embodiments of the present invention;

[0092] Figure 3 This is a schematic diagram illustrating the interchangeability of shot-receiver relationships in OBN seismic data provided in an embodiment of the present invention;

[0093] Figure 4 A schematic diagram of wavefield extension of a pre-stack common receiver point gather provided in an embodiment of the present invention;

[0094] Figure 5 This is a schematic diagram of lateral velocity interpolation in the velocity field provided in an embodiment of the present invention;

[0095] Figure 6A schematic diagram of a three-dimensional pre-stack common receiver gather for OBN seismic data suppression of multiple wavefronts provided in an embodiment of the present invention;

[0096] Figure 7 To Figure 6 The diagram shows the multiple wave model data obtained by multiple wave prediction from the three-dimensional pre-stack common receiver gather.

[0097] Figure 8 To Figure 6 The diagram shows a seismic gather obtained by multiple wave suppression of a three-dimensional pre-stack common receiver gather.

[0098] Figure 9 A schematic diagram of a two-dimensional pre-stack common receiver gather for OBN seismic data suppression of multiple wavefronts provided in an embodiment of the present invention;

[0099] Figure 10 To Figure 9 The diagram shows the multiple wave model data obtained by multiple wave prediction from a two-dimensional pre-stack common receiver gather.

[0100] Figure 11 To Figure 9 The diagram shows a seismic gather obtained by multiple wave suppression of a two-dimensional pre-stack common receiver gather.

[0101] Figure 12 To Figure 9 The diagram shows a seismic gather obtained by forward modeling a two-dimensional pre-stack common receiver gather.

[0102] Figure 13 A schematic flowchart of the multiple wave suppression method provided in an embodiment of the present invention;

[0103] Figure 14 This is a schematic diagram of the structure of the multiple wave prediction device provided in an embodiment of the present invention;

[0104] Figure 15 This is a schematic diagram of the structure of the multiple wave suppression device provided in an embodiment of the present invention;

[0105] Figure 16 A structural block diagram of a computer device provided in an embodiment of the present invention. Detailed Implementation

[0106] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0107] In this article, the term "and / or" is merely a description of the relationship between related objects, indicating that there can be three relationships. For example, A and / or B can represent three situations: A exists alone, A and B exist simultaneously, and B exists alone.

[0108] The terms "first" and "second," etc., used in the specification and drawings of this application are used to distinguish different objects or to distinguish different treatments of the same object, rather than to describe a specific order of objects.

[0109] Furthermore, the terms "comprising" and "having," and any variations thereof, used in the description of this application are intended to cover non-exclusive inclusion. For example, a process, method, system, product, or apparatus that includes a series of steps or units is not limited to the steps or units listed, but may optionally include other steps or units not listed, or may optionally include other steps or units inherent to such process, method, product, or apparatus.

[0110] It should be noted that in the embodiments of this application, the words "exemplary" or "for example" are used to indicate examples, illustrations, or explanations. Any embodiment or design scheme described as "exemplary" or "for example" in the embodiments of this application should not be construed as being more preferred or advantageous than other embodiments or design schemes. Specifically, the use of the words "exemplary" or "for example" is intended to present the relevant concepts in a specific manner.

[0111] In the description of this application, unless otherwise stated, "a plurality of" means two or more.

[0112] The embodiments of this application will now be described in detail with reference to the accompanying drawings.

[0113] This application provides a method for predicting multiple waves, referring to... Figure 1 As shown, the method includes:

[0114] S101: Perform a positive Fourier transform on the acquired pre-stack common receiver point seismic gathers to obtain the transformed seismic gathers in the three-dimensional frequency wavenumber domain.

[0115] S102: The transformed seismic gather is subjected to wavefield extension according to a preset phase shift method to obtain the downlink wavefield;

[0116] S103: Based on the acquired offset velocity field, perform wave field weighting on each point in the transverse grid corresponding to the downlink wave field at any depth to obtain the interpolated downlink wave field.

[0117] S104: Perform an inverse Fourier transform on the interpolated downlink wave field to obtain seismic data in the spatiotemporal domain;

[0118] S105: Based on the spatiotemporal seismic data and the pre-acquired post-stack depth domain migration profile, determine the seismic data triggered by the secondary source.

[0119] S106: Perform Fourier transform on the seismic data excited by the secondary source to obtain the wavefield in the three-dimensional frequency-wavenumber domain.

[0120] S107: The three-dimensional frequency wavenumber domain wavefield is extended in the frequency wavenumber domain according to a preset phase shift method to obtain the up-going wavefield.

[0121] S108: Perform an inverse Fourier transform on the ascending wave field to obtain multiple wave model data in the spatiotemporal domain.

[0122] The multiple wave prediction method provided in this embodiment of the invention can predict multiple waves of both OBN seismic data and offshore towed seismic data. If the offshore seismic data is OBN-acquired seismic data, the pre-stack common receiver gather described in step S101 is extracted from the OBN-acquired seismic data by exchanging the shot and receiver points. If the offshore seismic data is towed seismic data, the pre-stack common receiver gather described in step S101 is directly extracted from the towed seismic data.

[0123] Reference Figure 2 As shown in the embodiment of the present invention, by analyzing the acquired pre-stack common receiver point seismic gathers... Fourier transform and phase-shift wavefield extension are performed based on the obtained offset velocity field. and post-stack depth domain offset profile The method implements inverse offset operation to predict the corresponding multiple waves. In the multiple wave prediction process, a velocity interpolation weighting method is introduced, which can predict surface multiple waves of all orders at once, and expands the application scope of the multiple wave prediction method.

[0124] In this embodiment of the invention, for OBN-acquired seismic data, since the receiver point for OBN-acquired seismic data is located on the seabed while the shot point is located on the sea surface, in order to improve the accuracy of multiple wave prediction and computational efficiency, and to facilitate subsequent wavefield extrapolation, refer to... Figure 3 As shown, the positions of the shot point and the receiver point can be swapped according to the shot-receiver interchange principle. That is, the receiver point receives at the sea surface, and the shot point excites at the seabed. The pre-stack common receiver point gather can then be extracted from the swapped seismic data. .in, and These represent the horizontal and vertical coordinates within the seismic acquisition area, respectively. These are the travel time coordinates of earthquake propagation. This is due to the pre-stack common receiver gather. It can be either 3D or 2D seismic data; therefore, if the pre-stack common receiver point gather... If the data is two-dimensional seismic data, then its y-direction coordinate is 0.

[0125] In an optional embodiment, in step S101 above, the acquired pre-stack common receiver point seismic gathers are subjected to a positive Fourier transform to obtain the transformed seismic gathers in the three-dimensional frequency-wavenumber domain. Specifically, it includes:

[0126] The pre-stack common receiver point seismic gather Performing a one-dimensional fast Fourier transform along the time direction on each trace yields the seismic gathers in the frequency-spatial domain. ,in, The angular frequency of the seismic data;

[0127] The frequency-spatial domain seismic gathers Each trace in the dataset undergoes a two-dimensional fast Fourier transform along the lateral direction (x-direction) to obtain a two-dimensional frequency-wavenumber domain seismic gather. ,in, for Wave number in direction;

[0128] The seismic gathers in the two-dimensional frequency-wavenumber domain Each trace in the dataset undergoes a three-dimensional fast Fourier transform along the longitudinal direction (y-direction) to obtain the transformed seismic gathers in the three-dimensional frequency-wavenumber domain. ,in, for Wave number in direction.

[0129] Of course, it should be noted here that if the pre-stack common receiving point gather Since the data is two-dimensional seismic data, and its y-axis coordinate is 0, the step of performing a forward Fourier transform in the longitudinal direction is not required.

[0130] In an optional embodiment, the step S102 above, which involves performing wavefield extension on the transformed seismic gather using a preset phase-shift method to obtain the downlink wavefield, specifically includes:

[0131] Substituting the transformed seismic gathers into the following formula for downward wavefield extension, we obtain the downward wavefield:

[0132] , formula 1;

[0133] in, Indicates the downlink wave field. This represents the wavefield of the transformed seismic gather at depth z; The angular frequency of the seismic data. for Wave number in direction, for Wave number in direction, These are the depth coordinates of the underground strata. This represents the step size of the seismic wavefield extension, where, , Indicates the speed of seawater. This indicates the maximum frequency of the pre-stack common receiver point seismic gather; Denotes the continuation operator for downward propagation, where This represents the wave number when the wave field extends downwards to a depth of z. In the formula, the symbol "+" indicates that the wave field propagates downwards. In the offset velocity field Velocity value at depth position.

[0134] In this embodiment of the invention, due to the seismic gathers in the three-dimensional frequency wavenumber domain... The shot point and receiver point were interchanged; therefore, it can be understood that the wave field changes from the sea surface position according to... The step size is propagated downwards layer by layer, such as Figure 4 As shown, this represents the path of the seismic data received from the sea surface as it propagated downwards again.

[0135] In an optional embodiment, step S103 above involves performing wavefield weighting on each point in the transverse grid corresponding to the downlink wavefield at any depth based on the acquired offset velocity field, to obtain the interpolated downlink wavefield. Specifically, this includes:

[0136] Based on the acquired offset velocity field, the range of lateral velocity variation of the downflow wave field at any depth is determined, and the range of lateral velocity variation is divided into multiple velocity intervals.

[0137] Wavefield extension is performed using the endpoints of the multiple velocity intervals to obtain wavefield values ​​corresponding to the endpoints of the multiple velocity intervals, thus obtaining multiple wavefield intervals corresponding to the multiple velocity intervals;

[0138] The wave field interval corresponding to the velocity interval of each point in the horizontal direction of the extended grid is determined, and wave field weighting is performed based on the two endpoints of the wave field interval to obtain the interpolated downlink wave field.

[0139] In an optional embodiment, step S104 above determines the lateral velocity variation range of the downlink wave field at any depth based on the acquired offset velocity field, and divides the lateral velocity variation range into multiple velocity intervals, specifically including:

[0140] At any depth, the range of lateral velocity variation of the downlink wave field is obtained based on the lateral and longitudinal coordinates in the offset velocity field.

[0141] The velocity within the lateral velocity variation range is divided into multiple velocity end values ​​at equal intervals;

[0142] The range of speeds between any two adjacent speed values ​​is considered as a speed interval, resulting in multiple speed intervals.

[0143] In one specific embodiment, the wavefield is extended using the endpoints of the plurality of velocity intervals to obtain wavefield endpoints corresponding to the endpoints of the plurality of velocity intervals, thereby obtaining a plurality of wavefield intervals corresponding to the plurality of velocity intervals. Specifically, this includes:

[0144] Wavefield extension is performed using the multiple velocity endpoints to obtain wavefield endpoints corresponding to each velocity endpoint. The wavefield range between each two adjacent wavefield endpoints is taken as a wavefield interval, resulting in multiple wavefield intervals corresponding to the multiple velocity intervals.

[0145] In one specific embodiment, the wavefield interval corresponding to the velocity interval to which each point in the transverse direction of the extended grid belongs is determined. Based on the two endpoints of the wavefield interval, wavefield weighting is performed to obtain the interpolated downlink wavefield, specifically including:

[0146] Determine the velocity range in which the velocity value of each point in the horizontal direction of the extended grid lies;

[0147] Based on the correspondence between the multiple velocity intervals and the multiple wavefield intervals, the wavefield interval corresponding to the velocity value at each point in the lateral direction is determined. Wavefield weighting is then performed based on the following formula to determine the wavefield value at each point in the lateral direction, resulting in the interpolated downlink wavefield:

[0148] , formula 2;

[0149] In the formula, and Let be the endpoint value of the i-th wavefield interval. , These are the weighting coefficients.

[0150] The inventors discovered that in wavefield extrapolation based on the phase-shifting method, for a certain depth... Wavenumber calculation can only handle velocity fields that vary vertically with depth, and it is not sufficient for velocity fields that vary laterally. In order to improve the accuracy of wavefield extrapolation, this embodiment of the invention uses a wavefield weighting method to satisfy the situation of wavefield extrapolation with lateral velocity variation by executing the specific implementation process of the above step S103, which is more consistent with the actual geological structure.

[0151] Specifically, when the computation depth is When dealing with wave fields, one can determine the velocity field. In , The range of lateral velocity changes obtained from coordinates [ The velocity range is divided into equal intervals, and defined as follows: , , , ,like Figure 5 As shown, this can be divided into equal parts. Speed ​​range; Assumption The minimum speed, i.e. , The maximum speed, i.e. Then use each velocity value respectively , , , By performing a wavefield extension, the wavefield corresponding to the transverse velocity can be obtained. , , , Similarly, it is possible to divide the downlink wave field into The wave field interval of the portion.

[0152] Next, refer to Figure 5 As shown, the velocity value of each point on the extended mesh is extracted. Using the aforementioned n-1 speed intervals as the criterion, determine the interval in which the speed value lies, such as... Within, it is based on the interval where the speed value is located. The corresponding wave field range can be searched. Based on the above formula 2, wave field weighting is performed to obtain the interpolated wave field at each grid point.

[0153] In this embodiment of the invention, the weighting coefficients are obtained in the following manner. and The weighting coefficients are determined based on the coordinate distance between each point in the horizontal direction of the extended grid and the endpoint of the wave field interval. and Since the coordinates of any point are known, and the coordinates of the two nearest whole grid points are also known, we can calculate the distances from any point to these two whole grid points, such as... Figure 5 As shown, different weight values ​​are given as weighting coefficients. , .

[0154] In an optional embodiment, in step S105 above, the interpolated downlink wavefield is subjected to inverse Fourier transform to obtain seismic data in the spatiotemporal domain. Specifically, it includes:

[0155] The interpolated downlink wave field Along Perform an inverse fast Fourier transform on the direction to obtain wavefield data in the frequency and wavenumber domains. ;

[0156] The wavefield data in the frequency wavenumber domain Along The direction is subjected to inverse fast Fourier transform to obtain seismic data in the frequency spatial domain. ;

[0157] The seismic data in the frequency spatial domain Performing a fast inverse Fourier transform along the frequency direction yields seismic data in the spatiotemporal domain. .

[0158] In this embodiment of the invention, after obtaining the downlink wave field after the above interpolation, when the downlink wave field continues to propagate downwards and encounters a strong reflection interface at a certain depth underground, i.e., a post-stack depth domain offset profile is generated. The strata have a high reflection coefficient, and the wave field, under the influence of a strong reflection interface, possesses the conditions to form a secondary source, generating back-transmitted energy. Since the wave field propagates in the frequency-wavenumber domain, while the strata reflection coefficient is obtained in the spatial domain, for the sake of computational convenience and efficiency, this application transforms both into the spatiotemporal domain, where they interact to form an instantaneous response, triggering a secondary source.

[0159] In an optional embodiment, step S106 above, which determines the seismic data excited by the secondary source based on the spatiotemporal seismic data and the pre-acquired post-stack depth domain migration profile, specifically includes:

[0160] Based on the spatiotemporal seismic data and the pre-acquired post-stack depth-domain migration profile, the seismic data excited by the secondary source are obtained according to the following formula. :

[0161] , formula 3;

[0162] in, This represents the post-stack depth domain offset profile. This represents seismic data in the spatiotemporal domain.

[0163] In this embodiment of the invention, the obtained spatiotemporal domain seismic data The reflection points are formed by interacting with the pre-acquired post-stack depth domain migration profile. Using the post-stack depth domain migration profile as the reflection coefficient, the seismic data excited by the secondary source are obtained based on Equation 3. .

[0164] In an optional embodiment, step S107 above, which involves extending the three-dimensional frequency-wavenumber domain wavefield in the frequency-wavenumber domain using a preset phase-shift method to obtain an upward wavefield, specifically includes:

[0165] Substituting the three-dimensional frequency-wavenumber domain wavefield into the following formula for upward wavefield extension, we obtain the upward wavefield:

[0166] , formula 4;

[0167] in, Indicates the upward wave field. Represents the three-dimensional frequency-wavenumber domain wavefield; The angular frequency of the seismic data. for Wave number in direction, for Wave number in direction, These are the depth coordinates of the underground strata. This represents the step size of the seismic wavefield extension. Among them, , Indicates the speed of seawater. This indicates the maximum frequency of the pre-stack common receiver point seismic gather; Denotes the continuation operator for upward propagation, where This represents the wave number when the wave field extends upwards to a depth of z. In the formula, the symbol "-" indicates that the wave field propagates upwards. In the offset velocity field Velocity value at depth position.

[0168] In this embodiment of the invention, the excitation of the secondary source changes the propagation direction of the wave field, that is, the seismic wave field propagates in the opposite direction, forming an upward wave field. ,like Figure 4 As shown, under the influence of the secondary source, the wave field propagates towards the sea surface along the propagation path of the dashed line in the figure. The ascending wave field is obtained through the above formula 4.

[0169] In an optional embodiment, step S108 above, which involves performing an inverse Fourier transform on the ascending wave field to obtain multiple wave model data in the spatiotemporal domain, specifically includes:

[0170] Determine the total multiple wave field when the upward wave field ascends to a depth of 0;

[0171] Performing an inverse Fourier transform on the total multiple wave field yields the result at a depth of [missing information]. Multiple wave model data generated at the location .

[0172] In this embodiment of the invention, the ascending wave field will originate from a depth of [depth missing]. The location was transmitted back to the sea surface and recorded again by the detector, equivalent to the aforementioned pre-stack common receiver point seismic gather. Offset profile in post-stack depth domain and offset velocity field Under the influence of the sea surface and underground depths of Multiple reflections occur between the strongly reflective interfaces, causing the primary wave to become a first-order multiple wave, which in turn becomes a second-order multiple wave, and so on, forming the total multiple wave field generated at the sea surface location, i.e., at a subsurface depth Z=0. ; and then the multiple wave fields recorded again on the sea surface Performing a fast inverse Fourier transform maps the frequency-wavenumber domain to the spatiotemporal domain, yielding the multiple wave model data generated at a subsurface depth of Z=0, i.e., at the sea surface. .

[0173] In this embodiment of the invention, through the above steps S101 to S108, the multiple wave model data formed by secondary reflection points at all underground depths can be calculated. To achieve the goal of having a common receiver point seismic gather for the entire pre-stack process. Multiple wave predictions. The acquired pre-stack common receiver gathers. With multiple wave model data By performing a matching subtraction, the seismic data after suppressing multiple waves is obtained, thus achieving multiple wave suppression.

[0174] By acquiring the pre-stack common receiver point seismic gathers from all the seismic data to be suppressed, and using the aforementioned multiple prediction method, the multiple model data corresponding to the pre-stack common receiver point seismic gathers are obtained. Then, matching and subtraction are performed, which can achieve multiple suppression for all gathers in the entire seismic data.

[0175] In this embodiment of the invention, the pre-stack common receiver point gather With multiple wave model data The specific implementation method for matching and subtraction can be found in the detailed descriptions in the relevant existing technologies, and will not be repeated here.

[0176] To provide a clearer explanation of the embodiments of the present invention, in a specific example, see [reference needed]. Figures 6-8 As shown below, the following are examples of... Figure 6 Taking the example of the specific process of multiple wave prediction and suppression in the pre-stack common receiver gather of OBN seismic data, the multiple wave prediction and suppression process of this embodiment of the invention is described in detail below:

[0177] Through the above steps S101 to S108, Figure 6 Multiple wave prediction was performed on the 3D pre-stack common receiver gathers in the OBN seismic data shown. Figure 7 The multi-wave model data shown is from Figure 6 and Figure 7 The comparison reveals that, Figure 7 Mid-wave multiple model data and Figure 6 The original multiples have a good correspondence, indicating that the multiple prediction method provided by the embodiment of the present invention has high multiple prediction accuracy and can predict the multiples of each order in the three-dimensional pre-stack common receiver point set at one time, with good actual prediction effect.

[0178] Furthermore, by Figure 6 The three-dimensional pre-stack common receiver point gather shown is Figure 7 The multi-wave model data shown is obtained by matching and subtracting. Figure 8 The seismic gathers after suppression of multiple waves shown are compared with those shown. Figure 6 and Figure 8 You can see Figure 7 The predicted multiple waves, in Figure 8 The seismic traces shown above have all been suppressed, with good suppression of multiples, and are able to suppress surface multiples in OBN seismic data.

[0179] To provide a clearer explanation of the embodiments of the present invention, in another specific example, see [link to example]. Figures 9-12 As shown below, the following are examples of... Figure 9 Taking the example of the specific process of multiple wave prediction and suppression in the pre-stack common receiver gather of OBN seismic data, the multiple wave prediction and suppression process of this embodiment of the invention is described in detail below:

[0180] Through the above steps S101 to S108, Figure 9 Multiple wave prediction was performed on the two-dimensional pre-stack common receiver gathers in the OBN seismic data shown. Figure 10 The multi-wave model data shown is from Figure 9 and Figure 10 The comparison reveals that, Figure 10 Mid-wave multiple model data and Figure 9 The original multiples have a good correspondence, indicating that the multiple prediction method provided by the embodiment of the present invention has high multiple prediction accuracy and can predict the multiples of each order in the two-dimensional pre-stack common receiver point set at one time, with good actual prediction effect.

[0181] Furthermore, by Figure 9 The three-dimensional pre-stack common receiver point gather shown is Figure 10 The multi-wave model data shown is obtained by matching and subtracting. Figure 11 The seismic gathers after suppression of multiple waves shown are compared with those shown. Figure 9 and Figure 11 You can see Figure 10 The predicted multiple waves, in Figure 11 The seismic traces shown above have all been suppressed, with good suppression of multiples, effectively suppressing surface multiples in OBN seismic data. Furthermore, Figure 11 With through Figure 9 The earthquake gathers shown are obtained by forward modeling. Figure 12 By comparing the primary wave seismic data shown, we can see... Figure 11 The shape of the central axis and Figure 12 The similarity is very high, and, Figure 11 The absence of residual multiples indicates that the multiple prediction method provided by the embodiments of the present invention can effectively predict multiples of all orders, thereby effectively suppressing multiples.

[0182] In summary, the multiple wave prediction method provided by the embodiments of the present invention, compared with the SRME method used in conventional technology, is applicable to the prediction of multiple waves in marine seismic data acquired when the excitation point and the receiver point are not in the same plane and the seabed nodes are sparsely distributed. Therefore, it can realize the prediction of multiple waves in both OBN seismic data and marine towed seismic data, thereby enabling the suppression of surface multiple waves in marine seismic data. It has a wide range of applications, strong versatility, and practical application value.

[0183] The multiple wave prediction method provided in this invention employs a phase-shift wavefield extension approach to propagate seismic data downwards. Based on a pre-acquired post-stack depth-domain migration profile, it determines the secondary sources formed at strong reflection interfaces during downward propagation of the seismic data. It then acquires the multiple waves of various orders formed by the reflection of the secondary source wavefield between the sea surface and the strong reflection interface during upward propagation, enabling the prediction of all multiple wave orders in a single operation. Furthermore, to overcome the limitation that phase-shift wavefield extension can only propagate in horizontally layered media, a wavefield interpolation weighting method is introduced. Based on the acquired migration velocity field, wavefield weighting is applied to each point in the transverse direction of the extension grid corresponding to the downflow wavefield at any depth, resulting in the interpolated downflow wavefield. This effectively solves the problem of the invariance of the transverse velocity field and is applicable to multiple wave prediction of marine-acquired seismic data in strata with complex underground structures, thus improving the accuracy of multiple wave prediction.

[0184] The multiple wave prediction method provided in this invention performs multiple wave prediction on the acquired pre-stack common receiver point seismic gathers, realizing single-gather multiple wave prediction processing. It achieves multiple wave prediction in the frequency wavenumber domain. Since frequency domain seismic data is decoupled and each frequency is independent of each other, multi-core parallel computing can be realized, resulting in fast computing speed and improved computing efficiency. Moreover, it does not require a large amount of memory space and is suitable for multiple wave prediction processing of large amounts of marine-acquired seismic data.

[0185] The multiple wave prediction method provided in this embodiment of the invention, compared with conventional multiple wave prediction methods, does not require picking up the layer information of multiple wave generation during the multiple wave prediction process. It can predict the multiple wave information of all layers contained in a given profile, which is convenient for multiple wave prediction processing of large amounts of marine seismic data.

[0186] Example 2

[0187] Based on the same inventive concept, this application also provides a method for suppressing multiple waves, referring to... Figure 13 As shown, the method includes:

[0188] S101: Perform a positive Fourier transform on the acquired pre-stack common receiver point seismic gathers to obtain the transformed seismic gathers in the three-dimensional frequency wavenumber domain.

[0189] S102: The transformed seismic gather is subjected to wavefield extension according to a preset phase shift method to obtain the downlink wavefield;

[0190] S103: Based on the acquired offset velocity field, perform wave field weighting on each point in the transverse grid corresponding to the downlink wave field at any depth to obtain the interpolated downlink wave field.

[0191] S104: Perform an inverse Fourier transform on the interpolated downlink wave field to obtain seismic data in the spatiotemporal domain;

[0192] S105: Based on the spatiotemporal seismic data and the pre-acquired post-stack depth domain migration profile, determine the seismic data triggered by the secondary source.

[0193] S106: Perform Fourier transform on the seismic data excited by the secondary source to obtain the wavefield in the three-dimensional frequency-wavenumber domain.

[0194] S107: The three-dimensional frequency wavenumber domain wavefield is extended in the frequency wavenumber domain according to a preset phase shift method to obtain the up-going wavefield.

[0195] S108: Perform an inverse Fourier transform on the upward wave field to obtain multiple wave model data in the spatiotemporal domain;

[0196] S109: Match and subtract the acquired pre-stack common receiver gathers with the spatiotemporal multiple model data to obtain the seismic data after suppressing multiples.

[0197] In this embodiment of the invention, the specific implementation process of steps S101 to S108 can be referred to the detailed description of the multiple wave prediction method in Embodiment 1 above, and will not be repeated here. The specific implementation method of matching subtraction in step S109 can be referred to the detailed description in the relevant prior art, and will not be repeated here.

[0198] The multiple suppression method provided in this invention, compared with the SRME method used in conventional technology, is applicable to the prediction of multiple waves in marine seismic data acquired when the excitation point and the receiver point are not in the same plane and the seabed nodes are sparsely distributed. Therefore, it can realize the prediction of multiple waves in both OBN seismic data and marine towed seismic data, and thus can suppress surface multiple waves in marine seismic data. It has a wide range of applications, strong versatility, and practical application value.

[0199] The multiple wave suppression method provided in this invention employs a phase-shift wavefield extension approach to propagate seismic data downwards. Based on a pre-acquired post-stack depth-domain migration profile, it determines the secondary sources formed at strong reflection interfaces during downward propagation of the seismic data. It then acquires the multiple waves of various orders formed by reflections between the sea surface and strong reflection interfaces as the secondary source wavefield propagates upwards, enabling the prediction of all multiple wave orders in a single operation. Furthermore, to overcome the limitation that phase-shift wavefield extension can only propagate in horizontally layered media, a wavefield interpolation weighting method is introduced. Based on the acquired migration velocity field, wavefield weighting is applied to each point in the transverse direction of the extension grid corresponding to the downflow wavefield at any depth, resulting in the interpolated downflow wavefield. This effectively solves the problem of the invariance of the transverse velocity field and is applicable to multiple wave prediction of marine seismic data in strata with complex underground structures, thus improving the accuracy of multiple wave suppression.

[0200] The multiple wave suppression method provided in this invention performs multiple wave prediction on the acquired pre-stack common receiver point seismic gathers, realizing single-gather multiple wave prediction processing. It achieves multiple wave prediction in the frequency wavenumber domain. Since frequency domain seismic data is decoupled and each frequency is independent of each other, multi-core parallel computing can be realized, resulting in fast computing speed and improved computing efficiency. Moreover, it does not require a large amount of memory space and is suitable for multiple wave prediction processing of large amounts of marine-acquired seismic data.

[0201] The multiple wave suppression method provided in this embodiment of the invention, compared with conventional multiple wave prediction methods, does not require picking up the layer information of multiple wave generation during the multiple wave prediction process. It can predict the multiple wave information of all layers contained in a given profile, which is convenient for multiple wave suppression processing of large amounts of marine seismic data.

[0202] Example 3

[0203] Based on the same inventive concept, embodiments of this application also provide a multiple wave prediction device, referring to... Figure 14 As shown, it includes:

[0204] The first transformation module 101 is used to perform a forward Fourier transform on the acquired pre-stack common receiver point seismic gathers to obtain the transformed seismic gathers in the three-dimensional frequency wavenumber domain.

[0205] The first extension module 102 is used to extend the wavefield of the transformed seismic gather according to a preset phase shift method to obtain the downlink wavefield.

[0206] The interpolation module 103 is used to perform wave field weighting on each point in the horizontal direction of the extended grid corresponding to the downlink wave field at any depth based on the acquired offset velocity field, so as to obtain the interpolated downlink wave field.

[0207] The second transformation module 104 is used to perform an inverse Fourier transform on the interpolated downlink wave field to obtain seismic data in the spatiotemporal domain.

[0208] The determination module 105 is used to determine the seismic data excited by the secondary source based on the seismic data in the spatiotemporal domain and the pre-acquired post-stack depth domain migration profile.

[0209] The third transformation module 106 is used to perform Fourier transform on the seismic data excited by the secondary source to obtain the wavefield in the three-dimensional frequency wavenumber domain.

[0210] The second extension module 107 is used to extend the three-dimensional frequency wavenumber domain wavefield in the frequency wavenumber domain according to a preset phase shift method to obtain an upward wavefield.

[0211] The fourth transformation module 108 is used to perform an inverse Fourier transform on the ascending wave field to obtain multiple wave model data in the spatiotemporal domain.

[0212] Example 4

[0213] Based on the same inventive concept, this application also provides a multiple wave suppression device, referring to... Figure 15 As shown, it includes:

[0214] The first transformation module 101 is used to perform a forward Fourier transform on the acquired pre-stack common receiver point seismic gathers to obtain the transformed seismic gathers in the three-dimensional frequency wavenumber domain.

[0215] The first extension module 102 is used to extend the wavefield of the transformed seismic gather according to a preset phase shift method to obtain the downlink wavefield.

[0216] The interpolation module 103 is used to perform wave field weighting on each point in the horizontal direction of the extended grid corresponding to the downlink wave field at any depth based on the acquired offset velocity field, so as to obtain the interpolated downlink wave field.

[0217] The second transformation module 104 is used to perform an inverse Fourier transform on the interpolated downlink wave field to obtain seismic data in the spatiotemporal domain.

[0218] The determination module 105 is used to determine the seismic data excited by the secondary source based on the seismic data in the spatiotemporal domain and the pre-acquired post-stack depth domain migration profile.

[0219] The third transformation module 106 is used to perform Fourier transform on the seismic data excited by the secondary source to obtain the wavefield in the three-dimensional frequency wavenumber domain.

[0220] The second extension module 107 is used to extend the three-dimensional frequency wavenumber domain wavefield in the frequency wavenumber domain according to a preset phase shift method to obtain an upward wavefield.

[0221] The fourth transformation module 108 is used to perform an inverse Fourier transform on the upward wave field to obtain multiple wave model data in the spatiotemporal domain.

[0222] The suppression module 109 is used to match and subtract the pre-stack common receiver gather from the multiple wave model data to obtain the suppressed multiple wave seismic data.

[0223] Example 5

[0224] Based on the same inventive concept, embodiments of this application also provide a computer-readable storage medium storing instructions that, when executed on a terminal, cause the terminal to perform the multi-wave prediction method or the multi-wave suppression method described above.

[0225] The computer-readable storage medium may be, for example, but not limited to, an electrical, magnetic, optical, electromagnetic, infrared, or semiconductor system, apparatus, or device, or any combination thereof. More specific examples of computer-readable storage media (a non-exhaustive list) include: an electrical connection having one or more wires; a portable computer disk drive; a hard disk drive; a random access memory (RAM); a read-only memory (ROM); an erasable programmable read-only memory (EPROM); a register; a hard disk drive; an optical fiber; a portable compact disc read-only memory (CD-ROM); an optical storage device; a magnetic storage device; or any suitable combination thereof; or any other form of computer-readable storage medium known in the art. An exemplary storage medium is coupled to a processor, enabling the processor to read information from and write information to the storage medium. Of course, the storage medium may also be a component of the processor. The processor and the storage medium may reside in an application-specific integrated circuit (ASIC). In the embodiments of this application, the computer-readable storage medium can be any tangible medium that contains or stores a program that can be used by or in conjunction with an instruction execution system, apparatus, or device.

[0226] Example 6

[0227] Based on the same inventive concept, embodiments of this application also provide a computer program product containing instructions, which, when run on a computer, causes the computer to execute the multi-wave prediction method or the multi-wave suppression method described above.

[0228] Example 7

[0229] Based on the same inventive concept, this application also provides a computer device, characterized in that it includes a processor, a communication interface, a memory, and a communication bus, wherein the processor, the communication interface, and the memory communicate with each other through the communication bus;

[0230] Memory, used to store computer programs;

[0231] When the processor executes the program stored in the memory, it implements the aforementioned multiple wave prediction method or the aforementioned multiple wave suppression method.

[0232] Figure 16A possible structural diagram of the computer device involved in the above embodiments is shown. The computer device includes a processor 1002 and a communication interface 1003. The processor 1002 is used to control and manage the operation of the computer device, for example, executing the multi-wave prediction method or the multi-wave suppression method described above, and / or other processes of the techniques described herein. The communication interface 1003 is used to support communication between the computer device and other network entities, for example, executing the steps performed by the communication unit 902 described above. The computer device may also include a memory 1001 and a bus 1004, the memory 1001 being used to store the program code and data of the computer device.

[0233] The memory 1001 may be a memory in a computer device, and the memory may include volatile memory, such as random access memory; the memory may also include non-volatile memory, such as read-only memory, flash memory, hard disk or solid-state drive; the memory may also include a combination of the above types of memory.

[0234] The processor 1002 described above can implement or execute various exemplary logic blocks, modules, and circuits described in conjunction with the disclosure of this application. The processor can be a central processing unit, a general-purpose processor, a digital signal processor, an application-specific integrated circuit (ASIC), a field-programmable gate array (FPGA), or other programmable logic devices, transistor logic devices, hardware components, or any combination thereof. It can implement or execute various exemplary logic blocks, modules, and circuits described in conjunction with the disclosure of this application. The processor can also be a combination that implements computing functions, such as including one or more microprocessor combinations, a combination of a DSP and a microprocessor, etc.

[0235] Bus 1004 can be an Extended Industry Standard Architecture (EISA) bus, etc. Bus 1004 can be divided into address bus, data bus, control bus, etc. For ease of representation, Figure 16 The bus is represented by a single thick line, but this does not mean that there is only one bus or one type of bus.

[0236] Example 6

[0237] Based on the same inventive concept, embodiments of this application also provide a chip, which includes a processor and a communication interface, the communication interface and the processor being coupled together, the processor being used to run computer programs or instructions to implement the multiple wave prediction method or the multiple wave suppression method described above.

[0238] Specifically, the chip provided in this application also includes a memory for storing computer programs or instructions.

[0239] Figure 16 The computer device in the document may also be a chip. The chip includes one or more processors 1002 and a communication interface 1003.

[0240] In some embodiments, the chip further includes a memory 1001, which may include read-only memory and random access memory, and provides operation instructions and data to the processor 1002. A portion of the memory 1001 may also include non-volatile random access memory (NVRAM).

[0241] In some implementations, memory 1001 stores elements such as execution modules or data structures, or subsets thereof, or extended sets thereof.

[0242] In this embodiment of the application, the corresponding operation is executed by calling the operation instructions stored in the memory 1001 (which may be stored in the operating system).

[0243] Through the above description of the embodiments, those skilled in the art will clearly understand that, for the sake of convenience and brevity, only the division of the above functional modules is used as an example. In practical applications, the above functions can be assigned to different functional modules as needed, that is, the internal structure of the device can be divided into different functional modules to complete all or part of the functions described above. The specific working process of the system, device, and unit described above can be referred to the corresponding process in the foregoing method embodiments, and will not be repeated here.

[0244] Since the computer-readable storage medium, computer program product, and computer device in the embodiments of this application can be applied to the above methods, the technical effects that can be obtained can also be referred to the above method embodiments. The embodiments of this application will not be repeated here.

[0245] In the several embodiments provided in this application, it should be understood that the disclosed systems, devices, and methods can be implemented in other ways. For example, the device embodiments described above are merely illustrative; for instance, the division of units is only a logical functional division, and in actual implementation, there may be other division methods. For example, multiple units or components may be combined or integrated into another system, or some features may be ignored or not executed. Furthermore, the mutual coupling or direct coupling or communication connection shown or discussed may be through some interfaces; the indirect coupling or communication connection between devices or units may be electrical, mechanical, or other forms.

[0246] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment according to actual needs.

[0247] In addition, the functional units in the various embodiments of this application can be integrated into one processing unit, or each unit can exist physically separately, or two or more units can be integrated into one unit.

[0248] The above are merely specific embodiments of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions within the technical scope disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A method for predicting multiple waves, characterized in that, include: The pre-stack common receiver point seismic gathers were subjected to a positive Fourier transform to obtain the transformed seismic gathers in the three-dimensional frequency wavenumber domain. The transformed seismic gathers are subjected to wavefield extension using a preset phase-shift method to obtain the downlink wavefield; Based on the acquired offset velocity field, wavefield weighting is performed on each point in the transverse grid corresponding to the downlink wavefield at any depth to obtain the interpolated downlink wavefield, including: Based on the acquired offset velocity field, the range of lateral velocity variation of the downflow wave field at any depth is determined, and the range of lateral velocity variation is divided into multiple velocity intervals. Wavefield extension is performed using the endpoints of the multiple velocity intervals to obtain wavefield values ​​corresponding to the endpoints of the multiple velocity intervals, thus obtaining multiple wavefield intervals corresponding to the multiple velocity intervals; Determine the wave field interval corresponding to the velocity interval of each point in the horizontal direction of the extended grid, and perform wave field weighting based on the two endpoints of the wave field interval to obtain the interpolated downlink wave field; The interpolated downflow wave field is subjected to inverse Fourier transform to obtain seismic data in the spatiotemporal domain. Based on the aforementioned spatiotemporal seismic data and the pre-acquired post-stack depth-domain migration profile, the seismic data excited by the secondary source are determined, including: Based on the spatiotemporal seismic data and the pre-acquired post-stack depth-domain migration profile, the seismic data excited by the secondary source are obtained according to the following formula: ； in, This represents the post-stack depth domain offset profile. Represents seismic data in the spatiotemporal domain; Fourier transform is performed on the seismic data excited by the secondary source to obtain the three-dimensional frequency wavenumber domain wavefield. The three-dimensional frequency-wavenumber domain wavefield is extended in the frequency-wavenumber domain according to a preset phase-shift method to obtain an upward wavefield, including: Substituting the three-dimensional frequency-wavenumber domain wavefield into the following formula for upward wavefield extension, we obtain the upward wavefield: in, Indicates the upward wave field. Represents the three-dimensional frequency-wavenumber domain wavefield; The angular frequency of the seismic data. for Wave number in direction, for Wave number in direction, These are the depth coordinates of the underground strata. Indicates the step size of the seismic wavefield extension; Denotes the continuation operator for upward propagation, where This represents the wave number when the wave field extends upwards to a depth of z. In the formula, the symbol "-" indicates that the wave field propagates upwards. In the offset velocity field Velocity value at depth location; The upward wave field is subjected to inverse Fourier transform to obtain multiple wave model data in the spatiotemporal domain.

2. The method as described in claim 1, characterized in that, The process involves determining the lateral velocity variation range of the downflow wavefield at any depth based on the acquired offset velocity field, and dividing the lateral velocity variation range into multiple velocity intervals, including: At any depth, the range of lateral velocity variation of the downlink wave field is obtained based on the lateral and longitudinal coordinates in the offset velocity field. The velocity within the lateral velocity variation range is divided into multiple velocity end values ​​at equal intervals; The range of speeds between any two adjacent speed values ​​is considered as a speed interval, resulting in multiple speed intervals.

3. The method as described in claim 2, characterized in that, The step of performing wavefield extrapolation using the endpoints of the multiple velocity intervals to obtain wavefield endpoints corresponding to the endpoints of the multiple velocity intervals, and obtaining multiple wavefield intervals corresponding to the multiple velocity intervals, includes: Wavefield extension is performed using the multiple velocity endpoints to obtain wavefield endpoints corresponding to each velocity endpoint. The wavefield range between each two adjacent wavefield endpoints is taken as a wavefield interval, resulting in multiple wavefield intervals corresponding to the multiple velocity intervals.

4. The method as described in claim 3, characterized in that, The process of determining the wavefield interval corresponding to the velocity interval of each point in the horizontal direction of the extended grid, and performing wavefield weighting based on the two endpoints of the wavefield interval to obtain the interpolated downlink wavefield includes: Determine the velocity range in which the velocity value of each point in the horizontal direction of the extended grid lies; Based on the correspondence between the multiple velocity intervals and the multiple wavefield intervals, the wavefield interval corresponding to the velocity value at each point in the lateral direction is determined. Wavefield weighting is then performed based on the following formula to determine the wavefield value at each point in the lateral direction, resulting in the interpolated downlink wavefield: ； In the formula, The angular frequency of the seismic data. for Wave number in direction, for Wave number in direction, and Let be the endpoint value of the i-th wavefield interval. , These are the weighting coefficients.

5. The method as described in claim 4, characterized in that, Also includes: The weighting coefficients are obtained in the following manner. and : The weighting coefficients are determined based on the coordinate distance between each point in the horizontal direction of the extended grid and the endpoint of the wave field interval. and .

6. The method as described in claim 1, characterized in that, The step of performing a positive Fourier transform on the acquired pre-stack common receiver point seismic gathers to obtain the transformed seismic gathers in the three-dimensional frequency and wavenumber domain includes: Perform a one-dimensional fast Fourier transform along the time direction on each trace in the pre-stack common receiver point seismic trace set to obtain the frequency-spatial domain seismic trace set. Perform a two-dimensional fast Fourier transform on each trace in the frequency-space domain seismic trace set along the lateral direction to obtain a two-dimensional frequency-wavenumber domain seismic trace set. Perform a three-dimensional fast Fourier transform on each trace in the two-dimensional frequency-wavenumber domain seismic gather along the longitudinal direction to obtain the transformed seismic gather in the three-dimensional frequency-wavenumber domain.

7. The method as described in claim 1, characterized in that, The process of performing wavefield extension on the transformed seismic gathers according to a preset phase shift method to obtain the downlink wavefield includes: Substituting the transformed seismic gathers into the following formula for downward wavefield extension, we obtain the downward wavefield: ； in, Indicates the downlink wave field. Represents the transformed seismic gathers; The angular frequency of the seismic data. for Wave number in direction, for Wave number in direction, These are the depth coordinates of the underground strata. This represents the step size of the seismic wavefield extension, where, , Indicates the speed of seawater. This indicates the maximum frequency of the pre-stack common receiver point seismic gather; Denotes the continuation operator for downward propagation, where This represents the wave number when the wave field extends downwards to a depth of z. In the formula, the symbol "+" indicates that the wave field propagates downwards. In the offset velocity field Velocity value at depth position.

8. The method as described in claim 7, characterized in that, The step of performing an inverse Fourier transform on the interpolated downlink wavefield to obtain seismic data in the spatiotemporal domain includes: The interpolated downlink wave field Performing an inverse fast Fourier transform along the y-direction yields wavefield data in the frequency and wavenumber domains. ; The wavefield data in the frequency wavenumber domain Performing an inverse fast Fourier transform along the x-direction yields seismic data in the frequency spatial domain. ; The seismic data in the frequency spatial domain Performing a fast inverse Fourier transform along the frequency direction yields seismic data in the spatiotemporal domain. .

9. The method as described in claim 1, characterized in that, The step of performing an inverse Fourier transform on the ascending wave field to obtain multiple wave model data in the spatiotemporal domain includes: Determine the total multiple wave field when the upward wave field ascends to a depth of 0; Performing an inverse Fourier transform on the total multiple wave field yields the result at a depth of [missing information]. Multiple wave model data generated at the location .

10. The method as described in claim 1, characterized in that, Also includes: Pre-stack common receiver point seismic gathers were obtained using the following method: If the seismic data acquired at sea is OBN acquired seismic data, the shot-receiver points of the OBN acquired seismic data are interchanged, and the pre-stack common receiver point gather is extracted from the seismic data obtained by the interchange. If the seismic data acquired at sea is acquired by a towed cable, the pre-stack common receiver gather is extracted from the towed cable acquired seismic data.

11. A method for suppressing multiple waves, characterized in that, include: The pre-stack common receiver gathers are matched and subtracted from the multiple wave model data determined by the multiple wave prediction method according to any one of claims 1-10 to obtain the seismic data after suppressing multiple waves.

12. A multiple wave prediction device, characterized in that, include: The first transformation module is used to perform a forward Fourier transform on the acquired pre-stack common receiver point seismic gathers to obtain the transformed seismic gathers in the three-dimensional frequency wavenumber domain. The first extension module is used to extend the wavefield of the transformed seismic gather according to a preset phase shift method to obtain the downlink wavefield. The interpolation module is used to perform wavefield weighting processing on each point in the transverse grid corresponding to the downlink wavefield at any depth, based on the acquired offset velocity field, to obtain the interpolated downlink wavefield, including: Based on the acquired offset velocity field, the range of lateral velocity variation of the downflow wave field at any depth is determined, and the range of lateral velocity variation is divided into multiple velocity intervals. Wavefield extension is performed using the endpoints of the multiple velocity intervals to obtain wavefield values ​​corresponding to the endpoints of the multiple velocity intervals, thus obtaining multiple wavefield intervals corresponding to the multiple velocity intervals; Determine the wave field interval corresponding to the velocity interval of each point in the horizontal direction of the extended grid, and perform wave field weighting based on the two endpoints of the wave field interval to obtain the interpolated downlink wave field; The second transformation module is used to perform an inverse Fourier transform on the interpolated downlink wave field to obtain seismic data in the spatiotemporal domain. The determination module is used to determine the seismic data excited by the secondary source based on the spatiotemporal domain seismic data and the pre-acquired post-stack depth domain migration profile, including: Based on the spatiotemporal seismic data and the pre-acquired post-stack depth-domain migration profile, the seismic data excited by the secondary source are obtained according to the following formula: ； in, This represents the post-stack depth domain offset profile. Represents seismic data in the spatiotemporal domain; The third transformation module is used to perform Fourier transform on the seismic data excited by the secondary source to obtain the wavefield in the three-dimensional frequency wavenumber domain. The second extension module is used to extend the three-dimensional frequency wavenumber domain wavefield in the frequency wavenumber domain according to a preset phase shift method to obtain an upward wavefield, including: Substituting the three-dimensional frequency-wavenumber domain wavefield into the following formula for upward wavefield extension, we obtain the upward wavefield: in, Indicates the upward wave field. Represents the three-dimensional frequency-wavenumber domain wavefield; The angular frequency of the seismic data. for Wave number in direction, for Wave number in direction, These are the depth coordinates of the underground strata. Indicates the step size of the seismic wavefield extension; Denotes the continuation operator for upward propagation, where This represents the wave number when the wave field extends upwards to a depth of z. In the formula, the symbol "-" indicates that the wave field propagates upwards. In the offset velocity field Velocity value at depth location; The fourth transformation module is used to perform an inverse Fourier transform on the ascending wave field to obtain multiple wave model data in the spatiotemporal domain.

13. A multiple wave suppression device, characterized in that, include: The first transformation module is used to perform a forward Fourier transform on the acquired pre-stack common receiver point seismic gathers to obtain the transformed seismic gathers in the three-dimensional frequency wavenumber domain. The first extension module is used to extend the wavefield of the transformed seismic gather according to a preset phase shift method to obtain the downlink wavefield. The interpolation module is used to perform wavefield weighting processing on each point in the transverse grid corresponding to the downlink wavefield at any depth, based on the acquired offset velocity field, to obtain the interpolated downlink wavefield, including: Based on the acquired offset velocity field, the range of lateral velocity variation of the downflow wave field at any depth is determined, and the range of lateral velocity variation is divided into multiple velocity intervals. Wavefield extension is performed using the endpoints of the multiple velocity intervals to obtain wavefield values ​​corresponding to the endpoints of the multiple velocity intervals, thus obtaining multiple wavefield intervals corresponding to the multiple velocity intervals; Determine the wave field interval corresponding to the velocity interval of each point in the horizontal direction of the extended grid, and perform wave field weighting based on the two endpoints of the wave field interval to obtain the interpolated downlink wave field; The second transformation module is used to perform an inverse Fourier transform on the interpolated downlink wave field to obtain seismic data in the spatiotemporal domain. The determination module is used to determine the seismic data excited by the secondary source based on the spatiotemporal domain seismic data and the pre-acquired post-stack depth domain migration profile, including: Based on the spatiotemporal seismic data and the pre-acquired post-stack depth-domain migration profile, the seismic data excited by the secondary source are obtained according to the following formula: ； in, This represents the post-stack depth domain offset profile. Represents seismic data in the spatiotemporal domain; The third transformation module is used to perform Fourier transform on the seismic data excited by the secondary source to obtain the wavefield in the three-dimensional frequency wavenumber domain. The second extension module is used to extend the three-dimensional frequency wavenumber domain wavefield in the frequency wavenumber domain according to a preset phase shift method to obtain an upward wavefield, including: Substituting the three-dimensional frequency-wavenumber domain wavefield into the following formula for upward wavefield extension, we obtain the upward wavefield: in, Indicates the upward wave field. Represents the three-dimensional frequency-wavenumber domain wavefield; The angular frequency of the seismic data. for Wave number in direction, for Wave number in direction, These are the depth coordinates of the underground strata. Indicates the step size of the seismic wavefield extension; Denotes the continuation operator for upward propagation, where This represents the wave number when the wave field extends upwards to a depth of z. In the formula, the symbol "-" indicates that the wave field propagates upwards. In the offset velocity field Velocity value at depth location; The fourth transformation module is used to perform an inverse Fourier transform on the upward wave field to obtain multiple wave model data in the spatiotemporal domain. The suppression module is used to match and subtract the pre-stack common receiver gather from the multiple wave model data to obtain the suppressed multiple wave seismic data.

14. A computer-readable storage medium storing instructions that, when executed on a terminal, cause the terminal to perform the multiple wave prediction method as described in any one of claims 1-10 or the multiple wave suppression method as described in claim 11.

15. A computer device, characterized in that, It includes a processor, a communication interface, a memory, and a communication bus, wherein the processor, the communication interface, and the memory communicate with each other through the communication bus; Memory, used to store computer programs; The processor, when executing a program stored in memory, implements the multiple wave prediction method as described in any one of claims 1-10 or the multiple wave suppression method as described in claim 11.

Images