CN111242198B - Gather tuning optimization dimension increasing method based on sparse representation - Google Patents
Gather tuning optimization dimension increasing method based on sparse representation Download PDFInfo
- Publication number
- CN111242198B CN111242198B CN202010013149.9A CN202010013149A CN111242198B CN 111242198 B CN111242198 B CN 111242198B CN 202010013149 A CN202010013149 A CN 202010013149A CN 111242198 B CN111242198 B CN 111242198B
- Authority
- CN
- China
- Prior art keywords
- wavelet
- frequency
- signal
- matched
- time
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 46
- 238000005457 optimization Methods 0.000 title claims abstract description 13
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 20
- 238000004458 analytical method Methods 0.000 claims description 20
- 238000004422 calculation algorithm Methods 0.000 claims description 18
- 230000008569 process Effects 0.000 claims description 11
- 239000011159 matrix material Substances 0.000 claims description 10
- 238000013016 damping Methods 0.000 claims description 7
- 230000008859 change Effects 0.000 claims description 6
- 230000015572 biosynthetic process Effects 0.000 claims description 4
- 238000004364 calculation method Methods 0.000 claims description 4
- 238000005070 sampling Methods 0.000 claims description 4
- 238000003786 synthesis reaction Methods 0.000 claims description 4
- 238000012897 Levenberg–Marquardt algorithm Methods 0.000 claims description 3
- 238000012216 screening Methods 0.000 claims description 3
- 239000012141 concentrate Substances 0.000 claims description 2
- 230000000694 effects Effects 0.000 abstract description 10
- 238000010586 diagram Methods 0.000 description 9
- 238000012512 characterization method Methods 0.000 description 7
- 230000002776 aggregation Effects 0.000 description 2
- 238000004220 aggregation Methods 0.000 description 2
- 238000005530 etching Methods 0.000 description 2
- 238000001914 filtration Methods 0.000 description 2
- 238000003384 imaging method Methods 0.000 description 2
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 description 2
- 230000000737 periodic effect Effects 0.000 description 2
- 230000009466 transformation Effects 0.000 description 2
- 230000001052 transient effect Effects 0.000 description 2
- 238000010521 absorption reaction Methods 0.000 description 1
- 101150043283 ccdA gene Proteins 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000008021 deposition Effects 0.000 description 1
- 239000000284 extract Substances 0.000 description 1
- 239000003345 natural gas Substances 0.000 description 1
- 239000003209 petroleum derivative Substances 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 230000002441 reversible effect Effects 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/21—Design or setup of recognition systems or techniques; Extraction of features in feature space; Blind source separation
- G06F18/213—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods
- G06F18/2136—Feature extraction, e.g. by transforming the feature space; Summarisation; Mappings, e.g. subspace methods based on sparsity criteria, e.g. with an overcomplete basis
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/28—Determining representative reference patterns, e.g. by averaging or distorting; Generating dictionaries
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
- Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
- Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation
Landscapes
- Engineering & Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Artificial Intelligence (AREA)
- Bioinformatics & Computational Biology (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Evolutionary Biology (AREA)
- Evolutionary Computation (AREA)
- Physics & Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention provides a gather tuning optimization dimension-increasing method based on sparse representation, which comprises the following steps: step 1, constructing a time-frequency domain wavelet dictionary; step 2, performing Hilbert transform on the gather; step 3, determining matched wavelet control parameters; step 4, carrying out wavelet decomposition based on matching tracking; and 5, performing high-resolution time-frequency decomposition based on the matched wavelets. The gather tuning optimization dimension increasing method based on sparse representation avoids the seismic wave interference phenomenon generated by the superposition effect, eliminates the distortion influence caused by frequency aliasing, and can effectively improve the later inversion depiction precision.
Description
Technical Field
The invention relates to the technical field of geophysical prospecting of petroleum and natural gas resources, in particular to a gather tuning optimization dimension increasing method based on sparse representation.
Background
The land phase deposition reservoir layer in China has the characteristics of relatively thin thickness and large space variation, and can cause the propagation of seismic waves in the reservoir layer to generate the interference effect of thin interbed layers. In order to improve the identification accuracy of thin interbed reservoirs, the most important solution is to widen the frequency band of seismic data and improve the main frequency of the data. For imaging post-stack seismic data, the current common processing means include methods of expanding frequency bands and improving dominant frequency such as inverse Q filtering, frequency expansion and frequency division. The anti-Q filtering aims at the problems of high-frequency energy absorption attenuation and phase distortion caused by the increase of the propagation of seismic waves along with the travel time, and compensates the amplitude and the phase of the seismic waves, so that the main frequency of the seismic data is improved; the frequency expansion method is to compress the wavelength of the seismic signal, thereby achieving the effects of widening the frequency band of the seismic data and improving the main frequency of the seismic data; the frequency division method extracts high-frequency information through frequency spectrum decomposition to achieve the effect of improving the frequency of the seismic data. The methods play a certain role in improving the recognition precision of the thin interbed, but have yet to be studied in terms of fidelity of data.
Sparsity has become one of the important concepts in a wide range of signal processing fields, both theoretically and practically, and is always a very attractive signal feature. In the field of seismic signal processing, sparsity can be used as priori information to facilitate the solving of linear inversion problems such as seismic data denoising, deconvolution, seismic data reconstruction and the like. Since most signals in nature are usually not sparse, sparse representation is required. Sparse representation is the representation of a signal in the form of a linear combination of very few elements in a dictionary. As a feature of the signal, the sparsity of the seismic reflection coefficient sequence and the sparsity of the seismic data in the overcomplete dictionary can be reduced in dimension on the premise of guaranteeing effective information, and the solution of the seismic inversion problem is effectively promoted. From a frequency domain perspective, the data space dimension is smaller than the model space dimension, so the seismic inversion has multiple solutions. Under specific conditions, the solution of the sparse constraint inverse problem has unique and stable characteristics, and the optimal solution can be automatically selected.
The existing pre-stack inversion gather processing method mainly has the following problems: (1) because the prestack gather is a high-dimensional data volume, the data volume is large, the redundancy is high, and the interpretation efficiency is low one by one, the method adopted is a strategy of partial angle superposition, but the superposition can generate interference effect to reduce the seismic resolution; (2) the influence of frequency on the seismic attribute is not considered in the traditional pre-stack inversion, and meanwhile, the traditional time-frequency characterization method is unfavorable for fine characterization of thin interbings because the time-frequency resolution is not high. Therefore, the invention discloses a novel gather tuning optimization dimension increasing method based on sparse representation, and solves the technical problems.
Disclosure of Invention
The invention aims to provide a gather tuning optimization dimension increasing method based on sparse representation, which avoids the seismic wave interference phenomenon caused by superposition effect, eliminates the distortion influence caused by frequency aliasing and can effectively improve the later-stage reverse-modeling precision.
The aim of the invention can be achieved by the following technical measures: the gather tuning optimization dimension increasing method based on sparse representation comprises the following steps: step 1, constructing a time-frequency domain wavelet dictionary; step 2, performing Hilbert transform on the gather; step 3, determining matched wavelet control parameters; step 4, carrying out wavelet decomposition based on matching tracking; and 5, performing high-resolution time-frequency decomposition based on the matched wavelets.
The aim of the invention can be achieved by the following technical measures:
in step 1, a Ricker wavelet dictionary is optionally built when analyzing the seismic signals:
wherein ,ωR (t,f i ) Representing a dictionary of Ricker wavelets, t representing the center time of the Ricker wavelets, f i Representing the dominant frequency of the ith Ricker wavelet in the wavelet dictionary.
The Ricker wavelet has three control parameters, the convergence rate of sidelobe energy is variable, the delay time is controllable, the wavelet base has rich waveforms, and the matching with the seismic wavelet is flexible.
In step 2, for a non-stationary signal, its instantaneous frequency is:
wherein ,
is the instantaneous phase, t is the time;
each signal x (t) Hilbert within the gather is transformed into:
where P is the Cauchy principal value and τ represents the time shift.
The analytical signal for x (t) is thus obtained as:
in step 3, three control parameters of the matched wavelet are shifted by τ, frequency factor f and phase factor
The generated space discretization, any discretized point +.>
A matching wavelet can be generated according to the formula (1) to form a matching wavelet base phi= { omega δ (t)} δ∈Γ Γ is expressed as->
Is a collection of (3); determining a signal sampling rate and a signal length through a time sequence by using the obtained analytic signal; determining the instantaneous envelope and instantaneous phase of the signal by using the analysis signal, τ k The time t corresponding to the local maximum peak value of the instantaneous envelope can be roughly estimated; scanning frequency parameters of atoms in wavelet base in complex domain, determining several tau for a certain iteration m Obtaining the frequency a priori value f of the signal at these positions m Assuming that the matched wavelet is zero-phase wavelet, τ is used m and fm Determining a number of a priori information points (tau m ,f m 0); selecting a frequency scanning range by taking the prior information points as the center, and determining the peak frequency f of the matched wavelet k 。
In step 3, after the center time and the dominant frequency of the best matching wavelet omega are determined, according to the residual energy minimum principle, the amplitude and phase information of the matching wavelet are removed by utilizing a damping minimum flat method;
assuming that a seismic signal x (t) with limited bandwidth exists, iterating k-1 times through a matching pursuit algorithm, wherein a residual signal is R k-1 x; when iterating to the kth time, M matched wavelets omega are generated in total δm (t), m=1, 2, …, M; for residual signal R k- 1 x performs Hilbert transform structure analysis signal:
X＝R k-1 x+iR k-l x (5)
and constructing complex matched wavelets at the kth iteration in a similar way:
W δm (t)＝ω δm (t)+iω δm (t) (6)
residual signal R generated after k iterations k The energy of x is:
where X is the N non-column vector of the seismic signal, A represents the complex amplitude of the matched wavelet:
a m the amplitudes of the wavelets are matched for the real number domain,
to match the phase of the wavelet.
G δm (t)＝[ω δ (t,f 1 ),ω δ (t,f 2 ),ω δ (t,f 3 ),…,ω δ (t,f m )]For Nxm order matching wavelet momentsAn array, each column corresponding to a matching wavelet.
To match the seismic signal x (t) to the set of matched wavelets ωδ m The relative error between (t) is minimum, and the method is calculated by adopting a damping least squares method:
A＝[W T W+εI] -1 W T X (9)
wherein W is an N multiplied by M order complex matching wavelet matrix obtained after k iterations, epsilon is a damping factor, and I is an M multiplied by M order identity matrix; from equations (8) and (9), a control parameter a can be determined that matches wavelet amplitude and phase k And
in step 4, the matched wavelet is formed by time shifting, modulating and phase varying the basic wavelet ω (t); if the basic wavelet omega (t) is one basic wavelet meeting the condition in the matched wavelet base, making
As a control parameter of the matching wavelet, the matching wavelet is expressed as:
where τ represents the time shift, f represents the frequency factor,
expressed as a phase factor; the position of the matched wavelet is fixed by using time shift, the energy of the basic wavelet is concentrated near a few time-frequency points by using frequency factors, and the phase factors are used for controlling the shape change of the matched wavelet.
In step 4, after the overcomplete matching wavelet base Φ is constructed, although the matching wavelets in Φ are not all orthogonal, it is still required to satisfy ω i || 2 =1; selecting ω from the matched wavelet base Φ which most closely matches the given signal x (t) 1 The best matching requirement is their inner product<x,ω 1 >To match wavelet baseThe largest product of all atoms and x (t):
<x,ω 1 >＞<x,ω i >i≠1 ω i ∈Φ (11)
therefore, after one iteration, the signal x (t) is decomposed into:
x(t)＝<x,ω 1 >ω 1 +R 1 x (12)
R 1 x is the residual signal after the first iteration; r is R 1 x and omega 1 Orthogonalization, therefore:
<x,ω 1 >the larger means ω 1 Is in Hilbert space where the signal is located and
the direction is closest to the unit vector in the matched wavelet base; for residual signal R 1 x repeats the same matching process: screening from phi and R 1 x is the second atom which is most matched, the process is repeated continuously until the kth time (k is more than or equal to 0):
R k x＝<x,ω k-1 >ω k-1 +R k-1 x (14)
after k iterative decomposition calculations, the original signal x (t) is approximated by a synthesis of k matched wavelets, as residual R 1 And x is small enough, namely, when the sparse representation of the signal is good enough to approximate the signal or the residual signal meets a set threshold value, the iteration is stopped.
In step 5, according to the high-precision matching wavelets obtained in steps 3 and 4, using the independence of each matching wavelet in the wavelet set, selecting a needed time-frequency decomposition method, and respectively performing time-frequency domain expansion on the matching wavelets; obtaining a time-frequency spectrogram of the whole signal through linear superposition; the loop is repeated for each track in the gather to finally obtain the multi-resolution de-aliased tuning data volume.
The gather tuning optimization dimension increasing method based on sparse representation is mainly used for high-resolution prestack inversion gather processing, and further obtains a high-precision multi-resolution anti-aliasing tuning data body through a stable and reliable sparse representation method. The method fully utilizes rich seismic information contained in the pre-stack channel set, considers that the seismic attribute of the channel set is influenced by frequency while changing along with the change of offset, and the superposition of signals of the large-angle channel set can cause wave interference, influence the quality of seismic data and the signals can be influenced by frequency aliasing in the process of acquisition. Each track in the gather is sequentially unfolded, interference effects caused by superposition are avoided, and high-precision time-frequency unfolding is performed on each track by utilizing a matching pursuit algorithm, so that time-frequency resolution is improved, the influence of frequency aliasing is avoided, and the precision of thin reservoir characterization is improved. In the process of processing the gathers before the inversion work is carried out, the invention considers that the seismic attribute of the gathers is influenced by frequency while being changed along with the change of the offset. In addition, superposition of large-angle gather signals can cause wave interference, influence the quality of seismic data, and signals can be influenced by frequency aliasing in the acquisition process. Therefore, the method sequentially expands each track in the track set, avoids interference effect caused by superposition, utilizes a matching pursuit algorithm to carry out high-precision time-frequency expansion on each track, improves time-frequency resolution, avoids the influence of frequency aliasing, and improves the drawing precision of thin reservoir etching.
Drawings
FIG. 1 is a diagram of a Rick wavelet dictionary according to one embodiment of the present invention;
FIG. 2 is a diagram of a Hilbert transform signal according to an embodiment of the present invention;
FIG. 3 is a diagram illustrating the instantaneous envelope of a signal and the instantaneous frequency of the signal in accordance with one embodiment of the present invention;
FIG. 4 is a schematic diagram of a wavelet set decomposed by a matching pursuit algorithm in accordance with an embodiment of the present invention;
FIG. 5 is an exploded reconstruction contrast diagram of a matching pursuit algorithm according to an embodiment of the present invention;
FIG. 6 is a time-frequency resolution comparison chart of a conventional time-frequency analysis method and a time-frequency analysis method based on a matching pursuit algorithm according to an embodiment of the present invention;
FIG. 7 is a schematic diagram of a time-frequency characterization method based on a matching pursuit algorithm for distinguishing thin reservoirs in accordance with an embodiment of the present invention;
FIG. 8 is a diagram of a volume of tuning data for multi-resolution de-aliasing for a gather in an embodiment of the present invention;
FIG. 9 is a flowchart of an embodiment of the sparse representation based gather tuning optimization dimension increase method of the present invention.
Detailed Description
The foregoing and other objects, features and advantages of the invention will be apparent from the following more particular description of preferred embodiments, as illustrated in the accompanying drawings.
As shown in fig. 9, fig. 9 is a flowchart of the sparse representation-based gather tuning optimization dimension-increasing method of the present invention.
And step 101, constructing a time-frequency domain wavelet dictionary.
In order to enable the sparse representation result to better meet the requirement of analysis signals, the overcomplete atomic dictionary needs to meet certain conditions: the wavelet dictionary constructed should match as closely as possible the inherent structure or behavior of the signal; atoms in the dictionary should have a better time-frequency aggregation in order to better accommodate non-stationary signals. The Ricker wavelet has the characteristics of simple waveform, faster convergence and short delay time, and has higher time-frequency resolution in thin interbed analysis, so that the Ricker wavelet dictionary is selectively constructed when the seismic signals are analyzed.
In addition, the Ricker wavelet has three control parameters, the convergence rate of sidelobe energy is variable, the delay time is controllable, the wavelet base has rich waveforms, and the matching with the seismic wavelet is flexible.
The frequency domain of a signal is an efficient way to describe a signal, but it belongs to a global transformation and is completely non-localized in the time domain. In the physical and engineering fields, the instantaneous frequency is the most intuitive description of the phenomenon of periodic variation of complex signals. For a non-stationary signal, its instantaneous frequency is:
wherein ,
is the instantaneous phase. In order to find these transient properties, an analytical signal needs to be introduced.
The essence of resolving the signal is to convert a narrowband real signal into a complex signal, and the Hilbert transform can generate a unique complex signal for the narrowband real signal. Each signal x (t) Hilbert within the gather is transformed into:
wherein P is the Cauchy principal value. The analytical signal for x (t) is thus obtained as:
The three control parameters (time shift tau, frequency factor f, phase factor) of the matched wavelet
) The generated space discretization, any discretized point +.>
A matching wavelet can be generated according to the formula (1) to form a matching wavelet base phi= { omega δ (t)} δ∈Γ Γ is expressed as->
Is a collection of (a) to (b).
And determining the signal sampling rate and the signal length through a time sequence by using the obtained analysis signal. Determining the instantaneous envelope and instantaneous phase of the signal by using the analysis signal, τ k Can be roughly estimated as the corresponding time t at the instantaneous envelope local maximum peak. Scanning frequency parameters of atoms in wavelet base in complex domain, and determining several tau for a certain iteration m Obtaining the frequency a priori value f of the signal at these positions m Assuming that the matched wavelet is zero-phase wavelet, τ is used m and fm Determining a number of a priori information points (tau m ,f m ,0). Selecting a frequency scanning range by taking the prior information points as the center, and determining the peak frequency f of the matched wavelet k 。
Having determined the center time and dominant frequency of the best matching wavelet ω, the amplitude and phase information of these matching wavelets is next determined. And according to the principle of minimum residual energy, using a damping least squares method to calculate, and removing the amplitude and phase information of the matched wavelets.
Assuming that a seismic signal x (t) with limited bandwidth exists, iterating k-1 times through a matching pursuit algorithm, wherein a residual signal is R k-1 x. When iterating to the kth time, M matched wavelets omega are generated in total δm (t), m=1, 2, …, M. For residual signal R k- 1 x performs Hilbert transform structure analysis signal:
X＝R k-1 x+iR k-1 x (5)
and constructing complex matched wavelets at the kth iteration in a similar way:
W δm (t)＝ω δm (t)+iω δm (t) (6)
residual signals generated after k iterationsR k The energy of x is:
where A represents the complex amplitude of the matched wavelet:
a m the amplitudes of the wavelets are matched for the real number domain,
to match the phase of the wavelet.
G δm (t)＝[ω δ (t,f 1 ),ω δ (t,f 2 ),ω δ (t,f 3 ),…,ω δ (t,f m )]For an Nxm order matching wavelet matrix, each column corresponds to a matching wavelet.
To match the seismic signal x (t) to the set of matched wavelets ωδ m The relative error between (t) is the smallest (i.e. the residual signal energy Emin), calculated using the damped least squares method:
A＝[W T W+εI] -1 W T X (9)
where W is an n×m-order complex matched wavelet matrix obtained after k iterations, epsilon is a damping factor (usually 0.01), and I is an m×m-order identity matrix. From equations (8) and (9), a control parameter a can be determined that matches the amplitude and phase of the sub-wave k And
The overcomplete matched wavelet base (redundant dictionary) is a wavelet matrix made up of a series of matched wavelets formed by time shifting, modulating and phase changing the base wavelet ω (t). If omega (t) is a basic wavelet meeting the condition in the matched wavelet base, let
As a control parameter of the matching wavelet, the matching wavelet is expressed as:
where τ represents the time shift, f represents the frequency factor,
expressed as a phase factor. The position of the matched wavelet can be fixed by using time shift, the frequency factor concentrates the energy of the basic wavelet near a few time-frequency points, and the phase factor is used for controlling the shape change of the matched wavelet.
After the overcomplete matching wavelet base phi is built, although the matching wavelets in phi are not all orthogonal, the requirements of I omega are met i || 2 =1. Selecting ω from the matched wavelet base Φ which most closely matches the given signal x (t) 1 The best matching requirement is their inner product<x,ω 1 >To match all atoms in the wavelet base to the largest one of the x (t) inner products:
<x,ω 1 >＞<x,ω i > i≠1 ω i ∈Φ (11)
therefore, after one iteration, the signal x (t) is decomposed into:
x(t)＝<x,ω 1 >ω 1 +R 1 x (12)
R 1 x is the residual signal after the first iteration. R is R 1 x and omega 1 Orthogonalization, therefore:
<x,ω 1 >the larger means ω 1 Is in Hilbert space where the signal is located and
the direction is closest to the unit vector in the matching wavelet base. For residual signal R 1 x repeats the same matching process: screening from phi and R 1 x is the second atom which is most matched, the process is repeated continuously until the kth time (k is more than or equal to 0):
R k x＝<x,ω k-1 >ω k-1 +R k-1 x (14)
after k iterative decomposition calculations, the original signal x (t) can be approximated by a synthesis of k matched wavelets, as residual R 1 And x is small enough, namely, the sparse representation of the signal is good enough to approximate the signal or the residual signal meets a set threshold value, and the iteration is stopped.
According to the high-precision matching wavelets obtained in the steps 103 and 104, the independence of each matching wavelet in the wavelet set is utilized, a needed time-frequency decomposition method is selected, and the time-frequency domain expansion is carried out on each matching wavelet. And obtaining a time-frequency spectrogram of the whole signal through linear superposition. The above loop is repeated for each track in the gather to finally obtain the multi-resolution de-aliased tuning data volume.
In one embodiment, a specific technical solution is described by taking a pre-stack common imaging point gather as an example:
the first step: construction of time-frequency domain wavelet dictionary
In order to enable the sparse representation result to better meet the requirement of analysis signals, the overcomplete atomic dictionary needs to meet certain conditions: the structured wavelet dictionary should match as closely as possible the inherent structure or behavior of the signal; atoms in the dictionary should have a better time-frequency aggregation in order to better accommodate non-stationary signals. The Ricker wavelet is utilized to have the characteristics of simple waveform, faster convergence and short delay time, and has higher time-frequency resolution in thin interbed analysis, so that the Ricker wavelet dictionary is selected to be constructed when the seismic signals are analyzed. As shown in FIG. 1, a portion of the waveforms in the Ricker wavelet dictionary are truncated and displayed.
And a second step of: hilbert transform of a gather
The signal description in the frequency domain is an efficient method, but it belongs to a global transformation and is completely non-localized in the time domain. In the physical and engineering fields, the instantaneous frequency is the most intuitive description of the phenomenon of periodic variation of complex signals. In order to find these transient properties, it is necessary to introduce an analytical signal. The essence of resolving the signal is to convert a narrowband real signal into a complex signal, and the Hilbert transform can generate a unique complex signal for the narrowband real signal. We perform a Hilbert transform on each signal x (t) in the gather to obtain the resolved signal shown in fig. 2.
And a third step of: determination of matching wavelet control parameters
Discretizing the space generated by three control parameters of matched wavelets, any discrete point
A matching wavelet can be generated according to the formula (1) to form a matching wavelet base phi= { omega δ (t)} δ∈Γ Γ is expressed as
Is a set of (3).
And determining the signal sampling rate and the signal length through a time sequence by using the obtained analysis signal. Determining the instantaneous envelope and instantaneous phase of the signal by using the analysis signal, τ k Can be roughly estimated as the corresponding time t at the instantaneous envelope local maximum peak. Scanning frequency parameters of atoms in wavelet base in complex domain, and determining several tau for a certain iteration m Obtaining the frequency a priori value f of the signal at these positions m Assuming that the matched wavelet is zero-phase wavelet, τ is used m and fm Determining a number of a priori information points (tau m ,f m ,0). Selecting a frequency scanning range by taking the prior information points as the center, and determining the peak frequency f of the matched wavelet k 。
FIG. 3 shows that we fix the time t by solving for the instantaneous envelope of the constructed analytic signal, and we obtain the frequency a priori value f of the signal at these locations by solving for the instantaneous frequency at these points m 。
Having determined the center time and dominant frequency of the best matching wavelet ω, the amplitude and phase information of these matching wavelets is next determined. The amplitude and phase information of these matched wavelets can be calculated and obtained using the damped least squares method according to the residual energy minimization principle.
Fourth step: wavelet decomposition based on matching pursuits
The overcomplete matched wavelet base (redundant dictionary) is a wavelet matrix of a series of matched wavelets that are generated substantially by time shifting, modulating and phase changing the basic wavelet ω (t). After the overcomplete matching wavelet base phi is built, although the matching wavelets in phi are not fully orthogonal, the requirements of omega are met i || 2 =1. Selecting ω from the matched wavelet base Φ which most closely matches the given signal x (t) 1 The best matching requirement is their inner product<x,ω 1 >To match all atoms in the wavelet base to the largest one of the x (t) inner products. After k iterative decomposition calculations, the original signal x (t) can be approximated by a synthesis of k matched wavelets, as residual R 1 And x is small enough, namely, the sparse representation of the signal is good enough to approximate the signal or the residual signal meets a set threshold value, and the iteration is stopped.
FIG. 4 is a wavelet set decomposed by a matching pursuit algorithm showing individual matched Ricker wavelets resulting from signal decomposition; fig. 5 is a graph of a match tracking algorithm decomposition reconstruction contrast, and can be intuitively seen that a reconstruction signal obtained by linear superposition of matched wavelets decomposed by the match tracking algorithm is basically identical to an original signal, so that due information can be kept.
Fifth step: wavelet decomposition based on matching pursuits
And (3) according to the high-precision matched wavelets obtained in the steps (3) and (4), utilizing the independence of each matched wavelet in the wavelet set, selecting a needed time-frequency decomposition method, and respectively performing time-frequency domain expansion on the matched wavelets. And obtaining a time-frequency spectrogram of the whole signal through linear superposition. The loop is repeated on each track in the gather to finally obtain the multi-resolution de-aliased tuning data volume.
Fig. 6 is a time-frequency resolution comparison chart of a conventional time-frequency analysis method and a time-frequency analysis method based on a matching pursuit algorithm. The time-frequency diagrams obtained by using short-time Fourier transform and wiener distribution can be seen that STFT has high requirement on window selection, and if the time resolution is to be improved, the window must be as small as possible, which inevitably loses frequency resolution; conversely, if the frequency resolution is to be improved by using a large window, the time resolution must be sacrificed. The WVD analysis has better time-frequency resolution, a window function is not needed to be selected, the reverse sequence of the self signal is used as an analysis time window, the influence of the time window is reduced to the minimum, the product of the time-frequency distribution time length and the bandwidth reaches the minimum, but serious cross-term interference exists. The graphs (c) and (d) are STFT and WVD time-frequency graphs obtained after optimization by using a matching pursuit algorithm, so that the time-frequency resolution is greatly improved, and the frequency is more finely depicted.
FIG. 7 is a time-frequency characterization method based on a matching pursuit algorithm to distinguish thin reservoirs. The traditional time-frequency analysis method can not accurately identify the thin interbed, and the high-resolution time-frequency characterization means based on the matching pursuit algorithm greatly improves the identification capability.
Fig. 8 is a diagram of a multi-resolution de-aliased tuning data volume for a gather. By sequentially expanding each track in the gather, interference effect caused by superposition is avoided, high-precision time-frequency expansion is performed on each track by utilizing a matching pursuit algorithm, time-frequency resolution is improved, influence of frequency aliasing is avoided, and drawing precision of thin reservoir etching is improved.
In the method, the influence of frequency on the seismic attribute of the gather while changing along with the change of offset is considered in the gather processing process before the inversion work is carried out. In addition, superposition of the large-angle gather signals can cause wave interference, influence the quality of seismic data, and the signals can be influenced by frequency aliasing in the acquisition process. Therefore, the method sequentially expands each track in the track set, avoids interference effect caused by superposition, utilizes a matching pursuit algorithm to carry out high-precision time-frequency expansion on each track, improves time-frequency resolution, avoids the influence of frequency aliasing, and improves the precision of the thin reservoir characterization.
Claims (1)
1. The gather tuning optimization dimension increasing method based on sparse representation is characterized by comprising the following steps of:
step 1, constructing a time-frequency domain wavelet dictionary;
step 2, performing Hilbert transform on the gather;
step 3, determining matched wavelet control parameters;
step 4, carrying out wavelet decomposition based on matching tracking;
step 5, performing high-resolution time-frequency decomposition based on the matched wavelets;
in step 1, a Ricker wavelet dictionary is optionally built when analyzing the seismic signals:
wherein ,ωR (t,f i ) Representing a dictionary of Ricker wavelets, t' representing the center time of the Ricker wavelet, f i Representing the dominant frequency of the ith Ricker wavelet in the wavelet dictionary;
the Ricker wavelet has three control parameters, the convergence rate of sidelobe energy is variable, the delay time is controllable, the wavelet base waveform is rich, and the Ricker wavelet is flexibly matched with the seismic wavelet;
in step 2, for a non-stationary signal, its instantaneous frequency is:
wherein ,
is the instantaneous phase, t is the time;
the seismic signal x (t) Hilbert is transformed into:
wherein P is a Cauchy principal value, and τ represents a time shift;
the analytical signal for x (t) is thus obtained as:
in step 3, three control parameters of the matched wavelet are shifted by τ, frequency factor f and phase factor
The generated space discretization, any discretized point +.>
A matching wavelet can be generated according to the formula (1) so as to form a matching wavelet base phi; determining a signal sampling rate and a signal length through a time sequence by using the obtained analytic signal; determining the instantaneous envelope and instantaneous phase of the signal by using the analysis signal, τ k The corresponding time at the peak of the local maximum value of the instantaneous envelope; scanning frequency parameters of atoms in wavelet base in complex domain, determining several tau for a certain iteration m Obtaining the frequency a priori value f of the signal at these positions m Assuming that the matched wavelet is zero-phase wavelet, τ is used m and fm Determining a number of a priori information points (tau m ,f m 0); selecting a frequency scanning range by taking the prior information points as the center, and determining the peak frequency f of the matched wavelet k ；
In the step 3, after the center time and the main frequency of the best matching omega are determined, according to the residual energy minimum principle, the damping minimum level method is utilized to calculate, and the amplitude and phase information of the matching wavelet are removed;
there is a limited bandwidth seismic signal x (t), iterated k-1 times by a matching pursuit algorithm, residual signal R k- 1 x; when iterating to the kth time, M matched wavelets omega are generated in total δm (t), m=1, 2, …, M; for residual signal R k-1 x performs Hilbert transform structure analysis signal:
X＝R k-1 x+iR k-1 x (5)
and constructing complex matched wavelets at the kth iteration in a similar way:
W δm (t)＝ω δm (t)+iω δm (t) (6)
in the formula Wδm (t) is a complex matched wavelet;
residual signal R generated after k iterations k The energy of x is:
wherein X is the N-dimensional column vector of the seismic signal;
in the formula am Matching amplitudes of wavelets for a real number domain;
G δm (t)＝[ω δ (t,f 1 ),ω δ (t,f 2 ),ω δ (t,f 3 ),…,ω δ (t,f m )]for N multiplied by m order matching wavelet matrix, each column corresponds to a matching wavelet;
to minimize the relative error between the seismic signal x (t) and the matched wavelet base Φ, a damped least squares method is used to calculate:
A＝[W T W+εI] -1 W T X (9)
wherein W is an N multiplied by M order complex matching wavelet matrix obtained after k iterations, epsilon is a damping factor, and I is an M multiplied by M order identity matrix; from equations (8) and (9), a control parameter a is determined that matches the wavelet amplitude and phase k And
in step 4, the matched wavelet is formed by time shifting, modulating and phase changing the basic wavelet ω (t); if the basic wavelet omega (t) is one basic wavelet meeting the condition in the matched wavelet base, making
As a control parameter of the matching wavelet, the matching wavelet is expressed as:
where τ represents the time shift, f represents the frequency factor,
expressed as a phase factor; fixing the position of the matched wavelet by using time shift, wherein the frequency factor concentrates the energy of the basic wavelet near a few time-frequency points, and the phase factor is used for controlling the shape change of the matched wavelet;
in step 4, after the overcomplete matching wavelet base Φ is constructed, although the matching wavelets in Φ are not all orthogonal, the requirement is satisfied
Selecting ω from the matched wavelet base Φ which most closely matches a given seismic signal x (t) 1 The best matching requirement is their inner product<x,ω 1 >To match all atoms in the wavelet base to the largest one of the x (t) inner products:
<x,ω 1 >><x,ω i >i≠1 ω i ∈Φ (11)
therefore, after one iteration, the seismic signal x (t) is decomposed into:
x(t)＝<x,ω 1 >ω 1 +R 1 x (12)
R 1 x is the residual signal after the first iteration; r is R 1 x and omega 1 Orthogonalization, therefore:
<x,ω i >the larger means ω 1 Is the unit vector in the Hilbert space where the signal is located and m=1, 2, …, and the M direction is closest to the matched wavelet base; for residual signal R 1 x repeats the same matching process: screening from phi and R 1 x is the second atom which is most matched, the process is repeated continuously until the kth time, and k is more than or equal to 0:
R k x＝<x,ω k-1 >ω k-1 +R k-1 x (14)
after k times of iterative decomposition calculation, the original seismic signal x (t) is approximated by the synthesis of k matched wavelets, and when the residual signal meets a set threshold value, iteration is stopped;
in step 5, according to the matched wavelets obtained in steps 3 and 4, using the independence of each matched wavelet in the wavelet set, selecting a needed time-frequency decomposition method, and respectively performing time-frequency domain expansion on the matched wavelets; obtaining a time-frequency spectrogram of the whole signal through linear superposition; the loop is repeated on each track in the gather to finally obtain the multi-resolution de-aliased tuning data volume.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010013149.9A CN111242198B (en) | 2020-01-06 | 2020-01-06 | Gather tuning optimization dimension increasing method based on sparse representation |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010013149.9A CN111242198B (en) | 2020-01-06 | 2020-01-06 | Gather tuning optimization dimension increasing method based on sparse representation |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111242198A CN111242198A (en) | 2020-06-05 |
| CN111242198B true CN111242198B (en) | 2023-05-09 |
Family
ID=70875976
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010013149.9A Active CN111242198B (en) | 2020-01-06 | 2020-01-06 | Gather tuning optimization dimension increasing method based on sparse representation |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111242198B (en) |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104570107A (en) * | 2013-10-29 | 2015-04-29 | 中国石油化工股份有限公司 | Time-frequency analysis method based on improved matching pursuit algorithm |
| CN109752758A (en) * | 2019-01-16 | 2019-05-14 | 中国石油大学(华东) | A kind of seismic data decomposition method, system and storage medium and terminal |
Family Cites Families (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103293552B (en) * | 2013-05-24 | 2015-10-28 | 中国石油天然气集团公司 | A kind of inversion method of Prestack seismic data and system |
| CN106483561B (en) * | 2015-08-24 | 2019-02-01 | 中国石油化工股份有限公司 | A kind of wavelet decomposition method of frequency domain |
| CN107807390B (en) * | 2016-09-09 | 2019-08-23 | 中国石油化工股份有限公司 | The processing method and system of seismic data |
| CN106597532B (en) * | 2016-11-14 | 2020-06-30 | 中国石油化工股份有限公司 | Pre-stack seismic data frequency band expanding method combining well data and horizon data |
| US10295685B2 (en) * | 2017-04-06 | 2019-05-21 | Saudi Arabian Oil Company | Generating common image gather using wave-field separation |
-
2020
- 2020-01-06 CN CN202010013149.9A patent/CN111242198B/en active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104570107A (en) * | 2013-10-29 | 2015-04-29 | 中国石油化工股份有限公司 | Time-frequency analysis method based on improved matching pursuit algorithm |
| CN109752758A (en) * | 2019-01-16 | 2019-05-14 | 中国石油大学(华东) | A kind of seismic data decomposition method, system and storage medium and terminal |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111242198A (en) | 2020-06-05 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US8538702B2 (en) | Geologic features from curvelet based seismic attributes | |
| Zu et al. | A periodically varying code for improving deblending of simultaneous sources in marine acquisition | |
| US11852767B2 (en) | Robust full waveform inversion of seismic data method and device | |
| Li et al. | Investigation of generalized S-transform analysis windows for time-frequency analysis of seismic reflection data | |
| CN102221708B (en) | Fractional-Fourier-transform-based random noise suppression method | |
| Zhang et al. | Nonstretching NMO correction of prestack time-migrated gathers using a matching-pursuit algorithm | |
| CN102116868A (en) | Seismic wave decomposition method | |
| CN110895348B (en) | Method, system and storage medium for extracting low-frequency information of seismic elastic impedance | |
| CN103728660A (en) | Multi-channel matching tracking method based on seismic data | |
| CN112882099B (en) | Earthquake frequency band widening method and device, medium and electronic equipment | |
| CN105093315B (en) | A method of removal coal seam strong reflectance signal | |
| CN113608259B (en) | Seismic thin layer detection method based on ICEEMDAN constraint generalized S transformation | |
| Aeron et al. | Broadband dispersion extraction using simultaneous sparse penalization | |
| CN111242198B (en) | Gather tuning optimization dimension increasing method based on sparse representation | |
| CN112255690B (en) | Self-adaptive surrounding rock strong reflection separation method based on seismic phase decomposition | |
| Liu | Local time-frequency transform and its application to ground-roll noise attenuation | |
| CN111505709A (en) | Attenuation qualitative analysis method based on sparse spectral decomposition | |
| Xin et al. | Seismic high-resolution processing method based on spectral simulation and total variation regularization constraints | |
| Xu et al. | Enhancing the resolution of time–frequency spectrum using directional multichannel matching pursuit | |
| Zheng et al. | Vibrator data denoising based on fractional wavelet transform | |
| Chen et al. | Nonstationary spectral inversion of seismic data | |
| CN112666603A (en) | Lp norm constraint-based phase identification method and system | |
| Wang et al. | Seismic resolution enhancement with variational modal-based fast-matching pursuit decomposition | |
| Dmitriev et al. | Efficient four-dimensional supergrouping algorithm for enhancement of high-channel count seismic data | |
| CN114428282B (en) | Seismic signal time-frequency conversion method based on downscaling S conversion |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |