A multiple wave suppression method based on characteristic reflection nonlinear convolution model

Abstract

The invention discloses a multiple wave suppression method based on a characteristic reflection nonlinear convolution model, including collecting multi-component seafloor node data and preprocessing the data; constructing a plane wave basis function according to the preprocessed parameters, converting the frequency space domain data into time plane wave domain data; introducing an adaptive matching factor to achieve adaptive matching of the time plane wave domain data, separating upgoing and downgoing wave seafloor node data; respectively converting the separated upgoing and downgoing wave seafloor node data into the frequency domain, performing upgoing and downgoing wave deconvolution; converting the deconvolution result back into the time-space domain, obtaining a single reflection Green's function, that is, single reflection seafloor node data, to achieve suppression. The method can eliminate the influence of source wavelets, and solve the problem that the traditional free surface related multiple wave suppression processing technology cannot adapt to the seismic data collected by the seafloor node observation method.

Description

Technical Field

[0001] The invention relates to the technical field of geophysical exploration, and in particular to a multiple wave suppression method based on a characteristic reflection nonlinear convolution model. Background Art

[0002] Seabed node acquisition technology is conducive to high-precision seismic imaging due to its long offset, all-round observation, and rich low-frequency signals. It plays an important role in complex structural imaging, but the unique acquisition method of the seabed node observation system also brings challenges to subsequent data processing. In the seabed node observation system, the seismic sources are densely excited on the sea surface, and the observation nodes are located on the seabed, allowing sparse and irregular layout. The obvious depth difference of the shot check point position makes the seabed node data have the characteristics of dual reference planes. Due to the existence of two strong characteristic reflection surfaces, the seawater-air interface and the seawater-seabed interface, the seabed node data contains not only upgoing wave signals from the strata, but also downgoing wave signals from the sea surface. Multiple waves related to the water body and free surface are also widely developed.

[0003] Since the current seismic data processing process and imaging methods are mainly aimed at primary reflections, it is usually necessary to separate the uplink and downlink signals before performing velocity analysis and imaging. The presence of multiple waves will also affect the quality of velocity analysis and imaging. Multiple wave suppression is an important issue in seafloor node data processing. The traditional free surface related multiple wave suppression method requires that the collected data is regular and complete, but in the seafloor node observation system, the seafloor observation nodes are usually sparsely distributed, and the collected seismic data are insufficiently sampled in the common shot gathers. This makes the traditional free surface related multiple wave suppression method unable to adapt to the seismic data collected by the seafloor node observation method. Summary of the invention

[0004] The purpose of this section is to summarize some aspects of embodiments of the present invention and briefly introduce some preferred embodiments. Some simplifications or omissions may be made in this section and the specification abstract and the invention title of this application to avoid blurring the purpose of this section, the specification abstract and the invention title, and such simplifications or omissions cannot be used to limit the scope of the present invention.

[0005] In view of the above existing problems, the present invention is proposed.

[0006] Therefore, the present invention provides a multiple wave suppression method based on a characteristic reflection nonlinear convolution model, which can solve the problem that the traditional free surface related multiple wave suppression method cannot adapt to the seismic data collected by the seabed node observation method.

[0007] In order to solve the above technical problems, the present invention provides the following technical solution, a multiple wave suppression method based on a characteristic reflection nonlinear convolution model, comprising:

[0008] Collect multi-component seafloor node data and pre-process the data;

[0009] Constructing plane wave basis functions according to the preprocessed parameters to convert frequency space domain data into time plane wave domain data;

[0010] Adaptive matching factors are introduced to achieve adaptive matching of time plane wave domain data and separate upgoing and downgoing seafloor node data;

[0011] The separated upgoing and downgoing wave seafloor node data are converted into the frequency domain respectively, and the upgoing and downgoing wave deconvolution is performed;

[0012] The deconvolution result is converted back to the time-space domain to obtain the first-reflection Green's function, that is, the first-reflection seabed node data, to achieve suppression.

[0013] As a preferred solution of the multiple wave suppression method based on the characteristic reflection nonlinear convolution model described in the present invention, the preprocessing includes inputting multi-component seabed node data into the header keyword program, obtaining observation system information, and sorting out the common detection point gather data of the pressure component and the vertical velocity component.

[0014] As a preferred solution of the multiple wave suppression method based on the characteristic reflection nonlinear convolution model described in the present invention, the observation system information includes the number of time sampling points, the time sampling interval, the spatial position of the observation node relative to the observation system, the excitation position of the sea surface source relative to the observation system and the inclination sampling interval of the seismic data.

