A method for decoupling and separating multi-source aliasing seismic data

Abstract

The application relates to the technical field of oil exploration, and particularly discloses a multi-seismic-source aliasing seismic data decoupling separation method, which comprises the following steps: inputting aliasing seismic data and establishing a mathematical model, converting a separation problem into a signal recovery problem with constraints; reorganizing the data into a common offset domain, constructing a cascade shift embedded Hankel matrix to form a low-rank prior; designing an energy-aware modulation function based on spectral component values to adaptively distinguish effective signals and crosstalk; fusing a weighted spectral norm and a transform domain sparsity constraint to construct a joint structured penalty optimization model; and adopting a variable splitting iteration framework to decompose a sub-problem, alternately solving and adaptively updating weights. The application protects effective signals by relying on an adaptive weight strategy, separates crosstalk by combining low-rank and sparse double characteristics, overcomes signal leakage and amplitude distortion problems, improves waveform fidelity and separation signal-to-noise ratio, and takes into account convergence stability and calculation efficiency, so that the application has important application value in complex geological oil and gas exploration and high-precision seismic imaging.

Description

Technical Field

[0001] This invention relates to the field of petroleum exploration technology, specifically to a method for decoupling and separating multi-source superimposed seismic data. Background Technology

[0002] In oil and gas exploration, multi-source cascade acquisition is a highly innovative seismic data acquisition technology that can break through the time limitations of traditional independent excitation, significantly improve acquisition efficiency and reduce operating costs, and obtain massive amounts of high-density, wide-azimuth seismic data. However, this technology directly leads to the overlap of wavefields excited by different sources in the record, forming severe random or coherent crosstalk between sources, which seriously masks the effective reflection signal, causing a sharp drop in the data signal-to-noise ratio and affecting the accuracy of subsequent imaging.

[0003] Therefore, for the needs of oil and gas exploration, it is particularly important to develop a method that can efficiently extract independent source wavefields from aliased records. It is urgent to study a multi-source seismic data decoupling and reconstruction method that can effectively achieve high-fidelity data separation in complex interference environments and avoid effective wave amplitude damage and signal leakage during the denoising process. Summary of the Invention

[0004] The purpose of this invention is to solve the problem of overlapping seismic records acquired by multi-source mixed excitation, resulting in aliased data, and to provide a decoupling and separation method for multi-source aliased seismic data based on adaptive low-rank compact characterization constraints.

[0005] To achieve the above objectives, the technical solution provided by this invention is: a method for decoupling and separating multi-source cascading seismic data, comprising the following steps: S1. Input the seismic data acquired by multi-source mixed excitation, establish a mathematical model of multi-source mixed seismic data, and transform the aliasing separation problem into a constrained signal recovery problem. S2. Reassemble the aliased data from the common shot domain to the common offset domain, construct a cascaded translational embedded Hankel matrix for each common offset gather along the spatial direction, and form a low-rank prior structure by utilizing the spatial coherence of seismic data. S3. Design an energy-sensing modulation function based on the spectral component values ​​of the Hankel matrix, assign differentiated weights to different spectral components, and achieve adaptive differentiation between effective signals and crosstalk interference. S4. By integrating the weighted spectral norm and constraints with the sparsity of the transform domain, a low-rank-sparse joint structured penalty optimization model is constructed. S5. The original optimization problem is decomposed into multiple sub-problems using a variable splitting iterative framework, and the sub-problems are solved alternately and the weights are updated adaptively. Furthermore, in step S1, the mathematical model for the multi-source superimposed seismic data is established as follows: In simultaneous excitation and acquisition by multiple seismic sources, assuming a total of Multiple seismic sources participate in the excitation, with random time delays set between each source to avoid completely synchronous excitation; let the first... The record of conventional earthquakes generated by a single seismic source is as follows: This includes the number of time sampling points. and number of receiving channels Two dimensions of information; defining a corresponding wavefield coupling mapping operator for each seismic source. The wavefield coupling mapping operator is used to comprehensively describe the mapping relationship from single source data to aliased data, and specifically includes the following two basic operations: The first operation is a timing shift operation. According to the Excitation delay of each seismic source The seismic data is subjected to a cyclic shift along the time axis, expressed by the following formula: ;in Indicates the time sampling point. Indicates the channel number to be received.

[0006] The second operation is the spatial extraction operation. Based on the acquired geometric relationships, seismic data undergoes trace extraction or spatial interpolation processing, using the following formula: ;in For the first Trace mapping function corresponding to each earthquake source; Combining the above two operations, the first The wavefield coupling mapping operator corresponding to each earthquake source is defined as: That is, a temporal translation operation is performed first, followed by a spatial extraction operation; based on the above definition, the received aliased seismic data This can be represented as a linear superposition of all source data after being processed by their respective wavefield coupling mapping operators: ;in This is a stacked vector of all source data. The total wavefield coupling mapping operator matrix, This refers to random inter-source crosstalk clutter during the acquisition process; the goal of aliased data separation is to separate the known aliased data... Wavefield coupling mapping operator Under these conditions, the original seismic records corresponding to each seismic source are reconstructed from the observation data. .

