CN105467442A - A globally optimized time-varying sparse deconvolution method and apparatus - Google Patents

Abstract

The invention discloses a globally optimized time-varying sparse deconvolution method and apparatus. The method comprises the following steps: a frequency domain time-varying part spectrum deconvolution model of unstable state seismic data is constructed; a reflection coefficient pulse number is determined through utilization of the limited update rate principle and an unstable state seismic data effective bandwidth; through a very fast simulated annealing algorithm and the frequency domain time-varying part spectrum deconvolution model, the time position of reflection coefficient pulses of the frequency domain time-varying part spectrum deconvolution model is searched; with regard to the searched time positions, iteration is carried out through utilization of a least square method, and reflection coefficient pulse amplitudes of the time positions are determined; for each iteration, available time position and amplitudes are determined according to system energy; through utilization of a square error of part spectrum observation data and part spectrum synthesized data, whether convergence of the iteration exists is determined; if the convergence exists, iteration stops; and according to the reflection coefficient pulse number and the available time positions and amplitudes, a nonzero reflection coefficient sequence is reconstructured. Through utilization of time varying wavelets, the method is more realistic. The non-linear method and the linearity method are combined, and the calculating precision is raised.

Description

The time thinning thin the Method of Deconvolution and device of global optimization

Technical field

The present invention relates to geophysical survey inversion technique field, particularly relate to a kind of time thinning thin the Method of Deconvolution and device of global optimization.

Background technology

In field of geophysical exploration, the waveform data gathered can regard the addition of waveforms of wavelet when received after underground propagation and boundary reflection that ground excites as, be the result of seismic wavelet and reflection coefficient convolution, namely meet convolution model conventional in subsequently seismic data processing and inverting.Wherein, reflection coefficient is relevant with the elastic characteristic of underground medium, and contain the interfacial information of subsurface lithologic, reflection coefficient has more broadband and more high resolving power compared to waveform data.Therefore, in order to obtain high-resolution reflection coefficient, the conventional method such as deconvolution, inverting carries out asking for of reflection coefficient, to reaching compact wavelet, widening frequency band and putting forward high-resolution effect.

But, conventional the Method of Deconvolution hypothesis seismic data is stable state, underground medium is perfectly elastic, namely wavelet on space-time be do not have vicissitudinous, but (wavelet is in the process of underground propagation for this and actual conditions, the effects such as attenuation by absorption can transform the waveform of wavelet, frequency spectrum and phase place) be not inconsistent, directly affects stability and the accuracy of conventional deconvolution algorithms.Secondly, conventional the Method of Deconvolution processes poststack seismic data in time domain, be equivalent to the seismic data information having used Whole frequency band, but seismic data is lower than being relative inaccurate in about 5Hz and the frequency band higher than about 100Hz, therefore easily introduce high and low frequency noise, the precision of result is restricted.Again, it is a nonlinear process that the position of reflection coefficient is asked for, and use conventional linear to be similar to means and ask for, its result is easily absorbed in local solution, and depends on initial model, and therefore precision is restricted unavoidably.

Summary of the invention

The invention provides a kind of time thinning thin the Method of Deconvolution and device of global optimization, at least to solve in current geologic prospecting process, use conventional the Method of Deconvolution to obtain reflection coefficient, the problem that result precision is not high.

According to an aspect of the present invention, provide a kind of time thinning thin the Method of Deconvolution of global optimization, comprising: the partial spectrum deconvolution model become when building the frequency field of unstable state seismic data; Utilize the effective bandwidth determination reflection coefficient pulse number of limited turnover rate principle and described unstable state seismic data; Utilize very fast simulated annealing algorithm, in conjunction with described partial spectrum deconvolution model, the time location of reflection coefficient pulse is searched for, for each time location searched, utilize least square method to carry out iteration, determine the reflection coefficient pulse amplitude at each time location place; For each iteration, determine available reflection coefficient pulse time positions and amplitude according to system capacity; Utilize the square error of partial spectrum observation data and partial spectrum generated data, judge whether iteration restrains; If convergence, then stop iteration, according to described reflection coefficient pulse number, described available reflection coefficient pulse time positions and amplitude, reconstruct non-zero reflection coefficient sequence.

In one embodiment, the partial spectrum deconvolution model become when building the frequency field of unstable state seismic data, comprising:

When described unstable state seismic data extracts, become seismic wavelet, the expression becoming seismic wavelet time described is as follows: wherein, become the phase spectrum of seismic wavelet during expression, become the spectral amplitude of seismic wavelet when A (ω) represents, HT represents Hilbert transform, and ω represents angular frequency;

Build described partial spectrum deconvolution model in conjunction with becoming seismic wavelet time described, described partial spectrum deconvolution model is as follows:

$S\left(f_m\right)/W\left(f_m\right)=\underset{k=1}{\overset{K}{\Σ}}{a}_k\exp (-i2{\mathrm{\π\τ}}_k{f}_m)+N\left(f_m\right)/W\left(f_m\right),$

Wherein, S (f _m) represent that the Fourier transform of unstable state seismic data is in frequency f _mthe value at place; M=1,2 ..., M, M>=K; M represents at sampling frequency domain (f _low, f _high) frequency points of up-sampling; K represents the reflection coefficient pulse number of non-zero; a _krepresent the amplitude of a kth reflection coefficient pulse; τ _krepresent the time location of a kth reflection coefficient pulse; I represents imaginary unit; N (f _m) represent that the Fourier transform of noise is in frequency f _mthe value at place; W (f _m) become the Fourier transform of seismic wavelet in frequency f when representing _mthe value at place.

In one embodiment, utilize the effective bandwidth determination reflection coefficient pulse number of limited turnover rate principle and described unstable state seismic data, comprising: effective bandwidth is (f _low, f _high) the unstable state seismic data non-zero reflection coefficient pulse number K that can reconstruct be K=[| f _high-f _low| × L/2] _int, wherein, L represents the cycle of non-zero reflection coefficient pulse, [] _intrepresent and round downwards.

In one embodiment, utilize least square method to carry out iteration, determine the reflection coefficient pulse amplitude at each time location place, comprising:

If noiseless and the prior-constrained K of the sparse number of plies is less than the first preset value, adopt the amplitude of following formulae discovery reflection coefficient pulse:

$\widehat{\mathrm{\α}}={\left(G^{T}G\right)}^{-1}{G}^{T}d,$

Wherein, represent the reflection coefficient pulse amplitude vector obtained; $G={\left[\begin{array}{cc}{G}_{mk}^r& G_{mk}^i\end{array}\right]}^T,$ represent complex-exponential function exp (-i2 π τ _kf _m) real part, represent complex-exponential function exp (-i2 π τ _kf _m) imaginary part; τ _krepresent the time location of a kth reflection coefficient pulse; f _mrepresent m sample frequency; D represents known partial spectrum observation data; T representing matrix transposition;

If noise is stronger and the prior-constrained K of the sparse number of plies is greater than the second preset value, adopt the amplitude of following formulae discovery reflection coefficient pulse:

$\widehat{\mathrm{\α}}={(G^{T}G+\mathrm{\λ}I)}^{-1}{G}^{T}d,$

