CN105842732A - Inversion method of multichannel sparse reflection coefficient and system thereof - Google Patents

Abstract

The invention provides an inversion method of a multichannel sparse reflection coefficient and a system thereof. The method comprises the following steps of collecting post-stack seismic data; according to the post-stack seismic data, extracting a seismic wavelet; according to the seismic wavelet, generating a seismic wavelet convolution matrix; according to the seismic wavelet convolution matrix and the post-stack seismic data, establishing a multichannel reflection coefficient inversion equation; acquiring an initial value of a multichannel sparse reflection coefficient model parameter which is given in advance; through a massive Bayes learning theory framework, combining the multichannel reflection coefficient inversion equation to update the multichannel sparse reflection coefficient model parameter and finally acquiring the updated multichannel sparse reflection coefficient. A reflection coefficient profile acquired through inversion possesses high resolution and lateral continuity of the profile is maintained.

Description

The inversion method of the sparse reflection coefficient of multiple tracks and system

Technical field

The present invention is about oil-gas exploration and development technique field, especially with regard to the inverting of physical prospecting earthquake model parameter Technology, is concretely inversion method and the system of the sparse reflection coefficient of a kind of multiple tracks.

Background technology

During oil exploration and exploitation, precise seismic interpretation generally requires and effectively disappears during data processes Except the impact of band limit wavelet, the reflection coefficient of renwing stratum, improve the resolution of seismic data.At present, the estimation of reflection coefficient Mainly using inversion method based on convolution model, the quality of its result is by noise, seismic wavelet, the vacation of reflection coefficient And if affecting in terms of the hypothesis these four of convolution model.

Reflection coefficient based on convolution model of the prior art is estimated to limit wavelet, reflection coefficient statistically rule for band Carry out a series of it is assumed that and obtain corresponding reflection coefficient by stable the Method for Numerical Inversion.But, based on spectral factorization and Geology assumes that the thin bed reflection coefficient inverting of constraint belongs to the constraint inversion method of nonlinear function class, and its essence is expectation Solve from convolution formula and obtain the most sparse reflection coefficient pulse train.Recognized further by the priori of reflection coefficient pulse The convergent tendency solved in constraint refutation process, exports in the reflection coefficient result of each seismic channel except a few pulses is nonzero value Outward, other are all zero.Therefore, in terms of thin bed reflection coefficient inverting and tuning thickness thickness of thin layer identified below, such method Accurate inversion result can be obtained.

But, there is following limitation in this technology:

(1), the noise immunity problem of this technology itself.

In the case of seismic data signal to noise ratio is relatively low, improves inverting in an iterative process by the constraint of prior information and tie The ability of fruit is restricted.Although the reflection coefficient that final inverting obtains maintains openness, but exists bigger with true model Difference.

(2), the holding of lateral continuity.

Traditional is to reach by calculating single seismic channel overall situation sparse solution based on sparse constraint reflection coefficient inversion method To inverting reflection coefficient purpose, do not consider reflection coefficient seriality on adjacent seismic channel, and reflection coefficient sequence exists Actual formation has certain horizontal dependency, when carrying out single track reflection coefficient inverting, does not consider whether result protects Hold this character.When noise level is higher or reflection coefficient amplitude is more weak, traditional single track reflection coefficient inverting cannot root Accurately recover position and the amplitude of weak signal according to existing information, thus cause whole inversion result that deviation occurs.

Therefore, how to research and develop out a kind of new scheme, it can make full use of reflecting interface practically lowerly Plastid has the feature of certain lateral extension, reaches accurately to estimate, reflection coefficient position and amplitude especially in weak signal Identification aspect, and to be finally inversed by reflection coefficient in the section of low signal-to-noise ratio be this area technical barrier urgently to be resolved hurrily.

Summary of the invention

For the above-mentioned technical problem overcoming prior art to exist, the invention provides the sparse reflection coefficient of a kind of multiple tracks Inversion method and system, with poststack seismic data for input data, under block Bayesian Learning Theory framework, by multiple Seismic channel reflection coefficient sparse prior assume and seismic channel reflection coefficient between dependency it is assumed that utilize MAP estimation and EM algorithm solves the Sparse Pulse reflection coefficient of multiple tracks simultaneously, and final inverting obtains reflection coefficient section, has superelevation Resolution, and maintain the lateral continuity of section.

It is an object of the invention to provide the inversion method of the sparse reflection coefficient of a kind of multiple tracks, described method includes: adopt Collection post-stack seismic data；Seismic wavelet is extracted according to described post-stack seismic data；Earthquake is generated according to described seismic wavelet Wavelet convolution matrix；Multiple tracks reflection coefficient is set up according to described seismic wavelet convolution matrix and described post-stack seismic data Inversion equation；Obtain the initial value of previously given multiple tracks sparse Bayesian reflectivity model parameter；By block Bayes Theory of learning framework combines multiple tracks reflection coefficient inversion equation reflectivity model sparse to described multiple tracks parameter and is updated, Multiple tracks sparse reflectivity model parameter after renewal.

In a preferred embodiment of the invention, described multiple tracks sparse reflectivity model parameter include non-negative parameter, Covariance matrix, noise variance and average.

In a preferred embodiment of the invention, multiple tracks reflection coefficient is combined by block Bayesian Learning Theory framework anti- Drill equation reflectivity model sparse to described multiple tracks parameter to be updated, the sparse reflectivity model of multiple tracks after being updated Parameter includes: obtain maximum iteration time and model modification tolerable error；According to the previously given sparse reflection coefficient of multiple tracks The initial value of model parameter combines multiple tracks reflection coefficient inversion equation and determines that renewal is many by block Bayesian Learning Theory framework Road sparse reflectivity model parameter；The average of renewal is obtained out from the multiple tracks sparse reflectivity model parameter updated；From The initial value of previously given multiple tracks sparse reflectivity model parameter obtains out the initial value of average；Determine described renewal The difference of the initial value of average and average；Whether the difference described in judgement is less than model modification tolerable error；When being judged as NO, Initial by the previously given multiple tracks sparse reflectivity model parameter of multiple tracks sparse reflectivity model parameter iteration updated Value, after record iterations, returns and performs to pass through according to the initial value of previously given multiple tracks sparse reflectivity model parameter Bayesian Learning Theory framework determines the multiple tracks sparse reflectivity model parameter of renewal.

In a preferred embodiment of the invention, described multiple tracks reflection system is combined by block Bayesian Learning Theory framework Number inversion equation reflectivity model sparse to described multiple tracks parameter is updated, the sparse reflection coefficient of multiple tracks after being updated Model parameter also includes: when described difference is less than or equal to model modification tolerable error, the multiple tracks that output updates is sparse instead Penetrating Modulus Model parameter, the multiple tracks that the average in the multiple tracks sparse reflectivity model parameter of described renewal is after inverting is sparse Reflection coefficient.

In a preferred embodiment of the invention, described multiple tracks reflection system is combined by block Bayesian Learning Theory framework Number inversion equation reflectivity model sparse to described multiple tracks parameter is updated, the sparse reflection coefficient of multiple tracks after being updated Model parameter also includes: whether the iterations described in judgement is less than described maximum iteration time；When being judged as NO, output The multiple tracks sparse reflectivity model parameter updated, the average in the multiple tracks sparse reflectivity model parameter of described renewal is The sparse reflection coefficient of multiple tracks after inverting.

In a preferred embodiment of the invention, described method also includes: by MAP estimation and maximum phase Multiple tracks sparse reflectivity model parameter after hoping algorithm estimate described renewal.

The Inversion System of the sparse reflection coefficient of a kind of multiple tracks, described system bag are it is an object of the invention to provide Include post-stack seismic data harvester, be used for gathering post-stack seismic data；Seismic wavelet extraction device, for folding according to described Rear geological data extracts seismic wavelet；Convolution matrix generation device, for generating seismic wavelet pleat according to described seismic wavelet Product matrix；Inversion equation sets up device, for according to described seismic wavelet convolution matrix and described post-stack seismic data Set up multiple tracks reflection coefficient inversion equation；Initial value acquisition device, for obtaining previously given multiple tracks sparse reflection coefficient mould The initial value of shape parameter；Inverting device, for combining multiple tracks reflection coefficient inverting side by block Bayesian Learning Theory framework Journey reflectivity model sparse to described multiple tracks parameter is updated, the multiple tracks sparse reflectivity model ginseng after being updated Number.

In a preferred embodiment of the invention, described multiple tracks sparse reflectivity model parameter include non-negative parameter, Covariance matrix, noise variance and average.