[0015] As a preferred solution of the multiple wave suppression method based on the characteristic reflection nonlinear convolution model described in the present invention, wherein: the plane wave basis function includes:

[0016] In the plane wave domain, according to the physical relationship between the magnitudes of the pressure component and the vertical velocity component data, the amplitudes of the two are adjusted to the same magnitude. The magnitude adjustment relationship is as follows:

[0017]

[0018] Where P′ represents the pressure component after magnitude adjustment, k Z represents the vertical wave number component, ρ represents the seawater density, ω represents the signal wave number after Fourier transformation, P represents the pressure component before magnitude adjustment, and p X represents the lateral ray parameter component of the local plane wave signal, v sea_water Represents the speed of sound wave propagation in sea water.

[0019] As a preferred solution of the multiple wave suppression method based on the characteristic reflection nonlinear convolution model described in the present invention, the adaptive matching factor includes selecting the direct wave and obtaining it by the least square method, and the equation is constructed as follows:

[0020]

[0021] Among them, P represents the pressure component data, V Z represents the vertical velocity component data, w represents the best matched filter of the pressure component and the vertical velocity component, Represents the convolution matrix composed of vertical velocity component data, Representation Matrix The transpose of P Represents the matrix composed of pressure component data, and d P The matrix arrangement of refers to the definition of the general convolution matrix equation.

[0022] As a preferred solution of the multiple wave suppression method based on the characteristic reflection nonlinear convolution model described in the present invention, wherein: the direct wave includes, and the trajectory of the direct wave in the time and space domain is expressed as:

[0023]

[0024]

[0025] Where t represents the arrival time curve of the direct wave signal in the time and space domain, that is, the propagation time of the source signal in the seawater before reaching the detector, v sea_water represents the sound wave propagation speed of seawater, z represents the seawater depth where the geophone is placed, x0 represents the transverse spatial coordinate where the geophone is placed, x represents the transverse spatial coordinate where the source is excited, y represents the longitudinal spatial coordinate where the geophone is placed, y0 represents the longitudinal spatial coordinate where the source is excited, Indicates the lateral spatial distance between the earthquake source and the geophone. Indicates the longitudinal spatial distance between the earthquake source and the detector.

[0026] As a preferred solution of the multiple wave suppression method based on the characteristic reflection nonlinear convolution model described in the present invention, wherein: the direct wave also includes, and the direct wave trajectory in the time plane domain is expressed as follows:

[0027]

[0028] Where t0 represents the arrival time curve of the direct wave in the plane wave domain, p x represents the transverse ray parameter component of the direct wave, p y It represents the longitudinal ray parameter component of the direct wave. The definitions of the transverse and longitudinal ray parameters are as follows:

[0029]

[0030]

[0031] The physical meaning of each variable is the same as that defined in the above formula.

[0032] As a preferred solution of the multiple wave suppression method based on the characteristic reflection nonlinear convolution model described in the present invention, wherein: the upgoing wave and downgoing wave seabed node data include:

[0033] The formation process of the up and down waves of the seabed node data, that is, the free surface related multiple waves contained therein, is expressed in the plane wave domain as follows:

[0034]

[0035]

[0036] Where p represents the ray parameter vector of the signal (including the component information of the ray parameters in each direction), ω represents the wave number of the signal after Fourier transformation, and z r represents the seawater depth where the geophone is placed, z0 represents the water depth of the source excitation, exp() represents the natural exponential function, U(p,ω,z r ) represents the upgoing wave signal in the frequency plane wave domain, P(p,ω,z0) represents the full wave field signal observed at the sea surface, i represents the imaginary unit of the complex signal, k z represents the vertical wave number component of the wave field, g(p,ω) represents the primary reflection Green's function (i.e., the primary wave information to be extracted), a(ω) represents the wave impedance of the sea surface (usually taken as a constant -1), s(p,ω) represents the excitation signal of the earthquake source, and D(p,ω,z r ) represents the downlink wave signal in the frequency plane wave domain, P d (p,ω,z r ) represents the direct wave signal.

[0037] As a preferred solution of the multiple wave suppression method based on the characteristic reflection nonlinear convolution model described in the present invention, wherein: the reflection Green's function includes:

[0038]

[0039] Among them, P (U) (p,ω,z r ) represents the upgoing wave signal in the frequency plane wave domain (complex data points on a certain ray parameter), P (D) (p,ω,z r) represents the downlink wave signal in the frequency plane wave domain (complex data points on a certain ray parameter), g(p,ω,z r ) represents the final obtained first reflection Green's function (i.e. the first wave information that needs to be extracted).