Wherein, λ represents ratio of damping; I representation unit matrix.

In one embodiment, for each iteration, determine available reflection coefficient pulse time positions and amplitude according to system capacity, comprising:

Adopt system capacity J described in following formulae discovery:

J＝||Gα-d|| ²+λ||α|| ²，

Wherein, $G={\left[\begin{array}{cc}{G}_{mk}^r& G_{mk}^i\end{array}\right]}^T,G_{mk}^r=cos\left(2{\mathrm{\π\τ}}_k{f}_m\right)$ Represent complex-exponential function exp (-i2 π τ _kf _m) real part, represent complex-exponential function exp (-i2 π τ _kf _m) imaginary part; τ _krepresent the time location of a kth reflection coefficient pulse; f _mrepresent m sample frequency; α represents the amplitude vector of reflection coefficient pulse, α=[a ₁, a ₂..., a _k] ^t, a _krepresent the amplitude of a kth reflection coefficient pulse; G α expressed portion open score generated data; D represents known partial spectrum observation data; || || ²represent the quadratic sum asking for data, then open radical sign; λ represents ratio of damping;

If described system capacity J reaches preset range, then determine that corresponding reflection coefficient pulse time positions and amplitude are available.

In one embodiment, utilize the square error of partial spectrum observation data and partial spectrum generated data, judge whether iteration restrains, and comprising:

Adopt square error J described in following formulae discovery _a:

J _a＝||Gα-d|| ²，

Wherein, $G={\left[\begin{array}{cc}{G}_{mk}^r& G_{mk}^i\end{array}\right]}^T,G_{mk}^r=cos\left(2{\mathrm{\π\τ}}_k{f}_m\right)$ Represent complex-exponential function exp (-i2 π τ _kf _m) real part, represent complex-exponential function exp (-i2 π τ _kf _m) imaginary part; τ _krepresent the time location of a kth reflection coefficient pulse; f _mrepresent m sample frequency; α represents the amplitude vector of reflection coefficient pulse, α=[a ₁, a ₂..., a _k] ^t, a _krepresent the amplitude of a kth reflection coefficient pulse; G α expressed portion open score generated data; D represents known partial spectrum observation data; || || ²represent the quadratic sum asking for data, then open radical sign;

If described square error J _areach preset threshold range, then determine iteration convergence.

In one embodiment, according to described reflection coefficient pulse number, described available reflection coefficient pulse time positions and amplitude, reconstruct non-zero reflection coefficient sequence, comprising:

Described reflection coefficient pulse number, described available reflection coefficient pulse time positions and amplitude are brought in non-zero reflection coefficient sequence expression formula, the expression formula of described non-zero reflection coefficient sequence r (t) is as follows:

$r\left(t\right)=\underset{k=1}{\overset{K}{\Σ}}{a}_k\mathrm{\δ}(t-{\mathrm{\τ}}_k),$

Wherein, K represents the reflection coefficient pulse number of non-zero; τ _krepresent the time location of a kth reflection coefficient pulse; a _krepresent the amplitude of a kth reflection coefficient pulse; δ (t) indicating impulse function.

In one embodiment, described method also comprises: in an iterative process, if the maximum iteration time that iterations reaches default detected, then stops iteration.

According to another aspect of the present invention, provide a kind of global optimization time become sparse deconvolution device, comprising: construction unit, the partial spectrum deconvolution model become during for building the frequency field of unstable state seismic data; First determining unit, for utilizing the effective bandwidth determination reflection coefficient pulse number of limited turnover rate principle and described unstable state seismic data; Iteration unit, for utilizing very fast simulated annealing algorithm, in conjunction with described partial spectrum deconvolution model, the time location of reflection coefficient pulse is searched for, for each time location searched, utilize least square method to carry out iteration, determine the reflection coefficient pulse amplitude at each time location place; Second determining unit, for for each iteration, determines available reflection coefficient pulse time positions and amplitude according to system capacity; Judging unit, for utilizing the square error of partial spectrum observation data and partial spectrum generated data, judges whether iteration restrains; Reconfiguration unit, for when iteration convergence, stops iteration, according to described reflection coefficient pulse number, described available reflection coefficient pulse time positions and amplitude, and reconstruct non-zero reflection coefficient sequence.

In one embodiment, described construction unit comprises:

Extraction module, for becoming seismic wavelet when described unstable state seismic data extracts, the expression becoming seismic wavelet time described is as follows: wherein, become the phase spectrum of seismic wavelet during expression, become the spectral amplitude of seismic wavelet when A (ω) represents, HT represents Hilbert transform, and ω represents angular frequency;

Build module, for building described partial spectrum deconvolution model in conjunction with becoming seismic wavelet time described, described partial spectrum deconvolution model is as follows:

$S\left(f_m\right)/W\left(f_m\right)=\underset{k=1}{\overset{K}{\Σ}}{a}_k\exp (-i2{\mathrm{\π\τ}}_k{f}_m)+N\left(f_m\right)/W\left(f_m\right),$

Wherein, S (f _m) represent that the Fourier transform of unstable state seismic data is in frequency f _mthe value at place; M=1,2 ..., M, M>=K; M represents at sampling frequency domain (f _low, f _high) frequency points of up-sampling; K represents the reflection coefficient pulse number of non-zero; a _krepresent the amplitude of a kth reflection coefficient pulse; τ _krepresent the time location of a kth reflection coefficient pulse; I represents imaginary unit; N (f _m) represent that the Fourier transform of noise is in frequency f _mthe value at place; W (f _m) become the Fourier transform of seismic wavelet in frequency f when representing _mthe value at place.

By the time thinning thin the Method of Deconvolution and device of global optimization of the present invention, consider the deformation problems of wavelet in underground propagation process, include the change of amplitude and phase place, reflect actual conditions more accurately, improve the precision asking for reflection coefficient.Further, reflection coefficient asked for the optimization problem being divided into position and size, tackle nonlinear position optimization problem with non-linear method, tackle linear amplitude size with linear method and ask for problem, suit measures to local conditions, shoot the arrow at the target, not only improve precision, have also contemplated that the problem of calculated amount.In addition, by the process of simulating nature circle metal annealing crystallization, the searching method of acquisition more meets spontaneous phenomenon, can obtain globally optimal solution, improves the precision of reflection coefficient location estimation.Owing to make use of partial spectrum information, can effectively avoid out-of-band noise, particularly to the high frequency noise that deconvolution effect has the greatest impact, improve the degree of stability of deconvolution process.

Accompanying drawing explanation

Accompanying drawing described herein is used to provide a further understanding of the present invention, and form a application's part, schematic description and description of the present invention, for explaining the present invention, does not form limitation of the invention.In the accompanying drawings:

Fig. 1 be the global optimization of the embodiment of the present invention time thinning thin the Method of Deconvolution process flow diagram;

Fig. 2 be the global optimization of the embodiment of the present invention time become the structured flowchart of sparse deconvolution device;

Fig. 3 is the schematic diagram of the true reflection coefficient sequence model of the embodiment of the present invention;

Fig. 4 A is the waveform schematic diagram of the time-varying wavelet that the embodiment of the present invention adopts;

