Seismic sequence stratigraphy based on aurora-optimized two-dimensional stepwise variational mode decomposition

Abstract

The application belongs to the field of geophysical exploration, and particularly relates to a method for seismic sequence division based on aurora optimization two-dimensional step-by-step variation mode decomposition, comprising the following steps: obtaining an original two-dimensional seismic signal; based on the original two-dimensional seismic signal, an optimal penalty factor upper limit is solved by using an aurora optimization algorithm; based on the optimal penalty factor upper limit, the original two-dimensional seismic signal is subjected to two-dimensional step-by-step variation mode decomposition to obtain multi-order two-dimensional modes with different central frequencies; based on the multi-order two-dimensional modes, a frequency feature profile is calculated and generated, a cycle curve is extracted, and a sequence division of strata is performed. The method of the application extracts a real mode while ensuring the continuity of a phase axis, avoids a pseudo mode and mode aliasing, and obtains a sequence division result corresponding to a well logging curve cycle analysis result.

Description

Technical Field

[0001] This invention belongs to the field of geophysical exploration, specifically involving a seismic sequence decomposition method based on auroral-optimized two-dimensional progressive variational mode decomposition. Background Technology

[0002] Traditional methods often use time-frequency analysis for seismic sequence delineation. These methods extract specific frequency slices through spectral decomposition, effectively characterizing the lateral distribution of sedimentary bodies. However, their physical nature is still constrained by the Heisenberg uncertainty principle—there is a theoretical limit to time-frequency resolution, and the preset basis functions are difficult to match the non-stationary characteristics of seismic signals, resulting in a weakened physical correlation between frequency division energy and sedimentary cycle sequence.

[0003] To overcome the physical constraints of traditional time-frequency analysis methods, those skilled in the art have turned to mode decomposition methods. Successive variational mode decomposition (SVMD) is a commonly used mode decomposition method. By extracting modes sequentially, it avoids the need for a pre-defined number of modes and reduces inter-mode overlap, thus improving the accuracy and efficiency of signal analysis. SVMD, based on variational mode extraction, is a mode decomposition method with advantages such as high efficiency and accurate decomposition, and is widely used. One-dimensional SVMD performed channel-by-channel is a simple and effective analysis method in sequence partitioning; however, the characteristic of decomposing each channel separately leads to reduced decomposition efficiency, and the final two-dimensional profile is laterally discontinuous, making it unsuitable for sequence partitioning.

[0004] Therefore, there is an urgent need to introduce a seismic sequence delineation method that can solve the problem of inaccurate cycle identification and sequence delineation in two-dimensional seismic profiles using traditional time-frequency analysis and mode decomposition methods. Summary of the Invention

[0005] Based on this, the present invention proposes a seismic sequence decomposition method based on aurora-optimized two-dimensional progressive variational mode decomposition.

[0006] In a first aspect, embodiments of this application provide a seismic sequence decomposition method based on aurora-optimized two-dimensional stepwise variational mode decomposition, including:

[0007] S1, acquire the raw two-dimensional seismic signal;

[0008] S2, based on the original two-dimensional seismic signal, uses the aurora optimization algorithm to solve for the upper limit of the optimal penalty factor;

[0009] S3, based on the upper limit of the optimal penalty factor, performs two-dimensional stepwise variational mode decomposition on the original two-dimensional seismic signal to obtain multi-order two-dimensional modes with different center frequencies;

[0010] S4. Based on the multi-order two-dimensional modes, calculate and generate frequency characteristic profiles, extract cycle curves, and perform stratigraphic sequence division.

[0011] In one possible implementation, step S2 specifically involves:

[0012] S2.1, Set the number of candidate solutions for the upper limit of the optimal penalty factor and the maximum number of evaluations in the aurora optimization algorithm, and initialize the high-energy particle swarm;

[0013] S2.2, taking the original two-dimensional seismic signal as input, execute the aurora optimization algorithm, use any candidate solution of the upper limit of the optimal penalty factor in the high-energy particle swarm to perform a temporary two-dimensional stepwise variational mode decomposition on the original two-dimensional seismic signal, decompose to obtain two-dimensional modes of all orders, calculate the two-dimensional information entropy of each order of two-dimensional mode in turn, and take the minimum value of the two-dimensional information entropy as the fitness function value of this iteration.

