CN109375265B - Ideal seismic spectrum decomposing method based on phase-varying Rake wavelet matching pursuit - Google Patents

Ideal seismic spectrum decomposing method based on phase-varying Rake wavelet matching pursuit Download PDF

Info

Publication number
CN109375265B
CN109375265B CN201810962986.9A CN201810962986A CN109375265B CN 109375265 B CN109375265 B CN 109375265B CN 201810962986 A CN201810962986 A CN 201810962986A CN 109375265 B CN109375265 B CN 109375265B
Authority
CN
China
Prior art keywords
wavelet
phase
spectrum
signal
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.)
Expired - Fee Related
Application number
CN201810962986.9A
Other languages
Chinese (zh)
Other versions
CN109375265A (en
Inventor
薛姣
顾汉明
蔡成国
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Geosciences
Original Assignee
China University of Geosciences
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by China University of Geosciences filed Critical China University of Geosciences
Priority to CN201810962986.9A priority Critical patent/CN109375265B/en
Publication of CN109375265A publication Critical patent/CN109375265A/en
Application granted granted Critical
Publication of CN109375265B publication Critical patent/CN109375265B/en
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30—Analysis
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/32—Transforming one recording into another or one representation into another
    • G01V1/325—Transforming one representation into another
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/40—Transforming data representation
    • G01V2210/48—Other transforms

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

The invention discloses an ideal seismic spectrum decomposition method based on phase-varying Rake wavelet matching pursuit, which comprises the steps of firstly constructing a phase-varying Rake wavelet base; then, matching wavelet and corresponding amplitude coefficient are obtained according to a matching tracking method, and the amplitude coefficient is multiplied by wavelet atoms to be used as a sub-signal obtained through decomposition; and then, carrying out Fourier transform on the sub-signals corresponding to the single reflection interface, calculating an amplitude spectrum and a rotation phase spectrum of the interface reflection seismic wave, and placing the Fourier transform spectrum at a corresponding delay time to finally obtain an ideal time frequency spectrum decomposition result. The time-frequency decomposition result comprises an amplitude spectrum and a rotation phase spectrum, and the rotation phase spectrum can provide richer geological and reservoir information. The amplitude spectrum decomposition result accurately reflects the delay time of the seismic wavelet arriving at the interface and the distribution of the energy of the seismic wavelet along with the frequency change, is not limited by the limit of an inaccurate principle, has extremely high time resolution and frequency resolution, and is an ideal time spectrum for seismic signals.

Description

Ideal seismic spectrum decomposing method based on phase-varying Rake wavelet matching pursuit
Technical Field
The invention belongs to the field of oil and gas geophysical exploration, and particularly relates to an ideal earth spectrum decomposition method based on phase-varying Rake wavelet matching pursuit.
Background
The seismic spectrum decomposition technology converts time domain seismic signals into a time-frequency domain for analysis, can acquire richer stratum information from seismic data, and is widely applied to the aspects of determining stratum thickness, sedimentary facies interpretation, oil gas detection and the like.
Common seismic spectrum decomposition methods mainly comprise linear time-frequency analysis, bilinear time-frequency analysis and self-adaptive time-frequency analysis methods. A linear time-frequency analysis method, namely a kernel function decomposition method, realizes the time-frequency locality of signals by using a window function, and mainly comprises short-time Fourier transform (STFT), Continuous Wavelet Transform (CWT), S Transform (ST) and Generalized S Transform (GST). The linear time-frequency analysis method is limited by a Heisenberg uncertainty principle, the time resolution and the frequency resolution of the linear time-frequency analysis method are mutually restricted and contained, and the time resolution and the frequency resolution cannot be simultaneously reduced to any degree. The distribution of signal energy in a time-frequency domain is obtained by the bilinear time-frequency analysis method, and mainly comprises Cohen time-frequency distribution and affine time-frequency distribution. The common Wigner-Ville distribution (WVD) is a Cohen-like time frequency distribution, and compared with a linear time frequency analysis method, the WVD takes a self signal as a window, so that the self-adaptive matching is realized, the limitation of an inaccurate principle is avoided, and the time frequency resolution is higher. Cross terms exist in the WVD decomposition, serious interference is generated on a frequency spectrum decomposition result, and the cross interference terms can be reduced by applying a smooth function to carry out two-dimensional low-pass filtering on the WVD. The higher the smoothing degree of the time-frequency domain, the smaller the influence of the interference term, and the lower the resolution of the time-frequency spectrum, and a trade-off between the joint time-frequency resolution and the influence of the coherent term is needed.
A seismic spectrum decomposition method based on Matching Pursuit (MP) is a self-adaptive time-frequency analysis method, firstly an over-complete wavelet base is constructed, then sparse decomposition of signals is realized based on matching pursuit, and finally time-frequency distribution of each decomposed sub-component signal is obtained by utilizing WVD, and superposition and summation are carried out. The seismic spectrum decomposition result based on matching pursuit has high time-frequency resolution, and avoids the influence of cross terms, thereby being widely applied to the field of seismic data interpretation and inversion. Wavelet atoms commonly used in matching pursuit algorithms include Gabor wavelets, Morlet wavelets, and rake wavelets. The matching tracking algorithm based on the Rake wavelets usually selects the Rake wavelets with zero phase to establish a wavelet base, and calculates the amplitude spectrum and the phase spectrum of the matched wavelets as a whole, wherein the phase spectrum obtained by the method is the linear phase change caused by time delay. However, the zero-phase rake wavelet base has poor redundancy compared to the actual seismic wavelets, and phase-varying rake wavelets need to be constructed. Meanwhile, in the actual seismic recording, not only linear phase change caused by time delay exists, but also phase rotation caused by lamella, velocity dispersion, energy attenuation and the like exists, and conventional matching tracking (no matter zero-phase or non-zero-phase Rake wavelets, real domain or complex domain matching tracking) does not separately consider the rotating phase spectrum.
Disclosure of Invention
In order to overcome the defects of the prior art, the invention provides an ideal seismic spectrum decomposition method based on phase-varying Rake wavelet matching pursuit, which not only has high time-frequency resolution, but also can obtain a rotating phase spectrum, and the rotating phase spectrum contains richer geological information and is expected to play an important role in reservoir prediction and oil gas detection.
The technical scheme adopted by the invention for solving the technical problems is as follows: an ideal seismic spectrum decomposition method based on phase-varying Rake wavelet matching pursuit is constructed, and the method comprises the following steps:
(1) inputting a path of seismic data as an initial residual signal;
(2) constructing a complex analysis signal of a current residual signal, then calculating an instantaneous envelope of the complex analysis signal, and determining the delay time of a matched wavelet atom according to the instantaneous envelope;
(3) constructing a zero-phase Rake wavelet according to the delay time of the matched wavelet atoms, scanning the dominant frequency and the phase of the zero-phase Rake wavelet, and constructing a phase-variable Rake wavelet base;
(4) extracting a wavelet atom which is most matched with a current residual signal from a phase-changing Rake wavelet base, solving an amplitude coefficient corresponding to the most matched wavelet atom, multiplying the amplitude coefficient and the most matched wavelet atom to obtain a sub-signal obtained by decomposition, then subtracting the sub-signal obtained by decomposition from the current residual signal, and taking the obtained residual as a new residual signal;
(5) repeating the steps (2) to (4) until the iteration number reaches a preset value or the residual energy is less than a certain threshold value;
(6) after iteration is finished, obtaining a plurality of phase-variable Rake matching wavelets and corresponding amplitude coefficients, and respectively and correspondingly multiplying the matching wavelets and the amplitude coefficients to obtain single-interface reflection seismic signals serving as a sub-signal set;
(7) and performing Fourier transform on each sub-signal in the sub-signal set, and corresponding the obtained Fourier transform spectrum to the delay time corresponding to the interface to finally obtain an ideal time-frequency spectrum decomposition result.
Further, in the ideal seismic spectrum decomposition method based on phase-varying Rake wavelet matching pursuit of the invention, the step (2) specifically comprises:
constructing a complex analysis signal of the residual signal s (t):
S(t)＝s(t)+isH(t)，
wherein s isH(t) represents the Hilbert transform of the real signal s (t), the instantaneous amplitude envelope of the complex analysis signal being:
Figure BDA0001774267440000031
and taking the time corresponding to the maximum value of the instantaneous amplitude envelope as the delay time of the matched wavelet atom.
Further, in the ideal seismic spectrum decomposition method based on phase-varying Rake wavelet matching pursuit of the invention, the step (3) specifically comprises:
(31) constructing a zero-phase Rake wavelet according to the delay time of the matched wavelet atoms:
Figure BDA0001774267440000032
where τ is the delay time of the matched wavelet atom, fdDominant frequency of the Rake wavelet;
(32) by scanning the time-frequency parameter gamma ═ fdPhi, constructing a phase-variable Rake wavelet atom library, wherein the scanning range of the wavelet main frequency is 0 to Nyquist frequency, the scanning range of the phase is 0-180 degrees, and the constructed phase-variable Rake wavelet library is as follows:
wγ(t)＝w(t-τ,fd)cos(φ)-wH(t-τ,fd)sin(φ)，
where φ is the phase of the Rake wavelet, wH(t) is the Hilbert transform of the zero-phase Rake wavelet w (t).
Further, in the ideal seismic spectrum decomposition method based on phase-varying Rake wavelet matching pursuit of the invention, the step (4) specifically comprises:
for the first random iteration number k, all wavelet atoms in the phase-varying Rake wavelet base and the current residual signal sres k-1(t) performing correlation calculation to obtain a correlation coefficient corr:
Figure BDA0001774267440000041
wherein,
Figure BDA0001774267440000042
representing the current residual signal and the wavelet atom wγ(t) the inner product of the (t),
Figure BDA0001774267440000043
representing wavelet atoms wγ(t) a signal normalization operator;
then, the maximum value of the correlation coefficient and the wavelet atom w corresponding to the maximum value are obtainedγk(t)，
Figure BDA0001774267440000044
The wavelet atom wγk(t) is the wavelet atom which is most matched with the current residual signal, and then the wavelet atom w is obtained through calculationγk(t) amplitude coefficient C corresponding tok：
Figure BDA0001774267440000045
The sub-signal decomposed by the kth iteration is Ckwγk(t) subtracting the k-th decomposed sub-signal from the current residual to obtain a new residual signal sres k(t)：
sres k(t)＝sres k-1(t)-Ckwγk(t)。
Further, in the ideal seismic spectrum decomposition method based on phase-varying Rake wavelet matching pursuit of the invention, the step (6) specifically comprises:
will sk(t)＝Ckwγk(t) as the sub-signals generated in the kth iteration, all sub-signals are formed into a sub-signal set.
Further, in the ideal seismic spectrum decomposition method based on phase-varying Rake wavelet matching pursuit of the invention, the step (7) specifically comprises:
evaluating a subsignal sk(t) Fourier transform spectrum Sk(f)：
Figure BDA0001774267440000051
Wherein f represents the frequency of the fourier transform;
then Fourier transform spectrum Sk(f) Corresponding to the delay time corresponding to the interface to obtain the ideal time spectrum IDFTs(t,f)：
Obtaining the amplitude spectrum and the phase spectrum of the reflected seismic waves of the interface,
further, in the ideal seismic spectrum decomposition method based on phase-varying Rake wavelet matching pursuit, the ideal time frequency spectrum IDFTsThe amplitude spectrum a (t, f) and the phase spectrum Φ (t, f) of (t, f) are:
Figure BDA0001774267440000053
further, in the seismic spectrum decomposition method based on phase-varying Rake wavelet matching pursuit of the invention, the step (7) further comprises obtaining an ideal time frequency spectrum IDFTsRotational phase spectrum Θ (t, f) of (t, f):
Figure BDA0001774267440000055
wherein phik(f) Is a propagation time of τkBest match the phase of the rake wavelet.
Compared with the prior art, the invention has the beneficial effects that:
the amplitude spectrum is formed by Fourier transform of a single-interface fixed propagation time sub-signal, the energy of a reflection wavelet of a single interface is reflected in time, and matching pursuit is a sparse decomposition method, so that the amplitude spectrum has ideal time resolution; in terms of frequency, the Fourier transform is directly carried out on the sub-signals with the same length as the seismic record length, the frequency resolution cannot be reduced due to the influence of the window length, and therefore the ideal frequency resolution is achieved.
The phase spectrum contains two types of phase information, the first is linear phase change due to time delay, and the second is phase rotation due to thin layer, velocity dispersion and energy attenuation. The phase spectrum obtained by the conventional matching tracking based on the zero-phase Rake wavelet only reflects the linear phase change caused by time delay, and the phase spectrum obtained by the matching tracking based on the non-zero-phase Rake wavelet reflects the comprehensive effect of two influences of the linear phase change and the phase rotation. The linear phase change caused by the time delay has no direct relation with the lithology and the reservoir, so for the conventional time frequency analysis, only amplitude spectrum information is generally utilized, and phase spectrum information is less utilized. In fact, in seismic exploration, the rotating phase spectrum can play an important role in formation identification and reservoir prediction because the rotating phase contains information of thin layers, velocity dispersion and energy attenuation.
Drawings
The invention will be further described with reference to the accompanying drawings and examples, in which:
FIG. 1 is a flow chart of an ideal seismic spectrum decomposition method based on phase-varying Rake wavelet matching pursuit according to the present invention;
FIG. 2 is a seismic trace of the model of example 1;
FIG. 3 is a sub-signal set decomposed by the model seismic trace based on the matching pursuit of the phase-varying Rake wavelets of example 1;
FIG. 4(a) is an amplitude spectrum of an ideal seismic spectrum decomposition result obtained by using the ideal seismic spectrum decomposition method based on phase-varying Rake wavelet matching pursuit according to the invention for a model seismic trace of example 1;
FIG. 4(b) is a graphical representation of an ideal amplitude spectrum;
FIG. 4(c) is an ideal rotating phase spectrum;
FIGS. 5(a) and 5(b) are an amplitude spectrum and a phase spectrum obtained by the S-transform time-frequency analysis method for the model seismic traces of example 1;
FIG. 6 shows the result of the matching pursuit-WVD time-frequency decomposition of the model seismic traces of example 1.
Detailed Description
For a more clear understanding of the technical features, objects and effects of the present invention, embodiments of the present invention will now be described in detail with reference to the accompanying drawings.
Referring to fig. 1, fig. 1 is a flowchart of an embodiment of an ideal seismic spectrum decomposition method based on phase-varying raeket wavelet matching pursuit of the present invention, and the method includes the following steps:
and S1, inputting single-channel seismic data as an initial residual signal.
S2, constructing a complex analysis signal of the current residual signal S (t):
S(t)＝s(t)+isH(t)
wherein s isH(t) represents the Hilbert transform of the real signal s (t), then the instantaneous amplitude A (t) of the complex analysis signal is:
Figure BDA0001774267440000071
and taking the time corresponding to the maximum value of the instantaneous amplitude envelope as the delay time tau of the matched wavelet atom.
S3, constructing a zero-phase Rake wavelet according to the delay time tau of the matched wavelet atom calculated in the step S2
w(t)＝{1-2[πfd(t-τ)]2}exp{-[πfd(t-τ)]2}
Wherein, the delay time tau of the matching wavelet atom is the propagation time of the Rake wavelet in the underground medium, fdTo match the dominant frequencies of the rake wavelet. To obtain a Rake wavelet with any constant phase, an analytic signal expression of a Rake wavelet with a zero phase is obtained
W(t)＝w(t)+iwH(t)
Wherein wH(t) is a zero phaseHilbert transform of bit-Rake wavelets. After rotating its phase by phi, the new non-zero phase complex field Rake wavelet is
W'(t)＝W(t)exp(iφ)
Namely, it is
W'(t)＝[w(t)+iwH(t)][cos(φ)+isin(φ)]
Then the time-frequency atoms of the constant-phase Rake wavelet can be represented as
wγ(t)＝w(t-τ,fd)cos(φ)-wH(t-τ,fd)sin(φ)
If the sampling frequency of the signal is fsThe maximum frequency of the atoms of the Rake wavelet is f according to the Nyquist's laws2, in practice, the primary frequency of seismic waves is usually within 100 Hz; the phase interval of the Rake wavelets is 0-180 degrees. Scanning time-frequency parameter gamma ═ tau, fdPhi, then obtaining the Rake wavelet base. The smaller the scanning interval, the more wavelet atoms are generated, the more sparse the decomposition result is, but the calculation amount required by time-frequency decomposition is increased.
In the invention, the phase of the Rake wavelet is assumed to be phi, and the time-frequency parameter gamma is changed to { f through scanningdPhi scans the dominant frequency and phase of the zero-phase Rake wavelet to construct a phase-variable Rake wavelet atom library, wherein the scanning range of the dominant frequency of the wavelet is 0 to Nyquist frequency, the scanning range of the phase is 0 DEG to 180 DEG, and the constructed phase-variable Rake wavelet library wγ(t) is:
wγ(t)＝w(t-τ,fd)cos(φ)-wH(t-τ,fd)sin(φ)
wherein wH(t) is the Hilbert transform of the zero-phase Rake wavelet w (t). Generally, if the sampling frequency of the signal is fsThe frequency sweep interval may be the minimum frequency interval in the fourier transform, i.e., △ f-fsand/N, wherein N is the data length. The scan interval of the phase angle may take pi/36 radians. The smaller the scanning interval, the more wavelet atoms are generated, the more sparse the decomposition result is obtained, but the more time is required for calculation, and a compromise needs to be made between the calculation accuracy and the calculation time.
S4、Combining all wavelets in the Rake wavelet base with the current residual signal sres k-1(t) performing a correlation calculation with a correlation coefficient of
Figure BDA0001774267440000081
Calculating maximum value of correlation coefficient and corresponding wavelet atom
Figure BDA0001774267440000082
Figure BDA0001774267440000083
Where D is a library of matching Rake wavelet atoms, the wavelet atom wγk(t) is the wavelet atom that matches best with the current residual signal, and its corresponding amplitude coefficient is:
Figure BDA0001774267440000091
the signal decomposed by the kth iteration is
Figure BDA0001774267440000092
Subtracting the subsignal decomposed at the k time from the current residual to obtain a new residual signal:
Figure BDA0001774267440000093
and S5, repeating the steps S2 to S4 for the new residual signal until the energy of the residual signal is less than a certain threshold value or the iteration number exceeds the set maximum iteration number, and ending the iteration.
S6, multiplying the decomposed matching Rake wavelets by the corresponding amplitude coefficients to obtain sub signals:
the original seismic signal may be decomposed into sub-signals and decomposed residuals using matching pursuits:
Figure BDA0001774267440000095
where K is the total number of iterations performed, RK(t) is the matching pursuit decomposition residue, which can be ignored, and the in-situ seismic signal is approximately written as:
Figure BDA0001774267440000096
all sub-signals can be regarded as reflection signals of a single interface, and the actual seismic record can be regarded as superposition of a plurality of reflection interface reflection waves.
S7, for all the decomposed sub-signals, the amplitude spectrum and the phase spectrum are obtained by fourier transform. For any sub-signal, its propagation time τkIs calculated according to the extreme value of the instantaneous amplitude of the analytic signal, so that the propagation time of all the sub-signals in the underground medium is known. The conventional time-frequency decomposition method based on matching pursuit usually performs WVD on the decomposed sub-signals to obtain a high-resolution time-frequency spectrum. In fact, according to the stacking principle, the seismic record can be regarded as stacking of a plurality of single interface reflected waves, so that the decomposed wavelets are subjected to fourier transform:
Figure BDA0001774267440000101
and Fourier transform the result Sk(f) The ideal time frequency spectrum IDFT can be obtained by placing the IDFT on the corresponding propagation time of the interfaces(t,f)：
Figure BDA0001774267440000102
Ideally, the spectrum is a complex matrix, whose elements are written in exponential form,
Figure BDA0001774267440000103
the amplitude spectrum a (t, f), the phase spectrum Φ (t, f), and the rotational phase spectrum Θ (t, f) are:
Figure BDA0001774267440000105
Figure BDA0001774267440000106
in one embodiment of the invention, a model seismic trace (FIG. 2) is match tracked (FIG. 3) and its ideal amplitude spectrum and rotational phase spectrum are calculated (see FIGS. 4(a) - (c)). Comparing the spectrum decomposition result (fig. 4) based on the present invention with the spectrum decomposition results of other time-frequency analysis methods (fig. 5(a), 5(b), 6), it can be seen that: the time-frequency resolution of S transformation is relatively poor, and the phase spectrum of the S transformation comprises a linear phase and a rotating phase, so that the practicability is poor; the amplitude spectrum in the WVD decomposition result based on matching pursuit represents the energy of the signal but not the amplitude, and no phase spectrum exists; the spectrum decomposition result of the invention has higher time and spectrum resolution, can obtain the information of the rotation phase and has higher application value.
While the present invention has been described with reference to the embodiments shown in the drawings, the present invention is not limited to the embodiments, which are illustrative and not restrictive, and it will be apparent to those skilled in the art that various changes and modifications can be made therein without departing from the spirit and scope of the invention as defined in the appended claims.

Claims (8)

1. An ideal seismic spectrum decomposition method based on phase-varying Rake wavelet matching pursuit is characterized by comprising the following steps:
(1) inputting a path of seismic data as an initial residual signal;
(2) constructing a complex analysis signal of a current residual signal, then calculating an instantaneous envelope of the complex analysis signal, and determining the delay time of a matched wavelet atom according to the instantaneous envelope;
(3) constructing a zero-phase Rake wavelet according to the delay time of the matched wavelet atoms, scanning the dominant frequency and the phase of the zero-phase Rake wavelet, and constructing a phase-variable Rake wavelet base;
(4) extracting a wavelet atom which is most matched with a current residual signal from a phase-changing Rake wavelet base, solving an amplitude coefficient corresponding to the most matched wavelet atom, multiplying the amplitude coefficient and the most matched wavelet atom to obtain a sub-signal obtained by decomposition, then subtracting the sub-signal obtained by decomposition from the current residual signal, and taking the obtained residual as a new residual signal;
(5) repeating the steps (2) to (4) until the iteration number reaches a preset value or the residual energy is less than a certain threshold value;
(6) after iteration is finished, obtaining a plurality of phase-variable Rake matching wavelets and corresponding amplitude coefficients, and respectively and correspondingly multiplying the matching wavelets and the amplitude coefficients to obtain single-interface reflection seismic signals serving as a sub-signal set;
(7) and performing Fourier transform on each sub-signal in the sub-signal set, and corresponding the obtained Fourier transform spectrum to the delay time corresponding to the interface to finally obtain an ideal time-frequency spectrum decomposition result.
2. The method for decomposing an ideal seismic spectrum based on the phase-varying Ricker wavelet matching pursuit as claimed in claim 1, wherein the step (2) comprises:
constructing a complex analysis signal of the residual signal s (t):
S(t)＝s(t)+isH(t)，
wherein s isH(t) represents the Hilbert transform of the real signal s (t), the instantaneous amplitude envelope of the complex analysis signal being:
Figure FDA0002191051620000011
and taking the time corresponding to the maximum value of the instantaneous amplitude envelope as the delay time of the matched wavelet atom.
3. The method for decomposing an ideal seismic spectrum based on the phase-varying Ricker wavelet matching pursuit as claimed in claim 1, wherein the step (3) comprises:
(31) constructing a zero-phase Rake wavelet according to the delay time of the matched wavelet atoms:
w(t)＝{1-2[πfd(t-τ)]2}exp{-[πfd(t-τ)]2}，
where τ is the delay time of the matched wavelet atom, fd(t- τ) is the dominant frequency of the Rake wavelet at time t- τ;
(32) by scanning the time-frequency parameter gamma ═ fdPhi, constructing a phase-variable Rake wavelet atom library, wherein the scanning range of the wavelet main frequency is 0 to Nyquist frequency, the scanning range of the phase is 0-180 degrees, and the constructed phase-variable Rake wavelet library is as follows:
wγ(t)＝w(t-τ,fd)cos(φ)-wH(t-τ,fd)sin(φ)，
where φ is the phase of the Rake wavelet, wH(t-τ,fd) Is a zero phase Rake wavelet w (t-tau, f)d) Hilbert transform of w (t- τ, f)d) Indicating a dominant frequency of fdA zero-phase rake wavelet with a delay time τ.
4. The method for decomposing an ideal seismic spectrum based on the phase-varying Ricker wavelet matching pursuit as claimed in claim 3, wherein the step (4) comprises:
for the first random iteration number k, all wavelet atoms in the phase-varying Rake wavelet base and the current residual signal sres k-1(t) performing correlation calculation to obtain a correlation coefficient corr:
Figure FDA0002191051620000021
wherein,
Figure FDA0002191051620000022
representing the current residual signal and the wavelet atom wγ(t) the inner product of the (t),
Figure FDA0002191051620000023
representing wavelet atoms wγ(t) a signal normalization operator;
then, the maximum value of the correlation coefficient and the corresponding wavelet atom are obtained
Figure FDA0002191051620000024
Figure FDA0002191051620000025
Wherein D is a matching Rake wavelet atom library;
the wavelet atom
Figure FDA0002191051620000026
I.e. the wavelet atom which is most matched with the current residual signal, and then the wavelet atom is obtained by calculation
Figure FDA0002191051620000027
Corresponding amplitude coefficient Ck：
The sub-signal decomposed by the k iteration is
Figure FDA0002191051620000031
Subtracting the sub-signal decomposed at the k time from the current residual to obtain a new residual signal sres k(t)：
Figure FDA0002191051620000032
5. The method for decomposing an ideal seismic spectrum based on phase-varying Ricker wavelet matching pursuit according to claim 4, wherein the step (6) comprises:
will be provided with
Figure FDA0002191051620000033
As the sub-signals generated in the kth iteration, all the sub-signals constitute a set of sub-signals.
6. The method for decomposing an ideal seismic spectrum based on phase-varying Ricker wavelet matching pursuit according to claim 5, wherein the step (7) comprises:
evaluating a subsignal sk(t) Fourier transform spectrum Sk(f)：
Figure FDA0002191051620000034
Wherein f represents the frequency of the fourier transform;
then Fourier transform spectrum Sk(f) Corresponding to the delay time corresponding to the interface to obtain the ideal time frequency spectrum IDFTs(t,f)：
Figure FDA0002191051620000035
Obtaining the amplitude spectrum and phase spectrum, delta (t-tau), of the reflected seismic wave at the interfacek) Representing a delay time of taukThe unit pulse function of (2).
7. The method of claim 5, wherein the IDFT is an ideal time spectrum IDFTsThe amplitude spectrum a (t, f) and the phase spectrum Φ (t, f) of (t, f) are:
Figure FDA0002191051620000036
Figure FDA0002191051620000037
wherein, δ (t- τ)k) Representing a delay time of taukUnit pulse function of Sk(f) Is a sub-signal sk(t) a Fourier transform spectrum of the (t),
Figure FDA0002191051620000038
is a Fourier transform spectrum Sk(f) The phase of (c).
8. The method of claim 5, wherein the step (7) further comprises obtaining an IDFT (inverse discrete Fourier transform) of the ideal time-frequency spectrumsRotational phase spectrum Θ (t, f) of (t, f):
Figure FDA0002191051620000041
wherein phik(f) Is a propagation time of τkIs best matched to the phase, delta (t-tau), of the Rake waveletk) Representing a delay time of taukThe unit pulse function of (2).
CN201810962986.9A 2018-08-22 2018-08-22 Ideal seismic spectrum decomposing method based on phase-varying Rake wavelet matching pursuit Expired - Fee Related CN109375265B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810962986.9A CN109375265B (en) 2018-08-22 2018-08-22 Ideal seismic spectrum decomposing method based on phase-varying Rake wavelet matching pursuit

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810962986.9A CN109375265B (en) 2018-08-22 2018-08-22 Ideal seismic spectrum decomposing method based on phase-varying Rake wavelet matching pursuit