In a preferred embodiment of the invention, described inverting device includes: acquisition module, is used for obtaining greatest iteration Number of times and model modification tolerable error；Update coefficient determination module, for according to the previously given sparse reflection coefficient of multiple tracks The initial value of model parameter combines multiple tracks reflection coefficient inversion equation and determines renewal by block Bayesian Learning Theory framework Multiple tracks sparse reflectivity model parameter；Update average acquisition module, for from the multiple tracks sparse reflectivity model ginseng updated Number obtains out the average of renewal；Initial value acquisition module, for from previously given multiple tracks sparse reflectivity model parameter Initial value in obtain out the initial value of average；Difference determines module, for determining the average of described renewal and the initial of average The difference of value；First judge module, for judging that whether described difference is less than described model modification tolerable error；Iteration mould Block, for when the first described judge module is judged as NO, pre-with the multiple tracks sparse reflectivity model parameter iteration updated The initial value of first given multiple tracks sparse reflectivity model parameter, after record iterations, returns the renewal system described in performing Number determines module.

In a preferred embodiment of the invention, described inverting device also includes: the first reflection coefficient output module, uses In time being judged as YES when the first described judge module, the multiple tracks sparse reflectivity model parameter that output updates, described renewal Multiple tracks sparse reflectivity model parameter in average be the sparse reflection coefficient of multiple tracks after inverting.

In a preferred embodiment of the invention, described inverting device also includes: the second judge module, is used for judging institute Whether the iterations stated is less than described maximum iteration time；Second reflection coefficient output module, for when described second When judge module is judged as NO, the multiple tracks sparse reflectivity model parameter that output updates, the sparse reflection of multiple tracks of described renewal Average in Modulus Model parameter is the sparse reflection coefficient of the multiple tracks after inverting.

In a preferred embodiment of the invention, described system also includes: reflectivity model parameter estimation apparatus, uses In estimating the multiple tracks sparse reflectivity model parameter after described renewal by MAP estimation and EM algorithm.

In a preferred embodiment of the invention, described renewal coefficient determination module is carried out by equation below:

$\begin{array}{c}{x}^{*}\stackrel{\mathrm{\Δ}}{=}{\mathrm{\μ}}_x={(\mathrm{\λ}{\Σ}_0^{-1}+D^{T}D)}^{-1}{D}^{T}y\\ ={\Σ}_0{D}^T{(\mathrm{\λ}I+D{\Σ}_0{D}^T)}^{-1}y\end{array};$

$\begin{array}{cc}{\mathrm{\γ}}_i\←\frac{1}{L}{X}_i{B}^{-1}{X}_i^T+{\left({\mathrm{\Ξ}}_x\right)}_{ii}& {\∀}_i\end{array};$

$B\←\left(\frac{1}{M}\underset{i=1}{\overset{M}{\Σ}}\frac{{\left({\mathrm{\Ξ}}_x\right)}_{ii}}{{\mathrm{\γ}}_i}\right)B+\frac{1}{M}\underset{i=1}{\overset{M}{\Σ}}\frac{X_i^T{X}_i}{r_i};$

$D=A\⊗I_L$

${\mathrm{\Ξ}}_x={({\mathrm{\Γ}}^{-1}+\frac{1}{\mathrm{\λ}}{A}^{T}A)}^{-1}$

$\mathrm{\Γ}\stackrel{\mathrm{\Δ}}{=}diag({\mathrm{\γ}}_1,...,{\mathrm{\γ}}_M)$

Wherein, x^*For reflection coefficient inversion result, μ_xFor average, λ is noise variance, and I is unit battle array, γ_iJoin for non-negative Number, L is the number of channels of multiple tracks reflection coefficient, and B is covariance matrix, and X is multiple tracks reflection coefficient sequence,Represent and any i is taken Value, N is that earthquake record time sampling is counted, and M is that the time sampling of reflectivity model is counted, and Y is multichannel seismic data, and A is The convolution matrix of seismic wavelet,Represent vector not Luo Beini this norm computing of crow.

The beneficial effects of the present invention is, it is provided that a kind of multiple tracks sparse reflection coefficient inversion method and system, with poststack Seismic data is input data, under block Bayesian Learning Theory framework, by multiple sparse elder generations of seismic channel reflection coefficient Test the dependency between hypothesis and seismic channel reflection coefficient it is assumed that utilize MAP estimation and EM algorithm to solve simultaneously The model parameter of the Sparse Pulse reflection coefficient of multiple tracks, what this technology obtained is reflection coefficient section, has the resolution of superelevation, And maintaining the lateral continuity of section, its computational efficiency is higher, for identifying thin layer, the weak reflective information of renwing stratum, improves Seismic data resolution provides technical support.

For the above and other objects, features and advantages of the present invention can be become apparent, preferred embodiment cited below particularly, And coordinate institute's accompanying drawings, it is described in detail below.

Accompanying drawing explanation

In order to be illustrated more clearly that the embodiment of the present invention or technical scheme of the prior art, below will be to embodiment or existing In having technology to describe, the required accompanying drawing used is briefly described, it should be apparent that, the accompanying drawing in describing below is only this Some embodiments of invention, for those of ordinary skill in the art, on the premise of not paying creative work, it is also possible to Other accompanying drawing is obtained according to these accompanying drawings.

The stream of the embodiment one of the inversion method of the sparse reflection coefficient of a kind of multiple tracks that Fig. 1 provides for the embodiment of the present invention Cheng Tu；

Fig. 2 is the flow chart of the embodiment one of step S106 in Fig. 1；

Fig. 3 is the flow chart of the embodiment two of step S106 in Fig. 1；

Fig. 4 is the flow chart of the embodiment three of step S106 in Fig. 1；

The stream of the embodiment two of the inversion method of the sparse reflection coefficient of a kind of multiple tracks that Fig. 5 provides for the embodiment of the present invention Cheng Tu；

The knot of the embodiment one of the Inversion System of the sparse reflection coefficient of a kind of multiple tracks that Fig. 6 provides for the embodiment of the present invention Structure block diagram；

The enforcement of inverting device in the Inversion System of the sparse reflection coefficient of a kind of multiple tracks that Fig. 7 provides for the embodiment of the present invention The structured flowchart of mode one；

The enforcement of inverting device in the Inversion System of the sparse reflection coefficient of a kind of multiple tracks that Fig. 8 provides for the embodiment of the present invention The structured flowchart of mode two；

The enforcement of inverting device in the Inversion System of the sparse reflection coefficient of a kind of multiple tracks that Fig. 9 provides for the embodiment of the present invention The structured flowchart of mode three；

The embodiment two of the Inversion System of the sparse reflection coefficient of a kind of multiple tracks that Figure 10 provides for the embodiment of the present invention Structured flowchart；

The inverting flow chart of the sparse reflection coefficient of multiple tracks in the specific embodiment that Figure 11 provides for the present invention；

Figure 12 is the model schematic of true stratiform reflection coefficient in specific embodiment；

Figure 13 is the schematic diagram of the seismic wavelet selecting mixed phase wavelet in specific embodiment；

Figure 14 be in specific embodiment by the convolution operation of seismic wavelet and reflection coefficient generate without making an uproar earthquake record Schematic diagram；

Figure 15 is multiple tracks sparse reflection coefficient inversion result schematic diagram in specific embodiment；

Figure 16 is single track sparse Bayesian reflection coefficient inversion result schematic diagram；

Figure 17 is inversion result and the root-mean-square error schematic diagram of true model in specific embodiment；

Figure 18 is noisy composite traces schematic diagram in specific embodiment；

Figure 19 is multiple tracks sparse reflection coefficient inversion result schematic diagram；

Figure 20 is single track sparse Bayesian reflection coefficient inversion result schematic diagram；

Figure 21 inversion result and the root-mean-square error schematic diagram of true model；

Figure 22 is gas field, Sichuan province actual poststack seismic data schematic diagram；

Figure 23 is multiple tracks sparse reflection coefficient inversion result schematic diagram；

Figure 24 is the spectral contrast schematic diagram of raw data and inversion result.

Detailed description of the invention

Below in conjunction with the accompanying drawing in the embodiment of the present invention, the technical scheme in the embodiment of the present invention is carried out clear, complete Describe, it is clear that described embodiment is only a part of embodiment of the present invention rather than whole embodiments wholely.Based on Embodiment in the present invention, it is every other that those of ordinary skill in the art are obtained under not making creative work premise Embodiment, broadly falls into the scope of protection of the invention.

The invention belongs to geophysical prospecting for oil earthquake model parametric inversion field.Particularly relate to based on convolution model The implementation method of nonlinear reflection coefficient high-resolution inversion.

There is to make full use of reflecting interface in actual underground geologic bodies the feature of certain lateral extension, reach instead Penetrate coefficient positions and amplitude is accurately estimated, especially the identification aspect of weak signal, and be finally inversed by the section of low signal-to-noise ratio Reflection coefficient, the present invention proposes a kind of multiple tracks reflection coefficient inversion technique based on block Bayesian Learning Theory framework of knowing clearly.

