CN105549078A - Five-dimensional interpolation processing method and apparatus of irregular seismic data - Google Patents

Abstract

The invention discloses a five-dimensional interpolation processing method and apparatus of irregular seismic data. The method comprises the following steps: transforming the irregular seismic data to a frequency domain and a space domain; for frequency slices of each piece of seismic data transformed to the frequency domain and the space domain, according to a fast Fourier transform (FFT) algorithm, performing product operation on matrixes and vectors of the seismic data so as to solve frequency wave number domain data after five-dimensional interpolation; and transforming the frequency wave number domain data after the five-dimensional interpolation to the frequency domain and the space domain and then transforming to a time domain and the space domain. The technical scheme provided by the invention realizes rapid processing of the irregular seismic data, the vast majority of matrix and vector product operation in the whole processing process is realized by use of the FFT algorithm, and thus the seismic data processing efficiency is improved.

Description

Five dimension interpolation process method and devices of irregular geological data

Technical field

The present invention relates to seismic exploration technique field, particularly a kind of five dimension interpolation process method and devices of irregular geological data.

Background technology

Seismic prospecting is mainly divided into seismic data acquisition, process and explains three phases.In the seismic data acquisition stage, no matter be land exploration or seafari, by the impact of various complicated factor, the spatial sampling of actual acquisition data can only be approximate rule, or even irregular.Such as, for land exploration, in city, the area of highway, the surface conditions complexity such as the network of waterways, cannot the recording geometry of layout rules, irregular collection must be carried out; For seafari, by the impact of ocean current, there is the phenomenon of feathering in actual cable survey line, is difficult to the spatial sampling ensureing rule equally.In the seism processing stage, many data processing algorithms (as migration algorithm and multiple removal algorithm) all need or at least have benefited from the geological data of rule space sampling.In addition, carry out in flakes or time-lapse seismic data processing time, the recording geometry parameter of the geological data of different batches collection is generally different, and this also relates to the spatial sampling regularization problem of geological data.Therefore, it is very necessary for utilizing special interpolation algorithm to solve geological data spatial sampling irregular problem.

At present, most of geological data interpolation algorithm is also only confined to three-dimensional interpolation, but the geological data of field acquisition is the function of five dimension coordinates in essence, bidimensional is used for determining shot point locus, bidimensional is used for determining geophone station locus, also have one dimension to be used for determining the sampled point time, development five dimension interpolation algorithm can utilize the geological data of collection more fully, and then obtains better interpolation result.Compared with D interpolation algorithm, the primary difficult problem that five dimension interpolation algorithms face is exactly that googol is according to amount and calculated amount.The three-dimensional minimum weight norm interpolation algorithm of Liu and Sacchi (2004) is extended to five dimensions by Trad (2009), during algorithm realization, all matrixes and vector product computing can utilize Fast Fourier Transform (FFT) (FFT) to complete, thus ensure that counting yield, but algorithm requires that the spatial sampling of geological data is regular, irregular acquiring seismic data cannot be used for, there is significant limitation.Jin (2010) supposes based on irregular spatial sampling, propose the five dimension geological data interpolation algorithms based on decay minimum norm Fourier inverting, algorithm adopts unequal interval Fast Fourier Transform (FFT) (NFFT) to realize matrix and vector product computing frequently, improve counting yield to a certain extent, but because the counting yield of high dimensional data NFFT is far away from FFT, and NFFT itself is a kind of approximate data, the method for Jin is still to be improved in counting yield.

Summary of the invention

Embodiments provide a kind of five dimension interpolation process methods of irregular geological data, in order to improve the efficiency of seismic data process, the method comprises:

By irregular earthquake data transformation to frequency field and spatial domain;

Each is transformed to the frequency slice of the geological data of frequency field and spatial domain, according to fast fourier transform algorithm, product calculation is carried out to the matrix of geological data and vector, asks for the frequency-wavenumber domain data after five dimension interpolation;

By the frequency-wavenumber domain data transformation after five dimension interpolation to frequency field and spatial domain, then transform to time domain and spatial domain.

Embodiments provide a kind of five dimension interpolation processors of irregular geological data, in order to improve the efficiency of seismic data process, this device comprises:

Geological data conversion module, for by irregular earthquake data transformation to frequency field and spatial domain;

Five dimension interpolation processing modules, for transforming to the frequency slice of the geological data of frequency field and spatial domain for each, according to fast fourier transform algorithm, product calculation is carried out to the matrix of geological data and vector, asks for the frequency-wavenumber domain data after five dimension interpolation;

Five dimension interpolated data conversion modules, for tieing up the frequency-wavenumber domain data transformation after interpolation by five to frequency field and spatial domain, then transform to time domain and spatial domain.

Compared with prior art, the technical scheme that the embodiment of the present invention provides, by by irregular earthquake data transformation to frequency field and spatial domain; Each is transformed to the frequency slice of the geological data of frequency field and spatial domain, according to fast fourier transform algorithm, product calculation is carried out to the matrix of geological data and vector, asks for the frequency-wavenumber domain data after five dimension interpolation; By the frequency-wavenumber domain data transformation after five dimension interpolation to frequency field and spatial domain, transform to time domain and spatial domain again, achieve the fast processing to irregular geological data, whole computation process only relates to a small amount of NFFT computing, matrix and the vector product computing of the overwhelming majority can utilize FFT to realize, substantially increase counting yield, improve the efficiency of seismic data process, there is important actual application value.

Accompanying drawing explanation

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

Fig. 1 is the schematic flow sheet of five dimension interpolation process methods of irregular geological data in the embodiment of the present invention;

Fig. 2 is the geological data recording geometry schematic diagram in the embodiment of the present invention before interpolation;

Fig. 3 is the geological data recording geometry schematic diagram of the same CMP bin in the embodiment of the present invention before interpolation;

Fig. 4 be after interpolation corresponding with Fig. 3 in the embodiment of the present invention the geological data recording geometry schematic diagram of same CMP bin;

Fig. 5 adopts comparison diagram operation time of different fast algorithm to represent intention in the embodiment of the present invention;

Fig. 6 is interpolation Qian CMP road collection geological data schematic diagram corresponding with Fig. 3 in the embodiment of the present invention;

Fig. 7 is interpolation Hou CMP road collection geological data corresponding with Fig. 4 in the embodiment of the present invention;

Fig. 8 is the structural representation of five dimension interpolation processors of irregular geological data in the embodiment of the present invention.

Embodiment

For making the object, technical solutions and advantages of the present invention clearly understand, below in conjunction with embodiment and accompanying drawing, the present invention is described in further details.At this, exemplary embodiment of the present invention and illustrating for explaining the present invention, but not as a limitation of the invention.

The present invention proposes a kind of new iterative scheme, whole computation process only relates to a small amount of NFFT computing, and matrix and the vector product computing of the overwhelming majority can utilize FFT to realize, and substantially increase counting yield, have important actual application value.Describe in detail as follows below in conjunction with accompanying drawing 1 to 8 pair of program.

Fig. 1 is the schematic flow sheet of five dimension interpolation process methods of irregular geological data in the embodiment of the present invention, and as shown in Figure 1, this disposal route comprises the steps:

Step 101: by irregular earthquake data transformation to frequency field and spatial domain;

Step 102: the frequency slice each being transformed to the geological data of frequency field and spatial domain, according to fast fourier transform algorithm, carries out product calculation to the matrix of geological data and vector, asks for the frequency-wavenumber domain data after five dimension interpolation;

Step 103: by the frequency-wavenumber domain data transformation after five dimension interpolation to frequency field and spatial domain, then transform to time domain and spatial domain.

Compared with prior art, the technical scheme that the embodiment of the present invention provides, by by irregular earthquake data transformation to frequency field and spatial domain; Each is transformed to the frequency slice of the geological data of frequency field and spatial domain, according to fast fourier transform algorithm, product calculation is carried out to the matrix of geological data and vector, asks for the frequency-wavenumber domain data after five dimension interpolation; By the frequency-wavenumber domain data transformation after five dimension interpolation to frequency field and spatial domain, transform to time domain and spatial domain again, achieve the fast processing to irregular geological data, whole computation process only relates to a small amount of NFFT computing, matrix and the vector product computing of the overwhelming majority can utilize FFT to realize, substantially increase counting yield, improve the efficiency of seismic data process, there is important actual application value.

In one embodiment, can comprise before above-mentioned steps 101:

Irregular geological data carries out denoising, static correction and linear NMO process;

The irregular geological data treating five dimension interpolation processing is extracted from the irregular geological data after process;

During concrete enforcement, irregular geological data carries out denoising, static correction and linear NMO process, can comprise: carry out preceding processing operations to the geological data gathered, comprise denoising, static correction, linear NMO etc.

During concrete enforcement, from the irregular geological data after process, extract the irregular geological data treating five dimension interpolation processing, can comprise: extract the geological data of wanting to carry out five dimension interpolation

Above-mentioned t _itfor time-sampling coordinate, subscript it={0,1,2 ..., N _t-1}, namely time orientation has N _tindividual sampled point, if time sampling interval is Δ t, then t _it=it Δ t;

Above-mentioned for spatial sampling coordinate, subscript ix={0,1,2 ..., N _x-1}, namely direction in space has N _xindividual sampled point, ${\tilde{x}}_{ix}={({\tilde{x}}_0,{\tilde{x}}_1,{\tilde{x}}_2,{\tilde{x}}_3)}_{ix}=({\left({\tilde{x}}_0\right)}_{ix},{\left({\tilde{x}}_1\right)}_{ix},{\left({\tilde{x}}_2\right)}_{ix},{\left({\tilde{x}}_3\right)}_{ix})$ A point in four-dimentional space, wherein be a kind of general representation, according to actual needs, shot point x coordinate, shot point y coordinate, geophone station x coordinate, geophone station y coordinate can be represented respectively, also can represent central point x coordinate, central point y coordinate respectively, geophone offset, position angle.

Above-mentioned steps 101 can comprise: treat that the irregular earthquake data transformation of five dimension interpolation processing is to frequency field and spatial domain by what extract.

During concrete enforcement, treat that the irregular earthquake data transformation of five dimension interpolation processing can comprise to frequency field and spatial domain by what extract:

Right in every one-dimensional coordinate be normalized conversion, transformation results is designated as x _ix=((x ₀) _ix, (x ₁) _ix, (x ₂) _ix, (x ₃) _ix), the geological data after conversion is expressed as

Above-mentioned normalization transform method is wherein, η={ 0,1,2,3} represents different dimensions, for any one-dimensional coordinate ix={0,1,2 ..., N _x-1}, gets its maximal value and minimum value, is designated as vmax respectively _ηand vmin _η, particularly ${\mathrm{vmax}}_{\mathrm{\η}}=\underset{ix=0,1,2,...,N_x-1}{max}\left({\left({\tilde{x}}_{\mathrm{\η}}\right)}_{ix}\right)$ With ${\mathrm{vmin}}_{\mathrm{\η}}=\underset{ix=0,1,2,...,N_x-1}{min}\left({\left({\tilde{x}}_{\mathrm{\η}}\right)}_{ix}\right),$ N _ηbe according to actual needs and setting interpolation after coordinate count, after interpolation, the sampling interval of each space dimension is $\mathrm{\Δ}{\tilde{x}}_{\mathrm{\η}}=\frac{{\mathrm{vmax}}_{\mathrm{\η}}-{\mathrm{vmin}}_{\mathrm{\η}}}{n_{\mathrm{\η}}-1}.$

Utilize Fast Fourier Transform (FFT) by geological data d (x _ix, t _it) transform to frequency field d (x _ix, ω _{i ω}).

Above-mentioned d (x _ix, ω _{i ω}) and d (x _ix, t _it) meet following relational expression, wherein, ω _{i ω}=i ω Δ ω, i ω=0,1,2 ..., N _t-1}.Here need to ensure geological data d (x _ix, t _it) meet sampling thheorem, if geological data d is (x _ix, t _it) useful signal energy be just in frequency band in, wherein then demand fulfillment condition i ω _max≤ floor (N _t/ 2).

