CN105425301A - Frequency domain three-dimensional irregular earthquake data reconstruction method - Google Patents

Abstract

The invention brings forward a frequency domain three-dimensional irregular earthquake data reconstruction method. The method is characterized in that first of all, three-dimensional earthquake data in a time domain is converted to a frequency domain by use of Fourier transform, and then, a projection onto convex set (POCS) algorithm is employed and curvelet transform capable of describing localized features of the earthquake data is introduced; and in an iteration process, a new threshold parameter attenuating according to an index rule is brought forward, and each frequency slice is individually reconstructed by use of a soft threshold operator, such that the iteration frequency is reduced, the reconstruction precision is improved, and the purpose of reconstructing the three-dimensional earthquake data is realized. According to the invention, the new threshold parameter attenuating according to the index rule is brought forward and each frequency slice is reconstructed in individually by use of the soft threshold operator, such that the disadvantage of quite slow convergence speed of a conventional threshold parameter is overcome, the calculation complexity of an algorithm is reduced, the calculation efficiency is substantially improved, and the operation time is reduced.

Description

The three-dimensional irregular earthquake data re-establishing method of a kind of frequency field

Technical field

What the present invention relates to is a kind of method lacking Reconstruction of seismic data, specifically the three-dimensional irregular earthquake data re-establishing method of a kind of frequency field.

Technical background

In the wild in data acquisition, due to collecting device, the restriction of project cost and topographic condition, geological data carries out irregular sampling along direction in space being usually less than nyquist sampling theorem requirement, thus cause the geological data that collects irregular, imperfect, there is spatial aliasing, be difficult to the requirement meeting subsequent treatment, such as velocity analysis, Free Surface Multiple attenuation, migration etc., in this case, need to develop good Reconstruction of seismic data method, reconstruct the seismic trace of disappearance, be met the geological data of processing requirements, and the method also can instruct field data collection, compression earthquake data acquisition amount.

Data re-establishing method based on warp wavelet does not need the prior imformation of underground structure, and can reconstruct the seismic trace of rule disappearance and irregular disappearance, and computing velocity is fast, precision is high, is widely used.Xue Nian etc. adopt threshold value process of iteration, achieve two-dimension earthquake data reconstruction (the geological data interpolation converted based on Curvelet and the denoising [Master's thesis] based on warp wavelet, Southwest Jiaotong University, 2010), Liu Guochang etc. adopt convex set projection algorithm, bent wave zone achieve based on POCS algorithm two-dimension earthquake data reconstruction (based on curvelet conversion disappearance geological data interpolation method [J]. geophysical prospecting for oil, 2011,46 (2), 237-245.), the method more can rebuild non-linear lineups, more be applicable to the seismic wave field with anisotropic feature, but do not have to further investigate the three dimension data reconstruct method based on warp wavelet, and select also not study further to threshold parameter, only adopt index threshold parameter equation, its rebuild after effect or limited.And posttension China wait achieve on this basis based on the 3D seismic data of warp wavelet rebuilds (based on jitter sampling and warp wavelet 3D seismic data reconstruction [J]. Chinese Journal of Geophysics, 2013,56 (5): 1637-1649.), by successively rebuilding time slice, obtaining and rebuilding effect preferably, and inquired into threshold parameter selection strategy.Wang Benfeng etc. propose bent ripple method for reconstructing under different threshold value (POCS combines bent wave zone Reconstruction of seismic data method [J] of Jitter sampling theory of improvement. Chinese Journal of Geophysics, 2015,50 (1): 20-27.), but the method is rebuild in time domain, and only adopt hard-threshold parameter to process, when data reconstruction comprises high and low-yield amplitude inphase axle simultaneously, because POCS algorithm is in front iterative process several times, lower than the lineups coefficient below threshold parameter by filtering, therefore these methods are rebuilding the limited efficiency of short arc significant wave lineups, the low-yield lineups of seismic signal can not be recovered completely and improve the signal to noise ratio (S/N ratio) in local anomaly district, and iterations is many, reconstruction efficiency is low, constrains the further application of the method.

Summary of the invention

The object of the invention is in order to reconstructing of quick high accuracy road 3D seismic data can be lacked, and provide the three-dimensional irregular earthquake data re-establishing method of a kind of frequency field.

Technical solution of the present invention:

The present invention proposes the three-dimensional irregular earthquake data re-establishing method of frequency field, first utilize Fourier transform that time domain 3D seismic data is transformed to frequency field, then adopt convex set projection algorithm idea and introduce the warp wavelet can portraying geological data local characteristic.Propose the new threshold parameter of exponentially decaying in an iterative process, adopt soft-threshold operator, each frequency slice is rebuild separately, reduce iterations and improve reconstruction precision, thus reaching the object of the irregular geological data of reconstruction of three-dimensional.

The three-dimensional irregular earthquake data re-establishing method of a kind of frequency field, the Problems of Reconstruction of geological data is described as recovering complete data by one group of deficiency of data by the effect of linear operator, supposes as lower linear forward model:

y _obs＝Mf

Here y _obs∈ R ⁿthe geological data that representative gathers; F ∈ R ⁿ, and N>=n, represents to be reconstructed without alias data; M ∈ R ^{n × N}represent diagonal matrix, its element 1 and 0 represents known seismic trace and unknown seismic trace respectively, and suppose that coefficient x is the rarefaction representation of f in bent wave zone C, then above-mentioned equation is:

Y _obs=Ax and

Here subscript ^hrepresent associate matrix, as the data y from collection _obswhen middle reconstruction is without alias data f, because x is sparse, thus Corresponding Sparse Algorithm can be adopted to solve this underdetermined equation.

After sparse promotion inverting, reconstruction signal by determine, simultaneously

$\begin{array}{ccc}\stackrel{\‾}{x}=\arg \min \left|\right|x|{|}_1& subjectto& y_{obs}=Ax\end{array}$

In this expression formula, represent estimated value, L ₁norm is defined as x [i] is i-th element in vector x, and by solving above-mentioned equation, the original data without alias are just rebuild out;

In Reconstruction of seismic data process, adopt convex set projection algorithm, concrete reconstruction procedures is as follows:

Step 1: input has the 3D seismic data in disappearance road in the time domain, then adopts Fourier transform, data reconstruction is transformed to frequency field from time domain;

Step 2: adopt warp wavelet to carry out sparse transformation to frequency field 3D seismic data, obtain bent wave system number in bent wave zone, select suitable threshold parameter formula λ according to the size of bent wave system number _i(i=1,2,3, N, wherein N is iterations) and formula;

Step 3: in bent wave zone, adopts soft-threshold operator to process, is also namely greater than threshold value λ _ibent wave system number deduct the value of a threshold size, be less than threshold value-λ _ibent wave system number add the value of a threshold size, and other bent wave system number zero setting;

Step 4: the bent wave system number after thresholding is done contrary flexure wave conversion and obtains time domain geological data, and then with disappearance geological data do not lack seismic trace be filled into inverse transformation after geological data in go;

Step 5: finally the geological data that obtain is substituted into step (2), re-start iteration, until run N end, then does inverse-Fourier transform to the geological data after iteration N time and namely obtains final reconstructed results.

Further, described warp wavelet is defined as:

$C(j,l,k)=<f,{\mathrm{\&phi;}}_{j,l,k}>={\&Integral;}_{R^2}f\left(x\right)\stackrel{\&OverBar;}{{\mathrm{\&phi;}}_{j,l,k}\left(x\right)}dx$

In formula: φ _{j, l, k}represent bent wave function, j, l, k represent yardstick respectively, direction and location parameter, and f (x) is geological data, and its frequency field definition is:

$C(j,l,k)=\frac{1}{{\left(2\mathrm{\&pi;}\right)}^2}\&Integral;\hat{f}\left(\mathrm{\&omega;}\right)\stackrel{\&OverBar;}{{\widehat{\mathrm{\&phi;}}}_{j,l,k}\left(\mathrm{\&omega;}\right)}d\mathrm{\&omega;}=\frac{1}{{\left(2\mathrm{\&pi;}\right)}^2}\&Integral;\hat{f}\left(\mathrm{\&omega;}\right)U_j\left(R_{{\mathrm{\&theta;}}_l}\mathrm{\&omega;}\right)e^{i<x_k^{(j,l)},\mathrm{\&omega;}>}d\mathrm{\&omega;}$

The bent wave system number obtained after conversion, available C{j}{l} (k ₁, k ₂) represent its structure, wherein j represents yardstick, and l represents direction, (k ₁, k ₂) represent matrix coefficient on j yardstick l direction.

Further, described soft-threshold operator, its expression formula is as follows:

$T_k(i,j)=\left\{\begin{array}{cc}{S}_{k-1}(i,j)-{\mathrm{\&lambda;}}_k,& S_{k-1}(i,j)\&GreaterEqual;{\mathrm{\&lambda;}}_k\\ S_{k-1}(i,j)+{\mathrm{\&lambda;}}_k,& S_{k-1}(i,j)\&le;-{\mathrm{\&lambda;}}_k\\ 0,& \left|S_{k-1}(i,j)\right|<{\mathrm{\&lambda;}}_k\end{array}\right.,{\mathrm{\&lambda;}}_k\&Element;\mathrm{\&lambda;}$

T _krepresent soft-threshold operator, S _k-1represent the bent wave system number of the data reconstruction that kth-1 iteration obtains, meet f _tand F ^-1 _trepresent the positive inverse-Fourier transform about time variable t, represent the time domain data reconstruction that kth-1 iteration obtains, wherein λ represents that N ties up threshold value set, λ={ λ ₁, λ ₂, λ _n, and meet λ ₁> λ ₂> > λ _n, N represents maximum iteration time.

Further, threshold parameter formula, its expression formula is as follows:

${\mathrm{\&lambda;}}_k=Max\&CenterDot;e^{\left(\ln \right(\mathrm{\&epsiv;})-\ln (Max)\&rsqb;{\left(\frac{k-1}{N-1}\right)}^{1/3}}$

Wherein Max is | the maximal value of CFy|, the i.e. maximal value of warp wavelet absolute coefficient.N is total iterations, and ε is the little value close to zero, relevant with the energy of noise in data.

Advantage of the present invention: time domain 3D seismic data is transformed to frequency field by utilizing Fourier transform by the present invention, and adopt convex set projection algorithm to carry out data reconstruction with the warp wavelet can portraying geological data local characteristic, thus avoid low-yield significant wave lineups coefficient by filtering, improve disappearance Reconstruction of seismic data precision, the new threshold parameter of exponentially decaying is proposed in an iterative process simultaneously, adopt soft-threshold operator, each frequency slice is rebuild separately, iterations can be reduced, reduce the computation complexity of algorithm, improve counting yield significantly, save operation time.

Accompanying drawing explanation

Fig. 1 represents the three-dimensional irregular earthquake data re-establishing method process flow diagram in embodiment of the present invention medium frequency territory.

The different threshold parameter of Fig. 2 rebuilds signal to noise ratio (S/N ratio) curve comparison figure.

Fig. 3 is two-dimensional random lack sampling and reconstructed results comparison diagram thereof.

Fig. 4 is for rebuilding front and back 2-d spectrum comparison diagram.

Fig. 5 is signal to noise ratio (S/N ratio) and sampling rate relation curve comparison diagram.

Fig. 6 is noisy data reconstruction processes figure.

Embodiment

Following case study on implementation for illustration of the present invention, but is not used for limiting the scope of the invention.

Embodiment 1

The step realizing the method mainly comprises, the structure of Reconstructed equation, and time domain transforms to frequency field, convex set projection reconstruction algorithm, threshold parameter process etc.Concrete steps are as follows:

Step 1: the structure of Reconstructed equation.The Problems of Reconstruction of geological data can be described as recovering complete data by one group of deficiency of data by the effect of linear operator, supposes as lower linear forward model

y _obs＝Mf

Here y _obs∈ R ⁿthe geological data that representative gathers; F ∈ R ⁿ, and N>=n, represents to be reconstructed without alias data; M ∈ R ^{n × N}represent diagonal matrix, its element 1 and 0 represents known seismic trace and unknown seismic trace respectively, and suppose that coefficient x is the rarefaction representation of f in bent wave zone C, then above-mentioned equation can be write as

Y _obs=Ax and

Here subscript ^hrepresent associate matrix, when from the data y collected _obswhen middle reconstruction is without alias data f, because x is sparse, thus this underdetermined equation is made to have solution.

After sparse promotion inverting, reconstruction signal by determine, simultaneously

$\begin{array}{ccc}\stackrel{\&OverBar;}{x}=\arg \min \left|\right|x|{|}_1& subjectto& y_{obs}=Ax\end{array}$

In this expression formula, represent estimated value, L ₁norm is defined as x [i] is i-th element in vector x, and by solving above-mentioned underdetermined equation, the original complete and geological data of rule just can be rebuild out.

Step 2: time domain transforms to frequency field.The present invention adopts Fourier transform to calculate time domain geological data, and its discrete Fourier transform (DFT) formula is as follows:

$Y\left(m\mathrm{\&Delta;}f\right)=\mathrm{\&Delta;}t\underset{n=0}{\overset{N-1}{\&Sigma;}}y\left(n\mathrm{\&Delta;}t\right)e^{-i2\mathrm{\&pi;}m\mathrm{\&Delta;}f*n\mathrm{\&Delta;}t}$

In formula: m, n are integer.m＝0,1,…M-1；n＝0,1,…N-1。N △ t=T is the window length of spectrum analysis.

Step 3: convex set projection reconstruction algorithm.The shortcoming low for traditional reconstruction algorithm reconstruction precision and computing velocity is slow, propose the algorithm carrying out rebuilding in frequency field, and adopt soft-threshold operator and the new threshold parameter formula of exponentially decaying, its algorithm iteration expression formula is:

${\stackrel{\&OverBar;}{f}}_k={\stackrel{\&OverBar;}{y}}_{obs}+{\mathrm{MF}}_t^{-1}{C}^{-1}{T}_k{\mathrm{CF}}_t{\stackrel{\&OverBar;}{f}}_{k-1},k=1,2,\mathrm{...},N$

Wherein, f _krepresent the time domain data reconstruction that kth time iteration obtains, represent that acquired original is to data y _obs(t, x, y), meets f _tand F ^-1 _trepresent the positive inverse-Fourier transform about time variable t, C and C ^-1represent positive and negative warp wavelet, its warp wavelet is defined as:

$C(j,l,k)=<f,{\mathrm{\&phi;}}_{j,l,k}>={\&Integral;}_{R^2}f\left(x\right)\stackrel{\&OverBar;}{{\mathrm{\&phi;}}_{j,l,k}\left(x\right)}dx$

In formula: φ _{j, l, k}represent bent wave function, j, l, k represent yardstick respectively, direction and location parameter, and f (x) is geological data, and its frequency field definition is:

$C(j,l,k)=\frac{1}{{\left(2\mathrm{\&pi;}\right)}^2}\hat{f}\left(\mathrm{\&omega;}\right)\stackrel{\&OverBar;}{{\widehat{\mathrm{\&phi;}}}_{j,l,k}\left(\mathrm{\&omega;}\right)}d\mathrm{\&omega;}\frac{1}{{\left(2\mathrm{\&pi;}\right)}^2}\&Integral;\hat{f}\left(\mathrm{\&omega;}\right)U_j\left(R_{{\mathrm{\&theta;}}_l}\mathrm{\&omega;}\right)e^{i<x_k^{(j,l)},\mathrm{\&omega;}>}d\mathrm{\&omega;}$

The bent wave system number obtained after conversion, available C{j}{l} (k ₁, k ₂) represent its structure, wherein j represents yardstick, and l represents direction, (k ₁, k ₂) represent matrix coefficient on j yardstick l direction.

T _krepresent soft-threshold operator, its element meets:

$T_k(i,j)=\{\begin{array}{cc}{S}_{k-1}(i,j)-{\mathrm{\&lambda;}}_k,& S_{k-1}(i,j)\&GreaterEqual;{\mathrm{\&lambda;}}_k\\ S_{k-1}(i,j)+{\mathrm{\&lambda;}}_k,& S_{k-1}(i,j)\&le;-{\mathrm{\&lambda;}}_k\\ 0,& \left|S_{k-1}(i,j)\right|<{\mathrm{\&lambda;}}_k\end{array},{\mathrm{\&lambda;}}_k\&Element;\mathrm{\&lambda;}$

S _k-1represent the bent wave system number of the data reconstruction that kth-1 iteration obtains, meet λ represents that N ties up threshold value set, λ={ λ ₁, λ ₂, λ _n, and meet λ ₁> λ ₂> > λ _n, N represents maximum iteration time.

Step 4: threshold parameter formula.In an iterative process, its thresholding parameter value generally changes from big to small, also namely most of bent wave system number is removed in front iteration several times, only leave the lineups that energy is stronger, and in the end several times in iteration, the bent wave system number of the overwhelming majority retains, and only removes the random noise that some energy are less, i.e. threshold parameter λ _imeet:

||CFy|| _∞＝λ ₁＞λ ₂＞···＞ε

In formula, ε is the little value close to zero, relevant with the energy of noise in data, different pieces of information ε value difference to some extent.The present invention on the basis of linear threshold parameter equation and index threshold parameter equation in the past, in conjunction with frequency field warp wavelet feature, propose to select by the threshold parameter of rule decay, under the prerequisite ensureing precision, can improve speed of convergence sooner, this threshold parameter formula is:

${\mathrm{\&lambda;}}_k=Max\&CenterDot;e^{\left(\ln \right(\mathrm{\&epsiv;})-\ln (Max)\&rsqb;{\left(\frac{k-1}{N-1}\right)}^{1/3}}$

Wherein Max is | the maximal value of CFy|, the i.e. maximal value of warp wavelet absolute coefficient.N is total iterations, and ε is the minimal value close to zero, relevant with the energy of noise in data.

Realizing the method concrete operations is:

(1) threshold parameter comparative analysis

In order to correlation data rebuilds effect better, definition signal to noise ratio snr=20log ₁₀|| f ₀|| ₂/ || f-f ₀|| ₂, unit is dB, wherein f ₀represent primary model data, f represents reconstructed results, and signal to noise ratio (S/N ratio) is higher, represent reconstructed results and model data is more close, effect is unreasonable to be thought, the present invention is simultaneously in 3D seismic data process of reconstruction, the scale parameter of warp wavelet is 5, and the angle number on the thickest yardstick is 8.

Then the threshold parameter formula adopting conventional linear threshold value, index threshold and the present invention to propose carries out data reconstruction to the data (Fig. 3 (c) represents the time slice two-dimensional random lack sampling figure of 3D seismic data) after theoretical model 50% lack sampling.Reconstructed results is as shown in Figure 2 a (when Fig. 2 (a) represents iterations 50, each iterations and signal to noise ratio (S/N ratio) figure), Fig. 2 (a) for maximum iteration time be 50 times time, each iterations threshold parameter different from three kinds rebuilds rear Between Signal To Noise Ratio curve map, after can finding out the threshold parameter reconstruction that in each iterative process, the present invention proposes, signal to noise ratio (S/N ratio) is the highest, next is index threshold parameter, is finally only linear threshold parameter.Maximum iteration time is changed from 5 ~ 100 times simultaneously, Between Signal To Noise Ratio curve map after calculating each maximum iteration time and rebuilding, as Suo Shi Fig. 2 (b) (Fig. 2 (b) represents maximum iteration time and its signal to noise ratio (S/N ratio) figure of each computing), therefrom also can find out in different maximum iteration time, it is maximum that the threshold parameter that the present invention proposes rebuilds rear signal to noise ratio (S/N ratio), and can find out, wanting to obtain signal to noise ratio (S/N ratio) is the reconstructed results of 20dB, linear threshold needs iteration about 70 times, index threshold needs iteration about 17 times, and the threshold parameter that the present invention proposes only needs iteration about 13 times, certain computing time can be saved.And for the index square root threshold parameter that index threshold and the present invention propose, after iterations is greater than 30 times, signal to noise ratio (S/N ratio) recruitment is relatively less with iterations increase, and therefore, the three-dimensional irregular earthquake data re-establishing method of frequency field of the present invention all selects 30 iteration.

(2) data reconstruction processes

Adopt sound wave finite difference method, simulate the Seismic forward record under uniform grid sampling, and just drill geological data by wave detector to obtained, shot point and time carry out being arranged in 3-D data volume.Ideal data are (Fig. 3 (c) represents the time slice two-dimensional random lack sampling figure of 3D seismic data) as Suo Shi Fig. 3 (c), wherein time slice is 0.48s, shot record migration respective distances is 1512m (128 big gun), and the distance that common receiver is corresponding is 1512m (128 road).Then the three-dimensional irregular earthquake data re-establishing method of frequency field of the present invention is adopted to rebuild it, rebuild as shown in Figure 3 d (Fig. 3 (d) represents the two-dimensional random lack sampling reconstructed results of cutting into slices sometime), signal to noise ratio (S/N ratio) after reconstruction is 31.67dB, and reconstruction precision is higher.From the analysis of reconstruction surrounding time section 2-d spectrum, (the 2-d spectrum analysis chart of a certain section of Fig. 4 (e) original earthquake data as shown in Figure 4; Fig. 4 (f) adopts the 2-d spectrum analysis chart of the inventive method reconstructed results), 2-d spectrum after the inventive method reconstruction is closest to the 2-d spectrum of raw data, this also illustrates the inventive method can reflect seismic event before variation characteristic, remain weak significant wave lineups, improve the signal to noise ratio (S/N ratio) in local anomaly district.

(3) different sampling rate compares

In order to contrast the reconstruction effect after a peacekeeping two dimension lack sampling, random lack sampling rate is kept to increase progressively from 10 ~ 80%, and then rebuild, signal to noise ratio (S/N ratio) after record reconstruction and the relation between sampling rate, as Suo Shi Fig. 5 (g) (Fig. 5 (g) represents the sampling rate after adopting warp wavelet to rebuild and Between Signal To Noise Ratio figure), therefrom can find out that sampling rate is fewer, effect after reconstruction is poor, but comparatively speaking, the effect that two dimension lack sampling is rebuild rebuilds better effects if than one dimension lack sampling, therefore need to adopt three-dimensional method for reconstructing, increase the information in other directions, to make up the deficiency of low-dimensional method for reconstructing.Simultaneously in order to the superiority of outstanding warp wavelet, itself and the three dimension data reconstruct algorithm based on Fourier transform are compared, Fourier transform also adopts 30 iteration when rebuilding, as Suo Shi Fig. 5 (h) (Fig. 5 (h) represents the sampling rate after adopting Fourier transform to rebuild and Between Signal To Noise Ratio figure), can find out, when adopting identical sampling rate, the effect that warp wavelet is rebuild is far superior to Fourier basis, this also illustrates warp wavelet more can reflect seismic event before variation characteristic, improve the reconstruction precision in local anomaly district.Also can find out, no matter be warp wavelet or Fourier transform, the reconstruction effect of its 3D seismic data is all better than two-dimension earthquake data reconstruction effect simultaneously.

(4) noisy data reconstruction

Because geological data all contains noise usually, so the noise resisting ability of the inventive method under needing to check two-dimensional random lack sampling, gaussian random noise is added in Fig. 3 (c) time slice, as Suo Shi Fig. 6 (i) (a certain noisy slice map that Fig. 6 (i) is 3D seismic data), then carry out anti-noise and rebuild simulated experiment, two-dimensional random lack sampling result is (Fig. 6 (j) is two-dimensional random lack sampling slice map sometime) as Suo Shi Fig. 6 (j), then the inventive method is adopted to rebuild it, its reconstructed results is (Fig. 6 (k) adopts the inventive method reconstructed results figure) as Suo Shi Fig. 6 (k), signal to noise ratio (S/N ratio) 10.92dB, it rebuilds front and back error section (Fig. 6 (l) rebuilds front and back Error Graph) as Suo Shi Fig. 6 (l), can be found out by contrast, effective information after reconstruction is almost constant, data reconstruction effect is better, this also illustrates that the inventive method has good noise resisting ability in Reconstruction of seismic data, can apply in the process of real data.

Claims (6)

1. the three-dimensional irregular earthquake data re-establishing method of frequency field, it is characterized in that, first utilize Fourier transform that time domain 3D seismic data is transformed to frequency field, then adopt convex set projection algorithm idea and introduce the warp wavelet can portraying geological data local characteristic; Propose the new threshold parameter of exponentially decaying in an iterative process, adopt soft-threshold operator, each frequency slice is rebuild separately, thus reduce iterations and improve reconstruction precision, reach the object of reconstruction of three-dimensional geological data.

2. the three-dimensional irregular earthquake data re-establishing method of a kind of frequency field according to claim 1, it is characterized in that, the Problems of Reconstruction of geological data is described as recovering complete data by one group of deficiency of data by the effect of linear operator, supposes as lower linear forward model:

y _obs＝Mf

Here y _obs∈ R ⁿthe geological data that representative gathers; F ∈ R ⁿ, and N>=n, represents to be reconstructed without alias data; M ∈ R ^{n × N}represent diagonal matrix, its element 1 and 0 represents known seismic trace and unknown seismic trace respectively, and suppose that coefficient x is the rarefaction representation of f in bent wave zone C, then above-mentioned equation is:

Y _obs=Ax and

Here subscript ^hrepresent associate matrix, as the data y from collection _obswhen middle reconstruction is without alias data f, because x is sparse, thus Corresponding Sparse Algorithm can be adopted to solve this underdetermined equation;

After sparse promotion inverting, reconstruction signal by determine, simultaneously

$\begin{array}{ccc}\stackrel{\&OverBar;}{x}=\arg \min \left|\right|x|{|}_1& subjectto& y_{obs}\end{array}=Ax$

In this expression formula, represent estimated value, L ₁norm is defined as x [i] is i-th element in vector x, and by solving above-mentioned equation, the original data without alias are just rebuild out.

3. the three-dimensional irregular earthquake data re-establishing method of a kind of frequency field according to claim 1 and 2, is characterized in that, in Reconstruction of seismic data process, adopt convex set projection algorithm, concrete reconstruction procedures is as follows:

Step 1: input has the 3D seismic data in disappearance road in the time domain, then adopts Fourier transform, data reconstruction is transformed to frequency field from time domain;

Step 2: adopt warp wavelet to carry out sparse transformation to frequency field 3D seismic data, obtain bent wave system number in bent wave zone, select suitable threshold parameter formula λ according to the size of bent wave system number _i(i=1,2,3 ..., N, wherein N is iterations);

Step 3: in bent wave zone, adopts soft-threshold operator to process, is also namely greater than threshold value λ _ibent wave system number deduct the value of a threshold size, be less than threshold value-λ _ibent wave system number add the value of a threshold size, and other bent wave system number zero setting;

Step 4: the bent wave system number after thresholding is done contrary flexure wave conversion and obtains time domain geological data, and then with disappearance geological data do not lack seismic trace be filled into inverse transformation after geological data in go;

Step 5: finally the geological data that obtain is substituted into step (2), re-start iteration, until run N end, then does inverse-Fourier transform to the geological data after iteration N time and namely obtains final reconstructed results.

4. the three-dimensional irregular earthquake data re-establishing method of a kind of frequency field according to claim 3, it is characterized in that, described warp wavelet is defined as:

In formula: φ _{j, l, k}represent bent wave function, j, l, k represent yardstick respectively, direction and location parameter, and f (x) is geological data, and its frequency field definition is:

The bent wave system number obtained after conversion, available C{j}{l} (k ₁, k ₂) represent its structure, wherein j represents yardstick, and l represents direction, (k ₁, k ₂) represent matrix coefficient on j yardstick l direction.

5. the three-dimensional irregular earthquake data re-establishing method of a kind of frequency field according to claim 3, it is characterized in that, described soft-threshold operator, its expression formula is as follows:

$T_k(i,j)=\left\{\begin{array}{cc}{S}_{k-1}(i,j)-{\mathrm{\&lambda;}}_k,& S_{k-1}(i,j)\&GreaterEqual;{\mathrm{\&lambda;}}_k\\ S_{k-1}(i,j)+{\mathrm{\&lambda;}}_k,& S_{k-1}(i,j)\&le;-{\mathrm{\&lambda;}}_k\\ 0,& \left|S_{k-1}(i,j)\right|<{\mathrm{\&lambda;}}_k\end{array}\right.,{\mathrm{\&lambda;}}_k\&Element;\mathrm{\&lambda;}$

T _krepresent soft-threshold operator, S _k-1represent the bent wave system number of the data reconstruction that kth-1 iteration obtains, meet f _tand F ^-1 _trepresent the positive inverse-Fourier transform about time variable t, represent the time domain data reconstruction that kth-1 iteration obtains, wherein λ represents that N ties up threshold value set, λ={ λ ₁, λ ₂..., λ _n, and meet λ ₁> λ ₂> ... > λ _n, N represents maximum iteration time.

6. the three-dimensional irregular earthquake data re-establishing method of a kind of frequency field according to claim 3, it is characterized in that, threshold parameter formula, its expression formula is as follows:

${\mathrm{\&lambda;}}_k=Max\&CenterDot;e^{(ln\left(\mathrm{\&epsiv;}\right)-ln\left(Max\right)\&rsqb;{\left(\frac{k-1}{N-1}\right)}^{1/3}}$

Wherein Max is the maximal value of warp wavelet absolute coefficient; N is total iterations, and ε is the little value close to zero, relevant with the energy of noise in data.