[0007] Furthermore, step S2 is specifically implemented as follows: Assuming aliased data The Middle Cannon The data of Tao is Define offset distance ,in The spatial coordinates of the gun point, Spatial coordinates of the receiving point; common offset domain data By extracting all offsets To obtain through the Way: ;in, The trace number in the common offset trace set is used, and the center point coordinates are taken. As a spatial marker of the Dao; For each common offset gather, construct a concatenated translation embedding matrix along the spatial direction, and let the common offset gather... Include Channel data, each channel data is denoted as Each data point has a length of [length missing]. Column vectors; Cascaded translation embedding matrix The construction method is as follows: adjacent The data is stacked vertically to form one column of a matrix, and then a sliding window is used to construct each column of the matrix sequentially; specifically, the matrix... The The column is from the first Dao Dao Di The data is stacked vertically, and the matrix contains a total of The column, the formula is: ; Equivalently, the above construction process can be represented in operator form: ;in Construct operators for the Hankel matrix. The number of rows in a block controls the size of the Hankel matrix and the strength of its low-rank properties; The value of needs to be balanced between low rank and computational efficiency. That is, half of the number of channels, rounded down.

[0008] Furthermore, in step S3, the specific design of the energy-sensing modulation function is as follows: The Hankel matrix constructed in step S2 Spectral decomposition yields: ;in It is a left singular vector matrix. It is a right singular vector matrix. This is a diagonal matrix of spectral component values. The spectral component values ​​are arranged in descending order. Let be the rank of the matrix, and be the superscript. Indicates conjugate transpose; Design an energy-sensing modulation function, assigning different suppression weights to different spectral component values: ;in For the first The weights corresponding to each spectral component value For the first Spectral component values This is a numerical robustness parameter used to prevent numerical instability caused by the denominator approaching zero. ,in , This represents the maximum spectral component value. The weights are adaptively updated during the iteration process based on the current estimation results. Let the weights be the first... The estimated signal at the next iteration is The corresponding Hankel matrix is The weight update formula is: ;in Representation matrix The Each spectral component value; initial weights are taken as follows: As the iteration proceeds, the weights gradually adapt to the actual spectral component value distribution of the data, achieving adaptive adjustment.

[0009] Furthermore, step S4 is specifically as follows: The definition of the weighted spectral norm sum is based on the energy-sensing modulation function, and is applied to the matrix. Define its weighted spectral norm sum as: ; in, For matrix The Spectral component values The corresponding adaptive weights designed in step S3; Orthogonal frequency basis expansion transform is used as a sparse transform operator. Two-dimensional orthogonal frequency basis expansion transforms spatiotemporal signals to the frequency domain; The sparsity constraint uses the ℓ1 norm: The ℓ1 norm is a convex relaxation of the ℓ0 norm. By combining low-rank constraints and sparsity constraints, an optimization model with a weighted spectral norm and a joint structured penalty for sparsity is established: ; The For the weighted spectral norm and terms, where, For the transformation operator to the common offset domain, Construct operators for the Hankel matrix; For sparse regularization, It is an orthogonal frequency basis expansion transform operator; For data fidelity items, It is the Frobenius norm; This is a data fidelity parameter, with a value range of [0.5, 2].

[0010] Furthermore, step S5 is specifically as follows: Perform variable splitting and introduce auxiliary variables. and ,make The Hankel matrix representation of the signal in the common offset domain. This represents the sparse representation of the signal in the transform domain; it transforms the original unconstrained optimization problem into an equivalent constrained optimization problem. ; The constraints are: , ; Introducing Lagrange multipliers and An augmented Lagrangian function is constructed using a proportional approach, and then minimized by alternately updating the variables; let the current iteration number be... The suppression parameters are The update steps for each variable are as follows: Sub-problem update: Fix other variables, The update subproblem is a least-squares problem involving the weighted spectral norm and structured penalty; this subproblem has a closed-form solution, which is obtained through a flexible threshold filter based on the weighted spectral component values; the update formula is: ; in , From the matrix The spectral decomposition process performs threshold shrinkage on each spectral component value, with the threshold value being proportional to the corresponding weight. Sub-problem update: Fix other variables, The update subproblem is a least squares problem that solves the ℓ1 norm structured penalty; this subproblem is also solved using flexible threshold filtering. The updated formula is: ; This allows for independent soft thresholding of each element of the matrix, where... ; This is a sparse structured penalty parameter, with a value range of [0.01, 0.1]. Sub-problem update: Fix other variables, The update subproblem is a quadratic optimization problem, and its optimality condition corresponds to a linear system: ; in The adjoint operator for constructing the Hankel matrix corresponds physically to the Hankel averaging operation; the linear system is solved efficiently using the conjugate gradient method; when the wave field is coupled to the mapping operator... Dao Collection Recombination Operator Sparse transformation operators When all operators are orthogonal or quasi-orthogonal, the coefficient matrix has a good condition number, and the conjugate gradient method converges quickly. The linear system can be further simplified to: ; in Construct the adjoint operator for the Hankel matrix; Lagrange multiplier update: Update the Lagrange multipliers according to the degree of constraint violation. The update formula is as follows: ; ; Weight Update: The adaptive weights are updated based on the distribution of spectral component values ​​in the current iteration result. The update formula is as follows: ; The weights are gradually adapted to the actual structural characteristics of the data. During algorithm initialization, the estimated signal Initialization is performed using the adjoint operation of the wave-field coupling mapping operator, i.e. The Lagrange multipliers are initialized to zero, and the weights are initialized to 1; the suppression parameters... The value is 1.0.