[0014] S2.3, Repeat step S2.2 until the fitness function values ​​corresponding to all candidate solutions with the upper limit of the optimal penalty factor in the high-energy particle swarm are obtained;

[0015] S2.4, update the high-energy particle swarm based on all fitness function values;

[0016] S2.5, repeat steps S2.2-S2.4. After multiple iterations under the constraint of the maximum number of evaluations, the fitness function value converges to the minimum value. The candidate solution corresponding to the upper limit of the optimal penalty factor at this time is the upper limit of the optimal penalty factor.

[0017] In one possible implementation, the two-dimensional information entropy of the two-dimensional modality. The expression is:

[0018]

[0019] in, Indicates the first Two-dimensional modes; and They represent the first The longitudinal and lateral scales of the two-dimensional modality; Indicates the first Two-dimensional modes in spatial coordinates The amplitude at that point; Indicates the first The amplitude probability distribution of the two-dimensional mode.

[0020] In one possible implementation, step S3 specifically involves:

[0021] S3.1, perform two-dimensional stepwise variational mode decomposition on the original two-dimensional seismic signal, set the upper limit parameter of the penalty factor to the optimal upper limit of the penalty factor, and set the penalty factor... Set to the preset minimum value;

[0022] S3.2, based on penalty factor Alternately iteratively update the first The iteration continues until the convergence condition is met or the maximum number of iterations is reached, then the iteration stops, yielding the second-order two-dimensional modal components and their corresponding center frequencies and Lagrange multipliers. Two-dimensional modal components and their corresponding center frequencies and Lagrange multipliers;

[0023] S3.3, Regarding the penalty factor Perform exponential growth updates;

[0024] S3.4, Repeat steps S3.2-S3.3 until the penalty factor is reached. When the upper limit of the optimal penalty factor is exceeded, the loop terminates, and the two-dimensional modal components output in the last iteration are the first... Two-dimensional modal frequency domain analytic signal and its corresponding center frequency , will the Two-dimensional modal frequency domain analytic signal Perform completion and inverse Fourier transform to obtain the first... Two-dimensional modes;

[0025] Step S3.5, repeat steps S3.1-S3.4L times to obtain the first to Lth two-dimensional modes and their corresponding center frequencies.

[0026] In one possible implementation, the first The iterative formula for the two-dimensional modal components is:

[0027]

[0028] in, Indicates the first The iteration of the ... Two-dimensional modal components, , Indicates the maximum number of iterations; Indicates the first The next iteration of the two-dimensional planar Hilbert mask, , This represents the frequency matrix of the original two-dimensional seismic signal. Indicates the first The iteration of the ... The center frequencies corresponding to the two-dimensional modal components are express and The sign function of the inner product; Indicates an analytical signal; The frequency domain signal representing the original two-dimensional seismic signal; Indicates the penalty factor; Indicates the first The iteration of the ... Two-dimensional modal components; Indicates the first The iteration of the ... Lagrange multipliers corresponding to two-dimensional modal components; Indicates the first The center frequency corresponding to the two-dimensional mode.

[0029] In one possible implementation, the first The iterative formula for the center frequency corresponding to the two-dimensional modal component is:

[0030]

[0031] in, Indicates the first The iteration of the ... The center frequency corresponding to the two-dimensional modal components.

[0032] In one possible implementation, the first The iterative formula for the Lagrange multipliers corresponding to the two-dimensional modal components is:

[0033]

[0034] in, Indicates the first The iteration of the ... The Lagrange multipliers corresponding to the two-dimensional modal components are Indicates the first The iteration of the ... Lagrange multipliers corresponding to two-dimensional modal components; This represents the update parameters of the Lagrange multipliers; Indicates the first The sum of the two-dimensional modal frequency domain analytic signals of all orders prior to order 1; Indicates the first The iteration of the ... Two-dimensional modal components.

[0035] In one possible implementation, step S4 specifically involves: calculating the energy-weighted average frequency at each spatial location point in the original two-dimensional seismic signal, traversing all spatial points to obtain a frequency feature profile that can characterize the local frequency changes of the original two-dimensional seismic signal; extracting a single-channel signal from the frequency feature profile as a cyclic curve, performing sedimentary cyclic analysis, and completing the stratigraphic sequence division.

