Residual multiple suppression method, device, equipment, medium and program product

Abstract

The application discloses a residual multiple wave suppression method, device, equipment, medium and program product, and the method comprises the steps of: performing NMO correction processing on first common point gather data; traversing different zero offset distances and different parabolic curvatures of the corrected data to obtain a parabolic Radon record, and performing in-phase weighting and absolute value processing to obtain a parabolic Radon spectrum; performing event tracking on the parabolic Radon spectrum to obtain the time position of each multiple wave event in each seismic trace; extracting the first residual multiple wave event by using KL transform filtering, and performing inverse NMO correction processing to obtain the second residual multiple wave event; and eliminating the residual multiple wave by least square filtering based on the first common point gather data and the second residual multiple wave event. According to the technical scheme, the parabolic Radon spectrum is created, and the filtering transform is introduced, so that the multiple wave tracking and suppression are realized, the residual multiple wave is further eliminated, and the signal-to-noise ratio of the effective primary wave is improved.

Description

Technical Field

[0001] This invention relates to the field of seismic exploration technology, and in particular to a method, apparatus, equipment, medium, and program product for suppressing residual multiples. Background Technology

[0002] Currently, the main seismic exploration method used in oil and gas exploration is the reflection wave method, which treats the primary reflected wave as the effective signal, while direct waves, shallow refracted waves, and multiple waves are classified as interference waves. Among them, interference such as direct waves and shallow refracted waves can usually be directly removed from the seismic record; however, multiple waves and effective primary waves overlap in the record and cannot be suppressed by simple removal.

[0003] Current methods for multiple removal mainly include predictive deconvolution, multiple suppression based on apparent velocity differences, and multiple suppression based on wave theory. These methods separate and suppress multiples by utilizing differences in periodicity, apparent velocity, or propagation path between primary and multiple waves. Predictive deconvolution uses the periodicity of multiples to suppress them through deconvolution operations. Methods based on apparent velocity differences use the different distribution regions of primary and multiple waves in the velocity spectrum for filtering and separation. Methods based on wave theory predict multiples by simulating the wave field propagation process and subtract them from the original data. These methods each have their applicable conditions in practical seismic data processing and can all achieve a certain degree of suppression effect. However, due to limitations in complex seafloor morphology, data acquisition conditions, and the accuracy of the methods themselves, some multiple energy remains after processing, making complete removal difficult. Therefore, how to effectively suppress residual multiples is a key problem that urgently needs to be solved in seismic data processing. Summary of the Invention

[0004] This invention provides a method, apparatus, device, medium, and program product for suppressing residual multiples, in order to solve the problem of residual contamination and interference of the effective signal by the residual multiples after the initial suppression, and to improve the signal-to-noise ratio of the effective signal.

[0005] According to one aspect of the present invention, a method for suppressing residual multiples is provided, the method comprising: Dynamic correction processing is performed on the first common center point gather data to obtain the second common center point gather data; For the second common center point gather data, when traversing different zero offset distances and different parabolic curvatures, gathers are superimposed along the parabolic trajectory to obtain parabolic Radon domain records, and the parabolic Radon domain records are weighted in phase and the absolute value is taken to obtain the parabolic Radon spectrum. Parabolic phase axis tracing is performed on the parabolic Lardon spectrum to obtain the temporal position of each multiple phase axis in each seismic trace in the second common center point gather data; The first remaining multiple phase axis is extracted from each multiple phase axis using KL transform filtering. The first residual multiple phase axis is subjected to reaction correction processing to obtain the second residual multiple phase axis; Based on the first common center point gather data and the second residual multiple phase axis, the residual multiple is eliminated by least squares filtering to obtain the third common center point gather data.

[0006] According to another aspect of the present invention, a residual multiple suppression device is provided, the device comprising: The data correction module is used to perform dynamic correction processing on the first common center point gather data to obtain the second common center point gather data; The Radon spectrum creation module is used to perform a parabolic Radon domain record by superimposing the traces along the parabolic trajectory when the second common center point gather data traverses different zero offset distances and different parabolic curvatures, and to perform in-phase weighting on the parabolic Radon domain record and take the absolute value to obtain the parabolic Radon spectrum. The phase axis tracking module is used to perform parabolic phase axis tracking on the parabolic Lardon spectrum to obtain the time position of each multiple phase axis in each seismic trace in the second common center point gather data; The in-phase axis extraction module is used to extract the first remaining multiple in-phase axis from each multiple in-phase axis using KL transform filtering; The reaction correction processing module is used to perform reaction correction processing on the first residual multiple phase axis to obtain the second residual multiple phase axis. The multiple wave cancellation module is used to eliminate the remaining multiple waves by least square filtering based on the first common center point gather data and the second remaining multiple wave phase axis, so as to obtain the third common center point gather data.

[0007] According to another aspect of the present invention, an electronic device is provided, the electronic device comprising: At least one processor; and A memory communicatively connected to the at least one processor; wherein, The memory stores a computer program that can be executed by the at least one processor, the computer program being executed by the at least one processor to enable the at least one processor to perform the residual multiple wave suppression method according to any embodiment of the present invention.

[0008] According to another aspect of the present invention, a computer-readable storage medium is provided, the computer-readable storage medium storing computer instructions for causing a processor to execute and implement the residual multiple wave suppression method according to any embodiment of the present invention.

[0009] According to another aspect of the present invention, a computer program product is provided, the computer program product comprising a computer program that, when executed by a processor, implements the residual multiple suppression method described in any embodiment of the present invention.

[0010] The technical solution of this invention achieves the tracking and suppression of the residual multiple phase axis by creating a high-resolution parabolic Lardon spectrum and introducing filtering transformation. The energy of the residual multiple wave is further attenuated, improving the recognizability of the primary wave, reducing the adverse effects on the authenticity and reliability of seismic imaging, and avoiding misleading subsequent seismic geological interpretation.

[0011] It should be understood that the description in this section is not intended to identify key or essential features of the embodiments of the present invention, nor is it intended to limit the scope of the invention. Other features of the invention will become readily apparent from the following description. Attached Figure Description

[0012] To more clearly illustrate the technical solutions in the embodiments of the present invention, the accompanying drawings used in the description of the embodiments will be briefly introduced below. Obviously, the accompanying drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0013] Figure 1 This is a flowchart of a residual multiple suppression method provided in Embodiment 1 of this application; Figure 2 This is a flowchart of a residual multiple suppression method according to Embodiment 2 of this application; Figure 3 This is a schematic diagram of the first common center point gather data provided according to Embodiment 2 of this application; Figure 4 This is a schematic diagram of the second common center point gather data (a) and the parabolic Radon spectrum (b) provided according to Embodiment 2 of this application; Figure 5 This is a schematic diagram of the second residual multiple phase axis provided according to Embodiment 2 of this application; Figure 6 This is a schematic diagram of the third common center point gather data provided according to Embodiment 2 of this application; Figure 7 This is a schematic diagram of a residual multiple wave suppression device according to Embodiment 3 of this application; Figure 8 This is a schematic diagram of the structure of an electronic device that implements the residual multiple suppression method provided in Embodiment 4 of this application. Detailed Implementation

[0014] To enable those skilled in the art to better understand the present invention, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort should fall within the scope of protection of the present invention.

[0015] It should be noted that the terms "first," "second," etc., in the specification, claims, and accompanying drawings of this invention are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments of the invention described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover a non-exclusive inclusion; for example, a process, method, system, product, or apparatus that comprises a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0016] Example 1 Figure 1 This is a flowchart of a residual multiple suppression method according to Embodiment 1 of this application. This embodiment is applicable to situations where residual multiples affect the identification of effective signals. The method can be executed by a residual multiple suppression device, which can be implemented in hardware and / or software. This residual multiple suppression device can be configured in an electronic device with residual multiple suppression capabilities. For example... Figure 1 As shown, the method includes: S110. Perform dynamic correction processing on the first common center point gather data to obtain the second common center point gather data.

[0017] The first common center point gather data can be a collection of reflected waves from multiple seismic traces recorded when the seismic wave excitation and receiving points are symmetrically arranged about the same center point. Theoretically, reflected waves at different offsets in the reflected wave collection all correspond to the same underground reflection point, enabling the first common center point gather data to reflect the reflection characteristics of that reflection point at different observation distances, thus providing the original data basis for subsequent velocity analysis and stacking imaging.

[0018] For example, in seismic exploration using a multiple coverage observation system, a series of equally spaced excitation and reception points are arranged along the survey line. For a reflection point on a subsurface reflecting interface, multiple pairs of shot points and receiver points correspond to the reflection point, and the center point of each shot point and receiver point combination falls directly above the reflection point. The seismic traces recorded by these multiple shot point and receiver point combinations are extracted and arranged according to their offset, forming the first common center point gather data. Each seismic trace in this first common center point gather data contains reflected wave information from the same reflection point, but the offset and travel time of each trace are different. Each trace also records waveform characteristics such as amplitude, frequency, and phase.

[0019] The second common center point gather data can be understood as the dataset obtained by correcting the time difference of each seismic trace record in the first common center point gather according to the offset. In this second common center point gather data, the phase axis of the primary reflection wave has been flattened into a horizontal straight line, and the arrival times of each seismic trace at the same reflection interface tend to be consistent, giving the effective signal good linear coherence. However, due to inaccurate correction velocity, the multiple waves still retain residual time differences, thus forming time-shift characteristics different from the effective waves within the second common center point gather data, providing a basis for subsequent stacking and wavefield separation.

[0020] For example, the second common center point gather data may contain six seismic traces with offsets from near to far. The arrival times of the reflected waves from the same subsurface interface in each seismic trace are 1.0 seconds, 1.1 seconds, 1.3 seconds, 1.6 seconds, 2.0 seconds, and 2.5 seconds, respectively. After dynamic correction, the arrival times are all corrected to 1.0 seconds. The originally curved phase axis of the primary reflected wave is flattened into a horizontal line, while the phase axis of the multiple waves cannot be completely flattened due to velocity differences and still retains the tilted or curved time difference characteristics.

[0021] The dynamic correction process corrects for the travel time differences caused by varying offsets among seismic traces in a common center point ensemble, restoring them to the travel time of self-excitation and self-reception at the same reflection point. Specifically, since seismic traces with larger offsets have longer propagation paths and later arrival times of reflected waves, dynamic correction uses velocity parameters to adjust the time coordinates of each trace, aligning traces with different offsets in time and thus reflecting the true time location of the underground reflection point. After correction, the phase axes of the reflected waves from each trace are flattened, facilitating subsequent stacking to enhance the effective signal and providing a basis for velocity analysis.

[0022] For example, in the first common center point gather data, seismic traces with different offsets receive reflected waves from the same reflection interface, and the travel time of the seismic traces increases with the offset. The dynamic correction process can be to calculate a time correction amount that varies with the offset for each seismic trace. This correction amount is determined based on the speed of seismic wave propagation and the magnitude of the offset. Then, the sampling time of each seismic trace is subtracted from the corresponding correction amount. After the above trace-by-trace adjustment, the arrival time of the reflected waves, which originally varied with the offset, is uniformly corrected to the self-excitation and self-reception time when the offset is zero, and the phase axis of the reflected waves in each seismic trace changes from curved to straight.

[0023] In this embodiment, dynamic correction processing is performed on the first common center point gather data to correct the travel time of each seismic trace to the self-excitation and self-reception time, so that the phase axis of the reflected wave is straightened from curved to straight, and the seismic traces are consistent in time, eliminating the influence of offset variation on the arrival time of the reflected wave. At the same time, the straightening effect also provides an intuitive basis for velocity analysis.

[0024] S120. For the second common center point gather data, when traversing different zero offset distances and different parabolic curvatures, gathers are superimposed along the parabolic trajectory to obtain parabolic Radon domain records, and the parabolic Radon domain records are weighted in phase and the absolute value is taken to obtain the parabolic Radon spectrum.

[0025] Zero offset time refers to the total time taken for a seismic wave to travel from the excitation point to the underground reflecting interface and back to the receiving point when the excitation point and the receiving point are located at the same location. In this observation method, the seismic wave is incident vertically on the reflecting interface and returns along the same path. Therefore, zero offset time can represent the two-way vertical travel time of the underground reflecting interface and can reflect the depth and shape of the reflecting interface.

[0026] For example, in the second common center point gather data, the offset of the seismic trace is 500 meters, the travel time of the reflected wave is 1.2 seconds, while the zero offset time of the same reflection interface is 1.0 seconds.

[0027] In seismic data processing, the travel time of reflected waves from each seismic trace within the second common center point gather does not exhibit a linear relationship with offset, but rather approximates a parabolic shape. The curvature of this parabola can be understood as its curvature. When the offset is small, the travel time of each seismic trace increases slowly with increasing offset, resulting in a relatively gentle parabola. Conversely, when the offset is large, the travel time increases rapidly with increasing offset, leading to a more pronounced parabolic curvature. The magnitude of the parabolic curvature directly depends on the propagation speed of the seismic wave in the medium. Lower speeds result in greater curvature and a more pronounced bending of the reflected wave's phase axis; higher speeds result in smaller curvature and a straighter reflected wave's phase axis.

[0028] For example, in the second common center point gather data, there are three seismic trace records with offsets of 100 meters, 300 meters, and 500 meters, and the travel times of the reflected waves of each seismic trace can be 1.00 seconds, 1.08 seconds, and 1.20 seconds, respectively. When the data points corresponding to different offsets and different travel times are plotted on a coordinate system with the square of the offset as the horizontal axis and the square of the travel time as the vertical axis, the data points roughly fall on a straight line. However, when converted back to the time-offset coordinate system, the data points are connected to form an upward-curving parabola.

[0029] In seismic data processing, the Parabolic Radon domain record can be understood as a result obtained after the original seismic data has undergone the Parabolic Radon transform. The Parabolic Radon transform superimposes and rearranges the seismic wave energy along parabolic paths with different parabolic curvatures, so that the primary reflection waves and multiple waves that were originally mixed in the time and space domain are significantly distinguished in the Parabolic Radon domain according to the differences in the curvature of the parabolic waves.

[0030] For example, the second common center point gather data can contain seismic wave phase axes with different parabolic curvatures. After undergoing parabolic Ladon transformation, these phase axes are mapped to a new parameter space with zero offset and parabolic curvature as coordinates, forming a parabolic Ladon domain record. In the parabolic Ladon domain, energies with different parabolic curvatures are separated; energies with the same parabolic curvature value are clustered at the same location, while energies with different parabolic curvature values ​​are clustered at different locations. This transforms the overlapping information in the second common center point gather data into a record format that can be identified and processed according to parabolic curvature characteristics.

[0031] Among them, the parabolic Radon spectrum can be a two-dimensional energy distribution map obtained by performing a parabolic Radon transform on the second common center point gather data in seismic data processing. The parabolic Radon spectrum uses zero offset time and parabolic curvature as coordinates to rearrange and focus the seismic wave energy along the parabolic path in the second common center point gather data according to the curvature magnitude. The in-phase axes of different parabolic curvatures will form energy clusters at different curvature positions on the spectrum, thereby showing the curvature characteristics and relative energy intensity of various components in the seismic wave field, providing an intuitive basis for subsequent identification and separation of different wave field components.

[0032] For example, the second common center point gather data may contain three in-phase axes with different curvatures, with parabolic curvatures of 0.5, 1.0, and 2.0, respectively. After the parabolic Radon transform, the energy of the three in-phase axes will be concentrated at the corresponding parabolic curvature positions of 0.5, 1.0, and 2.0 on the parabolic Radon spectrum, forming three independent energy clusters. The amplitude of each energy cluster reflects the relative intensity of different curvature components in the second common center point gather data. By reading the position and amplitude of the energy clusters on the parabolic Radon spectrum, the distribution of different curvature components in the wave field and their energy intensity can be quantitatively determined.

[0033] In forming the parabolic Radon domain record, when traversing different zero offsets and different parabolic curvatures, for each possible curvature parameter value, all sampling points within the entire time range must be scanned along the parabolic path defined by the curvature parameter value. This scanning process needs to be executed one by one for all preset curvature values ​​and all time points to ensure that the energy under each curvature is fully calculated, thereby constructing a complete Radon domain record.

[0034] For example, the curvature values ​​to be scanned can be 0.5, 1.0, and 2.0, and there can be 1000 time sampling points. When forming the Radon domain record, first take curvature 0.5, and scan along the parabolic path corresponding to curvature 0.5 from the first time point to the 1000th time point, taking into account all the data points passed along the path; then take curvature 1.0 and repeat the same scanning process; finally take curvature 2.0 and scan the whole thing again.

[0035] In forming the parabolic Radon domain record, stacking gathers along the parabolic trajectory can be understood as finding all seismic traces and corresponding time positions in the second common center point gather data that satisfy the parabolic equation determined by the curvature parameter and the zero offset, given the curvature parameter and the zero offset. The amplitude values ​​at the time positions corresponding to the curvature parameter and the zero offset are extracted and added together to obtain the total energy value. This total energy value is the Radon transform output result at the position corresponding to the curvature parameter and the zero offset.

[0036] For example, the second common center point gather data can have 5 seismic traces with offsets of 100m, 200m, 300m, 400m, and 500m. For a parabolic curvature of 0.5 and a zero offset of 1.0 second, the time position corresponding to each offset can be calculated using the parabolic equation determined by 0.5 and 1.0. The 5 seismic traces can correspond to 1.05 seconds, 1.20 seconds, 1.45 seconds, 1.80 seconds, and 2.25 seconds, respectively. The amplitude values ​​at the corresponding time points are extracted, and the 5 values ​​are added together to obtain a sum. This sum is the superposition result at the position corresponding to 0.5 curvature and 1.0 second zero offset. The superposition process is repeated for each curvature and each zero offset to complete the construction of the entire parabolic Radon domain record.