Generate the discrete inverse Fourier transform matrix F of unequal interval (x _ix, k _ik), wherein, ix={0,1,2 ..., N _x-1} is the line index of matrix, ik={0,1,2 ..., n ₀n ₁n ₂n ₃-1} is matrix column index.

Above-mentioned wherein, k _ik=(k ₀(i ₀), k ₁(i ₁), k ₂(i ₂), k ₃(i ₃)) _ika point in four-dimentional space, k _η(i _η)=i _η-c _η, c _η=floor ((n _η-1)/2), i _η=0,1,2 ..., n _η-1}, η=0,1,2,3, function f loor () expression rounds downwards, k _ikx _ixrepresent four dimensional vector k _ikwith x _ixinner product.And ik and i ₀, i ₁, i ₂, i ₃keep relational expression ik=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃, wherein, n ₀, n ₁, n ₂, n ₃definition see step (3).

In one embodiment, above-mentioned steps 102 can comprise:

To each frequency slice i ω=0,1,2 ..., i ω _max, process according to following formula:

b _iω＝F ^Hd _iω；

y _iω＝CG(b _iω,w _iω)；

u _iω＝conj(y _iω)·*y _iω；

w _iω+1＝(sqrt(u _iω/sum(u _iω))+w _iω)/2；

Wherein, i ω _maxbe the upper frequency limit index chosen according to actual needs, selection principle makes geological data d (x _ix, t _it) useful signal energy be just in frequency band in;

W ₀be one and comprise n ₀n ₁n ₂n ₃the column vector of individual element, and w ₀each element equal wherein, w ₀w _{i ω}situation during middle subscript i ω=0;

F=F (x _ix, k _ik), be a N _xrow, n ₀n ₁n ₂n ₃the matrix of row; Wherein, F (x _ix, k _ik) be the discrete inverse Fourier transform matrix of unequal interval, ix={0,1,2 ..., N _x-1} is the line index of matrix, ik={0,1,2 ..., n ₀n ₁n ₂n ₃-1} is matrix column index, wherein, k _ik=(k ₀(i ₀), k ₁(i ₁), k ₂(i ₂), k ₃(i ₃)) _ika point in four-dimentional space, k _η(i _η)=i _η-c _η, c _η=floor ((n _η-1)/2), i _η=0,1,2 ..., n _η-1}, η=0,1,2,3, function f loor () expression rounds downwards, k _ikx _ixrepresent four dimensional vector k _ikwith x _ixinner product, ik and i ₀, i ₁, i ₂, i ₃keep relational expression ik=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃;

$d_{i\mathrm{\ω}}={\[d(x_0,{\mathrm{\ω}}_{i\mathrm{\ω}}),d(x_1,{\mathrm{\ω}}_{i\mathrm{\ω}}),d(x_2,{\mathrm{\ω}}_{i\mathrm{\ω}}),...,d(x_{N_x-1},{\mathrm{\ω}}_{i\mathrm{\ω}})\]}^H,$ Be one and comprise N _xthe column vector of individual element;

Function conj () expression asks complex conjugate to each element of input vector, and exports a vector; Operational symbol * represents that two vectorial corresponding elements are multiplied, and obtains a vector; Function sqrt () expression is extracted square root to each element of input vector, and exports a vector; Function sum () expression is sued for peace to all elements of input vector, and exports a scalar; Operational symbol ^hrepresent the conjugate transpose asking matrix;

Function y=CG (b, w) is expressed as follows preconditioned conjugate gradient method, 1. sets constant, comprising: iteration ends relative tolerance tol=10 ^-4, maximum iteration time maxit=30, A=F ^hf; 2. calculate iteration initializaing variable, comprising: r ₀=b, ρ _-1=1; 3. iteration ends absolute tolerance is calculated for i={0,1,2 ..., maxit-1}, performs and processes operation as follows:

If or i=maxit, jumps out circulation, returns y=y _i, otherwise continue circulation:

z _i＝w·*r _i；

ρ _i＝r _i ^Hz _i；

q _i＝Ap _i；

α _i＝ρ _i/p _i ^Hq _i；

y _i+1＝y _i+α _ip _i；

r _i+1＝r _i-α _iq _i；

p _i+1＝z _i+(ρ _i/ρ _i-1)p _i。

In one embodiment, in described function y=CG (b, w), each iteration needs calculating matrix and vector product computing q=Ap, wherein A=F ^hf has Multilevel Block Toeplitz matrix structure, adopts following fast algorithm to calculate q=Ap:

1. N is established ₀=n ₁n ₂n ₃, N ₁=n ₂n ₃, N ₂=n ₃, N ₃=1, m _η=2n _η-1, η=0,1,2,3, generating length is as follows m _ηvector with wherein η=0,1,2,3:

2. vectorial a={a (i) is established | i=i _a(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., m _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, wherein, and i _a(i ₀, i ₁, i2, i ₃)=((i ₀m ₁+ i ₁) m ₂+ i ₂) m ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing, as follows to the element assignment of a, initializing variable

For i ₀=0,1,2 ..., m ₀-1}, performs following circulation (1) to (12)

$\left(1\right)---t_1^r=t_0^r+i_0^r\left(i_0\right)$

$\left(2\right)---t_1^c=t_0^c+i_0^c\left(i_0\right)$

(3) for i ₁=0,1,2 ..., m ₁-1}, performs following circulation (4) to (12)

$\left(4\right)---t_2^r=t_1^r+i_1^r\left(i_1\right)$

$\left(5\right)---t_2^c=t_1^c+i_1^c\left(i_1\right)$

(6) for i ₂=0,1,2 ..., m ₂-1}, performs following circulation (7) to (12)

$\left(7\right)---t_3^r=t_2^r+i_2^r\left(i_2\right)$

$\left(8\right)---t_3^c=t_2^c+i_2^c\left(i_2\right)$

(9) for i ₃=0,1,2 ..., m ₃-1}, performs following circulation (10) to (12)

$\left(10\right)---t_4^r=t_3^r+i_3^r\left(i_3\right)$

$\left(11\right)---t_4^c=t_3^c+i_3^c\left(i_3\right)$

$\left(12\right)---a\left(i_a\right(i_0,i_1,i_2,i_3\left)\right)=A(t_4^r,t_4^c);$

Wherein, i (i) represents i-th element of vectorial i, and a (i) represents i-th element of vectorial a, A (i _r, i _c) representing matrix A i-th _rrow, i _cthe element of row;

3. four-dimensional Fast Fourier Transform (FFT) is done to the four-dimensional array represented by vectorial a, and return vector

Wherein, vector $\hat{a}=\left\{\hat{a}\left(i\right)\right|i=i_{\hat{a}}(i_0,i_1,i_2,i_3),i_{\mathrm{\η}}=0,1,2,...,m_{\mathrm{\η}}-1,\mathrm{\η}=0,1,2,3\}$ That the vectorization of a four-dimensional array represents, $i_{\hat{a}}(i_0,i_1,i_2,i_3)=\left(\right(i_0\·m_1+i_1)\·m_2+i_2)\·m_3+i_3,$ I is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing, function F FT ₄() expression does four-dimensional Fast Fourier Transform (FFT) to four-dimensional array, and returns a four-dimensional array;