[0036] In one possible implementation, the energy-weighted average frequency at each spatial location point The calculation formula is:

[0037]

[0038] in, Indicates the total number of two-dimensional modes; Indicates the first Two-dimensional mode at point The amplitude at that point; Indicates the first The center frequency corresponding to the two-dimensional mode.

[0039] This invention proposes and implements a seismic sequence decomposition method based on aurora-optimized two-dimensional stepwise variational mode decomposition, based on the application of stepwise variational mode decomposition in sequence stratigraphy. Unlike existing stepwise variational mode decomposition techniques, this method extends the one-dimensional analytic signal before mode decomposition, proposes a definition for the two-dimensional analytic signal, and adds lateral consistency constraints to the optimization conditions to ensure lateral consistency in the decomposition of the two-dimensional signal. Finally, aurora optimization is applied to the upper limit of the input penalty factor to avoid pseudo-modes or mode aliasing, further improving the decomposition effect on the two-dimensional signal. Experiments using synthetic images demonstrate the feasibility of this method.

[0040] This invention introduces the Polar Lights Optimizer (PLO) algorithm to adaptively optimize key parameters of two-dimensional stepwise variational mode decomposition (2D-SVMD). This method integrates rotational motion, elliptical motion, and particle collision perturbation strategies, possessing good global exploration capabilities and local convergence performance. By simulating the complex dynamic behavior of high-energy particles, it effectively avoids getting trapped in local optima, achieving efficient optimization of two-dimensional mode decomposition parameters, thereby enhancing the frequency focusing and spatial structure fidelity of modes, and further improving the decomposition quality in seismic stratigraphic pattern recognition.

[0041] The results obtained by using real seismic datasets demonstrate the effectiveness of aurora-optimized two-dimensional stepwise variational mode decomposition. While ensuring the continuity of the phase axis, it extracts the true modes, avoids pseudo-modes and mode aliasing, and obtains sequence division results that correspond well with the well logging curve cycle analysis results.

[0042] The seismic sequence decomposition method based on aurora-optimized two-dimensional stepwise variational mode decomposition proposed in this application is simple and direct. It is a signal decomposition of two-dimensional seismic data, which is only related to the properties of the data itself. It has the advantages of high computational efficiency and strong practicality. Attached Figure Description

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

[0044] Figure 1 A flowchart of a seismic sequence decomposition method based on aurora-optimized two-dimensional stepwise variational mode decomposition provided in an embodiment of the present invention;

[0045] Figure 2 The original two-dimensional seismic signal synthesized image provided in the embodiments of the present invention;

[0046] Figure 3 The fitness function value convergence curve of the synthesized image based on the original two-dimensional seismic signal provided in the embodiments of the present invention;

[0047] Figure 4 This is a schematic diagram of all effective modes of the synthesized image based on the original two-dimensional seismic signal, provided in the embodiments of this application.

[0048] Figure 5 This is a schematic diagram of the frequency feature profile of a synthesized image based on the original two-dimensional seismic signal, provided in an embodiment of this application.

[0049] Figure 6 Post-stack seismic dataset provided in the embodiments of this application;

[0050] Figure 7 The fitness function value convergence curve based on the post-stack seismic dataset provided in the embodiments of this application;

[0051] Figure 8 This is a schematic diagram of all valid modes based on the post-stack seismic dataset provided in the embodiments of this application;

[0052] Figure 9 This is a schematic diagram of the frequency feature profile based on a post-stack seismic dataset provided in an embodiment of this application;

[0053] Figure 10 This is a comparison chart of cycle curves extracted from post-stack seismic datasets and cycle curves obtained from well logging data, provided in an embodiment of the present invention. Figure 10 In the image, 'a' represents the cycle curve extracted from the post-stack seismic dataset. Figure 10 In the middle, b represents the cyclic curve obtained from well logging data. Detailed Implementation

[0054] To better understand the technical solution of this application, the embodiments of this application will be described in detail below with reference to the accompanying drawings.