This technology is with poststack seismic data for input data, under block Bayesian Learning Theory framework, by multiple Seismic channel reflection coefficient sparse prior assume and seismic channel reflection coefficient between dependency it is assumed that utilize MAP estimation and EM algorithm solves the Sparse Pulse reflection coefficient of multiple tracks simultaneously.What this technology obtained is reflection coefficient section, has super High resolution, and maintain the lateral continuity of section.Its computational efficiency is higher, and for identifying thin layer, renwing stratum is weak instead Penetrate information, improve seismic data resolution and provide technical support.

The stream of the embodiment one of the inversion method of the sparse reflection coefficient of a kind of multiple tracks that Fig. 1 provides for the embodiment of the present invention Cheng Tu, as shown in Figure 1, in embodiment one, the method includes:

S101: gathering post-stack seismic data, in a particular embodiment, post-stack seismic data such as can use x, and (t n) comes Representing, wherein, x represents multitrace seismogram value, t express time, and n represents Taoist monastic name.

S102: extracting seismic wavelet according to described post-stack seismic data, in a particular embodiment, seismic wavelet is all As available w (t) represents, wherein, t express time, w represents seismic wavelet.

S103: generate seismic wavelet convolution matrix, in a particular embodiment, earthquake according to described seismic wavelet Ripple convolution matrix such as can represent with A；

S104: set up multiple tracks reflection coefficient according to described seismic wavelet convolution matrix and described post-stack seismic data Inversion equation, in a particular embodiment, multitrace seismogram convolution model equation such as can represent with Y=AX+V.Its Middle Y represents that multitrace seismogram, A represent seismic wavelet convolution matrix, and X represents multiple tracks reflection coefficient sequence, and V represents that multiple tracks is observed Noise.

S105: obtain the initial value of previously given multiple tracks sparse reflectivity model parameter.In the present invention, described Multiple tracks sparse reflectivity model parameter includes non-negative parameter γ_i, covariance matrix B, noise variance λ and mean μ_x, i represents The index of time sampling.Its each self-corresponding initial value can be expressed as γ⁰、B⁰、λ⁰、

S106: combine multiple tracks reflection coefficient inversion equation by block Bayesian Learning Theory framework sparse to described multiple tracks Reflectivity model parameter is updated, the multiple tracks sparse reflectivity model parameter after being updated.

Fig. 2 is the flow chart of the embodiment one of step S106, as shown in Figure 2, in embodiment one, and step S106 bag Include:

S201: obtain maximum iteration time and model modification tolerable error, in a particular embodiment, greatest iteration Number of times such as can use it_maxRepresenting, model modification tolerable error such as can represent with E.

S202: combine multiple tracks reflection coefficient according to the initial value of previously given multiple tracks sparse reflectivity model parameter anti- Drill equation and determine renewal multiple tracks sparse reflectivity model parameter by block Bayesian Learning Theory framework.Introduce first below Principle.

1, reflection coefficient inverting ultimate principle based on convolution model

In field of seismic exploration, widely used convolution model describes the seismic wave process at underground propagation, its concrete table Reaching formula is:

S (t)=w (t) * r (t)+n (t) (1)

Wherein, t express time, s (t) represents earthquake record, and r (t) represents reflection coefficient, and w (t) represents seismic wavelet, n T () represents seismic noise, * represents convolution operation.In regular hour scope and spatial dimension, seismic wavelet is considered stable Constant.Therefore, multichannel deconvolution model can be written as:

S=WR+N (2)

Wherein, S represents that multitrace seismogram, W represent the convolution matrix of stable state seismic wavelet, and R represents multiple tracks reflection coefficient, N is noise.Owing to S is affected by W and N, especially seismic wavelet continuity on time orientation, causes directly using S Resolution reservoir is restricted, particularly the reservoir of micro-structure.Process of seismic data processing is eliminated by various mathematical methods Noise and the seismic wavelet impact on seismic data, it is therefore an objective to recover the reflection coefficient sequence R of subsurface formations.Solve reflection coefficient Most common method is predictive deconvolution.The method assumes that reflection coefficient sequence obeys white noise distribution, and seismic wavelet is Little phase place wavelet.The auto-correlation of earthquake record can be used to replace the auto-correlation of seismic wavelet by this hypothesis.At autocorrelation domain Neutron deficiency one is eliminated.When two kinds of assumed conditions are disagreed with practical situation, it was predicted that deconvolution result will be affected. Wiggins proposes minimum entropy deconvolution method, and the method introduces the variance mould of a linear operator and seismic data convolution result Maximum.Ask for the optimum linear operator maximum so that object function by the method for linear iteraction, finally this operator is acted on Seismic data obtains the reflection coefficient of simplest structure.But, owing to using linear system, the method can not suitably process band Limit data.When processing noisy data, result is the most unstable.

Process the limitation of data to solve linear system, some mathematical functions, such as entropy norm and Cauchy criterion etc. Control the sparse non-linear inversion of reflection coefficient or the Method of Deconvolution is suggested.These methods pass through material matches and a non-thread Property function between mutually constraint widen the frequency band of grandfather tape limit data.By nonlinear function as prior-constrained reflection Coefficient inversion method can obtain the sparse reflection coefficient inversion result with little ring.Compared to material matches item, non-thread When the weight of property function is relatively large, ringing will reduce, but inversion result is unsatisfactory for formula (1).Therefore material matches Relation balance between item and nonlinear function item is the key that this kind of method is successful.

Assume that the thin bed reflection coefficient inversion method of constraint can also regard nonlinear system as based on spectral factorization and geology Inverting, and be all that expectation obtains the most sparse reflection coefficient from convolution model.But such reflection coefficient inverting side at present Method there is also certain problem, mainly has and includes: the noise immunity of method itself, and signal to noise ratio is lower ground shaken data and will be reduced instead Penetrate the accuracy of coefficient inverting；Refutation process needs have certain priori to recognize reflection coefficient pulse number, particularly needs Pulse number to be determined in advance；Inversion method does not accounts for the inaccurate impact on result of seismic wavelet；Computational efficiency is tied with inverting The precision of fruit is not reaching to preferable compromise and processes；Certain lateral continuity is there is in reflection coefficient in actual formation, and this Class method does not accounts for keeping the geometric relativity on inversion result direction in space.

2, multiple tracks sparse reflection coefficient inversion method

In order to from poststack seismic data exactly inverting obtain reflection coefficient section, keep the horizontal of reflection coefficient section Seriality, and remain able to obtain relatively accurate inversion result in the case of noise level is higher.Therefore, with formula (2) implementation method of the present invention that derives based on.

Assuming that seismic wavelet is stable within regular hour and spatial dimension, reflectivity model meets openness vacation If condition, and spatially there is certain dependency.Accurate method is used to estimate seismic wavelet, structure by seismic profile Building the inversion equation about reflection coefficient, being write as matrix form has:

Y=AX+V (3)

Wherein,For multichannel seismic data, the row of Y represents that per pass geological data time sampling is counted, and list shows ground Shake number of channels, i.e. Convolution matrix for seismic wavelet； For multiple tracks reflectivity model；V represents noise, and N is that earthquake record time sampling is counted, and L represents the number of channels of multiple tracks reflection coefficient. M is that the time sampling of reflectivity model is counted.

Within regular hour and spatial dimension, it is believed that the propagation of seismic wavelet is stablized constant, and particular moment Reflection coefficient between each seismic channel, there is certain dependency.For reflectivity model, openness vacation need to be met If the non-zero line number of condition, i.e. X have to be lower than a certain threshold value, in order to guarantee that inversion method can converge to globally optimal solution.Instead It is identical for penetrating coefficient nonzero value index value in X matrix each column, i.e. non-zero reflection system on identical reflex time position Transversely there is seriality in number.Actual reflectivity model may be the most slowly varying, i.e. there is thin layer.Therefore, jointly Openness assumed condition just for the model of less L.

In order to make full use of reflection coefficient lateral continuity and feature of temporal correlation between adjacent seismic channel, Present invention introduces block management loading framework.In this framework, multiple tracks fractal model can be converted to Block unidirectional amount reflectivity model.Assume to be separate between non-zero reflection coefficient row in X matrix, i.e. different The non-zero reflection coefficient of time location is separate, and the reflection coefficient of identical time location takes in each seismic channel From Gauss distribution, i.e.

Wherein, p () represents probability density function, γ_iIt is non-negative parameter, for controlling the openness of reflection coefficient matrix X. Work as γ_iWhen=0, the i-th row reflectance value of homography X is zero.B_iIt is a positive definite matrix, for Description Matrix X the i-th row Middle reflection coefficient is worth correlation structure, owing to reflection coefficient is unknown quantity, and therefore B_iNeeds are estimated.

OrderV=vec (V^T), by multiple tracks Inversion equation is converted to the Bayesian model of bulk: y=Dx+v.WhereinRepresenting Kronecker product, i.e. tensor product, T represents square Battle array transposition computing, vec representing matrix presses rearrangement computing, I_LFor L rank unit matrix.In X, K non-zero row shows in x and is K Nonzero block.