[0011] The advantages of this invention compared to existing technologies are as follows: Based on adaptive weighted spectral norm and minimization theory combined with multi-domain constraints, this invention achieves accurate identification and protection of large spectral component values ​​representing effective signals by introducing an adaptive weighting strategy inversely proportional to the spectral component values. Furthermore, it leverages the dual characteristics of low rank in the co-offset domain and sparsity in the transform domain of seismic data to achieve efficient separation of crosstalk clutter from complex aliasing sources. This invention effectively overcomes the "signal leakage" and amplitude distortion problems caused by excessive suppression of effective components in traditional denoising processes. It not only significantly improves waveform fidelity and separation signal-to-noise ratio but also ensures the convergence stability and computational efficiency of the algorithm when processing large-scale data by relying on the ADMM solution framework. Therefore, this method for decoupling and reconstructing aliased data based on adaptive weights and dual prior constraints has significant research significance and application value in oil and gas exploration and high-precision seismic imaging processing under complex geological conditions. Attached Figure Description

[0012] Figure 1This is a flowchart of a method for decoupling and separating multi-source superimposed seismic data according to the present invention.

[0013] Figure 2 This is a raw single-shot record diagram without coupling crosstalk, based on the Marmousi model forward model.

[0014] Figure 3 This is a composite seismic record image based on the Marmousi model forward model and after random time delay stacking.

[0015] Figure 4 This is the result image after separation using the traditional method.

[0016] Figure 5 This is a diagram showing the result of separation using the method of the present invention.

[0017] Figure 6 It is a composite seismic record map obtained from multiple sources in the actual work area.

[0018] Figure 7 It is a seismic data map of the actual work area obtained by using the traditional decoupling and reconstruction method.

[0019] Figure 8 This is a seismic data map of the actual work area obtained using the method of this invention. Detailed Implementation

[0020] Various exemplary embodiments of the present invention will now be described in detail with reference to the accompanying drawings. It should be noted that, unless otherwise specifically stated, the relative arrangement, numerical expressions, and values ​​of the components and steps set forth in these embodiments do not limit the scope of the present invention.

[0021] The following description of at least one exemplary embodiment is merely illustrative and is in no way intended to limit the invention or its application or use.

[0022] Techniques, methods, and equipment known to those skilled in the art may not be discussed in detail, but where appropriate, such techniques, methods, and equipment should be considered part of the specification.

[0023] In all examples shown and discussed herein, any specific values ​​should be interpreted as merely exemplary and not as limitations. Therefore, other examples of exemplary embodiments may have different values.

[0024] The following detailed description of a multi-source superimposed seismic data decoupling and separation method of the present invention, with reference to the accompanying drawings, provides further insight.

[0025] Combined with appendix Figures 1-8 This invention will be described in detail below.

[0026] A method for decoupling and separating multi-source cascading seismic data includes the following steps: Step S1: Establish a mathematical model for multi-source superimposed seismic data or input seismic data acquired through multi-source mixed excitation; The method for decoupling and separating multi-source cascading seismic data based on adaptive low-rank compact representation constraints according to claim 1 is characterized in that the process of establishing the mathematical model of the multi-source cascading seismic data in step one is as follows: In simultaneous excitation and acquisition by multiple seismic sources, assuming a total of Multiple seismic sources participate in the excitation, with random time delays set between each source to avoid completely synchronous excitation; let the first... The record of conventional earthquakes generated by a single seismic source is as follows: This includes the number of time sampling points. and number of receiving channels Two dimensions of information; defining a corresponding wavefield coupling mapping operator for each seismic source. The wavefield coupling mapping operator is used to comprehensively describe the mapping relationship from single source data to aliased data, and specifically includes the following two basic operations: The first operation is a timing shift operation. Its physical meaning is based on the first Excitation delay of each seismic source The mathematical expression for performing a cyclic shift on the time axis on seismic data is: ;in Indicates the time sampling point. This indicates the received channel number. This operation reflects the time migration effect caused by the different excitation times of different seismic sources.

[0027] The second operation is the spatial extraction operation. Its physical meaning is to perform trace extraction or spatial interpolation processing on seismic data according to the acquisition geometry. The mathematical expression is: ;in For the first The trace mapping function corresponding to each seismic source describes the spatial sampling pattern at the time of excitation by that source. This operation reflects the geometric differences in acquisition caused by the different spatial locations of different seismic sources.