[0037] In this process, in-phase weighting of the parabolic Radon domain records can be achieved through adaptive filtering in the transform domain, which enhances the effective wave energy based on signal coherence. The adaptive filtering process involves first transforming the seismic gathers to the parabolic Radon domain, then scanning the energy distribution at each parabolic curvature parameter, and constructing weighting factors based on the phase similarity or amplitude consistency of amplitudes at adjacent points.

[0038] For example, in the parabolic latant domain of the second common center point gather data containing multiple wave interference, the energy of the primary reflected wave converges at a specific parabolic curvature to form a continuous cluster, while multiple waves and random noise are scattered at other locations; in-phase weighting can be achieved by scanning point by point along the parabolic curvature parameter, assigning high weights to sampling points with consistent amplitude and phase between adjacent channels and accumulating them, while assigning low weights to sampling points with disordered phases.

[0039] In the embodiments of this application, by traversing the zero offset time and parabolic curvature of each seismic trace and performing in-phase weighting processing, a high-resolution parabolic Lardon spectrum is obtained, which achieves more refined effective wave energy resolution and effectively avoids energy aliasing between effective waves and interfering waves, thereby better preserving the high-frequency details of the seismic profile while suppressing interference.

[0040] S130. Perform parabolic phase axis tracing on the parabolic Lardon spectrum to obtain the time position of each multiple phase axis in each seismic trace in the second common center point gather data.

[0041] In this context, the phase axis of multiple waves can be understood as the continuous phase line of the reflected waves recorded by the same receiver array after the seismic wave has undergone two or more reflections in the subsurface medium. Unlike primary reflections, multiple waves cannot be completely flattened after dynamic correction in the common midpoint gather data. They still retain varying degrees of tilt, curvature, or residual time difference, often appearing in the record in a layered or oblique pattern parallel to the primary wave, thus interfering with the effective signal.

[0042] For example, in the second common center point gather data, multiple waves from the same reflecting interface undergo two or more reflections. The equivalent velocity corresponding to the propagation path of the multiple waves is lower than that of the primary wave. After dynamic correction, the primary wave is flattened into a horizontal straight line, while the multiple waves retain a positive time difference and present an overall downward curved shape. They always appear later than the effective wave at the same reflecting interface in time.

[0043] The time position of the multiple phase axis can be the sampling time point corresponding to the multiple signal on each seismic trace record. Using the time of the zero-offset seismic trace as a reference, the time position corresponding to seismic traces with different offset distances can be determined by combining the parabolic curvature. The time positions distributed along the parabolic trajectory in the order of offset distance can be connected to form the parabolic phase axis of the multiple.

[0044] For example, in the second common center point gather data, the sampling time point of the multiple wave in the zero offset trace can be a fixed reference. As the offset gradually increases, the sampling time point of the multiple wave in each seismic trace increases synchronously. A series of discretely distributed sampling time points are fitted to obtain the parabolic phase axis of the multiple wave. The curvature of the parabola directly controls the increment of the sampling time point of each seismic trace.

[0045] In this context, in-phase axis tracing refers to identifying and continuously calibrating in-phase axes of seismic waves originating from the same subsurface reflection interface along a seismic profile. This tracing process, based on waveform similarity, phase consistency, and amplitude characteristics, uses manual interpretation or automated algorithms to sequentially connect wave peaks or troughs corresponding to the same reflection interface in discrete seismic traces, forming in-phase reflection axes that characterize the continuity of subsurface strata.

[0046] For example, on a seismic profile, the initial peak position of the reflection interface can be selected as a seed point. Then, a search is performed trace by trace along the time or spatial direction. Waveform characteristics are compared within a specific time window of adjacent seismic traces, and the peak or trough with the highest similarity to the seed point is selected as the tracking point on the seismic trace. This tracking point is then used as a new reference for the next seismic trace search, and this process is iterated until the entire survey line is covered. After each seismic trace is matched, the travel time coordinates of the corresponding tracking point are recorded. Finally, all tracking points are connected according to the seismic trace number to form a complete in-phase axis trajectory of the reflection interface.

[0047] In the embodiments of this application, by tracing the in-phase axis on the parabolic Radon spectrum, the wave fields with different parabolic curvatures are separated on the parabolic Radon spectrum. The primary wave and the multiple waves do not interfere with each other, and the random noise is suppressed because it cannot be coherently superimposed, thereby significantly reducing the tracking difficulty and improving the reliability and accuracy of in-phase axis identification.

[0048] S140. Use KL transform filtering to extract the first remaining multiple phase axis from each multiple phase axis.

[0049] KL transform filtering can be understood as calculating the covariance matrix of the input data and performing eigenvalue decomposition, projecting the original data onto the orthogonal transformation space spanned by eigenvectors, filtering the transformed components according to the magnitude of the eigenvalues, removing low eigenvalue components that represent interference information, retaining the main components that can represent the main information, and then performing an inverse transform to restore the processed data, thereby achieving the separation of effective information and interference information and completing the filtering process.

[0050] For example, KL transform filtering can treat seismic gathers as multi-dimensional samples, obtain uncorrelated principal components through covariance matrix eigenvalue decomposition, concentrate effective signals in principal components with high energy proportions, and distribute noise mainly in low-energy components. After discarding low-energy noise components, inverse transform is performed to reconstruct the gathers, thereby achieving noise suppression.

[0051] The first residual multiple phase axis refers to the residual phase axis formed by the energy of multiple waves that remain in the seismic record after multiple wave suppression processing. This is because the assumptions underlying the suppression algorithm do not fully match the complexity and inhomogeneity of the actual subsurface medium. This results in the energy of multiple waves not being completely attenuated. This residual phase axis is temporally distributed near or after the primary wave phase axis. Although its waveform is similar to the primary wave, its time difference relationship and velocity characteristics do not conform to the propagation law of the primary wave. Therefore, it interferes with the effective signal of the primary wave in subsequent processing, reduces the signal-to-noise ratio and resolution of the stacked profile, and affects the accuracy of velocity analysis and the reliability of geological interpretation.

[0052] For example, in seismic data processing, most of the multiple wave energy can be effectively eliminated after multiple wave suppression. However, in areas with complex geological structures, due to large strata dip angles or drastic changes in the velocity field, multiple waves no longer appear periodically, resulting in some multiple waves not being completely attenuated and remaining in the data. On the processed profile, the residual energy of the multiple waves forms a set of in-phase axes with weak amplitude but similar waveforms at a position approximately twice the travel time below the strong reflection interface.

[0053] In this context, phase axis extraction can be understood as the process of systematically collecting and outputting the travel time data, amplitude values, and phase information of each seismic trace corresponding to the tracked phase axis trajectory from the seismic data into a structured numerical sequence. The extracted data sequence includes the time-depth relationship, waveform characteristic parameters, and confidence index of each seismic trace, which provides quantitative input for velocity modeling, inversion calculation, and geological horizon interpretation. The essence of phase axis extraction is to realize the conversion from profile graphic information to a numerical format that can be used for analysis and calculation.

[0054] For example, after completing the in-phase axis tracking, the time coordinate and waveform amplitude value corresponding to each sampling point along the parabolic trajectory can be read. For each seismic trace, the sampling point with the maximum amplitude or zero phase within the time window of the parabolic trajectory can be selected as the representative value, and recorded one by one according to the original seismic trace number, finally outputting an in-phase axis information file containing three columns of data: trace number, time, and amplitude.

[0055] In the embodiments of this application, by using KL filtering transformation to extract the residual multiple phase axes, accurate travel time and spatial location information can be obtained, which facilitates targeted filtering in subsequent processing, thereby reducing interference to the primary wave signal and improving the clarity and reliability of seismic profile interpretation.

[0056] S150. Perform reaction correction processing on the first residual multiple phase axis to obtain the second residual multiple phase axis.

[0057] In this context, reaction correction can be an inverse transformation operation used in seismic data processing to recover the true propagation path of waves. The reaction correction process takes conventionally dynamically corrected horizontally stacked data as input and, based on velocity field and offset information, reverses the reflection time, corrected to zero offset, to restore the original travel time corresponding to each receiver seismic trace, thereby reconstructing the original kinematic characteristics of the wavefield.

[0058] For example, in the dynamic correction process, the arrival times of reflected waves from each seismic trace within the second common center point gather data are corrected to the arrival times at zero offset, and then stacked to form a horizontal profile. However, before performing pre-stack migration, it is necessary to restore the original time difference information of each seismic trace to meet the input requirements of the migration algorithm. At this time, based on the same velocity field and offset parameters of each seismic trace used in the dynamic correction, a time shift needs to be applied in reverse to the pre-stack corrected second common center point gather data, so that the arrival times of reflected waves from each seismic trace are restored from the zero offset state to the original non-hyperbolic shape before dynamic correction.