Each element assuming noise vector v is independent, and has Gauss distribution, i.e. p (υ_i)～Ν (0, λ), λ is The variance of Gauss distribution.Therefore for block model, Gauss likelihood model is set up:

Inversion problem is converted into the sparse linear regression problem of a known primary signal, prior information is given for x:Wherein:

Represent for any i value.

By Bayesian formula, can obtain the priori probability density about x, it meets Gauss distribution equally:

Its average and variance writing:

${\mathrm{\μ}}_x=\frac{1}{\mathrm{\λ}}{\Σ}_x{D}^{T}y---\left(8\right)$

$\begin{array}{c}{\Σ}_x={({\Σ}_0^{-1}+\frac{1}{\mathrm{\λ}}{D}^{T}D)}^T\\ ={\Σ}_0-{\Σ}_0{D}^T{(\mathrm{\λ}I+D{\Σ}_0{D}^T)}^{-1}D{\Σ}_0\end{array}---\left(9\right)$

Given hyper parameter λ, γ_i,B_i,For the MAP estimation of x, can direct the most worth from Posterior estimator Arrive, it may be assumed that

$\begin{array}{c}{x}^{*}\stackrel{\mathrm{\Δ}}{=}{\mathrm{\μ}}_x={(\mathrm{\λ}{\Σ}_0^{-1}+D^{T}D)}^{-1}{D}^{T}y\\ ={\Σ}_0{D}^T{(\mathrm{\λ}I+D{\Σ}_0{D}^T)}^{-1}y\end{array}---\left(10\right)$

For the estimation of hyper parameter, solved by Equations of The Second Kind maximum Likelihood or EM algorithm, But unlike that traditional Bayesian Learning Theory framework, reflect the covariance matrix of same time location reflection coefficient dependency B_iNeed to estimate as hyper parameter equally.During the reflection coefficient covariance matrix difference of different time position, it will make The reflection coefficient over-fitting solved.In order to avoid this problem, the present invention uses covariance matrix (positive definite matrix) B, represents The reflection coefficient of different time points position meets identical correlated characteristic, then formula (6) is reduced to:Whereindiag(γ₁,…γ_M) represent with γ₁,…γ_MIt it is the square formation of main diagonal element.Found by derivation Even if correlated characteristic is different, the present invention also can obtain preferable result, and the characteristic of overall situation sparse solution will not be produced shadow Ring.

Note Θ={ γ₁,…,γ_M, B, λ }, the present invention uses EM algorithm (EM) to estimate corresponding hyper parameter, passes through Corresponding algorithm obtains so that Probability p (y；Θ) maximum parameter value.Sparse x vector is carried out marginalisation integration, is closed In hyper parameter and the marginal likelihood function of matrix, the log expressions of equivalence is as follows:

$\begin{array}{c}L\left(\mathrm{\Θ}\right)\stackrel{\mathrm{\Δ}}{=}-2\log \∫p\left(y\right|x;\mathrm{\λ})p(x;{\mathrm{\γ}}_i,B_i,\∀i)dx\\ =y^T{\left({\Σ}_y\right)}^{-1}y+\log \left|{\Σ}_y\right|\end{array}---\left(11\right)$

Wherein,X is as implicit variable so that following object function reaches maximum:

$\begin{array}{c}Q\left(\mathrm{\Θ}\right)=E_{x|y;{\mathrm{\Θ}}^{\left(old\right)}}\[\log p(y,x;\mathrm{\Θ})\]\\ =E_{x|y;{\mathrm{\Θ}}^{\left(old\right)}}\[\log p(y,x;\mathrm{\λ})\]+E_{x|y;{\mathrm{\Θ}}^{\left(old\right)}}\[\log p(x;{\mathrm{\γ}}_1,\mathrm{...},{\mathrm{\γ}}_M,B)\]\end{array}---\left(12\right)$

In equation, right-hand member Section 1 is unrelated with hyper parameter γ and matrix B, in order to estimate hyper parameterAnd B, will Section 1 is omitted, and by replacement, expectation function derivation makes its null parameter value, obtains about hyper parameter γ's Learning method:

$\begin{array}{cc}{\mathrm{\γ}}_i\←\frac{Tr\[B^{-1}({\Σ}_x^i+{\mathrm{\μ}}_x^i{\left({\mathrm{\μ}}_x^i\right)}^T)\]}{L}& i=1,...,M\end{array}---\left(13\right)$

All elements sum on the wherein mark of Tr representing matrix, i.e. matrix leading diagonal.Represent vector μ_x(i-1) L+1 to iL index value.Equally,In representing matrix, ranks index is the square formation that (i-1) L+1 to iL extracts.

Equally, object function for the gradient of matrix B is:

$\frac{\∂Q}{\∂B}=-\frac{M}{2}{B}^{-1}+\frac{1}{2}\underset{i=1}{\overset{M}{\Σ}}\frac{1}{{\mathrm{\γ}}_i}{B}^{-1}({\Σ}_x^i+{\mathrm{\μ}}_x^i{\left({\mathrm{\μ}}_x^i\right)}^T)B^{-1}---\left(14\right)$

The learning method about matrix B is obtained equal to zero by gradient:

$B\←\frac{1}{M}\underset{i=1}{\overset{M}{\Σ}}\frac{{\Σ}_x^i+{\mathrm{\μ}}_x^i{\left({\mathrm{\μ}}_x^i\right)}^T}{{\mathrm{\γ}}_i}---\left(15\right)$

Finally need the variance of noise is estimated, object function is made corresponding simplification, obtains the letter with λ as variable Number:

$\begin{array}{c}Q\left(\mathrm{\λ}\right)=E_{x|y;{\mathrm{\Θ}}^{\left(old\right)}}\[\log p\left(y\right|x;\mathrm{\λ})\]\\ \∝-\frac{NL}{2}\log \mathrm{\λ}-\frac{1}{2\mathrm{\λ}}{E}_{x|y;{\mathrm{\Θ}}^{\left(old\right)}}\[{||y-Dx||}_2^2\]\\ =-\frac{NL}{2}\log \mathrm{\λ}-\frac{1}{2\mathrm{\λ}}\[{||y-Dx||}_2^2+E_{x|y;{\mathrm{\Θ}}^{\left(old\right)}}\[{||D(x-{\mathrm{\μ}}_x)||}_2^2\]\]\\ =-\frac{NL}{2}\log \mathrm{\λ}-\frac{1}{2\mathrm{\λ}}\[{||y-{\mathrm{D\μ}}_x||}_2^2+Tr\left({\Σ}_x{D}^{T}D\right)\]\\ =-\frac{NL}{2}\log \mathrm{\λ}-\frac{1}{2\mathrm{\λ}}\[{||y-{\mathrm{D\μ}}_x||}_2^2+\widehat{\mathrm{\λ}}Tr\left({\Σ}_x\right({{\Σ}_x}^{-1}-{{\Σ}_0}^{-1}\left)\right)\]\\ =-\frac{NL}{2}\log \mathrm{\λ}-\frac{1}{2\mathrm{\λ}}\[{||y-{\mathrm{D\μ}}_x||}_2^2+\widehat{\mathrm{\λ}}\[ML-Tr\left({\Σ}_x{{\Σ}_0}^{-1}\right)\]\]\end{array}---\left(16\right)$

WhereinRepresent λ value during front an iteration.To formula (16) about λ derivation, making derivative is zero, obtains about λ Learning method:

$\mathrm{\λ}\←\frac{{||y-{\mathrm{D\μ}}_x||}_2^2+\mathrm{\λ}\[ML-Tr\left({\Σ}_x{\Σ}_0^{-1}\right)\]}{NL}---\left(17\right)$

Wherein | | | |₂Represent the l of vector₂Norm computing.

Owing to employing Kronecker product (tensor product) during calculating, overall calculation amount is caused to increase.By introducing Bayesian learning method in many vector measurements, simplifies hyper parameter, positive definite matrix and the learning process of variance, finally gives corresponding Estimation formulas:

${\mathrm{\Ξ}}_x={({\mathrm{\Γ}}^{-1}+\frac{1}{\mathrm{\λ}}{A}^{T}A)}^{-1}---\left(18\right)$

$\begin{array}{cc}{\mathrm{\γ}}_i\←\frac{1}{L}{X}_i{B}^{-1}{X}_i^T+{\left({\mathrm{\Ξ}}_x\right)}_{ii}& {\∀}_i\end{array}---\left(19\right)$

$B\←\left(\frac{1}{M}\underset{i=1}{\overset{M}{\Σ}}\frac{{\left({\mathrm{\Ξ}}_x\right)}_{ii}}{{\mathrm{\γ}}_i}\right)B+\frac{1}{M}\underset{i=1}{\overset{M}{\Σ}}\frac{X_i^T{X}_i}{r_i}---\left(20\right)$