[0040] As a preferred solution of the multiple wave suppression method based on the characteristic reflection nonlinear convolution model described in the present invention, wherein: the first deconvolution includes:

[0041] In order to ensure the stability of the deconvolution process, a stability factor is introduced and the deconvolution process is improved as follows:

[0042]

[0043] in, Indicates P (D) (p,ω,z r ), and ε represents a small positive number for stabilizing the solution.

[0044] Beneficial effects of the present invention: The method of the present invention is based on a feedback loop model, and obtains a nonlinear convolution model of up and down waves of seabed node data. The model shows that the underground primary reflection Green's function can be extracted by deconvolution of the up and down waves. The method assumes that the multiple wave prediction relationships of plane wave components in different directions in the seabed node data are independent of each other, so deconvolution can be performed frequency by frequency and plane wave direction by plane wave in the plane wave domain, with high computational efficiency, and only local plane wave decomposition of the common detection point gathers of the seabed node data is required, thus avoiding the problem of sparse spatial sampling of the common shot gathers of the seabed node data. BRIEF DESCRIPTION OF THE DRAWINGS

[0045] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the following briefly introduces the drawings required for describing the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative labor. Among them:

[0046] Figure 1 A schematic flow chart of a multiple wave suppression method based on a characteristic reflection nonlinear convolution model provided by an embodiment of the present invention;

[0047] Figure 2 A schematic diagram of a feedback loop model of a multiple wave suppression method based on a characteristic reflection nonlinear convolution model provided by an embodiment of the present invention.

[0048] Figure 3 A schematic diagram of a seafloor node observation system for a multiple wave suppression method based on a characteristic reflection nonlinear convolution model provided by an embodiment of the present invention;

[0049] Figure 4 A schematic diagram of application data of a multiple wave suppression method based on a characteristic reflection nonlinear convolution model provided by one embodiment of the present invention;

[0050] Figure 5 A schematic diagram of uplink and downlink wave separation of a multiple wave suppression method based on a characteristic reflection nonlinear convolution model provided by an embodiment of the present invention;

[0051] Figure 6 A schematic diagram of extracting a primary reflection Green's function by deconvolution of up and down waves of a multiple wave suppression method based on a characteristic reflection nonlinear convolution model provided by an embodiment of the present invention;

[0052] Figure 7 A comparative analysis diagram of a multiple wave suppression method based on a characteristic reflection nonlinear convolution model provided by an embodiment of the present invention and the common detection point gathers of the original seabed node data; DETAILED DESCRIPTION

[0053] In order to make the above-mentioned purposes, features and advantages of the present invention more obvious and easy to understand, the specific implementation methods of the present invention are described in detail below in conjunction with the drawings of the specification. Obviously, the described embodiments are part of the embodiments of the present invention, but not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary persons in the art without creative work should fall within the scope of protection of the present invention.

[0054] In the following description, many specific details are set forth to facilitate a full understanding of the present invention, but the present invention may also be implemented in other ways different from those described herein, and those skilled in the art may make similar generalizations without violating the connotation of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0055] Secondly, the term "one embodiment" or "embodiment" as used herein refers to a specific feature, structure, or characteristic that may be included in at least one implementation of the present invention. The term "in one embodiment" that appears in different places in this specification does not necessarily refer to the same embodiment, nor does it refer to a separate or selective embodiment that is mutually exclusive with other embodiments.

[0056] The present invention is described in detail with reference to schematic diagrams. When describing the embodiments of the present invention, for the sake of convenience, the cross-sectional diagrams showing the device structure will not be partially enlarged according to the general scale, and the schematic diagrams are only examples, which should not limit the scope of protection of the present invention. In addition, in actual production, the three-dimensional dimensions of length, width and depth should be included.

[0057] At the same time, in the description of the present invention, it should be noted that the directions or positional relationships indicated by the terms "upper, lower, inner and outer" are based on the directions or positional relationships shown in the drawings, and are only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific direction, be constructed and operated in a specific direction, and therefore cannot be understood as limiting the present invention. In addition, the terms "first, second or third" are only used for descriptive purposes and cannot be understood as indicating or implying relative importance.