4. vectorial p={p (i) is established | i=i _p(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., n _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, wherein, and i _p(i ₀, i ₁, i ₂, i ₃)=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing; If vectorial p '=and p ' (i) | i=i _{p '}(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., m _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, wherein, and i _{p '}(i ₀, i ₁, i ₂, i ₃)=((i ₀m ₁+ i ₁) m ₂+ i ₂) m ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing; In the following manner to each element assignment of vectorial p ':

$p^{\′}\left(i_{p^{\′}}\right(i_0,i_1,i_2,i_3\left)\right)=\left\{\begin{array}{cc}p\left(i_p\right(i_0,i_1,i_2,i_3))& 0\≤i_{\mathrm{\η}}\≤n_{\mathrm{\η}}-1\\ 0& n_{\mathrm{\η}}\≤i_{\mathrm{\η}}\≤m_{\mathrm{\η}}-1\end{array}\right.;$

5. calculate wherein, function F FT ₄() and IFFT ₄() represents respectively and does four-dimensional Fast Fourier Transform (FFT) and inverse transformation to four-dimensional array, and returns a four-dimensional array; Operational symbol * represents that the corresponding element of two four-dimensional arrays is multiplied, and obtains a four-dimensional array;

Vector q '=and q ' (i) | i=i _{q '}(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., m _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, wherein, and i _{q '}(i ₀, i ₁, i ₂, i ₃)=((i ₀m ₁+ i ₁) m ₂+ i ₂) m ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing;

6. vectorial q={q (i) is established | i=i _q(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., n _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, wherein, and i _q(i ₀, i ₁, i ₂, i ₃)=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing, in the following manner to each element assignment of vectorial q, q (i _q(i ₀, i ₁, i ₂, i ₃))=q ' (i _{q '}(i ₀+ n ₀, i ₁+ n ₁, i ₂+ n ₂, i ₃+ n ₃)), wherein 0≤i _η≤ n _η-1, η=0,1,2,3.

During concrete enforcement, in whole method flow, above the 1. ~ 3. step only need to calculate once, operand can be ignored substantially, therefore, as the 5. as described in step, the calculated amount of matrix and vector product computing q=Ap mainly comprises once four-dimensional Fast Fourier Transform (FFT) and once four-dimensional inverse fast Fourier transform.As shown in Figure 5, the technical scheme that the present inventor's embodiment provides improves the efficiency of irregular seismic data process.

In one embodiment, also comprise: to y _{i ω}perform function wherein, i ω=0,1,2 ..., i ω _max;

Function represent and the element position of four-dimensional array y is adjusted, and return a four-dimensional array wherein:

Y={y (i) | i=i _y(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., n _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, i _y(i ₀, i ₁, i ₂, i ₃)=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing;

$\hat{y}=\left\{\hat{y}\left(i\right)\right|i=i_{\hat{y}}(i_0,i_1,i_2,i_3),i_{\mathrm{\η}}=0,1,2,...,n_{\mathrm{\η}}-1,\mathrm{\η}=0,1,2,3\}$ That the vectorization of a four-dimensional array represents, $i_{\hat{y}}(i_0,i_1,i_2,i_3)=\left(\right(i_0\·n_1+i_1)\·n_2+i_2)\·n_3+i_3,$ I is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing;

$\hat{y}\left(i_{\hat{y}}\right({\hat{i}}_0,{\hat{i}}_1,{\hat{i}}_2,{\hat{i}}_3))=y\left(i_y\right(i_0,i_1,i_2,i_3\left)\right),$ Wherein ${\hat{i}}_{\mathrm{\η}}=\mathrm{mod}(i_{\mathrm{\η}}-c_{\mathrm{\η}},n_{\mathrm{\η}}),$ C _η=floor ((n _η-1)/2), i _η=0,1,2 ..., n _η-1}, η=0,1,2,3, function mod (v ₁, v ₂) represent v ₁to v ₂ask 0 ~ v ₂remainder between-1, function f loor () expression rounds downwards.

In one embodiment, also comprise: right perform function wherein, i ω=0,1,2 ..., i ω _max;

Function represent four-dimensional array do four-dimensional inverse fast Fourier transform, and return a four-dimensional array s, wherein,

$\hat{y}=\left\{\hat{y}\left(i\right)\right|i=i_{\hat{y}}(i_0,i_1,i_2,i_3),i_{\mathrm{\η}}=0,1,2,...,n_{\mathrm{\η}}-1,\mathrm{\η}=0,1,2,3\}$ That the vectorization of a four-dimensional array represents, $i_{\hat{y}}(i_0,i_1,i_2,i_3)=\left(\right(i_0\·n_1+i_1)\·n_2+i_2)\·n_3+i_3,$ I is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing;

S={s (i) | i=i _s(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., n _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, i _s(i ₀, i ₁, i ₂, i ₃)=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing.

In one embodiment, above-mentioned steps 103 can comprise:

Utilize s _{i ω}generate frequency field five dimension data after interpolation wherein, for the spatial sampling coordinate after interpolation, ix={0,1,2 ..., n ₀n ₁n ₂n ₃-1}, i ω=0,1,2 ..., N _t-1},

${\hat{x}}_{ix}={\left(x_0\right(i_0),x_1(i_1),x_2(i_2),x_3(i_3\left)\right)}_{ix},$ Wherein, $x_{\mathrm{\η}}\left(i_{\mathrm{\η}}\right)={\mathrm{vmin}}_{\mathrm{\η}}+i_{\mathrm{\η}}\mathrm{\Δ}{\tilde{x}}_{\mathrm{\η}},$ I _η=0,1,2 ..., n _η-1}, η=0,1,2,3, for the sampling interval of each space dimension after interpolation, ix and i ₀, i ₁, i ₂, i ₃keep relational expression ix=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃; generating mode is as follows, with vector representative a part of data,

${\hat{d}}_{i\mathrm{\ω}}=\{\hat{d}({\hat{x}}_0,{\mathrm{\ω}}_{i\mathrm{\ω}}),\hat{d}({\hat{x}}_1,{\mathrm{\ω}}_{i\mathrm{\ω}}),\hat{d}({\hat{x}}_2,{\mathrm{\ω}}_{i\mathrm{\ω}}),...,\hat{d}({\hat{x}}_{n_0{n}_1{n}_2{n}_3-1},{\mathrm{\ω}}_{i\mathrm{\ω}})\};$

${\hat{d}}_{i\mathrm{\ω}}=\{\hat{d}({\hat{x}}_0,{\mathrm{\ω}}_{i\mathrm{\ω}}),\hat{d}({\hat{x}}_1,{\mathrm{\ω}}_{i\mathrm{\ω}}),\hat{d}({\hat{x}}_2,{\mathrm{\ω}}_{i\mathrm{\ω}}),...,\hat{d}({\hat{x}}_{n_0{n}_1{n}_2{n}_3-1},{\mathrm{\ω}}_{i\mathrm{\ω}})\},$

Then ${\hat{d}}_{i\mathrm{\&omega;}}=\left\{\begin{array}{cc}{s}_{i\mathrm{\&omega;}}& 0\&le;i\mathrm{\&omega;}\&le;{\mathrm{i\&omega;}}_{max}\\ conj\left(s_{N_t-i\mathrm{\&omega;}+1}\right)& N_t-{\mathrm{i\&omega;}}_{max}\&le;i\mathrm{\&omega;}\&le;N_t\\ 0& {\mathrm{i\&omega;}}_{max}<i\mathrm{\&omega;}<N_t-{\mathrm{i\&omega;}}_{\max }\end{array}\right.;$

Wherein, function conj () expression asks complex conjugate to each element of input vector, and 0 represents by n ₀n ₁n ₂n ₃the vector that individual 0 element is formed;

Inverse fast Fourier transform is utilized to incite somebody to action the time domain of conversion, obtains five dimension data after interpolation

Wherein, with meet following relational expression, $\hat{d}\left({\hat{x}}_{ix},t_{it}\right)=\underset{i\mathrm{\&omega;}=0}{\overset{N_t-1}{\mathrm{\&Sigma;}}}\hat{d}\left({\hat{x}}_{ix},{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}\right)e^{j2\mathrm{\&pi;}\frac{\left(i\mathrm{\&omega;}\right)\left(it\right)}{N_t}},$ Wherein, t _it=it Δ t, it={0,1,2 ..., N _t-1}, ix={0,1,2 ..., n ₀n ₁n ₂n ₃-1}.

Based on same inventive concept, additionally provide a kind of five dimension interpolation processors of irregular geological data in the embodiment of the present invention, as described in the following examples.It is similar that interpolation process method tieed up by principle and five of irregular geological data due to five dimension interpolation processor problems of irregular geological data, therefore the enforcement of five dimension interpolation processors of irregular geological data see the enforcement of five dimension interpolation process methods of irregular geological data, can repeat part and repeats no more.Following used, term " unit " or " module " can realize the software of predetermined function and/or the combination of hardware.Although the device described by following examples preferably realizes with software, hardware, or the realization of the combination of software and hardware also may and conceived.

Fig. 8 is the structural representation of five dimension interpolation processors of irregular geological data in the embodiment of the present invention, and as shown in Figure 8, this device comprises:

Geological data conversion module 10, for by irregular earthquake data transformation to frequency field and spatial domain;

Five dimension interpolation processing modules 20, for transforming to the frequency slice of the geological data of frequency field and spatial domain for each, according to fast fourier transform algorithm, product calculation is carried out to the matrix of geological data and vector, asks for the frequency-wavenumber domain data after five dimension interpolation;

Five dimension interpolated data conversion modules 30, for tieing up the frequency-wavenumber domain data transformation after interpolation by five to frequency field and spatial domain, then transform to time domain and spatial domain.

In an example, also comprise:

Geological data pretreatment module, for carrying out denoising, static correction and linear NMO process to irregular geological data;

Geological data abstraction module, for extracting the irregular geological data treating five dimension interpolation processing from the irregular geological data after process;

Described geological data conversion module specifically for: by extract treat that the irregular earthquake data transformation of five dimension interpolation processing is to frequency field and spatial domain.

In an example, five dimension interpolation processing modules specifically for:

To each frequency slice i ω=0,1,2 ..., i ω _max, process according to following formula:

b _iω＝F ^Hd _iω；

y _iω＝CG(b _iω,w _iω)；

u _iω＝conj(y _iω)·*y _iω；

w _iω+1＝(sqrt(u _iω/sum(u _iω))+w _iω)/2；

Wherein, i ω _maxbe the upper frequency limit index chosen according to actual needs, selection principle makes geological data d (x _ix, t _it) useful signal energy be just in frequency band in, wherein

W ₀be one and comprise n ₀n ₁n ₂n ₃the column vector of individual element, and w ₀each element equal wherein, w ₀w _{i ω}situation during middle subscript i ω=0;

F=F (x _ix, k _ik), be a N _xrow, n ₀n ₁n ₂n ₃the matrix of row; Wherein, F (x _ix, k _ik) be the discrete inverse Fourier transform matrix of unequal interval, ix={0,1,2 ..., N _x-1} is the line index of matrix, ik={0,1,2 ..., n ₀n ₁n ₂n ₃-1} is matrix column index, wherein, k _ik=(k ₀(i ₀), k ₁(i ₁), k ₂(i ₂), k ₃(i ₃)) _ika point in four-dimentional space, k _η(i _η)=i _η-c _η, c _η=floor ((n _η-1)/2), i _η=0,1,2 ..., n _η-1}, η=0,1,2,3, function f loor () expression rounds downwards, k _ikx _ixrepresent four dimensional vector k _ikwith x _ixinner product, ik and i ₀, i ₁, i ₂, i ₃keep relational expression ik=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃;

$d_{i\mathrm{\&omega;}}={\&lsqb;d\left(x_0,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}\right),d\left(x_1,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}\right),d\left(x_2,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}\right),\mathrm{...},d\left(x_{N_x-1},{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}\right)\&rsqb;}^H,$ Be one and comprise N _xthe column vector of individual element;

Function conj () expression asks complex conjugate to each element of input vector, and exports a vector; Operational symbol * represents that two vectorial corresponding elements are multiplied, and obtains a vector; Function sqrt () expression is extracted square root to each element of input vector, and exports a vector; Function sum () expression is sued for peace to all elements of input vector, and exports a scalar; Operational symbol ^hrepresent the conjugate transpose asking matrix;

Function y=CG (b, w) is expressed as follows preconditioned conjugate gradient method, 1. sets constant, comprising: iteration ends relative tolerance tol=10 ^-4, maximum iteration time maxit=30, A=F ^hf; 2. calculate iteration initializaing variable, comprising: r ₀=b, ρ _-1=1; 3. iteration ends absolute tolerance is calculated for i={0,1,2 ..., maxit-1}, performs and processes operation as follows:

If or i=maxit, jumps out circulation, returns y=y _i, otherwise continue circulation:

z _i＝w·*r _i；

ρ _i＝r _i ^Hz _i；

q _i＝Ap _i；

α _i＝ρ _i/p _i ^Hq _i；

y _i+1＝y _i+α _ip _i；

r _i+1＝r _i-α _iq _i；

p _i+1＝z _i+(ρ _i/ρ _i-1)p _i。

In an example, in described function y=CG (b, w), each iteration needs calculating matrix and vector product computing q=Ap, wherein A=F ^hf has Multilevel Block Toeplitz matrix structure, adopts following fast algorithm to calculate q=Ap:

1. N is established ₀=n ₁n ₂n ₃, N ₁=n ₂n ₃, N ₂=n ₃, N ₃=1, m _η=2n _η-1, η=0,1,2,3, generating length is as follows m _ηvector with wherein η=0,1,2,3:

2. vectorial a={a (i) is established | i=i _a(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., m _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, wherein, and i _a(i ₀, i ₁, i ₂, i ₃)=((i ₀m ₁+ i ₁) m ₂+ i ₂) m ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing, as follows to the element assignment of a, initializing variable

For i ₀=0,1,2 ..., m ₀-1}, performs following circulation (1) to (12)

$\left(1\right)---t_1^r=t_0^r+i_0^r\left(i_0\right)$

$\left(2\right)---t_1^c=t_0^c+i_0^c\left(i_0\right)$

(3) for i ₁=0,1,2 ..., m ₁-1}, performs following circulation (4) to (12)

$\left(4\right)---t_2^r=t_1^r+i_1^r\left(i_1\right)$

$\left(5\right)---t_2^c=t_1^c+i_1^c\left(i_1\right)$

(6) for i ₂=0,1,2 ..., m ₂-1}, performs following circulation (7) to (12)

$\left(7\right)---t_3^r=t_2^r+i_2^r\left(i_2\right)$

$\left(8\right)---t_3^c=t_2^c+i_2^c\left(i_2\right)$

(9) for i ₃=0,1,2 ..., m ₃-1}, performs following circulation (10) to (12)

$\left(10\right)---t_4^r=t_3^r+i_3^r\left(i_3\right)$

$\left(11\right)---t_4^c=t_3^c+i_3^c\left(i_3\right)$

$\left(12\right)---a\left(i_a\right(i_0,i_1,i_2,i_3\left)\right)=A(t_4^r,t_4^c);$

Wherein, i (i) represents i-th element of vectorial i, and a (i) represents i-th element of vectorial a, A (i _r, i _c) representing matrix A i-th _rrow, i _cthe element of row;

3. four-dimensional Fast Fourier Transform (FFT) is done to the four-dimensional array represented by vectorial a, and return vector

Wherein, vector $\hat{a}=\left\{\hat{a}\left(i\right)\right|i=i_{\hat{a}}(i_0,i_1,i_2,i_3),i_{\mathrm{\&eta;}}=0,1,2,...,m_{\mathrm{\&eta;}}-1,\mathrm{\&eta;}=0,1,2,3\}$ That the vectorization of a four-dimensional array represents, $i_{\hat{a}}(i_0,i_1,i_2,i_3)=\left(\right(i_0\&CenterDot;m_1+i_1)\&CenterDot;m_2+i_2)\&CenterDot;m_3+i_3,$ I is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing, function F FT ₄() expression does four-dimensional Fast Fourier Transform (FFT) to four-dimensional array, and returns a four-dimensional array;