[0059] Among them, the second residual multiple phase axis can refer to the residual multiple phase axis after reaction correction. This residual multiple phase axis no longer satisfies the hyperbolic time-distance relationship of normal dynamic correction, and the residual time difference has not been completely eliminated. In the second common center point gather data, it appears as a morphologically distorted phase axis. The apparent velocity of this morphologically distorted phase axis is lower than the apparent velocity of the effective primary reflection wave. It still retains the propagation characteristics of multiple waves and cannot correspond to the reflection of the actual underground strata. It is an interference signal that needs to be further removed in the multiple wave suppression process.

[0060] For example, after the second common center point gather data is subjected to reaction correction, some multiples cannot be completely flattened and still retain residual time differences on the order of tens of milliseconds. The phase axis is bent and distorted, and the apparent velocity is significantly lower than that of the first reflected wave. This distorted phase axis is the remaining multiple phase axis after reaction correction, and its energy needs to be suppressed in subsequent processes.

[0061] In this embodiment of the application, by performing reaction correction processing on the residual multiple phase axis, the true time-distance characteristics of the residual multiple can be preserved, and the residual multiple can be prevented from being incorrectly flattened and mixed into the effective reflection signal. This facilitates the subsequent algorithm to accurately identify the interfering phase axis, distinguish the residual multiple from the formation primary reflection wave, and provide a reliable data basis for subsequent multiple suppression.

[0062] S160. Based on the first common center point gather data and the second residual multiple phase axis, the residual multiple is eliminated by least squares filtering to obtain the third common center point gather data.

[0063] Least square filtering is a mathematical method that extracts the optimal signal from the observed data by minimizing the sum of squared errors. Its core is to minimize the sum of squared residuals between the model output and the actual observed values ​​in order to find the best fitting parameters. It is often used for noise reduction and system identification in signal processing.

[0064] For example, least squares filtering can be implemented by designing a set of filtering operators to minimize the sum of squared errors between the multiple wave estimates output by the filter and the multiple wave components in the actual record. This allows the filter to learn to identify and predict the morphology of the multiple waves. Subtracting the predicted multiple waves from the original seismic record preserves the primary reflection wave signal, which more closely resembles the actual strata.

[0065] Among them, residual multiples can refer to the multiple reflected waves that remain after the previous multiple wave suppression process. These residual multiples originate from the reciprocating reflection of seismic waves between underground interfaces and have not been completely eliminated. They still have the inherent time distance and apparent velocity characteristics of multiples and will be superimposed with the effective reflection signal of the actual underground strata, interfering with the identification of strata reflection information.

[0066] For example, after multiple wave suppression is applied to the first common center point gather data, some low-velocity multiple reflection energy with a low amplitude ratio is still not completely eliminated. This part of the energy is mixed with the primary reflection signal to form the remaining multiple waves.

[0067] Among them, the third common center point gather data refers to the common center point gather after suppressing the remaining multiple waves. The residual multiple reflection interference energy has been removed. The third common center point gather data mainly retains the effective reflection signal generated by the underground strata. The time distance relationship of the phase axis in the third common center point gather data is more in line with the actual geological reflection law in the underground. There is no longer aliasing interference caused by multiple waves between waveforms. The overall signal-to-noise ratio of the third common center point gather data is improved, which can provide more reliable basic gather data for subsequent velocity analysis, overlay imaging and geological interpretation.

[0068] For example, after interference suppression processing, the residual multiple energy with lower amplitude in the first common center point gather data is removed, leaving only the waveform corresponding to the first reflection of the formation, thus obtaining the third common center point gather data after eliminating the residual multiples.

[0069] Eliminating residual multiples can be a signal processing action targeting residual multiple reflection interference. Based on the differences in apparent velocity and time distance characteristics between multiples and primary reflections, the energy of residual multiples mixed in with the effective signal is separated and suppressed. While protecting the effective primary reflection waveform of the formation as much as possible, the interference component is weakened, thereby removing the waveform aliasing problem caused by residual multiple reflections and realizing the processing of residual multiple interference.

[0070] For example, based on the apparent velocity difference between residual multiples and primary reflections in the first common center point gather data, short-time window frequency wavenumber filtering is used to distinguish low-velocity residual multiples from high-velocity effective reflection signals in the frequency wavenumber domain. Energy suppression is applied to the wavenumber region representing the residual multiples, and then the data is transformed back to the seismic trace domain to preserve the primary reflection waveform of the strata as much as possible, thereby completing the elimination of residual multiples.

[0071] In this embodiment, least square filtering is introduced to eliminate residual multiples, effectively avoiding the adverse effects of dynamic correction stretching on common center point gather data, reducing the superposition pollution of effective reflection waveforms by low-velocity interference phase axes, improving the overall signal-to-noise ratio of the gather, and providing more reliable basic data for subsequent seismic stacking, imaging, and geological interpretation.

[0072] This application provides a method for suppressing residual multiples. The method involves performing dynamic correction on first common center point gather data to obtain second common center point gather data; for the second common center point gather data traversing different zero offsets and different parabolic curvatures, gathers are superimposed along the parabolic trajectory to obtain parabolic Ladon domain records; these records are then weighted in phase and their absolute values ​​are taken to obtain the parabolic Ladon spectrum; the parabolic Ladon spectrum is then tracked using parabolic phase axes to obtain the temporal position of each multiple phase axis in each seismic trace within the second common center point gather data; a first residual multiple phase axis is extracted from each multiple phase axis using KL transform filtering; the first residual multiple phase axis is subjected to reaction correction to obtain a second residual multiple phase axis; and based on the first common center point gather data and the second residual multiple phase axis, residual multiples are eliminated using least squares filtering to obtain third common center point gather data. This technical solution achieves multiple wave tracking and suppression by creating a parabolic Lardon spectrum and introducing filtering transformation, further eliminating residual multiple waves and laying a good foundation for subsequent velocity analysis and migration imaging.

[0073] Example 2 Figure 2 This is a flowchart of a residual multiple suppression method according to Embodiment 2 of this application. This embodiment is an optimization based on the above embodiment; schemes not described in detail in this embodiment are described in the above embodiment. Figure 2 As shown, the method in this embodiment of the application specifically includes the following steps: S210. Perform dynamic correction processing on the first common center point gather data to obtain the second common center point gather data.

[0074] For example, the first common center point gather data is ,use This represents the velocity curve of the first wave in the first common center point gather data. Perform dynamic correction to obtain the second common center point gather data after dynamic correction. At this point, the phase axis of the primary wave is corrected to be horizontal, while the phase axes of the remaining multiple waves are downward-curving parabolas. The dynamic correction process can be described as follows: ; in, Offset distance, representing the horizontal distance between the seismic source and the receiving detector; For travel time, it represents the propagation time of a seismic wave from the epicenter to the detector; For self-excitation and self-reception travel time, it represents the propagation time of seismic waves when the source and the detector are in the same position; Record length represents the total time range during which the detector records seismic signals during the seismic acquisition process. Figure 3 This is a schematic diagram of the first common center point gather data provided in Embodiment 2 of this application. Figure 4 (a) is a schematic diagram of the second common center point gather data provided according to Embodiment 2 of this application, as shown in the figure. Figure 3 and Figure 4 As shown in (a), Figure 3 and Figure 4 (a) With travel time as the vertical axis and seismic trace number as the horizontal axis, after dynamic correction processing, the effective reflection phase axis of the common center point gather data is flattened to zero offset time, presenting a horizontal straight line shape, while the multiple waves show downward bending or incomplete flattening due to the low correction speed.

[0075] S220. The zero offset time is cross-combined with the parabolic curvature to form multiple sets of parameters. The amplitude of the second common center point gather data is superimposed according to the parabolic trajectory corresponding to the parameters to obtain the superimposed energy corresponding to the parameters. The superimposed energy is assigned to the corresponding grid position of the Radon spectrum. After all the grid positions are assigned, the parabolic Radon domain record is obtained.

[0076] Wherein, the zero offset time is a preset discrete time point, and the parabolic curvature is a preset discrete sampling variable.

[0077] The superposition energy of the parameters can be obtained by selecting a set of zero-offset time and parabolic curvature parameters, extracting the seismic amplitude corresponding to the position of each seismic trace in the second common center point gather data according to the parabolic trajectory determined by the zero-offset time and parabolic curvature parameters, and accumulating the seismic amplitudes. This energy value can characterize the coherence of the seismic waveform along the parabolic trajectory. For example, when the actual time-distance trajectory of the underground reflected signal matches the current parabolic trajectory, the phase of the waveforms of each seismic trace tends to be consistent, and a high superposition energy will be obtained after accumulation. If the actual time-distance trajectory does not match the current parabolic trajectory, the waveform phases cancel each other out, and the obtained superposition energy will be significantly reduced.