Fig. 4 B be the embodiment of the present invention Fig. 4 A shown in spectral amplitude corresponding to time-varying wavelet;

Fig. 5 is the seismologic record data schematic diagram with decay characteristics of embodiment of the present invention synthesis;

Fig. 6 is the comparison diagram of the reflection coefficient that the reflection coefficient that obtains of the present invention and conventional method obtain;

Fig. 7 A is the waveform schematic diagram of the time-varying wavelet that the actual seismic data of the embodiment of the present invention obtains;

Fig. 7 B be the embodiment of the present invention Fig. 7 A shown in spectral amplitude corresponding to time-varying wavelet;

Fig. 8 is the design sketch that actual seismic data is obtained by conventional method;

Fig. 9 is the design sketch that actual seismic data is obtained by the inventive method.

Embodiment

Below in conjunction with the accompanying drawing in the embodiment of the present invention, be clearly and completely described the technical scheme in the embodiment of the present invention, obviously, described embodiment is only the present invention's part embodiment, instead of whole embodiments.Based on embodiments of the invention, those of ordinary skill in the art, not making the every other embodiment obtained under creative work prerequisite, belong to protection scope of the present invention.

Embodiments provide a kind of time thinning thin the Method of Deconvolution of global optimization, Fig. 1 be the global optimization of the embodiment of the present invention time thinning thin the Method of Deconvolution process flow diagram, as shown in Figure 1, the method comprises the following steps:

Step S101, the partial spectrum deconvolution model become when building the frequency field of unstable state seismic data.

Step S102, utilizes the effective bandwidth determination reflection coefficient pulse number of limited turnover rate principle and unstable state seismic data.

Step S103, utilizes very fast simulated annealing algorithm, joint portion open score deconvolution model, the time location of reflection coefficient pulse is searched for, for each time location searched, utilize least square method to carry out iteration, determine the reflection coefficient pulse amplitude at each time location place.

Step S104, for each iteration, determines available reflection coefficient pulse time positions and amplitude according to system capacity, such as, can judge whether system capacity meets pre-conditioned, if met, then represents that corresponding solution is available.

Step S105, utilizes the square error of partial spectrum observation data and partial spectrum generated data, judges whether iteration restrains.

Step S106, if convergence, then stops iteration, according to reflection coefficient pulse number, available reflection coefficient pulse time positions and amplitude, and reconstruct non-zero reflection coefficient sequence.

By said method, consider the deformation problems of wavelet in underground propagation process, include the change of amplitude and phase place, reflect actual conditions more accurately, improve the precision asking for reflection coefficient.Further, reflection coefficient asked for the optimization problem being divided into position and size, tackle nonlinear position optimization problem with non-linear method, tackle linear amplitude size with linear method and ask for problem, suit measures to local conditions, shoot the arrow at the target, not only improve precision, have also contemplated that the problem of calculated amount.In addition, by the process of simulating nature circle metal annealing crystallization, the searching method of acquisition more meets spontaneous phenomenon, can obtain globally optimal solution, improves the precision of reflection coefficient location estimation.Owing to make use of partial spectrum information, can effectively avoid out-of-band noise, particularly to the high frequency noise that deconvolution effect has the greatest impact, improve the degree of stability of deconvolution process.

The present invention can also obtain multiple inversion result, to carry out inversion result uncertainty analysis by the different initial model of input.

In one embodiment, step S101 can comprise: become seismic wavelet (can referred to as time-varying wavelet) when extracting from unstable state seismic data, in conjunction with time become seismic wavelet and build partial spectrum deconvolution model.

First the extraction becoming seismic wavelet below time pair is described.Seismic wavelet extraction mode is become theoretical based on time frequency analysis time conventional, composed by the means such as spectrum analog and horizontal weighting estimator wave amplitude, become seismic wavelet when being obtained by hypothesis sub-wave phase again, unlike this, when the following two kinds method can be adopted in the present embodiment to extract, become seismic wavelet.

(1) adopt segmentation to estimate the mode of wavelet, time dependent unstable state seismic data is regarded in the less time period as stable state, like this, in often a bit of, carry out wavelet extraction, time dependent wavelet storehouse can be obtained.Concrete, seismic data is windowed process, is divided into the time window that multiple size is identical, when selected, adopt statistically method to estimate wavelet in window, change seismic wavelet when finally obtaining.

When supposing each, in window, wavelet is minimum phase, and so the spectral amplitude of wavelet has relation one to one, the Hilbert transform each other of the intermediary heat spectrum of spectral amplitude and phase spectrum with its phase spectrum.Therefore, wavelet when can synthesize corresponding by Hilbert transform in window, that is:

Wherein, become the phase spectrum of seismic wavelet during expression, become the spectral amplitude of seismic wavelet when A (ω) represents, HT represents Hilbert transform, and ω represents angular frequency.

This method practicality is stronger, and more stable, the wavelet of acquisition can meet the deconvolution needs of this interval better.

(2) according to the phase spectrum of normal phasescan method determination time-varying wavelet.Deconvolution is generally assumed to be basis with Chang Xiangwei, so just can go the accuracy of measuring phase place by a given objective function with phase place change.By constantly revising phase place, obtain the phase place making the minimization of object function, namely best phase value.Finally, adopt when this phase place and the synthesis of seismic data spectral amplitude and become seismic wavelet.Wherein, this seismic data spectral amplitude can obtain by carrying out Fourier transform to seismic data.

Secondly in conjunction with time become the partial spectrum deconvolution model become when seismic wavelet builds frequency field and be described.

Suppose that the subsurface reflective boundary in current selected window has K, namely the reflection coefficient pulse of non-zero is K, and so reflection coefficient sequence just can be expressed as:

$r\left(t\right)=\underset{k=1}{\overset{K}{\Σ}}{a}_k\mathrm{\δ}(t-{\mathrm{\τ}}_k)---\left(2\right)$

Wherein, r (t) represents reflection coefficient sequence; K represents the reflection coefficient pulse number of non-zero; τ _krepresent the time location of a kth reflection coefficient pulse; a _krepresent the amplitude of a kth reflection coefficient pulse; δ (t) indicating impulse function, t represents the time.When K is much smaller than sampled point number, formula (2) can regard the sparse constraint to reflected impulse as.

Suppose that seismic data can be described by Robinson convolution model, so seismic data can be represented by following formula:

$s\left(t\right)=w\left(t\right)\×r\left(t\right)+n\left(t\right)=\underset{k=1}{\overset{K}{\Σ}}{a}_{k}w(t-{\mathrm{\τ}}_k)+n\left(t\right)---\left(3\right)$

Wherein, s (t) represents seismologic record (or seismic data); Become the convolution matrix of seismic wavelet storehouse composition when w (t) represents, need to use formula (1) herein; N (t) represents noise.

Formula (3) the right and left is carried out Fourier transform simultaneously, and the convolution model obtaining frequency field is shown below:

$S\left(f\right)=W\left(f\right)R\left(f\right)+N\left(f\right)=W\left(f\right)\underset{k=1}{\overset{K}{\Σ}}{a}_k\exp (-i2{\mathrm{\π\τ}}_{k}f)+N\left(f\right)---\left(4\right)$