4. vectorial p={p (i) is established | i=i _p(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., n _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, wherein, and i _p(i ₀, i ₁, i ₂, i ₃)=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing; If vectorial p '=and p ' (i) | i=i _{p '}(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., m _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents _,wherein, i _{p '}(i ₀, i ₁, i ₂, i ₃)=((i ₀m ₁+ i ₁) m ₂+ i ₂) m ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing; In the following manner to each element assignment of vectorial p ':

$p^{\&prime;}\left(i_{p^{\&prime;}}\right(i_0,i_1,i_2,i_3\left)\right)=\left\{\begin{array}{cc}p\left(i_p\right(i_0,i_1,i_2,i_3))& 0\&le;i_{\mathrm{\&eta;}}\&le;n_{\mathrm{\&eta;}}-1\\ 0& n_{\mathrm{\&eta;}}\&le;i_{\mathrm{\&eta;}}\&le;m_{\mathrm{\&eta;}}-1\end{array}\right.;$

5. calculate wherein, function F FT ₄() and IFFT ₄() represents respectively and does four-dimensional Fast Fourier Transform (FFT) and inverse transformation to four-dimensional array, and returns a four-dimensional array; Operational symbol * represents that the corresponding element of two four-dimensional arrays is multiplied, and obtains a four-dimensional array;

Vector q '=and q ' (i) | i=i _{q '}(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., m _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, wherein, and i _{q '}(i ₀, i ₁, i ₂, i ₃)=((i ₀m ₁+ i ₁) m ₂+ i ₂) m ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing;

6. vectorial q={q (i) is established | i=i _q(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., n _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, wherein, and i _q(i ₀, i ₁, i ₂, i ₃)=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing, in the following manner to each element assignment of vectorial q, q (i _q(i ₀, i ₁, i ₂, i ₃))=q ' (i _{q '}(i ₀+ n ₀, i ₁+ n ₁, i ₂+ n ₂, i ₃+ n ₃)), wherein 0≤i _η≤ n _η-1, η=0,1,2,3.

In an example, also comprise: to y _{i ω}perform function wherein, i ω=0,1,2 ..., i ω _max;

Function represent and the element position of four-dimensional array y is adjusted, and return a four-dimensional array wherein:

Y={y (i) | i=i _y(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., n _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, i _y(i ₀, i ₁, i ₂, i ₃)=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing;

$\hat{y}=\left\{\hat{y}\left(i\right)\right|i=i_{\hat{y}}(i_0,i_1,i_2,i_3),i_{\mathrm{\&eta;}}=0,1,2,...,n_{\mathrm{\&eta;}}-1,\mathrm{\&eta;}=0,1,2,3\}$ That the vectorization of a four-dimensional array represents, $i_{\hat{y}}(i_0,i_1,i_2,i_3)=\left(\right(i_0\&CenterDot;n_1+i_1)\&CenterDot;n_2+i_2)\&CenterDot;n_3+i_3,$ I is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing;

$\hat{y}\left(i_{\hat{y}}\right({\hat{i}}_0,{\hat{i}}_1,{\hat{i}}_2,{\hat{i}}_3\left)\right)=y\left(i_y\right(i_0,i_1,i_2,i_3\left)\right),$ Wherein ${\hat{i}}_{\mathrm{\&eta;}}=\mathrm{mod}(i_{\mathrm{\&eta;}}-c_{\mathrm{\&eta;}},n_{\mathrm{\&eta;}}),$ C _η=floor ((n _η-1)/2), i _η=0,1,2 ..., n _η-1}, η=0,1,2,3, function mod (v ₁, v ₂) represent v ₁to v ₂ask 0 ~ v ₂remainder between-1, function f loor () expression rounds downwards.

In an example, also comprise: right perform function wherein, i ω=0,1,2 ..., i ω _max;

Function represent four-dimensional array do four-dimensional inverse fast Fourier transform, and return a four-dimensional array s, wherein,

$\hat{y}=\left\{\hat{y}\left(i\right)\right|i=i_{\hat{y}}(i_0,i_1,i_2,i_3),i_{\mathrm{\&eta;}}=0,1,2,...,n_{\mathrm{\&eta;}}-1,\mathrm{\&eta;}=0,1,2,3\}$ That the vectorization of a four-dimensional array represents, $i_{\hat{y}}(i_0,i_1,i_2,i_3)=\left(\right(i_0\&CenterDot;n_1+i_1)\&CenterDot;n_2+i_2)\&CenterDot;n_3+i_3,$ I is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing;

S={s (i) | i=i _s(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., n _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, i _s(i ₀, i ₁, i ₂, i ₃)=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing.

In an example, described five dimension interpolated data conversion modules specifically for:

Utilize s _{i ω}generate frequency field five dimension data after interpolation wherein, for the spatial sampling coordinate after interpolation, ix={0,1,2 ..., n ₀n ₁n ₂n ₃-1}, i ω=0,1,2 ..., N _t-1},

${\hat{x}}_{ix}={\left(x_0\right(i_0),x_1(i_1),x_2(i_2),x_3(i_3\left)\right)}_{ix},$ Wherein, $x_{\mathrm{\&eta;}}\left(i_{\mathrm{\&eta;}}\right)={\mathrm{vmin}}_{\mathrm{\&eta;}}+i_{\mathrm{\&eta;}}\mathrm{\&Delta;}{\tilde{x}}_{\mathrm{\&eta;}},i_{\mathrm{\&eta;}}=\{0,1,2,...,n_{\mathrm{\&eta;}}-1\},$ η=0,1,2,3, for the sampling interval of each space dimension after interpolation, ix and i ₀, i ₁, i ₂, i ₃keep relational expression ix=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃; generating mode is as follows, with vector representative a part of data,

${\hat{d}}_{i\mathrm{\&omega;}}=\{\hat{d}({\hat{x}}_0,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}),\hat{d}({\hat{x}}_1,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}),\hat{d}({\hat{x}}_2,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}),...,\hat{d}({\hat{x}}_{n_0{n}_1{n}_2{n}_3-1},{\mathrm{\&omega;}}_{i\mathrm{\&omega;}})\};$

${\hat{d}}_{i\mathrm{\&omega;}}=\{\hat{d}({\hat{x}}_0,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}),\hat{d}({\hat{x}}_1,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}),\hat{d}({\hat{x}}_2,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}),...,\hat{d}({\hat{x}}_{n_0{n}_1{n}_2{n}_3-1},{\mathrm{\&omega;}}_{i\mathrm{\&omega;}})\},$

Then ${\hat{d}}_{i\mathrm{\&omega;}}=\left\{\begin{array}{cc}{s}_{i\mathrm{\&omega;}}& 0\&le;i\mathrm{\&omega;}\&le;{\mathrm{i\&omega;}}_{max}\\ conj\left(s_{N_t-i\mathrm{\&omega;}+1}\right)& N_t-{\mathrm{i\&omega;}}_{max}\&le;i\mathrm{\&omega;}\&le;N_t\\ 0& {\mathrm{i\&omega;}}_{max}<i\mathrm{\&omega;}<N_t-{\mathrm{i\&omega;}}_{\max }\end{array}\right.;$

Wherein, function conj () expression asks complex conjugate to each element of input vector, and 0 represents by n ₀n ₁n ₂n ₃the vector that individual 0 element is formed;

Inverse fast Fourier transform is utilized to incite somebody to action the time domain of conversion, obtains five dimension data after interpolation

Wherein, with meet following relational expression, $\hat{d}\left({\hat{x}}_{ix},t_{it}\right)=\underset{i\mathrm{\&omega;}=0}{\overset{N_t-1}{\mathrm{\&Sigma;}}}\hat{d}\left({\hat{x}}_{ix},{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}\right)e^{j2\mathrm{\&pi;}\frac{\left(i\mathrm{\&omega;}\right)\left(it\right)}{N_t}},$ Wherein, t _it=it Δ t, it={0,1,2 ..., N _t-1}, ix={0,1,2 ..., n ₀n ₁n ₂n ₃-1}.

Be described with example more below, so that understand how to implement the present invention.

The embodiment of the present invention provides five of irregular geological data dimension interpolation process methods mainly to comprise the steps:

(1), to the geological data gathered carry out preceding processing operations, comprise denoising, static correction, linear NMO etc.;

(2), the geological data carrying out five dimension interpolation is wanted in extraction fig. 2 is the geological data layout chart before the interpolation that provides of the embodiment of the present invention.

T described in step (2) _itfor time-sampling coordinate, subscript it={0,1,2 ..., N _t-1}, namely time orientation has N _tindividual sampled point, if time sampling interval is Δ t, then t _it=it Δ t.In this example, N _t=1200, Δ t=0.002 seconds;

Described in step (2) for spatial sampling coordinate, subscript ix={0,1,2 ..., N _x-1}, namely direction in space has N _xindividual sampled point, ${\tilde{x}}_{ix}={({\tilde{x}}_0,{\tilde{x}}_1,{\tilde{x}}_2,{\tilde{x}}_3)}_{ix}=({\left({\tilde{x}}_0\right)}_{ix},{\left({\tilde{x}}_1\right)}_{ix},{\left({\tilde{x}}_2\right)}_{ix},{\left({\tilde{x}}_3\right)}_{ix})$ A point in four-dimentional space, wherein be a kind of general representation, according to these routine needs, represent central point x coordinate, central point y coordinate respectively, geophone offset, position angle.In this example, N _x=7151.

(3), right in every one-dimensional coordinate be normalized conversion, transformation results is designated as x _ix=((x ₀) _ix, (x ₁) _ix, (x ₂) _ix, (x ₃) _ix), the geological data after conversion is expressed as

Normalization transform method described in step (3) is wherein, η={ 0,1,2,3} represents different dimensions, for any one-dimensional coordinate ix={0,1,2 ..., N _x-1}, gets its maximal value and minimum value, is designated as vmax respectively _ηand vmin _η, particularly with n _ηbe according to actual needs and setting interpolation after coordinate count, after interpolation, the sampling interval of each space dimension is in this example, vmin ₀=253, vmax ₀=263, vmin ₁=660, vmax ₁=670, vmin ₂=-160, vmax ₂=160, vmin ₃=0, vmax ₃=170, n ₀=11, n ₁=11, n ₂=65, n ₃=18, $\mathrm{\&Delta;}{\tilde{x}}_0=1,\mathrm{\&Delta;}{\tilde{x}}_1=1,\mathrm{\&Delta;}{\tilde{x}}_2=0.5,\mathrm{\&Delta;}{\tilde{x}}_3=10.$

(4), utilize Fast Fourier Transform (FFT) by geological data d (x _ix, t _it) transform to frequency field d (x _ix, ω _{i ω}).

D (x described in step (4) _ix, ω _{i ω}) and d (x _ix, t _it) meet following relational expression, $d(x_{ix},{\mathrm{\&omega;}}_{i\mathrm{\&omega;}})=\underset{it=0}{\overset{N_t-1}{\&Sigma;}}d(x_{ix},t_{it})e^{-j2\mathrm{\&pi;}\frac{\left(i\mathrm{\&omega;}\right)\left(it\right)}{N_t}},$ Wherein, t _it=it Δ t, ω _{i ω}=i ω Δ ω, i ω=0,1,2 ..., N _t-1}.Here need to ensure geological data d (x _ix, t _it) meet sampling thheorem, if geological data d is (x _ix, t _it) useful signal energy be just in frequency band in, wherein then demand fulfillment condition i ω _max≤ floor (N _t/ 2).In this example, Δ ω ≈ 2.618, i ω _max=175.

(5) the discrete inverse Fourier transform matrix F of unequal interval (x, is generated _ix, k _ik), wherein, ix={0,1,2 ..., N _x-1} is the line index of matrix, ik={0,1,2 ..., n ₀n ₁n ₂n ₃-1} is matrix column index.In this example, n ₀n ₁n ₂n ₃=141570.

Described in step (5) wherein, k _ik=(k ₀(i ₀), k ₁(i ₁), k ₂(i ₂), k ₃(i ₃)) _ika point in four-dimentional space, k _η(i _η)=i _η-c _η, c _η=floor ((n _η-1)/2), i _η=0,1,2 ..., n _η-1}, η=0,1,2,3, function f loor () expression rounds downwards, k _ikx _ixrepresent four dimensional vector k _ikwith x _ixinner product.And ik and i ₀, i ₁, i ₂, i ₃keep relational expression ik=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃, wherein, n ₀, n ₁, n ₂, n ₃definition see step (3).In this example, c ₀=5, c ₁=5, c ₂=32, c ₃=8.

(6), to each frequency slice i ω=0,1,2 ..., i ω _max, perform and circulate as follows, and store y _{i ω}.

b _iω＝F ^Hd _iω

y _iω＝CG(b _iω,w _iω)

u _iω＝conj(y _iω)·*y _iω。

w _iω+1＝(sqrt(u _iω/sum(u _iω))+w _iω)/2