[0055] It should be understood that the described embodiments are merely some, not all, of the embodiments in this application. All other embodiments obtained by those skilled in the art based on the embodiments in this application without inventive effort are within the scope of protection of this application.

[0056] The terminology used in the embodiments of this application is for the purpose of describing particular embodiments only and is not intended to be limiting of this application. The singular forms “a,” “the,” and “the” used in the embodiments of this application and the appended claims are also intended to include the plural forms unless the context clearly indicates otherwise.

[0057] It should be understood that the term "and / or" used in this article is merely a description of the relationship between related objects, indicating that three relationships can exist. For example, A and / or B can represent: A existing alone, A and B existing simultaneously, or B existing alone. Additionally, the character " / " in this article generally indicates that the preceding and following related objects have an "or" relationship.

[0058] See Figure 1 This is a flowchart of a seismic sequence decomposition method based on aurora-optimized two-dimensional stepwise variational mode decomposition provided in an embodiment of the present invention. Figure 1 As shown, the method includes:

[0059] Step S1: Obtain the original two-dimensional seismic signal.

[0060] See Figure 2 This is a synthesized image of the original two-dimensional seismic signal provided in an embodiment of the present invention. Figure 2 As shown, the horizontal and vertical resolutions of the synthetic image are both 256. The original two-dimensional seismic signal is composed of six ellipses with different center frequencies and spatial orientations and a rectangle superimposed. It specifically includes a DC (zero frequency) elliptical component and a DC rectangular component to simulate the complex frequency components and geological morphology in the seismic signal.

[0061] Step S2 involves using the aurora optimization algorithm, based on the original 2D seismic signal, to solve for the upper limit of the optimal penalty factor. Specifically:

[0062] Step S2.1: Set the input parameters for the aurora optimization algorithm and initialize the high-energy particle swarm optimization. The input parameters include the number of candidate solutions for the upper limit of the optimal penalty factor. E Maximum number of evaluations And the upper limit of the variable for the optimal penalty factor. and lower limit of variables .

[0063] High-energy particle swarm optimization in matrix The formal representation of a matrix The elements in the matrix represent candidate solutions for the upper bound of the optimal penalty factor. The expression is:

[0064]

[0065] in, The lower bound of the variable representing the upper limit of the optimal penalty factor. The upper limit of the variable representing the upper limit of the optimal penalty factor. and This constitutes the boundary of the variable space that forms the upper limit of the penalty factor; This represents a sequence of random numbers ranging from 0 to 1. E This represents the number of candidate solutions for the upper bound of the optimal penalty factor, which is the number of rows in the initial high-energy particle swarm in the Aurora optimization algorithm. E The larger the value, the longer the running time. In this embodiment, E The value is set to 50; This represents the expandable dimension of the solution space, which is the number of optimization variables. In this embodiment, only the upper limit of the penalty factor is optimized, therefore there is only one optimization variable. .

[0066] Therefore, the matrix of the high-energy particle swarm in this embodiment The expression is:

[0067]

[0068] Furthermore, when initializing the high-energy particle swarm, candidate solutions for the upper limit of the optimal penalty factor in the high-energy particle swarm are generated based on pseudo-random numbers in the solution space.

[0069] Step S2.2: Using the original two-dimensional seismic signal as input, execute the aurora optimization algorithm. Use any candidate solution with the upper limit of the optimal penalty factor in the high-energy particle swarm to perform a temporary two-dimensional progressive variational mode decomposition on the original two-dimensional seismic signal. Decompose to obtain two-dimensional modes of all orders. Calculate the two-dimensional information entropy of each order of two-dimensional mode in turn. Use the minimum value of the two-dimensional information entropy as the fitness function value of this iteration.

[0070] Two-dimensional information entropy of two-dimensional modes The expression is:

[0071]

[0072] in, Indicates the first Two-dimensional modes; and They represent the first The longitudinal and transverse scales of the second-order two-dimensional mode represent the first, second, third, and fourth dimensions, respectively. The number of rows and columns of the two-dimensional mode; Indicates the first Two-dimensional modes in spatial coordinates The amplitude at that point; Indicates the first The amplitude probability distribution of the two-dimensional mode.