[0058] In the present invention, unless otherwise clearly specified and limited, the terms "install, connect, connect" should be understood in a broad sense, for example: it can be a fixed connection, a detachable connection or an integral connection; it can also be a mechanical connection, an electrical connection or a direct connection, or it can be indirectly connected through an intermediate medium, or it can be the internal communication of two components. For ordinary technicians in this field, the specific meanings of the above terms in the present invention can be understood according to specific circumstances.

[0059] Example 1

[0060] Reference Figure 1-7 , which is the first embodiment of the present invention, provides a multiple wave suppression method based on a characteristic reflection nonlinear convolution model, comprising:

[0061] S1: Collect multi-component seafloor node data and pre-process the data;

[0062] Furthermore, the preprocessing includes inputting multi-component seafloor node data in standard SEG-Y format into the header keyword program, obtaining observation system information, and sorting out the common detection point gather data of the pressure component and the vertical velocity component according to the header keyword information in the SEG-Y format.

[0063] It should be noted that the observation system parameters are: the number of time sampling points, the time sampling interval, the spatial position of the observation node relative to the observation system, the excitation position of the sea surface source relative to the observation system and the dip sampling interval of the seismic data.

[0064] It should also be noted that Figure 2 It can be seen that the earthquake sources are densely excited at the sea surface, while the detectors (represented by triangles) are sparsely and irregularly arranged on the seabed.

[0065] S2: constructing plane wave basis functions according to the preprocessed parameters and converting the frequency space domain data into time plane wave domain data;

[0066] Furthermore, a one-dimensional fast Fourier transform is performed on the data channels to convert the pressure component and vertical velocity component of the co-detection point gather seismic data d(x, t) in the time (t)-space (x) domain to the frequency-space domain d(x, ω).

[0067] It should be noted that the plane wave basis function is constructed according to the observation system parameters, and the frequency-space domain data d(x,ω) is converted into plane wave domain data s(p,t). In the plane wave domain, the pressure component P(x,t) and the vertical velocity component V z The physical relationship between the (x, t) data magnitudes is to adjust the amplitude levels of the two to the same magnitude. The magnitude adjustment relationship is as follows:

[0068]

[0069] Where P′ represents the pressure component after magnitude adjustment, k Z represents the vertical wave number component, ρ represents the seawater density, ω represents the signal wave number after Fourier transformation, P represents the pressure component before magnitude adjustment, and p X represents the lateral ray parameter component of the local plane wave signal, v sea_water Represents the speed of sound wave propagation in sea water.

[0070] Furthermore, the dip range of the seabed node seismic data can be calculated based on the seawater velocity, and the plane wave basis function required for local plane wave decomposition can be constructed in combination with the dip sampling interval of the seismic data to perform local plane wave decomposition of the seismic data.

[0071] It should be noted that when performing fast Fourier transform on time-space domain data by channel, multi-threaded parallel computing is used to improve processing efficiency.

[0072] S3: Introduce adaptive matching factors to achieve adaptive matching of time plane wave domain data and separate upgoing and downgoing seafloor node data;

[0073] Furthermore, an adaptive matching factor w is introduced to achieve adaptive matching of pressure component data and vertical velocity component data, and separate the upgoing wave and downgoing wave seabed node data.

[0074] It should be noted that the upgoing wave component is U′=(P+w*V Z ), the downlink wave component is D′=(Pw*V Z ).

[0075] Furthermore, the adaptive matching factor w is obtained using the least squares method, and the equation is constructed as follows:

[0076]

[0077] Among them, P represents the pressure component data, V Z represents the vertical velocity component data, w represents the best matched filter of the pressure component and the vertical velocity component, Represents the convolution matrix composed of vertical velocity component data, Representation Matrix The transpose of P Represents the matrix composed of pressure component data, and d P The matrix arrangement of refers to the definition of the general convolution matrix equation.

[0078] It should be noted that the direct wave signal in the seafloor node data should be selected to construct the equation to solve the adaptive matching factor. This is because the direct wave in the seafloor node data belongs to the downlink wave. The adaptive matching factor obtained by solving can minimize the downlink wave energy in the data, thereby separating the uplink and downlink waves of the seafloor node data. At the same time, the trajectory of the direct wave signal in the seafloor node data is relatively fixed, which can be obtained by an analytical formula, and then the range of the direct wave signal data window is estimated according to the length of the seismic signal wavelet.

[0079] Furthermore, the trajectory of the direct wave in the seabed node data in the time and space domain is expressed as:

[0080]

[0081]