I ω described in step (6) _maxbe the upper frequency limit index chosen according to actual needs, selection principle makes geological data d (x _ix, t _it) useful signal energy be just in frequency band in, wherein

W described in step (6) ₀be one and comprise n ₀n ₁n ₂n ₃the column vector of individual element, and w ₀each element equal wherein, w ₀w _{i ω}situation during middle subscript i ω=0.

F=F (x described in step (6) _ix, k _ik), be a N _xrow, n ₀n ₁n ₂n ₃the matrix of row.

Described in step (6) $d_{i\mathrm{\&omega;}}={\&lsqb;d\left(x_0,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}\right),d\left(x_1,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}\right),d\left(x_2,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}\right),\mathrm{...},d\left(x_{N_x-1},{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}\right)\&rsqb;}^H,$ Be one and comprise N _xthe column vector of individual element.

Function conj () expression described in step (6) asks complex conjugate to each element of input vector, and exports a vector; Operational symbol * represents that two vectorial corresponding elements are multiplied, and obtains a vector; Function sqrt () expression is extracted square root to each element of input vector, and exports a vector; Function sum () expression is sued for peace to all elements of input vector, and exports a scalar; Operational symbol ^hrepresent the conjugate transpose asking matrix, lower same.

Function y=CG (b, w) described in step (6) is expressed as follows preconditioned conjugate gradient method, 1. sets constant.Comprise iteration ends relative tolerance tol=10 ^-4, maximum iteration time maxit=30, A=F ^hf; 2. iteration initializaing variable is calculated.Comprise r ₀=b, ρ _-1=1; 3. iteration ends absolute tolerance is calculated for i={0,1,2 ..., maxit-1}, performs and circulates as follows

If or i=maxit, jumps out circulation, returns y=y _i, otherwise continue circulation

z _i＝w·*r _i

ρ _i＝r _i ^Hz _i

q _i＝Ap _i

α _i＝ρ _i/p _i ^Hq _i

y _i+1＝y _i+α _ip _i

r _i+1＝r _i-α _iq _i

p _i+1＝z _i+(ρ _i/ρ _i-1)p _i

In function y=CG (b, w) described in step (6), each iteration needs calculating matrix and vector product computing q=Ap, wherein A=F ^hf has Multilevel Block Toeplitz matrix structure, and following fast algorithm can be adopted to calculate q=Ap.

1. N is established ₀=n ₁n ₂n ₃, N ₁=n ₂n ₃, N ₂=n ₃, N ₃=1.m _η＝2n _η-1,η＝0,1,2,3。Generating length is as follows m _ηvector with wherein η=0,1,2,3.In this example, N ₀=12870, N ₁=1170, N ₂=18, N ₃=1.m ₀＝21，m ₁＝21，m ₂＝129，m ₃＝35。

2. vectorial a={a (i) is established | i=i _a(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., m _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, wherein, and i _a(i ₀, i ₁, i ₂, i ₃)=((i ₀m ₁+ i ₁) m ₂+ i ₂) m ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing.As follows to the element assignment of a.Initializing variable

For i ₀=0,1,2 ..., m ₀-1}, performs following circulation (1) to (12)

$\left(1\right)---t_1^r=t_0^r+i_0^r\left(i_0\right)$

$\left(2\right)---t_1^c=t_0^c+i_0^c\left(i_0\right)$

(3) for i ₁=0,1,2 ..., m ₁-1}, performs following circulation (4) to (12)

$\left(4\right)---t_2^r=t_1^r+i_1^r\left(i_1\right)$

$\left(5\right)---t_2^c=t_1^c+i_1^c\left(i_1\right)$

(6) for i ₂=0,1,2 ..., m ₂-1}, performs following circulation (7) to (12)

$\left(7\right)---t_3^r=t_2^r+i_2^r\left(i_2\right)$

$\left(8\right)---t_3^c=t_2^c+i_2^c\left(i_2\right)$

(9) for i ₃=0,1,2 ..., m ₃-1}, performs following circulation (10) to (12)

$\left(10\right)---t_4^r=t_3^r+i_3^r\left(i_3\right)$

$\left(11\right)---t_4^c=t_3^c+i_3^c\left(i_3\right)$

$\left(12\right)---a\left(i_a\right(i_0,i_1,i_2,i_3\left)\right)=A(t_4^r,t_4^c);$

Wherein, i (i) represents i-th element of vectorial i, and a (i) represents i-th element of vectorial a, A (i _r, i _c) representing matrix A i-th _rrow, i _cthe element of row.

3. four-dimensional Fast Fourier Transform (FFT) is done to the four-dimensional array represented by vectorial a, and return vector wherein, vector $\hat{a}=\left\{\hat{a}\left(i\right)\right|i=i_{\hat{a}}(i_0,i_1,i_2,i_3),i_{\mathrm{\&eta;}}=0,1,2,...,m_{\mathrm{\&eta;}}-1,\mathrm{\&eta;}=0,1,2,3\}$ That the vectorization of a four-dimensional array represents, $i_{\hat{a}}(i_0,i_1,i_2,i_3)=\left(\right(i_0\&CenterDot;m_1+i_1)\&CenterDot;m_2+i_2)\&CenterDot;m_3+i_3,$ I is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing. function F FT ₄() expression does four-dimensional Fast Fourier Transform (FFT) to four-dimensional array, and returns a four-dimensional array.