[0028] Combining the two operations above, the first The wavefield coupling mapping operator corresponding to each earthquake source is defined as: That is, a temporal translation operation is performed first, followed by a spatial extraction operation. Based on the above definition, the received mixed seismic data... This can be represented as a linear superposition of all source data after being processed by their respective wavefield coupling mapping operators: ; in This is a stacked vector of all source data. The total wavefield coupling mapping operator matrix, This refers to random inter-source crosstalk clutter during the acquisition process. The goal of aliased data separation is to separate the known aliased data... Wavefield coupling mapping operator Under these conditions, the original seismic records corresponding to each seismic source are reconstructed from the observation data. Due to the wave-field coupling mapping operator The problem is usually underdetermined, and it has infinitely many solutions. Therefore, it is necessary to introduce additional prior information to impose structured penalty constraints on the solution space, which is the starting point for the design of subsequent steps.

[0029] Step S2: Collection reorganization and cascaded translation embedding matrix construction; The purpose of gather reassembly is to reassemble aliased data from common shot gathers into the common offset domain, thereby fully utilizing the spatial coherence of seismic data. In the common shot domain, the phase axes of seismic waves from different reflection interfaces within the same shot gather exhibit a hyperbolic distribution, indicating weak spatial coherence. In contrast, in the common offset domain, the phase axes of seismic waves from the same reflection interface appear as approximately horizontal or gradually varying curves, demonstrating good spatial coherence. This characteristic provides a physical basis for establishing intrinsic low-dimensional manifold constraints.

[0030] Assuming aliased data The Middle Cannon The data of Tao is Define offset distance ,in The spatial coordinates of the gun point, Spatial coordinates of the receiving point; common offset domain data By extracting all offsets To obtain through the Way: ;in, The trace number in the common offset trace set is used, and the center point coordinates are taken. As a spatial location marker of the Dao; through defining the Dao set recombination operator To achieve the above transformation process.

[0031] In the common offset domain, to further enhance the low-rank structure of the data and facilitate subsequent matrix factorization, a concatenated translation embedding matrix is ​​constructed along the spatial direction for each common offset gather. Let the common offset gather... Include Channel data, each channel data is denoted as Each data point has a length of [length missing]. Column vectors.

[0032] Cascaded translation embedding matrix The construction method is as follows: adjacent The data is stacked vertically to form one column of a matrix, and then a sliding window is used to construct each column of the matrix sequentially; specifically, the matrix... The The column is from the first Dao Dao Di The data is stacked vertically, and the matrix contains a total of The column, mathematical expression is: ; Equivalently, the above construction process can be represented in operator form: ;in Construct operators for the Hankel matrix. The number of rows in a block controls the size of the Hankel matrix and the strength of its low-rank properties; The value of needs to be balanced between low rank and computational efficiency. That is, half of the number of channels, rounded down.

[0033] Based on the physical characteristics of seismic data, when seismic phase axes exhibit linear or slowly varying spatial features, the constructed Hankel matrix possesses low-rank properties. Intuitively, if seismic data consists of K ideal linear phase axes, the theoretical rank of the Hankel matrix is ​​K. In reality, due to the slowly varying nature of phase axes, the finite bandwidth effect, and the presence of inter-source crosstalk clutter, the Hankel matrix exhibits approximately low rank. That is, the matrix's energy is mainly concentrated in a few large spectral components, while the smaller spectral components primarily correspond to inter-source crosstalk clutter and interference.

[0034] When coupling crosstalk exists, the spatial coherence of the interfering signal in the current common-offset channel set is disrupted because the interfering signal originates from seismic sources at different spatial locations, exhibiting characteristics similar to random inter-source crosstalk clutter. Therefore, the coupling crosstalk component is mainly distributed on the smaller spectral components of the Hankel matrix, while the effective signal component is mainly distributed on the larger spectral components. This distribution characteristic provides a theoretical basis for subsequent adaptive weighting processing based on spectral components.

[0035] Step S3: Introduce an energy-sensing modulation function for weighting; First, perform spectral decomposition on the Hankel matrix H constructed in step S2 to obtain: ;in It is a left singular vector matrix. It is a right singular vector matrix. This is a diagonal matrix of spectral component values. The spectral component values ​​are arranged in descending order. Let be the rank of the matrix, and be the superscript. This indicates the conjugate transpose.

[0036] The magnitude of the spectral component values ​​directly reflects the energy intensity of the signal component represented by the corresponding singular vector. Based on the analysis in step S2, for aliased seismic data: components with large spectral component values ​​correspond to effective signals with strong spatial coherence. These signals exhibit a low-rank structure in the Hankel matrix, with energy concentrated in the first few spectral component values; components with small spectral component values ​​correspond to coupling crosstalk and random inter-source crosstalk clutter whose spatial coherence is destroyed. Although coupling crosstalk has coherence in its original domain, its coherence is lost after transformation to the current co-offset domain, exhibiting characteristics similar to inter-source crosstalk clutter.

[0037] Based on the above physical analysis, this invention designs an energy-sensing modulation function that assigns different suppression weights to different spectral component values: ;in For the first The weights corresponding to each spectral component value For the first Spectral component values This is a numerical robustness parameter used to prevent numerical instability caused by the denominator approaching zero. ,in , This represents the maximum spectral component value.

[0038] The core design idea of ​​this weighting function is reflected in the following two aspects: First, when the spectral component values When it is large, the corresponding weight The smaller the spectral component value, the less suppression it experiences in subsequent optimizations, thus preserving the effective signal components and avoiding excessive compression of the effective signal; secondly, when the spectral component value... When it is small, the corresponding weight The magnitude is relatively large, and the suppression of this component is relatively large in subsequent optimization, thereby effectively suppressing coupling crosstalk and inter-source crosstalk clutter components.