Publications (2)

Publication Number Publication Date
CN109375265A CN109375265A (en) 2019-02-22
CN109375265B true CN109375265B (en) 2020-01-17

Family

ID=65404522

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810962986.9A Expired - Fee Related CN109375265B (en) 2018-08-22 2018-08-22 Ideal seismic spectrum decomposing method based on phase-varying Rake wavelet matching pursuit

Country Status (1)

Country Link
CN (1) CN109375265B (en)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110133717A (en) * 2019-04-15 2019-08-16 长江大学 Determine the method and apparatus of regional earthquake wave phase
CN111142156B (en) * 2020-01-06 2021-07-20 中国石油化工股份有限公司 Phase control-based seismic strong shielding time-frequency information extraction and stripping method
CN112213782B (en) * 2020-09-29 2022-03-04 中国石油大学(北京) Processing method and device for sub-phase seismic data and server
CN112666609B (en) * 2020-11-30 2022-06-17 中国石油大学(华东) Seismic data low-frequency compensation method and device based on sparse envelope
CN113703049B (en) * 2021-08-27 2024-04-23 中铁二十局集团安哥拉国际有限责任公司 Wavelet estimation method based on n times of Fourier spectrum
CN114384580B (en) * 2021-12-31 2023-05-02 同济大学 Ideal wavelet customizing method based on controllable focus
CN116626760B (en) * 2023-06-02 2024-05-10 四川省自然资源投资集团物探勘查院有限公司 Stratum discontinuity detection method and device based on self-adaptive high-order maximum entropy WVD

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2016063125A1 (en) * 2014-10-23 2016-04-28 Cgg Services Sa Imaging the near subsurface with surface consistent deconvolution operators
CN105022090A (en) * 2015-07-14 2015-11-04 北京博达瑞恒科技有限公司 Wavelet decomposition-based earthquake spectrum decomposition method
CN105353408B (en) * 2015-11-20 2017-10-27 电子科技大学 A kind of Wigner higher-order spectrum seismic signal spectral factorization methods based on match tracing
CN105388527B (en) * 2015-11-30 2018-02-06 中国石油大学(北京) A kind of gas-oil detecting method based on complex field matching pursuit algorithm
CN106855638A (en) * 2016-12-19 2017-06-16 中国石油天然气股份有限公司 Matching pursuit seismic spectrum decomposition method and device
CN107255831A (en) * 2017-06-14 2017-10-17 中国石油化工股份有限公司 A kind of extracting method of prestack frequency dispersion attribute