[0082] Where t represents the arrival time curve of the direct wave signal in the time and space domain (that is, the propagation time of the source signal in the seawater before it reaches the detector), v sea_water represents the sound wave propagation speed of seawater, z represents the seawater depth where the geophone is placed, x0 represents the transverse spatial coordinate where the geophone is placed, x represents the transverse spatial coordinate where the source is excited, y represents the longitudinal spatial coordinate where the geophone is placed, y0 represents the longitudinal spatial coordinate where the source is excited, Indicates the lateral spatial distance between the earthquake source and the geophone. Indicates the longitudinal spatial distance between the earthquake source and the detector.

[0083] It should be noted that, since the present invention performs data processing in the time plane wave domain, the above direct wave trajectory is converted into the time plane wave domain and expressed as:

[0084]

[0085] Where t0 represents the arrival time curve of the direct wave in the plane wave domain, p x represents the transverse ray parameter component of the direct wave, p y It represents the longitudinal ray parameter component of the direct wave. The definitions of the transverse and longitudinal ray parameters are as follows:

[0086]

[0087]

[0088] The physical meaning of each variable is the same as that defined in the above formula.

[0089] It should be noted that the direct wave trajectory in the time plane wave domain can be calculated from the direct wave trajectory t in the time domain. After separating the seafloor node data into upgoing and downgoing waves, the nonlinear convolution model of the upgoing and downgoing waves in the seafloor node data can be obtained by combining the characteristics of the feedback loop model and the seafloor node observation system. The feedback loop model is a classic model that describes the multiple waves related to the free surface of seismic data. It regards the upgoing reflection from the underground as a "secondary source" and inputs it into the underground again;

[0090] It should be noted that the above feedback loop model usually assumes that the seismic source and signal receiver in the seismic data observation system are both located on the free surface, but in the seabed node observation system, the signal receiver is deployed on the seabed, which does not meet the basic assumptions of the model.

[0091] Furthermore, the present invention derives a formation model of free surface related multiple waves under the seabed node observation system, which describes the formation process of the uplink and downlink waves of the seabed node data, i.e., the free surface related multiple waves contained therein, and is expressed in the plane wave domain as follows:

[0092]

[0093]

[0094] Where p represents the ray parameter vector of the signal (including the component information of the ray parameters in each direction), ω represents the wave number of the signal after Fourier transformation, and z r represents the seawater depth where the geophone is placed, z0 represents the water depth of the source excitation, exp() represents the natural exponential function, U(p,ω,z r ) represents the upgoing wave signal in the frequency plane wave domain, P(p,ω,z0) represents the full wave field signal observed at the sea surface, i represents the imaginary unit of the complex signal, k z represents the vertical wave number component of the wave field, g(p,ω) represents the primary reflection Green's function (i.e., the primary wave information to be extracted), a(ω) represents the wave impedance of the sea surface (usually taken as a constant -1), s(p,ω) represents the excitation signal of the earthquake source, and D(p,ω,z r ) represents the downlink wave signal in the frequency plane wave domain, P d (p,ω,z r ) represents the direct wave signal.

[0095] S4: Convert the separated upgoing and downgoing wave seafloor node data into the frequency domain respectively, and perform upgoing and downgoing wave deconvolution.

[0096] Furthermore, one-dimensional fast Fourier transform is performed on the data channel to convert the separated upgoing wave U′(p, t) and downgoing wave data D′(p, t) in the plane wave domain into the frequency domain s(p, ω), and deconvolution of the upgoing and downgoing waves is performed on a plane wave direction-by-plane wave and frequency-by-frequency basis.

[0097] Furthermore, according to the separated up and down wave data in the plane wave domain and the derived mathematical model of the up and down waves, the underground reflection Green's function can be obtained to include:

[0098]

[0099] Among them, P (U) (p,ω,z r ) represents the upgoing wave signal in the frequency plane wave domain (complex data points on a certain ray parameter), P (D) (p,ω,z r ) represents the downlink wave signal in the frequency plane wave domain (complex data points on a certain ray parameter), g(p,ω,z r ) represents the final obtained first reflection Green's function (i.e. the first wave information that needs to be extracted).

[0100] It should be noted that, according to the above formula, the separated uplink and downlink waves are deconvolved in the frequency plane wave domain to obtain the desired underground primary reflection Green's function. When processing, it is necessary to perform fast Fourier transform on the obtained uplink and downlink wave data one by one, and after converting to the frequency plane wave domain, perform the above deconvolution process one by one in plane wave direction and one by one in frequency.

[0101] Furthermore, in order to ensure the stability of the deconvolution process, a stability factor is introduced and the deconvolution process is approximated as follows:

[0102]

[0103] in, Indicates P (D) (p,ω,z r ), ε represents a small positive number for stable solution, usually 0.1, and the meanings of other variables have been mentioned before.

[0104] It should be noted that a stable deconvolution result can be obtained according to the above formula, where exp(2ik z z r ) represents a downward extension operator, which performs reference plane correction on the deconvolution result. This operator can be included in the deconvolution result.

[0105] S5: Convert the deconvolution result back to the time-space domain to obtain the first reflection Green's function, that is, the first reflection seabed node data, to achieve suppression.

[0106] The result obtained by the deconvolution is converted back to the time-space domain, which is the inverse process corresponding to the S2 local plane wave decomposition.

[0107] It should be noted that the deconvolution process eliminates the influence of the source wavelet and expands the frequency band, so when converting back to the time-space domain, a higher frequency band range needs to be calculated instead of only considering the frequency band range of the source wavelet as in S2.

[0108] Example 2

[0109] Reference Figure 3-7 , which is an embodiment of the present invention, provides a multiple wave suppression method based on a characteristic reflection nonlinear convolution model. In order to verify the beneficial effects of the present invention, scientific demonstration is carried out through experiments.

[0110] by Figure 4 The common detection point gather of the seabed node seismic data shown in the figure is taken as an example as input (the pressure and vertical velocity components are input at the same time, and the processing methods of the two are similar). The number of sampling points in the data space direction is 900, and the sampling interval is 5 meters; the number of time sampling points is 5000, and the sampling interval is 0.3 milliseconds; the range of the plane wave inclination direction contained in the data is determined according to the seawater velocity. The seawater velocity is 1500m / s, then the maximum ray parameter of the data plane wave decomposition is 1 / 1500=0.000667s / m, the sampling interval of the ray parameter of the plane wave decomposition is set to 0.000003s / m, and the number of samples is 601 directions.

[0111] Figure 4 The left side is the common detection point gather data, and the right side is the corresponding vertical well observation data (VSP record). Figure 3 The model shown in the figure is simulated, so only one event is the primary reflection from the underground interface we need. According to the VSP record, it can be analyzed that the third event from top to bottom is the primary reflection we want to extract, and the remaining events are mostly multiple waves related to the free surface. The first event from top to bottom is a direct wave. Although it is not a multiple wave, it is usually not used for migration imaging.

[0112] The fast Fourier transform algorithm is used to convert the seismic data d(x, t) in the time-space domain into the frequency-space domain d(x, ω). According to the acceleration characteristics of the fast Fourier transform algorithm, the number of time sampling points of the original data is filled to 8192 (the power of 2), and multi-threaded parallel calculation is used for each channel. According to the determined observation system and ray parameters, the common detection point channel data of the two components are plane wave decomposed respectively, and the data order is adjusted.

[0113] The least squares method is used to construct the equation to solve the adaptive filter factor. When solving the equation, the first 700 time sampling points of the plane wave domain data are selected to construct the equation. The filter length is about 50 time sampling points. For the seabed node data with different time sampling intervals, the corresponding filter length can be adjusted according to the time length of the seismic wavelet.

[0114] Figure 5 The schematic diagram of upgoing and downgoing wave separation of the present invention is shown. Upgoing and downgoing wave separation is a necessary prerequisite for subsequent upgoing and downgoing wave deconvolution. The upper left figure is the input seabed node time-space domain common detection point gather data, which contains both upgoing and downgoing waves; the upper right figure is the result of converting the time-space domain data to the plane wave domain; the lower left figure is the separated upgoing wave; and the lower right figure is the separated downgoing wave.

[0115] The separated uplink and downlink waves will be used for the subsequent uplink and downlink deconvolution. The separated plane wave domain uplink and downlink wave data are subjected to one-dimensional fast Fourier transform channel by channel, and the data are converted into the frequency plane wave domain. Then, the uplink and downlink deconvolution is performed data point by data point, and the stability factor of the deconvolution is about 0.001.

[0116] The first reflection Green's function in the plane wave domain obtained by deconvolution is converted back to the time-space domain. The inverse transformation is performed frequency by frequency and needs to cover the maximum frequency sampling range of the data (1 to 4096 frequency slices).

[0117] Figure 6 The left figure is the first reflection Green's function obtained after the up-and-downward wave deconvolution in the plane wave domain; the right figure is the multiple wave suppression result obtained after converting the left figure back to the time-space domain. It can be seen that the left figure is the first reflection Green's function obtained after the up-and-downward wave deconvolution in the plane wave domain. Figure 5 Compared with the original data in the upper left figure, the event axes in the processed results are significantly reduced, the multiple waves in the data are effectively suppressed, while the event axes corresponding to the primary wave are retained, and the resolution is also improved.