[0078] For example, by iterating through multiple sets of zero-offset time and parabolic curvature parameters, for each set of zero-offset time and parabolic curvature parameters, the seismic amplitude of each seismic trace in the second common center point gather data can be obtained along the corresponding parabolic trajectory and accumulated to calculate the superimposed energy. When the reflected signal trajectory matches the calculated trajectory, the superimposed energy reaches several times the seismic amplitude; when the reflected signal trajectory does not match the calculated trajectory, the superimposed energy drops to the level of background noise. The superimposed energy obtained by combining all zero-offset time and parabolic curvature parameters is arranged sequentially according to the parameter correspondence, thereby obtaining the parabolic Ladon domain record. The formula for calculating the parabolic field record can be described as follows: ; in, When the offset is zero, it represents the propagation time of the seismic wave under self-excitation and self-reception conditions; The parabolic curvature reflects how quickly the distance curve of the common center point gather changes with the increase of the offset distance. These are the lower and upper limits of the offset distance, respectively.

[0079] S230. Based on the first common center point gather data, calculate the in-phase weighting factor of the parabolic curvature using the first formula; calculate the product of the in-phase weighting factor and the parabolic Radon domain record, and take the absolute value to obtain the parabolic Radon spectrum.

[0080] The first formula is described by the following formula: ; in, This is the offset distance. For travel, For the curvature of the parabola, When the offset is zero, The time window length, The number of channels recorded in the offset field. The damping factor, These are the lower and upper limits of the offset distance, respectively. For the first common center point gather data, The in-phase weighting factor is defined as follows: the first common midpoint gather data includes the offset and the travel time; the parabolic Ladon domain record includes the parabolic curvature and the zero offset; the time window length is a preset time window length; the number of traces recorded in the offset domain is the sum of the number of seismic traces in the first common midpoint gather data; the damping factor is a preset damping factor; the lower limit value is a preset lower limit value; and the upper limit value is a preset upper limit value.

[0081] Specifically, when a phase axis with parabolic curvature exists in the second common center point gather data, the phase weighting factor corresponding to that parabolic curvature is larger; otherwise, it is smaller or zero. The time window length represents the time span of a segment of data extracted on the time axis. The damping factor is generally taken as 0.01 to 0.001 of the average amplitude of the second common center point gather data, and both are determined based on the actual application.

[0082] For example, the parabolic Radon domain record is weighted in phase based on the in-phase weighting factor, and the absolute value is taken to obtain the parabolic Radon spectrum. The calculation process can be described as follows: ; In order to avoid mistracking of the effective wave phase axis, the parabolic Ladon spectrum is used. Chinese correspondence All energy values ​​are set to zero. Figure 4 (b) A schematic diagram of the parabolic Radon spectrum provided in Embodiment 2 of this application, as shown in... Figure 4 As shown in (b), with zero offset as the vertical axis and parabolic curvature as the horizontal axis, the energy is focused into several discrete extreme points. The energy of each cluster structure corresponds to a phase axis. The more concentrated the energy, the better the fit of the phase axis.

[0083] S240. For the extreme point corresponding to the parabolic curvature value and the zero offset in the second common center point gather data, form the cluster structure energy of the extreme point in the parabolic Radon spectrum, search for the time position of the extreme point, and fit each multiple phase axis in the second common center point gather data according to the time position.

[0084] In the parabolic Ladon spectrum, an extremum point can be the location where the energy value on the spectral surface reaches a local maximum or minimum, representing the optimal match between the second common center point gather data and the parabolic trajectory determined by that extremum point. When the superposition energy of a set of parabolic curvature parameters and zero offset is the highest, it indicates that the parabola corresponding to that set of parabolic curvature parameters and zero offset parameters can accurately describe the linear or curved characteristics in the second common center point gather data. The extremum points corresponding to this set of parabolic curvature parameters and zero offset parameters directly reflect the implicit pattern and structural information in the data.

[0085] For example, in a parabolic Radon spectrum, an extremum can be a local peak on the spectral surface where the energy value is significantly higher than the surrounding values. This indicates that the parabolic trajectory determined by the extremum best matches the signal in the second common center point gather data. Each extremum corresponds to an actual existing in-phase axis. The coordinate position of the extremum precisely indicates the arrival time and curvature of the in-phase axis, while the energy amplitude of the extremum reflects the relative strength of the in-phase axis.

[0086] The clumping energy at extreme points can be understood as energy not concentrated in an isolated peak in the transform domain or parameter space, but rather clustered within a local region, forming clumps or patches. This clumping structure typically arises because the features of the second common center point gather data exhibit a certain range of variation or uncertainty in the parameter dimensions, causing the energy to diffuse along a specific direction. The centroid of the clumping structure usually represents the best estimated parameters of the second common center point gather data features, while the width of the clumping energy distribution on the spectral surface reflects the magnitude of variation in the feature parameters or the degree of ambiguity in the data.

[0087] For example, in a parabolic Ladon spectrum, the energy clusters at extreme points can manifest as small patches of energy with a continuous distribution within a local region, rather than a single, isolated peak. When the parabolic curvature of the in-phase axis varies slightly with spatial location, the energy of that in-phase axis will not concentrate on a single peak. Not a single point, but several adjacent points. Value and The energy diffuses between the values, forming an elliptical or strip-shaped high-energy region, where the energy amplitude at different locations within the region gradually decreases from the center outwards.

[0088] The search for extreme points can be carried out by using an algorithm to traverse all locations on a given data space or parameter plane to find the coordinates where the energy value reaches a local maximum or minimum. Searching for extreme points typically involves comparing the values ​​of each point with those of its neighbors, filtering out locations that are significantly higher or lower than the surrounding area, thereby locating the most prominent feature or the strongest response in the data.

[0089] For example, searching for extrema in a parabolic Ladon spectrum can involve traversing the entire two-dimensional spectral plane, calculating the energy difference between each location and its neighboring points, and identifying local peak locations where the energy is significantly higher than the neighborhood. When there is a single phase axis, a cluster of energy structures will appear on the spectral surface; the highest energy point at the center of this cluster is marked as the extremum. When multiple phase axes exist, several dispersed extrema need to be identified, each corresponding to the best-fit parameters of a single phase axis. Regions with lower energy or gradual changes will not be marked.

[0090] Among them, fitting the phase axis of multiple waves can be based on the propagation law of multiple waves in the recorded signal, using a preset trajectory model to match and approximate the energy trajectory of multiple waves, adjusting the trajectory model parameters to make the model trajectory as aligned as possible with the changing trend of the phase axis of multiple waves in the actual signal, characterizing the spatiotemporal distribution characteristics of multiple waves, and providing support for subsequent separation processing of multiple waves.

[0091] For example, in a parabolic Radon spectrum, the least squares method can be used to fit the phase axis of the multiple waves. For the established energy clusters of the multiple waves within the spectral domain, a model is fitted based on a standard parabola. By iteratively fine-tuning the curvature and zero offset of the parabola, the sum of squared residuals between the fitted trajectory nodes and the actual extreme points of the multiple wave energy is continuously calculated. The fitting process can use minimizing the sum of squared residuals as the optimal criterion, continuously correcting the model parameters until the standard parabolic trajectory accurately matches the distribution trend of the multiple wave energy, thus completing the fitting of the phase axis of the multiple waves and precisely locking the complete shape and characteristic parameters of the phase axis.

[0092] For example, a pre-defined sliding window can be used to perform sliding window averaging on the parabolic Radon spectrum, and then the parabolic curvature value in the second common center point gather data can be used. And at zero offset The in-phase axis forms a pattern in the parabolic Radon spectrum. For the energy of the clumping structure with the central extremum, the extreme point is searched within the sliding window. The time location, and then based on the coordinates of that extreme point. The corresponding in-phase axis in the second common center point gather data is fitted. A specific example of the calculation formula for fitting the in-phase axis can be described as follows: ; in, This indicates the trace number of each seismic trace in the second common center point gather data. This is the offset of the seismic trace data. Extreme point When traveling through various earthquake zones.

[0093] S250. For seismic traces with trace numbers greater than or equal to a trace number threshold, a multiple wave recording segment of a short time window length is extracted from the seismic trace centered on the travel time of the trace number; based on the multiple wave recording segment, the multiple wave phase axis recording corrected to horizontal is extracted by KL transform filtering; the multiple wave phase axis recording is reverse rearranged and put back into the seismic trace to obtain the first residual multiple wave phase axis.

[0094] Wherein, the channel number threshold is a preset channel number threshold, and the short window length is a preset short window length.

[0095] The travel time for each trace can refer to the time it takes for the seismic signal corresponding to different seismic traces to propagate from the excitation location to the receiving location. For example, the travel time for a near-offset trace can be 1.2 seconds, while the travel time for a far-offset trace, which has a longer propagation path, increases to 1.4 seconds.

[0096] In this context, a multiple wave recording segment can be a data interval in a seismic record occupied by signals formed by multiple reflected waves. For example, a single effective wave may appear around 1.2 seconds on the time axis, and the corresponding multiple wave recording segment is distributed within the data range of 2.4 seconds to 2.7 seconds.