[0073] Two-dimensional information entropy of two-dimensional modes Able to quantify the first The energy distribution sparsity of the two-dimensional modes.

[0074] Step S2.3, repeat step S2.2 until the fitness function values ​​corresponding to all candidate solutions with the upper limit of the optimal penalty factor in the high-energy particle swarm are obtained;

[0075] Step S2.4: Update the high-energy particle swarm based on all fitness function values.

[0076] It is particularly important to note that the update process for high-energy particle swarm optimization needs to consider both rotational and elliptical motion searches. The rotational motion search first performs a fine-grained search along the gradient direction of the current optimal solution, while the elliptical motion search introduces random disturbances. To avoid getting trapped in local optima, collisions and disturbances between particles need to be considered throughout the process.

[0077] Step S2.5, repeat steps S2.2-S2.4, up to the maximum number of evaluations. Under the constraint, after multiple iterations, the fitness function value converges to its minimum value. The candidate solution corresponding to this minimum value is the optimal upper bound of the penalty factor. See also... Figure 3 This is the convergence curve of the fitness function value of the synthesized image based on the original two-dimensional seismic signal provided in the embodiment of the present invention.

[0078] In this embodiment, the maximum number of evaluations .

[0079] It should be noted that, E The value is user-defined, that is, how many candidate solutions are tried in one iteration. E Defined as 50, 50 candidate solutions are randomly generated in the solution space based on pseudo-random numbers. Two-dimensional stepwise variational mode decomposition is performed on the original two-dimensional seismic signal using the 50 candidate solutions, and 50 minimum two-dimensional information entropy values ​​(i.e., fitness function values) are obtained. Then, the high-energy particle swarm optimization algorithm is used to update the rotational and elliptical motion search, and 50 better candidate solutions are regenerated. This process is repeated until the maximum number of evaluations is reached, and the upper limit of the penalty factor that minimizes the two-dimensional information entropy is found, which is the optimal penalty factor upper limit.

[0080] Step S3: Based on the upper limit of the optimal penalty factor, perform two-dimensional stepwise variational mode decomposition on the original two-dimensional seismic signal to obtain multiple two-dimensional modes with different center frequencies. Specifically:

[0081] Step S3.1: Perform two-dimensional stepwise variational mode decomposition on the original two-dimensional seismic signal, set the upper limit parameter of the penalty factor to the optimal upper limit of the penalty factor, and set the penalty factor... Set to the preset minimum value.

[0082] Step S3.2, based on the penalty factor Alternately iteratively update the first The iteration continues until the convergence condition is met or the maximum number of iterations is reached, then the iteration stops, yielding the second-order two-dimensional modal components and their corresponding center frequencies and Lagrange multipliers. Two-dimensional modal components and their corresponding center frequencies and Lagrange multipliers.

[0083] No. The iterative formula for the two-dimensional modal components is:

[0084]

[0085] in, Indicates the first The iteration of the ... Two-dimensional modal components, , This represents the maximum number of iterations, in this embodiment. ; Indicates the first The next iteration uses a two-dimensional planar Hilbert mask to construct a two-dimensional analytic signal in the frequency domain. , This represents the frequency matrix of the original two-dimensional seismic signal. Indicates the first The iteration of the ... The center frequencies corresponding to the two-dimensional modal components are A zero matrix indicates that modes are extracted starting from frequency 0. express and The sign function of the inner product; Indicates an analytical signal; The frequency domain signal representing the original two-dimensional seismic signal is obtained by two-dimensional Fourier transform of the original two-dimensional seismic signal; The penalty factor is a dynamically adjusted quantity. Indicates the first The iteration of the ... Two-dimensional modal components, A zero matrix; Indicates the first The iteration of the ... The Lagrange multipliers corresponding to the two-dimensional modal components are dynamically updated "error correction terms" during the iteration process, used to measure the current reconstruction error. In subsequent solution processes, this error correction term... It will also be iteratively updated to guide the optimization direction and ensure the fidelity of the final decomposition result; Indicates the first The center frequency corresponding to the two-dimensional mode.

[0086] No. The iterative formula for the center frequency corresponding to the two-dimensional modal component is:

[0087]

[0088] in, Indicates the first The iteration of the ... The center frequency corresponding to the two-dimensional modal components.