[0118] Figure 7 The figure below is a comparative analysis of the processing results of the present invention and the common detection point gathers of the original seafloor node data. The left figure is the common detection point gathers of the seafloor node data, and shows the waveform of the channel where the data receiver is located. The right figure is the deconvolution processing result of the present invention after reference plane correction, and shows the waveform of the channel where the receiver is located. Figure 3In the seafloor node observation system shown, the source and the receiver are located on two different depths of the datum plane, and the up and down wave deconvolution results eliminate this datum plane difference, so the datum plane of the processing result is corrected back to the seafloor node observation system to facilitate comparison with the original data. By comparing with the original data, it can be seen that the third phase axis (i.e., the first reflection) from top to bottom is retained, and the remaining free surface related multiple waves are effectively suppressed. In addition, by comparing the waveform diagram of the extraction channel, it is found that the reflection interface position in the processing result is a zero-phase sub-wave, indicating that the processing result can eliminate the influence of the source sub-wave to a certain extent.

[0119] The present invention is a free surface correlated multiple wave suppression technology for seafloor node data based on a characteristic reflection nonlinear convolution model, and is mainly used for free surface correlated multiple wave suppression of seafloor node observation data in exploration seismic. The method firstly sorts the pressure component and vertical velocity component of the seafloor node data into a common detection point gather, converts the common detection point gather data into a plane wave domain by local plane wave decomposition, merges and separates the double detection data in the plane wave domain to obtain uplink and downlink wave records, then converts the uplink and downlink wave records into the frequency domain, performs uplink and downlink wave deconvolution frequency by frequency and plane wave direction by plane wave direction, and then converts the deconvolution result back to the time-space domain for subsequent processing.

[0120] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solutions of the present invention may be modified or replaced by equivalents without departing from the spirit and scope of the technical solutions of the present invention, which should all be included in the scope of the claims of the present invention.

[0121] It should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that the technical solutions of the present invention may be modified or replaced by equivalents without departing from the spirit and scope of the technical solutions of the present invention, which should all be included in the scope of the claims of the present invention.

[0122] Those skilled in the art will appreciate that the embodiments of the present application can be provided as methods, systems, or computer program products. Therefore, the present application can adopt the form of complete hardware embodiments, complete software embodiments, or embodiments in combination with software and hardware. Moreover, the present application can adopt the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) that contain computer-usable program code. The scheme in the embodiments of the present application can be implemented in various computer languages, for example, object-oriented programming language Java and literal scripting language JavaScript, etc.

[0123] The present application is described with reference to the flowcharts and / or block diagrams of the methods, devices (systems), and computer program products according to the embodiments of the present application. It should be understood that each process and / or box in the flowchart and / or block diagram, as well as the combination of the processes and / or boxes in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing device to generate a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 A process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0124] These computer program instructions may also be stored in a computer-readable memory capable of directing a computer or other programmable data processing device to operate in a specific manner, so that the instructions stored in the computer-readable memory produce an article of manufacture comprising an instruction device, which implements the process Figure 1 A process or multiple processes and / or boxes Figure 1 A function specified in one or more boxes.

[0125] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, thereby providing instructions for implementing the process. Figure 1 A process or multiple processes and / or boxes Figure 1 The steps for the functions specified in one or more boxes.

[0126] Although the preferred embodiments of the present application have been described, those skilled in the art may make other changes and modifications to these embodiments once they have learned the basic creative concept. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments and all changes and modifications falling within the scope of the present application.

[0127] Obviously, those skilled in the art can make various changes and modifications to the present application without departing from the spirit and scope of the present application. Thus, if these modifications and variations of the present application fall within the scope of the claims of the present application and their equivalents, the present application is also intended to include these modifications and variations.

Claims