Wherein, f represents frequency, and S (f) represents the Fourier spectrum of seismologic record; W (f) represents the Fourier spectrum in wavelet storehouse; R (f) represents the Fourier spectrum of reflection coefficient; N (f) represents the Fourier spectrum of noise; I represents the imaginary unit of plural number; The exponential function that it is the end that exp represents with natural constant e.

If by frequency with sampling interval Δ f at frequency domain (f _low, f _high) on get the different frequency of M, then f _m=(m-1) × Δ f+f ₁, m=1,2 ..., M, M>=K, M represents at sampling frequency domain (f _low, f _high) frequency points of up-sampling.In order to make full use of the information of partial spectrum, usual frequency chooses the dominant frequency that need meet seismic data, and frequency sampling interval delta f need meet Δ f≤1/L, and L represents the cycle of non-zero reflection coefficient pulse.Obviously, [f ₁, f _m] be (f _low, f _high) subset, even if work as | f _m-f ₁| when × L/2 is greater than the real reflection coefficient pulse number in underground, also can reconstruct all stratum reflection coefficient pulses exactly.

Thus, eliminate the impact of wavelet, the partial spectrum deconvolution equation that can construct unstable state seismic data is as follows:

$S\left(f_m\right)/W\left(f_m\right)=\underset{k=1}{\overset{K}{\Σ}}{a}_k\exp (-i2{\mathrm{\π\τ}}_k{f}_m)+N\left(f_m\right)/W\left(f_m\right)---\left(5\right)$

Wherein, N (f _m) represent that the Fourier transform of noise is in frequency f _mthe value at place; W (f _m) become the Fourier transform of seismic wavelet in frequency f when representing _mthe value at place; S (f _m) represent that the Fourier transform of unstable state seismic data is in frequency f _mthe value at place; $\tilde{N}\left(f_m\right)=N\left(f_m\right)/W\left(f_m\right).$

Write formula (5) as matrix form as follows:

$\tilde{R}\left(f_m\right)=G\mathrm{\α}+\tilde{N}={\left[\begin{array}{cc}{R}^r& R^i\end{array}\right]}^T={\left[\begin{array}{cc}{G}_{mk}^r& G_{mk}^i\end{array}\right]}^T\mathrm{\α}+{\left[\begin{array}{cc}{\tilde{N}}^r& {\tilde{N}}^i\end{array}\right]}^T---\left(6\right)$

Wherein, $\tilde{R}\left(f_m\right)={\left[\begin{array}{cc}{R}^r& R^i\end{array}\right]}^T$ Represent the partial spectrum observation data of reflection coefficient; R ^r=real (S (f _m)/W (f _m)) real part of expressed portion open score; R ⁱ=imag (S (f _m)/W (f _m)) imaginary part of expressed portion open score; $G={\left[\begin{array}{cc}{G}_{mk}^r& G_{mk}^i\end{array}\right]}^T;G_{mk}^r=cos\left(2{\mathrm{\π\τ}}_k{f}_m\right)$ Represent complex-exponential function exp (-i2 π τ _kf _m) real part; represent complex-exponential function exp (-i2 π τ _kf _m) imaginary part, the value of G is subject to the time location τ of reflection coefficient pulse _kcontrol; α represents the amplitude vector of reflection coefficient pulse, α=[a ₁, a ₂..., a _k] ^t; $\tilde{N}={\left[\begin{array}{cc}{\tilde{N}}^r& {\tilde{N}}^i\end{array}\right]}^T$ Represent the partial spectrum observation data of noise; ${\tilde{N}}^r=real(N\left(f_m\right)/W\left(f_m\right))$ Represent the real part of noise section spectrum vector; represent the imaginary part of noise section spectrum vector; Real part is got in real () expression; Imaginary part is got in imag () expression; T represents matrix transpose.

In one embodiment, step S102 can comprise: according to limited turnover rate sampling and the reconstruction theory of non-band-limit signal, effective bandwidth is (f _low, f _high) the unstable state seismic data non-zero reflection coefficient pulse number K that can reconstruct at most be K=[| f _high-f _low| × L/2] _int, wherein, L represents the cycle of non-zero reflection coefficient pulse, [] _intrepresent and round downwards.

In the present embodiment, based on the limited turnover rate principle (Vetterli in Bandlimited Signals sampling and reconstruction theory, 2002), the prior estimate to reflection coefficient pulse number can be obtained, the situation selecting the time variant deconvolution result that accidentally causes inaccurate and unstable due to pulse number can be effectively reduced like this.If the effective bandwidth (f of seismic signal _low, f _high) as the sampling bandwidth of sparse stratum pulse-echo coefficient.When processing actual seismic data, seismic signal in assumption period L is needed to be the result that pulse-echo coefficient is steadily sampled.Usually the approximate value of the seismic data length effective writing time in this cycle as period L can be got.