[0089] No. The iterative formula for the Lagrange multipliers corresponding to the two-dimensional modal components is:

[0090]

[0091] in, Indicates the first The iteration of the ... The Lagrange multipliers corresponding to the two-dimensional modal components are Indicates the first The iteration of the ... The Lagrange multipliers corresponding to the two-dimensional modal components are A zero matrix; The update parameters of the Lagrange multipliers are usually called the time-step or dual ascending step size. Their core function is to control the update speed of the Lagrange multipliers in each iteration. Indicates the first The sum of the two-dimensional modal frequency domain analytic signals of all orders prior to order 1; Indicates the first The iteration of the ... Two-dimensional modal components.

[0092] Step S3.3, regarding the penalty factor Perform exponential growth updates.

[0093] Step S3.4: Repeat steps S3.2-S3.3 until the penalty factor is reached. When the upper limit of the optimal penalty factor is exceeded, the loop terminates, and the two-dimensional modal components output in the last iteration are the first... Two-dimensional modal frequency domain analytic signal and its corresponding center frequency , will the Two-dimensional modal frequency domain analytic signal Perform completion and inverse Fourier transform to obtain the first... Two-dimensional modes;

[0094] Step S3.5, repeat steps S3.1-S3.4L times to obtain the first to Lth two-dimensional modes and their corresponding center frequencies.

[0095] See Figure 4 This diagram illustrates all effective modes of the synthesized image based on the original two-dimensional seismic signal, as provided in this embodiment of the application. Two-dimensional successive variational mode decomposition (2D-SVMD) decomposes the input signal into six identical two-dimensional modes (or simply modes), each representing a different frequency band component of the signal. In the mode decomposition results of the synthesized image, the first-order mode represents the boundary of the DC component, the second-order mode represents the DC component in the synthesized image, and the third to sixth-order modes represent ellipses of different frequencies that make up the synthesized image. Experimental results clearly demonstrate that the third to sixth-order modes successfully separated the elliptical components of different spatial frequencies in the original image. This verifies that each mode decomposed by 2D-SVMD has a clear physical meaning, and its center frequency exhibits a monotonically increasing trend with the increase of the mode number.

[0096] Step S4: Based on the multi-order two-dimensional modes, calculate and generate frequency characteristic profiles, extract cyclic curves, and perform stratigraphic sequence division.

[0097] Specifically, the energy-weighted average frequency at each spatial location point in the original two-dimensional seismic signal is calculated, and all spatial points are traversed to obtain a frequency characteristic profile that can characterize the local frequency changes of the original two-dimensional seismic signal; single-channel signals are extracted from the frequency characteristic profile as cyclic curves, sedimentary cyclic analysis is performed, and stratigraphic sequence division is completed.

[0098] Energy-weighted average frequency at each spatial location point The calculation formula is:

[0099]

[0100] in, Indicates the total number of two-dimensional modes; Indicates the first Two-dimensional mode at point The amplitude at that point; Indicates the first The center frequency corresponding to the two-dimensional mode.

[0101] See Figure 5 This is a schematic diagram of the frequency feature profile of a synthesized image based on the original two-dimensional seismic signal, provided in an embodiment of this application.

[0102] Based on the above, a real post-stack seismic dataset is used as the input of the original two-dimensional seismic signal to verify the practicality and effectiveness of the method described in this invention when processing real seismic exploration data. See also... Figure 6 The post-stack seismic dataset provided in the embodiments of this application is as follows: Figure 6 As shown, the trace spacing in this dataset is 6.25 meters, and the number of stacks is 24. The profile reveals complex stratigraphic contact relationships and reflection characteristics, posing a challenge to sequence stratigraphy.

[0103] Parameter optimization: First, for the post-stack seismic dataset, the PLO (Propagation-Optimization-Lower-Low) algorithm is used, with the minimum two-dimensional information entropy of the decomposed modes as the objective function, to automatically find the optimal upper bound of the penalty factor suitable for this specific dataset. See also Figure 7 This is the convergence curve of the fitness function value based on the post-stack seismic dataset provided in the embodiments of this application.