[0039] Compared to traditional standard spectral norm summing methods, which assign the same weight to all spectral components, standard spectral norm summing fails to distinguish between valid signals and interference components, easily leading to over-compression of large spectral components and causing amplitude loss and waveform distortion of the valid signal. The adaptive weighting mechanism of this invention achieves differentiated processing of signals and interference, effectively protecting the amplitude and waveform characteristics of the valid signal while suppressing coupling crosstalk.

[0040] Since the distribution of spectral component values ​​of the separated data cannot be accurately obtained in the initial stage of the algorithm, the weights need to be adaptively updated based on the current estimation results during the iteration process. Let the... The estimated signal at the next iteration is The corresponding Hankel matrix is The weight update formula is: ;in Representation matrix The Each spectral component value; initial weights are taken as follows: As the iteration proceeds, the weights gradually adapt to the actual spectral component value distribution of the data, achieving adaptive adjustment.

[0041] Step S4: Construct a weighted spectral norm and sparsity joint structured penalty optimization model; The definition of the weighted spectral norm sum is based on the energy-sensing modulation function, and is applied to the matrix. Define its weighted spectral norm sum as: ;in, For matrix The Spectral component values The corresponding adaptive weights designed in step S3. And the standard spectral norm and In comparison, the weighted spectral norm and the degree of suppression of each spectral component value are different, which can better balance signal preservation and interference suppression.

[0042] The weighted spectral norm sum has the following important properties: First, since the weights depend on the spectral component values ​​themselves, the weighted spectral norm sum is a non-convex function, which is different from the convexity of the standard spectral norm sum; second, the weighted spectral norm sum is a tighter approximation of the rank function, and can more accurately characterize the low-rank properties of the matrix compared to the standard spectral norm sum; third, the weights are automatically adjusted with iteration, without the need for manual preset, and have good adaptability.

[0043] Besides exhibiting low-rank properties in the spatial domain, seismic data also possesses sparse representation characteristics in certain transform domains. This invention employs orthogonal frequency basis expansion transform as a sparse transform operator. Two-dimensional orthogonal frequency basis expansion transforms a spatiotemporal signal into the frequency domain. For a seismic signal composed of a finite number of frequency components, its transformation coefficients are sparsity, meaning that most coefficients are close to zero and only a few coefficients have large amplitude values.

[0044] The sparsity constraint uses the ℓ1 norm: The ℓ1 norm is a convex relaxation of the ℓ0 norm (the number of non-zero elements), which can effectively promote the recovery of sparse solutions and has computational advantages.

[0045] By combining low-rank constraints and sparsity constraints, an optimization model with weighted spectral norm and joint structured penalty for sparsity is established: ; The optimization model consists of three components: the first term For the weighted spectral norm and terms, where For the transformation operator to the common offset domain, The operator for constructing the Hankel matrix is ​​introduced, which promotes the low-rank nature of the solution in the Hankel matrix representation in the common offset domain and protects the effective signal components through adaptive weighting; the second term For sparse regularization, For orthogonal frequency basis expansion transform operators, this term promotes the sparsity of the solution in the transform domain; the third term For data fidelity items, It is the Frobenius norm, which guarantees the consistency between the solution and the observed aliased data.

[0046] In the model The sparse structure penalty parameter controls the strength of the sparsity constraint. Too large a value will result in over-smoothing. If the value is too small, the compact characterization constraint will be insufficient; the recommended value range is [0.01, 0.1]. This parameter is for data fidelity and controls the strength of data fitting. It is relatively insensitive to the results and the recommended value range is [0.5, 2].

[0047] The physical meaning of this joint structured penalty model is: under the premise of ensuring consistency with the observation data, it seeks a solution that simultaneously satisfies the dual structural characteristics of low rank and sparsity, and the low rank constraint adaptively protects the effective signal components through a weighting mechanism, fully mining the inherent structural information of the seismic data, thereby achieving high-precision separation of aliased data.

[0048] Step S5: Solve iteratively using a variable splitting iterative framework; Since the objective function of the optimization model established in step S4 contains multiple non-smooth terms, it is difficult to solve directly. This invention adopts a variable splitting iterative framework, which decomposes the original problem into several subproblems by introducing auxiliary variables, and then iteratively solves each subproblem until convergence.

[0049] First, perform variable splitting and introduce auxiliary variables. and ,make The Hankel matrix representation of the signal in the common offset domain. This represents the sparse representation of the signal in the transform domain. The original unconstrained optimization problem is transformed into an equivalent constrained optimization problem: ; The constraints are: , ; Introducing Lagrange multipliers and An augmented Lagrangian function is constructed using a proportional approach, and then minimized by alternately updating the variables. Let the current iteration number be . The suppression parameters are The update steps for each variable are as follows: Sub-problem update: Fix other variables, The update subproblem is a least-squares problem involving the weighted spectral norm and structured penalty. This subproblem has a closed-form solution, which is obtained through a flexible threshold filter based on the weighted spectral component values. The update formula is: ; in , From the matrix The spectral decomposition performs threshold shrinkage on each spectral component value, with the threshold value being proportional to the corresponding weight.