4. vectorial p={p (i) is established | i=i _p(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., n _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, wherein, and i _p(i ₀, i ₁, i ₂, i ₃)=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing.If vectorial p '=and p ' (i) | i=i _{p '}(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., m _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, wherein, and i _{p '}(i ₀, i ₁, i ₂, i ₃)=((i ₀m ₁+ i ₁) m ₂+ i ₂) m ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing.In the following manner to each element assignment of vectorial p '

$p^{\&prime;}\left(i_{p^{\&prime;}}\right(i_0,i_1,i_2,i_3\left)\right)=\left\{\begin{array}{cc}p\left(i_p\right(i_0,i_1,i_2,i_3))& 0\&le;i_{\mathrm{\&eta;}}\&le;n_{\mathrm{\&eta;}}-1\\ 0& n_{\mathrm{\&eta;}}\&le;i_{\mathrm{\&eta;}}\&le;m_{\mathrm{\&eta;}}-1\end{array}\right..$

5. calculate wherein, function F FT ₄() and IFFT ₄() represents respectively and does four-dimensional Fast Fourier Transform (FFT) and inverse transformation to four-dimensional array, and returns a four-dimensional array; Operational symbol * represents that the corresponding element of two four-dimensional arrays is multiplied, and obtains a four-dimensional array; Vector q '=and q ' (i) | i=i _{q '}(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., m _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, wherein, and i _{q '}(i ₀, i ₁, i ₂, i ₃)=((i ₀m ₁+ i ₁) m ₂+ i ₂) m ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing.

6. vectorial q={q (i) is established | i=i _q(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., n _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, wherein, and i _q(i ₀, i ₁, i ₂, i ₃)=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing.In the following manner to each element assignment of vectorial q, q (i _q(i ₀, i ₁, i ₂, i ₃))=q ' (i _{q '}(i ₀+ n ₀, i ₁+ n ₁, i ₂+ n ₂, i ₃+ n ₃)), wherein 0≤i _η≤ n _η-1, η=0,1,2,3.

7. in whole method flow, above the 1. ~ 3. step only need to calculate once, operand can be ignored substantially, therefore, as the 5. as described in step, the calculated amount of matrix and vector product computing q=Ap mainly comprises once four-dimensional Fast Fourier Transform (FFT) and once four-dimensional inverse fast Fourier transform.

(7), to y _{i ω}perform function and store wherein

Function described in step (7) represent and the element position of four-dimensional array y is adjusted, and return a four-dimensional array wherein,

Y={y (i) | i=i _y(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., n _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, i _y(i ₀, i ₁, i ₂, i ₃)=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing;

$\hat{y}=\left\{\hat{y}\left(i\right)\right|i=i_{\hat{y}}(i_0,i_1,i_2,i_3),i_{\mathrm{\&eta;}}=0,1,2,...,n_{\mathrm{\&eta;}}-1,\mathrm{\&eta;}=0,1,2,3\}$ That the vectorization of a four-dimensional array represents, $i_{\hat{y}}(i_0,i_1,i_2,i_3)=\left(\right(i_0\&CenterDot;n_1+i_1)\&CenterDot;n_2+i_2)\&CenterDot;n_3+i_3,$ I is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing.

Concrete method of adjustment is as follows, $\hat{y}\left(i_{\hat{y}}\right({\hat{i}}_0,{\hat{i}}_1,{\hat{i}}_2,{\hat{i}}_3\left)\right)=y\left(i_y\right(i_0,i_1,i_2,i_3\left)\right),$ Wherein ${\hat{i}}_{\mathrm{\&eta;}}=\mathrm{mod}(i_{\mathrm{\&eta;}}-c_{\mathrm{\&eta;}},n_{\mathrm{\&eta;}}),$ C _η=floor ((n _η-1)/2), i _η=0,1,2 ..., n _η-1}, η=0,1,2,3, function mod (v ₁, v ₂) represent v ₁to v ₂ask 0 ~ v ₂remainder between-1, function f loor () expression rounds downwards.

(8), right perform function and store s _{i ω}, wherein i ω=0,1,2 ..., i ω _max.

Function described in step (8) represent four-dimensional array do four-dimensional inverse fast Fourier transform, and return a four-dimensional array s.Wherein,

$\hat{y}=\left\{\hat{y}\left(i\right)\right|i=i_{\hat{y}}(i_0,i_1,i_2,i_3),i_{\mathrm{\&eta;}}=0,1,2,...,n_{\mathrm{\&eta;}}-1,\mathrm{\&eta;}=0,1,2,3\}$ That the vectorization of a four-dimensional array represents, $i_{\hat{y}}(i_0,i_1,i_2,i_3)=\left(\right(i_0\&CenterDot;n_1+i_1)\&CenterDot;n_2+i_2)\&CenterDot;n_3+i_3,$ I is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing;

S={s (i) | i=i _s(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., n _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, i _s(i ₀, i ₁, i ₂, i ₃)=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing.

(9), s is utilized _{i ω}generate frequency field five dimension data after interpolation wherein, for the spatial sampling coordinate after interpolation, ix={0,1,2 ..., n ₀n ₁n ₂n ₃-1}, i ω=0,1,2 ..., N _t-1}.

Described in step (9) ${\hat{x}}_{ix}={\left(x_0\right(i_0),x_1(i_1),x_2(i_2),x_3(i_3\left)\right)}_{ix},$ Wherein, $x_{\mathrm{\&eta;}}\left(i_{\mathrm{\&eta;}}\right)={\mathrm{vmin}}_{\mathrm{\&eta;}}+i_{\mathrm{\&eta;}}\mathrm{\&Delta;}{\tilde{x}}_{\mathrm{\&eta;}},$ I _η=0,1,2 ..., n _η-1}, η=0,1,2,3, definition see step (3).And ix and i ₀, i ₁, i ₂, i ₃keep relational expression ix=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃.

Described in step (9) generating mode is as follows, conveniently states, with vector representative a part of data, particularly, ${\hat{d}}_{i\mathrm{\&omega;}}=\{\hat{d}({\hat{x}}_0,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}),\hat{d}({\hat{x}}_1,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}),\hat{d}({\hat{x}}_2,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}),...,\hat{d}({\hat{x}}_{n_0{n}_1{n}_2{n}_3-1},{\mathrm{\&omega;}}_{i\mathrm{\&omega;}})\},$ Then

${\hat{d}}_{i\mathrm{\&omega;}}=\left\{\begin{array}{cc}{s}_{i\mathrm{\&omega;}}& 0\&le;i\mathrm{\&omega;}\&le;{\mathrm{i\&omega;}}_{max}\\ conj\left(s_{N_t-i\mathrm{\&omega;}+1}\right)& N_t-{\mathrm{i\&omega;}}_{max}\&le;i\mathrm{\&omega;}\&le;N_t\\ 0& {\mathrm{i\&omega;}}_{max}<i\mathrm{\&omega;}<N_t-{\mathrm{i\&omega;}}_{\max }\end{array}\right.$

Wherein, function conj () expression asks complex conjugate to each element of input vector.0 represents by n ₀n ₁n ₂n ₃the vector that individual 0 element is formed.

(10) inverse fast Fourier transform, is utilized to incite somebody to action the time domain of conversion, obtains five dimension data after interpolation

Described in step (10) with meet following relational expression, $\hat{d}\left({\hat{x}}_{ix},t_{it}\right)=\underset{i\mathrm{\&omega;}=0}{\overset{N_t-1}{\mathrm{\&Sigma;}}}\hat{d}\left({\hat{x}}_{ix},{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}\right)e^{j2\mathrm{\&pi;}\frac{\left(i\mathrm{\&omega;}\right)\left(it\right)}{N_t}},$ Wherein, t _it=it Δ t, it={0,1,2 ..., N _t-1}, ix={0,1,2 ..., n ₀n ₁n ₂n ₃-1}.

Fig. 2 is the geological data recording geometry schematic diagram in the embodiment of the present invention before interpolation; Fig. 3 is the geological data recording geometry schematic diagram of the same CMP bin in the embodiment of the present invention before interpolation; Fig. 4 be after interpolation corresponding with Fig. 3 in the embodiment of the present invention the geological data recording geometry schematic diagram of same CMP bin; Fig. 5 adopts comparison diagram operation time of different fast algorithm to represent intention in the embodiment of the present invention; Fig. 6 is interpolation Qian CMP road collection geological data schematic diagram corresponding with Fig. 3 in the embodiment of the present invention; Fig. 7 is interpolation Hou CMP road collection geological data corresponding with Fig. 4 in the embodiment of the present invention.

Compared with prior art, by known with introducing of above-described embodiment shown in Fig. 2 to Fig. 7, five dimension interpolation fast algorithms of the irregular acquiring seismic data that the embodiment of the present invention provides, core proposes a kind of new iterative scheme, substantially increases counting yield.Whole computation process only relates to a small amount of NFFT computing, and matrix and the vector product computing of the overwhelming majority can utilize FFT to realize, and substantially increase counting yield, have important actual application value.

Those skilled in the art should understand, embodiments of the invention can be provided as method, system or computer program.Therefore, the present invention can adopt the form of complete hardware embodiment, completely software implementation or the embodiment in conjunction with software and hardware aspect.And the present invention can adopt in one or more form wherein including the upper computer program implemented of computer-usable storage medium (including but not limited to magnetic disk memory, CD-ROM, optical memory etc.) of computer usable program code.

The present invention describes with reference to according to the process flow diagram of the method for the embodiment of the present invention, equipment (system) and computer program and/or block scheme.Should understand can by the combination of the flow process in each flow process in computer program instructions realization flow figure and/or block scheme and/or square frame and process flow diagram and/or block scheme and/or square frame.These computer program instructions can being provided to the processor of multi-purpose computer, special purpose computer, Embedded Processor or other programmable data processing device to produce a machine, making the instruction performed by the processor of computing machine or other programmable data processing device produce device for realizing the function of specifying in process flow diagram flow process or multiple flow process and/or block scheme square frame or multiple square frame.

These computer program instructions also can be stored in can in the computer-readable memory that works in a specific way of vectoring computer or other programmable data processing device, the instruction making to be stored in this computer-readable memory produces the manufacture comprising command device, and this command device realizes the function of specifying in process flow diagram flow process or multiple flow process and/or block scheme square frame or multiple square frame.

These computer program instructions also can be loaded in computing machine or other programmable data processing device, make on computing machine or other programmable devices, to perform sequence of operations step to produce computer implemented process, thus the instruction performed on computing machine or other programmable devices is provided for the step realizing the function of specifying in process flow diagram flow process or multiple flow process and/or block scheme square frame or multiple square frame.

The foregoing is only the preferred embodiments of the present invention, be not limited to the present invention, for a person skilled in the art, the embodiment of the present invention can have various modifications and variations.Within the spirit and principles in the present invention all, any amendment done, equivalent replacement, improvement etc., all should be included within protection scope of the present invention.

Claims (14)

1. five dimension interpolation process methods of irregular geological data, is characterized in that, comprising:

By irregular earthquake data transformation to frequency field and spatial domain;

Each is transformed to the frequency slice of the geological data of frequency field and spatial domain, according to fast fourier transform algorithm, product calculation is carried out to the matrix of geological data and vector, asks for the frequency-wavenumber domain data after five dimension interpolation;

By the frequency-wavenumber domain data transformation after five dimension interpolation to frequency field and spatial domain, then transform to time domain and spatial domain.

2. five dimension interpolation process methods of irregular geological data as claimed in claim 1, is characterized in that, described by before irregular earthquake data transformation to frequency field and spatial domain, comprising:

Denoising, static correction and linear NMO process are carried out to irregular geological data;

The irregular geological data treating five dimension interpolation processing is extracted from the irregular geological data after process;

Described by irregular earthquake data transformation to frequency field and spatial domain, comprising: by extract treat that the irregular earthquake data transformation of five dimension interpolation processing is to frequency field and spatial domain.

3. five dimension interpolation process methods of irregular geological data as claimed in claim 1, it is characterized in that, each is transformed to the frequency slice of the geological data of frequency field and spatial domain, according to fast fourier transform algorithm, product calculation is carried out to the matrix of geological data and vector, ask for the frequency-wavenumber domain data after five dimension interpolation, comprising:

To each frequency slice i ω=0,1,2 ..., i ω _max, process according to following formula:

b _iω＝F ^Hd _iω；

y _iω＝CG(b _iω,w _iω)；

u _iω＝conj(y _iω)·*y _iω；

w _iω+1＝(sqrt(u _iω/sum(u _iω))+w _iω)/2；

Wherein, i ω _maxbe the upper frequency limit index chosen according to actual needs, selection principle makes geological data d (x _ix, t _it) useful signal energy be just in frequency band in;

W ₀be one and comprise n ₀n ₁n ₂n ₃the column vector of individual element, and w ₀each element equal wherein, w ₀w _{i ω}situation during middle subscript i ω=0;

F=F (x _ix, k _ik), be a N _xrow, n ₀n ₁n ₂n ₃the matrix of row; Wherein, F (x _ix, k _ik) be the discrete inverse Fourier transform matrix of unequal interval, ix={0,1,2 ..., N _x-1} is the line index of matrix, ik={0,1,2 ..., n ₀n ₁n ₂n ₃-1} is matrix column index, wherein, k _ik=(k ₀(i ₀), k ₁(i ₁), k ₂(i ₂), k ₃(i ₃)) _ika point in four-dimentional space, k _η(i _η)=i _η-c _η, c _η=floor ((n _η-1)/2), i _η=0,1,2 ..., n _η-1}, η=0,1,2,3, function f loor () expression rounds downwards, k _ikx _ixrepresent four dimensional vector k _ikwith x _ixinner product, ik and i ₀, i ₁, i ₂, i ₃keep relational expression ik=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃;

$d_{\mathrm{i\&omega;}}=[d(x_0,{\mathrm{\&omega;}}_{\mathrm{i\&omega;}}),d(x_1,{\mathrm{\&omega;}}_{\mathrm{i\&omega;}}),d(x_2,{\mathrm{\&omega;}}_{\mathrm{i\&omega;}}),...,d{(x_{N_{x^{-1}}},{\mathrm{\&omega;}}_{\mathrm{i\&omega;}})]}^H,$ Be one and comprise N _xthe column vector of individual element;

Function conj () expression asks complex conjugate to each element of input vector, and exports a vector; Operational symbol * represents that two vectorial corresponding elements are multiplied, and obtains a vector; Function sqrt () expression is extracted square root to each element of input vector, and exports a vector; Function sum () expression is sued for peace to all elements of input vector, and exports a scalar; Operational symbol ^hrepresent the conjugate transpose asking matrix;

Function y=CG (b, w) is expressed as follows preconditioned conjugate gradient method, 1. sets constant, comprising: iteration ends relative tolerance tol=10 ^-4, maximum iteration time maxit=30, A=F ^hf; 2. calculate iteration initializaing variable, comprising: r ₀=b, ρ _-1=1; 3. iteration ends absolute tolerance is calculated for i={0,1,2 ..., maxit-1}, performs and processes operation as follows:

If or i=maxit, jumps out circulation, returns y=y _i, otherwise continue circulation:

z _i＝w·*r _i；

ρ _i＝r _i ^Hz _i；

q _i＝Ap _i；

α _i＝ρ _i/p _i ^Hq _i；

y _i+1＝y _i+α _ip _i；

r _i+1＝r _i-α _iq _i；

p _i+1＝z _i+(ρ _i/ρ _i-1)p _i。

4. five dimension interpolation process methods of irregular geological data as claimed in claim 3, it is characterized in that, in described function y=CG (b, w), each iteration needs calculating matrix and vector product computing q=Ap, wherein A=F ^hf has Multilevel Block Toeplitz matrix structure, adopts following fast algorithm to calculate q=Ap:

1. N is established ₀=n ₁n ₂n ₃, N ₁=n ₂n ₃, N ₂=n ₃, N ₃=1, m _η=2n _η-1, η=0,1,2,3, generating length is as follows m _ηvector with wherein η=0,1,2,3:

2. vectorial a={a (i) is established | i=i _a(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., m _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, wherein, and i _a(i ₀, i ₁, i ₂, i ₃)=((i ₀m ₁+ i ₁) m ₂+ i ₂) m ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing, as follows to the element assignment of a, initializing variable

For i ₀=0,1,2 ..., m ₀-1}, performs following circulation (1) to (12)

(1) $t_1^r=t_0^r+i_0^r\left(i_0\right)$

(2) $t_1^c=t_0^c+i_0^c\left(i_0\right)$

(3) for i ₁=0,1,2 ..., m ₁-1}, performs following circulation (4) to (12)

(4) $t_2^r=t_1^r+i_1^r\left(i_1\right)$

(5) $t_2^c=t_1^c+i_1^c\left(i_1\right)$

(6) for i ₂=0,1,2 ..., m ₂-1}, performs following circulation (7) to (12)

(7) $t_3^r=t_2^r+i_2^r\left(i_2\right)$

(8) $t_3^c=t_2^c+i_2^c\left(i_2\right)$

(9) for i ₃=0,1,2 ..., m ₃-1}, performs following circulation (10) to (12)

(10) $t_4^r=t_3^r+i_3^r\left(i_3\right)$

(11) $t_4^c=t_3^c+i_3^c\left(i_3\right)$

(12) $a\left(i_a\right(i_0,i_1,i_2,i_3\left)\right)=A(t_4^r,t_4^c);$

Wherein, i (i) represents i-th element of vectorial i, and a (i) represents i-th element of vectorial a, A (i _r, i _c) representing matrix A i-th _rrow, i _cthe element of row;

3. four-dimensional Fast Fourier Transform (FFT) is done to the four-dimensional array represented by vectorial a, and return vector

Wherein, vector $\hat{a}=\left\{\hat{a}\left(i\right)\right|i=i_{\hat{a}}(i_0,i_1,i_2,i_3),i_{\mathrm{\&eta;}}=0,1,2,...,m_{\mathrm{\&eta;}}-1,\mathrm{\&eta;}=0,1,2,3\}$ That the vectorization of a four-dimensional array represents, $i_{\hat{a}}(i_0,i_1,i_2,i_3)=\left(\right(i_0\&CenterDot;m_1+i_1)\&CenterDot;m_2+i_2)\&CenterDot;m_3+i_3,$ I is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing, function F FT ₄() expression does four-dimensional Fast Fourier Transform (FFT) to four-dimensional array, and returns a four-dimensional array;