[0104] The post-stack seismic dataset is input into the 2D-SVMD algorithm, and high-fidelity mode decomposition is performed using the upper limit of the optimal penalty factor found in the previous step. This yields a series of two-dimensional modes that can reflect information about geological bodies at different scales. (See [link to relevant documentation]). Figure 8 This diagram illustrates all effective modes based on the post-stack seismic dataset provided in this embodiment. Among the modes decomposed by 2D-SVMD, the first and second modes exhibit strong energy, steep dip angles, and overlapping linear or arc-shaped characteristics. This is not geological reflection but rather typical acquisition imprints or migration noise, severely interfering with the identification of effective signals. In stark contrast to these strong interferences, the third mode, as the core result of this decomposition, clearly reveals the main layered reflection structures in the study area with extremely high energy. Its phase axis is gentle and continuous, and its structural morphology conforms to geological laws, serving as a key basis for structural and stratigraphic interpretation. Meanwhile, the fifth mode, as a secondary signal component, reflects weaker, higher-frequency stratigraphic details. Furthermore, the algorithm successfully separates high-frequency random noise into the sixth mode, which exhibits typical spatial disorder distribution characteristics.

[0105] Based on the modal components obtained from the decomposition, the frequency characteristic profile of this seismic profile is calculated and generated. (See [link]) Figure 9 This is a schematic diagram of the frequency characteristic profile based on a post-stack seismic dataset provided in an embodiment of this application. From... Figure 9By extracting single-channel signals from the frequency characteristic profile as cyclic curves, sedimentary cycle analysis can be performed, ultimately completing the stratigraphic sequence division. See also... Figure 10 This is a comparison chart of the cycle curves extracted from post-stack seismic datasets and the cycle curves obtained from well logging data, provided in an embodiment of the present invention. Figure 10 In the image, 'a' represents the cycle curve extracted from the post-stack seismic dataset. Figure 10 In the middle, b represents the cycle curve obtained from well logging data, for comparison. Figure 10 a and Figure 10 As can be seen from Figure b, the peaks and troughs of the two curves show a high correlation, verifying the effectiveness of the seismic sequence decomposition method based on aurora-optimized two-dimensional stepwise variational mode decomposition in this invention for cyclic analysis and sequence decomposition.

[0106] The above description is merely a specific embodiment of the present invention. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this invention should be included within the protection scope of this invention. The protection scope of this invention should be determined by the scope of the claims.

Claims

1. A method for seismic sequence division based on aurora optimized two-dimensional progressive variational modal decomposition, characterized in that, The method comprises the following steps: S1, obtaining an original two-dimensional seismic signal; S2, based on the original two-dimensional seismic signal, using the aurora optimization algorithm to solve the optimal penalty factor upper limit; S3, based on the optimal penalty factor upper limit, performing two-dimensional step-by-step variational mode decomposition on the original two-dimensional seismic signal to obtain multi-order two-dimensional modes with different center frequencies; S4, based on the multi-order two-dimensional modes, calculating and generating a frequency feature profile, extracting a cycle curve, and performing stratigraphic sequence division; Step S2 is specifically: S2.1, set the number of optimal penalty factor upper limit candidate solutions in the aurora optimization algorithm and the maximum evaluation times, and initialize the high-energy particle swarm; S2.2, taking the original two-dimensional seismic signal as input, executing the aurora optimization algorithm, using any optimal penalty factor upper limit candidate solution in the high-energy particle swarm to perform a temporary two-dimensional step-by-step variational mode decomposition on the original two-dimensional seismic signal, and decomposing to obtain two-dimensional modes of all orders, and calculating the two-dimensional information entropy of each order two-dimensional mode in turn, and taking the minimum value of the two-dimensional information entropy as the fitness function value of this iteration; S2.3, repeat step S2.2 until the fitness function values corresponding to all optimal penalty factor upper limit candidate solutions in the high-energy particle swarm are obtained; S2.4, update the high-energy particle swarm based on all fitness function values; S2.5, repeat steps S2.2-S2.4, and after multiple iterations under the limitation of the maximum evaluation times, the fitness function value converges to the minimum value, at this time the optimal penalty factor upper limit candidate solution corresponding to the minimum value is the optimal penalty factor upper limit.