[0097] Among them, multiple phase axis records can be signal records composed of the continuous arrangement of energy points formed by multiple waves in a seismic record. For example, within the time interval of 2.4 seconds to 2.7 seconds on the time axis, the energy points of multiple waves from each seismic trace are arranged in an orderly and continuous manner, forming a segment of identifiable multiple phase axis records.

[0098] In this context, de-rearrangement can be the reverse processing operation that restores the rearranged signal data to its original data arrangement. For example, the original data can be numbered 1, 2, 3, 4 in sequence according to the channel numbers. After rearrangement, the data order is adjusted to 1, 3, 2, 4. Performing de-rearrangement processing can restore the rearranged sequence back to the original data arrangement with the channel numbers 1, 2, 3, 4 in sequence.

[0099] For example, given a channel number threshold Regarding the Taoist name Earthquake routes, with travel times corresponding to each route number. A pre-defined short time window length of multiple waveform recordings is extracted from the center, and the multiple waveform recordings are aligned along the starting position to correct the target phase axis to horizontal. Using the extracted multiple waveform recordings as input, the horizontally corrected multiple waveform phase axis recordings are extracted using the KL transform filtering method. The extracted multiple waveform phase axis recordings are then stitched together to obtain the first remaining multiple waveform phase axis. .

[0100] S260. Perform reaction correction processing on the first residual multiple phase axis to obtain the second residual multiple phase axis.

[0101] For example, the reaction correction process can be described as follows: ; in, It is the residual multiple phase axis after reaction correction, i.e., the second residual multiple phase axis. Figure 5 This is a schematic diagram of the second residual multiple phase axis provided in Embodiment 2 of this application, as shown below. Figure 5 As shown, Figure 5 With travel time as the vertical axis and seismic trace number as the horizontal axis, after reaction correction, the remaining multiple wave phase axes will appear to be bent downwards, forming a significant difference from the flattened primary wave, and will appear as tilted or twisted interference energy on the profile.

[0102] S270, Based on the first common center point gather data and the first residual multiple phase axis, through The adaptive subtraction operation of the norm determines the filter factor.

[0103] in, The norm can characterize the overall magnitude of a vector. It can be calculated by taking the square of each element of the vector, summing the sums, and then taking the square root of the sum. For example, In seismic data processing, the norm can be used to quantitatively measure the overall deviation between two sets of data. The seismic residual vector can be composed of four elements: 2, -4, 3, and -1. The corresponding norm for this seismic residual vector can then be calculated. The norm is approximately 5.48. The larger the norm, the higher the overall residual deviation between the two sets of data.

[0104] Adaptive subtraction can be a processing method that automatically adjusts the subtraction weights based on the characteristics of the data itself, subtracting the reference signal from the input signal to suppress unwanted interference signals. For example, for seismic records, the amplitude and phase of each seismic trace are automatically adapted, and the estimated multiple wave signals are subtracted from the seismic records to reduce multiple wave interference.

[0105] Here, the filter factor can be a set of coefficients used to perform filtering operations. By performing operations with the input data, the processing objective of suppressing useless signals and retaining the target signal is achieved. For example, a set of filter factors is randomly selected and convolved with the seismic trace data. This set of filter factors is used to reduce noise components and retain effective seismic reflection signals.

[0106] The adaptive subtraction operation includes: subtracting the product of the first residual multiple phase axis and the filter factor from the first common center point gather data to obtain a difference, and then performing... Norm calculation yields the sum of squared error energies.

[0107] For example, using The process of determining the filter factor using norm adaptive subtraction can be described as follows: ; in, The sum of squares of error energy; The filter factor is *, which indicates convolution operation. Indicates based on The least squares constraint process for norms.

[0108] S280. When the sum of squared error energy is less than the preset sum of squared error energy and the length of the filter factor is 1, the filter factor is calculated using the second formula.

[0109] The second formula is described using the following formula: ; in, This is the offset distance. For travel, This represents a filter factor of length 1. This represents the matrix transpose operation. This represents the first common center point gather data. This indicates the first remaining multiple wave phase axis.

[0110] The preset sum of squared error energy is determined based on the actual application. Considering the computational cost in real-world scenarios, only a filter factor of length 1 is used.

[0111] S290. Based on the first common center point gather data, the first residual multiple phase axis, and the filter factor, the third common center point gather data is calculated using the third formula.

[0112] The third formula is described using the following formula: ; in, For the third common center point gather data, For the first common center point gather data, is the in-phase axis of the first remaining multiple waves, and * indicates convolution operation. Figure 6 This is a schematic diagram of the third common center point gather data provided in Embodiment 2 of this application, as shown below. Figure 6 As shown, Figure 6 Using travel time as the vertical axis and earthquake track number as the horizontal axis, compare... Figure 3 and Figure 6 It can be observed Figure 6 Multiple waves in the middle and near channels were basically eliminated.

[0113] This application provides a residual multiple suppression method, which involves performing dynamic correction processing on first common center point gather data to obtain second common center point gather data; cross-combining the zero offset time with the parabolic curvature to form multiple sets of parameters; performing amplitude superposition on the second common center point gather data according to the parabolic trajectory corresponding to the parameters to obtain the superposition energy corresponding to the parameters; assigning the superposition energy to the corresponding grid position of the Radon spectrum; obtaining the parabolic Radon domain record after all grid positions are assigned; calculating the in-phase weighting factor of the parabolic curvature using a first formula based on the first common center point gather data; calculating the product of the in-phase weighting factor and the parabolic Radon domain record, and taking the absolute value to obtain the parabolic Radon spectrum; and targeting the parabolic curvature value and the zero offset time in the second common center point gather data. The extreme point corresponding to the offset is identified in the parabolic Lardon spectrum, forming a cluster structure energy at the extreme point. The time position of the extreme point is searched, and each multiple phase axis in the second common center point gather data is fitted based on the time position. For seismic traces with trace numbers greater than or equal to a trace number threshold, a short time window length of multiple recording segment is extracted from the seismic trace centered on the travel time of the trace number. Based on the multiple recording segment, the corrected horizontal multiple phase axis record is extracted through KL transform filtering. The multiple phase axis record is reverse rearranged and placed back into the seismic trace to obtain the first residual multiple phase axis. The first residual multiple phase axis is subjected to reaction correction processing to obtain the second residual multiple phase axis. Based on the first common center point gather data and the first residual multiple phase axis, through... The adaptive subtraction operation of the norm determines the filtering factor; when the sum of squared error energy is less than the preset sum of squared error energy and the length of the filtering factor is 1, the filtering factor is calculated using the second formula; based on the first common center point gather data, the first residual multiple phase axis, and the filtering factor, the third common center point gather data is calculated using the third formula. This technical solution creates a high-resolution parabolic Lardon spectrum by introducing a phase weighting factor, enabling the tracking of the residual multiple phase axis. It also uses KL transform filtering to extract the residual multiple phase axis, and utilizes the filtering factor to eliminate the residual multiple phase axis, thereby improving the signal-to-noise ratio and fidelity of the primary reflected wave, making the effective signal clearer and more prominent, and reducing interference and distortion of the effective wave's phase axis by the multiple waves.

[0114] Example 3 Figure 7This is a schematic diagram of a residual multiple suppression device according to Embodiment 3 of this application. This embodiment is applicable to situations where residual multiples affect the identification of effective signals. The residual multiple suppression device can be implemented in hardware and / or software, and can be configured in electronic devices with residual multiple suppression capabilities. Figure 7 As shown, the device includes: Data correction module 310 is used to perform dynamic correction processing on the first common center point gather data to obtain the second common center point gather data; The Radon spectrum creation module 320 is used to perform a parabolic Radon domain record by superimposing the traces along the parabolic trajectory when the second common center point gather data traverses different zero offset distances and different parabolic curvatures, and to perform in-phase weighting on the parabolic Radon domain record and take the absolute value to obtain the parabolic Radon spectrum. The phase axis tracking module 330 is used to perform parabolic phase axis tracking on the parabolic Radon spectrum to obtain the time position of each multiple phase axis in each seismic trace in the second common center point gather data; The in-phase axis extraction module 340 is used to extract the first residual multiple in-phase axis from each multiple in-phase axis using KL transform filtering; The reaction correction processing module 350 is used to perform reaction correction processing on the first residual multiple phase axis to obtain the second residual multiple phase axis. The multiple cancellation module 360 ​​is used to eliminate the remaining multiples by least square filtering based on the first common center point gather data and the second remaining multiple phase axis to obtain the third common center point gather data.

