Multi-scale deposition cycle phase space fine characterization method, system and application

Abstract

The invention discloses a multi-scale deposition cycle phase space fine characterization method and system and application, the method combines a variational mode decomposition method and synchronous extrusion optimal basis wavelet transform, and provides a working process for multi-scale deposition cycle characterization. The method comprises the following steps: firstly, giving a relationship between geological structures with different scales and intrinsic mode functions through theoretical derivation, and proposing a method for determining the number of the intrinsic mode functions; in consideration of the problem of transverse continuity of calculation efficiency and a decomposition result, a variational mode decomposition algorithm is improved, namely, the number of intrinsic mode functions and the center frequency of the intrinsic mode functions are determined according to a smooth reflection coefficient amplitude spectrum. On this basis, multi-scale deposition cycle characterization is realized by combining synchronous extrusion optimal basis wavelet transform and ridge extraction methods; the key parameters in the processing process are all extracted from the seismic data, operation is easy, the calculation speed is high, the method is suitable for massive seismic data processing, and an important tool is provided for subsequent sequence stratigraphy research.

Description

Technical Field

[0001] The present invention belongs to the technical field of petroleum geophysical exploration, and is mainly a multi-scale sedimentary cycle phase space fine characterization method, system and application for auxiliary sequence stratigraphy research. Background Art

[0002] The basic argument of sequence stratigraphy is that the geometry and lithology of stratigraphic units are controlled by four major factors: sea level rise and fall, tectonic subsidence, sediment supply, and climate (McLaughlin et al., 2005). These factors act at different scales and show periodic changes, so sedimentary sequences also have multi-scale and periodic characteristics. Sedimentary cycles are the basic building blocks of sedimentary sequences and are the result of the superposition of periodic sedimentary events at different scales. Regional or global sedimentary events affect the superposition pattern of sedimentary sequences at different scales. Sedimentary cycles are generally composed of some basic units, including transgressive cycle units and regressive cycle units. For transgressive cycles, the thickness of sandstone increases with depth, while the grain size of sediment changes from fine to coarse, indicating a transgressive sedimentary sequence. For regressive cycles, the thickness of sandstone decreases with depth, while the grain size of sediment changes from coarse to fine, indicating a regressive sedimentary sequence (Melani et al., 2020; Rossi et al., 2020; Zhang et al., 2017). Typically, a large sequence is composed of several small-scale cycles, recording evidence of transgressive and regressive sequences in the sedimentary environment (Goldhammer et al., 1990; Melani et al., 2020). Therefore, the detailed characterization and division of multi-scale sedimentary cycles is the basis for subsequent sedimentary environment research, sequence stratigraphic analysis and reservoir prediction, and is also an important step in establishing a structural and stratigraphic framework.

[0003] Well logging technology is the most common geophysical tool for characterizing and dividing multi-scale sedimentary cycles by making fine measurements of underground stratigraphic media in the wellbore (Vermeer et al., 1992; Coconi-Morales et al., 2010). However, well logging data is only a view of one well and its spatial coverage is limited. In contrast, seismic reflection data has better spatial coverage, carries lithology and physical property information of underground strata, and contains a large amount of information about multi-scale sedimentary cycles (Zeng, 2013). Therefore, effectively extracting sedimentary cycle information from seismic data is of great significance for the identification and interpretation of sedimentary cycles.

[0004] Data-driven decomposition methods have been used to extract sedimentary cycle information from seismic data. Liu et al. (2015) used the empirical mode decomposition method to decompose seismic signals into intrinsic mode functions representing sedimentary cycles of different scales to assist in the characterization of multi-scale sedimentary cycles. However, the empirical mode decomposition method cannot artificially determine the optimal number of intrinsic mode functions. To address this problem, Li et al. (2017) introduced a method based on variational mode decomposition, which can artificially define the number of intrinsic mode functions. Although this method provides more flexibility, the interpretation of its results is still highly subjective and susceptible to the experience of the interpreter. Therefore, a key issue that these data-driven decomposition methods need to solve is to establish a deterministic relationship between intrinsic mode functions and geological structures of different scales. In addition, more reasonable criteria are needed to determine the optimal number of intrinsic mode functions (Tian et al., 2022).

[0005] In order to address the limitations of previous work, the applicant previously proposed an enhanced multi-channel variational modal decomposition method and a workflow for seismic data decomposition to assist in multi-scale sequence stratigraphic interpretation (Tian et al., 2022). The method first clarified the relationship between intrinsic mode functions and geological structures at different scales. Subsequently, a method for determining the number of intrinsic mode functions based on scale space representation and compression mapping operators was proposed. Finally, a multi-channel variational modal decomposition method was used to decompose seismic data to assist in the interpretation of sedimentary sequence boundaries at different scales. However, the decomposition results of these modal decomposition methods, i.e., seismic reflection waveforms at different scales, still need to be interpreted by experienced interpreters to infer sedimentary sequences, which hinders the direct division of sedimentary cycle units.

[0006] The time-frequency analysis method can intuitively reflect the instantaneous main frequency change characteristics of the seismic signal caused by the change of thin layer thickness within the sedimentary cycle (Tian et al., 2022). Therefore, the time-frequency analysis method is also one of the most commonly used tools for characterizing sedimentary cycles. Accurately characterizing the main frequency change characteristics caused by the change of thin layer thickness is the key issue in using the time-frequency analysis method to assist in the division of sedimentary cycles. The pioneering research in this field can be traced back to 1980, when the sedimentary cycle division method based on linear time-frequency transformation was proposed by the Russian Institute of Exploration Geophysics (Mushin et al., 2000).