Ξ_xIllustrate that γ_iAbove be Mahalanobis distance, after easily represent and the symbol that introduces.

Wherein,The Frobenius of representing matrix not Luo Beini this norm of crow.Relative to many vector Bayesian learning sides Method, the method uses Mahalanobis Mahalanobis generalised distance to replace L2 norm to carry out the similarity between description vectors.Iteration Owing to being updated positive definite matrix B in learning process, therefore between multiple tracks, reflection coefficient sequence is estimated tool by spatial coherence Constrained effect.Multiple tracks reflection coefficient inversion result is obtained eventually through formula (10).

That is, step S202 can determine the mean μ of renewal by (10)_x, determine that renewal is many by formula (19) Road model hyper parameter γ_i, determined the covariance matrix B of renewal by formula (20), determine renewal by formula (21) Noise variance λ.In a particular embodiment, the sparse reflection coefficient of multiple tracks after renewal can be expressed as γ¹、B1、λ¹、

S203: obtain out the average of renewal from the multiple tracks sparse reflectivity model parameter updated；

S204: obtain out the initial of average from the initial value of previously given multiple tracks sparse reflectivity model parameter Value；

S205: determine the difference of the average of described renewal and the initial value of average；

Whether S206: the difference described in judgement, less than model modification tolerable error, i.e. judges

S207: when being judged as NO, by the previously given multiple tracks of the multiple tracks sparse reflectivity model parameter iteration updated The initial value of sparse reflectivity model parameter, after record iterations, returns and performs step S202, according to previously given many By Bayesian Learning Theory framework, the initial value of road sparse reflectivity model parameter determines that the multiple tracks of renewal is sparse and reflects system Digital-to-analogue shape parameter.In a particular embodiment, iteration can represent with j herein.

Namely in embodiment one, it may be judged whether when jumping out loop iteration condition, when being unsatisfactory forAnd Iterations now is 1, not up to maximum iteration time it_max, return the multiple tracks sparse reflection coefficient mould redefining renewal Shape parameter.

Follow-up when judging whether to jump out loop iteration condition, need to be unsatisfactory forAnd iterations is not up to Maximum iteration time it_maxTime, return the multiple tracks sparse reflectivity model parameter redefining renewal.

Fig. 3 is the flow chart of the embodiment two of step S106, from the figure 3, it may be seen that in embodiment two, described method Including:

S301: obtain maximum iteration time and model modification tolerable error, in a particular embodiment, greatest iteration Number of times such as can use it_maxRepresenting, model modification tolerable error such as can represent with E.

S302: combine multiple tracks reflection coefficient according to the initial value of previously given multiple tracks sparse reflectivity model parameter anti- Drill equation and determine renewal multiple tracks sparse reflectivity model parameter by block Bayesian Learning Theory framework.Formula can be passed through (10) mean μ of renewal is determined_x, multiple tracks model non-negative parameter γ of renewal is determined by formula (19)_i, pass through formula (20) determine the covariance matrix B of renewal, determined the noise variance λ of renewal by formula (21).

S303: obtain out the average of renewal from the multiple tracks sparse reflectivity model parameter updated；

S304: obtain out the initial of average from the initial value of previously given multiple tracks sparse reflectivity model parameter Value；

S305: determine the difference of the average of described renewal and the initial value of average；

Whether S306: the difference described in judgement, less than model modification tolerable error, i.e. judges

S307: when being judged as NO, by the previously given multiple tracks of the multiple tracks sparse reflectivity model parameter iteration updated The initial value of sparse reflectivity model parameter, returns and performs step S302, redefine the sparse reflection coefficient of multiple tracks of renewal Multiple tracks sparse reflectivity model parameter, and record iterations, in a particular embodiment, iteration can carry out table with j herein Show.

S308: when being judged as YES, the multiple tracks sparse reflectivity model parameter that output updates, the multiple tracks of described renewal is dilute Dredge the sparse reflection coefficient of multiple tracks after the average in reflectivity model parameter is inverting.

Namely in embodiment two, it may be judged whether when jumping out loop iteration condition, when meetingTime, jump Go out circulation, the multiple tracks sparse reflectivity model parameter that output updates.

Fig. 4 is the flow chart of the embodiment three of step S106, as shown in Figure 4, in embodiment three, and described method Including:

S401: obtain maximum iteration time and model modification tolerable error, in a particular embodiment, greatest iteration Number of times such as can use it_maxRepresenting, model modification tolerable error such as can represent with E.

S402: combine multiple tracks reflection coefficient according to the initial value of previously given multiple tracks sparse reflectivity model parameter anti- Drill equation and determine renewal multiple tracks sparse reflectivity model parameter by block Bayesian Learning Theory framework.Formula can be passed through (10) mean μ of renewal is determined_x, the multiple tracks model hyper parameter γ of renewal is determined by formula (19)_i, by formula (20) Determine the covariance matrix B of renewal, determined the noise variance λ of renewal by formula (21).

S403: obtain out the average of renewal from the multiple tracks sparse reflectivity model parameter updated；

S404: obtain out the initial of average from the initial value of previously given multiple tracks sparse reflectivity model parameter Value；

S405: determine the difference of the average of described renewal and the initial value of average；

Whether S406: the difference described in judgement, less than model modification tolerable error, i.e. judges

When being judged as NO, the multiple tracks previously given with the multiple tracks sparse reflectivity model parameter iteration updated is sparse instead Penetrate the initial value of Modulus Model parameter, return and perform step S402, redefine the multiple tracks sparse reflectivity model ginseng of renewal Number, and record iterations, in a particular embodiment, iteration can represent with j herein.

S407: whether continue the iterations described in judging less than described maximum iteration time；

S408: when being judged as NO, the multiple tracks sparse reflectivity model parameter that output updates, the multiple tracks of described renewal is dilute Dredge the sparse reflection coefficient of multiple tracks after average is inverting in reflectivity model parameter.

Namely in embodiment three, it may be judged whether when jumping out loop iteration condition, as j < it_maxBe false i.e. iteration time When number is not less than described maximum iteration time, jump out circulation, the multiple tracks sparse reflectivity model parameter that output updates.

A kind of multiple tracks sparse reflection coefficient multiple tracks sparse reflectivity model parameter that Fig. 5 provides for the embodiment of the present invention The flow chart of the embodiment one of inversion method, as shown in Figure 5, in embodiment two, described method also includes:

S501: gathering post-stack seismic data, in a particular embodiment, post-stack seismic data such as can use x, and (t n) comes Represent.

S502: extracting seismic wavelet according to described post-stack seismic data, in a particular embodiment, seismic wavelet is all As available w (t) represents.

S503: generate seismic wavelet convolution matrix according to described seismic wavelet；

S504: set up multiple tracks reflection coefficient according to described seismic wavelet convolution matrix and described post-stack seismic data Inversion equation.

S505: obtain the initial value of previously given multiple tracks sparse reflectivity model parameter.In the present invention, described Multiple tracks sparse reflectivity model parameter includes non-negative parameter γ_i, covariance matrix B, noise variance λ and mean μ_x。

S506: combine multiple tracks reflection coefficient inversion equation by block Bayesian Learning Theory framework sparse to described multiple tracks Reflectivity model parameter is updated, the multiple tracks sparse reflectivity model parameter after being updated.

S507: estimate the sparse reflection coefficient of multiple tracks after described renewal by MAP estimation and EM algorithm Model parameter.

As it has been described above, be the inversion method of a kind of multiple tracks sparse reflectivity model parameter that the present invention provides, be given Multiple tracks based on convolution model sparse reflection coefficient inverting framework, derived based on convolution model time domain Bayes in detail Study reflection coefficient inversion formula, describes poststack single track reflection coefficient inversion method based on Method Using Relevance Vector Machine method, is given The detailed derivation of management loading reflection coefficient inverting based on multiple tracks temporal correlation, passes through Mahalanobis distance characterizes the dependency between multiple tracks reflection coefficient, utilizes this information to reach the reflection of multiple seismic channels and is Numerical digit puts the exact inversion with amplitude.Use the multiple tracks reflection coefficient inverting framework that the present invention provides, select the multiple tracks of the present invention Sparse reflection coefficient inversion equation, inverting has the reflection coefficient seismic channel of dependency.With traditional sparse constraint reflection coefficient Comparison of Inversions, the advantage that this method exists in terms of weak signal identification and noise immunity, it is finally reached subsurface reflective boundary The accurate inverting of corresponding parameter.

A kind of multiple tracks sparse reflection coefficient multiple tracks sparse reflectivity model parameter that Fig. 6 provides for the embodiment of the present invention The structured flowchart of the embodiment one of Inversion System, it will be appreciated from fig. 6 that in embodiment one, the sparse reflectivity model of multiple tracks The Inversion System of parameter includes:

Post-stack seismic data harvester 101, is used for gathering post-stack seismic data, in a particular embodiment, and poststack Geological data such as can use x, and (t, n) represents, wherein, x represents multitrace seismogram, t express time, and n represents Taoist monastic name.

Seismic wavelet extraction device 102, for extracting seismic wavelet according to described post-stack seismic data, in concrete reality Executing in mode, seismic wavelet such as can represent with w (t), and wherein, t express time, w represents seismic wavelet.

Convolution matrix generation device 103, for generating seismic wavelet convolution matrix according to described seismic wavelet, specifically Embodiment in, seismic wavelet convolution matrix such as can represent with A；

Inversion equation sets up device 104, for according to described seismic wavelet convolution matrix and described poststack earthquake Data set up multiple tracks reflection coefficient inversion equation, and in a particular embodiment, multitrace seismogram convolution model equation is such as Useful Y=AX+V represents.Wherein Y represents that multitrace seismogram, A represent seismic wavelet convolution matrix, and X represents that multiple tracks reflection is Number Sequence, V represents multiple tracks observation noise.

Initial value acquisition device 105, for obtaining the initial value of previously given multiple tracks sparse reflectivity model parameter. In the present invention, described multiple tracks sparse reflectivity model parameter includes non-negative parameter γ_i, covariance matrix B, noise variance λ and mean μ_x.Its each self-corresponding initial value can be expressed as γ⁰、B⁰、λ⁰、

Inverting device 106, combines multiple tracks reflection coefficient inversion equation to described by block Bayesian Learning Theory framework Multiple tracks sparse reflectivity model parameter is updated, the multiple tracks sparse reflectivity model parameter after being updated.Fig. 7 is anti- Drilling the structured flowchart of the embodiment one of device 106, as shown in Figure 7, in embodiment one, inverting device 106 includes:

Acquisition module 201, is used for obtaining maximum iteration time and model modification tolerable error, in specific embodiment In, maximum iteration time such as can use it_maxRepresenting, model modification tolerable error such as can represent with E.

Update coefficient determination module 202, initial for according to previously given multiple tracks sparse reflectivity model parameter Value combines multiple tracks reflection coefficient inversion equation and determines the renewal sparse reflection coefficient of multiple tracks by block Bayesian Learning Theory framework Model parameter.The mean μ of renewal can be determined by (10)_x, the multiple tracks model hyper parameter of renewal is determined by formula (19) γ_i, determined the covariance matrix B of renewal by formula (20), determined the noise variance λ of renewal by formula (21).? In specific embodiment, the multiple tracks sparse reflectivity model parameter after renewal can be expressed as γ¹、B¹、λ¹、

Update average acquisition module 203, for obtaining out renewal from the multiple tracks sparse reflectivity model parameter updated Average；

Initial value acquisition module 204, for from the initial value of previously given multiple tracks sparse reflectivity model parameter Obtain out the initial value of average；

Difference determines module 205, is used for the difference of the average determining described renewal and the initial value of average；

First judge module 206, for judging that described difference, whether less than model modification tolerable error, i.e. judges

Iteration module 207, for when the first described judge module is judged as NO, with the sparse reflection of the multiple tracks updated be The initial value of the multiple tracks sparse reflectivity model parameter that digital-to-analogue shape parameter iteration is previously given, after record iterations, returns Perform step S202, pass through Bayesian Learning Theory according to the initial value of previously given multiple tracks sparse reflectivity model parameter Framework determines the multiple tracks sparse reflectivity model parameter of renewal.In a particular embodiment, iteration can carry out table with j herein Show.

Namely in embodiment one, it may be judged whether when jumping out loop iteration condition, when being unsatisfactory forAnd Iterations now is 1, not up to maximum iteration time it_max, return the multiple tracks sparse reflection coefficient mould redefining renewal Shape parameter.

Follow-up when judging whether to jump out loop iteration condition, need to be unsatisfactory forAnd iterations is not up to Maximum iteration time it_maxTime, return the multiple tracks sparse reflectivity model parameter redefining renewal.

Fig. 8 is the structured flowchart of the embodiment two of inverting device, and as shown in Figure 8, in embodiment two, described is anti- Drill device 106 also to include:

First reflection coefficient output module 208, for when the first described judge module is judged as YES, output updates Multiple tracks sparse reflectivity model parameter, after the average in the multiple tracks sparse reflectivity model parameter of described renewal is inverting The sparse reflection coefficient of multiple tracks.

Namely in embodiment two, it may be judged whether when jumping out loop iteration condition, when meetingTime, jump Go out circulation, the multiple tracks sparse reflectivity model parameter that output updates.

Fig. 9 is the structured flowchart of the embodiment three of inverting device, and as shown in Figure 9, in embodiment three, described is anti- Drill device also to include:

Second judge module 209, for judging that whether described iterations is less than described maximum iteration time；

Second reflection coefficient output module 210, for when the second described judge module is judged as NO, output updates Multiple tracks sparse reflectivity model parameter, after in the multiple tracks sparse reflectivity model parameter of described renewal, average is inverting The sparse reflection coefficient of multiple tracks.

Namely in embodiment three, it may be judged whether when jumping out loop iteration condition, as j < it_maxBe false i.e. iteration time When number is not less than described maximum iteration time, jump out circulation, the multiple tracks sparse reflectivity model parameter that output updates.

A kind of multiple tracks sparse reflection coefficient multiple tracks sparse reflectivity model parameter that Figure 10 provides for the embodiment of the present invention The structured flowchart of embodiment two of Inversion System, as shown in Figure 10, in embodiment two, described system also includes:

Reflectivity model parameter estimation apparatus 107, for estimating by MAP estimation and EM algorithm Multiple tracks sparse reflectivity model parameter after described renewal.

As it has been described above, be the Inversion System of a kind of multiple tracks sparse reflectivity model parameter that the present invention provides, be given Multiple tracks based on convolution model sparse reflection coefficient inverting framework, derived based on convolution model time domain Bayes in detail Study reflection coefficient inversion formula, describes poststack single track reflection coefficient inversion method based on Method Using Relevance Vector Machine method, is given The detailed derivation of management loading reflection coefficient inverting based on multiple tracks temporal correlation, passes through Mahalanobis distance characterizes the dependency between multiple tracks reflection coefficient, utilizes this information to reach the reflection of multiple seismic channels and is Numerical digit puts the exact inversion with amplitude.Use the multiple tracks reflection coefficient inverting framework that the present invention provides, select the multiple tracks of the present invention Sparse reflection coefficient inversion equation, inverting has the reflection coefficient seismic channel of dependency.With traditional sparse constraint reflection coefficient Comparison of Inversions, the advantage that this method exists in terms of weak signal identification and noise immunity, it is finally reached subsurface reflective boundary The accurate inverting of corresponding parameter.

Below in conjunction with specific embodiment, technical scheme is discussed in detail.Figure 11 for the present invention provide concrete The inverting flow chart of multiple tracks sparse reflectivity model parameter in embodiment.As shown in Figure 11, the sparse reflection of the multiple tracks of the present invention Inverting reflection coefficient is converted into reflectivity model probability by Modulus Model parametric inversion scheme by Bayesian frame Hyper parameter estimation problem, and introduce the covariance matrix between multiple tracks reflection coefficient to retrain reflection coefficient refutation process, Finally give the reflection coefficient inversion schemes of overall situation sparse solution.By MAP estimation and EM algorithm to hyper parameter and Correlation matrix learns.It is utilized respectively formula (19), (20) and (21) in an iterative process and calculates γ_i, B and λ, pass through formula (10) mean μ is calculated_xAnd replace the reflection coefficient that need to estimate, until front and back the difference of the object function of twice iteration reaches tolerance by mistake Difference, the μ of final output_xBeing multiple tracks reflection coefficient impulse magnitude simultaneously, the call number at the nonzero value place of output γ is same The position of Shi Duodao reflection coefficient pulse, and output matrix B reaction is that the mutual of reflection coefficient impulse magnitude between multiple tracks is closed System.

Therefore, the flow process of multiple tracks sparse reflectivity model parametric inversion technology is as follows:

1) adjacent seismic channel number L for inverting reflection coefficient, maximum iteration time it are given_max, model modification is tolerated Error E；2) a given relatively accurate seismic wavelet；

3) generate seismic wavelet convolution matrix, set up multiple tracks reflection coefficient inversion equation；

4) by Bayesian frame inverting carries out for the reflection coefficient in given adjacent seismic channel:

A) non-negative parameter γ, the initial value of covariance matrix B are given for the noise level in earthquake record；

B) given mean μ_xInitial value；

C) γ, B and λ are updated respectively by formula (19), (20) and (21)；

D) mean μ is updated by formula (10)_x；

E) judging to jump out loop iteration condition, iteration j, if j is < it_maxAndThen return to 3c)；Otherwise jump out circulation, export γ, μ_x；