[0115] In this embodiment of the application, the Radon spectrum creation module 320 includes: The Radon domain record creation unit is used to cross-combine the zero offset time with the parabolic curvature to form multiple sets of parameters, and to perform amplitude superposition on the second common center point gather data according to the parabolic trajectory corresponding to the parameters to obtain the superposition energy corresponding to the parameters. The superposition energy is then assigned to the corresponding grid position of the Radon spectrum. After all grid positions are assigned, the parabolic Radon domain record is obtained. Wherein, the zero offset time is a preset discrete time point, and the parabolic curvature is a preset discrete sampling variable. The in-phase weighting factor calculation unit is used to calculate the in-phase weighting factor of the parabolic curvature based on the first common center point gather data and using a first formula; wherein, the first formula is described by the following formula: ; in, This is the offset distance. For travel, For the curvature of the parabola, When the offset is zero, The time window length, The number of channels recorded in the offset field. The damping factor, These are the lower and upper limits of the offset distance, respectively. For the first common center point gather data, The phase weighting factor is defined as follows: the first common midpoint gather data includes the offset and the travel time; the parabolic Ladon domain record includes the parabolic curvature and the zero offset; the time window length is a preset time window length; the number of traces recorded in the offset domain is the sum of the number of seismic traces in the first common midpoint gather data; the damping factor is a preset damping factor; the lower limit value is a preset lower limit value; and the upper limit value is a preset upper limit value. The Radon spectrum calculation unit is used to calculate the product of the in-phase weighting factor and the parabolic Radon domain record, and take the absolute value to obtain the parabolic Radon spectrum.

[0116] In this embodiment, the co-phase axis tracking module 330 includes: The multiple phase axis fitting unit is used to search for the extremum point corresponding to the parabolic curvature value and the zero offset in the second common center point gather data, to form the cluster structure energy of the extremum point in the parabolic Radon spectrum, to find the time position of the extremum point, and to fit each multiple phase axis in the second common center point gather data according to the time position.

[0117] In this embodiment of the application, the co-phase axis extraction module 340 includes: The multiple wave recording segment extraction unit is used to extract a short-time window length of multiple wave recording segment from the seismic trace with the trace number greater than or equal to the trace number threshold, centered on the travel time of the trace number; wherein, the trace number threshold is a preset trace number threshold, and the short-time window length is a preset short-time window length. The multiple wave in-phase axis recording extraction unit is used to extract the corrected horizontal multiple wave in-phase axis recording based on the multiple wave recording segment through KL transform filtering. The first residual multiple phase axis acquisition unit is used to reverse rearrange the multiple phase axis records and put them back into the seismic trace to obtain the first residual multiple phase axis.

[0118] In this embodiment of the application, the multiple wave cancellation module 360 ​​includes: The filter factor determination unit is used to determine the filter factor based on the first common center point gather data and the first residual multiple phase axis, through... An adaptive subtraction operation of the norm determines the filter factor; wherein, the adaptive subtraction operation includes: subtracting the product of the first residual multiple phase axis and the filter factor from the first common center point gather data to obtain a difference, and performing a step on the difference. Norm calculation yields the sum of squared error energies; The filter factor calculation unit is used to calculate the filter factor using a second formula when the sum of squared error energy is less than a preset sum of squared error energy and the length of the filter factor is 1; wherein the second formula is described by the following formula: ; in, This is the offset distance. For travel, This represents a filter factor of length 1. This represents the matrix transpose operation. This represents the first common center point gather data. Indicates the first residual multiple wave phase axis; The third common center point gather data calculation unit is used to calculate the third common center point gather data based on the first common center point gather data, the first residual multiple phase axis, and the filter factor using a third formula; wherein, the third formula is described by the following formula: ; in, For the third common center point gather data, For the first common center point gather data, is the in-phase axis of the first remaining multiple waves, and * indicates convolution operation.

[0119] The residual multiple suppression device provided in this embodiment of the invention can execute the residual multiple suppression method provided in any embodiment of the invention, and has the corresponding functional modules and beneficial effects of the method.

[0120] Example 4 Figure 8 This is a structural diagram of an electronic device that implements the residual multiple suppression method provided in Embodiment 4 of this application, as shown below. Figure 8The diagram illustrates a schematic representation of an electronic device 10 that can be used to implement embodiments of the present invention. The electronic device is intended to represent various forms of digital computers, such as laptop computers, desktop computers, workstations, personal digital assistants, servers, blade servers, mainframe computers, and other suitable computers. The electronic device can also represent various forms of mobile devices, such as personal digital processors, cellular phones, smartphones, wearable devices (e.g., helmets, glasses, watches, etc.), and other similar computing devices. The components shown herein, their connections and relationships, and their functions are merely illustrative and are not intended to limit the implementation of the invention described and / or claimed herein.

[0121] like Figure 8 As shown, the electronic device 10 includes at least one processor 11 and a memory, such as a read-only memory (ROM) 12, a random access memory (RAM) 13, etc., which is communicatively connected to the at least one processor 11. The memory stores a computer program that can be executed by the at least one processor 11, and the computer program is executed by the at least one processor 11 to enable the at least one processor 11 to perform the method provided by the present invention.

[0122] The processor 11 can perform various appropriate actions and processes based on a computer program stored in the read-only memory (ROM) 12 or a computer program loaded from the storage unit 18 into the random access memory (RAM) 13. The RAM 13 can also store various programs and data required for the operation of the electronic device 10. The processor 11, ROM 12, and RAM 13 are interconnected via a bus 14. An input / output (I / O) interface 15 is also connected to the bus 14.

[0123] Multiple components in electronic device 10 are connected to I / O interface 15, including: input unit 16, such as keyboard, mouse, etc.; output unit 17, such as various types of displays, speakers, etc.; storage unit 18, such as disk, optical disk, etc.; and communication unit 19, such as network card, modem, wireless transceiver, etc. Communication unit 19 allows electronic device 10 to exchange information / data with other devices through computer networks such as the Internet and / or various telecommunications networks.

[0124] Processor 11 can be a variety of general-purpose and / or special-purpose processing components with processing and computing capabilities. Some examples of processor 11 include, but are not limited to, a central processing unit (CPU), a graphics processing unit (GPU), various special-purpose artificial intelligence (AI) computing chips, various processors running machine learning model algorithms, digital signal processors (DSPs), and any suitable processor, controller, microcontroller, etc. Processor 11 performs the various methods and processes described above, such as the methods provided in this invention.

[0125] In some embodiments, the methods provided herein may be implemented as a computer program tangibly contained in a computer-readable storage medium, such as storage unit 18. In some embodiments, part or all of the computer program may be loaded and / or installed on electronic device 10 via ROM 12 and / or communication unit 19. When the computer program is loaded into RAM 13 and executed by processor 11, one or more steps of the methods described above may be performed. Alternatively, in other embodiments, processor 11 may be configured to execute the methods by any other suitable means (e.g., by means of firmware).

[0126] Various embodiments of the systems and techniques described above herein can be implemented in digital electronic circuit systems, integrated circuit systems, field programmable gate arrays (FPGAs), application-specific integrated circuits (ASICs), application-specific standard parts (ASSPs), systems-on-chip (SoCs), complex programmable logic devices (CPLDs), computer hardware, firmware, software, and / or combinations thereof. These various embodiments may include implementations in one or more computer programs that can be executed and / or interpreted on a programmable system including at least one programmable processor, which may be a dedicated or general-purpose programmable processor, capable of receiving data and instructions from a storage system, at least one input device, and at least one output device, and transmitting data and instructions to the storage system, the at least one input device, and the at least one output device.

[0127] Computer programs used to implement the methods of the present invention may be written in any combination of one or more programming languages. These computer programs may be provided to a processor of a general-purpose computer, a special-purpose computer, or other programmable data processing device, such that when executed by the processor, the computer programs cause the functions / operations specified in the flowcharts and / or block diagrams to be performed. The computer programs may be executed entirely on a machine, partially on a machine, or as a standalone software package, partially on a machine and partially on a remote machine, or entirely on a remote machine or server.

[0128] In the context of this invention, a computer-readable storage medium stores computer instructions that are used to cause a processor to execute and implement the method provided by this invention.

[0129] The present invention also provides a computer program product comprising a computer program that, when executed by a processor, implements the method provided according to embodiments of the present invention. A computer-readable storage medium may be a tangible medium that may contain or store a computer program for use by or in conjunction with an instruction execution system, apparatus, or device. The computer-readable storage medium may include, but is not limited to, electronic, magnetic, optical, electromagnetic, infrared, or semiconductor systems, apparatus, or devices, or any suitable combination thereof. Alternatively, the computer-readable storage medium may be a machine-readable signal medium. More specific examples of machine-readable storage media include electrical connections based on one or more wires, portable computer disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM or flash memory), optical fibers, portable compact disk read-only memory (CD-ROM), optical storage devices, magnetic storage devices, or any suitable combination thereof.

[0130] To provide interaction with a user, the systems and techniques described herein can be implemented on an electronic device having: a display device (e.g., a CRT (cathode ray tube) or LCD (liquid crystal display) monitor) for displaying information to the user; and a keyboard and pointing device (e.g., a mouse or trackball) through which the user provides input to the electronic device. Other types of devices can also be used to provide interaction with the user; for example, feedback provided to the user can be any form of sensory feedback (e.g., visual feedback, auditory feedback, or tactile feedback); and input from the user can be received in any form (including sound input, voice input, or tactile input).