[0007] However, the time-frequency spectrum of the linear time-frequency analysis method diffuses along the frequency axis on the time-frequency plane (Huang et al., 2016), resulting in low time-frequency resolution and hindering the identification of sedimentary cycles. To address this defect, Daubechies et al. (2011) introduced the synchronous squeezing transform. The synchronous squeezing transform obtains a high-resolution time-frequency spectrum by synchronously squeezing the time-frequency spectrum along the frequency axis. Although this method provides a powerful tool for sedimentary cycle analysis, the choice of time-frequency atoms determines the interpretability of the obtained results. Gao et al. (2001) believed that when the base wavelet is very similar to the seismic wavelet, the time-frequency spectrum of the continuous wavelet transform has stronger noise robustness and readability. In order to better match the seismic wavelet, Gao et al. (2006) proposed a base wavelet called the three-parameter wavelet. By adjusting its three parameters, its waveform can be matched with the seismic wavelet of the seismic data to be analyzed. With matching the seismic wavelet as the criterion, we proposed an optimal base wavelet construction method combining the three-parameter wavelet and deep neural network. The method first estimates the seismic wavelet using an alternating iterative deep neural network, and then uses a three-parameter wavelet to fit the seismic wavelet to obtain the optimal base wavelet. Although the optimal base wavelet transform achieves high-resolution spectral decomposition of seismic data (Tian et al., 2021), its frequency resolution is low, which is not conducive to the characterization of the main frequency characteristics of the sedimentary cycle.

[0008] In order to solve the above limitations, we have previously combined the synchronous extrusion transform with the optimal base wavelet and proposed a synchronous extrusion optimal basic wavelet transform to achieve the optimal characterization of sedimentary cycles. This method achieves the optimal representation of the main frequency characteristics of seismic data by defining the main frequency positioning conditions of the synchronous extrusion transform and the similarity coefficient conditions between the base wavelet and the seismic wavelet, and then achieves the optimal characterization of sedimentary cycles (Tian et al., 2022). However, this method directly processes seismic data containing multi-scale sequence stratigraphic information, resulting in the superposition of information of sedimentary cycles of different scales in the time-frequency spectrum. The use of this time-frequency spectrum that simultaneously contains multi-scale sedimentary sequence information will reduce the accuracy of subsequent multi-scale sedimentary cycle unit division.

[0009] The above prior art has the following disadvantages:

[0010] (1) The existing modal decomposition method is used to decompose seismic data, which results in seismic reflection waveforms of different scales containing geological structure information of different scales. The results still need to be interpreted by experienced interpreters to infer the sedimentary laws, which will introduce human experience errors and is not convenient for the direct division of sedimentary cycle units.

[0011] (2) The existing sedimentary cycle characterization method based on synchronous extrusion optimal basis wavelet transform directly processes seismic data containing multi-scale sedimentary sequence information, resulting in the time-frequency spectrum containing sedimentary cycle information of different scales at the same time. Using this mixed time-frequency spectrum of multi-scale sedimentary cycle information will reduce the accuracy of subsequent multi-scale sedimentary cycle unit division. Summary of the invention

[0012] In order to overcome the shortcomings of the above-mentioned prior art, the purpose of the present invention is to provide a method for fine characterization of the phase space of multi-scale sedimentary cycles, which combines the variational mode decomposition method with the synchronous extrusion optimal basis wavelet transform, and proposes a workflow for multi-scale sedimentary cycle characterization. First, the relationship between geological structures of different scales and intrinsic mode functions is given through theoretical derivation, and a method for determining the number of intrinsic mode functions is proposed. Taking into account the computational efficiency and the lateral continuity of the decomposition results, we improved the variational mode decomposition algorithm, that is, determining the number of intrinsic mode functions and their center frequencies based on the smoothed reflection coefficient amplitude spectrum. On this basis, a new method for extracting the instantaneous main frequency attributes of seismic data is proposed in combination with the synchronous extrusion optimal basis wavelet transform and the ridge extraction method, thereby realizing the characterization of multi-scale sedimentary cycles.

[0013] In order to achieve the above object, the present invention adopts the following technical scheme to achieve it: a multi-scale sedimentary cycle phase space fine characterization method, comprising the following steps:

[0014] Calculate the average Fourier amplitude spectrum of 2D seismic data;

[0015] Calculate the average reflection coefficient amplitude spectrum;

[0016] The scale space representation method is used to calculate the smoothed reflection coefficient amplitude spectrum. According to the relationship between the reflection coefficient amplitude spectrum and the empirical mode, the number of extreme values ​​of the smoothed reflection coefficient spectrum is set to the number of intrinsic mode function components, and the extreme value frequency is set to the center frequency of each intrinsic mode function component.

[0017] The variational mode decomposition method is used to decompose the input two-dimensional seismic data to obtain the intrinsic mode function components that can characterize the seismic reflection characteristics of geological structures of different scales;

[0018] The synchronous squeezed optimal basis wavelet transform is used to calculate each seismic trace of each two-dimensional intrinsic mode function component to obtain the time-frequency spectrum of each seismic trace.

[0019] The ridge extraction method is used to extract the instantaneous main frequency curve from the time-frequency spectrum of each seismic trace to obtain a two-dimensional instantaneous main frequency attribute profile.

[0020] The instantaneous main frequency attribute profile is used to interpret the lateral distribution of sedimentary cycle units, that is, the detailed characterization of the multi-scale sedimentary cycle phase space.

[0021] Further, calculating the average reflection coefficient amplitude spectrum includes extracting the seismic wavelet amplitude spectrum from the average Fourier amplitude spectrum of the seismic data using a compression mapping operator, and obtaining the reflection coefficient amplitude spectrum by dividing the average Fourier amplitude spectrum of the seismic data by the seismic wavelet amplitude spectrum.

[0022] Furthermore, the compression mapping operator is used to extract the seismic wavelet amplitude spectrum from the seismic data amplitude spectrum;

[0023] The compression mapping operator for calculating the seismic sub-spectrum is introduced. Assume that the seismic sub-spectrum is within the specified frequency band [f l ,f h ] appears as a smooth and unimodal function:

[0024]

[0025] q∈(0,1), 0≤F q (|W|; f)≤1, let c q >0, α, β>1, the compression mapping operator P is

[0026]

[0027] Set the initial value of the iteration to:

[0028]

[0029] By iterating the compression mapping algorithm |W| k =P(|W| k-1 ) so that it converges to a fixed point |W| * , the fixed point corresponds to the estimated seismic sub-spectrum (|W| * ) 1 / q ;