[0050] Sub-problem update: Fix other variables, The update subproblem is a least-squares problem involving the ℓ1 norm structured penalty. This subproblem is also solved using flexible threshold filtering.

[0051] The updated formula is: ; This allows for independent soft thresholding of each element in the matrix, where... .

[0052] Sub-problem update: Fix other variables, The update subproblem is a quadratic optimization problem, and its optimality condition corresponds to a linear system: ; in The adjoint operator for constructing the Hankel matrix corresponds physically to the Hankel averaging operation. This linear system can be solved efficiently using the conjugate gradient method. When the wave-field coupling mapping operator... Dao Collection Recombination Operator Sparse transformation operators When all operators are orthogonal or quasi-orthogonal, the coefficient matrix has a good condition number, and the conjugate gradient method converges quickly.

[0053] The linear system can be further simplified to: ;in Construct the adjoint operator for the Hankel matrix; Lagrange multiplier update: Update the Lagrange multipliers according to the degree of constraint violation. The update formula is as follows: ; ; Weight Update: The adaptive weights are updated based on the distribution of spectral component values ​​in the current iteration result. The update formula is as follows: ; The weights are gradually adapted to the actual structural characteristics of the data.

[0054] During algorithm initialization, the estimated signal Initialization is performed using the adjoint operation of the wave-field coupling mapping operator, i.e. The Lagrange multipliers are initialized to zero, and the weights are initialized to 1. Suppression parameters. The typical value is 1.0.

[0055] Step S6: Iterative convergence judgment.

[0056] Step S7: Output the separation results.

[0057] The specific implementation process of the multi-source superimposed seismic data decoupling and separation method of the present invention is as follows: 1. Marmousi model testing To comprehensively verify the effectiveness and amplitude preservation performance of the adaptive weighted joint constraint method proposed in this invention when processing seismic data with complex geological structures, this embodiment selects the Marmousi velocity model, which has strong lateral velocity variation and complex fault block characteristics, as the test benchmark. The excitation source is a zero-phase Ricker wavelet with a dominant frequency of 25Hz, placed at a depth of 2m. The first shot is located 50m horizontally, with a shot spacing of 25m, and a total of 320 shots are fired. The receiver arrangement adopts a fixed receiving mode with full surface coverage, with 320 geophones deployed, also with a channel spacing of 25m. The time sampling interval is set to 2ms, and the number of sampling points recorded per shot is 2000, i.e., the total recording time is 4.0s. Based on this, in order to realistically simulate the aliasing acquisition operation scenario, an aliasing data model was constructed by simultaneously exciting two seismic sources from a ship. The spatial interval between the two seismic sources was set to 50m, and the excitation delay time of the second seismic source relative to the first seismic source strictly followed a uniform random distribution in the interval [100ms, 500ms]. Through this random delay superposition, an aliasing seismic record with irregular and strong interference characteristics was generated.

[0058] The data separation and processing procedure strictly follows Figure 1 The technical approach shown involves using an adaptive weight matrix inversely proportional to the spectral component values ​​to iteratively update and separate the aliased wave field. The algorithm is set such that when the relative error between two adjacent iterations is less than [a certain value], the aliased wave field is separated. The calculation terminates when the number of iterations reaches a preset maximum value. To visually demonstrate the superiority of the method of this invention, an experiment was conducted comparing the method of this invention with traditional standard spectral norm and minimization methods. The test results are as follows: Figures 2 to 5 As shown. Specifically, Figure 2The original single-shot record without coupling crosstalk is shown, with clear wavefield and well-defined stratigraphic levels, serving as the standard reference for this test; Figure 3 The image shows the superimposed seismic records after random time delay. It can be seen that the wave fields from different sources are severely intertwined and overlapped in the spatiotemporal domain. The effective reflection signals from deep layers are almost completely masked by the strong-energy crosstalk clutter between sources. Figure 4 The results show the separation achieved using a traditional method. While this method suppresses background inter-source crosstalk to some extent, significant artifacts remain in regions of intense wavefield crossover, and the amplitude of some effective signals is excessively smoothed. In contrast, Figure 5 The results of separation using the method of this invention are presented. Thanks to the introduction of dual prior constraints and an adaptive weighting strategy, this result not only completely eliminates random aliasing inter-source crosstalk clutter, but also perfectly restores the continuity and detail features of the in-phase axis, with extremely high waveform fidelity. Figure 2 The actual records shown are highly consistent, which fully demonstrates the robustness and advancement of the method of the present invention in complex structures and strong aliasing backgrounds.