[0131] The systems and technologies described herein can be implemented in computing systems that include backend components (e.g., as data servers), or middleware components (e.g., application servers), or frontend components (e.g., user computers with graphical user interfaces or web browsers through which users can interact with implementations of the systems and technologies described herein), or any combination of such backend, middleware, or frontend components. The components of the system can be interconnected via digital data communication of any form or medium (e.g., communication networks). Examples of communication networks include local area networks (LANs), wide area networks (WANs), blockchain networks, and the Internet.

[0132] A computing system can include clients and servers. Clients and servers are generally located far apart and typically interact through communication networks. The client-server relationship is created by computer programs running on the respective computers and having a client-server relationship with each other. The server can be a cloud server, also known as a cloud computing server or cloud host, which is a hosting product within the cloud computing service system to address the shortcomings of traditional physical hosts and VPS services, such as high management difficulty and weak business scalability.

[0133] This invention also provides a computer program product, including a computer program that, when executed by a processor, can implement the methods provided in any embodiment of this application.

[0134] In the implementation of the computer program product, computer program code for performing the operations of this application can be written in one or more programming languages ​​or a combination thereof. Programming languages ​​include object-oriented programming languages ​​such as Java, Smalltalk, and C++, as well as conventional procedural programming languages ​​such as C or similar languages. The program code can be executed entirely on the user's computer, partially on the user's computer, as a standalone software package, partially on the user's computer and partially on a remote computer, or entirely on a remote computer or server. In cases involving remote computers, the remote computer can be connected to the user's computer via any type of network—including a local area network (LAN) or a wide area network (WAN)—or can be connected to an external computer (e.g., via the Internet using an Internet service provider).

[0135] It should be understood that the various forms of processes shown above can be used, with steps reordered, added, or deleted. For example, the steps described in this invention can be executed in parallel, sequentially, or in different orders, as long as the desired result of the technical solution of this invention can be achieved, and this is not limited herein.

[0136] The specific embodiments described above do not constitute a limitation on the scope of protection of this invention. Those skilled in the art should understand that various modifications, combinations, sub-combinations, and substitutions can be made according to design requirements and other factors. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of this invention should be included within the scope of protection of this invention.

Claims

1. A method for suppressing residual multiples, characterized in that, The method includes: Dynamic correction processing is performed on the first common center point gather data to obtain the second common center point gather data; For the second common center point gather data, when traversing different zero offset distances and different parabolic curvatures, gathers are superimposed along the parabolic trajectory to obtain parabolic Radon domain records, and the parabolic Radon domain records are weighted in phase and the absolute value is taken to obtain the parabolic Radon spectrum. Parabolic phase axis tracing is performed on the parabolic Lardon spectrum to obtain the temporal position of each multiple phase axis in each seismic trace in the second common center point gather data; The first remaining multiple phase axis is extracted from each multiple phase axis using KL transform filtering. The first residual multiple phase axis is subjected to reaction correction processing to obtain the second residual multiple phase axis; Based on the first common center point gather data and the second residual multiple phase axis, the residual multiple is eliminated by least squares filtering to obtain the third common center point gather data.

2. The method according to claim 1, characterized in that, For the second common center point gather data traversing different zero offset distances and different parabolic curvatures, gather superposition along the parabolic trajectory yields the parabolic Ladon domain record, including: The zero offset time is cross-combined with the parabolic curvature to form multiple sets of parameters. The amplitude of the second common center point gather data is superimposed according to the parabolic trajectory corresponding to the parameters to obtain the superimposed energy corresponding to the parameters. The superimposed energy is assigned to the corresponding grid position of the Radon spectrum. After all grid positions are assigned, the parabolic Radon domain record is obtained. Wherein, the zero offset time is a preset discrete time point, and the parabolic curvature is a preset discrete sampling variable.

3. The method according to claim 2, characterized in that, The parabolic Radon domain record is weighted in phase, and the absolute value is taken to obtain the parabolic Radon spectrum, including: Based on the first common center point gather data, the in-phase weighting factor of the parabolic curvature is calculated using the first formula; Calculate the product of the in-phase weighting factor and the parabolic Radon domain record, and take the absolute value to obtain the parabolic Radon spectrum; The first formula is described by the following formula: ； in, This is the offset distance. For travel, For the curvature of the parabola, When the offset is zero, The time window length, The number of channels recorded in the offset field. The damping factor, These are the lower and upper limits of the offset distance, respectively. For the first common center point gather data, The in-phase weighting factor is defined as follows: the first common midpoint gather data includes the offset and the travel time; the parabolic Ladon domain record includes the parabolic curvature and the zero offset; the time window length is a preset time window length; the number of traces recorded in the offset domain is the sum of the number of seismic traces in the first common midpoint gather data; the damping factor is a preset damping factor; the lower limit value is a preset lower limit value; and the upper limit value is a preset upper limit value.

4. The method according to claim 1, characterized in that, Parabolic phase axis tracing is performed on the parabolic Lardon spectrum to obtain the temporal position of each multiple phase axis in each seismic trace in the second common center point gather data, including: For the extreme point corresponding to the parabolic curvature value and the zero offset in the second common center point gather data, the cluster structure energy of the extreme point is formed in the parabolic Radon spectrum, the time position of the extreme point is searched, and the phase axis of each multiple wave in the second common center point gather data is fitted according to the time position.

5. The method according to claim 1, characterized in that, Extracting the first residual multiple phase axis from each multiple phase axis using KL transform filtering includes: For seismic traces with trace numbers greater than or equal to a trace number threshold, a multiple wave recording segment of a short time window length is extracted from the seismic trace centered on the travel time of the trace number; wherein, the trace number threshold is a preset trace number threshold, and the short time window length is a preset short time window length. Based on the aforementioned multiple wave recording segments, the multi-wave in-phase axis recordings that have been corrected to be horizontal are extracted by KL transform filtering. The multiple phase axis records are reverse rearranged and placed back into the seismic trace to obtain the first residual multiple phase axis.

6. The method according to claim 1, characterized in that, Based on the first common center point gather data and the second residual multiple phase axis, the residual multiples are eliminated by least squares filtering to obtain the third common center point gather data, including: Based on the first common center point gather data and the first residual multiple phase axis, through An adaptive subtraction operation of the norm determines the filter factor; wherein, the adaptive subtraction operation includes: subtracting the product of the first residual multiple phase axis and the filter factor from the first common center point gather data to obtain a difference, and performing a step on the difference. Norm calculation yields the sum of squared error energies; When the sum of squared error energies is less than a preset sum of squared error energies and the length of the filter factor is 1, the filter factor is calculated using the second formula; wherein, the second formula is described by the following formula: ； in, This is the offset distance. For travel, This represents a filter factor of length 1. This represents the matrix transpose operation. This represents the first common center point gather data. Indicates the first residual multiple wave phase axis; Based on the first common center point gather data, the first residual multiple phase axis, and the filter factor, the third common center point gather data is calculated using the third formula; wherein, the third formula is described by the following formula: ； in, For the third common center point gather data, For the first common center point gather data, is the in-phase axis of the first remaining multiple waves, and * indicates convolution operation.

7. A residual multiple wave suppression device, characterized in that, The device includes: The data correction module is used to perform dynamic correction processing on the first common center point gather data to obtain the second common center point gather data; The Radon spectrum creation module is used to perform a parabolic Radon domain record by superimposing the traces along the parabolic trajectory when the second common center point gather data traverses different zero offset distances and different parabolic curvatures, and to perform in-phase weighting on the parabolic Radon domain record and take the absolute value to obtain the parabolic Radon spectrum. The phase axis fitting module is used to perform parabolic phase axis tracking on the parabolic Radon spectrum to obtain the time position of each multiple phase axis in each seismic trace in the second common center point gather data; The in-phase axis extraction module is used to extract the first remaining multiple in-phase axis from each multiple in-phase axis using KL transform filtering; The data reaction correction module is used to perform reaction correction processing on the first residual multiple phase axis to obtain the second residual multiple phase axis. The multiple wave cancellation module is used to eliminate the remaining multiple waves by least square filtering based on the first common center point gather data and the second remaining multiple wave phase axis, so as to obtain the third common center point gather data.

8. An electronic device, characterized in that, The electronic device includes: At least one processor; and A memory communicatively connected to the at least one processor; wherein, The memory stores a computer program that can be executed by the at least one processor to enable the at least one processor to perform the method of any one of claims 1-6.

9. A computer-readable storage medium, characterized in that, The computer-readable storage medium stores computer instructions that cause a processor to execute the method of any one of claims 1-6.

10. A computer program product, characterized in that, The computer program product includes a computer program that, when executed by a processor, implements the method according to any one of claims 1-6.

Images