[0030] The reflection coefficient amplitude spectrum |R is obtained by dividing the seismic data amplitude spectrum |S(f)| by the seismic wavelet amplitude spectrum |W(f| k (f)|:

[0031]

[0032] Here ε≠0, preventing |W(f)|=0.

[0033] Furthermore, according to the relationship between the reflection coefficient spectrum and the empirical mode, the number of extreme values ​​of the smoothed reflection coefficient spectrum is set to the number of intrinsic mode function components, and the extreme value frequency is the center frequency of each intrinsic mode function component. The relationship between the reflection coefficient amplitude spectrum and geological structures of different scales is obtained based on theoretical derivation, as shown in the following formula:

[0034]

[0035] The seismic reflection of structures of different scales is defined as the empirical mode function, where |R k (f) | represents the amplitude spectrum of the reflection coefficient corresponding to the geological structure at the kth scale, c kj represents the reflection coefficient value of the jth sampling point in the kth empirical mode function, (t km -t kn ) represents the time difference between two reflection coefficient points m and n at the kth scale. It is defined that a smaller absolute value k represents a large-scale structure, and a larger value k represents a small-scale structure. Equation (1) shows that |R k The frequency of the first extreme point of (f)| and (t km -t kn ) is inversely proportional to.

[0036] Furthermore, the input two-dimensional seismic data is decomposed by using the variational mode decomposition method to obtain the intrinsic mode function components that can characterize the seismic reflection characteristics of geological structures of different scales, including:

[0037] The signal separation and extraction is achieved by decomposing the signal into eigenmode functions with finite bandwidth, which is achieved by solving the following optimization problem:

[0038]

[0039] Among them, s(t) represents single-channel seismic data, u m (t) represents the mth intrinsic mode function, ω m represents the center frequency of the mth eigenmode function, represents the Hilbert transform of the mth intrinsic mode function. Further, the Lagrange multiplier method shown in the following formula is used to solve the above optimization problem:

[0040]

[0041] The value of the center frequency is determined by the frequency of the extreme points of the smoothed reflection coefficient amplitude spectrum.

[0042] Furthermore, the synchronous squeezed optimal basis wavelet transform is used to calculate each seismic trace of each two-dimensional intrinsic mode function component, and the time-frequency spectrum of each seismic trace is obtained, including:

[0043] The optimal basic wavelet is used as the base wavelet in the synchronous squeezing transform framework. The optimal basic wavelet is obtained by fitting the seismic wavelet with a three-parameter wavelet. The three-parameter wavelet is defined in the time domain as:

[0044]

[0045] Where σ, τ and β represent the modulation frequency, energy attenuation factor and energy delay factor respectively. Λ=(σ,τ,β) represents the set of σ, τ and β.

[0046]

[0047] Based on the concept of the above method, the present invention provides a multi-scale sedimentary cycle phase space fine characterization system, including a first calculation module, a second calculation module, a mapping module, a time-frequency spectrum calculation module and a characterization module;

[0048] The first calculation module is used to calculate the average Fourier amplitude spectrum and the average reflection coefficient amplitude spectrum of the two-dimensional seismic data;

[0049] The second calculation module calculates the smoothed reflection coefficient amplitude spectrum using a scale space representation method, and according to the relationship between the reflection coefficient amplitude spectrum and the empirical mode, sets the number of extreme values ​​of the smoothed reflection coefficient spectrum to the number of intrinsic mode function components, and the extreme value frequency is the center frequency of each intrinsic mode function component;

[0050] The mapping module uses the variational mode decomposition method to decompose the input two-dimensional seismic data to obtain the intrinsic mode function components that can characterize the seismic reflection characteristics of geological structures at different scales;

[0051] The time-frequency spectrum calculation module uses synchronous extrusion optimal basis wavelet transform to calculate each seismic trace of each two-dimensional intrinsic mode function component to obtain the time-frequency spectrum of each seismic trace.

[0052] The characterization module uses the ridge extraction method to extract the instantaneous main frequency curve from the time-frequency spectrum of each seismic trace to obtain a two-dimensional instantaneous main frequency attribute profile. The instantaneous main frequency attribute profile is used to explain the lateral distribution law of the sedimentary cycle unit, that is, a detailed characterization of the multi-scale sedimentary cycle phase space.

[0053] The present invention can also provide a computer device, including a processor and a memory, the memory is used to store a computer executable program, the processor reads the computer executable program from the memory and executes it, and when the processor executes the computer executable program, the multi-scale sedimentary cycle phase space fine characterization method described in the present invention can be implemented.

[0054] A computer-readable storage medium is also provided, in which a computer program is stored. When the computer program is executed by a processor, the multi-scale sedimentary cycle phase space fine characterization method of the present invention can be implemented.

[0055] In addition, the multi-scale sedimentary cycle phase space fine characterization method described in the present invention can also be used to identify and interpret sedimentary cycle information extracted from seismic data.

[0056] Compared with the prior art, the present invention has the following beneficial effects:

[0057] The present application takes into account that seismic reflection data simultaneously contains seismic reflection information of sedimentary sequences of different scales. It overcomes the problem that traditional methods cannot effectively separate and characterize the multi-scale sedimentary cycle information contained in seismic data. The present application can effectively separate the seismic reflection wave characteristics of different scales from the seismic data by improving the variational modal decomposition method, and on this basis, the synchronous extrusion optimal basis wavelet transform is used to accurately characterize the changing laws of sedimentary cycles of different scales. The instantaneous main frequency attribute proposed in combination with the ridge extraction method can effectively characterize the lateral distribution range of sedimentary cycles of different scales. The key parameters in the processing process of the present invention are all extracted from the data, the operation is simple, the calculation speed is fast, and it is suitable for the processing of massive seismic data, which provides an important tool for subsequent sequence stratigraphic research. BRIEF DESCRIPTION OF THE DRAWINGS

[0058] Figure 1 It is a flow chart of multi-scale sedimentary cycle characterization of the present invention.

[0059] Figure 2 is an example of model data; wherein, (a) is the established geological model containing two sedimentary cycle scales: large scale (one transgressive cycle unit, three regressive cycle units) and small scale (two transgressive cycle units, five regressive cycle units); (b) is the reflection coefficient corresponding to the model; (c) is the synthetic seismic record; (d) is the intrinsic mode function 1 calculated using the variational mode decomposition method; (e) is the intrinsic mode function 2 calculated using the variational mode decomposition method; (f) is the time-frequency spectrum of the synthetic seismic record; (g) is the time-frequency spectrum of the intrinsic mode function 1; (h) is the instantaneous main frequency curve of the intrinsic mode function 1; (i) is the time-frequency spectrum of the intrinsic mode function 2; (j) is the instantaneous main frequency curve of the intrinsic mode function 2.

[0060] Figure 3 is the amplitude spectrum of the model data; (a) shows the amplitude spectrum of the model earthquake record (dashed line), the reflection coefficient amplitude spectrum (solid line) and its smoothed spectrum (dotted line). (b) shows the reflection coefficient amplitude spectrum (solid line), the reflection coefficient smoothed spectrum (dotted line), the reflection coefficient amplitude spectrum of the intrinsic mode function 1 (dotted line) and the intrinsic mode function 2 (dashed line).

[0061] Figure 4 is an actual data example; where (a) is the seismic profile of the target work area; (c) and (e) are the intrinsic mode function 1 and intrinsic mode function 2 of the seismic profile decomposed by the variational mode decomposition method; (b), (d) and (f) are the instantaneous main frequency attributes of the seismic profile to be analyzed, the intrinsic mode function 1 and the intrinsic mode function 2, respectively; the curve at the well point represents the logging acoustic impedance curve.

[0062] Figure 5 The results of characterizing and dividing the two-scale sedimentary cycles at Well 2 using the proposed method clearly demonstrate the effectiveness of the proposed method in decomposing the characteristics of two-scale seismic signals. The synchronous extrusion optimal basis wavelet transform can effectively characterize the main frequency variation characteristics of sedimentary cycles of different scales, and the sedimentary cycle unit division results based on this can well match the logging lithology interpretation results. DETAILED DESCRIPTION

[0063] In order to enable people in the field of the present invention to better understand the scheme of the present invention, the technical scheme in the embodiment of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiment of the present invention. Obviously, the described embodiment is only a part of the embodiment of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work should fall within the scope of protection of the present invention.

[0064] Sedimentary cycle is a sequence concept that describes the repeated deposition of sedimentary processes and sedimentary conditions in the same order. It is manifested as the alternation of lithology and lithofacies, and is the way in which rock deposition and erosion interact in geological history. The formation and evolution of sedimentary cycles are affected by global climate change, plate tectonic movement, sea level change, and the speed of sound of river erosion. Seismic data contain rich information about underground sedimentary bodies. Using time-frequency analysis technology, one-dimensional time-domain seismic records can be transformed into two-dimensional time-frequency domains to extract more information. In the time-frequency diagram, the distribution and movement trend of amplitude energy groups are observed. Low-frequency information corresponds to coarse-grained rock components or thicker strata, and high-frequency information corresponds to finer-grained rock components or thinner strata. The type of sedimentary cycle can be determined based on the movement direction and shape of the amplitude energy group. Combined with geological background knowledge, the identified sedimentary cycle can be interpreted, the formation and evolution process of the sedimentary cycle can be analyzed, and the sedimentary cycle information can be used to guide reservoir prediction and oil and gas exploration.

[0065] The present invention is further described in detail below in conjunction with the accompanying drawings:

[0066] The multi-scale sedimentary cycle phase space fine characterization method of the present invention has the following process: Figure 1 As shown; specifically including the following steps:

[0067] 1) Obtain two-dimensional seismic data and calculate its average Fourier spectrum S(f) using Fourier transform.