1. A multiple wave suppression method based on a characteristic reflection nonlinear convolution model, characterized in that: include, Collect multi-component seafloor node data and pre-process the data; Constructing plane wave basis functions according to the preprocessed parameters to convert frequency space domain data into time plane wave domain data; Adaptive matching factors are introduced to achieve adaptive matching of time plane wave domain data and separate upgoing and downgoing seafloor node data; The separated upgoing and downgoing wave seafloor node data are converted into the frequency domain respectively, and the upgoing and downgoing wave deconvolution is performed; The deconvolution result is converted back to the time-space domain to obtain the first-reflection Green's function, that is, the first-reflection seafloor node data, to achieve suppression; The plane wave basis functions include: In the plane wave domain, according to the physical relationship between the magnitudes of the pressure component and the vertical velocity component data, the amplitudes of the two are adjusted to the same magnitude. The magnitude adjustment relationship is as follows: Where P′ represents the pressure component after magnitude adjustment, k Z represents the vertical wave number component, ρ represents the seawater density, ω represents the signal wave number after Fourier transformation, P represents the pressure component before magnitude adjustment, and p X represents the lateral ray parameter component of the local plane wave signal, v sea_water Indicates the speed of sound wave propagation in seawater; The adaptive matching factor includes selecting the direct wave and obtaining it by the least square method. The equation is constructed as follows: Among them, P represents the pressure component data, V Z represents the vertical velocity component data, w represents the best matched filter of the pressure component and the vertical velocity component, Represents the convolution matrix composed of vertical velocity component data, Representation Matrix The transpose of P Represents the matrix composed of pressure component data; The direct wave includes: the trajectory of the direct wave in the time and space domain is expressed as: Where t represents the arrival time curve of the direct wave signal in the time and space domain, that is, the propagation time of the source signal in the seawater before reaching the detector, v sea_water represents the sound wave propagation speed of seawater, z represents the seawater depth where the geophone is placed, x0 represents the transverse spatial coordinate where the geophone is placed, x represents the transverse spatial coordinate where the source is excited, y represents the longitudinal spatial coordinate where the geophone is placed, y0 represents the longitudinal spatial coordinate where the source is excited, Indicates the lateral spatial distance between the earthquake source and the geophone. Indicates the longitudinal spatial distance between the earthquake source and the geophone; The direct wave also includes a direct wave trajectory in the time plane domain expressed as follows: Where t0 represents the arrival time curve of the direct wave in the plane wave domain, p x represents the transverse ray parameter component of the direct wave, p y It represents the longitudinal ray parameter component of the direct wave. The definitions of the transverse and longitudinal ray parameters are as follows: Among them, the physical meaning of each variable is the same as that in the above formula; The upgoing and downgoing seafloor node data include: The formation process of the up and down waves of the seabed node data, that is, the free surface related multiple waves contained therein, is expressed in the plane wave domain as follows: Among them, p represents the ray parameter vector of the signal, ω represents the wave number of the signal after Fourier transformation, and z r represents the seawater depth where the geophone is placed, z0 represents the water depth of the source excitation, exp() represents the natural exponential function, U(p,ω,z r ) represents the upgoing wave signal in the frequency plane wave domain, P(p,ω,z0) represents the full wave field signal observed at the sea surface, i represents the imaginary unit of the complex signal, k z represents the vertical wave number component of the wave field, g(p,ω) represents the primary reflection Green's function, a(ω) represents the wave impedance of the sea surface, s(p,ω) represents the excitation signal of the earthquake source, and D(p,ω,z r ) represents the downlink wave signal in the frequency plane wave domain, P d (p,ω,z r ) represents the direct wave signal.

2. The method for suppressing multiple waves based on the characteristic reflection nonlinear convolution model according to claim 1, characterized in that: The pretreatment includes: Input multi-component seafloor node data into the header keyword program, obtain observation system information, and sort out the common detection point gather data of pressure component and vertical velocity component.

3. The method for suppressing multiple waves based on the characteristic reflection nonlinear convolution model according to claim 2, characterized in that: The observation system information includes the number of time sampling points, the time sampling interval, the spatial position of the observation node relative to the observation system, the excitation position of the sea surface source relative to the observation system and the dip sampling interval of the seismic data.

4. The method for suppressing multiple waves based on a characteristic reflection nonlinear convolution model according to claim 1, characterized in that: The reflection Green's function includes, Among them, P (U) (p,ω,z r ) represents the upgoing wave signal in the frequency plane wave domain, P (D) (p,ω,z r ) represents the downlink signal in the frequency plane wave domain, g(p,ω,z r ) represents the final first reflection Green's function.

5. The method for suppressing multiple waves based on the characteristic reflection nonlinear convolution model according to claim 4, characterized in that: The deconvolution comprises: In order to ensure the stability of the deconvolution process, a stability factor is introduced and the deconvolution process is improved as follows: in, Indicates P (D) (p,ω,z r ), and ε represents a small positive number for stabilizing the solution.

Images