[0059] 2. Actual work area data testing To verify the effectiveness and adaptability of the method of the present invention in actual industrial production, actual seismic data collected in the exploration area were selected for testing. Figure 6 The study presents the aliasing record of actual data, revealing severe wavefield aliasing in the seismic record due to the combined excitation from multiple sources. The phase axes from different sources intersect and overlap in time and space, making it extremely difficult to directly identify effective reflection horizons. This poses a significant challenge to subsequent velocity analysis and imaging processing. Figure 7 For the separated seismic data obtained by the traditional decoupling reconstruction method, although the coupling crosstalk is suppressed to a certain extent after processing by the traditional decoupling reconstruction method, residual inter-source crosstalk clutter still exists in some areas, and the continuity of the effective signal is somewhat lacking. Figure 8 To separate the shot record using the method of this invention, a comparison is made. Figure 7 and Figure 8As can be seen, after applying the decoupling and reconstruction method based on adaptive weights and dual constraints proposed in this invention, the wavefields of the mixed seismic sources were successfully and clearly separated. The comparison shows that the separated recording background is clean, the phase axis continuity is good, and the discontinuities and fine structural morphology are clearly visible, with no obvious residual crosstalk, inter-source crosstalk clutter, or artificial artifacts. The above test results demonstrate that the method of this invention has significant advantages in processing real-world complex seismic data. The algorithm of this invention can intelligently identify data characteristics, applying smaller suppression weights to the large spectral components representing the effective signal, while applying larger suppression weights to the small spectral components representing inter-source crosstalk clutter. This mechanism efficiently suppresses coupling crosstalk while preserving the amplitude and phase information of the actual seismic data to the greatest extent, avoiding the "signal leakage" and "signal distortion" problems common in traditional methods. It has extremely high waveform fidelity and can meet the actual needs of high-precision oil and gas exploration.

[0060] The present invention and its embodiments have been described above. This description is not restrictive, and the accompanying drawings are only one embodiment of the present invention; the actual structure is not limited thereto. In conclusion, if those skilled in the art are inspired by this description and design similar structures and embodiments without departing from the spirit of the invention, such designs should fall within the protection scope of the present invention.

Claims

1. A method for decoupling and separating multi-source superimposed seismic data, characterized in that, Includes the following steps: S1. Input the seismic data acquired by multi-source mixed excitation, establish a mathematical model of multi-source mixed seismic data, and transform the aliasing separation problem into a constrained signal recovery problem. S2. Reassemble the aliased data from the common shot domain to the common offset domain, construct a cascaded translational embedding Hankel matrix for each common offset gather along the spatial direction, and form a low-rank prior structure by utilizing the spatial coherence of seismic data. S3. Design an energy-sensing modulation function based on the spectral component values ​​of the Hankel matrix, assign differentiated weights to different spectral components, and achieve adaptive differentiation between effective signals and crosstalk interference; S4. By integrating the weighted spectral norm and constraints with the sparsity of the transform domain, a low-rank-sparse joint structured penalty optimization model is constructed. S5. The original optimization problem is decomposed into multiple subproblems using a variable splitting iterative framework, and the subproblems are solved iteratively and the weights are updated adaptively.

2. The method for decoupling and separating multi-source superimposed seismic data according to claim 1, characterized in that: In step S1, the mathematical model for the multi-source superimposed seismic data is established as follows: In simultaneous excitation and acquisition by multiple seismic sources, assuming a total of Multiple seismic sources participate in the excitation, with random time delays set between each source to avoid completely synchronous excitation; let the first... The record of conventional earthquakes generated by a single seismic source is as follows: This includes the number of time sampling points. and number of receiving channels Two dimensions of information; defining a corresponding wavefield coupling mapping operator for each seismic source. The wavefield coupling mapping operator is used to comprehensively describe the mapping relationship from single source data to aliased data, and specifically includes the following two basic operations: The first operation is a timing shift operation. According to the Excitation delay of each seismic source The seismic data is subjected to a cyclic shift along the time axis, expressed by the following formula: ;in Indicates the time sampling point. Indicates the channel number to be received; The second operation is spatial extraction. Based on the acquired geometric relationships, seismic data undergoes trace extraction or spatial interpolation processing, using the following formula: ;in For the first Trace mapping function corresponding to each earthquake source; Combining the above two operations, the first The wavefield coupling mapping operator corresponding to each earthquake source is defined as: That is, a temporal translation operation is performed first, followed by a spatial extraction operation; based on the above definition, the received aliased seismic data This can be represented as a linear superposition of all source data after being processed by their respective wavefield coupling mapping operators: ;in This is a stacked vector of all source data. The total wavefield coupling mapping operator matrix, This refers to random inter-source crosstalk clutter during the acquisition process; the goal of aliased data separation is to separate the known aliased data... Coupled with wave field mapping operator Under these conditions, the original seismic records corresponding to each seismic source are reconstructed from the observation data. .

3. The method for decoupling and separating multi-source superimposed seismic data according to claim 2, characterized in that: The specific implementation method of step S2 is as follows: Assuming aliased data The Middle Cannon The data of Tao is Define offset distance ,in The spatial coordinates of the gun point, Spatial coordinates of the receiving point; common offset domain data By extracting all offsets To obtain through the Way: ;in, The trace number in the common offset trace set is used, and the center point coordinates are taken. As a spatial marker of the Dao; For each common offset gather, construct a concatenated translation embedding matrix along the spatial direction, and let the common offset gather... Include Channel data, each channel data is denoted as Each data point has a length of [length missing]. Column vectors; Cascaded translation embedding matrix The construction method is as follows: adjacent The data is stacked vertically to form one column of a matrix, and then a sliding window is used to construct each column of the matrix sequentially; specifically, the matrix... The The column is from the first Dao Dao Di The data is stacked vertically, and the matrix contains a total of The column, the formula is: ; Equivalently, the above construction process can be represented in operator form: ;in Construct operators for the Hankel matrix. The number of rows in a block controls the size of the Hankel matrix and the strength of its low-rank properties; The value of needs to be balanced between low rank and computational efficiency. That is, half of the number of channels, rounded down.