[0068] 2) Calculate the average reflection coefficient amplitude spectrum of the original two-dimensional seismic data; including two steps:

[0069] ①Use the compression mapping operator to extract the seismic wavelet amplitude spectrum from the seismic data amplitude spectrum;

[0070] A compression mapping operator is introduced to calculate the seismic sub-spectrum, which assumes that the seismic sub-spectrum is within the specified frequency band [f l ,f h ] is a smooth and unimodal function. Definition

[0071]

[0072] q∈(0,1), obviously 0≤F q (|W|; f)≤1. Let c q >0, α, β>1, the compression mapping operator P can be defined as

[0073] P[|W|;f]=c q F q α (|W|; f)[1-F q (|W|; f)] β (3)

[0074] Set the initial value of the iteration to:

[0075]

[0076] By iterating the compression mapping algorithm |W| k =P(|W| k-1 ) so that it converges to a fixed point |W| * This fixed point corresponds to the estimated seismic wavelet spectrum (|W| * ) 1 / q .

[0077] ② Use the seismic data amplitude spectrum |S(f)| to divide the seismic wavelet amplitude spectrum |W(f)| to obtain the reflection coefficient amplitude spectrum |R k (f)|:

[0078]

[0079] Here ε≠0, preventing |W(f)|=0.

[0080] 3) Calculate the smoothed reflection coefficient amplitude spectrum using the scale space representation method:

[0081] The kernel function used in the scale space representation framework is:

[0082]

[0083] Here n represents the scale parameter, the reflection coefficient amplitude spectrum |R k The scale space representation of (f)| is defined as:

[0084] L(f;n)=|R k (f)|*G(f;n) (7)

[0085] Here * represents the convolution symbol, and the scale space representation result is the original reflection coefficient amplitude spectrum |R k (f)| global trend (i.e. smoothing result). As the scale parameter n increases, L(f;n) will become smoother and the number of maximum frequencies will decrease. Here, we use the central frequency f of the seismic data c To determine the parameter n:

[0086] n=μf c (8)

[0087] Where μ = 0.2 ~ 0.4, is an empirical parameter. c Determined by the following formula:

[0088]

[0089] For each seismic data set, the optimal value of n is determined, which directly affects the number of extreme values ​​of the scale space representation results. This dependence stems from the inherent characteristics of the seismic data itself.

[0090] According to the relationship between the reflection coefficient amplitude spectrum and the empirical mode, the number of extreme values ​​of the smoothed reflection coefficient amplitude spectrum is set to the number of intrinsic mode function components, and the extreme value frequency is set to the center frequency of each intrinsic mode function component.

[0091] 4) The variational mode decomposition method is used to decompose the input two-dimensional seismic data to obtain the intrinsic mode function components that can characterize the seismic reflection characteristics of geological structures of different scales.