5) to X assignment, the index of nonzero element, reflection coefficient width in the position correspondence vector γ of multiple tracks non-zero reflection coefficient The corresponding μ of value_x。

Owing to adjacent multiple tracks reflection coefficient exists certain scope at the dependency of identical time location, input the most every time L needs to test.Although the present invention considers the change on time location, but for the stratigraphic model of big rise and fall, should Method needs to reduce number of channels L carrying out inverting reflection coefficient.

Figure 12 is the model schematic of true stratiform reflection coefficient in specific embodiment, as shown in Figure 12, meta time identical The reflection coefficient put has certain dependency.Reflection coefficient amplitude is gradually reduced, and increases over time reflection Coefficient positions moves closer to, 20ms reduce and put 2ms.At 120ms and at 121ms, two reflection coefficient are set, and relatively Amplitude is the least.

Figure 13 is the schematic diagram of the seismic wavelet selecting mixed phase wavelet in specific embodiment, and dominant frequency is 40Hz.

Figure 14 be in specific embodiment by the convolution operation of seismic wavelet and reflection coefficient generate without making an uproar earthquake record Schematic diagram, Direct Recognition cannot go out reflectance signature from earthquake record.The seismic channel that reflection coefficient entirety amplitude is relatively small None-identified goes out the reflection coefficient of thin layer.

Figure 15 is multiple tracks sparse reflection coefficient inversion result schematic diagram in specific embodiment, and Figure 16 is single track sparse Bayesian Reflection coefficient inversion result schematic diagram, Figure 17 is the root-mean-square error signal of inversion result and true model in specific embodiment Figure.From Figure 15-Figure 17, all reflection coefficient time locations in the reflection coefficient inversion method result that the present invention proposes With amplitude all can accurately inverting, and maintain horizontal dependency, reached to identify the high-resolution effect of thin layer.Right Ratio single track sparse Bayesian inversion method, in terms of weak signal identification, the present invention is with the obvious advantage.By to inversion result with true The root-mean-square error of model is it can be seen that multiple tracks sparse reflection coefficient inversion method accuracy of identification is high, and error is little, and inverting thin layer is anti- Penetrate coefficient ability higher.

Figure 18 is noisy composite traces schematic diagram in specific embodiment, overall section signal to noise ratio snr=10.Due to reflection letter Number amplitude is gradually reduced, and causes right-hand member seismic channel signal to noise ratio significantly lower than 10.

Figure 19 is multiple tracks sparse reflection coefficient inversion result schematic diagram, and Figure 20 is single track sparse Bayesian reflection coefficient inverting The root-mean-square error schematic diagram of result schematic diagram, Figure 21 inversion result and true model.From Figure 19-Figure 21, due to noise Interference, the result that the present invention obtains can not be completely superposed with true reflectivity model, but horizontal owing to considering reflection coefficient Dependency, inversion result is better than single track sparse Bayesian inversion result.By logarithm root-mean-square error figure, this conclusion also can be described.

Figure 22 is gas field, Sichuan province actual poststack seismic data schematic diagram, from section it can be seen that this areal geology Constructing relatively easy, reflection line-ups is clear, without bigger structural relief.Demarcated by well shake and obtain zero phase that matching degree is higher Position seismic wavelet.

Figure 23 is multiple tracks sparse reflection coefficient inversion result schematic diagram, and Figure 24 is the frequency spectrum pair of raw data and inversion result Compare schematic diagram.From Figure 23-Figure 24, results of real data shows, the inverting that the method using the present invention to propose obtains Result resolution significantly improves, and maintains the lateral continuity of original seismic profile.The reflection coefficient width of adjacent seismic channel Value change is rationally.By process before and after data frequency spectrum contrast it can be seen that in the range of seismic data effective band inverting tie Really amplitude spectrum keeps consistent with initial data amplitude spectrum.Reasonably widened frequency band, it is ensured that reflection coefficient inversion result simultaneously Accuracy.

In sum, the inversion method of the sparse reflection coefficient of a kind of multiple tracks that the present invention proposes and system, with poststack earthquake Data is input data, under block Bayesian Learning Theory framework, by false to multiple seismic channel reflection coefficient sparse prior And if the dependency between seismic channel reflection coefficient is it is assumed that utilize MAP estimation and EM algorithm to solve multiple tracks simultaneously Sparse Pulse reflection coefficient, what this technology obtained is reflection coefficient section, has the resolution of superelevation, and maintains section Lateral continuity, its computational efficiency is higher, for identify thin layer, the weak reflective information of renwing stratum, improve seismic data resolution Provide technical support.

The key problem in technology point of the present invention is: ask for multiple tracks model hyper parameter γ_i, covariance matrix B, noise variance λ and Mean μ_x。

Compared with existing technical scheme, the invention have the advantages that

The present invention considers stratum reflection coefficient seriality on locus, special by this dependency of numerical representation method Point, utilizes nonlinear inversion, improves the precision asking for reflection coefficient.Reflection coefficient position will be asked for and amplitude problem will be led to Crossing block Bayesian learning frame transformations is the estimation to model hyper parameter, simplifies refutation process, reduces corresponding calculating Amount.The parameter adaptive derivation algorithm of iteration makes inversion result not affected by initial given parameters value.In refutation process The positive definite matrix adding the horizontal correlation constraint of reflection coefficient innovatively is estimated, thus ensures the space of inversion result section Seriality.There is more preferable advantage in the method in terms of identifying weak reflective information and noise immunity, improves reflection coefficient inverting steady Qualitative.Better characteristics is shown in poststack multiple tracks fractal inversion method.

One of ordinary skill in the art will appreciate that all or part of flow process realizing in above-described embodiment method, Ke Yitong Crossing computer program and complete to instruct relevant hardware, described program can be stored in general computer read/write memory medium In, this program is upon execution, it may include such as the flow process of the embodiment of above-mentioned each method.Wherein, described storage medium can be magnetic Dish, CD, read-only store-memory body (Read-Only Memory, ROM) or random store-memory body (Random Access Memory, RAM) etc..

Those skilled in the art are it will also be appreciated that various functions that the embodiment of the present invention is listed are by hardware or soft Part realizes depending on specifically applying the design requirement with whole system.Those skilled in the art can be specific for every kind Application, it is possible to use the function described in the realization of various methods, but this realization is understood not to protect beyond the embodiment of the present invention The scope protected.

The present invention applies specific embodiment principle and the embodiment of the present invention are set forth, above example Explanation be only intended to help to understand method and the core concept thereof of the present invention；Simultaneously for one of ordinary skill in the art, According to the thought of the present invention, the most all will change, in sum, in this specification Hold and should not be construed as limitation of the present invention.

Claims (14)

1. an inversion method for the sparse reflection coefficient of multiple tracks, is characterized in that, described method includes:

Gather post-stack seismic data；

Seismic wavelet is extracted according to described post-stack seismic data；

Seismic wavelet convolution matrix is generated according to described seismic wavelet；

Multiple tracks reflection coefficient inversion equation is set up according to described seismic wavelet convolution matrix and described post-stack seismic data；

Obtain the initial value of previously given multiple tracks sparse Bayesian reflectivity model parameter；

Multiple tracks reflection coefficient inversion equation reflection coefficient sparse to described multiple tracks is combined by block Bayesian Learning Theory framework Model parameter is updated, the multiple tracks sparse reflectivity model parameter after being updated.

Method the most according to claim 1, is characterized in that, described multiple tracks sparse reflectivity model parameter includes non-negative Parameter, covariance matrix, noise variance and average.

Method the most according to claim 2, is characterized in that, described combines multiple tracks by block Bayesian Learning Theory framework Reflection coefficient inversion equation reflectivity model sparse to described multiple tracks parameter is updated, and the multiple tracks after being updated is sparse instead Penetrate Modulus Model parameter to include:

Obtain maximum iteration time and model modification tolerable error；

Initial value according to previously given multiple tracks sparse reflectivity model parameter combines multiple tracks reflection coefficient inversion equation and leads to Cross block Bayesian Learning Theory framework and determine the multiple tracks sparse reflectivity model parameter of renewal；

The average of renewal is obtained out from the multiple tracks sparse reflectivity model parameter updated；

The initial value of average is obtained out from the initial value of previously given multiple tracks sparse reflectivity model parameter；

Determine the difference of the average of described renewal and the initial value of average；

Whether the difference described in judgement is less than model modification tolerable error；

When being judged as NO, with the sparse reflection of the multiple tracks that the multiple tracks sparse reflectivity model parameter iteration updated is previously given it is The initial value of digital-to-analogue shape parameter, after record iterations, returns and performs according to the previously given sparse reflectivity model of multiple tracks The initial value of parameter determines the multiple tracks sparse reflectivity model parameter of renewal by Bayesian Learning Theory framework.