4. vectorial p={p (i) is established | i=i _p(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., n _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, wherein, and i _p(i ₀, i ₁, i ₂, i ₃)=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing; If vectorial p '=and p ' (i) | i=i _{p '}(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., m _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, wherein, and i _{p '}(i ₀, i ₁, i ₂, i ₃)=((i ₀m ₁+ i ₁) m ₂+ i ₂) m ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing; In the following manner to each element assignment of vectorial p ':

$p^{\&prime;}\left(i_{p^{\&prime;}}\right(i_0,i_1,i_2,i_3\left)\right)=\left\{\begin{array}{cc}p\left(i_p\right(i_0,i_1,i_2,i_3))& 0\&le;i_{\mathrm{\&eta;}}\&le;n_{\mathrm{\&eta;}}-1\\ 0& n_{\mathrm{\&eta;}}\&le;i_{\mathrm{\&eta;}}\&le;m_{\mathrm{\&eta;}}-1\end{array}\right.;$

5. calculate wherein, function F FT ₄() and IFFT ₄() represents respectively and does four-dimensional Fast Fourier Transform (FFT) and inverse transformation to four-dimensional array, and returns a four-dimensional array; Operational symbol * represents that the corresponding element of two four-dimensional arrays is multiplied, and obtains a four-dimensional array;

Vector q '=and q ' (i) | i=i _{q '}(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., m _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, wherein, and i _{q '}(i ₀, i ₁, i ₂, i ₃)=((i ₀m ₁+ i ₁) m ₂+ i ₂) m ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing;

6. vectorial q={q (i) is established | i=i _q(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., n _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, wherein, and i _q(i ₀, i ₁, i ₂, i ₃)=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing, in the following manner to each element assignment of vectorial q, q (i _q(i ₀, i ₁, i ₂, i ₃))=q ' (i _{q '}(i ₀+ n ₀, i ₁+ n ₁, i ₂+ n ₂, i ₃+ n ₃)), wherein 0≤i _η≤ n _η-1, η=0,1,2,3.

5. five dimension interpolation process methods of irregular geological data as claimed in claim 3, is characterized in that, also comprise: to y _{i ω}perform function wherein, i ω=0,1,2 ..., i ω _max;

Function represent and the element position of four-dimensional array y is adjusted, and return a four-dimensional array wherein:

Y={y (i) | i=i _y(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., n _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, i _y(i ₀, i ₁, i ₂, i ₃)=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing;

$\hat{y}=\left\{\hat{y}\left(i\right)\right|i=i_{\hat{y}}(i_0,i_1,i_2,i_3),i_{\mathrm{\&eta;}}=0,1,2,...,n_{\mathrm{\&eta;}}-1,\mathrm{\&eta;}=0,1,2,3\}$ That the vectorization of a four-dimensional array represents, $i_{\hat{y}}(i_0,i_1,i_2,i_3)=\left(\right(i_0\&CenterDot;n_1+i_1)\&CenterDot;n_2+i_2)\&CenterDot;n_3+i_3,$ I is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing;

$\hat{y}\left(i_{\hat{y}}\right({\hat{i}}_0,{\hat{i}}_1,{\hat{i}}_2,{\hat{i}}_3\left)\right)=y\left(i_y\right(i_0,i_1,i_2,i_3\left)\right),$ Wherein ${\hat{i}}_{\mathrm{\&eta;}}=\mathrm{mod}(i_{\mathrm{\&eta;}}-c_{\mathrm{\&eta;}},n_{\mathrm{\&eta;}}),$ C _η=floor ((n _η-1)/2), i _η=0,1,2 ..., n _η-1}, η=0,1,2,3, function mod (v ₁, v ₂) represent v ₁to v ₂ask 0 ~ v ₂remainder between-1, function f loor () expression rounds downwards.

6. five dimension interpolation process methods of irregular geological data as claimed in claim 5, is characterized in that, also comprise: be right perform function wherein, i ω=0,1,2 ..., i ω _max;

Function represent four-dimensional array do four-dimensional inverse fast Fourier transform, and return a four-dimensional array s, wherein,

$\hat{y}=\left\{\hat{y}\left(i\right)\right|i=i_{\hat{y}}(i_0,i_1,i_2,i_3),i_{\mathrm{\&eta;}}=0,1,2,...,n_{\mathrm{\&eta;}}-1,\mathrm{\&eta;}=0,1,2,3\}$ That the vectorization of a four-dimensional array represents, $i_{\hat{y}}(i_0,i_1,i_2,i_3)=\left(\right(i_0\&CenterDot;n_1+i_1)\&CenterDot;n_2+i_2)\&CenterDot;n_3+i_3,$ I is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing;

S={s (i) | i=i _s(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., n _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, i _s(i ₀, i ₁, i ₂, i ₃)=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing.

7. five dimension interpolation process methods of irregular geological data as claimed in claim 6, it is characterized in that, by the frequency-wavenumber domain data transformation after interpolation to frequency and spatial domain, then transform to Time and place territory, complete the five dimension interpolation processing to irregular geological data, comprising:

Utilize s _{i ω}generate frequency field five dimension data after interpolation wherein, for the spatial sampling coordinate after interpolation, ix={0,1,2 ..., n ₀n ₁n ₂n ₃-1}, i ω=0,1,2 ..., N _t-1},

${\hat{x}}_{ix}={\left(x_0\right(i_0),x_1(i_1),x_2(i_2),x_3(i_3\left)\right)}_{ix},$ Wherein, $x_{\mathrm{\&eta;}}\left(i_{\mathrm{\&eta;}}\right)={\mathrm{vmin}}_{\mathrm{\&eta;}}+i_{\mathrm{\&eta;}}\mathrm{\&Delta;}{\tilde{x}}_{\mathrm{\&eta;}},$ I _η=0,1,2 ..., n _η-1}, η=0,1,2,3, for the sampling interval of each space dimension after interpolation, ix and i ₀, i ₁, i ₂, i ₃keep relational expression $ix=\left(\right(i_0\&CenterDot;n_1+i_1)\&CenterDot;n_2+i_2)\&CenterDot;n_3+i_3;$ generating mode is as follows, with vector representative a part of data,

${\hat{d}}_{i\mathrm{\&omega;}}=\{\hat{d}({\hat{x}}_0,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}),\hat{d}({\hat{x}}_1,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}),\hat{d}({\hat{x}}_2,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}),...,\hat{d}({\hat{x}}_{n_0{n}_1{n}_2{n}_3-1},{\mathrm{\&omega;}}_{i\mathrm{\&omega;}})\};$

${\hat{d}}_{i\mathrm{\&omega;}}=\{\hat{d}({\hat{x}}_0,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}),\hat{d}({\hat{x}}_1,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}),\hat{d}({\hat{x}}_2,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}),...,\hat{d}({\hat{x}}_{n_0{n}_1{n}_2{n}_3-1},{\mathrm{\&omega;}}_{i\mathrm{\&omega;}})\},$

Then ${\hat{d}}_{i\mathrm{\&omega;}}=\left\{\begin{array}{cc}{s}_{i\mathrm{\&omega;}}& 0\&le;i\mathrm{\&omega;}\&le;{\mathrm{i\&omega;}}_{max}\\ conj\left(s_{N_t-i\mathrm{\&omega;}+1}\right)& N_t-{\mathrm{i\&omega;}}_{\max }\&le;i\mathrm{\&omega;}\&le;N_t\\ 0& {\mathrm{i\&omega;}}_{max}<i\mathrm{\&omega;}<N_t-{\mathrm{i\&omega;}}_{\max }\end{array}\right.;$

Wherein, function conj () expression asks complex conjugate to each element of input vector, and 0 represents by n ₀n ₁n ₂n ₃the vector that individual 0 element is formed;

Inverse fast Fourier transform is utilized to incite somebody to action the time domain of conversion, obtains five dimension data after interpolation

Wherein, with meet following relational expression, $\hat{d}\left({\hat{x}}_{ix},t_{it}\right)=\underset{i\mathrm{\&omega;}=0}{\overset{N_t-1}{\&Sigma;}}\hat{d}({\hat{x}}_{ix},{\mathrm{\&omega;}}_{i\mathrm{\&omega;}})e^{j2\mathrm{\&pi;}\frac{\left(i\mathrm{\&omega;}\right)\left(it\right)}{N_t}},$ Wherein, it={0,1,2 ..., N _t-1}, ix={0,1,2 ..., n ₀n ₁n ₂n ₃-1}.

8. five dimension interpolation processors of irregular geological data, is characterized in that, comprising:

Geological data conversion module, for by irregular earthquake data transformation to frequency field and spatial domain;

Five dimension interpolation processing modules, for transforming to the frequency slice of the geological data of frequency field and spatial domain for each, according to fast fourier transform algorithm, product calculation is carried out to the matrix of geological data and vector, asks for the frequency-wavenumber domain data after five dimension interpolation;

Five dimension interpolated data conversion modules, for tieing up the frequency-wavenumber domain data transformation after interpolation by five to frequency field and spatial domain, then transform to time domain and spatial domain.

9. five dimension interpolation processors of irregular geological data as claimed in claim 8, is characterized in that, also comprise:

Geological data pretreatment module, for carrying out denoising, static correction and linear NMO process to irregular geological data;

Geological data abstraction module, for extracting the irregular geological data treating five dimension interpolation processing from the irregular geological data after process;

Described geological data conversion module specifically for: by extract treat that the irregular earthquake data transformation of five dimension interpolation processing is to frequency field and spatial domain.

10. five dimension interpolation processors of irregular geological data as claimed in claim 8, is characterized in that, described five dimension interpolation processing modules specifically for:

To each frequency slice i ω=0,1,2 ..., i ω _max, process according to following formula:

b _iω＝F ^Hd _iω；

y _iω＝CG(b _iω,w _iω)；

u _iω＝conj(y _iω)·*y _iω；

w _iω+1＝(sqrt(u _iω/sum(u _iω))+w _iω)/2；

Wherein, i ω _maxbe the upper frequency limit index chosen according to actual needs, selection principle makes geological data d (x _ix, t _it) useful signal energy be just in frequency band in, wherein

W ₀be one and comprise n ₀n ₁n ₂n ₃the column vector of individual element, and w ₀each element equal wherein, w ₀w _{i ω}situation during middle subscript i ω=0;

F=F (x _ix, k _ik), be a N _xrow, n ₀n ₁n ₂n ₃the matrix of row; Wherein, F (x _ix, k _ik) be the discrete inverse Fourier transform matrix of unequal interval, ix={0,1,2 ..., N _x-1} is the line index of matrix, ik={0,1,2 ..., n ₀n ₁n ₂n ₃-1} is matrix column index, wherein, k _ik=(k ₀(i ₀), k ₁(i ₁), k ₂(i ₂), k ₃(i ₃)) _ika point in four-dimentional space, k _η(i _η)=i _η-c _η, c _η=floor ((n _η-1)/2), i _η=0,1,2 ..., n _η-1}, η=0,1,2,3, function f loor () expression rounds downwards, k _ikx _ixrepresent four dimensional vector k _ikwith x _ixinner product, ik and i ₀, i ₁, i ₂, i ₃keep relational expression ik=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃;

$d_{\mathrm{i\&omega;}}=[d(x_0,{\mathrm{\&omega;}}_{\mathrm{i\&omega;}}),d(x_1,{\mathrm{\&omega;}}_{\mathrm{i\&omega;}}),d(x_2,{\mathrm{\&omega;}}_{\mathrm{i\&omega;}}),...,d{(x_{N_{x^{-1}}},{\mathrm{\&omega;}}_{\mathrm{i\&omega;}})]}^H,$ Be one and comprise N _xthe column vector of individual element;

Function conj () expression asks complex conjugate to each element of input vector, and exports a vector; Operational symbol * represents that two vectorial corresponding elements are multiplied, and obtains a vector; Function sqrt () expression is extracted square root to each element of input vector, and exports a vector; Function sum () expression is sued for peace to all elements of input vector, and exports a scalar; Operational symbol ^hrepresent the conjugate transpose asking matrix;

Function y=CG (b, w) is expressed as follows preconditioned conjugate gradient method, 1. sets constant, comprising: iteration ends relative tolerance tol=10 ^-4, maximum iteration time maxit=30, A=F ^hf; 2. calculate iteration initializaing variable, comprising: 3. iteration ends absolute tolerance is calculated for i={0,1,2 ..., maxit-1}, performs and processes operation as follows:

If or i=maxit, jumps out circulation, returns y=y _i, otherwise continue circulation:

z _i＝w·*r _i；

ρ _i＝r _i ^Hz _i；

q _i＝Ap _i；

α _i＝ρ _i/p _i ^Hq _i；

y _i+1＝y _i+α _ip _i；

r _i+1＝r _i-α _iq _i；

p _i+1＝z _i+(ρ _i/ρ _i-1)p _i。

Five dimension interpolation processors of 11. irregular geological datas as claimed in claim 10, is characterized in that,

In described function y=CG (b, w), each iteration needs calculating matrix and vector product computing q=Ap, wherein A=F ^hf has Multilevel Block Toeplitz matrix structure, adopts following fast algorithm to calculate q=Ap:

1. N is established ₀=n ₁n ₂n ₃, N ₁=n ₂n ₃, N ₂=n ₃, N ₃=1, m _η=2n _η-1, η=0,1,2,3, generating length is as follows m _ηvector with wherein η=0,1,2,3:

2. vectorial a={a (i) is established | i=i _a(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., m _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, wherein, and i _a(i ₀, i ₁, i ₂, i ₃)=((i ₀m ₁+ i ₁) m ₂+ i ₂) m ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing, as follows to the element assignment of a, initializing variable

For i ₀=0,1,2 ..., m ₀-1}, performs following circulation (1) to (12)

(1) $t_1^r=t_0^r+i_0^r\left(i_0\right)$

(2) $t_1^c=t_0^c+i_0^c\left(i_0\right)$

(3) for i ₁=0,1,2 ..., m ₁-1}, performs following circulation (4) to (12)

(4) $t_2^r=t_1^r+i_1^r\left(i_1\right)$

(5) $t_2^c=t_1^c+i_1^c\left(i_1\right)$

(6) for i ₂=0,1,2 ..., m ₂-1}, performs following circulation (7) to (12)

(7) $t_3^r=t_2^r+i_2^r\left(i_2\right)$

(8) $t_3^c=t_2^c+i_2^c\left(i_2\right)$

(9) for i ₃=0,1,2 ..., m ₃-1}, performs following circulation (10) to (12)

(10) $t_4^r=t_3^r+i_3^r\left(i_3\right)$

(11) $t_4^c=t_3^c+i_3^c\left(i_3\right)$

(12) $a\left(i_a\right(i_0,i_1,i_2,i_3\left)\right)=A(t_4^r,t_4^c);$

Wherein, i (i) represents i-th element of vectorial i, and a (i) represents i-th element of vectorial a, A (i _r, i _c) representing matrix A i-th _rrow, i _cthe element of row;

3. four-dimensional Fast Fourier Transform (FFT) is done to the four-dimensional array represented by vectorial a, and return vector

Wherein, vector $\hat{a}=\left\{\hat{a}\left(i\right)\right|i=i_{\hat{a}}(i_0,i_1,i_2,i_3),i_{\mathrm{\&eta;}}=0,1,2,...,m_{\mathrm{\&eta;}}-1,\mathrm{\&eta;}=0,1,2,3\}$ That the vectorization of a four-dimensional array represents, $i_{\hat{a}}(i_0,i_1,i_2,i_3)=\left(\right(i_0\&CenterDot;m_1+i_1)\&CenterDot;m_2+i_2)\&CenterDot;m_3+i_3,$ I is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing, function F FT ₄() expression does four-dimensional Fast Fourier Transform (FFT) to four-dimensional array, and returns a four-dimensional array;

4. vectorial p={p (i) is established | i=i _p(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., n _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, wherein, and i _p(i ₀, i ₁, i ₂, i ₃)=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing; If vectorial p '=and p ' (i) | i=i _{p '}(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., m _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, wherein, and i _{p '}(i ₀, i ₁, i ₂, i ₃)=((i ₀m ₁+ i ₁) m ₂+ i ₂) m ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing; In the following manner to each element assignment of vectorial p ':

$p^{\&prime;}\left(i_{p^{\&prime;}}\right(i_0,i_1,i_2,i_3\left)\right)=\left\{\begin{array}{cc}p\left(i_p\right(i_0,i_1,i_2,i_3))& 0\&le;i_{\mathrm{\&eta;}}\&le;n_{\mathrm{\&eta;}}-1\\ 0& n_{\mathrm{\&eta;}}\&le;i_{\mathrm{\&eta;}}\&le;m_{\mathrm{\&eta;}}-1\end{array}\right.;$

5. calculate wherein, function F FT ₄() and IFFT ₄() represents respectively and does four-dimensional Fast Fourier Transform (FFT) and inverse transformation to four-dimensional array, and returns a four-dimensional array; Operational symbol * represents that the corresponding element of two four-dimensional arrays is multiplied, and obtains a four-dimensional array;

Vector q '=and q ' (i) | i=i _{q '}(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., m _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, wherein, and i _{q '}(i ₀, i ₁, i ₂, i ₃)=((i ₀m ₁+ i ₁) m ₂+ i ₂) m ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing;

6. vectorial q={q (i) is established | i=i _q(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., n _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, wherein, and i _q(i ₀, i ₁, i ₂, i ₃)=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing, in the following manner to each element assignment of vectorial q, q (i _q(i ₀, i ₁, i ₂, i ₃))=q ' (i _{q '}(i ₀+ n ₀, i ₁+ n ₁, i ₂+ n ₂, i ₃+ n ₃)), wherein 0≤i _η≤ n _η-1, η=0,1,2,3.

Five dimension interpolation processors of 12. irregular geological datas as claimed in claim 11, is characterized in that, also comprise: to y _{i ω}perform function wherein, i ω=0,1,2 ..., i ω _max;

Function represent and the element position of four-dimensional array y is adjusted, and return a four-dimensional array wherein:

Y={y (i) | i=i _y(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., n _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, i _y(i ₀, i ₁, i ₂, i ₃)=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing;

$\hat{y}=\left\{\hat{y}\left(i\right)\right|i=i_{\hat{y}}(i_0,i_1,i_2,i_3),i_{\mathrm{\&eta;}}=0,1,2,...,n_{\mathrm{\&eta;}}-1,\mathrm{\&eta;}=0,1,2,3\}$ That the vectorization of a four-dimensional array represents, $i_{\hat{y}}(i_0,i_1,i_2,i_3)=\left(\right(i_0\&CenterDot;n_1+i_1)\&CenterDot;n_2+i_2)\&CenterDot;n_3+i_3,$ I is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing;

$\hat{y}\left(i_{\hat{y}}\right({\hat{i}}_0,{\hat{i}}_1,{\hat{i}}_2,{\hat{i}}_3\left)\right)=y\left(i_y\right(i_0,i_1,i_2,i_3\left)\right),$ Wherein ${\hat{i}}_{\mathrm{\&eta;}}=\mathrm{mod}(i_{\mathrm{\&eta;}}-c_{\mathrm{\&eta;}},n_{\mathrm{\&eta;}}),$ C _η=floor ((n _η-1)/2), i _η=0,1,2 ..., n _η-1}, η=0,1,2,3, function mod (v ₁, v ₂) represent v ₁to v ₂ask 0 ~ v ₂remainder between-1, function f loor () expression rounds downwards.

Five dimension interpolation processors of 13. irregular geological datas as claimed in claim 12, is characterized in that, also comprise:

Right perform function wherein, i ω=0,1,2 ..., i ω _max;

Function represent four-dimensional array do four-dimensional inverse fast Fourier transform, and return a four-dimensional array s, wherein,

$\hat{y}=\left\{\hat{y}\left(i\right)\right|i=i_{\hat{y}}(i_0,i_1,i_2,i_3),i_{\mathrm{\&eta;}}=0,1,2,...,n_{\mathrm{\&eta;}}-1,\mathrm{\&eta;}=0,1,2,3\}$ That the vectorization of a four-dimensional array represents, $i_{\hat{y}}(i_0,i_1,i_2,i_3)=\left(\right(i_0\&CenterDot;n_1+i_1)\&CenterDot;n_2+i_2)\&CenterDot;n_3+i_3,$ I is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing;

S={s (i) | i=i _s(i ₀, i ₁, i ₂, i ₃), i _η=0,1,2 ..., n _η-1, η=0,1,2,3} is that the vectorization of a four-dimensional array represents, i _s(i ₀, i ₁, i ₂, i ₃)=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃, i is vector index, i ₀, i ₁, i ₂, i ₃for four-dimensional array indexing.

Five dimension interpolation processors of 14. irregular geological datas as claimed in claim 13, is characterized in that, described five dimension interpolated data conversion modules specifically for:

Utilize s _{i ω}generate frequency field five dimension data after interpolation wherein, for the spatial sampling coordinate after interpolation, ix={0,1,2 ..., n ₀n ₁n ₂n ₃-1}, i ω=0,1,2 ..., N _t-1},

${\hat{x}}_{ix}={\left(x_0\right(i_0),x_1(i_1),x_2(i_2),x_3(i_3\left)\right)}_{ix},$ Wherein, $x_{\mathrm{\&eta;}}\left(i_{\mathrm{\&eta;}}\right)={\mathrm{vmin}}_{\mathrm{\&eta;}}+i_{\mathrm{\&eta;}}\mathrm{\&Delta;}{\tilde{x}}_{\mathrm{\&eta;}},$ I _η=0,1,2 ..., n _η-1}, η=0,1,2,3, for the sampling interval of each space dimension after interpolation, ix and i ₀, i ₁, i ₂, i ₃keep relational expression ix=((i ₀n ₁+ i ₁) n ₂+ i ₂) n ₃+ i ₃; generating mode is as follows, with vector representative a part of data,

${\hat{d}}_{i\mathrm{\&omega;}}=\{\hat{d}({\hat{x}}_0,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}),\hat{d}({\hat{x}}_1,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}),\hat{d}({\hat{x}}_2,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}),...,\hat{d}({\hat{x}}_{n_0{n}_1{n}_2{n}_3-1},{\mathrm{\&omega;}}_{i\mathrm{\&omega;}})\};$

${\hat{d}}_{i\mathrm{\&omega;}}=\{\hat{d}({\hat{x}}_0,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}),\hat{d}({\hat{x}}_1,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}),\hat{d}({\hat{x}}_2,{\mathrm{\&omega;}}_{i\mathrm{\&omega;}}),...,\hat{d}({\hat{x}}_{n_0{n}_1{n}_2{n}_3-1},{\mathrm{\&omega;}}_{i\mathrm{\&omega;}})\},$

Then ${\hat{d}}_{i\mathrm{\&omega;}}=\left\{\begin{array}{cc}{s}_{i\mathrm{\&omega;}}& 0\&le;i\mathrm{\&omega;}\&le;{\mathrm{i\&omega;}}_{max}\\ conj\left(s_{N_t-i\mathrm{\&omega;}+1}\right)& N_t-{\mathrm{i\&omega;}}_{\max }\&le;i\mathrm{\&omega;}\&le;N_t\\ 0& {\mathrm{i\&omega;}}_{max}<i\mathrm{\&omega;}<N_t-{\mathrm{i\&omega;}}_{\max }\end{array}\right.;$

Wherein, function conj () expression asks complex conjugate to each element of input vector, and 0 represents by n ₀n ₁n ₂n ₃the vector that individual 0 element is formed;

Inverse fast Fourier transform is utilized to incite somebody to action the time domain of conversion, obtains five dimension data after interpolation

Wherein, with meet following relational expression, $\hat{d}\left({\hat{x}}_{ix},t_{it}\right)=\underset{i\mathrm{\&omega;}=0}{\overset{N_t-1}{\&Sigma;}}\hat{d}({\hat{x}}_{ix},{\mathrm{\&omega;}}_{i\mathrm{\&omega;}})e^{j2\mathrm{\&pi;}\frac{\left(i\mathrm{\&omega;}\right)\left(it\right)}{N_t}},$ Wherein, $j=\sqrt{-1},t_{it}=it\&CenterDot;\mathrm{\&Delta;}t,$ it＝{0,1,2,…,N _t-1}，ix＝{0,1,2,…,n ₀n ₁n ₂n ₃-1}。