4. The method for decoupling and separating multi-source superimposed seismic data according to claim 3, characterized in that: In step S3, the specific design method of the energy-sensing modulation function is as follows: The Hankel matrix constructed in step S2 Spectral decomposition yields: ;in It is a left singular vector matrix. It is a right singular vector matrix. This is a diagonal matrix of spectral component values. The spectral component values ​​are arranged in descending order. Let be the rank of the matrix, and be the superscript. Indicates conjugate transpose; Design an energy-sensing modulation function, assigning different suppression weights to different spectral component values: ;in For the first The weights corresponding to each spectral component value For the first Spectral component values This is a numerical robustness parameter used to prevent numerical instability caused by the denominator approaching zero. ,in , This represents the maximum spectral component value. The weights are adaptively updated during the iteration process based on the current estimation results. Let the weights be the first... The estimated signal at the next iteration is The corresponding Hankel matrix is The weight update formula is: ;in Representation matrix The Each spectral component value; initial weights are taken as follows: As the iteration proceeds, the weights gradually adapt to the actual spectral component value distribution of the data, achieving adaptive adjustment.

5. The method for decoupling and separating multi-source superimposed seismic data according to claim 4, characterized in that: Step S4 is as follows: The definition of the weighted spectral norm sum is based on the energy-sensing modulation function, and is applied to the matrix. Define its weighted spectral norm sum as: ； in, For matrix The Spectral component values The corresponding adaptive weights designed in step S3; Orthogonal frequency basis expansion transform is used as a sparse transform operator. Two-dimensional orthogonal frequency basis expansion transforms spatiotemporal signals to the frequency domain; The sparsity constraint uses the ℓ1 norm: The ℓ1 norm is a convex relaxation of the ℓ0 norm. By combining low-rank constraints and sparsity constraints, an optimization model with a weighted spectral norm and a joint structured penalty for sparsity is established: ； The For the weighted spectral norm and terms, where, For the transformation operator to the common offset domain, Construct operators for the Hankel matrix; For sparse regularization, It is an orthogonal frequency basis expansion transform operator; For data fidelity items, It is the Frobenius norm; This is a data fidelity parameter, with a value range of [0.5, 2].

6. The method for decoupling and separating multi-source superimposed seismic data according to claim 5, characterized in that: Step S5 is as follows: Perform variable splitting and introduce auxiliary variables. and ,make The Hankel matrix representation of the signal in the common offset domain. This represents the sparse representation of a signal in the transform domain. The original unconstrained optimization problem is transformed into an equivalent constrained optimization problem: ； The constraints are: , ; Introducing Lagrange multipliers and An augmented Lagrangian function is constructed using a proportional approach, and then minimized by alternately updating the variables; let the current iteration number be... The suppression parameters are The update steps for each variable are as follows: Sub-problem update: Fix other variables, The update subproblem is a least-squares problem involving the weighted spectral norm and structured penalty; this subproblem has a closed-form solution, which is obtained through a flexible threshold filter based on the weighted spectral component values; the update formula is: ； in , From the matrix The spectral decomposition process performs threshold shrinkage on each spectral component value, with the threshold value being proportional to the corresponding weight. Sub-problem update: Fix other variables, The update subproblem is a least squares problem that solves the ℓ1 norm structured penalty; this subproblem is also solved using flexible threshold filtering. The updated formula is: ; This allows for independent soft thresholding of each element of the matrix, where... ; This is a sparse structured penalty parameter, with a value range of [0.01, 0.1]. Sub-problem update: Fix other variables, The update subproblem is a quadratic optimization problem, and its optimality condition corresponds to a linear system: ； in The adjoint operator for constructing the Hankel matrix corresponds physically to the Hankel averaging operation; the linear system is solved efficiently using the conjugate gradient method; when the wave field is coupled to the mapping operator... Dao Collection Recombination Operator Sparse transformation operators When all operators are orthogonal or quasi-orthogonal, the coefficient matrix has a good condition number, and the conjugate gradient method converges quickly. The linear system can be further simplified to: ； in Constructing adjoint operators for operators of Hankel matrices; Lagrange multiplier update: Update the Lagrange multipliers according to the degree of constraint violation. The update formula is as follows: ； ； Weight Update: The adaptive weights are updated based on the distribution of spectral component values ​​in the current iteration result. The update formula is as follows: ； The weights are gradually adapted to the actual structural characteristics of the data. During algorithm initialization, the estimated signal Initialization is performed using the adjoint operation of the wave-field coupling mapping operator, i.e. The Lagrange multipliers are initialized to zero, and the weights are initialized to 1; the suppression parameters... The value is 1.0.

Images