Also Published As

Publication number Publication date
CN109375265A (en) 2019-02-22

Similar Documents

Publication Publication Date Title
CN109375265B (en) Ideal seismic spectrum decomposing method based on phase-varying Rake wavelet matching pursuit
US11333783B2 (en) Integrated method for estimation of seismic wavelets and synthesis of seismic records in depth domain
Puryear et al. Layer-thickness determination and stratigraphic interpretation using spectral inversion: Theory and application
US8265875B2 (en) Interpolation of periodic data
AU2023202187B2 (en) Use nuos technology to acquire optimized 2d data
Wang Multichannel matching pursuit for seismic trace decomposition
Velis Stochastic sparse-spike deconvolution
Ji et al. Source description of the 1999 Hector Mine, California, earthquake, part I: Wavelet domain inversion theory and resolution analysis
Liu et al. Time-frequency analysis of seismic data using local attributes
US9310502B2 (en) Method of representing seismic signals
US10169652B2 (en) Spatial expansion seismic data processing method and apparatus
EP2689271B1 (en) Simultaneous wavelet extraction and deconvolution in the time domain
Hale A method for estimating apparent displacement vectors from time-lapse seismic images
Li et al. Plane-wave least-squares reverse time migration with a preconditioned stochastic conjugate gradient method
CA2672426A1 (en) Regularisation of irregularly sampled seismic data
van der Baan et al. Nonstationary phase estimation using regularized local kurtosis maximization
EP3273274A1 (en) Device and method for estimating pre-stack wavelet model from seismic gathers
de Matos et al. Detecting stratigraphic discontinuities using time-frequency seismic phase residues
Zhang et al. Depth-domain seismic reflectivity inversion with compressed sensing technique
Cao et al. Joint deblending and data reconstruction with focal transformation
Bonomi et al. Wavefield migration plus Monte Carlo imaging of 3D prestack seismic data
Turco et al. Geostatistical interpolation of non-stationary seismic data
Garabito Prestack seismic data interpolation and enhancement with common‐reflection‐surface–based migration and demigration
CN114779332A (en) Seismic data deposition background removing method and device and electronic equipment
Wu et al. Iterative deblending based on the modified singular spectrum analysis

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
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20200117

CF01 Termination of patent right due to non-payment of annual fee