[0092] Variational mode decomposition is an adaptive and non-recursive signal processing method that separates and extracts signals by decomposing them into a series of intrinsic mode functions with finite bandwidths. This is achieved by solving the following optimization problem:

[0093]

[0094] Among them, s(t) represents single-channel seismic data, u m (t) represents the mth intrinsic mode function, ω m represents the center frequency of the mth eigenmode function, represents the Hilbert transform of the mth intrinsic mode function. Further, the Lagrange multiplier method shown in the following formula is used to solve the above optimization problem:

[0095]

[0096] It should be noted that in the present invention, when solving the above formula, the center frequency of each eigenmode function is fixed and will not change with the iterative process. The value of the center frequency is determined by the frequency of the extreme point of the smoothed reflection coefficient amplitude spectrum. There are two reasons for this treatment: first, considering that the waveforms generated by structures of the same scale are similar, the different frequency components of the amplitude spectrum can represent the characteristics of geological structures of different scales. Secondly, in order to ensure the calculation efficiency and the lateral continuity of the calculation results, since variational mode decomposition is a single-channel calculation method, not fixing the center frequency when calculating two-dimensional seismic data will result in poor lateral continuity of the decomposition results.

[0097] 5) Using synchronous squeezing optimal basis wavelet transform, each seismic trace of each two-dimensional intrinsic mode function component is calculated to obtain the time-frequency spectrum.

[0098] Synchronous Squeezing Optimal Basis Wavelet Transform T S (t,η) is defined as:

[0099]

[0100] in, is the continuous wavelet transform of the signal; t, a and η represent time, scale and frequency respectively; is the instantaneous frequency of the signal, defined as:

[0101]

[0102] The present invention adopts the optimal basic wavelet as the base wavelet in the synchronous squeezing transformation framework. The optimal basic wavelet fits the seismic wavelet through the three-parameter wavelet to achieve the best match for the seismic wavelet. The three-parameter wavelet is defined in the time domain as:

[0103]

[0104] Wherein, σ, τ and β represent the modulation frequency, energy attenuation factor and energy delay factor respectively. Λ=(σ,τ,β) represents the set of σ, τ and β. In addition,

[0105]

[0106] Wavelet transform can transform one-dimensional time-domain seismic records into two-dimensional time-frequency domain, thereby extracting more information to reveal the changes in rock physical properties inside the sedimentary body and provide important phase markers for sedimentary phase research.

[0107] 6) The ridge extraction method is used to extract the instantaneous main frequency curve from the time-frequency spectrum of each seismic trace to obtain a two-dimensional instantaneous main frequency attribute profile.

[0108] 7) Use the instantaneous main frequency attributes to explain the lateral distribution of sedimentary cyclic units.

[0109] In order to verify the effectiveness of the proposed method, a Figure 2 (a) shows a sedimentary cycle model with two scales. From a large-scale perspective, the model includes one transgressive cycle unit and three regressive cycle units. From a small-scale perspective, the model includes two transgressive cycle units and five regressive cycle units. Generally speaking, large-scale sedimentary sequence interfaces correspond to strong seismic reflection amplitudes and have large reflection coefficient values; small-scale sequence interfaces correspond to weak seismic reflection amplitudes and have small reflection coefficient values. Therefore, the present invention is based on Figure 2 The sedimentary cycle geological model shown in (a) is designed Figure 2 (b) shows the reflection coefficient value of the geological model. The reflection coefficient value of the large-scale sedimentary cycle unit boundary is 0.08, and the reflection coefficient value of the small-scale sedimentary cycle unit boundary is 0.04. The reflection coefficient is convolved with a 30Hz main frequency zero phase Ricker wavelet to obtain Figure 2 (c) shows the synthetic seismic record. Figure 2 (c) The characteristics of the medium-strong reflection waveform can be used to approximately determine the location of the boundary of the large-scale cyclonic unit, but it is not accurate. Figure 2 (f) shows the time-frequency spectrum of the synthetic seismic record calculated by synchronous extrusion optimal basis wavelet transform; it is not difficult to see that the time-frequency spectrum contains two frequency change components at the same time, and these two components can reflect the instantaneous main frequency change trend caused by two scale sedimentary cycles. The low-frequency component here represents the thickness change trend of the large-scale sedimentary cycle, and the high-frequency component represents the thickness change trend of the small-scale sedimentary cycle. Directly using Figure 2 The time-frequency spectrum shown in (f) can approximately assist in the division of sedimentary cycle units of the two scales, but the overlap of the two scales can easily lead to division errors, especially at the location where the frequency characteristics of the two scales overlap, it is impossible to determine which scale of sedimentary cycle the frequency characteristics at that location are caused by.

[0110] In order to better utilize seismic data to assist in the division of sedimentary cycle units of different scales, the present invention further uses the variational mode decomposition method to decompose the seismic waveform characteristics of different scales contained in the seismic data. Then, the synchronous squeezed optimal basis wavelet transform is used to extract sedimentary cycle information of different scales from the intrinsic mode function. The method first uses Fourier transform to calculate the amplitude spectrum of the synthetic seismic trace, such as Figure 3 (a) is shown by the dotted line. Then, the compression mapping method is used to extract the seismic wavelet amplitude spectrum from the synthetic seismic trace amplitude spectrum, and the two are divided to obtain the reflection coefficient amplitude spectrum, as shown in Figure 3 (a) is shown by the solid line. Further, the scale space representation is used to calculate the smoothed reflection coefficient spectrum, as Figure 3(a) dotted line. It can be seen from the smoothed reflection coefficient spectrum that it contains two extreme values, and the frequencies corresponding to the extreme value positions are 21Hz and 52Hz respectively. Therefore, the number of intrinsic mode functions is set to 2, and the center frequencies of eigenmode function 1 and eigenmode function 2 are fixed to 21Hz and 52Hz respectively. With the above variational mode decomposition parameters, the variational mode decomposition method is used to decompose the synthetic seismic record trace to obtain Figure 2 (d) and (e) show the eigenmode function 1 and eigenmode function 2. It is not difficult to see that the eigenmode function 1 reflects the strong amplitude and low frequency seismic reflection waveform characteristics corresponding to the large-scale sedimentary cycle, and the eigenmode function 2 reflects the weak amplitude and high frequency seismic waveform characteristics corresponding to the small-scale sedimentary cycle. Figure 3 The dotted line and dashed line in (b) are the amplitude spectra of eigenmode function 1 and eigenmode function 2, respectively. It can be seen that eigenmode function 1 and eigenmode function 2 reflect the low-frequency large-scale and high-frequency small-scale characteristics of seismic data, respectively. In order to better characterize the sedimentary cycle, the synchronous extrusion optimal basis wavelet transform is further used to calculate the eigenmode function 1 and eigenmode function 2 to obtain Figure 2 (g) and (i) are the time-frequency spectra. It is not difficult to see that Figure 2 (g) clearly depicts the instantaneous main frequency variation characteristics caused by the thickness change of large-scale sedimentary cycles. In order to more clearly characterize the sedimentary cycles and assist in the division of sedimentary cycle units, the time-frequency spectrum ridge extraction method is used to obtain Figure 2 (h) shows the instantaneous main frequency curve and smoothes it. According to the change trend of the smoothed instantaneous main frequency curve, one transgressive and three regressive sedimentary cycle units in the large-scale sedimentary cycle can be accurately divided. This result is consistent with Figure 2 The geological model shown in (a) has a good match. Figure 2 (i) Clearly depicts the instantaneous main frequency characteristics caused by the thickness change of small-scale sedimentary cycles. The time-frequency spectrum ridge extraction method is used to calculate Figure 2 (j) The instantaneous main frequency curve is smoothed. According to the change trend of the smoothed instantaneous main frequency curve, two transgressive and five regressive sedimentary cycle units in the small-scale sedimentary cycle can be accurately divided. This result is consistent with Figure 2 The geological model shown in (a) has a good match.

[0111] According to the model example, the variational mode decomposition method can well separate the seismic reflection waveform characteristics caused by sedimentary cycles of different scales. The time-frequency spectrum calculated by using synchronous extrusion optimal basis wavelet transform for different intrinsic mode functions based on variational mode decomposition only contains the instantaneous frequency characteristics of sedimentary cycles of one scale, which is very helpful to assist the subsequent division of multi-scale sedimentary cycle units.

[0112] Figure 4 (a) shows an actual 2D seismic profile used to verify the workflow proposed in the present invention. The data comes from a basin rich in oil and gas, a craton basin with an area of ​​250,000 square kilometers. Figure 4 The T3 layer in (a) is the main source rock formed during the maximum flooding period, corresponding to a set of oil shales with a thickness of more than 100 meters. Between T2 and T3 are thin interbeds of sandstone and mudstone formed in the alternating stages of deep lake shales and turbidity current sediments, which are the main reservoirs. It should be noted that the thin interbeds between T2 and T3 belong to the quaternary sequence within two tertiary sequences (i.e., the high-water system domain and the low-water system domain). The purpose of the present invention is to use the proposed method to assist in the characterization and division of sedimentary cycles of two scale sequences (i.e., tertiary and quaternary sequences).

[0113] The smoothed reflection coefficient spectrum calculated by the technical process proposed by the present invention contains two extreme points, and the corresponding extreme frequencies are 22 Hz and 48 Hz respectively. The two extreme frequencies represent the frequency characteristics of the third-order sequence and the fourth-order sequence on the seismic section. Therefore, the present invention sets the number of intrinsic mode functions to 2, and the center frequencies of intrinsic mode function 1 and intrinsic mode function 2 are set to 22 Hz and 48 Hz respectively. Figure 4 (c) and (e) are the intrinsic mode function 1 and intrinsic mode function 2 respectively obtained by using the empirical mode decomposition method. Figure 4 (c) It can better characterize the low-frequency seismic reflection characteristics of the third-order sequence. Figure 4 (e) It can better characterize the high-frequency seismic reflection characteristics of the fourth-order sequence.

[0114] Furthermore, the synchronous squeezed optimal basis wavelet transform and ridge extraction method are used to Figure 4 The seismic profiles shown in (a), (c) and (e) are calculated to obtain Figure 4 (b), (d) and (f) show the instantaneous main frequency property profiles. Figure 4 Compared with (d) and (f), Figure 4 (b) has poor lateral continuity. This is because the original seismic data contains seismic reflection characteristics of two scale sequences. The time-frequency spectrum directly calculated using the synchronous squeezing optimal basis wavelet transform method has two components. It is difficult to accurately extract the instantaneous main frequency characteristics of the time-frequency spectrum using the ridge extraction method based on this. Figure 4 The instantaneous main frequency attribute profiles shown in (d) and (f) have good lateral continuity, and respectively characterize large-scale (third-order sequence) sedimentary cycle units and small-scale (fourth-order sequence) sedimentary cycle units. At five well points, the present invention divides the sedimentary cycle units according to the instantaneous main frequency attributes, and the results of the division of the two scale sedimentary cycle units have a good matching relationship with the formation thickness variation law characterized by the logging impedance curve.

[0115] Figure 5 The figure shows the characterization results of two scale sedimentary cycles at Well 2. The spectrum of the original wellside seismic trace directly calculated by synchronous extrusion optimal basis wavelet transform contains two frequency components, which respectively characterize the large-scale sedimentary cycle (third-order sequence) and the small-scale sedimentary cycle (fourth-order sequence). Due to the existence of multiple frequency components, the ridge extraction method cannot accurately calculate the instantaneous main frequency point from the time-frequency spectrum of the original seismic trace, which is also Figure 4 (b) is the reason for the poor lateral continuity of the instantaneous main frequency profile shown. This time-frequency spectrum with multi-scale information mixed together is not conducive to the interpretation of multi-scale sedimentary cycles. The intrinsic mode function 1 and the intrinsic mode function 2 obtained by the proposed VMD decomposition process respectively reflect the seismic reflection waveform characteristics of the two scale sedimentary sequences. The time-frequency spectrum calculated by the synchronous extrusion optimal basis wavelet transform for the intrinsic mode function 1 and the intrinsic mode function 2 respectively does not contain other frequency components and can effectively characterize the changing laws of large-scale and small-scale sedimentary cycles. The instantaneous main frequency curve extracted by the ridge extraction method can effectively characterize the changing trend of the stratum thickness in sedimentary cycles of different scales, and the sedimentary cycle units divided accordingly can well match the logging lithology interpretation results. The above results further verify the effectiveness of the method proposed in the present invention in dividing multi-scale sedimentary cycles.