2. The method of claim 1, wherein the method is based on the aurora-optimized two-dimensional stepwise variational modal decomposition. Two-dimensional information entropy of two-dimensional modal The expression is: ， wherein, represents the th two-dimensional mode; and represent the longitudinal and transverse scale of the th two-dimensional mode, respectively; represents the amplitude of the th two-dimensional mode at the spatial coordinate ; represents the amplitude probability distribution of the th two-dimensional mode.

3. The method of claim 1, wherein the method is based on the aurora-optimized two-dimensional stepwise variational modal decomposition. Step S3 is specifically: S3.1, performing two-dimensional step-by-step variational modal decomposition on the original two-dimensional seismic signal, setting the upper limit parameter of the penalty factor to the optimal upper limit of the penalty factor, and setting the lower limit parameter of the penalty factor to the preset minimum value; setting the lower limit parameter of the penalty factor to the preset minimum value; S3.2, based on the penalty factor , alternately iteratively updating the first order two-dimensional modal components and the corresponding center frequencies and Lagrange multipliers until a convergence condition is met or a maximum number of iterations is reached, stopping the iteration, and obtaining the first order two-dimensional modal components and the corresponding center frequencies and Lagrange multipliers; S3.3, on the penalty factor exponential growth update; S3.4, repeat steps S3.2-S3.3 until the penalty factor When the optimal penalty factor upper limit is exceeded, the loop terminates, and the two-dimensional modal component output in the last iteration is the n order two-dimensional modal frequency-domain analytical signal and its corresponding center frequency , the n order two-dimensional modal frequency-domain analytical signal is completed and Fourier inverse transformed to obtain the n order two-dimensional modal Step S3.5, repeat steps S3.1-S3.4 for L times, and obtain the first-order to L-order two-dimensional modes and the corresponding center frequencies in turn.

4. The method of claim 3, wherein the method is based on the aurora-optimized two-dimensional stepwise variational modal decomposition. No. The iterative formula for the two-dimensional modal components is: ， wherein, represents the order two-dimensional modal component of the th iteration, , represents the maximum number of iterations; represents the two-dimensional plane Hilbert mask of the , represents the frequency matrix of the original two-dimensional seismic signal, represents the order two-dimensional modal component of the th iteration, represents the sign function of the inner product with ; represents the analytic signal; represents the frequency domain signal of the original two-dimensional seismic signal; represents the penalty factor; represents the order two-dimensional modal component of the th iteration; represents the order two-dimensional modal component of the th iteration; represents the center frequency of the order two-dimensional modal.

5. The method of claim 4, wherein the method is based on the aurora-optimized two-dimensional stepwise variational modal decomposition. No. The iterative formula for the center frequency corresponding to the two-dimensional modal component is: ， wherein, represents the center frequency corresponding to the second-order two-dimensional modal component of the nth iteration.

6. The method of claim 5, wherein the method is based on the aurora-optimized two-dimensional stepwise variational modal decomposition. No. The iterative formula for the Lagrange multipliers corresponding to the two-dimensional modal components is: ， wherein, represents the Lagrange multiplier corresponding to the dimensional modal component of the dimensional modal component of the dimensional modal component of the dimensional modal component of the dimensional modal component of the represents an update parameter of the Lagrange multiplier; represents the sum of the dimensional modal component of the represents the sum of the dimensional modal component of the dimensional modal component of the 7. The method of claim 1, wherein the method is based on the aurora-optimized two-dimensional stepwise variational modal decomposition. Step S4 is specifically: calculate the energy-weighted average frequency of each spatial position point in the original two-dimensional seismic signal, traverse all spatial points to obtain a frequency feature profile capable of representing the local frequency variation of the original two-dimensional seismic signal; extract a single-channel signal from the frequency feature profile as a cycle curve, perform sedimentary cycle analysis, and complete stratigraphic sequence division.

8. The method of claim 7, wherein the method is based on the aurora-optimized two-dimensional stepwise variational modal decomposition. Energy weighted average frequency at each spatial location point The formula for calculating the energy weighted average frequency is: ， wherein, denotes the total number of two-dimensional modes; denotes the amplitude of the two-dimensional mode of order at the point denotes the center frequency corresponding to the two-dimensional mode of order .​

Images