Method the most according to claim 3, is characterized in that, described combines multiple tracks by block Bayesian Learning Theory framework Reflection coefficient inversion equation reflectivity model sparse to described multiple tracks parameter is updated, and the multiple tracks after being updated is sparse instead Penetrate Modulus Model parameter also to include:

When described difference is less than or equal to model modification tolerable error, the multiple tracks sparse reflectivity model ginseng that output updates Number, the average in the multiple tracks sparse reflectivity model parameter of described renewal is the sparse reflection coefficient of multiple tracks after inverting.

Method the most according to claim 3, is characterized in that, described combines multiple tracks by block Bayesian Learning Theory framework Reflection coefficient inversion equation reflectivity model sparse to described multiple tracks parameter is updated, and the multiple tracks after being updated is sparse instead Penetrate Modulus Model parameter also to include:

Whether the iterations described in judgement is less than described maximum iteration time；

When being judged as NO, the multiple tracks sparse reflectivity model parameter that output updates, the sparse reflection of multiple tracks of described renewal is Average in digital-to-analogue shape parameter is the sparse reflection coefficient of the multiple tracks after inverting.

6. according to the method described in claim 3 or 4 or 5, it is characterized in that, according to previously given multiple tracks sparse reflection coefficient mould The initial value of shape parameter combines multiple tracks reflection coefficient inversion equation and determines that renewal is many by block Bayesian Learning Theory framework Road sparse reflectivity model parameter is carried out by equation below:

$\begin{array}{c}{x}^{*}\stackrel{\mathrm{\&Delta;}}{=}{\mathrm{\&mu;}}_x={({\mathrm{\&lambda;\&Sigma;}}_0^{-1}+D^{T}D)}^{-1}{D}^{T}y\\ ={\mathrm{\&Sigma;}}_0{D}^T{(\mathrm{\&lambda;}I+{\mathrm{D\&Sigma;}}_0{D}^T)}^{-1}y\end{array};$

$\begin{array}{cc}{\mathrm{\&gamma;}}_i\&LeftArrow;\frac{1}{L}{X}_i{B}^{-1}{X}_i^T+{\left({\mathrm{\&Xi;}}_x\right)}_{ii}& {\&ForAll;}_i\end{array};$

$B\&LeftArrow;\left(\frac{1}{M}\underset{i=1}{\overset{M}{\&Sigma;}}\frac{{\left({\mathrm{\&Xi;}}_x\right)}_{ii}}{{\mathrm{\&gamma;}}_i}\right)B+\frac{1}{M}\underset{i=1}{\overset{M}{\&Sigma;}}\frac{X_i^T{X}_i}{r_i};$

$D=A\&CircleTimes;I_L$

${\mathrm{\&Xi;}}_x={({\mathrm{\&Gamma;}}^{-1}+\frac{1}{\mathrm{\&lambda;}}{A}^{T}A)}^{-1}$

$\mathrm{\&Gamma;}\stackrel{\mathrm{\&Delta;}}{=}diag({\mathrm{\&gamma;}}_1,...,{\mathrm{\&gamma;}}_M)$

Wherein, x^*For reflection coefficient inversion result, μ_xFor average, λ is noise variance, and I is unit battle array, γ_iFor non-negative parameter, L is The number of channels of multiple tracks reflection coefficient, B is covariance matrix, and X is multiple tracks reflection coefficient sequence,Representing for any i value, N is Earthquake record time sampling is counted, and M is that the time sampling of reflectivity model is counted, and Y is multichannel seismic data, and A is earthquake The convolution matrix of ripple,Represent not Luo Beini this F norm computing of crow of vector.

Method the most according to claim 6, is characterized in that, described method also includes:

The multiple tracks sparse reflectivity model parameter after described renewal is estimated by MAP estimation and EM algorithm.

8. an Inversion System for the sparse reflection coefficient of multiple tracks, is characterized in that, described system includes:

Post-stack seismic data harvester, is used for gathering post-stack seismic data；

Seismic wavelet extraction device, for extracting seismic wavelet according to described post-stack seismic data；

Convolution matrix generation device, for generating seismic wavelet convolution matrix according to described seismic wavelet；

Inversion equation sets up device, for setting up according to described seismic wavelet convolution matrix and described post-stack seismic data Multiple tracks reflection coefficient inversion equation；

Initial value acquisition device, for obtaining the initial value of previously given multiple tracks sparse reflectivity model parameter；

Inverting device, for combining multiple tracks reflection coefficient inversion equation to described multiple tracks by block Bayesian Learning Theory framework Sparse reflectivity model parameter is updated, the multiple tracks sparse reflectivity model parameter after being updated.

System the most according to claim 8, is characterized in that, described multiple tracks sparse reflectivity model parameter includes non-negative Parameter, covariance matrix, noise variance and average.

System the most according to claim 9, is characterized in that, described inverting device includes:

Acquisition module, is used for obtaining maximum iteration time and model modification tolerable error；

Update coefficient determination module, combine many for the initial value according to previously given multiple tracks sparse reflectivity model parameter By block Bayesian Learning Theory framework, road reflection coefficient inversion equation determines that the sparse reflectivity model of multiple tracks of renewal is joined Number；

Update average acquisition module, for obtaining out the average of renewal from the multiple tracks sparse reflectivity model parameter updated；

Initial value acquisition module, for obtaining out all from the initial value of previously given multiple tracks sparse reflectivity model parameter The initial value of value；

Difference determines module, is used for the difference of the average determining described renewal and the initial value of average；

First judge module, for judging that whether described difference is less than described model modification tolerable error；

Iteration module, for when the first described judge module is judged as NO, with the sparse reflectivity model of multiple tracks updated The initial value of the multiple tracks sparse reflectivity model parameter that parameter iteration is previously given, after record iterations, returns and performs institute The renewal coefficient determination module stated.

11. systems according to claim 10, is characterized in that, described inverting device also includes:

First reflection coefficient output module, for when the first described judge module is judged as YES, the multiple tracks that output updates is dilute Dredging reflectivity model parameter, the average in the multiple tracks sparse reflectivity model parameter of described renewal is the multiple tracks after inverting Sparse reflection coefficient.

12. systems according to claim 10, is characterized in that, described inverting device also includes:

Second judge module, for judging that whether described iterations is less than described maximum iteration time；

Second reflection coefficient output module, for when the second described judge module is judged as NO, the multiple tracks that output updates is dilute Dredging reflectivity model parameter, the average in the multiple tracks sparse reflectivity model parameter of described renewal is the multiple tracks after inverting Sparse reflection coefficient.

13., according to the system described in claim 10 or 11 or 12, is characterized in that, described renewal coefficient determination module is by such as Lower formula is carried out:

$\begin{array}{c}{x}^{*}\stackrel{\mathrm{\&Delta;}}{=}{\mathrm{\&mu;}}_x={({\mathrm{\&lambda;\&Sigma;}}_0^{-1}+D^{T}D)}^{-1}{D}^{T}y\\ ={\mathrm{\&Sigma;}}_0{D}^T{(\mathrm{\&lambda;}I+{\mathrm{D\&Sigma;}}_0{D}^T)}^{-1}y\end{array};$

$\begin{array}{cc}{\mathrm{\&gamma;}}_i\&LeftArrow;\frac{1}{L}{X}_i{B}^{-1}{X}_i^T+{\left({\mathrm{\&Xi;}}_x\right)}_{ii}& {\&ForAll;}_i\end{array};$

$B\&LeftArrow;\left(\frac{1}{M}\underset{i=1}{\overset{M}{\&Sigma;}}\frac{{\left({\mathrm{\&Xi;}}_x\right)}_{ii}}{{\mathrm{\&gamma;}}_i}\right)B+\frac{1}{M}\underset{i=1}{\overset{M}{\&Sigma;}}\frac{X_i^T{X}_i}{r_i};$

$D=A\&CircleTimes;I_L$

${\mathrm{\&Xi;}}_x={({\mathrm{\&Gamma;}}^{-1}+\frac{1}{\mathrm{\&lambda;}}{A}^{T}A)}^{-1}$

$\mathrm{\&Gamma;}\stackrel{\mathrm{\&Delta;}}{=}diag({\mathrm{\&gamma;}}_1,...,{\mathrm{\&gamma;}}_M)$

Wherein, x^*For reflection coefficient inversion result, μ_xFor average, λ is noise variance, and I is unit battle array, γ_iFor non-negative parameter, L is The number of channels of multiple tracks reflection coefficient, B is covariance matrix, and X is multiple tracks reflection coefficient sequence,Representing for any i value, N is Earthquake record time sampling is counted, and M is that the time sampling of reflectivity model is counted, and Y is multichannel seismic data, and A is earthquake The convolution matrix of ripple,Represent vector not Luo Beini this norm computing of crow.

14. systems according to claim 13, is characterized in that, described system also includes:

Reflectivity model parameter estimation apparatus, for estimating described renewal by MAP estimation and EM algorithm After multiple tracks sparse reflectivity model parameter.