[0116] On the other hand, the present invention can also provide a multi-scale sedimentary cycle phase space fine characterization system, including a first calculation module, a second calculation module, a mapping module, a time-frequency spectrum calculation module and a characterization module;

[0117] The first calculation module is used to calculate the average Fourier amplitude spectrum and the average reflection coefficient amplitude spectrum of the two-dimensional seismic data;

[0118] The second calculation module calculates the smoothed reflection coefficient amplitude spectrum using a scale space representation method, and according to the relationship between the reflection coefficient amplitude spectrum and the empirical mode, sets the number of extreme values ​​of the smoothed reflection coefficient spectrum to the number of intrinsic mode function components, and the extreme value frequency is the center frequency of each intrinsic mode function component;

[0119] The mapping module uses the variational mode decomposition method to decompose the input two-dimensional seismic data to obtain the intrinsic mode function components that can characterize the seismic reflection characteristics of geological structures at different scales;

[0120] The time-frequency spectrum calculation module uses synchronous extrusion optimal basis wavelet transform to calculate each seismic trace of each two-dimensional intrinsic mode function component to obtain the time-frequency spectrum of each seismic trace.

[0121] The characterization module uses the ridge extraction method to extract the instantaneous main frequency curve from the time-frequency spectrum of each seismic trace to obtain a two-dimensional instantaneous main frequency attribute profile. The instantaneous main frequency attribute profile is used to explain the lateral distribution law of the sedimentary cycle unit, that is, a detailed characterization of the multi-scale sedimentary cycle phase space.

[0122] In summary, the present invention provides a method for fine characterization of the phase space of multi-scale sedimentary cycles. The method takes into account that seismic reflection data simultaneously contains seismic reflection information of sedimentary sequences of different scales. It overcomes the problem that traditional methods cannot effectively separate and characterize the multi-scale sedimentary cycle information contained in seismic data. This technology can effectively separate the seismic reflection wave characteristics of different scales from seismic data by improving the variational mode decomposition method. On this basis, the time-frequency spectrum calculated by the synchronous extrusion optimal basis wavelet transform can accurately characterize the changing laws of sedimentary cycles of different scales. The instantaneous main frequency attribute proposed in combination with the ridge extraction method can effectively characterize the lateral distribution range of sedimentary cycles of different scales.

[0123] On the other hand, the present invention provides a computer-readable storage medium, in which a computer program is stored. When the computer program is executed by a processor, the multi-scale sedimentary cycle phase space fine characterization method described in the present invention can be implemented.

[0124] The computer device may be a laptop computer, a desktop computer or a workstation.

[0125] The present invention can also provide a computer device, including a processor and a memory, the memory is used to store a computer executable program, the processor reads the computer executable program from the memory and executes it, and when the processor executes the computer executable program, the multi-scale sedimentary cycle phase space fine characterization method described in the present invention can be implemented.

[0126] The processor may be a central processing unit (CPU), a digital signal processor (DSP), an application specific integrated circuit (ASIC), or an off-the-shelf field programmable gate array (FPGA).

[0127] The memory of the present invention may be an internal storage unit of a laptop computer, a desktop computer or a workstation, such as a memory or a hard disk; or an external storage unit, such as a mobile hard disk or a flash memory card.

[0128] Computer-readable storage media may include computer storage media and communication media. Computer storage media include volatile and non-volatile, removable and non-removable media implemented in any method or technology for storing information such as computer-readable instructions, data structures, program modules or other data. Computer-readable storage media may include: read-only memory (ROM), random access memory (RAM), solid-state drives (SSD) or optical disks, etc. Among them, random access memory may include resistance random access memory (ReRAM) and dynamic random access memory (DRAM)

[0129] In summary, the present invention discloses a method for fine characterization of multi-scale sedimentary cycle phase space. The method takes into account that seismic reflection data simultaneously contains seismic reflection information of sedimentary sequences of different scales. The present invention can effectively separate seismic reflection wave characteristics of different scales from seismic data by improving the variational modal decomposition method. On this basis, the time-frequency spectrum calculated by synchronous extrusion optimal basis wavelet transform can accurately characterize the changing law of sedimentary cycles of different scales. The instantaneous main frequency attribute proposed in combination with the ridge extraction method can effectively characterize the lateral distribution range of sedimentary cycles of different scales, overcoming the problem that traditional methods cannot effectively separate and characterize the multi-scale sedimentary cycle information contained in seismic data. The key parameters in the processing process of the present invention are all extracted from seismic data, which is simple to operate and fast in calculation. It is suitable for processing massive seismic data, which provides an important tool for subsequent sequence stratigraphic research.

[0130] The above contents are only for explaining the technical idea of ​​the present invention and cannot be used to limit the protection scope of the present invention. Any changes made on the basis of the technical solution in accordance with the technical idea proposed by the present invention shall fall within the protection scope of the claims of the present invention.

Claims

1. A multi-scale sedimentary cycle phase space fine characterization method, characterized in that: The following steps are involved: Calculate the average Fourier amplitude spectrum of 2D seismic data; Calculate the average reflection coefficient amplitude spectrum of 2D seismic data; The scale space representation method is used to calculate the smoothed reflection coefficient amplitude spectrum. According to the relationship between the reflection coefficient amplitude spectrum and the empirical mode, the number of extreme values ​​of the smoothed reflection coefficient spectrum is set to the number of intrinsic mode function components, and the extreme value frequency is set to the center frequency of each intrinsic mode function component. The variational mode decomposition method is used to decompose the input two-dimensional seismic data to obtain the intrinsic mode function components that can characterize the seismic reflection characteristics of geological structures of different scales; The synchronous squeezed optimal basis wavelet transform is used to calculate each seismic trace of each two-dimensional intrinsic mode function component to obtain the time-frequency spectrum of each seismic trace. The ridge extraction method is used to extract the instantaneous main frequency curve from the time-frequency spectrum of each seismic trace to obtain a two-dimensional instantaneous main frequency attribute profile. The instantaneous main frequency attribute profile is used to interpret the lateral distribution of sedimentary cycle units, that is, the detailed characterization of the multi-scale sedimentary cycle phase space.