The partial spectrum deconvolution equation of the matrix form according to formula (6), G is the matrix of coefficients of 2M × K; D is the observation data vector of 2M × 1; α is the reflection coefficient amplitude vector of K × 1.If known observation data d, then formula (6) is containing the time location { τ by reflection coefficient pulse _kand amplitude { a _k2K unknown number forming, and unknown number { τ _kin be contained in matrix G, this illustrate Matrix division (6) be a Nonlinear System of Equations.

By the square matches error J of known partial spectrum observation data d and partial spectrum generated data G α _asetting up objective function is:

J _a＝||Gα-d|| ²(7)

Usually, the time location of reflection coefficient pulse and amplitude can by making J _aobtain minimum way to solve.If the partial spectrum observation data of noise meet zero-mean gaussian distribution and variances sigma ²known, then make the square error of Data Matching expect to meet E (J _a)≤2M σ ²time, the time location { τ of K the reflection coefficient pulse asked for _kand amplitude { a _kcan as of an equation feasible solution.

In fact, the time location { τ of reflection coefficient pulse is asked for _kthe process of a non-linear optimizing, and ask the process of amplitude to be linear at the time location place of given reflection coefficient pulse.Therefore, non-linear optimization method and linear least square method combine by the present invention, estimate best reflection coefficient amplitude while searching out reflection coefficient optimum position.

The present invention adopts the very fast simulated annealing algorithm (VFSA) of global optimization method to carry out the time location { τ of reflection coefficient pulse _ksearch, at the time location that each simulated annealing searches, estimate the amplitude { a of corresponding position with linear least square method _k.Loop iteration, until search best reflection coefficient pulse time positions and amplitude.

In one embodiment, the reflection coefficient pulse amplitude at each time location place can be determined in step S103 as follows, comprising:

If noiseless and the prior-constrained K of the sparse number of plies is less than the first preset value, adopt with the amplitude of following formula (8) computational reflect coefficient pulse.

By $\frac{\∂J_a}{\∂\widehat{\mathrm{\α}}}=2G^T(G\widehat{\mathrm{\α}}-d)=0$ Derive following formula:

$\widehat{\mathrm{\α}}={\left(G^{T}G\right)}^{-1}{G}^{T}d---\left(8\right)$

Wherein, represent the reflection coefficient pulse amplitude vector obtained; $G={\left[\begin{array}{cc}{G}_{mk}^r& G_{mk}^i\end{array}\right]}^T,$ $G_{mk}^r=cos\left(2{\mathrm{\π\τ}}_k{f}_m\right)$ Represent complex-exponential function exp (-i2 π τ _kf _m) real part, $G_{mk}^i=-sin\left(2{\mathrm{\π\τ}}_k{f}_m\right)$ Represent complex-exponential function exp (-i2 π τ _kf _m) imaginary part; τ _krepresent the time location of a kth reflection coefficient pulse; f _mrepresent m sample frequency; D represents known partial spectrum observation data; T representing matrix transposition;

If noise is comparatively strong (such as, noise is greater than a certain preset value) and the prior-constrained K of the sparse number of plies is greater than the second preset value (when namely K is larger), the further stabilization result of damped least squares method can be adopted, namely according to the amplitude of following formulae discovery reflection coefficient pulse:

$\widehat{\mathrm{\α}}={(G^{T}G+\mathrm{\λ}I)}^{-1}{G}^{T}d---\left(9\right)$

Wherein, λ represents ratio of damping (quadratic form regularization factors); I representation unit matrix.

Obviously, due to the finiteness of seismic data bandwidth, the frequency points chosen and ask for sparsely layer by layer count prior-constrained K be little relative to seismologic record sampling number.Therefore, the calculated amount of formula (9) is also little.

In one embodiment, step S104 can comprise: according to the system capacity J of formula (7) and (9) calculating simulation annealing algorithm, if system capacity J reaches preset range, then determine that corresponding reflection coefficient pulse time positions and amplitude are available.

Concrete, following formulae discovery system capacity J can be adopted:

J＝||Gα-d|| ²+λ||α|| ²(10)

Wherein, $G={\left[\begin{array}{cc}{G}_{mk}^r& G_{mk}^i\end{array}\right]}^T,G_{mk}^r=cos\left(2{\mathrm{\π\τ}}_k{f}_m\right)$ Represent complex-exponential function exp (-i2 π τ _kf _m) real part, represent complex-exponential function exp (-i2 π τ _kf _m) imaginary part; τ _krepresent the time location of a kth reflection coefficient pulse; f _mrepresent m sample frequency; α represents the amplitude vector of reflection coefficient pulse, α=[a ₁, a ₂..., a _k] ^t, a _krepresent the amplitude of a kth reflection coefficient pulse; G α expressed portion open score generated data; D represents known partial spectrum observation data; || || ²represent the quadratic sum asking for data, then open radical sign; λ represents ratio of damping.

In one embodiment, step S105 can comprise: calculate square error J _aif, square error J _areach preset threshold range, then determine iteration convergence.

Concrete, following formulae discovery square error J can be adopted _a:

J _a＝||Gα-d|| ²(11)

Wherein, $G={\left[\begin{array}{cc}{G}_{mk}^r& G_{mk}^i\end{array}\right]}^T,G_{mk}^r=cos\left(2{\mathrm{\π\τ}}_k{f}_m\right)$ Represent complex-exponential function exp (-i2 π τ _kf _m) real part, represent complex-exponential function exp (-i2 π τ _kf _m) imaginary part; τ _krepresent the time location of a kth reflection coefficient pulse; f _mrepresent m sample frequency; α represents the amplitude vector of reflection coefficient pulse, α=[a ₁, a ₂..., a _k] ^t, a _krepresent the amplitude of a kth reflection coefficient pulse; G α expressed portion open score generated data; D represents known partial spectrum observation data; || || ²represent the quadratic sum asking for data, then open radical sign.

In one embodiment, step S106 can comprise: be brought in non-zero reflection coefficient sequence expression formula by reflection coefficient pulse number, available reflection coefficient pulse time positions and amplitude, i.e. formula shown in formula (2).

In one embodiment, said method can also comprise: in an iterative process, if the maximum iteration time that iterations reaches default detected, then stops iteration.

Obviously, based on the stochastic inverse of the Sparse Pulse reflection coefficient of partial spectrum be the process of an iterative.In every single-step iteration, need the time location searching for K reflection coefficient upgrade matrix of coefficients G, then solving equation (8) or (9) obtain the amplitude of reflection coefficient calculate current system energy J, then differentiate whether receive current solution by simulated annealing criterion.Square error can be calculated, to judge whether current iteration restrains according to formula (11).If square error J _areach given limits of error ε or iterations reaches given maximum iteration time, iteration stopping, export the non-zero pulses reflection coefficient of reconstruct.

Method flow of the present invention is as follows: choose time variant deconvolution frequency band, extract time-varying wavelet, the partial spectrum deconvolution model become when obtaining frequency field according to time-varying wavelet, the pulse number become when utilizing limited turnover rate to determine, the reflecting interface position become when providing, the amplitude size at initialization TVR time variable reflection interface, interactive refreshing reflecting interface position and optimization amplitude size, the equal computing system energy of each iteration, finally carries out convergence judgement, if convergence, export optimum time variant deconvolution result; If do not restrained, return and continue interactive refreshing reflecting interface position and optimize amplitude size.

Based on same inventive concept, the embodiment of the present invention additionally provide a kind of global optimization time become sparse deconvolution device, may be used for the method realized described by above-described embodiment.The principle of dealing with problems due to this device is similar to said method, and therefore the enforcement of this device see the enforcement of said method, can repeat part and repeat no more.Following used, term " unit " or " module " can realize the software of predetermined function and/or the combination of hardware.Although the system described by following examples preferably realizes with software, hardware, or the realization of the combination of software and hardware also may and conceived.

Fig. 2 be the global optimization of the embodiment of the present invention time become the structured flowchart of sparse deconvolution device, as shown in Figure 2, this device comprises: construction unit 21, first determining unit 22, iteration unit 23, second determining unit 24, judging unit 25 and reconfiguration unit 26, be specifically described this structure below.

Construction unit 21, the partial spectrum deconvolution model become during for building the frequency field of unstable state seismic data;

First determining unit 22, for utilizing the effective bandwidth determination reflection coefficient pulse number of limited turnover rate principle and unstable state seismic data;

Iteration unit 23, for utilizing very fast simulated annealing algorithm, joint portion open score deconvolution model, the time location of reflection coefficient pulse is searched for, for each time location searched, utilize least square method to carry out iteration, determine the reflection coefficient pulse amplitude at each time location place;

Second determining unit 24, for for each iteration, determines available reflection coefficient pulse time positions and amplitude according to system capacity;

Judging unit 25, for utilizing the square error of partial spectrum observation data and partial spectrum generated data, judges whether iteration restrains;

Reconfiguration unit 26, for when iteration convergence, stops iteration, according to reflection coefficient pulse number, available reflection coefficient pulse time positions and amplitude, and reconstruct non-zero reflection coefficient sequence.

By said apparatus, consider the deformation problems of wavelet in underground propagation process, include the change of amplitude and phase place, reflect actual conditions more accurately, improve the precision asking for reflection coefficient.Further, reflection coefficient asked for the optimization problem being divided into position and size, tackle nonlinear position optimization problem with non-linear method, tackle linear amplitude size with linear method and ask for problem, suit measures to local conditions, shoot the arrow at the target, not only improve precision, have also contemplated that the problem of calculated amount.In addition, by the process of simulating nature circle metal annealing crystallization, the searching method of acquisition more meets spontaneous phenomenon, can obtain globally optimal solution, improves the precision of reflection coefficient location estimation.Owing to make use of partial spectrum information, can effectively avoid out-of-band noise, particularly to the high frequency noise that deconvolution effect has the greatest impact, improve the degree of stability of deconvolution process.

In one embodiment, construction unit 21 can comprise: extraction module and structure module.

Extraction module, becomes seismic wavelet during for extracting from unstable state seismic data, time to become the expression of seismic wavelet as follows: wherein, become the phase spectrum of seismic wavelet during expression, become the spectral amplitude of seismic wavelet when A (ω) represents, HT represents Hilbert transform, and ω represents angular frequency;

Build module, in conjunction with time become seismic wavelet and build partial spectrum deconvolution model, partial spectrum deconvolution model is as follows:

$S\left(f_m\right)/W\left(f_m\right)=\underset{k=1}{\overset{K}{\Σ}}{a}_k\exp (-i2{\mathrm{\π\τ}}_k{f}_m)+N\left(f_m\right)/W\left(f_m\right),$

Wherein, S (f _m) represent that the Fourier transform of unstable state seismic data is in frequency f _mthe value at place; M=1,2 ..., M, M>=K; M represents at sampling frequency domain (f _low, f _high) frequency points of up-sampling; K represents the reflection coefficient pulse number of non-zero; a _krepresent the amplitude of a kth reflection coefficient pulse; τ _krepresent the time location of a kth reflection coefficient pulse; I represents imaginary unit; N (f _m) represent that the Fourier transform of noise is in frequency f _mthe value at place; W (f _m) become the Fourier transform of seismic wavelet in frequency f when representing _mthe value at place.

In one embodiment, the first determining unit 22 specifically for: effective bandwidth is (f _low, f _high) the unstable state seismic data non-zero reflection coefficient pulse number K that can reconstruct be K=[| f _high-f _low| × L/2] _int, wherein, L represents the cycle of non-zero reflection coefficient pulse, [] _intrepresent and round downwards.

In one embodiment, iteration unit 23 can comprise: the first computing module and the second computing module.

First computing module, in noiseless and the prior-constrained K of the sparse number of plies is less than the first preset value, adopt the amplitude of following formulae discovery reflection coefficient pulse:

$\widehat{\mathrm{\α}}={\left(G^{T}G\right)}^{-1}{G}^{T}d,$

Wherein, represent the reflection coefficient pulse amplitude vector obtained; $G={\left[\begin{array}{cc}{G}_{mk}^r& G_{mk}^i\end{array}\right]}^T,$ $G_{mk}^r=cos\left(2{\mathrm{\π\τ}}_k{f}_m\right)$ Represent complex-exponential function exp (-i2 π τ _kf _m) real part, $G_{mk}^i=-sin\left(2{\mathrm{\π\τ}}_k{f}_m\right)$ Represent complex-exponential function exp (-i2 π τ _kf _m) imaginary part; τ _krepresent the time location of a kth reflection coefficient pulse; f _mrepresent m sample frequency; D represents known partial spectrum observation data; T representing matrix transposition;

Second computing module, for when noise is stronger and the prior-constrained K of the sparse number of plies is greater than the second preset value, adopts the amplitude of following formulae discovery reflection coefficient pulse:

$\widehat{\mathrm{\α}}={(G^{T}G+\mathrm{\λ}I)}^{-1}{G}^{T}d,$

Wherein, λ represents ratio of damping; I representation unit matrix.

In one embodiment, the second determining unit 24 can comprise: the 3rd computing module and the first judge module.

3rd computing module, specifically for adopting following formulae discovery system capacity J:

J＝||Gα-d|| ²+λ||α|| ²，

Wherein, $G={\left[\begin{array}{cc}{G}_{mk}^r& G_{mk}^i\end{array}\right]}^T,G_{mk}^r=cos\left(2{\mathrm{\π\τ}}_k{f}_m\right)$ Represent complex-exponential function exp (-i2 π τ _kf _m) real part, represent complex-exponential function exp (-i2 π τ _kf _m) imaginary part; τ _krepresent the time location of a kth reflection coefficient pulse; f _mrepresent m sample frequency; α represents the amplitude vector of reflection coefficient pulse, α=[a ₁, a ₂..., a _k] ^t, a _krepresent the amplitude of a kth reflection coefficient pulse; G α expressed portion open score generated data; D represents known partial spectrum observation data; || || ²represent the quadratic sum asking for data, then open radical sign; λ represents ratio of damping;

First judge module, for judging whether system capacity J reaches preset range, and when judged result is for being, determines that corresponding reflection coefficient pulse time positions and amplitude are available.

In one embodiment, judging unit 25 can comprise: the 4th computing module and the second judge module.

4th computing module, specifically for adopting following formulae discovery square error J _a:

J _a＝||Gα-d|| ²，

Wherein, $G={\left[\begin{array}{cc}{G}_{mk}^r& G_{mk}^i\end{array}\right]}^T,G_{mk}^r=cos\left(2{\mathrm{\π\τ}}_k{f}_m\right)$ Represent complex-exponential function exp (-i2 π τ _kf _m) real part, represent complex-exponential function exp (-i2 π τ _kf _m) imaginary part; τ _krepresent the time location of a kth reflection coefficient pulse; f _mrepresent m sample frequency; α represents the amplitude vector of reflection coefficient pulse, α=[a ₁, a ₂..., a _k] ^t, a _krepresent the amplitude of a kth reflection coefficient pulse; G α expressed portion open score generated data; D represents known partial spectrum observation data; || || ²represent the quadratic sum asking for data, then open radical sign;

Second judge module, for judging square error J _awhether reach preset threshold range, and when judged result is for being, determine iteration convergence.

In one embodiment, reconfiguration unit 26 specifically for: reflection coefficient pulse number, available reflection coefficient pulse time positions and amplitude are brought in non-zero reflection coefficient sequence expression formula, the expression formula of non-zero reflection coefficient sequence r (t) is as follows:

$r\left(t\right)=\underset{k=1}{\overset{K}{\Σ}}{a}_k\mathrm{\δ}(t-{\mathrm{\τ}}_k),$

Wherein, K represents the reflection coefficient pulse number of non-zero; τ _krepresent the time location of a kth reflection coefficient pulse; a _krepresent the amplitude of a kth reflection coefficient pulse; δ (t) indicating impulse function.

In one embodiment, method also comprises: in an iterative process, if the maximum iteration time that iterations reaches default detected, then stops iteration.

Certainly, above-mentioned Module Division just a kind of signal divides, and the present invention is not limited thereto.As long as the Module Division of object of the present invention can be realized, protection scope of the present invention all should be belonged to.

In order to above-mentioned global optimization time thinning thin the Method of Deconvolution and device more clearly explain, be described below in conjunction with specific embodiment and accompanying drawing, but it should be noted that this embodiment is only to better the present invention is described, do not form and the present invention is limited improperly.

Fig. 3 is the schematic diagram of the true reflection coefficient sequence model of the embodiment of the present invention, and in figure 3, abscissa representing time, ordinate represents amplitude (or being called amplitude).As shown in Figure 3, reflection coefficient is pulse.Fig. 4 A is the waveform schematic diagram of the time-varying wavelet that the embodiment of the present invention adopts, Fig. 4 B be the embodiment of the present invention Fig. 4 A shown in spectral amplitude corresponding to time-varying wavelet.In Fig. 4 A and Fig. 4 B, solid line represents the superficial part of time-varying wavelet, and dotted line represents the deep of time-varying wavelet.Fig. 5 is the seismologic record data schematic diagram with decay characteristics of embodiment of the present invention synthesis.

Fig. 6 is the comparison diagram of the reflection coefficient that the reflection coefficient that obtains of the present invention and conventional method obtain.As shown in Figure 6, solid dot represents the reflection coefficient sequence (namely adopting method of the present invention) adopting time-varying wavelet inverting to obtain, asterisk represents the reflection coefficient sequence (namely adopting conventional method) adopting stable state wavelet to obtain, and solid line represents true reflection coefficient sequence model.Can find out, when not adopting time-varying wavelet, reflection coefficient position and size all create very large error, influenced little compared with deep when superficial part (time is less), this is because the wavelet that superficial part energy obtains by force meets superficial part seismic data feature more, and the reflection coefficient in deep cannot accurately recover.

Fig. 7 A is the waveform schematic diagram of time-varying wavelet that the actual seismic data of the embodiment of the present invention obtains, Fig. 7 B be the embodiment of the present invention Fig. 7 A shown in spectral amplitude corresponding to time-varying wavelet.Fig. 8 is the design sketch that actual seismic data is obtained by conventional method, and Fig. 9 is the design sketch that actual seismic data is obtained by the inventive method.By comparing conventional method deconvolution result and global optimization time-varying wavelet deconvolution result can clearly be found out, result has had very large improvement.The data division of superficial part (time is less) concentration of energy, be adopt stable state wavelet or time-varying wavelet all obtains good effect, reflecting interface is clear, and continuity is good, demonstrates stability and the accuracy of method; Deep energy is more weak, the data region that wavelet attenuation degree is large, it is undesirable that employing stable state wavelet carries out deconvolution effect, most reflecting interface is not all identified, and lineups are interrupted, poor continuity, and the sparse spectrum the Method of Deconvolution of the time-varying wavelet adopting the present invention to propose, result has clear improvement, and lineups are clear, reflecting interface continous-stable, and structural attitude obtains extraordinary embodiment, for subsequent treatment and explanation provide good data base, also to sequence dividing, petroleum-gas prediction is very helpful.

In sum, when Songliao basin improves day by day, utilizing available data to realize the maximum using of the hydrocarbon resources developed or be about to exploitation is the key issue that geophysical method needs to solve, address this problem and need more accurate simulation subsurface features, the result that choosing method is higher with obtaining resolution more targetedly.Embodiments provide a kind of time thinning thin the Method of Deconvolution and device of global optimization, it is frequency field time variant deconvolution, time-varying wavelet is adopted more coincidently to descend actual conditions, utilize the non-linear means combined with linear method can ask for high resolving power reflection coefficient sequence, to reaching the object improving non-stationary data resolution.The present invention is based on global optimization approach and ask for reflection coefficient pulse position, asked for the amplitude of reflection coefficient pulse simultaneously by least-squares algorithm, in frequency field, adopt the accurate estimation of time-varying wavelet realization to reflection coefficient sequence.

Describe and can be understood in process flow diagram or in this any process otherwise described or method, represent and comprise one or more for realizing the module of the code of the executable instruction of the step of specific logical function or process, fragment or part, and the scope of the preferred embodiment of the present invention comprises other realization, wherein can not according to order that is shown or that discuss, comprise according to involved function by the mode while of basic or by contrary order, carry out n-back test, this should understand by embodiments of the invention person of ordinary skill in the field.

In the description of this instructions, specific features, structure, material or feature that the description of reference term " embodiment ", " some embodiments ", " example ", " concrete example " or " some examples " etc. means to describe in conjunction with this embodiment or example are contained at least one embodiment of the present invention or example.In this manual, identical embodiment or example are not necessarily referred to the schematic representation of above-mentioned term.And the specific features of description, structure, material or feature can combine in an appropriate manner in any one or more embodiment or example.

Above-described specific embodiment; object of the present invention, technical scheme and beneficial effect are further described; be understood that; the foregoing is only specific embodiments of the invention; the protection domain be not intended to limit the present invention; within the spirit and principles in the present invention all, any amendment made, equivalent replacement, improvement etc., all should be included within protection scope of the present invention.

Claims (10)

1. a time thinning thin the Method of Deconvolution for global optimization, is characterized in that, comprising:

The partial spectrum deconvolution model become when building the frequency field of unstable state seismic data;

Utilize the effective bandwidth determination reflection coefficient pulse number of limited turnover rate principle and described unstable state seismic data;

Utilize very fast simulated annealing algorithm, in conjunction with described partial spectrum deconvolution model, the time location of reflection coefficient pulse is searched for, for each time location searched, utilize least square method to carry out iteration, determine the reflection coefficient pulse amplitude at each time location place;

For each iteration, determine available reflection coefficient pulse time positions and amplitude according to system capacity;

Utilize the square error of partial spectrum observation data and partial spectrum generated data, judge whether iteration restrains;

If convergence, then stop iteration, according to described reflection coefficient pulse number, described available reflection coefficient pulse time positions and amplitude, reconstruct non-zero reflection coefficient sequence.

2. method according to claim 1, is characterized in that, the partial spectrum deconvolution model become when building the frequency field of unstable state seismic data, comprising:

When described unstable state seismic data extracts, become seismic wavelet, the expression becoming seismic wavelet time described is as follows: wherein, become the phase spectrum of seismic wavelet during expression, become the spectral amplitude of seismic wavelet when A (ω) represents, HT represents Hilbert transform, and ω represents angular frequency;

Build described partial spectrum deconvolution model in conjunction with becoming seismic wavelet time described, described partial spectrum deconvolution model is as follows:

$S\left(f_m\right)/W\left(f_m\right)=\underset{k=1}{\overset{K}{\Σ}}{a}_k\exp (-i2{\mathrm{\π\τ}}_k{f}_m)+N\left(f_m\right)/W\left(f_m\right),$

Wherein, S (f _m) represent that the Fourier transform of unstable state seismic data is in frequency f _mthe value at place; M=1,2 ..., M, M>=K; M represents at sampling frequency domain (f _low, f _high) frequency points of up-sampling; K represents the reflection coefficient pulse number of non-zero; a _krepresent the amplitude of a kth reflection coefficient pulse; τ _krepresent the time location of a kth reflection coefficient pulse; I represents imaginary unit; N (f _m) represent that the Fourier transform of noise is in frequency f _mthe value at place; W (f _m) become the Fourier transform of seismic wavelet in frequency f when representing _mthe value at place.

3. method according to claim 1, is characterized in that, utilizes the effective bandwidth determination reflection coefficient pulse number of limited turnover rate principle and described unstable state seismic data, comprising:

Effective bandwidth is (f _low, f _high) the unstable state seismic data non-zero reflection coefficient pulse number K that can reconstruct be K=[| f _high-f _low| × L/2] _int, wherein, L represents the cycle of non-zero reflection coefficient pulse, [] _intrepresent and round downwards.

4. method according to claim 1, is characterized in that, utilizes least square method to carry out iteration, determines the reflection coefficient pulse amplitude at each time location place, comprising:

If noiseless and the prior-constrained K of the sparse number of plies is less than the first preset value, adopt the amplitude of following formulae discovery reflection coefficient pulse:

$\widehat{\mathrm{\α}}={\left(G^{T}G\right)}^{-1}{G}^{T}d,$

Wherein, represent the reflection coefficient pulse amplitude vector obtained; $G={\left[\begin{array}{cc}{G}_{mk}^r& G_{mk}^i\end{array}\right]}^T,$ represent complex-exponential function exp (-i2 π τ _kf _m) real part, represent complex-exponential function exp (-i2 π τ _kf _m) imaginary part; τ _krepresent the time location of a kth reflection coefficient pulse; f _mrepresent m sample frequency; D represents known partial spectrum observation data; T representing matrix transposition;

If noise is stronger and the prior-constrained K of the sparse number of plies is greater than the second preset value, adopt the amplitude of following formulae discovery reflection coefficient pulse:

$\widehat{\mathrm{\α}}={(G^{T}G+\mathrm{\λ}I)}^{-1}{G}^{T}d,$

Wherein, λ represents ratio of damping; I representation unit matrix.

5. method according to claim 1, is characterized in that, for each iteration, determines available reflection coefficient pulse time positions and amplitude, comprising according to system capacity:

Adopt system capacity J described in following formulae discovery:

J＝||Gα-d|| ²+λ||α|| ²，

Wherein, $G={\left[\begin{array}{cc}{G}_{mk}^r& G_{mk}^i\end{array}\right]}^T,$ represent complex-exponential function exp (-i2 π τ _kf _m) real part, represent complex-exponential function exp (-i2 π τ _kf _m) imaginary part; τ _krepresent the time location of a kth reflection coefficient pulse; f _mrepresent m sample frequency; α represents the amplitude vector of reflection coefficient pulse, α=[a ₁, a ₂..., a _k] ^t, a _krepresent the amplitude of a kth reflection coefficient pulse; G α expressed portion open score generated data; D represents known partial spectrum observation data; || || ²represent the quadratic sum asking for data, then open radical sign; λ represents ratio of damping;

If described system capacity J reaches preset range, then determine that corresponding reflection coefficient pulse time positions and amplitude are available.

6. method according to claim 1, is characterized in that, utilizes the square error of partial spectrum observation data and partial spectrum generated data, judges whether iteration restrains, and comprising:

Adopt square error J described in following formulae discovery _a:

J _a＝||Gα-d|| ²，

Wherein, $G={\left[\begin{array}{cc}{G}_{mk}^r& G_{mk}^i\end{array}\right]}^T,$ represent complex-exponential function exp (-i2 π τ _kf _m) real part, represent complex-exponential function exp (-i2 π τ _kf _m) imaginary part; τ _krepresent the time location of a kth reflection coefficient pulse; f _mrepresent m sample frequency; α represents the amplitude vector of reflection coefficient pulse, α=[a ₁, a ₂..., a _k] ^t, a _krepresent the amplitude of a kth reflection coefficient pulse; G α expressed portion open score generated data; D represents known partial spectrum observation data; || || ²represent the quadratic sum asking for data, then open radical sign;

If described square error J _areach preset threshold range, then determine iteration convergence.

7. method according to claim 1, is characterized in that, according to described reflection coefficient pulse number, described available reflection coefficient pulse time positions and amplitude, reconstruct non-zero reflection coefficient sequence, comprising:

Described reflection coefficient pulse number, described available reflection coefficient pulse time positions and amplitude are brought in non-zero reflection coefficient sequence expression formula, the expression formula of described non-zero reflection coefficient sequence r (t) is as follows:

$r\left(t\right)=\underset{k=1}{\overset{K}{\Σ}}{a}_k\mathrm{\δ}(t-{\mathrm{\τ}}_k),$

Wherein, K represents the reflection coefficient pulse number of non-zero; τ _krepresent the time location of a kth reflection coefficient pulse; a _krepresent the amplitude of a kth reflection coefficient pulse; δ (t) indicating impulse function.

8. method according to any one of claim 1 to 7, is characterized in that, described method also comprises: in an iterative process, if the maximum iteration time that iterations reaches default detected, then stops iteration.

9. global optimization time become a sparse deconvolution device, it is characterized in that, comprising:

Construction unit, the partial spectrum deconvolution model become during for building the frequency field of unstable state seismic data;

First determining unit, for utilizing the effective bandwidth determination reflection coefficient pulse number of limited turnover rate principle and described unstable state seismic data;

Iteration unit, for utilizing very fast simulated annealing algorithm, in conjunction with described partial spectrum deconvolution model, the time location of reflection coefficient pulse is searched for, for each time location searched, utilize least square method to carry out iteration, determine the reflection coefficient pulse amplitude at each time location place;

Second determining unit, for for each iteration, determines available reflection coefficient pulse time positions and amplitude according to system capacity;

Judging unit, for utilizing the square error of partial spectrum observation data and partial spectrum generated data, judges whether iteration restrains;

Reconfiguration unit, for when iteration convergence, stops iteration, according to described reflection coefficient pulse number, described available reflection coefficient pulse time positions and amplitude, and reconstruct non-zero reflection coefficient sequence.

10. device according to claim 9, is characterized in that, described construction unit comprises:

Extraction module, for becoming seismic wavelet when described unstable state seismic data extracts, the expression becoming seismic wavelet time described is as follows: wherein, become the phase spectrum of seismic wavelet during expression, become the spectral amplitude of seismic wavelet when A (ω) represents, HT represents Hilbert transform, and ω represents angular frequency;

Build module, for building described partial spectrum deconvolution model in conjunction with becoming seismic wavelet time described, described partial spectrum deconvolution model is as follows:

$S\left(f_m\right)/W\left(f_m\right)=\underset{k=1}{\overset{K}{\Σ}}{a}_k\exp (-i2{\mathrm{\π\τ}}_k{f}_m)+N\left(f_m\right)/W\left(f_m\right),$

Wherein, S (f _m) represent that the Fourier transform of unstable state seismic data is in frequency f _mthe value at place; M=1,2 ..., M, M>=K; M represents at sampling frequency domain (f _low, f _high) frequency points of up-sampling; K represents the reflection coefficient pulse number of non-zero; a _krepresent the amplitude of a kth reflection coefficient pulse; τ _krepresent the time location of a kth reflection coefficient pulse; I represents imaginary unit; N (f _m) represent that the Fourier transform of noise is in frequency f _mthe value at place; W (f _m) become the Fourier transform of seismic wavelet in frequency f when representing _mthe value at place.