2. The multi-scale sedimentary cycle phase space fine characterization method according to claim 1 is characterized in that: Calculating the average reflection coefficient amplitude spectrum of two-dimensional seismic data includes extracting the seismic wavelet amplitude spectrum from the average Fourier amplitude spectrum of the seismic data using a compression mapping operator, and obtaining the reflection coefficient amplitude spectrum by dividing the average Fourier amplitude spectrum of the seismic data by the seismic wavelet amplitude spectrum.

3. The multi-scale sedimentary cycle phase space fine characterization method according to claim 2 is characterized in that: The compression mapping operator is used to extract the seismic wavelet amplitude spectrum from the seismic data amplitude spectrum; The compression mapping operator for calculating the seismic sub-spectrum is introduced. Assume that the seismic sub-spectrum is within the specified frequency band [f l ,f h ] appears as a smooth and unimodal function: q∈(0,1), 0≤F q (|W|; f)≤1, let c q >0, α, β>1, the compression mapping operator P is Set the initial value of the iteration to: By iterating the compression mapping algorithm |W| k =P(|W| k-1 ) so that it converges to a fixed point |W| * , the fixed point corresponds to the estimated seismic sub-spectrum (|W| * ) 1 / q ; The reflection coefficient amplitude spectrum |R is obtained by dividing the seismic data amplitude spectrum |S(f)| by the seismic wavelet amplitude spectrum |W(f| k (f)|: Among them, ε≠0, preventing |W(f)|=0.

4. The multi-scale sedimentary cycle phase space fine characterization method according to claim 1, characterized in that: According to the relationship between the reflection coefficient spectrum and the empirical mode, the number of extreme values ​​of the smoothed reflection coefficient spectrum is set to the number of intrinsic mode function components, and the extreme value frequency is the center frequency of each intrinsic mode function component. The relationship between the reflection coefficient amplitude spectrum and geological structures of different scales is obtained based on theoretical derivation, as shown in the following formula: The seismic reflection of structures of different scales is defined as the empirical mode function, where |R k (f) | represents the amplitude spectrum of the reflection coefficient corresponding to the geological structure at the kth scale, c kj represents the reflection coefficient value of the jth sampling point in the kth empirical mode function, (t km -t kn ) represents the time difference between two reflection coefficient points m and n at the kth scale. It is defined that a smaller absolute value k represents a large-scale structure, and a larger value k represents a small-scale structure. Equation (1) shows that |R k The frequency of the first extreme point of (f)| and (t km -t kn ) is inversely proportional to.

5. The multi-scale sedimentary cycle phase space fine characterization method according to claim 1 is characterized in that: The variational mode decomposition method is used to decompose the input two-dimensional seismic data, and the intrinsic mode function components that can characterize the seismic reflection characteristics of geological structures of different scales are obtained, including: The signal separation and extraction is achieved by decomposing the signal into eigenmode functions with finite bandwidth, which is achieved by solving the following optimization problem: Among them, s(t) represents single-channel seismic data, u m (t) represents the mth intrinsic mode function, ω m represents the center frequency of the mth eigenmode function, represents the Hilbert transform of the mth intrinsic mode function. Further, the Lagrange multiplier method shown in the following formula is used to solve the above optimization problem: The value of the center frequency is determined by the frequency of the extreme points of the smoothed reflection coefficient amplitude spectrum.

6. The multi-scale sedimentary cycle phase space fine characterization method according to claim 1, characterized in that: The synchronous squeezed optimal basis wavelet transform is used to calculate each seismic trace of each two-dimensional intrinsic mode function component, and the time-frequency spectrum of each seismic trace is obtained, including: The optimal basic wavelet is used as the base wavelet in the synchronous squeezing transform framework. The optimal basic wavelet is obtained by fitting the seismic wavelet with a three-parameter wavelet. The three-parameter wavelet is defined in the time domain as: Among them, σ, τ and β represent the modulation frequency, energy attenuation factor and energy delay factor respectively, Λ=(σ,τ,β) represents the set of σ, τ and β, 7. A multi-scale sedimentary cycle phase space fine characterization system, characterized in that: It includes a first calculation module, a second calculation module, a mapping module, a time-frequency spectrum calculation module and a characterization module; The first calculation module is used to calculate the average Fourier amplitude spectrum and the average reflection coefficient amplitude spectrum of the two-dimensional seismic data; The second calculation module calculates the smoothed reflection coefficient amplitude spectrum using a scale space representation method, and according to the relationship between the reflection coefficient amplitude spectrum and the empirical mode, sets the number of extreme values ​​of the smoothed reflection coefficient spectrum to the number of intrinsic mode function components, and the extreme value frequency is the center frequency of each intrinsic mode function component; The mapping module uses the variational mode decomposition method to decompose the input two-dimensional seismic data to obtain the intrinsic mode function components that can characterize the seismic reflection characteristics of geological structures at different scales; The time-frequency spectrum calculation module uses synchronous extrusion optimal basis wavelet transform to calculate each seismic trace of each two-dimensional intrinsic mode function component to obtain the time-frequency spectrum of each seismic trace. The characterization module uses the ridge extraction method to extract the instantaneous main frequency curve from the time-frequency spectrum of each seismic trace to obtain a two-dimensional instantaneous main frequency attribute profile. The instantaneous main frequency attribute profile is used to explain the lateral distribution law of the sedimentary cycle unit, that is, a detailed characterization of the multi-scale sedimentary cycle phase space.

8. A computer device, characterized in that: It includes a processor and a memory, the memory is used to store a computer executable program, the processor reads part or all of the computer executable program from the memory and executes it, and when the processor executes part or all of the computer executable program, it can implement the multi-scale sedimentary cycle phase space fine characterization method described in any one of claims 1-6.

9. A computer-readable storage medium, characterized in that: A computer program is stored in a computer-readable storage medium, and when the computer program is executed by a processor, the method for fine characterization of the multi-scale sedimentary cycle phase space according to any one of claims 1 to 6 can be implemented.

10. An application of a multi-scale sedimentary cycle phase space fine characterization method, characterized in that: Used to identify and interpret sedimentary cycle information extracted from seismic data.

Images